💡 OBJETIVOS:

💡 MATERIALES PARA LA ACTIVIDAD:

Los materiales para el desarrollo de esta actividad son los siguientes.

► Imágenes satélite.

► script_laboratorio_7_clasificacion_manual

INTRODUCCIÓN

La clasificación manual de imágenes de satélite es un enfoque dentro del análisis de imágenes en el campo de la Teledetección, en el que el analista interpreta y clasifica los píxeles basándose en su experiencia, conocimiento del terreno y observación directa de las imágenes.

Los criterios utilizados para distinguir clases (agua, vegetación, suelo urbano…) son diversos:

Existen diferentes variantes, pero una de las más utilizadas es la separación de los niveles de densidad o intensidad de una imagen (sobre todo, de un determinado un índice espectral, pero también valores de reflectancia de una banda). Se basa en una idea simple: dividir los valores de esa imagen en dos o más intervalos y asignar cada intervalo a una categoría. Este procedimiento transforma un ráster continuo (pe. valores de NDVI) en un raster categórico (mapa temático). El procedimiento comprende varias acciones:

Una vez obtenidas las categorías, se puede analizar cada una de manera independiente o de manera conjunta. Por ejemplo, si se desea calcular la superficie de un embalse en distintos momentos del tiempo para comparar con el volumen de agua embalsada, el análisis resulta más sencillo si se eliminan los píxeles que no representan agua, conservando únicamente aquellos que sí lo hacen.

Es un procedimiento muy usado en las fases iniciales de cualquier investigación. Como resumen, los métodos de clasificación manual ofrecen ventajas y limitaciones:

Por esta razón, se usan

Tipología de métodos de clasificación:

Dentro de este enfoque, técnicas como la umbralización y el density slicing son métodos sencillos pero muy útiles para separar clases de información. Dependiendo del número de categorías a seleccionar, se puede diferenciar:

  • Umbralización (thresholding)

Consiste en separar píxeles según si su valor está por encima o por debajo de un umbral. Ejemplo :

  • NDVI > 0.3 → vegetación

  • NDVI ≤ 0.3 → no vegetación

Este es el caso más sencillo, y en realidad, es un ejemplo de creación de máscaras binarias. Ejemplo: el valor 0 de una imagen con el índice NDWI puede utilizarse como umbral para distinguir superficies acuáticas de superficies continentales. Los valores por encima de ese umbral se les asigna el valor de 1 (o TRUE), mientras que los inferiores pueden clasificarse como valores nulos (NA o FALSE). A partir del recuento del número de valores 1 es posible identificar masas de agua, establecer su perímetro y, en caso de disponer de imágenes obtenidas en diferentes momentos del tiempo, analizar su variación temporal.

Es un procedimiento fácil de implementar, rápido y funciona bien cuando hay contraste claro (agua vs tierra) o cuando sólo se necesita identificar una clase frente al resto. Pero, al mismo tiempo, es sensible a iluminación y ruido, con clases complejas o superpuestas.

  • Density slicing o clasificación multiumbral, en cuyo caso se establecen más de dos categorías. En este caso también existe una variante adaptativa por la que el umbral cambia de valor según la zona. En esencia, es una extensión del método anterior con la aparición de varias categorías.
Rango NDVI Clase Color
< 0 Agua Azul
0–0.2 Suelo desnudo Marrón
0.2–0.5 Vegetación baja Amarillo
> 0.5 Vegetación densa Verde

Como en el caso anterior, la clave de este procedimiento radica en la correcta elección de los valores o umbrales que separan las categorías. Su definición se puede realizar de forma manual o mediante herramientas estadísticas. Dependiendo del grado de control del proceso de elección de umbrales se pueden señalar:

  • Manual: el analista define los umbrales, por lo que es muy dependiente de la experiencia

  • Asistido: basado en histogramas (por ejemplo, cuantiles o métodos estadísticos), siendo más objetivo pero menos interpretativo.

Aspecto Umbralización Density slicing
Complejidad Baja Baja-media
Salida Binaria Multiclase
Uso Clasificación directa Interpretación visual
Flexibilidad Limitada Mayor

PREPARACIÓN DE LOS DATOS

Los paquetes básicos para realizar este ejercicio son los siguientes:

