# 
# If "?????????" then File - Reopen with Encoding... windows-1251
#
# МГУ имени М.В. Ломоносова, географический факультет
# кафедра физической географии и ландшафтоведения
# Д.Н. Козлов, Н.И.Лозбенев
# дисциплина 4 курса по выбору "Анализ данных в среде R"
# http://www.landscape.edu.ru/edu_help4_R.shtml
# 
# упражнение "Методы и алгоритмы цифровой почвенной картографии –
# модели почвенно-ландшафтных связей для категорий номинальной шкалы" 
# 
# объект исследования - почвенный покров Стрелецкой степи Центрально-черноземного заповедника
# 
# СОДЕРЖАНИЕ ЗАДАНИЯ
# 1. подготовка данных и рабочей среды R
# 2. анализ признакового пространства
# 3. моделирование почвенно-ландшафтных связей (ПЛС)
#  3.1. линейный дискриминантный анализ
#  3.2. ансамбль деревьев решений
#  3.3. метод опорных векторов
# 4. экспорт результатов моделирования
# 5. контрольные вопросы и задания для самостоятельной работы

# ПОДГОТОВКА ДАННЫХ И РАБОЧЕЙ СРЕДЫ R
##########################################################################################
# задача 1.1: скачать набор данных
##########################################################################################
# 
# распаковать в папку на диске архив с учебным набором данных 
# www.landscape.edu.ru/files/DSM_METdataset_v0.99.zip 
# либо https://yadi.sk/i/xEJB1-9DyeqUC
#
# актуальная версия скрипта R - ...
#
##########################################################################################
# задача 1.2: установить и активировать необходимые библиотеки R
##########################################################################################

library(raster)
library(car)
library(leaflet)
library(rasterVis)
library(randomForest) 
library(ithir)
library(ellipse)
library(MASS)
library(relaimpo)
library(e1071)

# при первом обращении установить пакеты
install.packages('car', dependencies = TRUE)
install.packages('leaflet', dependencies = TRUE)
install.packages('raster', dependencies = TRUE)
install.packages('ellipse', dependencies = TRUE)
install.packages('randomForest', dependencies = TRUE)
install.packages("relaimpo", dependencies = TRUE)
install.packages("rasterVis", dependencies = TRUE)
install.packages("relaimpo", dependencies = TRUE)
install.packages("MASS", dependencies = TRUE)
install.packages("e1071", dependencies = TRUE)

install.packages("devtools") # необходим для установки пакета "ithir" из репозитария Брендана Мелони
library(devtools)
install_bitbucket("brendo1001/ithir/pkg", dependencies = TRUE)
#
##########################################################################################
# задача 1.3: загрузить и визуализировать данные
##########################################################################################
#
setwd('C:/Dropbox/0_EDU/R/METsid') # назначить рабочую директорию 
setwd('D:/Dropbox/0_EDU/R/METsid') # назначить рабочую директорию 
setwd("C:/Users/User/YandexDisk/R/MET")

# загрузить растры формата SAGA-grid и определить их проекцию
covStack = stack(list.files(path = "...", pattern='*.sdat'))
proj4string(covStack) = CRS('+init=epsg:32637') # (UTM (WGS84), зона 37N)
# присвоить ковариатам компактные имена
names(covStack) = c("z","slp","cdep","facc","twi","tpi10","tpi60","plncurv","prfcurv","gnrcurv", "cscurv","aacn", "simwe") 

