Clasificación supervisada en R
1 Clasificación Supervisada
Vamos a ver cómo realizar tres tipos de clasificación de imágenes utilizando los siguientes algoritmos:
- Máxima verosimilitud
- Árboles de decisión
- Random Forests
- Evaluación del modelo y precisión
- Estadísticas básicas del resultado
1.1 Tomar áreas de entrenamiento
Para eso vamos a utilizar QGIS
1.2 Extracción de valores de raster por polígonos
Cargar los archivos que vamos a utilizar
## Warning: package 'raster' was built under R version 3.6.3## Loading required package: sp## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3library(curl)
#Ubicación para guardar la imagen descargada
location <- "D:/Descargas/Training.zip"
#Copiar el link para compartir desde google drive
GD_share_URL = "https://drive.google.com/open?id=173AcNuclHuF5Jh39riiNVeg81IbMGtzB"
# Reemplazar "open?" con "us?export=download&"
dwnld_URL <- gsub("open\\?", "uc\\?export=download\\&", GD_share_URL )
dl <- curl::curl_download(dwnld_URL,
destfile = location)
#Extraer el archivo descargado
unzip(location,
exdir = gsub("/Training.zip","",location))
#Substituir .zip por .shp
location <- gsub(".zip",".shp",location)
#Se pueden utilizar cualquiera de las dos funciones, pero optaremos por st_read
training<-shapefile(location)
training<-st_read(location)## Reading layer `Training' from data source `D:\Descargas\Training.shp' using driver `ESRI Shapefile'
## Simple feature collection with 9 features and 3 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -101.2346 ymin: 19.64483 xmax: -101.2202 ymax: 19.65714
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs#Ubicación para guardar la imagen descargada
location <- "D:/Descargas/imagen1.tif"
#Copiar el link para compartir desde google drive
GD_share_URL = "https://drive.google.com/open?id=17PyQnEIICpjNtPP3v59gLEaoY3zzzEcR"
# Reemplazar "open?" con "us?export=download&"
dwnld_URL <- gsub("open\\?", "uc\\?export=download\\&", GD_share_URL )
dl <- curl_download(dwnld_URL,
destfile = location)
S2<-stack(location)
names(S2)<-c("B","G","R","NIR","RE1","RE2","RE3","NNIR","SWIR1","SWIR2")
training## Simple feature collection with 9 features and 3 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -101.2346 ymin: 19.64483 xmax: -101.2202 ymax: 19.65714
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## id ClaseID DescID geometry
## 1 1 1 Vegetacion POLYGON ((-101.2292 19.6526...
## 2 2 1 Vegetacion POLYGON ((-101.229 19.65152...
## 3 3 1 Vegetacion POLYGON ((-101.2257 19.6571...
## 4 4 2 Urbano POLYGON ((-101.2346 19.6456...
## 5 5 2 Urbano POLYGON ((-101.2326 19.6504...
## 6 6 2 Urbano POLYGON ((-101.2332 19.6522...
## 7 7 3 Suelo POLYGON ((-101.2231 19.6455...
## 8 8 3 Suelo POLYGON ((-101.2213 19.6557...
## 9 9 3 Suelo POLYGON ((-101.221 19.65242...## [1] "id" "ClaseID" "DescID" "geometry"## [1] "+proj=longlat +datum=WGS84 +no_defs"
Extraer los valores de las áreas de entrenamiento ¿Qué es esto?
## [[1]]
## B G R NIR RE1 RE2 RE3 NNIR SWIR1 SWIR2
## [1,] 86 236 233 519 965 1393 1475 1530 866 533
## [2,] 67 247 210 459 1167 1399 1448 1495 936 481
## [3,] 94 215 220 482 1154 1529 1606 1682 892 473
## [4,] 76 215 217 540 1400 1363 1449 1616 947 551Convertir la lista de valores a un data.frame
#Se usa la notación :: para llamar funciones dentro de determinados paquetes sin necesidad de cargar el paquete mediante library(paquete)
#Aquí hay que hacer un pequeño truco para poder pasar de list a data.frame
extracted_info<-purrr::map(1:length(extracted_info), function(i) as.data.frame(extracted_info[[i]]))
#Ponerle nombres a cada lista de valores de acuerdo a la ClaseID
names(extracted_info)<-training$DescID
head(extracted_info, n = 1)## $Vegetacion
## B G R NIR RE1 RE2 RE3 NNIR SWIR1 SWIR2
## 1 86 236 233 519 965 1393 1475 1530 866 533
## 2 67 247 210 459 1167 1399 1448 1495 936 481
## 3 94 215 220 482 1154 1529 1606 1682 892 473
## 4 76 215 217 540 1400 1363 1449 1616 947 551#Ya tenemos la información extraida
extracted_info<-dplyr::bind_rows(extracted_info, .id = "id")
head(extracted_info, n = 1)## id B G R NIR RE1 RE2 RE3 NNIR SWIR1 SWIR2
## 1 Vegetacion 86 236 233 519 965 1393 1475 1530 866 533#Cambiar los nombres de las columnas
colnames(extracted_info)<-c("id","B","G","R","NIR","RE1","RE2","RE3","NNIR","SWIR1","SWIR2")
head(extracted_info, n = 1)## id B G R NIR RE1 RE2 RE3 NNIR SWIR1 SWIR2
## 1 Vegetacion 86 236 233 519 965 1393 1475 1530 866 533Veamos como se ven los datos por clase. Vamos a usar ggplot2 para graficar
library(ggplot2)
ggplot(extracted_info, aes(x = id,
y = NIR,
col = id)) +
geom_point(size = 3, alpha = 0.5) +
labs(x = "DescID", col = "DescID") +
theme_bw()
Podríamos ver cómo se ven estos valores en cada una de las bandas
##
## Attaching package: 'tidyr'## The following object is masked from 'package:raster':
##
## extractlibrary(ggplot2)
extracted_info %>%
pivot_longer(cols = -id,
names_to = "bands",
values_to = "reflectance") %>%
ggplot( aes(x = bands,
y = reflectance,
col = id,
group = id)) +
geom_point(size = 3, alpha = 0.5) +
stat_summary(aes(y = reflectance)
, fun.y=mean,
geom="line") +
labs(x = "DescID",
y = "Reflectancia",
col = "DescID") +
theme_bw()
1.3 Clasificación máxima verosimilitud (Maximum likelihood)
library(RStoolbox)
ClassMLC<-superClass(img = S2,
trainData = as(training, 'Spatial'),
responseCol = "DescID",
model = "mlc")## Loading required package: latticelibrary(tmap)
tm_shape(ClassMLC$map)+
tm_raster()+
tm_layout(scale = .8,
legend.position = c("right","bottom"),
legend.frame = T,
legend.bg.color = "white")
Si queremos exportar el resultado para abrirlo en QGIS y compararlo con la imagen original
1.4 Clasificación Árboles de decisión
A partir de la información de entrenamiento se genera un agrupamiento de las categorías.
library(rpart)
library(rpart.plot)
cartmodel <- rpart(as.factor(id)~.,
data = extracted_info,
method = 'class')
print(cartmodel)## n= 102
##
## node), split, n, loss, yval, (yprob)
## * denotes terminal node
##
## 1) root 102 57 Suelo (0.44117647 0.36274510 0.19607843)
## 2) B< 532.5 66 21 Suelo (0.68181818 0.01515152 0.30303030)
## 4) R>=720 45 0 Suelo (1.00000000 0.00000000 0.00000000) *
## 5) R< 720 21 1 Vegetacion (0.00000000 0.04761905 0.95238095) *
## 3) B>=532.5 36 0 Urbano (0.00000000 1.00000000 0.00000000) *
Ahora utilizamos la información obtenida de los datos de entrenamiento para extrapolarlos a toda la imagen y obtener la clasificación.
ClassCART <- predict(S2,
cartmodel,
type = "class")
cols <- c("light blue","orange","dark green")
plot(ClassCART, col = cols)
legend("bottomright",
legend=as.vector(t(ClassCART@data@attributes[[1]][2])),
fill=cols, bg="white")
1.4.1 Clasificación Random Forests
Otra manera de clasificar utilizando random forest. Yo recomiendo esta opción.
Fase de entrenamiento
## randomForest 4.6-14## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:ggplot2':
##
## marginrfmodel <- randomForest(as.factor(id) ~ .,
method = "rf",
data = extracted_info)
#Ver help(train_model_list) para ver todos los modelos disponibles
#Ver help(randomForest) para ver todos los argumentos disponiblesYa obtuvimos el bosque modelos (ya se ajustaron los árboles) Ver el Out of the Bag Error rate de toda la clasificación
##
## Call:
## randomForest(formula = as.factor(id) ~ ., data = extracted_info, method = "rf")
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 3
##
## OOB estimate of error rate: 0.98%
## Confusion matrix:
## Suelo Urbano Vegetacion class.error
## Suelo 45 0 0 0.00000000
## Urbano 0 36 1 0.02702703
## Vegetacion 0 0 20 0.00000000Ver los valores de importancia variable de acuerdo al decremente medio de la precisión. Las variables que tengan un mayor valor serán las que son más importantes para realizar la clasificación.
## MeanDecreaseGini
## B 13.8636821
## G 12.0863197
## R 7.9227143
## NIR 7.5493199
## RE1 1.3290514
## RE2 0.3323111
## RE3 0.3444938
## NNIR 0.8942046
## SWIR1 12.6556069
## SWIR2 7.2739039
Aquí se puede ver el comportamiento del error por clase conforme se aumenta el número de árboles construidos (el tamaño del bosque). Vemos que hay un número en el cual el error se minimiza. Este parámetro se puede indicar a la hora de que construye el modelo (revisar el argumento de “ntree”). Los colores de cada clase están en el mismo orden en el que aparecen en el resumen del modelo; por lo tanto, suelo: se indica en rojo, urbano: azul y vegetación: verde.
1.4.2 Evaluación del modelo
#Ubicación para guardar la imagen descargada
location <- "D:/Descargas/Testing.zip"
#Copiar el link para compartir desde google drive
GD_share_URL = "https://drive.google.com/open?id=1fFX9lepmSWEm8wYDv2CEUnPH4LwMy19H"
# Reemplazar "open?" con "us?export=download&"
dwnld_URL <- gsub("open\\?", "uc\\?export=download\\&", GD_share_URL )
dl <- curl::curl_download(dwnld_URL,
destfile = location)
#Extraer el archivo descargado
unzip(location,
exdir = gsub("/Testing.zip","",location))
#Substituir .zip por .shp
location <- gsub(".zip",".shp",location)
testing<-st_read(location)## Reading layer `Testing' from data source `D:\Descargas\Testing.shp' using driver `ESRI Shapefile'
## Simple feature collection with 6 features and 3 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -101.2352 ymin: 19.64641 xmax: -101.2196 ymax: 19.65617
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defsYa tenemos los datos de validación, así que vamos a correr el modelo para estos datos
extracted_test <- raster::extract(S2,testing)
#Se usa la notación :: para llamar funciones dentro de determinados paquetes sin necesidad de cargar el paquete mediante library(paquete)
#Aquí hay que hacer un pequeño truco para poder pasar de list a data.frame
extracted_test<-purrr::map(1:length(extracted_test), function(i) as.data.frame(extracted_test[[i]]))
#Ponerle nombres a cada lista de valores de acuerdo a la ClaseID
names(extracted_test)<-testing$DescID
head(extracted_test, n = 1)## $Vegetacion
## B G R NIR RE1 RE2 RE3 NNIR SWIR1 SWIR2
## 1 103 297 248 544 1191 1459 1528 1648 1137 667
## 2 156 311 284 535 1146 1336 1331 1384 1029 660
## 3 114 227 235 391 684 784 794 843 832 573
## 4 194 320 351 502 765 910 954 1056 887 637#Ya tenemos la información extraida
extracted_test<-dplyr::bind_rows(extracted_test, .id = "id")
head(extracted_test, n = 1)## id B G R NIR RE1 RE2 RE3 NNIR SWIR1 SWIR2
## 1 Vegetacion 103 297 248 544 1191 1459 1528 1648 1137 667#Cambiar los nombres de las columnas
colnames(extracted_test)<-c("id","B","G","R","NIR","RE1","RE2","RE3","NNIR","SWIR1","SWIR2")
head(extracted_test, n = 1)## id B G R NIR RE1 RE2 RE3 NNIR SWIR1 SWIR2
## 1 Vegetacion 103 297 248 544 1191 1459 1528 1648 1137 667Ya tenemos los datos extraidos del raster para la verificación, así que procedemos a correr el modelo para esos datos. Y vemos su matriz de error
ClassRF <- predict(rfmodel, extracted_test)
confusion_matrix<-table(ClassRF, extracted_test$id)
confusion_matrix##
## ClassRF Suelo Urbano Vegetacion
## Suelo 19 0 0
## Urbano 10 20 0
## Vegetacion 0 0 111.4.3 Clasificación de la imagen
Aquí ya obtenemos la imagen clasificada
## ID value
## 1 1 Suelo
## 2 2 Urbano
## 3 3 Vegetacion1.4.4 Calculo de precisión total y kappa
Ahora ya que tenemos todo deberíamos conocer cuál es la precisión del método
#Obtenemos el número de observaciones en los datos de validación
n_test <- nrow(extracted_test)
#Obtenemos el número de clases en los datos de validación
n_clases <- length(levels(extracted_test$id))
overall_accuracy <- sum(diag(confusion_matrix)) / n_test
overall_accuracy## [1] 0.8333333#Realizar sumas por renglones en la matriz de error
#Obtener la precisión del usuario
users_accuracy <- sapply(1:nrow(confusion_matrix), function(i){
confusion_matrix[i,i] / sum( confusion_matrix[i,])
})
names(users_accuracy) <- colnames(confusion_matrix)
users_accuracy## Suelo Urbano Vegetacion
## 1.0000000 0.6666667 1.0000000#Obtener la precisión del producto
producers_accuracy <- sapply(1:nrow(confusion_matrix), function(i){
confusion_matrix[i,i] / sum( confusion_matrix[,i])
})
names(producers_accuracy) <- row.names(confusion_matrix)
producers_accuracy## Suelo Urbano Vegetacion
## 0.6551724 1.0000000 1.0000000Una forma que quizás sea más sencilla
## Confusion Matrix and Statistics
##
## Reference
## Prediction Suelo Urbano Vegetacion
## Suelo 19 0 0
## Urbano 10 20 0
## Vegetacion 0 0 11
##
## Overall Statistics
##
## Accuracy : 0.8333
## 95% CI : (0.7148, 0.9171)
## No Information Rate : 0.4833
## P-Value [Acc > NIR] : 2.029e-08
##
## Kappa : 0.7423
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: Suelo Class: Urbano Class: Vegetacion
## Sensitivity 0.6552 1.0000 1.0000
## Specificity 1.0000 0.7500 1.0000
## Pos Pred Value 1.0000 0.6667 1.0000
## Neg Pred Value 0.7561 1.0000 1.0000
## Prevalence 0.4833 0.3333 0.1833
## Detection Rate 0.3167 0.3333 0.1833
## Detection Prevalence 0.3167 0.5000 0.1833
## Balanced Accuracy 0.8276 0.8750 1.00001.4.5 Estadísticas básicas de resultado
Se pueden calcular estadísticas de las imágenes clasificadas.
## value count
## [1,] 1 1845
## [2,] 2 4529
## [3,] 3 1969
Si quisieramos conocer el área ocupada por una categoría
#Area en km2
ar<-area(ClassRF_raster)
#Area promedio por pixel en km2
area_pix<-mean(ar@data@values)
area_pix## [1] 0.0003743235## [1] 0.01934744## value count
## [1,] 1 1845
## [2,] 2 4529
## [3,] 3 1969## num [1:3, 1:2] 1 2 3 1845 4529 ...
## - attr(*, "dimnames")=List of 2
## ..$ : NULL
## ..$ : chr [1:2] "value" "count"#Convertir la lista a data.frame
freq_table<-data.frame(matrix(unlist(resul_areas),
nrow=3,
byrow=F))
colnames(freq_table)<-c("id","num_pix")
freq_table$area_ha <- freq_table$num_pix * (area_pix * 100)
freq_table## id num_pix area_ha
## 1 1 1845 69.06269
## 2 2 4529 169.53113
## 3 3 1969 73.70430</div>
</div>