Cartografía de Suelos No. 1
Ivan Lizarazo
25 Febrero 2017
Algoritmos de aprendizaje de máquina para cartografía de suelos
Este tutorial está basado en el documento elaborado por T. Hengl para el instituto ISRIC - World Soil Information. En este documento se aplican algunos algoritmos de aprendizaje de máquina (ML en inglés) que pueden ser útiles para proyectos de cartografía de suelos. En particular, se usan algoritmos de alto desempeño tales como bosques aleatorios (random forests), máquinas de soporte vectorial (svm) y redes neuronales artificiales de múltiples capas (deep learning).
Tenga en cuenta que este tutorial no revisa los conceptos teóricos asociados a estos algoritmos y que ellos deben ser estudiados usando otras referencias.
Instalación de los paquetes de R que se requieren
La siguiente función, escrita por Juan Antonio Caro, permite descargar las librerías requeridas de una manera “automágica”:
InstalledPackage <- function(package)
{
available <- suppressMessages(suppressWarnings(sapply(package, require, quietly = TRUE, character.only = TRUE, warn.conflicts = FALSE)))
missing <- package[!available]
if (length(missing) > 0) return(FALSE)
return(TRUE)
}
CRANChoosen <- function()
{
return(getOption("repos")["CRAN"] != "@CRAN@")
}
UsePackage <- function(package, defaultCRANmirror = "http://cran.at.r-project.org")
{
if(!InstalledPackage(package))
{
if(!CRANChoosen())
{
chooseCRANmirror()
if(!CRANChoosen())
{
options(repos = c(CRAN = defaultCRANmirror))
}
}
suppressMessages(suppressWarnings(install.packages(package)))
if(!InstalledPackage(package)) return(FALSE)
}
return(TRUE)
}Ahora usemos el script mágico para usar nuestras librerías:
libraries <- c("plotKML", "sp", "randomForest", "nnet", "e1071", "GSIF",
"plyr","raster", "caret","Cubist", "GSIF", "xgboost")
for(library in libraries)
{
if(!UsePackage(library))
{
stop("Error!", library)
}
}Zona de estudio
Ebergötzen es una zona de 10 km x 10 km que está muy cerca de Göttingen (Alemania). Esta área ha sido estudiada desde hace muchos años, especialmente para el desarrollo de técnicas de cartografía digital del suelo (Gehrt and Böhner, 2001).
La variable eberg contiene 3670 observations (obtenidas con barreno) de textura del suelo a cinco profundidades (0-10, 10-30, 30-50, 50-70, and 70-90), y registros de clases de suelo de acuerdo con el sistema de clasificación de suelos de Alemania.
La variable eberg_grid contiene datos raster con resolución de 100 m que pueden ser usados como covariables para realizar la predicción de variables del suelo.
La variable eberg_grid25 contiene datos raster con resolución de 25 m.
La variable eberg_zones contiene datos vectoriales de tipo poligonal que representan la distribución de material parental en la zona de estudio (Silt and sand, Sandy material, Clayey derivats, Clay and loess).
La variable eberg_contours contine curvas de nivel derivadas del modelo de elevación digital (DEM) con resolución de 25 m, usando un intervalo entre curvas igual a 10 m.
Cargue de datos
Enseguida vamos a cargar en la sesión en R el conjunto de datos de Ebergotzen descrito anteriormente:
data(eberg)
data(eberg_grid)
coordinates(eberg) <- ~X+Y
proj4string(eberg) <- CRS("+init=epsg:31467")
gridded(eberg_grid) <- ~x+y
proj4string(eberg_grid) <- CRS("+init=epsg:31467")Las covariables (variables predictoras) originales se pueden sintetizar en nuevas covariables usando la transformación de componentes principales:
eberg_spc <- spc(eberg_grid, ~ PRMGEO6+DEMSRT6+TWISRT6+TIRAST6)## Converting PRMGEO6 to indicators...## Warning: In prcomp.default(formula = formulaString, x) :
## extra argument 'formula' will be disregarded## Converting covariates to principal components...eberg_grid@data <- cbind(eberg_grid@data, eberg_spc@predicted@data)Todos los análisis que siguen, usan una matriz de regresión (obtenida por medio de una operación espacial de superposición entre los puntos de muestreo y los datos raster) que contiene los valores de la variable de interés lo mismo que los valores de todas las covariables en todos los puntos de entrenamiento:
ov <- over(eberg, eberg_grid)
m <- cbind(ov, eberg@data)
dim(m)## [1] 3670 44El resultado anterior nos indica que la matriz de regresión consiste de 3670 filas (es decir, observaciones) y 44 columnas (valores de variables conocidos en cada observación).
Predicción espacial de las clases de suelos usando ML
En el primer ejemplo nos vamos a enfocar en obtener una clasificación de suelos usando los datos muestreados con barreno. En primer lugar, nos vamos a enfocar solamente en aquellas clases que ocurren con alta frecuencia y que, por ello, permiten realizar un modelamiento estadístico. Una regla aproximada para realizar esa selección es considerar aquellas clases que tengan por lo menos 5 observaciones:
xg = summary(m$TAXGRSC, maxsum=length(levels(m$TAXGRSC)))
str(xg)## Named int [1:13] 790 704 487 376 252 215 186 86 71 23 ...
## - attr(*, "names")= chr [1:13] "Braunerde" "Parabraunerde" "Pseudogley" "Regosol" ...selg.levs = attr(xg, "names")[xg > 5]
m$soiltype <- m$TAXGRSC
m$soiltype[which(!m$TAXGRSC %in% selg.levs)] <- NA
m$soiltype <- droplevels(m$soiltype)
str(summary(m$soiltype, maxsum=length(levels(m$soiltype))))## Named int [1:11] 790 704 487 376 252 215 186 86 71 43 ...
## - attr(*, "names")= chr [1:11] "Braunerde" "Parabraunerde" "Pseudogley" "Regosol" ...Como puede verse, existen 2 clases que tienen menos de 5 observaciones y ellas se han removido para realizar el análisis. Igualmente, también es conveniente remover aquellos puntos que contienen valores nulos (o que no tienen valores) para alguna combinación entre las covariables y la variable de interés:
m <- m[complete.cases(m[,1:(ncol(eberg_grid)+2)]),]
m$soiltype <- as.factor(m$soiltype)Ahora podemos ajustar un modelo predictivo usando la técnica de bosques aleatorios, usando cuatro covariables (material parental, elevación, índice de humedad topográfica -TWI- y banda termal Aster):
s <- sample.int(nrow(m), 500)
TAXGRSC.rf <- randomForest(x=m[-s,paste0("PC",1:10)], y=m$soiltype[-s], xtest=m[s,paste0("PC",1:10)], ytest=m$soiltype[s])
TAXGRSC.rf##
## Call:
## randomForest(x = m[-s, paste0("PC", 1:10)], y = m$soiltype[-s], xtest = m[s, paste0("PC", 1:10)], ytest = m$soiltype[s])
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 3
##
## OOB estimate of error rate: 59.5%
## Confusion matrix:
## Auenboden Braunerde Gley Kolluvisol Parabraunerde
## Auenboden 11 2 6 3 10
## Braunerde 2 290 6 13 59
## Gley 6 7 15 4 13
## Kolluvisol 4 20 5 22 37
## Parabraunerde 6 79 7 23 209
## Pararendzina 0 32 0 1 3
## Pelosol 0 32 0 3 3
## Pseudogley 5 64 5 7 83
## Ranker 0 6 0 0 4
## Regosol 0 66 2 3 42
## Rendzina 0 0 0 0 0
## Pararendzina Pelosol Pseudogley Ranker Regosol Rendzina
## Auenboden 0 2 3 0 0 0
## Braunerde 27 28 58 4 47 4
## Gley 0 0 8 0 4 0
## Kolluvisol 2 2 13 0 3 0
## Parabraunerde 5 5 61 2 27 0
## Pararendzina 61 35 5 2 2 1
## Pelosol 36 45 14 0 3 1
## Pseudogley 7 10 102 0 44 1
## Ranker 1 1 1 0 2 0
## Regosol 6 4 53 2 67 1
## Rendzina 3 7 0 0 1 9
## class.error
## Auenboden 0.7027027
## Braunerde 0.4609665
## Gley 0.7368421
## Kolluvisol 0.7962963
## Parabraunerde 0.5070755
## Pararendzina 0.5704225
## Pelosol 0.6715328
## Pseudogley 0.6890244
## Ranker 1.0000000
## Regosol 0.7276423
## Rendzina 0.5500000
## Test set error rate: 63.2%
## Confusion matrix:
## Auenboden Braunerde Gley Kolluvisol Parabraunerde
## Auenboden 3 1 3 1 2
## Braunerde 0 67 2 1 24
## Gley 0 2 3 0 2
## Kolluvisol 1 5 2 7 11
## Parabraunerde 2 9 4 6 40
## Pararendzina 0 7 0 0 0
## Pelosol 0 6 0 0 1
## Pseudogley 0 27 1 3 19
## Ranker 0 2 0 0 0
## Regosol 0 19 0 2 19
## Rendzina 0 1 0 0 0
## Pararendzina Pelosol Pseudogley Ranker Regosol Rendzina
## Auenboden 0 0 1 0 0 0
## Braunerde 6 4 10 0 17 0
## Gley 0 0 3 0 1 0
## Kolluvisol 0 0 4 0 0 0
## Parabraunerde 2 1 17 0 8 0
## Pararendzina 15 10 1 0 1 0
## Pelosol 11 19 1 1 1 0
## Pseudogley 0 2 19 0 12 0
## Ranker 0 0 0 0 0 0
## Regosol 0 6 10 0 11 0
## Rendzina 1 0 0 0 0 0
## class.error
## Auenboden 0.7272727
## Braunerde 0.4885496
## Gley 0.7272727
## Kolluvisol 0.7666667
## Parabraunerde 0.5505618
## Pararendzina 0.5588235
## Pelosol 0.5250000
## Pseudogley 0.7710843
## Ranker 1.0000000
## Regosol 0.8358209
## Rendzina 1.0000000TAXGRSC.rf$test$confusion[,"class.error"]## Auenboden Braunerde Gley Kolluvisol Parabraunerde
## 0.7272727 0.4885496 0.7272727 0.7666667 0.5505618
## Pararendzina Pelosol Pseudogley Ranker Regosol
## 0.5588235 0.5250000 0.7710843 1.0000000 0.8358209
## Rendzina
## 1.0000000Las variables xtest y ytest contienen una selección de los datos excluyendo 500 puntos. Esas variables permitem ajustar el modelo usando una técnica de validación cruzada. Los resultados muestran un error de predicción relativamente alto (60%). Es decir que la exactitud de clasificación es apenas del 40%.
Adicionalmente, podemos evaluar otros algoritmos ML tales como multinom (librería nnet) y máquinas de soporte vectorial (support vector machines) (librería e1071):
TAXGRSC.rf <- randomForest(x=m[,paste0("PC",1:10)], y=m$soiltype)
fm = as.formula(paste("soiltype~", paste(paste0("PC",1:10), collapse="+")))
fm## soiltype ~ PC1 + PC2 + PC3 + PC4 + PC5 + PC6 + PC7 + PC8 + PC9 +
## PC10TAXGRSC.mn <- multinom(fm, m)## # weights: 132 (110 variable)
## initial value 6119.428736
## iter 10 value 4161.338634
## iter 20 value 4118.296050
## iter 30 value 4054.454486
## iter 40 value 4020.653949
## iter 50 value 3995.113270
## iter 60 value 3980.172669
## iter 70 value 3975.188371
## iter 80 value 3973.743572
## iter 90 value 3973.073564
## iter 100 value 3973.064186
## final value 3973.064186
## stopped after 100 iterationsTAXGRSC.mn## Call:
## multinom(formula = fm, data = m)
##
## Coefficients:
## (Intercept) PC1 PC2 PC3 PC4
## Braunerde 5.45385949 -1.5523725 -1.1414899 3.8709194 1.4091389
## Gley 0.81237203 -0.8441025 -0.9312981 -0.8226797 -1.3941956
## Kolluvisol 2.07610351 -0.5182167 -0.8295361 0.5692995 0.1252458
## Parabraunerde 4.79697987 -0.7988520 -0.9690134 2.8326057 1.2534355
## Pararendzina 0.02502688 -3.9859766 4.3969612 5.3946981 1.0568509
## Pelosol 1.74206393 -4.6008224 -2.6307442 4.5101860 2.3914096
## Pseudogley 5.54640062 -1.6876641 -0.8466514 2.8029280 1.6803954
## Ranker -0.96272445 -4.2793437 -3.6187842 4.3917782 1.4555614
## Regosol 4.75933466 -1.8791610 -1.4315041 3.3350548 1.4528831
## Rendzina -3.79009984 -7.1819850 2.7137859 2.4190389 5.3561776
## PC5 PC6 PC7 PC8 PC9
## Braunerde 1.8267474 -0.53585896 -1.0202426 0.4130277 -2.1117407
## Gley 3.2509685 0.95970334 0.1373794 -1.3089189 -1.1892245
## Kolluvisol 0.1753722 -1.11878357 -3.0079465 -0.2429255 -0.8361207
## Parabraunerde 1.8110834 -0.23800375 -0.3752836 0.3904447 -1.5288602
## Pararendzina 2.1783580 -0.16180762 0.6286308 1.3313932 -1.8354678
## Pelosol -0.5658518 0.55226113 1.1437075 -0.3310972 -4.6781098
## Pseudogley 1.9471449 -0.02731435 -0.1180446 0.4029453 -2.0395526
## Ranker -0.4347812 -0.61573963 -2.3780706 -0.8746946 -6.1753129
## Regosol 2.0765431 -0.49133114 -1.4234241 -0.1413113 -2.2472025
## Rendzina 1.0357045 0.87248537 1.3751377 0.3783946 -5.4131087
## PC10
## Braunerde 4.0354695
## Gley 3.1921390
## Kolluvisol 2.3187566
## Parabraunerde 3.4629867
## Pararendzina 1.5015180
## Pelosol 2.4189837
## Pseudogley 5.5537219
## Ranker -0.6739856
## Regosol 4.2216540
## Rendzina 4.2512969
##
## Residual Deviance: 7946.128
## AIC: 8126.128TAXGRSC.svm <- svm(fm, m, probability=TRUE, cross=5)
TAXGRSC.svm##
## Call:
## svm(formula = fm, data = m, probability = TRUE, cross = 5)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 1
## gamma: 0.1
##
## Number of Support Vectors: 2375TAXGRSC.svm$tot.accuracy## [1] 40.90909Los resultados anteriores muestran que los tres métodos utilizados (random forest, multinom y svm) producen niveles de exactitud similares. Por ello, podemos “fusionar” las predicciones usando un promedio simple (otra alternativa sería usar el paquete caretEnsemble que permite optimizar la integración de varios modelos):
probs1 <- predict(TAXGRSC.mn, eberg_grid@data, type="probs", na.action = na.pass)
probs2 <- predict(TAXGRSC.rf, eberg_grid@data, type="prob", na.action = na.pass)
probs3 <- attr(predict(TAXGRSC.svm, eberg_grid@data, probability=TRUE, na.action = na.pass), "probabilities")Ahora, procedemos a la obtención del promedio:
leg <- levels(m$soiltype)
lt <- list(probs1[,leg], probs2[,leg], probs3[,leg])
probs <- Reduce("+", lt) / length(lt)
eberg_soiltype = eberg_grid
eberg_soiltype@data <- data.frame(probs)Es buena idea verificar que la suma de todas las predicciones es igual al 100%:
ch <- rowSums(eberg_soiltype@data)
summary(ch)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1 1 1 1 1 1Podemos plotear el resultado usando el paquete raster:
plot(raster::stack(eberg_soiltype), col=SAGA_pal[[1]], zlim=c(0,1))
La figura de arriba muestra la clasificación de suelos que se ha obtenido para la zona de Ebergotzen usando los tres modelos.
También es posible obtener un índice de confusión, que representa la incertidumbre de la clasificación, y establecer si es posible agregar algunas clases. También se puede generar una clasificación “cualitativa” si se selecciona para cada pixel aquella clase que tenga la mayor probabilidad:
eberg_soiltype$class <- NULL
summary(eberg_soiltype@data)## Auenboden Braunerde Gley
## Min. :0.001126 Min. :0.02568 Min. :0.001163
## 1st Qu.:0.003243 1st Qu.:0.13300 1st Qu.:0.008016
## Median :0.006650 Median :0.17772 Median :0.014450
## Mean :0.022304 Mean :0.22935 Mean :0.034343
## 3rd Qu.:0.022081 3rd Qu.:0.27390 3rd Qu.:0.035207
## Max. :0.461989 Max. :0.69116 Max. :0.710280
## Kolluvisol Parabraunerde Pararendzina
## Min. :0.002873 Min. :0.004612 Min. :0.001165
## 1st Qu.:0.023987 1st Qu.:0.083156 1st Qu.:0.002385
## Median :0.035491 Median :0.122953 Median :0.011258
## Mean :0.050266 Mean :0.184962 Mean :0.064523
## 3rd Qu.:0.051241 3rd Qu.:0.286666 3rd Qu.:0.040572
## Max. :0.450358 Max. :0.625822 Max. :0.662323
## Pelosol Pseudogley Ranker
## Min. :0.002002 Min. :0.01782 Min. :0.001376
## 1st Qu.:0.006346 1st Qu.:0.07926 1st Qu.:0.002947
## Median :0.013109 Median :0.12427 Median :0.004253
## Mean :0.059944 Mean :0.15993 Mean :0.006627
## 3rd Qu.:0.033575 3rd Qu.:0.18642 3rd Qu.:0.007135
## Max. :0.537763 Max. :0.61842 Max. :0.234746
## Regosol Rendzina
## Min. :0.008785 Min. :0.0003582
## 1st Qu.:0.034156 1st Qu.:0.0006643
## Median :0.073967 Median :0.0010215
## Mean :0.101642 Mean :0.0861144
## 3rd Qu.:0.147030 3rd Qu.:0.1496042
## Max. :0.572929 Max. :0.4378668DF <- eberg_soiltype@data
eberg_soiltype$class <- colnames(DF)[max.col(DF, ties.method = 'first')]
unique(eberg_soiltype$class) ## [1] "Pararendzina" "Pelosol" "Rendzina" "Braunerde"
## [5] "Parabraunerde" "Regosol" "Auenboden" "Kolluvisol"
## [9] "Pseudogley" "Gley"colors <- heat.colors(10)
plot(raster(eberg_soiltype[12]), col=colors, zlim=c(1,10))
Modelación numérica de las propiedades del suelo usando h2o
El algoritmo de bosques aleatorios (random forest) es útil para resolver problemas de clasificación y de regresión. Por ello, se puede aplicar en el modelamiento de propiedades del suelo (es decir, para generar predicciones de valores cuantitativos). El paquete randomForest de R puede no ser la mejor alternativa cuando se usan datos “muy pesados”. Una opción más interesante es la utilización del paquete h2o (Richter et al. 2015). h2o es una implementación basada en Java y todos los cálculos se realizan fuera de R, es decir dentro de JVM (Java Virtual Machine).
La parte más complicada es la instalación de h2o. Una guía para realizar ese proceso se encuentra en este link: http://h2o-release.s3.amazonaws.com/h2o/master/1247/docs-website/Ruser/R_studio.html
En el ejemplo que sigue vamos a cartografiar el contentido de arena para el horizonte superior. Empecemos por inicializar h2o:
#install.packages("h2o")
library(h2o)##
## ----------------------------------------------------------------------
##
## Your next step is to start H2O:
## > h2o.init()
##
## For H2O package documentation, ask for help:
## > ??h2o
##
## After starting H2O, you can use the Web UI at http://localhost:54321
## For more information visit http://docs.h2o.ai
##
## ----------------------------------------------------------------------##
## Attaching package: 'h2o'## The following objects are masked from 'package:raster':
##
## %in%, as.factor, is.factor## The following objects are masked from 'package:stats':
##
## cor, sd, var## The following objects are masked from 'package:base':
##
## %*%, %in%, &&, ||, apply, as.factor, as.numeric, colnames,
## colnames<-, ifelse, is.character, is.factor, is.numeric, log,
## log10, log1p, log2, round, signif, trunclocalH2O = h2o.init()## Connection successful!
##
## R is connected to the H2O cluster:
## H2O cluster uptime: 17 hours 42 minutes
## H2O cluster version: 3.10.3.6
## H2O cluster version age: 3 days
## H2O cluster name: H2O_started_from_R_Ivan_ldv605
## H2O cluster total nodes: 1
## H2O cluster total memory: 3.43 GB
## H2O cluster total cores: 8
## H2O cluster allowed cores: 2
## H2O cluster healthy: TRUE
## H2O Connection ip: localhost
## H2O Connection port: 54321
## H2O Connection proxy: NA
## R Version: R version 3.3.2 (2016-10-31)#h2o.shutdown()
#h2o.init(nthreads = -1)Enseguida, vamos a preparar la matriz de regresión y los sitios de las predicciones usando la función as.h2o de manera que sean visibles a la librería h2o:
eberg.hex = as.h2o(m, destination_frame = "eberg.hex")##
|
| | 0%
|
|=================================================================| 100%eberg.grid = as.h2o(eberg_grid@data, destination_frame = "eberg.grid")##
|
| | 0%
|
|=================================================================| 100%Ahora podemos ajustar un modelo random forest:
RF.m = h2o.randomForest(y = which(names(m)=="SNDMHT_A"), x = which(names(m) %in% paste0("PC",1:10)), training_frame = eberg.hex, ntree = 50)##
|
| | 0%
|
|==== | 6%
|
|============================================ | 68%
|
|=================================================================| 100%RF.m## Model Details:
## ==============
##
## H2ORegressionModel: drf
## Model ID: DRF_model_R_1487981499883_28
## Model Summary:
## number_of_trees number_of_internal_trees model_size_in_bytes min_depth
## 1 50 50 640844 20
## max_depth mean_depth min_leaves max_leaves mean_leaves
## 1 20 20.00000 923 1063 1016.36000
##
##
## H2ORegressionMetrics: drf
## ** Reported on training data. **
## ** Metrics reported on Out-Of-Bag training samples **
##
## MSE: 220.5606
## RMSE: 14.85128
## MAE: 10.05702
## RMSLE: 0.4263087
## Mean Residual Deviance : 220.5606Los resultados muestran que el valor R2 del modelo ajustado es cercano al 50%. Ello también se puede observar en un gráfico de valores observados vs valores obtenidos en la predicción:
plot(m$SNDMHT_A, as.data.frame(h2o.predict(RF.m, eberg.hex, na.action=na.pass))$predict, asp=1, xlab="measured", ylab="predicted")##
|
| | 0%
|
|=================================================================| 100%
Para obtener un mapa de contenido de arena basado en las predicciones se puede usar:
eberg_grid$RFx <- as.data.frame(h2o.predict(RF.m, eberg.grid, na.action=na.pass))$predict##
|
| | 0%
|
|=================================================================| 100%plot(raster(eberg_grid["RFx"]), col=SAGA_pal[[1]], zlim=c(10,90))
points(eberg, pch="+")
La librería h2o tiene otro algoritmo ML de interés que se conoce como “deep learning” y que no es otra cosa que es una red neuronal artificial de múltiples capas. El ajuste del modelo es similar al uso random forest:
DL.m <- h2o.deeplearning(y = which(names(m)=="SNDMHT_A"), x = which(names(m) %in% paste0("PC",1:10)), training_frame = eberg.hex)##
|
| | 0%
|
|====== | 10%
|
|==================== | 30%
|
|============================================== | 70%
|
|=================================================================| 100%DL.m## Model Details:
## ==============
##
## H2ORegressionModel: deeplearning
## Model ID: DeepLearning_model_R_1487981499883_29
## Status of Neuron Layers: predicting SNDMHT_A, regression, gaussian distribution, Quadratic loss, 42.601 weights/biases, 508,2 KB, 25.520 training samples, mini-batch size 1
## layer units type dropout l1 l2 mean_rate rate_rms
## 1 1 10 Input 0.00 %
## 2 2 200 Rectifier 0.00 % 0.000000 0.000000 0.016178 0.011101
## 3 3 200 Rectifier 0.00 % 0.000000 0.000000 0.161467 0.193774
## 4 4 1 Linear 0.000000 0.000000 0.001413 0.001068
## momentum mean_weight weight_rms mean_bias bias_rms
## 1
## 2 0.000000 0.000157 0.102569 0.334967 0.069185
## 3 0.000000 -0.018338 0.071091 0.952872 0.020253
## 4 0.000000 -0.002110 0.042361 0.062960 0.000000
##
##
## H2ORegressionMetrics: deeplearning
## ** Reported on training data. **
## ** Metrics reported on full training frame **
##
## MSE: 230.3774
## RMSE: 15.17819
## MAE: 10.99303
## RMSLE: 0.4484678
## Mean Residual Deviance : 230.3774Los resultados obtenidos mediante deep learning son comparables a los obtenidos mediante random forest. Sin embargo, los resultados de las predicciones (map) lucen diferente de las predicciones obtenidas mediante random forest:
eberg_grid$DLx <- as.data.frame(h2o.predict(DL.m, eberg.grid, na.action=na.pass))$predict##
|
| | 0%
|
|=================================================================| 100%plot(raster(eberg_grid["DLx"]), col=SAGA_pal[[1]], zlim=c(10,90))
points(eberg, pch="+")
La pregunta obvia es cuál de los dos métodos deberíamos usar? Teniendo en cuenta que los dos tienen un desempeño similar, la opción lógica es fusionar (emsamblar) las predicciones, de manera que cada método sea utilizado de acuerdo con algún sistema de ponderación. En este ejemplo, se puede usar el valor R2 para realizar esa ponderación:
rf.R2 = RF.m@model$training_metrics@metrics$r2
dl.R2 = DL.m@model$training_metrics@metrics$r2
eberg_grid$SNDMHT_A <- rowSums(cbind(eberg_grid$RFx*rf.R2, eberg_grid$DLx*dl.R2), na.rm=TRUE)/(rf.R2+dl.R2)
plot(raster(eberg_grid["SNDMHT_A"]), col=SAGA_pal[[1]], zlim=c(10,90))
points(eberg, pch="+")
La figura de arriba muestra la predicción del contenido de arena basada en los modelos integrados. De hecho, los resultados muestran patrones asociados a los métodos individuales pero tienen una mejor exactitud que cualquiera de ellos.
Resumen
Este documento ha ilustrado el uso de algunos algoritmos de aprendizaje de máquina (ML) aplicado a la cartografía de suelos, tanto de propiedades cualitativas como de propiedades cuantitativas. Los algoritmos ML son muy atractivos en cuanto permiten representar de una manera relaciones que son eminentemente locales, como ocurre a menudo con las propiedades de suelos. Los algoritmos ML no dependen en supuestos estadísticos ideales sino que explotan las características particulares de los datos disponibles. Sin embargo, un problema con estos algoritmos es que requieren una gran capacidad computacional y que los modelos obtenidos son eminentemente empíricos y, por consiguiente, únicamente válidos en las zonas geográficas específicas en que han sido estudiados.
Referencias
Hengl, T. (n.d.) Machine Learning Algorithms for Soil Mapping. Available at http://gsif.isric.org/doku.php?id=wiki:soilmapping_using_mla