par(mar=c(2, 2, 1.5, 0.2))  # поля графика снизу по часовой стрелке, см
plot(covStack[[1]],  main = "Абс. высота, м", col=terrain.colors(255))
plot(covStack[[2]],  main = "Крутизна, °", col=grey(seq(0,1,len = 25)))
plot(covStack[[3]],  main = "Глубина замкнутых понижений, м", col=colorRampPalette(c("lightgrey", "red"))( 10) )
plot(covStack[[4]],  main = "Водосборная площадь, м2", col=colorRampPalette(c("lightgrey", "darkblue"))( 50))
plot(covStack[[5]],  main = "ТопоИндекс влажности, log(CA/slope)")
plot(covStack[[6]],  main = "Отн. превышения в окрестности 10 м, м", col=colorRampPalette(c("dodgerblue4", "lightgrey", "darkorange3")) (50), zlim = c(-0.3,0.3))
plot(covStack[[7]],  main = "Отн. превышения в окрестности 60 м, м", col=colorRampPalette(c("dodgerblue4", "lightgrey", "darkorange3")) (50), zlim = c(-1.0,1.0))
plot(covStack[[8]],  main = "Плановая кривизна", col=colorRampPalette(c("dodgerblue4", "lightgrey", "darkorange3")) (50), zlim = c(-1.0,1.0))
plot(covStack[[9]],  main = "Профильная кривизна", col=colorRampPalette(c("darkblue", "lightgrey", "red"))( 50))
plot(covStack[[10]], main = "Общая кривизна", col=colorRampPalette(c("darkblue", "lightgrey", "red"))( 50))
plot(covStack[[11]], main = "Кросс-секционная кривизна", col=colorRampPalette(c("darkblue", "lightgrey", "red"))( 50))
plot(covStack[[12]], main = "Превышение над водотоком, м", col=colorRampPalette(c("lightgrey", "red"))( 30))
plot(covStack[[13]], main = "Слой перераспределенных осадков, м", col=colorRampPalette(c("lightgrey", "darkblue"))( 50))

# добавить горизонтали поверх изображения (комбинировать по желанию)
contour(covStack[[1]], drawlabels = TRUE, axes = FALSE, frame = TRUE, add=TRUE, col="black")

# загрузить таблицу с точками почвенного опробования
DSM_data = read.csv2("met_soilpoints.txt", header = TRUE, sep = "\t", 
                     quote = "\"", dec = ".", stringsAsFactors = FALSE) 
str(DSM_data) # показать структуру таблицы
DSM_data$sid = factor(DSM_data$sid,ordered=TRUE, levels=c("Лч","Чв","Чт","Чтк")) # индекс почв в порядковую шкалу
table(DSM_data$sid) # представленность почв (число точек описаний)

# преобразовать строки таблицы в объекты-точки с координатами Xm, Ym 
coordinates(DSM_data) = ~ Xm + Ym # функция пакета "raster"
proj4string(DSM_data) = CRS('+init=epsg:32637') # определить проекцию точек (UTM (WGS84), зона 37N)

# добавить точки (элементы обучающей выборки) к активной карте 
points(DSM_data, pch = 21:21, cex=0.5, col="darkslateblue") 
# и подписать индекс почвы
text(DSM_data$Xm, DSM_data$Ym, labels=DSM_data$sid, cex = 0.7, pos = 3, offset = 0.5)

# используем Leaflet map widget для визуализации точек и ковариат (необходимо подключение к Интернету)
leaflet() %>% 
  addMarkers(lng = 36.09170, lat = 51.57200, popup = "Метеостанция") %>% 
  addProviderTiles("Esri.WorldImagery") 

leaflet() %>% 
  addTiles() %>% 
  addProviderTiles("Esri.WorldImagery") %>% 
  addCircleMarkers(data = spTransform(DSM_data, CRS('+init=epsg:4326')), 
                  popup = paste0("<strong>sid: </strong>", DSM_data$sid), radius = 2, color = "red")

# добавим в виджет растр ковариат (выбор по номеру слоя из covStack)
covRST = covStack[[6]] # [6] - относительные превышения в окрестности 10 м
covRST = projectRaster(covRST, crs = '+init=epsg:4326') # перепроецировать растр в проекцию дол/шир (WGS84)
pal = colorNumeric(palette=rev(brewer.pal(3,"RdBu")), domain = values(covRST), na.color = "transparent") # сформировать палитру
leaflet() %>% 
  addProviderTiles("Esri.WorldImagery") %>% 
  addRasterImage(covRST, colors = pal, opacity = 0.99) %>%
  addLegend(pal = pal, values = values(covRST), title = "topographic wetness index") %>%
  addCircleMarkers(data = spTransform(DSM_data, CRS('+init=epsg:4326')), popup = paste0("<strong>sid: </strong>", DSM_data$sid), radius = 2, color = "black")
# повторите для других свойств рельефа, поигравшись с палитрой и прозрачностью


