1 INTRODUCCIÓN

Los incendios forestales son uno de los principales factores de perturbación en la mayoría de las zonas de vegetación de todo el mundo. Constituyen un desafío para la gestión de los ecosistemas, porque tienen el potencial de ser a la vez beneficiosos y perjudiciales. Por un lado, los incendios forestales son un proceso natural en algunos ecosistemas, que mantienen la salud y la biodiversidad, regulando la sucesión natural de especies y la acumulación de combustible, lo que a su vez controla la edad, estructura y composición de la vegetación. Por la misma razón, influye en los ciclos de nutrientes y de energía, regulando la productividad biótica, y la diversidad y estabilidad de la vida silvestre. Por otro lado, los incendios forestales constituyen una amenaza para la propiedad, las actividades económicas y hasta la vida humana.

A pesar de la importancia de los incendios, las estimaciones actuales del impacto de los incendios de vegetación a nivel mundial siguen siendo un desafío. Se estima que varios cientos de millones de hectáreas de bosques y otros tipos de vegetación se queman anualmente en todo el mundo, consumiendo varios miles de millones de toneladas de materia seca y liberando compuestos de emisión que afectan la composición y el funcionamiento de la atmósfera global y la salud humana. Según la FAO (FAO 2012), los incendios forestales son importantes factores de forzamiento del clima, ya que liberan aerosoles entre el 25 y el 35 % de las emisiones netas totales de CO2 a la atmósfera. Durante la última década en Canadá, los incendios forestales han consumido un promedio de 1,9 millones de ha/año y han inducido costos de supresión de incendios que van desde alrededor de $ 500 millones a $ 1 mil millones al año (Canadian Forest Service, 2012). En Europa, los incendios forestales queman más de medio millón de hectáreas de áreas boscosas cada año. Más del 95% de las áreas quemadas se encuentran en la región mediterránea, en la que se han producido incendios críticos en los últimos años (http://effis.jrc.it).

Debido a la amenaza que representan los incendios, se han desarrollado sistemas operativos para su uso en el manejo de incendios que incluyen la predicción del peligro de incendios, la detección y el control de incendios. Dados los aumentos esperados en los incendios en todo el mundo debido a los cambios climáticos, una mejor predicción del peligro de incendios y la detección de incendios tendrán beneficios significativos tanto desde el punto de vista económico como desde el punto de vista de la seguridad humana en todo el mundo. También existe la necesidad de una evaluación precisa de las áreas quemadas porque están relacionadas con las emisiones de gases de efecto invernadero a la atmósfera que deben tenerse en cuenta siguiendo los requisitos del protocolo de Kioto, así como para gestionar los impactos ambientales posteriores al incendio, como la regeneración y la erosión. Las imágenes de sensores remotos transportadas por el espacio pueden desempeñar un papel importante en estos sistemas. De hecho, las imágenes satelitales ofrecen las ventajas de una amplia cobertura regional, cero perturbaciones del área a ser vista, así como un método para adquirir datos en áreas menos accesibles de manera regular y rentable.

Otro campo en el que el uso de la teledetección es habitual hace referencia al análisis y cartografía de los incendios forestales. Un fuego activo se detecta con relativa facilidad porque genera calor, que detectan los sensores térmicos, y crea humo, visible en las imágenes de satélite. La detección de incendios activos con satélites está tan desarrollada que existen diferentes productos referentes a “fuegos activos” (http://modis-fire.umd.edu/index.php), registrados cada vez que un sensor MODIS o VIIRS pasa sobre un área.

Figure 1. Ejemplo de detección activa de incendios desde el espacio durante los incendios forestales australianos de 2019. Imagen satelital del humo de los incendios forestales sobre el este de Australia por la Agencia Espacial Europea, Flickr, CC BY-SA 2.0.

Una vez que se apaga un incendio, es importante evaluar cuánta área se quema y cuál ha sido la intensidad de los fuegos. Ambas circunstancias tienen importantes implicaciones para la recuperación de la vegetación y la biota del área quemada, y también puede ayudar avanzar en el análisis de los impactos del cambio climático antropogénico, como si la frecuencia, duración, tamaño e intensidad de los incendios forestales están cambiando a medida que la Tierra incrementa su temperatura. Si bien hay diferentes formas de evaluar el impacto de los incendios a través de imágenes de satélite, el enfoque más popular implica el uso de un índice denominado Normalized Burn Ratio (NBR), que no es el único índice ni necesariamente el mejor (http://giscenter.isu.edu/research/techpg/nasa_tlcc/PDF/Ch6.pdf). El NBR es calculado como:

\[NDVI= \frac{NIR-SWIR}{NIR+SWIR}\]

\[NBRI (Landsat 8) = \frac{B5 – B7} {B5 + B7}\]

\[NBRI (Landsat 4-7) = \frac{B4 – B7} {B4 + B7}\]

Para cuantificar la gravedad de un incendio, es necesario comparar los valores de NBR en dos imágenes, una de antes y otra de después del incendio, lo que lleva al cálculo de δNBR:

\[δNBR = NBR_{antes} – NBR_{despues}\]

Para minimizar los efectos de otros factores ambientales, es ideal utilizar fechas próximas entre sí para calcular δNBR.

Figure 2. Ejemplo de una imagen previa al incendio (arriba a la izquierda), una imagen posterior al incendio (arriba en el centro), un NBR previo al incendio (arriba a la derecha), un NBR posterior al incendio (abajo a la izquierda) y δNBR (abajo en el centro). Si bien el impacto del incendio es evidente en las imágenes originales, el uso de δNBR ayuda a resaltar la diferencia entre las dos imágenes.. Example of normalized burn ratio from a fire at Camp Gurnsey, Wyoming in 2006 by Jason Karl, The Landscape Toolbox, CC0 1.0.

Debido a que NBR se calcula como una diferencia normalizada, sus valores oscilan entre -1 y +1, por lo que δNBR oscila, en teoría, de -2 a +2, con valores más altos que indican una mayor gravedad del fuego. Sin embargo, al igual que con otros índices de vegetación, los valores extremos son raros, por lo que incluso los valores de 0.2-0.3 indican impactos de fuego relativamente severos.

2 DATOS

En este ejercicio se identificarán y categorizarán las áreas quemadas después de un incendio forestal. Se utilizaron imágenes landsat 8, a partir de los que se calcularán 2 índices de vegetación para detectar las áreas quemadas: el índice NDVI y el índice NBR.

3 IDENTIFICACIÓN DE SUPERFICIES QUEMADAS MEDIANTE EL INDICE NBR (Normalized burn ratio)

La manera más simple para delimitar el área quemada y evaluar la severidad del incendio es calcular la imagen del índice NBR y asociar valores bajos del índice con alta severidad.

El índice NBR es el resultado de normalizar la diferencia de reflectancias entre una banda del infrarrojo (IR) cercano y una del infrarrojo medio. El contraste entre las reflectancias de estas bandas está relacionado con la presencia de vegetación viva. La vegetación viva es muy reflectiva en la banda del IR cercano y no tanto en la del IR medio, mientras que la que ha sufrido daños tiene cada vez valores de reflectividad más bajos en la banda del IR cercano y más altos en la del IR medio.

\[NBRi= \frac{NIR – SWIR}{NIR+ SWIR}\]

Una manera más elaborada es comparar el NBR de dos imágenes: una anterior y otra posterior al incendio. De esta manera limitamos la perturbación que ejercen superficies no quemadas pero con un bajo NBR (suelos desnudos, caminos y ciertas superficies brillantes), resaltando las superficies que han cambiado tras el incendio. A este índice se le denomina dNBR (dNBR = NBRpre-NBRpos) y varía aproximadamente entre -0.5 y +1.3 (donde los valores más negativos indican vegetación regenerada, los valores entre -1 y +0.99 zonas no quemadas y los valores más positivos zonas afectadas, aumentando con las severidad de los daños).

3.1 Cálculo del índice NBR

Para iniciar la práctica se renombrará a los ficheros que sea necesario (en este caso la banda SWIR2, la banda NIR ya ha sido renombrada antes).

swir2_pre <- pre_st_cr$pre_st_cr.3
swir2_post <- post_st_cr$post_st_cr.3

El índice NBR se desarrolló originalmente para su uso con las bandas 4 y 7 de Landsat TM y ETM+, pero funciona con cualquier sensor multiespectral con una banda NIR entre 760 - 900 nm y una banda SWIR entre 2080 - 2350 nm. Por tanto, para Landsat 8, la banda utilizada es la SWIR2.

\[NBRi (Landsat 8) = \frac{B5 – B7} {B5 + B7}\]

A continuación se calcula el índice en la imagen previa y en la imagen posterior al incendio.

nbr_pre <- (nir_pre - swir2_pre)/(nir_pre + swir2_pre) 
nbr_post <- (nir_post - swir2_post)/(nir_post + swir2_post)

Se puede exportar el fichero con el índice NBR eh formato ráster

writeRaster(nbr_pre, datatype="FLT4S", filename = "./nbr_pre.tif", format = "GTiff", overwrite=TRUE)
writeRaster(nbr_post, datatype="FLT4S", filename = "./nbr_post.tif", format = "GTiff", overwrite=TRUE)

3.2 Visualización del índice NBR en las imágenes previas y posteriores al incendio.

El primer paso es definir apropiadamente la paleta de colores (para más información debe consultarse http://colorbrewer2.org/)

nbr_colors <- colorRampPalette(brewer.pal(11, "RdYlGn"))(100)

Representación gráfica de las dos imágenes

par(mfrow=c(1,2))

plot(nbr_pre,
     main = "Índice NBR pre-incendio",
     axes = FALSE,
     box = FALSE,
     col = nbr_colors,
     zlim = c(-1, 1))

plot(nbr_post,
     main = "Índice NBR post-incendio",
     axes = FALSE,
     box = FALSE,
     col = nbr_colors,
     zlim = c(-1, 1))

Figure 4. Índice NBR previos y posteriores a un incendio

3.3 Obtención del índice dNBR (Difference Normalized Burn Index)

Como se ha señalado previamente, la extensión y gravedad de un incendio se evidencia mediante el cálculo de la diferencia del índice NBR antes y después del incendio.

\[\Delta NBR=(NBR_{prefire}−NBR_{fire})\]

ΔNBR <- (nbr_pre) - (nbr_post)

3.4 Reclasificación del valor de dNBR según intervalos de severidad

Existen varios criterios para clasificar la severidad a partir del índice ΔNBR. A continuación se expone uno de ellos:

Nivel de severidad Rango ΔNBR (reescalado \(10^3\))
Recrecimiento vigoroso -500 a -251
Recrecimiento moderado -250 a -101
No quemado -100 a 99
Severidad baja 100 a 269
Severidad baja-moderada 270 a 439
Severidad moderada-alta 440 a 659
Alta severidad 660 a 1300

Los valores originales pueden ser multiplicados por 1000 para una mejor apreciación de las diferencias.

dNBR_scaled <- 1000*dNBR

De acuerdo con esa tabla, se establecerán los rangos utilizados para clasificar el dNBR en diferentes categorías.

NBR_ranges <- c(-Inf, -500, -1, -500, -251, 1, -251, -101, 2, -101, 99, 3, 99, 269, 4, 
                269, 439, 5, 439, 659, 6, 659, 1300, 7, 1300, +Inf, -1)

A partir de esos rangos se crea una matriz de clasificación

class.matrix <- matrix(NBR_ranges, ncol = 3, byrow = TRUE)

Esa matriz de clasificación es usada para reclasificar los valores del índice dNBR en las distintas categorías.

dNBR_reclass <- reclassify(dNBR_scaled, class.matrix, right=NA)

3.5 Visualización del dNBR

Esta fase supone crear la leyenda que se utilizará en el mapa, para lo cual es necesario construir la tabla de atributos de la leyenda.

dNBR_reclass <- ratify(dNBR_reclass)                       # Establece los valores del raster como una variable categórica
rat <- levels(dNBR_reclass)[[1]]                           # Crea una rat (tabla de atributos raster)

A continuación se crea el texto que aparecerá en la leyenda.

rat$legend  <- c("NA", "Recrecimiento vigoroso", "Recrecimiento moderado", "No quemado", 
                 "Severidad baja", "Severidad baja-moderada", "Severidad moderada-alta", "Severidad extrema") 
levels(dNBR_reclass) <- rat

Los colores correspondientes a cada categoría serán los siguientes:

mis_colores <- c("white", "darkolivegreen","darkolivegreen4","limegreen", "yellow2", "orange2", "red", "purple")

El mapa de severidad contendrá un título con dos líneas: en la priemra “Severidad de incendios” y en la segunda “Kalamos, Grecia”.

plot(dNBR_reclass,
     col=mis_colores,
     legend=F,
     box=F,
     axes=F, 
     main="Severidad de incendios \n Kalamos, Greece")

La leyenda se situará al lado derecho del mapa

legend(x='right', 
       legend = rat$legend, 
       fill = my_col, y='right')

Figure 6. Mapa de categorías de severidad de un incendio

Este mapa también puede exportarse como fichero raster:

writeRaster(dNBR_reclass, datatype="FLT4S", filename = "./dNBR_reclass.tif", format = "GTiff", overwrite=TRUE)

3.6 Análisis estadístico de la distribución de categorías

Además de un mapa, se puede analizar la distribución de los valores de dNBR según categorías

barplot(dNBR_reclass,
        main = "Distribución de las categorías de severidad de incendios",
        col = my_col,
        names.arg = c("A", "B", "C", "D", "E", "F", "G", "H"), 
        las=0)

legend(x='right', y='top', legend =rat$legend, fill = my_col)

Figure 7. Distribución de los valores de los píxeles pertenecientes a cada categoría. Niveles de severidad: A = NA, B = Recrecimiento vigoroso, C = Recrecimiento moderado, D = No quemado, E = Severidad baja, F = Severidad moderada-baja, G = Severidad moderada-alta, H = Severidad alta”

Al gráfico se le puede añadir las etiquetas con el número de píxeles de cada categoría.

vals <- getValues(dNBR_reclass)

pix_na <- length(subset(vals, vals == -1))
pix_1 <- length(subset(vals, vals == 1))
pix_2 <- length(subset(vals, vals == 2))
pix_3 <- length(subset(vals, vals == 3))
pix_4 <- length(subset(vals, vals == 4))
pix_5 <- length(subset(vals, vals == 5))
pix_6 <- length(subset(vals, vals == 6))
pix_7 <- length(subset(vals, vals == 7))

pix_values <- c(pix_na,pix_1,pix_2,pix_3,pix_4,pix_5,pix_6,pix_7)

bp <- barplot(dNBR_reclass,
        main = "Número de píxeles según categoría de severidad",
        col = my_col,
        names.arg = c("A", "B", "C", "D", "E", "F", "G", "H"), 
        las=0)

legend(x='right', y='top', legend =rat$legend, fill = my_col)

text(bp, y=2500, labels=pix_values)                          # El valor de y debe ser seleccionado manualmente

Figure 8. Distribución de valores de píxeles de cada categoría para el AOI seleccionado con número de píxeles por clase

También en forma de porcentajes sobre el total

pix_total <- length(vals)

area_na_perc <- round(pix_na/pix_total*100, digits = 5)
area_1_perc <- round(pix_1/pix_total*100, digits = 5)
area_2_perc <- round(pix_2/pix_total*100, digits = 5)
area_3_perc <- round(pix_3/pix_total*100, digits = 5)
area_4_perc <- round(pix_4/pix_total*100, digits = 5)
area_5_perc <- round(pix_5/pix_total*100, digits = 5)
area_6_perc <- round(pix_6/pix_total*100, digits = 5)
area_7_perc <- round(pix_7/pix_total*100, digits = 5)

porc_sup_quemada <- c(area_na_perc,area_1_perc,area_2_perc,area_3_perc,area_4_perc,area_5_perc,area_6_perc,area_7_perc)

bp <- barplot(dNBR_reclass,
        main = "Porcentaje de cada categoría de severidad",
        col = my_col,
        names.arg = c("A", "B", "C", "D", "E", "F", "G", "H"), 
        las=0)

legend(x='right', y='top', legend =rat$legend, fill = my_col)

text(bp, y=2500, labels = porc_sup_quemada)                          # El valor de y debe ser elegido manualmente

Figure 9. Distribución de valores de píxeles de cada categoría para el AOI seleccionado con porcentaje de área cubierta

ex <- extent(dNBR_reclass)

# Cálculo de la extensión en km2
x <- (ex[2]-ex[1])/1000
y <- (ex[4]-ex[3])/1000
area_total <- x*y

area_na <- round((area_na_perc/100)*area_total, digits = 5)
area_1 <- round((area_1_perc/100)*area_total, digits = 5)
area_2 <- round((area_2_perc/100)*area_total, digits = 5)
area_3 <- round((area_3_perc/100)*area_total, digits = 5)
area_4 <- round((area_4_perc/100)*area_total, digits = 5)
area_5 <- round((area_5_perc/100)*area_total, digits = 5)
area_6 <- round((area_6_perc/100)*area_total, digits = 5)
area_7 <- round((area_7_perc/100)*area_total, digits = 5)

area_cubierta <- c(area_na,area_1,area_2,area_3,area_4,area_5,area_6,area_7)

bp <- barplot(dNBR_reclass,
        main = "Superficie por categoría de NBR",
        col = my_col,
        names.arg = c("A", "B", "C", "D", "E", "F", "G", "H"), 
        las=0)

legend(x='right', y='top', legend =rat$legend, fill = my_col)

text(bp, y=2500, labels = area_cubierta) #y chosen manually

Figure 10.Distribución de valores de píxeles de cada categoría para el AOI seleccionado con área cubierta en kilómetros cuadrados

En los siguientes pasos construimos un rasterBrick con las bandas y nombres organizados para facilitar los pasos siguientes.

nombres <- c("AC","B","G","R","NIR","SWIR1","SWIR2","TERMAL")

Extraer bandas, construir un rasterStack y asignar nombres para la escena pre incendio:

bandas_pre <- file.path(destino_[2], tifs_pre[c(3:9,12)])
imagen_pre <- stack(bandas_pre)
names(imagen_pre) <- nombres

Y la del incendio:

bandas_fire <- file.path(destino_[3], tifs_fire[c(3:9,12)])
imagen_fire <- stack(bandas_fire)
names(imagen_fire) <- nombres

Para cada una, convertimos de DN o valores digitales a reflectancia (ρ) utilizando las constantes definidas en el cuadro 8.2 para las imágenes OLI y TIRS. Ver figura 9.61 para una composición RGB de ambas imágenes procesadas.

imagen_pre_fire <- (imagen_pre[[1:7]] * 0.0000275) + (-0.2)
names(imagen_pre_fire) <- nombres[1:7]
imagen_pre_fire$TERMAL <- (imagen_pre[[8]] * 0.00341802 + 149)
imagen_in_fire <- (imagen_fire[[1:7]] * 0.0000275) + (-0.2)
names(imagen_in_fire) <- nombres[1:7]
imagen_in_fire$TERMAL <- (imagen_fire[[8]] * 0.00341802 + 149)

Y utilizando una sesión interactiva se define un área más ajustada entorno al incendio forestal:

ext <- drawExtent()

Y recortamos nuestra imagen para trabajar en la zona de interés:

area_incendio_pre <- crop(x = imagen_pre_fire, y = ext)
area_incendio <- crop(x = imagen_in_fire, y = ext)

Un interesante ejercicio es comparar combinaciones de banda para conocer el comportamiento de un fuego activo y la pluma de humo (Figura 9.62).

4 OTROS ÍNDICES APROPIADOS PARA EL ANÁLISIS DE LA SEVERIDAD DE LOS INCENDIOS FORESTALES

4.1 Índice Área Quemada BAI

Este índice resalta la tierra quemada en la porción del espectro que va del rojo al infrarrojo cercano, al enfatizar la señal de carbón en las imágenes posteriores al incendio. El índice se calcula a partir de la distancia espectral de cada píxel a un punto espectral de referencia, donde convergen las áreas quemadas recientemente (Chuvievo [1998]).

\[BAI = \frac 1 {(0.1 - ρRED)^2 + (0.06 - ρNIR)^2}\]

bai <- 1 / ((0.1 - area_incendio$R)^2 + (0.06- area_incendio$NIR)^2)
bai[bai > 1000] <- NA

4.2 Índice NBRI

Su fórmula es similar a NDVI, excepto que utiliza longitudes de onda de infrarrojo cercano (NIR) e infrarrojo de onda corta (SWIR) (Key et al., 2002).

\[NBRI = \frac{NIR-SWIR}{NIR+SWIR}\]

nbri <- spectralIndices(area_incendio, 
                        nir ="NIR", 
                        swir3="SWIR2", 
                        indices = "NBRI")

4.3 NBRI Diferenciado, DNBRI

Dentro del programa FIREMON (Key and Benson, 2005), del servicio geológico de EE.UU. (USGS) se publica una categorización de la severidad de incendios basado en este índice:

\[ΔNBRI = NBRI_{pre} − NBRI_{post}\]

Donde \(NBRI_{pre}\) es el índice calculado en la misma zona geográfica previamente al incendio, y \(NBRI_{post}\) el índice calculado después del incendio.

Cuadro 9.4: Categorías de gravedad de quemaduras según proyecto FIREMON (USGS).

Límite Inferior Límite Superior Severidad
-1 -0.25 Alto rebrote post fuego
-0.25 -0.1 Bajo rebrote post fuego
-0.10 1 Sin quemar
0.10 0.27 Baja
0.27 0.44 Moderada a baja
0.44 0.66 Moderada-alta
0.66 1 Alta

Cálculo del índice NBRI en la imagen pre-incendio:

nbri_pre <- spectralIndices(area_incendio_pre, 
                            nir ="NIR", 
                            swir3="SWIR2", 
                            indices = "NBRI")

De la misma manera se puede calcular el índice en la imagen post-incendio, y a partir de ambas se calcula la diferencia:

dnbri <- nbri_pre - nbri

Se puede generar un mapa de severidad usando este índice, al que se aplica un proceso de reclasificación y conversión en factores de las categorías para su posterior tratamiento estadístico. Los pasos sería los siguientes:

  1. Creación de una matriz con rangos y nuevos valores

  2. Reclasificación mediante el uso de la función reclassify(l) donde el ráster (re)clasifica grupos de valores por otros valores definidos en una matriz rcl.

  3. Conversión en factores:

categoria <- ratify(firemon, 
                    count= T, 
                    ordered=T)
  1. Creación de un vector con un texto descriptivo de cada categoría.
clases <- c("Alto rebrote post fuego",
            "Bajo rebrote post fuego",
            "Sin quemar",
            "Quemadura baja gravedad",
            "Quemadura de gravedad moderada a baja",
            "Quemadura de gravedad moderada a alta",
            "Quemadura de alta gravedad")

Se asocia el vector como atributo a la tabla del ráster:

levels(categoria)[[1]] <- cbind(levels(categoria)[[1]], 
                                clase=clases)

A continuación, revisón de la frecuencia por categoría:

levels(categoria)

Y consultar por coordenadas (previa conversión a número de celda):

id <- cellFromXY(categoria,c(617397, 4132006))

Como IDds de celdas (300, 500, 4) se usa la función factorValues(), donde x corresponde a un rasterLayer y v un vector de identificación de celdas:

factorValues(x = categoria, v = categoria[c(id, 300, 500, 4)])

4.4 Índice NBRT1

Este índice utiliza una banda térmica para mejorar el NBR (Holden et al.2005), dado como resultado una mejor separabilidad entre tierra quemada y no quemada.

\[NBRT1 = \frac{ρNIR − ρSWIR (termal / 1000)}{ρNIR + ρSWIR (termal / 1000)}\]

La banda térmica debe calibrarse en relación a la temperatura de brillo (en º Kelvins).

nir <- area_incendio$NIR
swir <- area_incendio$SWIR2

termica <- area_incendio$TERMAL

nbrt1 <- (nir-swir*(termica/1000))/(nir+swir*(termica/1000))