Clasificacion Supervisada
Babotellos
11 de abril de 2021
Librerias
library(terra)
library(raster)entorno y procesamiento de datos
input <- function(name,num='',type,foldata=F,folder){
ruta <- getwd()
if (foldata == T){
carpeta <- folder
file <- paste0(ruta,'/',carpeta,'/',name, num, type)
}else{
file <- paste0(ruta,'/',name, num, type)
}
}
bruto_landsat <- rast(input('LC08_L1TP_009057_20140820_20200911_02_T1_B',1:7,'.tif',T,'Datos'))
names(bruto_landsat) <- c('ultra-blue','blue','green','red','NIR','SWIR1','SWIR2')Recorte
xmx <- 374838
xmn <- 358682
ymx <- 467578
ymn <- 454847
e <- ext(xmn, xmx, ymn, ymx)
# crop landsat by the extent
landsat <- crop(bruto_landsat, e)
LSbrick <- brick(input('math',type = '.tif'))datos de referencia
valle1 <- shapefile(input('valle1','','.shp',T,'shapes'))
valler <- valle1[98:99,]
valle <- rbind(valle1[10:97,],valle1[100:102,])
plotRGB(LSbrick,r=5,g=5,b=5, axes = T)
plot(valle, add= T)
Clasificacion
sampling sobre los poligonos
library(sp)
set.seed(1)
samples <- spsample(valle, 10000, type='random')
samples$estado <- over(samples,valle)$estadoExtraer valores espectrales
xy <- as.matrix(geom(samples)[,c('x','y')])
df <- extract(landsat, xy)[,-1]
head(df)sampdata <- data.frame(estado = samples$estado, df)entrenar al clasificador
library(rpart)
cartmodel <- rpart(as.factor(estado)~., data = sampdata, method = 'class',
minsplit = 5)
print(cartmodel)## n= 10000
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 10000 6666 c (0.33 0.23 0.15 0.072 0.21)
## 2) SWIR2< 8677.5 4872 1771 c (0.64 0.045 0.15 0.15 0.024)
## 4) SWIR1>=9335.5 4188 1097 c (0.74 0.052 0.17 0.011 0.026)
## 8) green< 9289.5 3187 405 c (0.87 0.055 0.035 0.013 0.024)
## 16) NIR>=17255.5 2954 216 c (0.93 0.01 0.038 0.0064 0.019) *
## 17) NIR< 17255.5 233 89 n (0.19 0.62 0 0.099 0.094) *
## 9) green>=9289.5 1001 387 p (0.31 0.044 0.61 0.004 0.03)
## 18) blue>=9657.5 196 77 c (0.61 0.2 0.02 0.02 0.15) *
## 19) blue< 9657.5 805 195 p (0.24 0.0062 0.76 0 0) *
## 5) SWIR1< 9335.5 684 19 rio (0.015 0 0 0.97 0.013) *
## 3) SWIR2>=8677.5 5128 3034 n (0.045 0.41 0.15 0.0014 0.39)
## 6) NIR< 16367.5 3206 1216 n (0.019 0.62 0.0075 0.0022 0.35)
## 12) SWIR2< 12171.5 2200 413 n (0.028 0.81 0.011 0.0032 0.15) *
## 13) SWIR2>=12171.5 1006 203 urbano (0 0.2 0 0 0.8)
## 26) red>=10999 283 113 n (0 0.6 0 0 0.4)
## 52) SWIR2< 13885 179 20 n (0 0.89 0 0 0.11) *
## 53) SWIR2>=13885 104 11 urbano (0 0.11 0 0 0.89) *
## 27) red< 10999 723 33 urbano (0 0.046 0 0 0.95) *
## 7) NIR>=16367.5 1922 1035 urbano (0.089 0.054 0.39 0 0.46)
## 14) SWIR2< 10213.5 1176 435 p (0.15 0.048 0.63 0 0.18)
## 28) green>=9423 890 172 p (0.12 0.048 0.81 0 0.029) *
## 29) green< 9423 286 105 urbano (0.24 0.045 0.08 0 0.63) *
## 15) SWIR2>=10213.5 746 66 urbano (0 0.064 0.024 0 0.91) *Arbol de clasificacion
plot(cartmodel, uniform=TRUE, main="Classification Tree")
text(cartmodel, cex = 0.6)
clasficar
modelo predictivo con criterio probabilistico
classified <- predict(landsat, cartmodel, na.rm = TRUE)
classified## class : SpatRaster
## dimensions : 424, 538, 5 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 358695, 374835, 454845, 467565 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : c, n, p, rio, urbano
## min values : 0, 0, 0, 0, 0
## max values : 0.9268788, 0.8882682, 0.8067416, 0.9722222, 0.9543568plot(classified)
Diseño combinado
class <- c("cutivo","suelo","pastos","rio","urbano")
mycolor <- c('darkgreen', 'Brown', 'chartreuse', 'blue', 'gray')
classdf <- data.frame(classvalue = c(1,2,3,4,5),
classnames = class,
color = mycolor)
lulc <- app(classified, fun = which.max)
lulcc <- as.factor(lulc)
levels(lulcc) <- as.character(classdf$classnames)
plot(lulcc,col=mycolor,main='Supervised clasification cover')
plot(valle,add=T,border='red',lwd=1)
plot(valler,col='white',border=F,add=T)
evaluacion del modelo
generacion de 5 muestras
set.seed(99)
k <- 5
j <- sample(rep(1:k, each = round(nrow(sampdata))/k))
table(j)## j
## 1 2 3 4 5
## 2000 2000 2000 2000 2000prueba del modelo
x <- list()
for (k in 1:5) {
train <- sampdata[j!= k, ]
test <- sampdata[j == k, ]
cart <- rpart(as.factor(estado)~., data=train, method = 'class',
minsplit = 5)
pclass <- predict(cart, test, na.rm = TRUE)
pclass <- apply(pclass, 1, which.max)
x[[k]] <- cbind(test$estado, as.integer(pclass))
}Convinar los 5 marco de datos
y <- do.call(rbind, x)
y <- data.frame(y)
colnames(y) <- c('observed', 'predicted')
# confusion matrix
conmat <- table(y)
# change the name of the classes
colnames(conmat) <- classdf$classnames
rownames(conmat) <- classdf$classnames
print(conmat)## predicted
## observed cutivo suelo pastos rio urbano
## cutivo 2841 132 310 10 41
## suelo 85 2005 49 3 170
## pastos 111 32 1311 0 54
## rio 22 35 0 661 0
## urbano 103 406 41 8 1570Precision general y "kappa"
n <- sum(conmat)
paste0('n= ',n)## [1] "n= 10000"##Precision general
diag <- diag(conmat)
OA <- sum(diag) / n
paste0('precision general ',OA)## [1] "precision general 0.8388"## Kappa
rowsums <- apply(conmat, 1, sum)
p <- rowsums / n
colsums <- apply(conmat, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
kappa <- (OA - expAccuracy) / (1 - expAccuracy)
paste0('Kappa = ',kappa)## [1] "Kappa = 0.789139998326724"Precision de usuario
PA <- diag / colsums
UA <- diag / rowsums
outAcc <- data.frame(producerAccuracy = PA, userAccuracy = UA)
outAccReplanteamiento del modelo
## remover urbano del modelo
valle2 <- rbind(valle1[10:97,],valle1[100,])
set.seed(1)
samples2 <- spsample(valle2, 10000, type='random')
samples2$estado <- over(samples2,valle2)$estadoxy2 <- as.matrix(geom(samples2)[,c('x','y')])
df2 <- extract(landsat, xy2)[,-1]
head(df2)sampdata2 <- data.frame(estado = samples2$estado, df2)cartmodel2 <- rpart(as.factor(estado)~., data = sampdata2, method = 'class',
minsplit = 4)
print(cartmodel2)## n= 10000
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 10000 5800 c (0.420000000 0.285900000 0.202400000 0.091700000)
## 2) NIR>=17052 5984 2124 c (0.645053476 0.024064171 0.326537433 0.004344920)
## 4) green< 9324.5 3879 332 c (0.914410931 0.014436710 0.064449600 0.006702758) *
## 5) green>=9324.5 2105 401 p (0.148693587 0.041805226 0.809501188 0.000000000) *
## 3) NIR< 17052 4016 1301 n (0.084661355 0.676045817 0.017430279 0.221862550)
## 6) SWIR1>=9872.5 3136 426 n (0.106505102 0.864158163 0.022321429 0.007015306) *
## 7) SWIR1< 9872.5 880 11 rio (0.006818182 0.005681818 0.000000000 0.987500000) *nuevo arbol de clasificacion
plot(cartmodel2, uniform=TRUE, main="Classification Tree")
text(cartmodel2, cex = 0.7)
### nuevo modelo
classified2 <- predict(landsat, cartmodel2, na.rm = TRUE)
classified2## class : SpatRaster
## dimensions : 424, 538, 4 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 358695, 374835, 454845, 467565 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : c, n, p, rio
## min values : 0.006818182, 0.005681818, 0.000000000, 0.000000000
## max values : 0.9144109, 0.8641582, 0.8095012, 0.9875000plot(classified2)
### Mapa
class2 <- c("cutivo","suelo","pastos","rio")
mycolor2 <- c('darkgreen', 'Brown', 'chartreuse', 'blue')
classdf2 <- data.frame(classvalue = c(1,2,3,4),
classnames = class2,
color = mycolor2)
lulc2 <- app(classified2, fun = which.max)
lulcc2 <- as.factor(lulc2)
levels(lulcc2) <- as.character(classdf2$classnames)
plot(lulcc2,col=mycolor2,main='Supervised clasification cover')
plot(valle2,add=T,border='black',lwd=1.2)
plot(valler,col='white',border=F,add=T)
evaluacion del modelo
set.seed(99)
k <- 5
j <- sample(rep(1:k, each = round(nrow(sampdata2))/k))
table(j)## j
## 1 2 3 4 5
## 2000 2000 2000 2000 2000x2 <- list()
for (k in 1:5) {
train <- sampdata2[j!= k, ]
test <- sampdata2[j == k, ]
cart <- rpart(as.factor(estado)~., data=train, method = 'class',
minsplit = 4)
pclass <- predict(cart, test, na.rm = TRUE)
pclass <- apply(pclass, 1, which.max)
x2[[k]] <- cbind(test$estado, as.integer(pclass))
}y2 <- do.call(rbind, x2)
y2 <- data.frame(y2)
colnames(y2) <- c('observed', 'predicted')
# confusion matrix
conmat2 <- table(y2)
print(conmat2)## predicted
## observed 1 2 3 4
## c 3572 308 314 6
## n 76 2683 95 5
## p 251 66 1707 0
## rio 28 22 0 867Precision general y "kappa"
n <- sum(conmat2)
paste0('n= ',n)## [1] "n= 10000"##Precision general
diag <- diag(conmat2)
OA2 <- sum(diag) / n
paste0('precision general ',OA2)## [1] "precision general 0.8829"## Kappa
rowsums <- apply(conmat2, 1, sum)
p <- rowsums / n
colsums <- apply(conmat2, 2, sum)
q <- colsums / n
expAccuracy <- sum(p*q)
kappa2 <- (OA2 - expAccuracy) / (1 - expAccuracy)
paste0('Kappa = ',kappa2)## [1] "Kappa = 0.831791128997401"Precision de usuario
PA <- diag / colsums
UA <- diag / rowsums
outAcc2 <- data.frame(producerAccuracy = PA, userAccuracy = UA)
outAcc2mejora considerablemente la precision general sin embargo sigue siendo bajo
sessionInfo()## R version 4.0.4 (2021-02-15)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.1252 LC_CTYPE=Spanish_Colombia.1252
## [3] LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] rpart_4.1-15 raster_3.3-7 sp_1.4-2 terra_1.0-10
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.5 codetools_0.2-18 lattice_0.20-41 digest_0.6.25
## [5] grid_4.0.4 jsonlite_1.7.2 magrittr_1.5 evaluate_0.14
## [9] highr_0.8 rlang_0.4.10 stringi_1.4.6 vctrs_0.3.6
## [13] rmarkdown_2.7 rgdal_1.5-23 tools_4.0.4 stringr_1.4.0
## [17] xfun_0.21 yaml_2.2.1 compiler_4.0.4 htmltools_0.5.1.1
## [21] knitr_1.31