################################################################################################################
# 2. АНАЛИЗ ПРИЗНАКОВОГО ПРОСТРАНСТВА (обосновать наиболее информативные признаки)
################################################################################################################
#
##########################################################################################
# задача 2.1: извлечь значения ковариат для точек почвенных описаний
##########################################################################################
#
DSM_data = data.frame(extract(covStack, DSM_data, method="simple", sp=1))
DSM_data = data.frame(DSM_data[ ,c("sid","z","slp","cdep","facc","twi","tpi10", "tpi60","plncurv","prfcurv","gnrcurv","cscurv","aacn", "simwe")])
format(head(DSM_data), digits=1, scientific=FALSE) # первые 6 строк
#
##########################################################################################
# задача 2.2: проверить независимость переменных
##########################################################################################
#
par(mar=c(0, 0, 0, 0))  # поля графика снизу по часовой стрелке, см
round(cor(DSM_data[ ,-1]), 2)   # рассчитать матрицу корреляций
symnum(cor(DSM_data[ ,-1]))  # компактный образ матрицы корреляций
plotcorr(cor(DSM_data[ ,-1]), type="lower")# образ в виде эллипсов
#
##########################################################################################
# задача 2.3: оценить значимость переменных в разделении почв
##########################################################################################
#
# визуальный анализ различимости категорий почв в пространстве признаков library(car)
scatterplotMatrix(~ slp + simwe + log(facc) + twi|sid, data=DSM_data, cex=0.8) 
scatterplotMatrix(~ tpi10 + twi + log(facc)|sid, data=DSM_data, cex=0.8)
# комбинируйте ковариаты в поиске мультипликативного эффекта
scatterplotMatrix(~ twi + log(facc*(tpi10+1)) + log(facc)| sid, data=DSM_data, cex=0.8)

# формальная оценка значимости переменных
lmMod = lm(as.numeric(sid) ~ ., data = DSM_data)  # подобрать линейную модель, разделяющую подтипы почв в факторном пространстве
relImportance = calc.relimp(lmMod, type = "lmg", rela = TRUE)  # рассчитать значимость переменных в долях от 100
sort(round(relImportance$lmg, 2), decreasing=TRUE)  # относительная значимость переменных (relative importance)

# Повторить для избранных переменных
lmMod = lm(as.numeric(sid) ~ twi + tpi10 + cdep + cscurv + aacn, data = DSM_data)  # подобрать линейную модель
relImportance = calc.relimp(lmMod, type = "lmg", rela = TRUE)  # рассчитать значимость переменных в долях от 100
sort(round(relImportance$lmg, 2), decreasing=TRUE)  # относительная значимость переменных (relative importance)

# Диаграмма рассеяния (точки опробования в пространстве признаков)
# исследуйте проекции многомерного признакового пространства почв на двухмерную плоскость, меняя переменные Х и У
par(mar=c(4, 4, 1, 0.2))  # поля графика снизу по часовой стрелке, см
# сформировать палитру (HEX-цвета почв "Лч", "Чв", "Чт", "Чтк")
sid_col = c("#92DC99", "#D39829", "#B7621C", "#793405") 
palette(sid_col)
# все пиксели в пространстве двух ковариат
plot(as.vector(covStack$tpi10), as.vector(covStack$twi), col="grey", pch=20, cex=0.6, 
     xlim = c(-0.2,0.2), ylim = c(4,16),  xlab="отн. превышения в окрестности 10 м, м", ylab="топоиндекс влажности") 
# отметить пиксели, опробованные при почвенном обследовании
points(DSM_data$tpi10, DSM_data$twi, col=DSM_data$sid, lwd=2, cex=1.2) 
# добавить эллипс доверительного интервала таксонов почв 50 и 90%
dataEllipse(DSM_data$tpi10, DSM_data$twi, DSM_data$sid, c(0.5,0.95), add=TRUE, 
            plot.points=TRUE, lwd=1, fill=TRUE, fill.alpha=0.15, pch=c(21,21,21,21), col=sid_col)
# добавить легенду
legend("topright",bg="white",pch=c(21,21,21,21), legend=c("Лч","Чв","Чт", "Чтк"), 
       inset = .02, bty="n", col= sid_col, cex=1.2, pt.lwd=2, pt.bg="white")

# ВЫВОД: области разных почв в пространстве факторов пересекаются

