Clasificación de coberturas aplicando Machine Learning
Universidad Distrital Francisco José de Caldas, Facultad de ingenierÃa, Programa de IngenierÃa Catastral y Geodesia
Jhonathan Aponte Saravia
2023-12-13
Introducción
En este ejercicio se intentó mostrar los proceso de clasificación de las oberturas, aplicando fundamentos conceptuales y algoritmos de machine learning, para ello, se utilizaron imágenes satelitales del programa CABERS 4A, que se encuentra disponible en la siguiente: Hacer clic aquÃ, ejecutándose los procesos de conversión de nivel digital a valores de reflectancia al tope de la atmosfera; Luego, se ejecutaron los procesos de aplicación de maching learning, a partir un conjunto de valores de pÃxeles que corresponden a coberturas conocidas, con ello se generó los proceso de entrenamiento para generar el modelo y con el otro grupo de pÃxeles se realizó el proceso de test de evaluación para conocer la confiabilidad, aplicando la metodologÃa de validación cruzada, luego se aplicó el proceso de predicción espacial en todo el conjunto de los datos para generar el proceso de agrupación de pÃxeles, logrando agrupar pÃxeles que corresponde a los tipos de cobertura del paisaje.Â
Por otro lado, fue comprobar la efectividad de la guÃa publicada por Kwesie Benjamin, (2021); y además, se esclarece la funcionalidad de los códigos, para que otras personas del idioma hispano puedan replicarlos con mayor facilidad, utilizando el programa R Package y RStudio.
LibrerÃas a cargar para el ejercicio.
library(sf)## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUElibrary(raster)## Loading required package: splibrary(ggplot2)
library(rgdal)## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-6, (SVN revision 1201)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.1, released 2021/12/27
## Path to GDAL shared files: C:/Users/57311/Documents/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/57311/Documents/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.library(RStoolbox)setwd("D:/.....")
files <- list.files(pattern='.sdat')
files## character(0)b1 <- raster("D:/A_manuscrito/Para ponencia/Datos/Band 1.sdat")
b2 <- raster("D:/A_manuscrito/Para ponencia/Datos/Band 2.sdat")
b3 <- raster("D:/A_manuscrito/Para ponencia/Datos/Band 3.sdat")
b4 <- raster("D:/A_manuscrito/Para ponencia/Datos/Band 4.sdat")Fase de preprocesamiento:
CaracterÃsticas de los datos de las bandas multiespectrales, con de tamaño de pixel 8 metros.
b1## class : RasterLayer
## dimensions : 1355, 1642, 2224910 (nrow, ncol, ncell)
## resolution : 8, 8 (x, y)
## extent : 338258, 351394, 8745258, 8756098 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs
## source : Band 1.sdat
## names : Band_1
## values : 114, 996 (min, max)El proceso de conversión de nivel digital a radiancia, nos valemos en el siguiente manuscrito: Hacer clic aquÃ
rad <- function(x1,x2,x3,x4){
rb1 <- 1.68*x1+(0)
rb2 <- 1.62*x2+(0)
rb3 <- 1.59*x3+(0)
rb4 <- 1.24*x4+(0)
return(list(rb1,rb2,rb3,rb4))
}
rabs<- stack(rad(b1,b2,b3,b4))Cálculo de distancia sol a la tierra, en unidades astronómicas.
dist<- function(j=228){
da <- (1-0.01672*cos(0.9856*(j-4)*(pi/180)))
return(da)}# Ingresamos el valor del dÃa juliano, en este caso es dÃa 228.
d<-dist()
d## [1] 1.012662Cálculo de ángulo cenital.
ac <- function(t=58.4011){
a <- cos((pi/180)*(90-t))
return(a)
}
z<- ac()
z## [1] 0.851737Cálculo del valor de reflectancia
Para determinar los valores de Esun, consultamos la siguiente url: Hacer clic aquÃ
ref <- function(){
rf1<- ((pi)*rabs$Band_1*d^2)/(1958*z)
rf2<- ((pi)*rabs$Band_2*d^2)/(1852*z)
rf3<- ((pi)*rabs$Band_3*d^2)/(1559*z)
rf4<- ((pi)*rabs$Band_4*d^2)/(1091*z)
return(list(rf1,rf2,rf3,rf4))
}
rfbs<-stack(ref())
rfbs## class : RasterStack
## dimensions : 1355, 1642, 2224910, 4 (nrow, ncol, ncell, nlayers)
## resolution : 8, 8 (x, y)
## extent : 338258, 351394, 8745258, 8756098 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs
## names : Band_1, Band_2, Band_3, Band_4
## min values : 0.2791054, 0.2051346, 0.2044561, 0.2106523
## max values : 3.320057, 3.384721, 3.946389, 4.397905Generando el proceso de rescalamiento.
res <- function(){
res1<-((rfbs$Band_1-0.2791054)/(3.320057-0.2791054))
res2<-((rfbs$Band_2-0.2051346)/(3.384721-0.2051346))
res3<-((rfbs$Band_3-0.2044561)/(3.946389-0.2044561))
res4<-((rfbs$Band_4-0.2106523)/(4.397905-0.2106523))
return(list(res1,res2,res3,res4))
}
resc<-stack(res())
resc## class : RasterStack
## dimensions : 1355, 1642, 2224910, 4 (nrow, ncol, ncell, nlayers)
## resolution : 8, 8 (x, y)
## extent : 338258, 351394, 8745258, 8756098 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs
## names : Band_1, Band_2, Band_3, Band_4
## min values : 1.185213e-08, 7.796832e-09, 7.050049e-09, 6.829703e-09
## max values : 0.9999998, 1.0000001, 1.0000000, 0.9999999summary(resc$Band_1)## Band_1
## Min. 1.185213e-08
## 1st Qu. 1.056563e-01
## Median 1.259338e-01
## 3rd Qu. 1.579509e-01
## Max. 9.999998e-01
## NA's 1.961200e+04hist(resc$Band_4)
plot(resc$Band_2)
Explorando los datos preprocesados
Activamos las siguientes librerÃas previamente:
library(ggplot2)
library(RStoolbox)Aplicamos el siguiente código para mostrar la combinación de bandas.
ggRGB(resc,r=3,g=2,b=1,stretch = "lin") + ggtitle("Combinación bandas 321") 
Procesamiento de datos aplicando algorÃtmos de maching Learning
A continuación, los códigos para la clasificación de coberturas, ha sido adaptado a partir de la guÃa, publicado por Kwesie Benjamin, (2021); para ver la guÃa Hacer clic aquÃ:.
Activamos las siguientes librerÃas:
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:raster':
##
## intersect, select, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(raster)
library(tmap)## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')library(caret)## Loading required package: latticeSe importan las ubicaciones de los objetos conocido, en formato vectorial de tipo polÃgono.
files <- list.files(pattern='.shp')
files## character(0)roiTraining <- sf::read_sf("D:/A_manuscrito/Para ponencia/Datos/Muestra.shp")
roiTraining## Simple feature collection with 5 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 338808 ymin: 8748904 xmax: 348323 ymax: 8753329
## Projected CRS: WGS 84 / UTM zone 18S
## # A tibble: 5 x 3
## id Cobertura geometry
## <dbl> <chr> <POLYGON [m]>
## 1 1 Glaciares ((339097.8 8750360, 338936 8750380, 338909.1 8750414,~
## 2 2 Lagunas ((342094.2 8750321, 342151.5 8750441, 342259.4 875041~
## 3 3 Humedales ((348201.6 8753272, 348299.4 8753329, 348323 8753238,~
## 4 4 Otros ((347790.4 8751947, 347891.5 8751995, 347942.1 875195~
## 5 5 Afloramiento roca ((339728.1 8749268, 339997.8 8749146, 340092.1 874903~Identificamos los contenidos de la tabla de atributos de dato vectorial.
roiTraining$Cobertura <- as.factor(roiTraining$Cobertura)Seleccionamos los atributos de los tipo de cobertura, zona de entrenamiento.
unique(roiTraining$Cobertura)## [1] Glaciares Lagunas Humedales Otros
## [5] Afloramiento roca
## Levels: Afloramiento roca Glaciares Humedales Lagunas OtrosPreviamente, los datos espaciales del dato raster y del dato vectorial debe estar en el mismo sistema de referencia espacial.
crs(roiTraining)## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["WGS 84 / UTM zone 18S",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 18S",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-75,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",10000000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## ID["EPSG",32718]]crs(resc)## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=18 +south +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["WGS 84 / UTM zone 18S",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 18S",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-75,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",10000000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["easting",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## ID["EPSG",32718]]Se extrae los valores de los pixeles de los datos raster que corresponde a la zona de entrenamiento.
extract <- raster::extract(resc,roiTraining,df=TRUE)
head(extract)## ID Band_1 Band_2 Band_3 Band_4
## 1 1 0.7886872 0.6732571 0.7515464 0.6889117
## 2 1 0.8751333 0.7918835 0.8948454 0.8193018
## 3 1 0.9519742 0.9146723 1.0000000 0.8993839
## 4 1 0.9199572 0.8418315 0.8958763 0.7515400
## 5 1 0.8644609 0.7679501 0.8969072 0.6878850
## 6 1 0.7694769 0.6441208 0.7639175 0.6334702Se genera el jion entre los coordenads de los puntos de entrenamiento con las coordenadas de los pixeles del dato raster, con (ID y id), para obtener una data frame con valores conjuntas.
extractMerged <- dplyr::inner_join(extract,roiTraining,by= c("ID"="id"))Suprimimos los valores que tienen datos nulos.
extractMerged$geometry <- NULLDividimos el data data frame, en proprorciones de 70% y 30% para aplicar los procesos de entrenamiento.
set.seed(123)
trainIndex <- caret::createDataPartition(extractMerged$ID,list = FALSE,p=0.7)
trainData <- extractMerged[trainIndex,] # 70% for training Data
testData <- extractMerged[-trainIndex,] # 30% for testing DataVisualizamos los datos de tipos de cobertura en la dataframe de entrenamiento.
classCount <- trainData %>%
dplyr::group_by(Cobertura) %>%
count()
print(classCount)## # A tibble: 5 x 2
## # Groups: Cobertura [5]
## Cobertura n
## <fct> <int>
## 1 Afloramiento roca 1390
## 2 Glaciares 1029
## 3 Humedales 639
## 4 Lagunas 699
## 5 Otros 349Mostrando la proporción de los pixeles de entremaniento.
(prop.table(table(trainData$Cobertura))*100)##
## Afloramiento roca Glaciares Humedales Lagunas
## 33.852898 25.060887 15.562591 17.023868
## Otros
## 8.499756Presentamos el contenido de la data frame del entrenamiento.
print(head(trainData))## ID Band_1 Band_2 Band_3 Band_4 Cobertura
## 1 1 0.7886872 0.6732571 0.7515464 0.6889117 Glaciares
## 2 1 0.8751333 0.7918835 0.8948454 0.8193018 Glaciares
## 4 1 0.9199572 0.8418315 0.8958763 0.7515400 Glaciares
## 5 1 0.8644609 0.7679501 0.8969072 0.6878850 Glaciares
## 7 1 0.7299892 0.6732571 0.7546392 0.5831622 Glaciares
## 8 1 0.7054428 0.6056192 0.7144330 0.6088295 GlaciaresSe prepara los datos para construir el modelo de entrenamiento de Random forest.
Variable de los objetos de la muestra.
respVar <- c("Cobertura")Predicción de variables con los datos de los pixeles del datos raster.
predVar <- c("Band_1","Band_2","Band_3","Band_4")Luego de construir el modelo utilizamos los datos de prueba.
set.seed(124)
cvControl <- caret::trainControl(method = 'cv',
number = 10,
savePredictions = TRUE,
verboseIter = FALSE)Aplicamos de algorÃtmo de random forest y avaluamos la exactitud de coeficiente de Kappa.
library(randomForest)## randomForest 4.7-1.1## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:dplyr':
##
## combine## The following object is masked from 'package:ggplot2':
##
## marginset.seed(124)
rfModel <- caret::train(trainData[,predVar],
trainData[,respVar],
method="rf",
metric = "Kappa",
ntree= 500,
trControl= cvControl,
tuneLength=6,
importance=TRUE)## note: only 3 unique complexity parameters in default grid. Truncating the grid to 3 .Mostramos el modelo de predicción, aplicando el método de validación cruzada rfModel.
print(rfModel)## Random Forest
##
## 4106 samples
## 4 predictor
## 5 classes: 'Afloramiento roca', 'Glaciares', 'Humedales', 'Lagunas', 'Otros'
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 3696, 3695, 3695, 3695, 3697, 3695, ...
## Resampling results across tuning parameters:
##
## mtry Accuracy Kappa
## 2 0.9929393 0.9907415
## 3 0.9934259 0.9913794
## 4 0.9931814 0.9910595
##
## Kappa was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 3.Evaluamos el desempeño del modelo utilizando los datos de prueba.
rfPredict <- predict(rfModel,testData)Utilizamos la matriz de confusión, teniendo en cuenta la matriz de Kappa.
confusionMatrix(rfPredict,testData$Cobertura)## Confusion Matrix and Statistics
##
## Reference
## Prediction Afloramiento roca Glaciares Humedales Lagunas Otros
## Afloramiento roca 581 0 0 0 0
## Glaciares 0 407 0 0 0
## Humedales 0 0 268 0 6
## Lagunas 1 0 0 333 0
## Otros 0 0 5 0 156
##
## Overall Statistics
##
## Accuracy : 0.9932
## 95% CI : (0.9881, 0.9965)
## No Information Rate : 0.3312
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9911
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: Afloramiento roca Class: Glaciares Class: Humedales
## Sensitivity 0.9983 1.0000 0.9817
## Specificity 1.0000 1.0000 0.9960
## Pos Pred Value 1.0000 1.0000 0.9781
## Neg Pred Value 0.9991 1.0000 0.9966
## Prevalence 0.3312 0.2316 0.1554
## Detection Rate 0.3307 0.2316 0.1525
## Detection Prevalence 0.3307 0.2316 0.1559
## Balanced Accuracy 0.9991 1.0000 0.9888
## Class: Lagunas Class: Otros
## Sensitivity 1.0000 0.96296
## Specificity 0.9993 0.99687
## Pos Pred Value 0.9970 0.96894
## Neg Pred Value 1.0000 0.99624
## Prevalence 0.1895 0.09220
## Detection Rate 0.1895 0.08879
## Detection Prevalence 0.1901 0.09163
## Balanced Accuracy 0.9996 0.97991Selección de caracterÃsticas (eliminación de caracterÃsticas recursivas).
Intenta eliminar las caracterÃsticas no predictoras para mejorar el rendimiento.
set.seed(124)
indexRfe <- createMultiFolds(trainData$Cobertura, times=5)Se indica una secuencia de cómo queremos que se subdividan las variables predictoras.
predSeq <-seq(from=1, to =length(predVar),by=2)Se genera la función de la eliminación de la caracterÃsticas recursivas(ECR).
rfeCont <- rfeControl(method="cv",
number = 10,
verbose=FALSE,
functions=rfFuncs,
index=indexRfe)Eliminación de las caracterÃsticas recursivas (ntree=500)
Previamente se procesamos nuestros datos de entrenamiento antes de introducirlos en el algorÃtmo ECR.
set.seed(124)
rfModelRfe <- caret::rfe(trainData[,predVar],
trainData[,respVar],
sizes = predSeq,
metric = "Kappa",
ntree=500,
method="rf",
rfeControl = rfeCont)Mostramos el resultados del modelo, para verificar las variables con las respectivas exactitudes y valores de coeficiente de Kappa.
print(rfModelRfe)##
## Recursive feature selection
##
## Outer resampling method: Cross-Validated (10 fold)
##
## Resampling performance over subset size:
##
## Variables Accuracy Kappa AccuracySD KappaSD Selected
## 1 0.9637 0.9523 0.020537 0.027229
## 3 0.9930 0.9908 0.003571 0.004682
## 4 0.9933 0.9912 0.003733 0.004895 *
##
## The top 4 variables (out of 4):
## Band_1, Band_4, Band_3, Band_2En seguida se evalúa con los datos que corresponde al test de datos aplicando el modelo de ECR.
rfElemPredict <- predict(rfModelRfe,testData)
confusionMatrix(rfElemPredict$pred,testData$Cobertura)## Confusion Matrix and Statistics
##
## Reference
## Prediction Afloramiento roca Glaciares Humedales Lagunas Otros
## Afloramiento roca 580 0 0 0 0
## Glaciares 0 407 0 0 0
## Humedales 0 0 268 0 5
## Lagunas 2 0 0 333 0
## Otros 0 0 5 0 157
##
## Overall Statistics
##
## Accuracy : 0.9932
## 95% CI : (0.9881, 0.9965)
## No Information Rate : 0.3312
## P-Value [Acc > NIR] : < 2.2e-16
##
## Kappa : 0.9911
##
## Mcnemar's Test P-Value : NA
##
## Statistics by Class:
##
## Class: Afloramiento roca Class: Glaciares Class: Humedales
## Sensitivity 0.9966 1.0000 0.9817
## Specificity 1.0000 1.0000 0.9966
## Pos Pred Value 1.0000 1.0000 0.9817
## Neg Pred Value 0.9983 1.0000 0.9966
## Prevalence 0.3312 0.2316 0.1554
## Detection Rate 0.3301 0.2316 0.1525
## Detection Prevalence 0.3301 0.2316 0.1554
## Balanced Accuracy 0.9983 1.0000 0.9892
## Class: Lagunas Class: Otros
## Sensitivity 1.0000 0.96914
## Specificity 0.9986 0.99687
## Pos Pred Value 0.9940 0.96914
## Neg Pred Value 1.0000 0.99687
## Prevalence 0.1895 0.09220
## Detection Rate 0.1895 0.08936
## Detection Prevalence 0.1907 0.09220
## Balanced Accuracy 0.9993 0.98300Aplicación de predicción espacial
Aplicamos la predicción espacial utilizando el modelo con los datos de las imágenes.
bawkuPredictions<- raster::predict(resc,rfModel)Mostramos los tipos de objetos
plot(bawkuPredictions)
Luego elegimos patrones para asociar los tipos de cobertura.
pal <-c('#ffdab9','#fffafa','#6b8e23','#87ceff','#b3ee3a')A continuación presentamos el mapa de coberturas.
tm_shape(bawkuPredictions)+
tm_raster(style = "cat",labels = c("Afloramiento roca",
"Glaciares","Humedales",
"Lagunas","Otros tipos de cobertura"
),
palette = pal,
title = "Leyenda")+
tm_compass(position = c("left", "top"))+
tm_scale_bar(position = c("right", "bottom"))+
tm_layout(main.title= "Tipos de cobertura imágenes CBERS 4A",
main.title.size =.9 )+
tm_layout(legend.outside = TRUE,
legend.outside.position = "right",
legend.title.size = 1) ## stars object downsampled to 1101 by 908 cells. See tm_shape manual (argument raster.downsample)