# Instalar y cargar paquetes
# install.packages("RStoolbox")
library(RStoolbox)
# install.packages("raster")
library(raster)Clusterización de variables ambientales
Introducción
La evaluación y planificación de actuaciones de un área afectada por un incendio forestal puede requerir la integración de multitud de variables en una extensión amplia del territorio. La mera observación de múltiples capas de información puede no ser suficiente a la hora de tomar decisiones. Por ello, es habitual realizar análisis de la información que nos permitan combinarla. Esta combinación puede hacerse de distintas maneras, desde técnicas estadísticas a análisis multicriterio, dependiendo de las fuentes de información que consideremos y el objetivo planteado. En esta sesión exploraremos una de estas opciones, la creación de áreas homogéneas mediante clusterización.
La creación de mapas de áreas homogéneas en virtud de distintas variables ambientales es una técnica muy utilizada para planificar actuaciones en el territorio. Esto es especialmente relevante en zonas muy extensas que no pueden ser manejadas sólo aplicando el conocimiento experto.
La idea clave de estas áreas homogéneas es que, con su creación, disponemos de una espacialización de ciertas variables, que reflejamos en un mapa. En lugar de contemplar el mundo con toda su complejidad, lo vemos sólo con la visión de las variables que hemos definido para generar dicho mapa. Es decir, podemos “leer” el territorio usando esas variables.
Por ejemplo, en una zona quemada, en lugar de tener un gradiente continuo de multitud de variables descriptoras (capacidad de regeneración, profundidad del suelo, etc) tendremos una diferenciación en zonas con características acotadas, como podrían ser “lugares con suelos poco profundos y alta insolación”. De esta forma podemos plantear actuaciones de manejo particularizadas a cada unidad homogénea.
Flujo de trabajo
En esta práctica aplicaremos una clasificación no supervisada basada en el algoritmo k-means, del que puedes leer más aquí. Para ello trabajaremos desde R con las capas que hemos generado durante las prácticas con qGIS. El flujo de trabajo será el siguiente:
- Selección de variables
- Preparación de los datos
- Elección del número óptimo de áreas homogéneas
- Generación del mapa
- Evaluación del mapa
Selección de variables
Para la creación del mapa de áreas homogéneas en primer lugar necesitaremos definir qué variables consideramos relevantes para nuestro objetivo. Una vez definidas estas variables podremos combinarlas para generar el mapa. Es habitual que tengamos que repetir el proceso varias veces, iterando sobre distintos conjuntos de variables hasta que demos con una combinación que nos satisfaga.
¿Qué variables crees que pueden ser útiles para nuestro caso de uso?
Preparación de los datos
En este ejercicio utilizaremos el paquete RStoolbox para realizar una clasificación no supervisada. Además utilizaremos el paquete raster para trabajar con las capas ya que RStoolbox utiliza los objetos espaciales característicos de este paquete. Para otros análisis espaciales con R el paquete terra es más adecuado (más moderno, sintáxis más sencillam más intuitivo). Es importante tener cuidado al trabajar con estos dos paquetes de forma simultánea ya que existe cierto solapamiento entre funciones. Por ello si trabajamos con ambos paquetes a la vez se recomienda utilizar la sintaxis paquete::funcion cuando llamamos a una función, asegurándonos de estar usando la correcta.
En primer lugar cargaremos las capas con las que vamos a trabajar y las exploraremos. Entonces las apilaremos en un mismo stack. La función stack del paquete raster crea una pila de rásters proveniente de distintos archivos. El algoritmo k-means requiere que todas las capas “encajen” entre sí, por lo que necesitaremos apilarlas en un stack.
# Cargar máscara
mask <- raster("capas_resultantes/mascara_23030.tif")
# Cargar capas con las que vamos a trabajar: temperaturas, precipitaciones, aridez, aspecto, pendiente
tmax <- raster("capas_resultantes/tmax_1971_2000_23030_clipbyextent.tif")
tmin <- raster("capas_resultantes/tmin_1971_2000_23030_clipbyextent.tif")
prec <- raster("capas_resultantes/p_1971_2000_23030_clipbyextent.tif")
aridez <- raster("capas_resultantes/i_aridez_23030.tif")
aspecto <- raster("capas_resultantes/aspecto_23030.tif")
pendiente <- raster("capas_resultantes/pendiente_en_grados_23030.tif")
# Explorar capas
tmax
tmin
prec
aridez
aspecto
pendiente
# Apilar capas:
my_stack <- raster::stack(tmax, tmin, prec, aridez, aspecto, pendiente) Hasta ahora hemos trabajado en qGIS con capas que tenían distintos distintas resoluciones y sistemas de referencia, sin haber tenido problemas por ello. qGIS proyecta al vuelo las capas que no están en el mismo sistema de referencia que nuestro proyecto, aunque como dijimos es más adecuado reproyectar las capas para asegurarnos de que todos los cálculos se hacen correctamente. R, sin embargo, no hace casi nada “under the hood”, sino que nos avisa cuando algo da problemas.
Al explorar las capas hemos visto que éstas tienen distintas propiedades. Aunque sabemos que todas están en el mismo sistema de referencia, no todas tienen el sistema asignado y, además, no tienen la misma resolución y extensión. Esto hace que sus celdas no coincidan, por lo que necesitaremos alinearlas. Dos capas ráster están alineadas cuando están en el mismo sistema de referencia, tienen la misma resolución y el mismo origen. Es decir, que encajan en todos los aspectos. Por tanto, tendremos que asignar el sistema de coordenadas a aquellas capas que no lo tienen (OJO: ¡sólo cuando estamos seguros de ello!) y remuestrear las capas para que tengan la misma resolución y origen.
# Asignar sistema de coordenadas:
crs(aspecto) <- CRS('+init=EPSG:23030')
crs(pendiente) <- CRS('+init=EPSG:23030')
# Reproyectar (no necesario en este caso)
# capa_EPSG_Y <- projectRaster(capa_EPSG_X, mask_EPSG_Y)
# Remuestrear: generar nueva capa con la resolución deseada
tmax_resampled <- resample(tmax, mask, method = "bilinear")
tmin_resampled <- resample(tmin, mask, method = "bilinear")
prec_resampled <- resample(prec, mask, method = "bilinear")
aridez_resampled <- resample(aridez, mask, method = "bilinear")
aspecto_resampled <- resample(aspecto, mask, method = "bilinear")
pendiente_resampled <- resample(pendiente, mask, method = "bilinear")Es importante prestar atención al método utilizado en el remuestreado, ya que el método más adecuado dependerá de la naturaleza de nuestra variable. Una vez reproyectadas y remuestreadas nuestras capas no deberíamos tener problemas para apilarlas.
my_stack <- raster::stack(tmax_resampled, tmin_resampled, prec_resampled, aridez_resampled, aspecto_resampled, pendiente_resampled)
# Explorar stack
my_stack
plot(my_stack)
# Asignar nombres a capas
names(my_stack) <- c("tmax", "tmin", "prec", "aridez", "aspecto", "pendiente")Por último, ya que el objetivo es generar un mapa de áreas homogéneas de la zona afectada por el incendio no tiene sentido considerar las zonas fuera del perímetro del mismo. Por ello, antes de proseguir aplicaremos una máscara a nuestro stack.
# Enmascarar: mi máscara tiene valores 1 en la zona del incendio y no valores en la zona exterior
stack_masked <- my_stack * mask
names(stack_masked) <- c("tmax", "tmin", "prec", "aridez", "aspecto", "pendiente")
plot(stack_masked)Elección del número óptimo de áreas homogéneas
Una vez construida nuestra pila de variables tenemos que decidir cuántas áreas son adecuadas. Para ello es habitual realizar la clusterización con distinto número de áreas y evaluar los resultados basándonos en la variabilidad de las áreas resultantes. Para ello calcularemos la suma de cuadrados dentro de cada grupo
## Choose appropriate number of classes
set.seed(25) # Definir semilla para obtener resultados reproducibles en procesos aleatorios, si no ejecutamos esta línea cada vez que lancemos el código obtendremos resultados ligeramente distintos.
wss <- numeric(10) # Inicializar el error de cuadrados totales (total within sum of squares error, wss)
for (i in 1:10) { # Ejecutar la clusterización con de 1 a 10 número de áreas
out <- unsuperClass(stack_masked, nClasses=i, nStarts=i*4, nIter = 2000, norm=T, clusterMap=F, algorithm = "Hartigan-Wong")
wss[i] <- out$model$tot.withinss # Guardar el wss de cada iteración
}
plot(1:10, wss, type = "b", xlab = "Number of Clusters", ylab = "Within groups sum of squares")Basándonos en esta gráfica ¿cuál sería el número de clases óptimo?
nr_optimo <- 4Generación del mapa
Ya podemos generar el mapa de áreas homogéneas deseado ejecutando la clasificación no supervisada.
clust <- unsuperClass(stack_masked,
nClasses = nr_optimo, nStarts = nr_optimo*4,
nIter=2000, norm=T, clusterMap=F, algorithm = "Hartigan-Wong")Evaluación del mapa
Por último evaluamos el mapa obtenido y lo exportamos a un archivo ráster que podremos guardar, compartir y abrir en qGIS o cualquier otro SIG.
plot(clust$map)
# Explorar características de las zonas
clust$model
writeRaster(clust$map, filename = "capas_resultantes/areas_homogeneas_4_23030.tif")