####################################################################################################
# 3. Моделирование почвенно-ландшафтных связей (ПЛС)
####################################################################################################
#
##########################################################################################
# задача 3.1: подготовка обучающей выборки
##########################################################################################
# разбиение массива точек на 2 части
# 70% - обучение и проверка точности; 30% - верификация (независимая проверка точности)
set.seed(123) # способ отбора случайной выборки
training <- sample(nrow(DSM_data), 0.70 * nrow(DSM_data)) # выбор 70% элементов выборки
training # выводит номера элементов обучающей выборки
DSM.train = DSM_data[training, ]  # обучающая выборка
DSM.verif = DSM_data[-training, ] # выборка для верификации
#
##########################################################################################
# задача 3.2: дискриминантная модель ПЛС (Linear Discriminant Analysis (lda)
##########################################################################################
# 
# построение дискриминантной модели (библиотека MASS) 
mod.lda = lda(sid ~ tpi10 + twi, data=DSM_data)  # с априорными вероятностями пропорционально размеру группы
mod.lda # структура модели

# модель с экспертно оцененными априорными вероятностями
mod.lda = lda(sid ~ tpi10 + twi, data=DSM_data, prior=c(0.15,0.30,0.40,0.15)) 
mod.lda

# наиболее вероятная категория для точек обучающей выборки
pred.lda = as.data.frame(predict(mod.lda)) # 
format(head(pred.lda), digits=1, scientific=FALSE) # наиболее вероятный класс, вероятности, дискриминантные функции

# оценка качества модели: матрица соответствия (строки - предсказанные, колонки - исходные), пакет (ithir)
goofcat(observed = DSM_data$sid, predicted = pred.lda$class)  

# наиболее вероятная категория для ячеек (пикселей) растра
map.lda.c = predict(covStack, mod.lda) # наиболее вероятный класс для каждого пикселя
# Вероятность принадлежности пикселя к данному классу 
map.lda.p1 = predict(covStack, mod.lda, type = "probs", index=2) # Лч
map.lda.p2 = predict(covStack, mod.lda, type = "probs", index=3) # Чв
map.lda.p3 = predict(covStack, mod.lda, type = "probs", index=4) # Чт
map.lda.p4 = predict(covStack, mod.lda, type = "probs", index=5) # Чтк
map.lda.p4 = predict(covStack, mod.lda, type = "probs", index=5) # Чтк
# объединим вероятности в общий стек
map.lda.p   = stack(map.lda.p1, map.lda.p2, map.lda.p3, map.lda.p4) 
names (map.lda.p) = c('Лч', 'Чв', 'Чт', 'Чтк') # переименуем

# модель ПЛС в признаковом пространстве (пакет 'car')
par(mar=c(4, 4, 1, 0.2))  # поля графика снизу по часовой стрелке, см
# сформировать палитру (HEX-цвета почв по порядку "Лч", "Чв", "Чт", "Чтк")
sid_col = c("#5AA75F", "#B77E29", "#A25817", "#793405") # цвета эллипсов
palette(c("#79C17E", "#E0B462", "#D89156", "#934D08"))  # цвета областей
# все пиксели в пространстве двух ковариат
plot(as.vector(covStack$tpi10), as.vector(covStack$twi), col=as.vector(map.lda.c), pch=20, cex=0.6, xlim = c(-0.2,0.2),
     ylim = c(4,16),  xlab="отн. превышения в окрестности 10 м, м", ylab="топоиндекс влажности") 
# добавить эллипс доверительного интервала таксонов почв 50 и 90%
dataEllipse(DSM_data$tpi10, DSM_data$twi, DSM_data$sid, c(0.5,0.95), add=TRUE, plot.points=TRUE, lwd=1,
            fill=TRUE, fill.alpha=0.15, pch=c(21,22,24,25), col=sid_col)
# добавить легенду
legend("topright",bg="white", legend=c("Лч","Чв","Чт", "Чтк"), inset = .02, 
       bty="n", col= sid_col, pch=c(21,22,24,25), cex=1.2, pt.lwd=2, pt.bg="white")

# модель ПЛС в территориальном пространстве (прогнозная карта)
cat = levels(as.factor(map.lda.c))[[1]]    # список числовых кодов атрибутивной таблицы растра
cat[["sid"]] = c("Лч","Чв", "Чт","Чтк")  # их расшифровка (легенда)
levels(map.lda.c) = cat      # добавить легенду к аттрибутам растра (формируем легенду)
sid_col = c("#92DC99", "#E6BE75", "#D89156", "#934D08") # цвета HEX для каждой категории по порядку
levelplot(map.lda.c, col.regions=sid_col, xlab="", ylab="") # построить изображение (пакет rasterVis)

