Clasificación de las Coberturas de la Tierra para la Subzona Hidrológica del Río Sogamoso
Diego Cortés Ramos
4/12/2021
1 Introducción
Este cuaderno de R (R Core Team 2020) contiene la clasificación de las coberturas de la tierra de la subzona hidrológica (SZH) del río Sogamoso. Apoyado en imágenes del satélite Landsat 8 del 04 de enero de 2015, escogida por su baja presencia de nubes, y utilizando la información proporcionada por el IDEAM para la cobertura de la tierra (IDEAM 2014) en el periodo 2010-2012 para el entrenamiento y validación del algoritmo de clasificación.
Se escogió el algoritmo Random Forest para el entrenamiento de la clasificación, y se utilizó un muestreo aleatorio estratificado de la información de referencia, la cual corresponde con las coberturas clasificadas por el IDEAM, mencionadas anteriormente.
Para la validación se calculó la matriz de confusión, con la cual se obtuvo la exactitud de productor y de usuario, además de la exactitud global y el coeficiente Kappa.
rm(list = ls())
library(sp)
library(terra)
library(raster)
library(dplyr)
library(ranger)
library(corrplot)
library(bookdown)
library(plotrix)
library(leaflet)2 Metodos
A continuación, se hará una breve descripción de la zona de estudio, la información utilizada y el los métodos empleados para el desarrollo del presente informe.
2.1 Zona de estudio
El IDEAM realizó la zonificación hidrográfica el país (2013) en tres tipos de unidades, las áreas hidrograficas, correspondientes a las 5 macrocuencas que componen el país (Magdalena-Cauca, Orinoco, Pacifico, Caribe y Amazonas), de las cuales la Magdalena-Cauca (McMC) es la que concentra la mayor parte de la población y economía de Colombia. Dentro de esta clasificación se ubican las zonas y subzonas hidrográficas, siendo estas ultimas las de menor extensión.
La zona hidrográfica del Medio Magdalena, comprende el tramo del río Magdalena entre los municipios de Honda (Tolima) y El Banco (Magdalena), abarcando los primeros pozos de explotación de hidrocarburos (Semana 2018). Entre los municipios de Barrancabermeja y Puerto Wilches se produce la entrega del río Sogamoso el cual es formado por la confluencia de los ríos Chicamocha y Suárez, generando una cuenca de mas de 23 000 km2. La SZH del río sogamoso se delimita desde este punto de confluencia y hasta la desembocadura en el Magdalena (Figura 2.1) abarcando un área de 3408 km2.
Actualmente, el país se encuentra debatiendo sobre la viabilidad de la implementación del fracking en el territorio. Para lo cual, se ha planteado la contrucción de “Proyectos Pilotos de Investigación,” el primero de ellos en Puerto Wilches (El Tiempo 2020).
Es por esto que conocer las dinámicas del recurso hídrico en cercanias del proyecto es de suma importancia para estudiar los verdaderos impactos de este tipo de proyectos. Y siendo las coberturas de la tierra uno de los principales insumos para el estudio del ciclo hidrológico mediante modelaciones numéricas, se ha desarrollado la presente clasificación para ser usada como insumo en futuras modelaciones hidrológicas.
workdir <- "D:/Documentos/Diego/03.Academico/MAESTRIA/PERCEPCION REMOTA"
sog <- vect(sprintf("%s/GIS/SHP/SZH_Sogamoso_32618.shp", workdir))
sog32618 <- spTransform(as(sog, "Spatial"),
CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
ext.sog <- ext(sog32618)
leaflet(sog32618) %>% addTiles() %>%
setView(lat = mean(ext.sog[3:4]), lng = mean(ext.sog[1:2]), zoom = 09) %>%
addPolygons(color = "#830000", fillOpacity = 0.1, popup = "SZH: Río Sogamoso")Figure 2.1: Área de estudio: Subzona hidrográfica del río Sogamoso
2.2 Información de entrada
Para desarrollar la presente clasificación, se cuenta con información del satelite Landsat 8 y la información de referencia aportada por el IDEAM (2014). Para el procesamiento de esta información se va a contar con el apoyo de la librería terra (Hijmans 2021) y se va a trabajar en el sistema coordenado WGS 84 / UTM zone 18N (EPSG: 32618), al cual fueron proyectadas las capas vectoriales de manera previa. A continuación, se muestran las variables en las cuales se guarda la información mencionada en, y en la Figura 2.2 se muestra una composición en color verdadero de las imágenes Landsat junto con la SZH del río Sogamoso. Si se desea conocer la información, se puede consultar en el link a continuación: Información de entrada
files <- list.files(sprintf("%s/GIS/LANDSAT_8_C1/LC08_L1TP_008055_20150104_20170415_01_T1/", workdir),
"B[1-7]\\.TIF$", full.names = T)
landsat <- rast(files)
names(landsat) <- c('ultra-blue', 'blue', 'green', 'red', 'NIR', 'SWIR1', 'SWIR2')
ref.data <- vect(sprintf("%s/GIS/SHP/COBERTURA_Sogamoso.shp", workdir))message("landsat")## landsatlandsat## class : SpatRaster
## dimensions : 7731, 7571, 7 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 514185, 741315, 683085, 915015 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## sources : LC08_L1TP_008055_20150104_20170415_01_T1_B1.TIF
## LC08_L1TP_008055_20150104_20170415_01_T1_B2.TIF
## LC08_L1TP_008055_20150104_20170415_01_T1_B3.TIF
## ... and 4 more source(s)
## names : ultra-blue, blue, green, red, NIR, SWIR1, ...message("\nsog")##
## sogsog## class : SpatVector
## geometry : polygons
## dimensions : 1, 11 (geometries, attributes)
## extent : 618777, 718699.2, 734238.5, 808470.1 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## names : OBJECTID_1 COD_AH COD_ZH COD_SZH NOM_AH NOM_ZH
## type : <num> <num> <num> <num> <chr> <chr>
## values : 75 2 24 2405 Magdalena Cauca Sogamoso
## NOM_SZH Shape_Leng Shape_Area RULEID Area
## <chr> <num> <num> <num> <num>
## RÃo Sogamoso 3.324 0.2789 3 3.408e+09message("\nref.data")##
## ref.dataref.data## class : SpatVector
## geometry : polygons
## dimensions : 1754, 5 (geometries, attributes)
## extent : 618777, 718699.2, 734238.5, 808470.1 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## names : OBJECTID_1 APOYO INSUMO CODIGO
## type : <num> <chr> <chr> <num>
## values : 1.413e+05 Rapideye 126213 201~ 411
## 1.42e+05 Rapideye 126214 201~ 3221
## 1.42e+05 Rapideye 126214 201~ 3222
## LEYENDA3N
## <chr>
## 4.1.1. Zonas Pantan~
## 3.2.2. Arbustal
## 3.2.2. ArbustalplotRGB(subset(landsat, c("red", "green", "blue")), axes = TRUE, stretch = "lin",
main = "Composición en color verdadero - Niveles digitales Landsat 8")
plot(sog, add = T, border = "#830000", lwd = 2)
Figure 2.2: Composición RGB432 con los niveles digitales de Landsat. Se presenta el área de interés de la SZH del río Sogamoso.
2.3 Preproceso
Inicialmente, se recortaron y enmascararon las imagenes a la extensión de la SZH, se pasaron los niveles digitales de las imagenes satelitales a reflectancia TOA y Se calcularon los indices normalized difference (NDVI), normalized difference water index (NDWI) y Enhanced Vegetation Index (EVI) (Figura 2.3), mediante la definicion de la función spec.trans.
spec.trans <- function(x, ind){
x.ret <- x
for(i in ind){
r <- switch (i,
"NDVI" = (x[["NIR"]] - x[["red"]]) / (x[["NIR"]] + x[["red"]]),
"NDWI" = (x[["NIR"]] - x[["SWIR1"]]) / (x[["NIR"]] + x[["SWIR1"]]),
"EVI" = 2.5 * (x[["NIR"]] - x[["red"]]) / (x[["NIR"]] + 6 * x[["red"]] - 7.5 * x[["blue"]] + 1)
)
nms <- c(names(x.ret), i)
x.ret <- c(x.ret, r)
names(x.ret) <- nms
}
return(x.ret)
}
spec.trans## function(x, ind){
## x.ret <- x
## for(i in ind){
## r <- switch (i,
## "NDVI" = (x[["NIR"]] - x[["red"]]) / (x[["NIR"]] + x[["red"]]),
## "NDWI" = (x[["NIR"]] - x[["SWIR1"]]) / (x[["NIR"]] + x[["SWIR1"]]),
## "EVI" = 2.5 * (x[["NIR"]] - x[["red"]]) / (x[["NIR"]] + 6 * x[["red"]] - 7.5 * x[["blue"]] + 1)
## )
## nms <- c(names(x.ret), i)
## x.ret <- c(x.ret, r)
## names(x.ret) <- nms
## }
## return(x.ret)
## }landsat <- terra::crop(landsat, sog); landsat <- mask(landsat, sog)
toa <- landsat * 2e-5 - 0.1
toa <- spec.trans(toa, c("NDVI", "EVI", "NDWI"))
message("toa")## toatoa## class : SpatRaster
## dimensions : 2475, 3331, 10 (nrow, ncol, nlyr)
## resolution : 30, 30 (x, y)
## extent : 618765, 718695, 734235, 808485 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## sources : memory (7 layers)
## memory
## memory
## ... and 1 more source(s)
## names : ultra-blue, blue, green, red, NIR, SWIR1, ...
## min values : -0.1000000, -0.1000000, -0.1000000, -0.1000000, -0.1000000, -0.1000000, ...
## max values : 0.676700, 0.705060, 0.714840, 0.759120, 0.839180, 1.210700, ...message("\nCapas contenidas en toa")##
## Capas contenidas en toanames(toa)## [1] "ultra-blue" "blue" "green" "red" "NIR"
## [6] "SWIR1" "SWIR2" "NDVI" "EVI" "NDWI"# Se guarda dado que se considera una información de interés por sí misma
writeRaster(toa, sprintf("%s/GIS/RASTER/toa_sogamoso.tif", workdir), overwrite = T)
plotRGB(subset(toa, c("red", "green", "blue")), axes = TRUE, stretch = "lin",
main = "Composición en color verdadero - Reflectancia TOA")
plot(sog, add = T, border = "#830000", lwd = 2)
Figure 2.3: Composicion RGB432 de la reflectancia TOA de las imagenes Landsat recortadas a la SZH de estudio.
2.4 Metodos
Para el desarrollo del presente informe se realizaron dos etapas posteriores a la carga y el preprocesamiento de la información ya realizada. Estás etápas fueron guiadas principalmente por las notas del curso de Percepción Remota, de Universidad Nacional de Colombia.
Inicialmente se desarrolló un análisis estadístico (Lizarazo 2021d), en el cual se calcularon los valores de desviación estándar medios, mínimos y máximos para todas las capas raster en análisis (7 bandas del sensor Landsat 8 y 3 indices espectrales). Posteriormente se desarrolló la clasificación de coberturas (Lizarazo 2021b, Lizarazo2021b) mediante el uso del algoritmo random forest Wright and Ziegler (2017) y se evaluó su exactitud (Lizarazo 2021a) tanto en la etapa de entrenamiento como en la de validación.
A su vez, en la clasificación se pueden identificar tres etapas. La primera, consiste en hacer un muestreo aleatorio estratificado sobre los datos de referencia para entrenar el algoritmo de clasificación; la segunda, es la clasificación en si misma, en donde con el algoritmo entrenado se clasifica cada uno de los píxeles de las imágenes recortadas; y la tercera, consiste en realizar el análisis de la exactitud del algoritmo, haciendo un nuevo muestreo y comparando los resultados de la clasificación con los de la información de referencia.
3 Resultados
3.1 Análisis estadístico
Se realizó el análisis estadístico de las imagenes TOA recortadas, utilizando la guia del cuaderno de R Uniband and multiband image statistics (Lizarazo 2021d) . En algunos casos, el análisis se hizo de manera diferenciada para las bandas del sensor y para los indices calculados, esto debido a que estos últimos no tienen los mismos rangos.
Inicialmente, se calcularon algunos parámetros estadísticos, como lo son los valores mínimos y máximos, la media y desviación estándar, los cuales se muestran en la Tabla 3.1. Es notable como la banda del infrarrojo cercano es la que presenta una mayor desviación estándar, por lo cual se espera que sea una banda con un alto contenido de información.
sts <- do.call(rbind, lapply(c("min", "max", "mean", "sd"), function(x)r <- global(toa, x, na.rm = T) %>% t))
rownames(sts) <- c("Minimo", "Máximo", "Medio", "Desv. Estandar")
knitr::kable(sts, digits = 3, caption = 'Estadística básicas')| ultra-blue | blue | green | red | NIR | SWIR1 | SWIR2 | NDVI | EVI | NDWI | |
|---|---|---|---|---|---|---|---|---|---|---|
| Minimo | -0.100 | -0.100 | -0.100 | -0.100 | -0.100 | -0.100 | -0.100 | -0.827 | -0.278 | -0.764 |
| Máximo | 0.677 | 0.705 | 0.715 | 0.759 | 0.839 | 1.211 | 1.211 | 3.709 | 10.154 | 8.837 |
| Medio | 0.084 | 0.067 | 0.058 | 0.043 | 0.217 | 0.121 | 0.055 | 0.643 | 0.443 | 0.300 |
| Desv. Estandar | 0.017 | 0.017 | 0.020 | 0.024 | 0.064 | 0.050 | 0.034 | 0.202 | 0.157 | 0.159 |
Adicionalmente, se calcularon algunas estadísticas multibanda. En este punto se difereció entre las bandas de los sensores Landast y los indices calculados posteriormente. En la Figura 3.1 se puede ver una alta correlación entre las bandas del visible y una aún mayor entre las bandas ultra blue y blue. En las bandas dentro del infrarojo se observa una menor correlación con las del visible, incluso entre ellas mismas. Las bandas NIR y SWIR2 se presenta una de las menores correlaciones, solamente superada por la que se dan entre las bandas red y NIR.
Para los indices se realizó un análisis similar (Figuras 3.2 y 3.3), se encontró que no una alta correlación entre NDVI y EVI, naturalmente, al tratarse de indices enfocados en la vegetación. Entre estos y el indice NDWI no se encontraron mayores correlaciones, como era esperado. Mediante los histogramas (Figura 3.2) se puede observar la distribución de los valores de cada uno de los indices.
pairs(toa[[1:7]], hist = T, cor = T)
Figure 3.1: Correlaciones entre las bandas espectrales de las imagenes de refelctancia TOA.
par(mfrow = c(1, 3))
for (i in 8:10){
hist(toa[[i]], plot=TRUE, ylab="Frecuencia", col="#09A1E8", main = names(toa)[i])
}
Figure 3.2: Histogramas de los indices espectrales calculados.
cmat <- cor(as.matrix(toa[[8:10]]), use ="pairwise.complete.obs")
ncmat <- round(cmat, 2)
corrplot(ncmat, type = "upper", order = "hclust",
tl.col = "black", tl.srt = 45)
Figure 3.3: Correlaciones entre los indices espectrales calculados.
3.2 Clasificación de las coberturas de la tierra
El proceso de clasificación se realizó en dos etapas. La primera corresponde con el entrenamiento del algortimo random forest (Breiman 2001) a partir de un muestro aleatorio estratificado de la capa de coberturas del IDEAM. Posteriormente se realizó la validación y el análisis de la exactitud. Para esta sección se utilizaron las guías A tutorial on land cover classification from satellite imagery using R y A tutorial on pixel-based land cover classification using random forests in R (Lizarazo 2021b, 2021c)
Random forest es un algoritmo de machine learning que enlaza multiples arboles de decisión mediante un vector aleatorio para luego combinar la clasificación de cada uno de estos y ofrecer una clasificación mas certera. Es un algoritmo que suele mostrar muy buenos resultados, además de ser eficiciente ante grandes bases de datos.
Para el proceso de clasificación se tuvo en cuenta el primer nivel de las coberturas, el cual corresponde a las categorias de: 1) zonas urbanas y mineras, 2) pastos y cultivos, 3) bosques y áreas naturales, 4) zonas humedas y pantanosas, 5) superficies de agua y 9) nubes.
code <- as.character(ref.data$CODIGO) %>% substr(1,1) %>% as.numeric
ref.data$code <- code
head(ref.data)## OBJECTID_1 APOYO INSUMO CODIGO
## 1 141287 Rapideye 126213 2011/08/08 411
## 2 141976 Rapideye 126214 2011/09/05 3221
## 3 142005 Rapideye 126214 2011/09/05 3222
## 4 142138 Rapideye 126214 2011/09/05 3221
## 5 142177 Rapideye 126214 2011/09/05 321111
## 6 142179 Rapideye 97183 2010/02/04 Rapideye 126214 2011/09/05 321111
## LEYENDA3N code
## 1 4.1.1. Zonas Pantanosas 4
## 2 3.2.2. Arbustal 3
## 3 3.2.2. Arbustal 3
## 4 3.2.2. Arbustal 3
## 5 3.2.1. Herbazal 3
## 6 3.2.1. Herbazal 33.2.1 Entrenamiento
El entrenamiento se hizo medainte un muestreo aleatorio estratificado, con 350 puntos por cada uno. A partir de este se creó una base de datos (train.sampdata) que contiene la clasificación de referencia y los valores en cada una de las bandas (incluyendo los indices espectrales) para el total de puntos generados. Posteriormente, se procede al entrenamiento del algoritmo random forest para lo cual se usa la función ranger (Wright and Ziegler 2017). El modelo de predicción generado se almacena en la variable fit.
set.seed(2303)
train.ptsamp <- spatSample(ref.data, size = 350, method = "random", strata = "code", chess = "")
train.nptsamp <- as(train.ptsamp, "Spatial")
train.nref.data <- as(ref.data, "Spatial")
train.ptsamp$class <- over(train.nptsamp, train.nref.data)$code
xy <- as.matrix(geom(train.ptsamp)[,c('x','y')])
df <- extract(toa, xy)
train.sampdata <- data.frame(class = as.factor(train.ptsamp$class), df)
train.sampdata <- na.omit(train.sampdata)
message("train.sampdata")## train.sampdatahead(train.sampdata)## class ultra.blue blue green red NIR SWIR1 SWIR2 NDVI
## 1 1 0.11906 0.10794 0.10418 0.11210 0.15248 0.16854 0.13574 0.1526192
## 2 1 0.09642 0.08104 0.07292 0.06558 0.22310 0.14184 0.08478 0.5456561
## 3 1 0.12768 0.11804 0.11070 0.11554 0.17200 0.21034 0.17440 0.1963553
## 4 1 0.12338 0.11130 0.10922 0.10798 0.18982 0.17684 0.12508 0.2748153
## 5 1 0.10802 0.09114 0.08424 0.06840 0.27940 0.16574 0.10138 0.6066705
## 6 1 0.10404 0.08868 0.07470 0.06382 0.19020 0.14786 0.08858 0.4975199
## EVI NDWI
## 1 0.09940622 -0.05002804
## 2 0.39037253 0.22266674
## 3 0.14403943 -0.10027724
## 4 0.20399821 0.03540064
## 5 0.52422360 0.25533540
## 6 0.34795489 0.12524404fit <- ranger(class~., data = train.sampdata, importance="impurity")
message("fit")## fitfit## Ranger result
##
## Call:
## ranger(class ~ ., data = train.sampdata, importance = "impurity")
##
## Type: Classification
## Number of trees: 500
## Sample size: 2094
## Number of independent variables: 10
## Mtry: 3
## Target node size: 1
## Variable importance mode: impurity
## Splitrule: gini
## OOB prediction error: 35.48 %En este punto es posible realizar un análisis de exactitud del entrenamiento del modelo (Lizarazo 2021c, 2021a). Para esto es importante primero definir la funcion “ovAcc,” encargada de evaluar la exactitud de las predicciones. Posteriormente se determina la matrix de confusión (Tabla 3.2) y se evalua la exactitud global, la exactitud de usuario y de productor y el coeficiente kappa (Tabla 3.3).
ovAcc <- function(conmat)
{
# number of total cases/samples
n = sum(conmat)
# number of correctly classified cases per class
diag = diag(conmat)
# Overall Accuracy
OA = sum(diag) / n
#
# observed (true) cases per class
rowsums = apply(conmat, 1, sum)
p = rowsums / n
# predicted cases per class
colsums = apply(conmat, 2, sum)
q = colsums / n
expAccuracy = sum(p*q)
kappa = (OA - expAccuracy) / (1 - expAccuracy)
#
# Producer accuracy
PA <- diag / colsums
# User accuracy
UA <- diag / rowsums
outAcc <- data.frame(producerAccuracy = PA, userAccuracy = UA)
#
global_acc = data.frame(overallAccuracy=OA, overallKappa=kappa)
return(list(global_acc, outAcc))
}pred <- predict(fit, data = train.sampdata , type = "response")
conmat <- table(pred = pred$predictions, obs = train.sampdata$class)
conmat.df <- as.data.frame.matrix(conmat)
knitr::kable(conmat.df, digits = 3, caption = 'Matríz de confusión del entrenamiento')| 1 | 2 | 3 | 4 | 5 | 9 | |
|---|---|---|---|---|---|---|
| 1 | 344 | 0 | 0 | 0 | 0 | 0 |
| 2 | 0 | 350 | 0 | 0 | 0 | 0 |
| 3 | 0 | 0 | 350 | 0 | 0 | 0 |
| 4 | 0 | 0 | 0 | 350 | 0 | 0 |
| 5 | 0 | 0 | 0 | 0 | 350 | 0 |
| 9 | 0 | 0 | 0 | 0 | 0 | 350 |
train.cart.acc <- ovAcc(conmat)
message(sprintf("Exactitud global: %.3g\nCoeficiente Kappa: %.3g",
train.cart.acc[[1]][1], train.cart.acc[[1]][2]))## Exactitud global: 1
## Coeficiente Kappa: 1colnames(train.cart.acc[[2]]) <- c("Productor", "Usuario")
knitr::kable(train.cart.acc[[2]], digits = 3, caption = 'Exactitud de productor y usuario para el entrenamiento')| Productor | Usuario | |
|---|---|---|
| 1 | 1 | 1 |
| 2 | 1 | 1 |
| 3 | 1 | 1 |
| 4 | 1 | 1 |
| 5 | 1 | 1 |
| 9 | 1 | 1 |
Pese a que se llega a una exactitud del 100% en el entrenamiento, puede que esto no sea un valor confiable y por el contrario indique la existencia de problemas de sobre-entrenamiento. Es por esto que el resultado de exactitud mas importante se obtendrá en la validación.
3.2.2 Clasificación
Finalmente, se procede a la clasificación, en donde a partir del modelo generado se obtiene el raster clasificado de las coberturas de la tierra (Figura 3.4). Adicionalmente, se pueden analizar las frecuencias de las coberturas mediante la Figura 3.5, en donde se ven como predominan los pastos y áreas de cultivo, seguido por bosques y áreas naturales.
classified <- predict(toa, fit, na.rm = T,
fun = function(model, ...) predict(model, ...)$predictions)## Aggregating predictions.. Progress: 86%. Estimated remaining time: 4 seconds.x.labels <- c(sprintf("Zonas urbanas y\nmineras"), "Pastos y cultivos", sprintf("Bosques y áreas\nnaturáles"), sprintf("Zonas humedas y\npantanosas"),
"Superficies de agua", "Nubes")
par(omi = c(0, 0, 0, 2.3))
plot(classified, legend = FALSE, col = c("#da0000ff", "#e6e431ff", "#0e7a11ff", "#60deb1cc", "#5da4fcff", "#dbeeffff"))
par(omi = c(0, 0, 0, 0), mai = c(0,0,0,0))
legend("right", x.labels,
fill = c("#da0000ff", "#e6e431ff", "#0e7a11ff", "#60deb1cc", "#5da4fcff", "#dbeeffff"))
Figure 3.4: Clasificación de las coberturas para la SZH del río Sogamoso
hist(classified$lyr1, plot = T, main = "Histograma de las coberturas clasificadas.",
col = c("#da0000ff", "#e6e431ff", "#0e7a11ff", "#60deb1cc", "#5da4fcff", "#dbeeffff"),
breaks = 0:6, ylab = "Frecuencia", xlab = "Clase")
Figure 3.5: Histograma de las coberturas clasificadas
3.2.3 Evaluación de la exactitud
El proceso de evaluación de la exactitud para la validación del modelo de clasificación se hizo seleccionando una nueva nube de puntos, en donde se compara la predicción realizada con la información de referencia. La nube de puntos se tomó con un nuevo muestreo estratificado, tomando 100 muestras por cada uno. Luego de realizar la extracción de los valores predichos y de referencia para cada muestra, se calculó la matriz de confusión (Tabla 3.4) y la exactitud global, el coeficiente capa y la exactitud de usuario y productor (Tabla 3.5).
set.seed(2908)
val.ptsamp <- spatSample(ref.data, size = 100, method = "random", strata = "code", chess = "")
val.ptsamp$reference <- as.factor(val.ptsamp$pol.id)
prediction <- extract(classified, val.ptsamp)
val.ptsamp$prediction <- as.factor(prediction$lyr1)
conmat <- table(val.ptsamp$prediction, val.ptsamp$reference)conmat.df <- as.data.frame.matrix(conmat[-nrow(conmat),])
knitr::kable(conmat.df, digits = 3, caption = 'Matríz de confusión de la validación',
row.names = T)| 1 | 2 | 3 | 4 | 5 | 6 | |
|---|---|---|---|---|---|---|
| 1 | 42 | 11 | 2 | 12 | 12 | 0 |
| 2 | 8 | 49 | 27 | 4 | 2 | 8 |
| 3 | 1 | 15 | 52 | 3 | 1 | 4 |
| 4 | 19 | 19 | 6 | 75 | 16 | 0 |
| 5 | 6 | 4 | 1 | 6 | 68 | 0 |
| 6 | 1 | 2 | 11 | 0 | 0 | 88 |
val.cart.acc <- ovAcc(conmat[-nrow(conmat),])
message(sprintf("Exactitud global: %.3g\nCoeficiente Kappa: %.3g",
val.cart.acc[[1]][1], val.cart.acc[[1]][2]))## Exactitud global: 0.65
## Coeficiente Kappa: 0.58colnames(val.cart.acc[[2]]) <- c("Productor", "Usuario")
knitr::kable(val.cart.acc[[2]], digits = 3, caption = 'Exactitud de productor y usuario para la validación')| Productor | Usuario |
|---|---|
| 0.545 | 0.532 |
| 0.490 | 0.500 |
| 0.525 | 0.684 |
| 0.750 | 0.556 |
| 0.687 | 0.800 |
| 0.880 | 0.863 |
Se obtuvieron valores de exactitud global de 0.65 y el coeficiente de kappa obtenido fue de 0.58. En cuanto a la exactitud global y de usuario se puede ver de mejor manera en la Figura 3.6, en donde se ve que la mayór exactitud se presenta en la cobertura 6, correspondiente a las nubes. En el caso de lacategoría 4, correspondiente con las zonas humedas y pantanosas, se destaca la gran diferencia entre la exactitud de usuario y de productor, esta categoría tiene un bajo error por omisión (0.25), sin embargo, el error de exceso si significativamente mayor (0.444). Caso contrario a lo que sucede con la categoría 3, bosques y áreas naturales, en donde el error por omisión es el mayor entre ambos.
barp(t(val.cart.acc[[2]]), col = c("#ff8008ff", "#0880ffff"), main = "Exactitud de usuario y de productor",
xlab = "Clase", ylab = "Exactitud", legend.lab = c("Productor", "Usuario"), legend.pos = c(0.5,0.8))
Figure 3.6: Exactitud de usuario y productor por categoría clasificada
4 Discusión
Como primer acercamiento a una clasificación de las coberturas de la tierra en la zona de estudio, se considera como satifactorio el presente ejercicio, en el cual se obtuvo una exactitud global de 0.65 y un coefficiente kappa de 0.58.
Durante el proceso se identificó que la clasificación de nubes fue errada, al comparar las Figuras 2.3 y 3.4, la razón de esto puede radicar en que las localización de las nubes no es estacionaria, y en la información de referencia se localizaron en una zona diferente a la que tenian el 4 de enero del 2015. Por lo cual el algoritmo erró en su clasificación, aún cuando su exactitud fue bastante buena.
Existe otro caso donde la cobertura de referencia fue significativamente diferente a la cobertura registrada por las imagenes satelitales, es el caso de la represa de hidrosogamoso, que fue llenada en el lapso de tiempo entre la información de referencia (2010-2012) y las imagenes Landsat (2015). Si bien en este caso, visualmente parece haber hecho una buena clasificación el algoritmo, es posible que los puntos allí ubicados durante la validación de la exactitud muestren un error por exceso de las “superficies de agua.”
Uno de los últimos aspectos notables fue la existencia de múltiples pixeles que fueron identificados como zonas urbanas y mineras en la parte baja de la SZH, la razón de este comportamiento es poco clara, pero se dió con una gran cantidad de pixeles aislados, por lo que es posible que una clasificación basada en objetos arrojara mejores resultados en este sentido.
5 Conclusiones
El algoritmo de randon forest mostró buenos resultados en este primer ejercicio. Sin embargo, las limitaciones y oportunidades de mejora identificadas en la sección anterior muestran el trabajo futuro que puede ser realizado para el mejoramiento de los resultados.
Las nubes en información de referencia probablemente va a traducirse en una clasificación errónea. Una posibilidad es hacer la remoción de estos poligonos de la capa de coberturas del IDEAM. Así mismo, puede considerarse hacer un pretratamiendo de las imagenes para enmascarar las zonas nubladas. En este mismo sentido, es notable que la categoría con mayor exactitud fue esta, por lo cual se realza la importancia de la verificación cualitativa y del criterio del profesional que está desarrollando la clasificación.
Finalmente, se resalta que la fecha de adquisición de los datos de referencia y de las imagenes satelitales sea lo mas cercana posible para evitar que cambios en la cobertura genere errores o confusiones en el algortimo de clasificación.
6 Reproducibilidad
Cite as: Cortés-Ramos, Diego, 2021. Clasificación de las Coberturas de la Tierra para la Subzona Hidrológica del río Sogamoso. Disponible en: https://rpubs.com/diacortesra/LandClass_Sogamoso
sessionInfo()## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] leaflet_2.0.4.1 plotrix_3.8-1 bookdown_0.21 corrplot_0.84
## [5] ranger_0.12.1 dplyr_1.0.5 raster_3.4-5 terra_1.0-10
## [9] sp_1.4-5
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.6 highr_0.8 bslib_0.2.4 compiler_4.0.3
## [5] pillar_1.5.1 jquerylib_0.1.3 tools_4.0.3 digest_0.6.27
## [9] jsonlite_1.7.2 evaluate_0.14 lifecycle_1.0.0 tibble_3.1.0
## [13] lattice_0.20-41 pkgconfig_2.0.3 rlang_0.4.10 Matrix_1.2-18
## [17] DBI_1.1.1 crosstalk_1.1.1 rgdal_1.5-23 yaml_2.2.1
## [21] xfun_0.22 stringr_1.4.0 knitr_1.31 htmlwidgets_1.5.3
## [25] generics_0.1.0 sass_0.3.1 vctrs_0.3.6 tidyselect_1.1.0
## [29] grid_4.0.3 glue_1.4.2 R6_2.5.0 fansi_0.4.2
## [33] rmarkdown_2.7 purrr_0.3.4 magrittr_2.0.1 codetools_0.2-16
## [37] htmltools_0.5.1.1 ellipsis_0.3.1 assertthat_0.2.1 utf8_1.2.1
## [41] stringi_1.5.3 crayon_1.4.1