Дерево решений (Decision Tree)

[Igg]

Здравствуйте !
Не могли бы уважаемые участники поделится опытом создания decision tree для растровых данных в R, по аналогии с тем что реализовано в Erdas или ENVI ?

зы: гугление особых успехов не принесло)

[Максим Дубинин]

обрисовываю в общих чертах, без углубления в детали, это выборка из более сложного моего кода, поэтому какие-то детали могут выпасть, но суть должна быть понятна, комментарии по коду

код использует куски из статьи про то как в R работать с растрами

Код: Выделить всё

#понадобится два модуля, их надо будет перед выполнением установить
library(rgdal)
library(tree)

#тренировочные данные
traindata = matrix(as.integer(scan(file=trainname)),ncol=numbands+1,byrow=T)
#название класса (у меня - один из параметров тренировочных данных)
class=factor(traindata[,numbands+1])
#немного переделаем, поставим класс первым параметром
trainframe=data.frame(class,traindata[,1:numbands])

#если нужно, делаем выборку с возвращением, скажем 50% из тренировочного набора
numrows=length(trainframe[,1])
indexAll=1:numrows
numtrees=30
step=0.5
sampledata=sample(indexAll,numrows*step,replace=T)
samplesarray = cbind(samplesarray,sampledata)
sampledata=samplesarray[,i]
testdata=trainframe[sampledata,]

#наконец получили класс есть, данные есть, строим дерево
atree=tree(class ~.,data=testdata,mindev=0,minsize=1)

#часть 2, дерево построили, теперь надо натравить его на сам растр и его отклассифицировать
image = ... #тут указываем само изображение
#получаем его размерности на будущее
x = new("GDALReadOnlyDataset", image)
width = dim(x)[2]
height = dim(x)[1]
#пропускаем первые два параметра, это lat,long, их классифицировать не надо
imagedata = data.frame(getRasterTable(x)[3:(numbands+2)])

#наконец классифицируем
B=predict.tree(atree,newdata=imagedata,type="vector")

#сохраняем результат обратно в новое изображение
tif_driver <- new("GDALDriver", "GTiff")
tif2 <- new("GDALTransientDataset", tif_driver, height, width, 1, 'Int8')
resmtx = matrix(B,width,height)
bnd1 <- putRasterData(tif2, resmtx)
tif_file = "D:\\Programming\\R\\trees\\simple-tree-from-raster\\s_2000res.tif"
saveDataset(tif2, tif_file)
GDAL.close(tif2)
GDAL.close(tif_driver)

[Анна]

а можно разобрать вот эту строку подробнее:

Код: Выделить всё

#название класса (у меня - один из параметров тренировочных данных)
class=factor(traindata[,numbands+1])

где здесь название класса? почему нужно использовать factor? что этим действием достигается?

[Максим Дубинин]

название класса берется в данном случае из последнего столбца

Код: Выделить всё

traindata[,numbands+1]

а factor, так как у меня в примере классы были 1 и 2, фактор переведет их в строки "1" и "2", иначе получится дерево регрессий, а не дерево решений, правда результат все равно будет один и тот же, я подозреваю, но не факт

[Анна]

спасибо большое!
а еще вопрос - в каком формате находится загружаемая матрица? это tif, img, grid или еще что-то (текстовый файл, таблица)?
и почему если это многоканальный растр, нельзя просто сказать numbands, а надо numbands+1?

[Максим Дубинин]

можно tif, можно img, что засунешь, то rgdal и съест
numbands+1 это потому что в тренингах поле class - следующее после собственно каналов, а numbands берется из самих растров

[Анна]

а что class не является еще одним каналом? как же он представлен в этом многоканальном растре?
я себе представляла матрицу тренингов чем-то вроде:

Код: Выделить всё

#пикселя   1b       2b    3b   class
1          234      123   56     1
2          230      131   46     1
3          134      223   16     2

где класс - это еще один канал, просто содержащий либо значения 1 либо 2. или что-то я недопонимаю?

[Максим Дубинин]

да, class это еще один "канал", а точнее поле-атрибут, но канал набора треннингов, а не растра
он никак в растре не представлен, растр - это растр, тренинги это данные с растра + класс

[Анна]

ок, тогда следующий вопрос - как представлены тренинги? какой-то особый тип данных\формат? и как их получить и привести в такой формат в R?

[Максим Дубинин]

Код: Выделить всё

234      123   56     1
230      131   46     1
134      223   16     2

Например так (канал1, канал2, канал3, класс). Ничего особенного. Я получаю их вот так.

[Анна]

все поняла спасибо. Это все просто в текстовом формате, а не в растре:) спасибо огромное!!

[Анна]

продолжаю мучать вопросами

Код: Выделить всё

samplesarray = cbind(samplesarray,sampledata)

в данном случае, как может в cbind использоваться samplesarray если она только инициализируется на этом шаге?

[Максим Дубинин]

а это мой глюк, возник когда я кусок кода вырезал из скрипта, дело в том, что для бэггинга, делается популяция выборок и все они запихиваются в массив

Код: Выделить всё

samplesarray = array()
#create samples array
for(i in 1:numtrees) {
	sampledata=sample(indexAll,numrows*step,replace=T)
	samplesarray = cbind(samplesarray,sampledata)
}
samplesarray = samplesarray[,-1]

А потом по одному берутся из samplesarray для создания моделей.