# Построить карты постериорных вероятностей
par(mfrow=c(1,4)) # разделить по 4 карты в один ряд, отменить - dev.off() 
par(mar=c(2.5, 0, 1, 2))  # поля графика снизу по часовой стрелке, см
for (i in 1:4) plot(map.lda.p[[i]], col = gray.colors(100, start = 0.9, end = 0.1))  

# Неопределенность прогноза наиболее вероятного таксона 
map.lda.pmax = stackApply(map.lda.p, indices=c(1,1,1,1), fun=max) 
cellStats(map.lda.pmax, mean) # средняя определенность прогноза
par(mfrow=c(1,1)) # восстановим отрисовку карт по одной в окне 
plot(map.lda.pmax, col = grey(seq(0, 1, len = 25)))

# обучение по части обучающей выборки (70%)
mod.lda.trn = lda(sid ~ tpi10 + twi, data=DSM.train) 
# расчет модельных значений
pred.lda.trn = as.data.frame(predict(mod.lda.trn)) # 70%
pred.lda.ver = as.data.frame(predict(mod.lda.trn, DSM.verif)) # 30%
# оценка точности и верификация модели 
goofcat(observed = DSM.train$sid, predicted = pred.lda.trn$class) # 70% 
goofcat(observed = DSM.verif$sid, predicted = pred.lda.ver$class) # 30%
#
##########################################################################################
# задача 3.3: модель ПЛС на основе ансамбля деревьев решений
##########################################################################################
# 
# модель с настройками по умолчанию
mod.rf  = randomForest(sid~tpi10 + twi, data=DSM_data) 
# Определение оптимального числа деревьев по графику
mod.rf  = randomForest(sid~tpi10 + twi, data=DSM_data, ntree=300, mtry=1, sampsize=50, nodesize=5, replace=TRUE) 

# Определение оптимального числа деревьев по графику
plot(mod.rf, col=c("black", "#92DC99", "#E6BE75", "#D89156", "#934D08"))

print(mod.rf) # показать основные результаты
importance(mod.rf) # относительная важность каждой ковариаты
varImpPlot(mod.rf) # то же в графической форме

# Предсказание подтипа почв почв для исходных точек на основе свойств рельефа
pred.rf = as.data.frame(predict(mod.rf)) # Наиболее вероятный класс
head(pred.rf)

# оценка качества модели
goofcat(observed = DSM_data$sid, predicted = pred.rf) 

# наиболее вероятная категория для элементов ОБУЧАЕМОЙ выборки (растра)
map.rf.c = predict(covStack, mod.rf) # наиболее вероятный класс
map.rf.p1 = predict(covStack, mod.rf, type = "prob", index=1) # Лч
map.rf.p2 = predict(covStack, mod.rf, type = "prob", index=2) # Чв
map.rf.p3 = predict(covStack, mod.rf, type = "prob", index=3) # Чт
map.rf.p4 = predict(covStack, mod.rf, type = "prob", index=4) # Чтк
# объединим вероятности в общий стек
map.rf.p  = stack(map.rf.p1, map.rf.p2, map.rf.p3, map.rf.p4) 
names (map.rf.p) = c('Лч', 'Чв', 'Чт', 'Чтк') # переименуем
par(mfrow=c(1,4)) # разместить карты по четыре в один ряд
par(mar=c(2.5, 0, 1, 2))  # поля графика снизу по часовой стрелке, см
for (i in 1:4) plot(map.rf.p[[i]], col = gray.colors(100, start = 0.9, end = 0.1))
# построить диаграммы рассеяния постериорных вероятностей lda и rf моделей
for (i in 1:4) plot(map.lda.p[[i]], map.rf.p[[i]])

# Неопределенность прогноза наиболее вероятного таксона 
par(mfrow=c(1,1)) # восстановим отрисовку карт по одной в окне 
map.rf.pmax = stackApply(map.rf.p, indices=c(1,1,1,1), fun=max) 
cellStats(map.rf.pmax, mean) # средняя определенность прогноза
plot(map.rf.pmax, col = grey(seq(0, 1, len = 25)))