library(terra)                             # Importación y transformación de imágenes ráster
## terra 1.8.93
library(ggplot2)                           # Visualización de datos ráster con RStoolbox
library(sf)                                # Manejo de ficheros en formato vectorial.
## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
library(RStoolbox)                        # Transformación y visualización de datos ráster.
## This is version 1.0.2.2 of RStoolbox
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.2
## Warning: package 'readr' was built under R version 4.3.2
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ readr     2.1.5
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ lubridate 1.9.5     ✔ tibble    3.3.1
## ✔ purrr     1.2.1     ✔ tidyr     1.3.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ tidyr::extract() masks terra::extract()
## ✖ dplyr::filter()  masks stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

En el presente apartado trabajaremos con el objeto ndvi, que contiene los valores del índice homónimo, y el objeto imagen, con el que se ha trabajado en unidades anteriores. Para ello, se establece la carpeta de trabajo

setwd("D:/G174_2026/LABORATORIO_7_Clasificacion_imagenes_PIXEL/datos/manual/")

A continuación se importan ambos objetos ráster

imagen <- rast("D:/G174_2026/LABORATORIO_7_Clasificacion_imagenes_PIXEL/datos/manual/imagen.tif")
names(imagen) <- c("coastal aerosol", "blue", "green", "red", "NIR", "SWIR1", "SWIR2")

A continuación, se calculará el índice NDVI

ndvi <- (imagen[["NIR"]] - imagen[["red"]]) / (imagen[["NIR"]] + imagen[["red"]])
ndvi
## class       : SpatRaster 
## size        : 1245, 1497, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 594090, 639000, 4190190, 4227540  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 10N (EPSG:32610) 
## source(s)   : memory
## varname     : imagen 
## name        :       NIR 
## min value   : -0.959918 
## max value   :  0.883263

Se realiza una primera representación gráfica del índice NDVI y se compara con la imagen en falso color.

par(mfrow = c(1,2))
#imagenFCC <- imagen[[c(5,4,3)]]
plotRGB(imagen, 
        r=5, g=4, b=3, 
        stretch = "lin",
        main="Imagen en falso color", 
        mar=c(3.5, 3.5, 2.5, 2.5))

plot(ndvi, main = "Índice NDVI")

par(mfrow = c(1,1))

UMBRALIZACIÓN

Método manual

Un enfoque tradicional es utiliza el valor 0 como umbral. Esto significa que los pixeles que tengan un valor superior a 0 se clasificarán como perteneciente a una determinada clase, por ejemplo, vegetación o agua, asignando NA a los píxeles con valores inferiores.

