1 Preparar sesión de trabajo.

1.1 Instalar paquetes.

Instalamos y cargamos paquetes necesarios siempre que no se encuentren instalados previamente

if(!require(rgdal)) suppressMessages( suppressWarnings( install.packages("rgdal") ))
if(!require(raster)) suppressMessages( suppressWarnings( install.packages("raster") ))
if(!require(randomForest)) suppressMessages( suppressWarnings( install.packages("randomForest") ))
if(!require(caret)) suppressMessages( suppressWarnings( install.packages("caret") ))
if(!require(foreach)) suppressMessages( suppressWarnings( install.packages("foreach") ))
if(!require(snow)) suppressMessages( suppressWarnings( install.packages("snow") ))
if(!require(import)) suppressMessages( suppressWarnings( install.packages("import") ))
if(!require(doParallel)) suppressMessages( suppressWarnings( install.packages("doParallel") ))

1.2 Definir directorios.

Podemos seleccionar el directorio de trabaja manualmente o usando la interfaz gráfica

### setwd("Z:/OneDrive - Universidad Mayor/documentos/application_work/mayor/cursos/Magister/Teledeteccion/Clase_04_Random_Forest")
ruta <- getwd()
ruta_imagenes <- file.path(ruta,"imagenes")
getwd()

## [1] "/home/matbox/Documents/TrabajosExtra/nelson.mattie/Taller_GEE_R/Clase_03_Modelamiento"

2 Manipulación de datos.

2.1 Cargar rasters.

Creamos una lista de las imágenes que están en el directorio imágenes

grids <- list.files(file.path(ruta_imagenes), pattern = "B\\d.tif$|ndvi.tif$", all.files=TRUE, full.name=TRUE)

Revisamos la lista

print(grids)

## [1] "/home/matbox/Documents/TrabajosExtra/nelson.mattie/Taller_GEE_R/Clase_03_Modelamiento/imagenes/B2.tif"   "/home/matbox/Documents/TrabajosExtra/nelson.mattie/Taller_GEE_R/Clase_03_Modelamiento/imagenes/B3.tif"  
## [3] "/home/matbox/Documents/TrabajosExtra/nelson.mattie/Taller_GEE_R/Clase_03_Modelamiento/imagenes/B4.tif"   "/home/matbox/Documents/TrabajosExtra/nelson.mattie/Taller_GEE_R/Clase_03_Modelamiento/imagenes/B5.tif"  
## [5] "/home/matbox/Documents/TrabajosExtra/nelson.mattie/Taller_GEE_R/Clase_03_Modelamiento/imagenes/B6.tif"   "/home/matbox/Documents/TrabajosExtra/nelson.mattie/Taller_GEE_R/Clase_03_Modelamiento/imagenes/ndvi.tif"

Construimos un stack de los rasters en la carpeta

stack_img <- stack(grids)

Comprobamos que es correcto

print(stack_img)

## class      : RasterStack 
## dimensions : 941, 1121, 1054861, 6  (nrow, ncol, ncell, nlayers)
## resolution : 0.0002694946, 0.0002694946  (x, y)
## extent     : -71.59339, -71.29129, -33.35022, -33.09663  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## names      :        B2,        B3,        B4,        B5,        B6,       ndvi 
## min values :         ?,         ?,         ?,         ?,         ?, -0.3210981 
## max values :         ?,         ?,         ?,         ?,         ?,  0.8064534

Ahora construimos un Brick. Esto es similar a un stack pero se carga directo en la memoria y es mucho mas rápido.

img <- brick(stack_img)

Comprobamos que todo está bien.

print(img)