# модель ПЛС в признаковом пространстве 
par(mfrow=c(1,1), mar=c(4, 4, 1, 0.2)) # отменить деление графического окна на части и поля графика снизу по часовой стрелке, см
# сформировать палитру (HEX-цвета почв по порядку "Лч", "Чв", "Чт", "Чтк")
sid_col = c("#5AA75F", "#B77E29", "#A25817", "#793405") # цвета эллипсов
palette(c("#79C17E", "#E0B462", "#D89156", "#934D08")) # цвета областей
plot(as.vector(covStack$tpi10), as.vector(covStack$twi), col=as.vector(map.rf.c), pch=20, cex=0.6, xlim = c(-0.2,0.2), ylim = c(4,16), 
     xlab="отн. превышения в окрестности 10 м, м", ylab="топоиндекс влажности")
dataEllipse(DSM_data$tpi10, DSM_data$twi, DSM_data$sid, c(0.5, 0.95), add = TRUE, plot.points = TRUE, lwd=1, 
            fill=TRUE, fill.alpha=0.15, pch=c(21,22,24,25), col=sid_col)
legend("topright",bg="white",pch=c(21,22,24,25), legend=c("Лч","Чв","Чт", "Чтк"), inset = .02, bty="n",
       col=sid_col,pt.bg=c("white","white","white","white"))

# Визуализируем карту
cat <- levels(as.factor(map.rf.c))[[1]] # список числовых кодов атрибутивной таблицы растра
cat[["sid"]] <- c("Лч","Чв", "Чт","Чтк")  # их расшифровка (легенда)
levels(map.rf.c) <- cat                 # добавить легенду к аттрибутам растра
area_colors <- c("#92DC99", "#E6BE75", "#D89156", "#934D08")      # цвета HEX для каждой категории по порядку
levelplot(map.rf.c, col.regions=area_colors, xlab="", ylab="")  # построить изображение

# обучение по части обучающей выборки (70%)
mod.rf.trn = randomForest(sid~tpi10+twi, data=DSM.train, ntree=100,mtry=2) 
# расчет модельных значений
pred.rf.trn = as.data.frame(predict(mod.rf.trn)) # 70%
pred.rf.ver = as.data.frame(predict(mod.rf.trn, DSM.verif)) # 30%
# оценка точности и верификация модели 
goofcat(observed = DSM.train$sid, predicted = pred.rf.trn$predict) # 70% 
goofcat(observed = DSM.verif$sid, predicted = pred.rf.ver$predict) # 30%
#
##########################################################################################
# задача 3.4: модель ПЛС методом опорных векторов (Support Vector Machine или SVM) 
##########################################################################################
#
# модель с начальными настройками (по умолчанию) 
mod.svm = svm(sid ~ tpi10 + twi, data = DSM_data) # библиотека е1071
summary(mod.svm)
goofcat(observed=DSM_data$sid,predicted=predict(mod.svm,newdata=DSM_data)) # оценка точности модели 

# подбор оптимальных настроек
svmtunepar = tune.svm(sid ~ tpi10 + twi, data = DSM_data, gamma = seq(0.25,2,0.25), cost = seq(1,501,25)) 
svmtunepar = tune.svm(sid ~ tpi10 + twi, data = DSM_data, gamma = seq(0.25,2,0.25), cost = seq(1,51,5), random=FALSE) 
summary(svmtunepar)
plot (svmtunepar,  xlab="gamma", ylab="cost")
svmtunepar$ # $best.parameters

# модель с экспертно заданными априорными вероятностями
costs = table(DSM_data$sid) 
costs['Лч']   = 10.20
costs['Чв']   = 1.25 
costs['Чт']   = 1.35
costs['Чтк']  = 1.20 
costs # показать априорные вероятности
mod.svm = svm(sid ~ tpi10 + twi,data = DSM_data, class.weights=costs, cost = 6, gamma=1) 
summary(mod.svm)
mod.svm
# оценка точности модели 
goofcat(observed=DSM_data$sid,predicted=predict(mod.svm,newdata=DSM_data))

# наиболее вероятный класс (прогнозная карта)
map.svm.c=predict(covStack, mod.svm, type = "class") 