En la función terra:: classify()[https://rspatial.github.io/terra/reference/classify.html] los argumentos principales son la imagen ráster de entrada y la matriz que sirve de guía para la reclasificación. Esta matriz contiene tres columnas.

  • La primera y la segunda columnas contienen los umbrales (límites inferior y superior de un rango).

  • La tercera columna contiene el valor que se asignará a ese rango.

A continuación se muestra un ejemplo:

umbral_mascara1 <- c(-Inf, 0, NA,  
                 0, Inf, 1)

matriz_mascara1 <- matrix(umbral_mascara1, 
                           ncol=3, 
                           byrow=TRUE)
matriz_mascara1
##      [,1] [,2] [,3]
## [1,] -Inf    0   NA
## [2,]    0  Inf    1

A continuación se clasifica el objeto ndvi en las categorías señaladas, y se representa gráficamente.

mascara1 <- classify(ndvi, 
                     matriz_mascara1)
plot(mascara1)

📝 ACTIVIDAD DE EVALUACIÓN 1:

A continuación se realizará un ejercicio en el que se va a utilizar el objeto ráster ndvi generado anteriormente. Si consideramos que los píxeles con valores superiores a 0.6 probablemente correspondan a bosque cerrado, utilícese este umbral para determinar si existe este tipo de cobertura del suelo.

► Crea la matriz con las clases por debajo y por encima.

► Clasifica el ráster según los umbrales definidos

Hasta este momento se ha clasificado el objeto ráster en 2 categorías, por encima y por debajo de un determinado valor. Sin embargo, también se pueden obtener categorías intermedias. Por ejemplo, la literatura asume que valores entre 0.2 y 0.4 corresponden a pastizales o áreas forestales degradadas. Entonces, se asignará el valor ‘1’ a todos los píxeles dentro de ese rango de valores y NA a los píxeles por encima y por debajo de esos límites

umbral_pasto <- c(-Inf, 0, NA,
                0, 0.3, 1, 
                0.3, Inf, NA)

matriz_pasto <- matrix(umbral_pasto, 
                      ncol=3, 
                      byrow=TRUE)

mascara_pasto <- classify(ndvi, 
                         matriz_pasto)
plot(mascara_pasto, col = "#d5d685")

Esta última imagen puede superponerse a una imagen en color natural para mejorar la interpretación del tipo de superficie

plotRGB(imagen, r=4, g=3, b=2, axes=TRUE, stretch="lin")
plot(mascara_pasto, 
     main = 'Mapa pasto', 
     col = "#d5d685", add=TRUE, legend=FALSE)

Identificación de umbrales usando un histograma de frecuencias.

Una manera de definir uno (o varios) umbral(es) cuando se carece de valores predeterminados es explorar la distribución de esos valores en el objeto ráster mediante la función hist(), que produce un histograma a partir del cual se pueden identificar tanto valores erróneos, valores extremos pero también cambios en la distribución de frecuencias que pudieran sugerir discontinuidades (umbrales).

hist(ndvi,
     main = "Distribución de valores del NDVI según píxeles",
     xlab = "NDVI",
     ylab= "Frecuencia",
     col = "green",
     xlim = c(-0.5, 1),
     breaks = 30,
     xaxt = 'n')
## Warning: [hist] a sample of 54% of the cells was used
axis(side=1, at = seq(-0.5,1, 0.05), labels = seq(-0.5,1, 0.05))

En este caso, se aprecia una primera discuntinuidad entorno a -0.30. A su vez, a partir de -0.05 se observa otra discontinuidad, que podría corrresponder a la aparición de agua. Obsérvese que también hay un corte que coincide con el valor 0.15, que podría representar el tránsito entre superficies no vegetadas y las vegetadas. Por la derecha de la distribución la frecuencia de valores cae rápidamente por encima de 0.35, que podría representar la aparición de superficies arboladas.

Selección del umbral mediante el método de Otsu

El procedimiento manual puede ser refinado utilizando el método de discretización de clases de Otsu, llamado así en honor a Nobuyuki Otsu en 1979 (Otsu, 1979), utiliza técnicas estadísticas, para resolver el problema. En su forma clásica está diseñado para separar en dos clases (es decir, una sola frontera o umbral):

o Clase 1: píxeles por debajo del umbral

o Clase 2: píxeles por encima del umbral

El objetivo es encontrar el valor óptimo del umbral que maximice la separación entre ambas clases (maximizando la varianza entre clases, una medida de la dispersión de valores). Su método calcula el valor umbral de forma que la dispersión dentro de cada segmento sea lo más pequeña posible, pero al mismo tiempo la dispersión sea lo más alta posible entre segmentos diferentes. Para ello se calcula el cociente entre ambas varianzas y se busca un valor umbral para el que este cociente sea máximo. Es, por lo tanto, muy útil para clasificaciones dicotómicas, como separar agua vs tierra, o vegetación vs no vegetación (en índices como NDVI). Existe una generalización llamada Otsu multiclase que busca varios umbrales, pero que no es muy utilizado ya que el coste computacional crece mucho al aumentar el número de clases; por ello, a menudo se prefieren métodos más avanzados (k-means, clasificación supervisada, etc.)

A continuación se mostrará una aplicación práctica de esta metodología, aplicándose a un raster en el que se ha calculado el índice NDWI: los valores negativos corresponden a terreno continental y los positivos a agua.

\[NDWI = (G - NIR)/(G + NIR)\]

El primer procedimiento consiste en calcular el citado índice

ndwi <- (imagen[["green"]] - imagen[["NIR"]]) / (imagen[["green"]] + imagen[["NIR"]])
plot(ndwi, col=rev(topo.colors(100)), main="NDWI")

Para destacar y extraer los cuerpos de agua de mayor tamaño utilizando el método de Otsu se recurre a una función del paquete EBImage. EBImage proporciona una funcionalidad de propósito general para el procesamiento y análisis de imágenes (Pau et al. [2010]) y es un paquete creado y mantenido por la iniciativa Bioconductor. Para la mantención e instalación de sus propios paquetes instalamos en primer lugar el paquete para gestión, BiocManager (Morgan [2021]).

#if (!require("BiocManager", quietly = TRUE))
#install.packages("BiocManager")

# BiocManager::install("EBImage", force = TRUE)

Y luego cargar en la sesión actual:

library(EBImage)
## 
## Attaching package: 'EBImage'
## The following object is masked from 'package:purrr':
## 
##     transpose
## The following objects are masked from 'package:terra':
## 
##     flip, rotate, thresh, watershed

Las funciones del paquete requieren transformación del objeto ráster a un formato numérico (array o matriz). La función terra::as.array() [https://rspatial.github.io/terra/reference/coerce.html] que toma un raster y devuelve como array.

array <- as.array(ndwi)
str(array)
##  num [1:1245, 1:1497, 1] -0.465 -0.469 -0.497 -0.545 -0.564 ...

La función EBIimage::otsu() devuelve un valor umbral o de corte para binarizar la imagen

umbral <- otsu(array, 
              range = c(-1,1))                              # Número de niveles de gris
umbral
## [1] -0.08984375

Aplicamos el valor de corte, pero antes duplicaremos la imagen original para trabajar sobre ella.

agua <- ndwi                                           
agua[agua >= umbral] <- 1
agua[agua < umbral] <- NA

El resultado del proceso anterior es una imagen binaria, en el que el valor 1 para los píxeles asignados como cuerpo de agua y el resto NA.

plot(agua, 
     col= "blue", 
     axes=F,
     box=T, 
     legend=F)

Esta imagen se puede comparar con la proporcionada por el ráster con el índice original

par(mfrow=c(1,2))
plot(ndwi, col = rev(viridis::viridis(100)), main="NDWI")
plot(agua, main="Superficies acuáticas (Otsu)", col = "blue")

par(mfrow=c(1,1))

Una vez creada la máscara, se podría aplicar al fichero original

ndvi_enmascarado <- mask(ndvi, agua)
plot(ndvi_enmascarado, col = "blue")

El método de Otsu es adecuado para la clasificación de índices espectrales sencillos, como NDVI (para separar las superficies con vegetación de las que no lo tienen), NDBI (extraer las zonas construidas) o el NDWI (extraer las superficies con agua). Como ventajas adicionales se encuentra el que es un procedimiento automático, sin intervención de la subjetividad del usuario, que funciona especialmente bien cuando los valores de esos índices muestran una distribución bimodal (histograma con dos clases dominantes), pero, por contra, no suele funcionar bien si la distribución es unimodal o presenta diversas clases dominantes.

📝 ACTIVIDAD DE EVALUACIÓN 2:

El fichero sentinel_santander.tif es una imagen multibanda correspondiente a Santanader, tomada por el satélite Sentinel 2.

Carga de las imágenes y representación gráfica

► Carga esta imagen en R y represéntalas gráficamente.

Las bandas son las siguientes

Banda Nombre
B1 Coastal aerosol
B2 Blue
B3 Green
B4 Red
B5 Vegetation red edge 1
B6 Vegetation red edge 2
B7 Vegetation red edge 3
B11 SWIR1
B12 SWIR2
B8A Narrow NIR
B8 NIR

► Cambia los nombres de las bandas y ordénalas

Cálculo del índice NDVI

► Calcular el índice NDVI, teniendo en cuenta que para Sentinel 2 la fórmula es la siguiente:

\[NDVI = \frac{B8 - B4}{B8 + B4}\]

► Combina en el mismo gráfico una imagen en color verdadero y la imagen del índice NDVI.

► A qué valores crees que corresponden:

  • Marismas (Bahía de Santander)

  • El agua (Bahía de Santander)

  • Zonas urbanas.

  • Zonas urbanas.

  • Vegetación.

Definición de un umbral para la reclasificación.

► Elabora un histograma de frecuencias con los valores del índice NDVI. ¿Proporciona alguna indicación útil?

► Tras inspeccionar el histograma, propon un umbral para identificar la zona de marismas.

Reclasificación y representación gráfica de las superficies seleccionadas

► A continuación se clasifica el objeto ndvi en las categorías señaladas, y se representa gráficamente.

► Representa gráficamente la zona de marismas superpuesta a la imagen en color natural

► A continuación, representa gráficamente las superficies acuáticas profundas.

► Representa gráficamente la zona de marismas superpuesta a la imagen en color natural

Aplicación del método de Otsu

► Dado que existen dificultades para identificar los límites de las superficies acuátivas de la bahía de Santander, se procede al cálculo y representación del índice NDWI.

► Aplica el procedimiento de Otsu para ver si es posible obtener un valor que actúe de umbral.

► Representa gráficamente la zona de marismas superpuesta a la imagen en color natural

► Compara esta imagen con la del índice NDWI original

📝 ACTIVIDAD DE EVALUACIÓN 3:

El monitoreo de glaciares mediante imágenes de satélite es una herramienta fundamental en la glaciología moderna, ya que permite observar de forma continua y sistemática zonas que suelen ser inaccesibles o peligrosas para el trabajo de campo. Gracias a sensores como los de Landsat 8, es posible analizar grandes extensiones de terreno sin necesidad de presencia física.

Una de las principales ventajas de este enfoque es la capacidad de seguimiento temporal. Las series de imágenes satelitales, disponibles desde hace décadas, permiten estudiar la evolución de los glaciares a largo plazo, observando tendencias como el retroceso frontal, la reducción de superficie o los cambios en la dinámica del hielo.

Además, el uso de técnicas de teledetección facilita la cartografía automática de nieve y hielo mediante índices espectrales como el NDSI o el NDVI. Estas herramientas permiten distinguir entre nieve, hielo, vegetación y roca, generando mapas temáticos que representan de forma clara la distribución de las distintas cubiertas del suelo.

El monitoreo satelital también permite estimar variables clave como la altitud de la línea de nieve, el área glaciar o incluso aproximaciones al balance de masa. La combinación de imágenes multiespectrales con modelos digitales del terreno mejora significativamente la precisión de estos análisis.

A escala espacial, los satélites ofrecen una cobertura regional y global, lo que permite comparar glaciares entre distintas cordilleras y construir inventarios completos a nivel continental o mundial. Esto supone una ventaja clara frente a los métodos tradicionales de observación local.

Sin embargo, este tipo de análisis presenta algunas limitaciones, como la interferencia de nubes, la confusión entre nieve y nubes brillantes o las dificultades para detectar hielo cubierto por detritos. Por ello, es habitual combinar distintos índices espectrales, modelos digitales del terreno y validaciones en campo para mejorar la fiabilidad de los resultados.

En conjunto, el monitoreo satelital de glaciares se ha convertido en una herramienta esencial para la investigación del cambio climático, la gestión de recursos hídricos y la evaluación de riesgos naturales en zonas de montaña, proporcionando información clave para comprender la evolución del sistema criosférico a escala global.

A continuación se desarrolla un ejercicio consistente en la identificación de los glaciares de la Cordillera Huascarán, en los Andes de Perú.

Definir el entorno

► Definir el directorio de trabajo con los datos.

► Cargar los archivos individuales en un único objeto ráster.

► Los valores de reflectancia superfial de las bandas aparecen como números digitales. Deben transformarse en valores reales de reflectancia superficial

► La reflectancia superficial oscila entre 0 y 1, pero algunas bandas poseen valores inferiores a 0 y superiores a 1. Pueden ser considerados resultado de errores y omitirse. Para ello, se deben ajustar los valores máximos y mínimos de todas las bandas a un rango entre 0 y 1 con la función terra::clamp.

Representación gráfica

► Representación gráfica de todas las bandas.

► Representa el área de trabajo combinando en el mismo gráfico una combinación de bandas en color natural y otra en falso color (RGB = SWIR1 – NIR – Red).

¿A qué colores en falso color corresponden las siguientes superficies:

  • Nieve/Hielo glaciar.

  • Terreno sin vegetación.

  • Terreno con vegetación.

  • Agua (lagos)

Obtención del índice NDSI

El índice más utilizado para la identificación y delimitación de superficies con nieve y hielo glaciar es el NDSI (Normalized Difference Snow Index). Su fórmula es la siguiente:

\[NDSI = \frac{Green - SWIR1}{Green + SWIR1}\]

► Calcula el índice NDSI para la zona de estudio.

Puede habar algún problema con divisiones en la que ambas bandas tienen un valor de 0. Para ello aplica la siguiente fórmula `ndsi[(green + swir1)] == 0] <- NA

► Si es necesario, cambia el nombre de la capa a ndsi

La literatura científica considera los siguientes umbrales

Clase NDSI
Nieve >0.4
Hielo glaciar limpio >0.5
Roca <0
Vegetación bajo

► Analiza un histograma con la frecuencia de valores del índice NDSI para determinar si los umbrales propuestos son aplicables el ejemplo señalado.

► Aplica el procedimiento de Otsu para confirmar estos resultados

► Aplica el valor recomendado de 0.4 en la literatura científica.

► Representa gráficamente la imagen, en el que el valor 1 para los píxeles asignados como cuerpo de agua y el resto NA.

► Esta última imagen puede superponerse a una imagen en color natural para mejorar la interpretación del tipo de superficie

Cálculo de la altitud del campo glaciar

► Importamos el fichero con el modelo digital del terreno

► Calculo de una máscara para aislar los píxeles del MDT dónde hay nieve

► A partir de esa máscara, calcular altura mínima y media a la que aparece la nieve.

► Representa con un histograma la distribución de altitudes con nieve.

EXTRACCIÓN DE VARIAS CATEGORÍAS.

Un proceso algo más complejo es la definición de sucesivos umbrales, que conducen a la aparición de más clases, técnica que en el ámbito del análisis de imágenes se conoce con el nombre de Density slicing. Este procedimiento divide los valores de los píxeles en diferentes rangos y asigna un valor (y se imprime también un color específico) a cada rango en el ráster de salida.

Mediante el “Density slicing” es posible generar diferentes categorías (que se interpretan como diferentes coberturas o usos del territorio) extendiendo el procedimiento anterior. Primero se crea la matriz de clases aprovechando los valores de corte obtenidos de la visualización del histograma.

ndvi_clases <- matrix(c(-Inf,  -0.3, 1,  
                        -0.3,  -0.05, 2,  
                        -0.05,  0.15, 3,  
                         0.15,  0.35, 4,  
                         0.35,  Inf, 5), 
                        ncol = 3, 
                       byrow = TRUE)

A continuación se clasifica el objeto NDVI y se representa gráficamente

ndvi_clasificado <- classify(ndvi, 
                            rcl = ndvi_clases)

plot(ndvi_clasificado, 
     main = "NDVI clasificado")

El objeto ráster obtenido a transformado una variable continua (el índice ndvi) en una variable discreta (las 5 categorías de vegetación). Este ráster con valores discretos es un ráster numérico, que puede ser transformado en un ráster categórico. Para ello, primero es necesario definir los nombres de las categorías mediante un dataframe.

levels <- data.frame(value = 1:5,
                     label = c("Manglares", "Agua profunda", "Urbano", "Suelo desnudo", "Vegetación"))

A continuación se aplicarán estos niveles al objeto ndvi_clasificado.

ndvi_factor <- as.factor(ndvi_clasificado)
levels(ndvi_factor) <- levels

Finalmente esto se puede representar gráficamente.

colores_ndvi<- c("#2604b1", "#15b60c",  "#b960cb", "#e92010", "#c1eb09") 
plot(ndvi_factor, 
     col = colores_ndvi, 
     main = "Usos según NDVI")

Igualmente, se puede usar la función terra::segregate() se puede descomponer el objeto ráster en sus diferentes categorías y representarlas gráficamente

ndvi_factor_segun_categorias <- segregate(ndvi_factor, 
                                          keep = TRUE, 
                                          other = NA)
my_colors <- c("#2604b1", "#15b60c",    "#b960cb", "#e92010", "#c1eb09")
par(mfrow=c(1, 5))  # set plotting layout
for (i in 1:nlyr(ndvi_factor_segun_categorias)) {
  plot(ndvi_factor_segun_categorias[[i]], col=my_colors[i], main=paste("Layer", i))
}

En la siguiente imagen se superponen todas las categorías otra vez

plot(ndvi_factor_segun_categorias[[1]], col=my_colors[1], legend=FALSE, main="Overlay of Categories")
for (i in 2:nlyr(ndvi_factor_segun_categorias)) {
  plot(ndvi_factor_segun_categorias[[i]], col=my_colors[i], legend=FALSE, add=TRUE)
}

A continuación se eliminarán del *Global Environment` la mayoría de los objetos creados anteriormente.

rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1", "imagen_veg2", "levels", "matriz_clas_suelo", "matriz_clas_veg1", "matriz_clas_veg2", "my_colors", "ndvi_clases", "ndvi_clasificado", "ndvi_factor", "ndvi_factor_segun_categorias", "umbral_suelo", "umbral_veg1", "umbral_veg2")
## Warning in rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1",
## "imagen_veg2", : objeto 'imagen_suelo' no encontrado
## Warning in rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1",
## "imagen_veg2", : objeto 'imagen_veg1' no encontrado
## Warning in rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1",
## "imagen_veg2", : objeto 'imagen_veg2' no encontrado
## Warning in rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1",
## "imagen_veg2", : objeto 'matriz_clas_suelo' no encontrado
## Warning in rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1",
## "imagen_veg2", : objeto 'matriz_clas_veg1' no encontrado
## Warning in rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1",
## "imagen_veg2", : objeto 'matriz_clas_veg2' no encontrado
## Warning in rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1",
## "imagen_veg2", : objeto 'umbral_suelo' no encontrado
## Warning in rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1",
## "imagen_veg2", : objeto 'umbral_veg1' no encontrado
## Warning in rm("colores_ndvi", "i", "imagen_suelo", "imagen_veg1",
## "imagen_veg2", : objeto 'umbral_veg2' no encontrado

📝 ACTIVIDAD DE EVALUACIÓN 4:

El 15 de agosto de 2023 se produjo un incendio forestal en la isla de Tenerife que es considerado de «sexta generación» debido a su gran intensidad y poder destructivo; en pocos días se quemaron casi 14000 has, convirtiéndole en el peor incendio que ha sufrido Canarias en los últimos 40 años.

El índice \(\Delta NBR\) (dNBR; Differenced Normalized Burn Ratio) es un índice utilizado para evaluar la severidad de incendios forestales comparando el estado de la vegetación antes y después del fuego. Mide cuánto ha cambiado la vegetación debido al fuego, siendo por ello uno de los indicadores más usados para obtener el perímetro de un incendios, clasificar y cartografiar la severidad de un incendio y evaluar sus impactos ecológicos.

Se basa en el índice Normalized Burn Ratio, que usa:

\[ \mathrm{NBR} = \frac{NIR - SWIR}{NIR + SWIR}\]

  • NIR (infrarrojo cercano) → sensible a vegetación sana

  • SWIR (infrarrojo de onda corta) → sensible a humedad y material quemado

El procedimiento es muy sencillo. Primero se calcula el NBR en dos momentos, uno antes del incendio y otro después del incendio. Posteriormente, a partir de ambas imágenes se calcula el

\[\Delta NBR = NBR_{pre} - NBR_{post}\]

como resta del valor, pixel a pixel, de la imagen anterior y posterior.

Su interpretación es sencilla

  • Valores altos positivos → mayor severidad (vegetación muy dañada).

  • Valores cercanos a 0 → poco o ningún cambio.

  • Valores negativos → regeneración o aumento de vegetación.

► Importa los ficheros ráster necesarios para el cálculo del índice NBR

► Calcula los valores del índice antes del incendio

Y después del incendio.

► A partir de los valores previos y posteriores del índice NBR, calcula y representa gráficamente los valores de deltaNBR.

dNBR Clase Colores
< -0.1 Regeneración alta darkgreen
-0.1 – 0.1 Sin cambio gray
0.1 – 0.27 Baja severidad yellow
0.27 – 0.44 Moderada-baja orange
0.44 – 0.66 Moderada-alta red
> 0.66 Alta severidad darkred

► Elabora la matriz de reclasificación siguiendo los valores del cuadro anterior

► Aplica el procedimiento de reclasificación

► Asigna los colores y su significado al ráster

► Asigna las etiquetas al ráster y cartografía el incendio.

► Cartografía la severidad del incendio.

► Crear y aplicar una máscara seleccionando sólo las categorías 4, 5 y 6 (el resto será NA)

► Representa gráficamente este objeto ráster.