## class      : RasterBrick 
## dimensions : 941, 1121, 1054861, 6  (nrow, ncol, ncell, nlayers)
## resolution : 0.0002694946, 0.0002694946  (x, y)
## extent     : -71.59339, -71.29129, -33.35022, -33.09663  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /tmp/Rtmp8QHgVd/raster/r_tmp_2021-02-18_205458_15779_84790.grd 
## names      :          B2,          B3,          B4,          B5,          B6,        ndvi 
## min values :  0.08366728,  0.05385134,  0.03522509,  0.01876169,  0.00793536, -0.32109806 
## max values :   0.4530444,   0.4990548,   0.5368390,   0.5475211,   0.6128311,   0.8064534

Ahora podemos graficar y analizar Ejemplo de una imagen color real

RGB<-plotRGB(img, r = 3, g = 2, b = 1,stretch="hist")

Todos los pasos anteriores nos permitieron preparar los datos

3 Random Forest

3.1 Cargamos las áreas que nos servirán para entrenar el modelo

Categorías:

Abierto=1
Agricola=2
Agua=3
Construido=4
Nativo=5
Plantacion=6

trainData <- shapefile("poligonos_clasificacion_numerico.shp", verbose=TRUE)

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/matbox/Documents/TrabajosExtra/nelson.mattie/Taller_GEE_R/Clase_03_Modelamiento", layer: "poligonos_clasificacion_numerico"
## with 60 features
## It has 1 fields
## Integer64 fields read as strings:  CLASS

## Warning in rgdal::readOGR(dirname(x), fn, stringsAsFactors = stringsAsFactors, : Z-dimension discarded

Graficamos para corroborar que todo esta correcto

plot(trainData)

Revisamos además los metadatos

head(trainData)

##   CLASS
## 0     1
## 1     1
## 2     1
## 3     1
## 4     1
## 5     1

Definimos una variable con el nombre de la columna que contiene nuestras categorías

responseCol <- "CLASS"

Ahora debemos extraer todos los valores de pixel correspondiente a cada categoría para cada polígono

dfAll = data.frame(matrix(vector(), nrow = 0, ncol = length(names(img)) + 1))
for (i in 1:length(unique(trainData[[responseCol]]))){
  category <- unique(trainData[[responseCol]])[i]
  categorymap <- trainData[trainData[[responseCol]] == category,]
  dataSet <- extract(img, categorymap)
  if(is(trainData, "SpatialPointsDataFrame")){
    dataSet <- cbind(dataSet, CLASS = as.numeric(rep(category, nrow(dataSet))))
    dfAll <- rbind(dfAll, dataSet[complete.cases(dataSet),])
  }
  if(is(trainData, "SpatialPolygonsDataFrame")){
    dataSet <- dataSet[!unlist(lapply(dataSet, is.null))]
    dataSet <- lapply(dataSet, function(x){cbind(x, CLASS = as.numeric(rep(category, nrow(x))))})
    df <- do.call("rbind", dataSet)
    dfAll <- rbind(dfAll, df)
  }
}

Las primeras observaciones son.

head(dfAll,10)

##           B2        B3        B4        B5        B6      ndvi CLASS
## 1  0.1312589 0.1231767 0.1450761 0.2544180 0.3892568 0.2737010     1
## 2  0.1315018 0.1239561 0.1474614 0.2550045 0.3924881 0.2672103     1
## 3  0.1337016 0.1263023 0.1477698 0.2584519 0.3896789 0.2724671     1
## 4  0.1315292 0.1226795 0.1463656 0.2581257 0.3910405 0.2762980     1
## 5  0.1291223 0.1221435 0.1473554 0.2555434 0.3956028 0.2685239     1
## 6  0.1291223 0.1221435 0.1473554 0.2555434 0.3956028 0.2685239     1
## 7  0.1285011 0.1213913 0.1453276 0.2526425 0.3912966 0.2696556     1
## 8  0.1279366 0.1197529 0.1424878 0.2490218 0.3905855 0.2721110     1
## 9  0.1277785 0.1202034 0.1411352 0.2473336 0.3908309 0.2733768     1
## 10 0.1278626 0.1194967 0.1418382 0.2510220 0.3860075 0.2779202     1