# модель ПЛС в признаковом пространстве
par(mfrow=c(1,1), mar=c(4, 4, 1, 0.2)) # отменить деление графического окна на части и поля графика снизу по часовой стрелке, см
# сформировать палитру (HEX-цвета почв по порядку "Лч", "Чв", "Чт", "Чтк")
sid_col = c("#5AA75F", "#B77E29", "#A25817", "#793405") # цвета эллипсов
palette(c("#79C17E", "#E0B462", "#D89156", "#934D08"))  # цвета областей
plot(as.vector(covStack$tpi10), as.vector(covStack$twi), col=as.vector(map.svm.c), pch=20, cex=0.6, xlim = c(-0.2,0.2), ylim = c(4,16), 
     xlab="отн. превышения в окрестности 10 м, м", ylab="топоиндекс влажности")
dataEllipse(DSM_data$tpi10, DSM_data$twi, DSM_data$sid, c(0.5, 0.95), add = TRUE, plot.points = TRUE, lwd=1, fill=TRUE, fill.alpha=0.15,
            pch=c(21,22,24,25), col=c("#92DC99", "#D39829", "#B7621C", "#793405"))
legend("topright",bg="white",pch=c(21,22,24,25), legend=c("Лч","Чв","Чт", "Чтк"), inset = .02, bty="n",
       col=c("#92DC99", "#E6BE75", "#B7621C", "#793405"),pt.bg=c("white","white","white","white"))

# Визуализируем карту, построенную методом опорных векторов
cat = levels(as.factor(map.svm.c))[[1]] # список числовых кодов атрибутивной таблицы растра
cat[["sid"]] = c("Лч","Чв", "Чт","Чтк") # их расшифровка (легенда)
levels(map.svm.c) = cat                 # добавить легенду к аттрибутам растра
area_colors = c("#92DC99", "#E6BE75", "#D89156", "#934D08")      # цвета HEX для каждой категории по порядку
levelplot(map.svm.c, col.regions=area_colors, xlab="", ylab="")  # построить изображение

# обучение по части обучающей выборки (70%)
mod.svm.trn = randomForest(sid~tpi10+twi, data=DSM.train, ntree=100,mtry=2) 
# расчет модельных значений
pred.svm.trn = as.data.frame(predict(mod.svm.trn)) # 70%
pred.svm.ver = as.data.frame(predict(mod.svm.trn, DSM.verif)) # 30%
# оценка точности и верификация модели 
goofcat(observed = DSM.train$sid, predicted = pred.svm.trn$predict) # 70% 
goofcat(observed = DSM.verif$sid, predicted = pred.svm.ver$predict) # 30%

##########################################################################################
# Задача 4: сохранить результаты моделирования
##########################################################################################
# варианты сохранения в файл
# 1) Одновременно предсказание и сохранение в оперативную память и в файл (.tiff, .grd, .sdat и др.: см. writeRaster {raster})
map.lda.c = predict(covStack, mod.lda, type = "class", filename = "map_lda_c.tiff", 
                    format = "GTiff", overwrite = TRUE, datatype = "FLT4S") 
# 2) Запись готовой карты в файл
writeRaster(map.lda.c, filename = "map_lda_с.tiff", 
            format = "GTiff", overwrite = TRUE, datatype = "FLT4S")

############################################################################
# Задача 5: Контрольные вопросы и задания для самостоятельной работы
############################################################################
# 1. самостоятельно постройте матрицы соответсвия (матрицы корреляций) между прогнозными 
# значениями четырех моделей по примеру:
goofcat(as.vector(map.lda.c), as.vector(map.rf.c))

# 2. измените априорные вероятности в моделях lda и svm и оцените изменение результата

# 3. Рассчитайте индекс спутанности (CI)
cellStats(map.lda.p1, mean) # Лч
cellStats(map.lda.p2, mean) # Чв
cellStats(map.lda.p3, mean) # Чт
cellStats(map.lda.p4, mean) # Чтк
# CI = 1 - (map.lda.p3 - map.lda.p2)
plot(CI, col = grey(seq(0, 1, len = 25)))

# 4.	Исследуйте влияние на результаты моделирования ПЛС расширения набора ковариат
mod.rf  = randomForest(sid ~ ., data=DSM_data, ntree = 300, mtry = 2) 
importance(mod.rf) # относительная важность каждой ковариаты
varImpPlot(mod.rf)

# 5.	Попробуйте предложить новый метод разделения признакового пространства на подобласти на основе обучающей выборки.
###########################################################################################
# END


#library(nnet) 
#library(C50)
#library(svmpath)
install.packages('C50', dependencies = TRUE)
#install.packages("svmpath", dependencies = TRUE)
install.packages('nnet', dependencies = TRUE)