Ahora seleccionamos un número de puntos aleatorios de la lista creada en el paso anterior Generalmente se generan bases de datos con miles de puntos por lo que debemos reducir el número Definimos cuantos puntos queremos usar. Ojo que si tenemos capacidad de proceso ilimitada podemos usarlos todos

nsamples <- 2000

Creamos una nueva lista de n puntos al azar

sdfAll <- dfAll[sample(1:nrow(dfAll), nsamples),]

3.2 Creación modelo

Usamos esa lista para entrenar el modelo

registerDoParallel(cores=4)
ptm <- proc.time()
modFit_rf_1 <- train(as.factor(CLASS) ~  B2 + B3 + B4 + B5 + ndvi , method = "parRF"
                     , data = sdfAll,ntree=50)
proc.time() - ptm

##    user  system elapsed 
##  22.892  15.388  30.051

Evaluamos el modelo

modFit_rf_1

## Parallel Random Forest 
## 
## 2000 samples
##    5 predictor
##    6 classes: '1', '2', '3', '4', '5', '6' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 2000, 2000, 2000, 2000, 2000, 2000, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##   2     0.8578323  0.7816358
##   3     0.8612052  0.7870863
##   5     0.8586398  0.7831923
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 3.

Variables mas importantes

imp_1<-varImp(modFit_rf_1)
plot(imp_1)

Ajusto mi modelo inicial Las variables mas importantes son B5, NDVI y B3

ptm <- proc.time()
modFit_rf_2 <- train(as.factor(CLASS) ~  B3 + B5 + ndvi , method = "parRF"
                     , data = sdfAll,ntree=50)

## note: only 2 unique complexity parameters in default grid. Truncating the grid to 2 .

proc.time() - ptm

##    user  system elapsed 
##  13.748   9.592   6.976

Evaluamos el modelo

modFit_rf_2

## Parallel Random Forest 
## 
## 2000 samples
##    3 predictor
##    6 classes: '1', '2', '3', '4', '5', '6' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 2000, 2000, 2000, 2000, 2000, 2000, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##   2     0.8175654  0.7213768
##   3     0.8165600  0.7197228
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.

Variables mas importantes

imp_2<-varImp(modFit_rf_2)
plot(imp_2, main="B3 + B5 + NDVI")

Que pasa si creamos un modelo sin ndvi

ptm <- proc.time()
modFit_rf_3 <- train(as.factor(CLASS) ~  B2 + B3 + B4 + B5, method = "parRF"
                     , data = sdfAll,ntree=50)
proc.time() - ptm

##    user  system elapsed 
##  21.416  13.784   9.924

Evaluamos el modelo

modFit_rf_3

## Parallel Random Forest 
## 
## 2000 samples
##    4 predictor
##    6 classes: '1', '2', '3', '4', '5', '6' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 2000, 2000, 2000, 2000, 2000, 2000, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##   2     0.8634085  0.7910895
##   3     0.8613495  0.7880867
##   4     0.8554984  0.7791201
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.

Variables mas importantes

imp_3<-varImp(modFit_rf_3)
plot(imp_3)

Modelo Sin capas que componen ndvi

ptm <- proc.time()
modFit_rf_4 <- train(as.factor(CLASS) ~  B2 + B3 + ndvi, method = "parRF"
                     , data = sdfAll,ntree=50)

## note: only 2 unique complexity parameters in default grid. Truncating the grid to 2 .

proc.time() - ptm

##    user  system elapsed 
##  13.920   9.608   7.523

Evaluamos el modelo

modFit_rf_4

## Parallel Random Forest 
## 
## 2000 samples
##    3 predictor
##    6 classes: '1', '2', '3', '4', '5', '6' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 2000, 2000, 2000, 2000, 2000, 2000, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##   2     0.8469978  0.7665002
##   3     0.8439184  0.7620161
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.

Variables mas importantes

imp_4<-varImp(modFit_rf_4)
plot(imp_4)

Caso especial

ptm <- proc.time()
modFit_rf_5 <- train(as.factor(CLASS) ~ B3+ B5 + B6 + ndvi, method = "parRF"
                     , data = sdfAll,ntree=50)
proc.time() - ptm

##    user  system elapsed 
##  21.304  14.104  10.612

Evaluamos el modelo

modFit_rf_5

## Parallel Random Forest 
## 
## 2000 samples
##    4 predictor
##    6 classes: '1', '2', '3', '4', '5', '6' 
## 
## No pre-processing
## Resampling: Bootstrapped (25 reps) 
## Summary of sample sizes: 2000, 2000, 2000, 2000, 2000, 2000, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##   2     0.8350094  0.7466027
##   3     0.8349808  0.7466777
##   4     0.8296604  0.7385960
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 2.

Variables mas importantes

imp_5<-varImp(modFit_rf_5)
plot(imp_5)

3.3 Predicción

Una vez generado el modelo podemos generar nuestra clasificación Este proceso puede ser bastante lento por lo que usaremos un cluster para aprovechar todos los núcleos de proceso que cuenta nuestro computador para acelerar el proceso

Definimos paleta de colores para los gráficos.

cpal <- c('white','darkorange','deepskyblue','bisque4','lightgreen','darkgreen')

Realizamos una predicción con el modelo completo

beginCluster()

## 4 cores detected, using 3

preds_rf_1 <- clusterR(img, raster::predict, args = list(model = modFit_rf_1))
endCluster()

Una vez terminado podemos graficar nuestra clasificación

plot(preds_rf_1,col=cpal, main="Modelo completo (B2 + B3 + B4 + B5 + NDVI)")

Predicción modelo con B5, NDVI y B3 solamente

beginCluster()

## 4 cores detected, using 3

preds_rf_2 <- clusterR(img, raster::predict, args = list(model = modFit_rf_2))
endCluster()

Una vez terminado podemos graficar nuestra clasificación

plot(preds_rf_2,col=cpal, main="B3 + B5 + NDVI")

Predicción modelo sin ndvi ( B2 + B3 + b4 + b5)

beginCluster()

## 4 cores detected, using 3

preds_rf_3 <- clusterR(img, raster::predict, args = list(model = modFit_rf_3))
endCluster()

Una vez terminado podemos graficar nuestra clasificación

plot(preds_rf_3,col=cpal, main="B2 + B3 + B4 + B5")

Predicción modelo con ndvi sin capas que lo componen ( B2 + B3 + ndvi)

beginCluster()

## 4 cores detected, using 3

preds_rf_4 <- clusterR(img, raster::predict, args = list(model = modFit_rf_4))
endCluster()

Una vez terminado podemos graficar nuestra clasificación

plot(preds_rf_4,col=cpal, main= "B2 + B3 + NDVI")

Predicción modelo con nuevas capas

beginCluster()

## 4 cores detected, using 3

preds_rf_5 <- clusterR(img, raster::predict, args = list(model = modFit_rf_5))
endCluster()

Una vez terminado podemos graficar nuestra clasificación

plot(preds_rf_5,col=cpal, main= "B3 + B5 + B6 + NDVI")

Revisamos todos los modelos para resumir su información.

lista_modelos <- ls(pattern="modFit")
resumen_modelos <- data.frame(matrix(ncol=7))
colnames(resumen_modelos) <- c( "Modelo", "Fórmula", "mtry", "Accuracy", "Kappa", "AccuracySD", "KappaSD" )
for (inc in 1:length(lista_modelos)){
  resumen_modelos[inc,1] <- lista_modelos[inc]
  resumen_modelos[inc,2] <-                                    as.character(eval(parse(text=lista_modelos[inc]))$call)[2]
  resumen_modelos[inc,3] <- eval(parse(text=lista_modelos[inc]))$results[nrow(eval(parse(text=lista_modelos[inc]))$results),1]
  resumen_modelos[inc,4] <- round(eval(parse(text=lista_modelos[inc]))$results[nrow(eval(parse(text=lista_modelos[inc]))$results),2],4)
  resumen_modelos[inc,5] <- round(eval(parse(text=lista_modelos[inc]))$results[nrow(eval(parse(text=lista_modelos[inc]))$results),3],4)
  resumen_modelos[inc,6] <- round(eval(parse(text=lista_modelos[inc]))$results[nrow(eval(parse(text=lista_modelos[inc]))$results),4],4)
  resumen_modelos[inc,7] <- round(eval(parse(text=lista_modelos[inc]))$results[nrow(eval(parse(text=lista_modelos[inc]))$results),5],4)
}
print(resumen_modelos)

##        Modelo                                     Fórmula mtry Accuracy  Kappa AccuracySD KappaSD
## 1 modFit_rf_1 as.factor(CLASS) ~ B2 + B3 + B4 + B5 + ndvi    5   0.8586 0.7832     0.0132  0.0196
## 2 modFit_rf_2           as.factor(CLASS) ~ B3 + B5 + ndvi    3   0.8166 0.7197     0.0093  0.0149
## 3 modFit_rf_3        as.factor(CLASS) ~ B2 + B3 + B4 + B5    4   0.8555 0.7791     0.0129  0.0207
## 4 modFit_rf_4           as.factor(CLASS) ~ B2 + B3 + ndvi    3   0.8439 0.7620     0.0080  0.0120
## 5 modFit_rf_5      as.factor(CLASS) ~ B3 + B5 + B6 + ndvi    4   0.8297 0.7386     0.0138  0.0212

3.4 Exportar la clasificación

Exportamos nuestra clasificación para manipularla en nuestro software favorito

rf_1 <- writeRaster(preds_rf_1, filename=file.path(ruta_imagenes,"RF_clasif_1.tif"), overwrite=TRUE)
rf_2 <- writeRaster(preds_rf_2, filename=file.path(ruta_imagenes,"RF_clasif_2.tif"), overwrite=TRUE)
rf_3 <- writeRaster(preds_rf_3, filename=file.path(ruta_imagenes,"RF_clasif_3.tif"), overwrite=TRUE)
rf_4 <- writeRaster(preds_rf_4, filename=file.path(ruta_imagenes,"RF_clasif_4.tif"), overwrite=TRUE)
rf_5 <- writeRaster(preds_rf_5, filename=file.path(ruta_imagenes,"RF_clasif_5.tif"), overwrite=TRUE)

Taller de especialización
Herramientas SIG:
Google Earth Engine y R.
Clase 3.
Modelamiento.

Autor: Nelson Mattie

Fecha: 23 de febrero de 2021

1 Preparar sesión de trabajo.

1.1 Instalar paquetes.

1.2 Definir directorios.

2 Manipulación de datos.

2.1 Cargar rasters.

3 Random Forest

3.1 Cargamos las áreas que nos servirán para entrenar el modelo

3.2 Creación modelo

3.3 Predicción

3.4 Exportar la clasificación

Taller de especialización Herramientas SIG: Google Earth Engine y R. Clase 3. Modelamiento.

Autor: Nelson Mattie

Fecha: 23 de febrero de 2021

1 Preparar sesión de trabajo.

1.1 Instalar paquetes.

1.2 Definir directorios.

2 Manipulación de datos.

2.1 Cargar rasters.

3 Random Forest

3.1 Cargamos las áreas que nos servirán para entrenar el modelo

3.2 Creación modelo

3.3 Predicción

3.4 Exportar la clasificación

Taller de especialización
Herramientas SIG:
Google Earth Engine y R.
Clase 3.
Modelamiento.