Veinte metros. A ese tamaño de píxel tuve que resignar toda la imagen antes de poder restar dos bandas de Sentinel-2 y ver, recién ahí, dónde había pasado el fuego de verdad. No era el número que esperaba encontrarme cuando empecé a hacer preguntas sobre el índice NBR en RStudio: pensé que el problema iba a ser el fuego, o el algoritmo, o mi poco manejo de R, y resultó que el primer obstáculo fue algo tan tonto como que dos bandas del mismo satélite no vienen ni siquiera del mismo tamaño.
La necesidad de medir esto no salió de la nada. El verano pasado el humo llegó hasta Corral, frente a la bahía de Valdivia, y se quedó ahí varios días, como una nube baja que no se iba ni con el viento del mar. Llevo metida en la teledetección forestal los fines de semana casi dos años, aprendiendo sola a punta de error y de volver a intentar, y cuando vi ese humo lo primero que pensé fue que quería un análisis de incendios propio, con mis datos, no solo con lo que salía en las noticias.
Desde que empecé a escribir sobre esto me han llegado hartas preguntas parecidas, la mayoría de gente que como yo partió de cero con RStudio y el paquete terra. Una de las últimas me la mandó Aldo Pinto, un lector desde Temuco que aparece cada cierto tiempo con una duda nueva, como si no hubiera pasado el tiempo entre un mensaje suyo y el siguiente. Así que en vez de contar esto como un diario, esta vez lo voy a ordenar como las preguntas que de verdad me han hecho, empezando por la más obvia.
¿Qué es el NBR y por qué sirve después de un incendio?
El NBR (Normalized Burn Ratio) compara dos partes de la luz que el ojo no alcanza a ver: el infrarrojo cercano (NIR), que es donde la vegetación sana refleja muchísimo, y el infrarrojo de onda corta (SWIR), que es donde el suelo quemado y la ceniza destacan más. Un bosque sano refleja harto NIR y poco SWIR; después de un incendio esa relación se da vuelta, porque ya no hay hojas verdes reflejando y el suelo queda expuesto. Restando y dividiendo esas dos bandas se arma un mapa donde lo quemado queda marcado con bastante claridad frente a lo que sigue en pie. Esa es la respuesta corta, y es la que doy cada vez que alguien me pregunta para qué sirve este índice en un análisis de incendios: para separar, banda contra banda, lo que el fuego dejó desnudo de lo que el bosque todavía conserva. No es lo mismo que el NDVI, que ya expliqué con calma en otra entrada y que usa otro par de bandas para medir qué tan verde está la vegetación en general; el NBR mira específicamente la cicatriz del fuego, no el verdor.
Para esto uso Sentinel-2, porque las imágenes son gratuitas y porque pasa seguido sobre esta parte del sur de Chile, algo que no es menor cuando la nubosidad de acá se pone pesada media semana sí y media no.
Las dos bandas que de verdad importan (y por qué no calzan a la primera)
Ya entré en detalle sobre qué son las bandas espectrales y cómo leerlas en otra entrada, así que aquí solo dejo las dos que hacen falta para este cálculo puntual. Para Sentinel-2, el NBR se arma con la Banda 8 (NIR), con una longitud de onda central de 833 nm, y la Banda 12 (SWIR), centrada en 2202 nm. El problema, el que me hizo perder tiempo la primera vez, es que estas bandas no comparten el tamaño de píxel: la Banda 8 viene a 10 metros y la Banda 12 a 20 metros. Si intentas restar una banda de 10 metros con una de 20 sin arreglar nada antes, R directamente no te deja avanzar, o peor, te deja avanzar con un resultado que no significa nada. Hay que remuestrear la de mayor detalle para que ambas queden a 20 metros y los píxeles calcen exactamente uno sobre otro. Ese fue el veinte metros del principio: no una elección mía, sino la resolución a la que me tuve que ajustar para que la resta tuviera sentido.
No hace falta procesar la escena completa de Sentinel-2
Una pregunta que me hacen seguido es si hay que cargar la escena entera, con todo lo que pesa, para calcular esto. No. Ya toqué en otra entrada cómo leer archivos raster en RStudio sin ahogarse en el tamaño del archivo, y la idea aquí es la misma: conviene recortar la imagen al área que te interesa antes de hacer cualquier cuenta, algo que también dejé por escrito en una entrada sobre recortar un raster con un shapefile. Para mi caso recorté todo alrededor de Corral y la bahía, que es la zona que realmente quería revisar, y el procesamiento pasó de tardar una eternidad a andar liviano. Si todavía no tienes las imágenes descargadas, escribí antes sobre cómo descargar imágenes Sentinel para procesar en RStudio fácilmente, que fue el paso obligado antes de llegar siquiera a la parte del NBR.
¿Qué pasa si las dos capas no están en la misma proyección?
Otra pregunta típica: qué pasa si una banda quedó en un sistema de coordenadas distinto al de la otra. Ya hay una entrada aparte donde explico con calma qué es una proyección y cómo corregirla en R, así que no me voy a extender aquí, solo dejo el dato puntual para este caso. Para esta zona del sur de Chile el sistema que corresponde es EPSG:32718, zona 18 sur, y antes de restar las bandas conviene revisar que ambas imágenes estén en ese mismo EPSG. Si no lo están, el resultado casi nunca avisa con un error claro: simplemente entrega un mapa de colores que no corresponde a nada real, y ahí es fácil perder una tarde entera sin entender por qué el rojo apareció donde no debería.
Por qué Google Earth no me sirvió para medir nada
Antes de meterme en RStudio probé lo obvio: abrir Google Earth y comparar visualmente las fotos de antes y después del incendio, año contra año. Sirve para hacerse una idea rápida, pero ahí se acaba; no hay forma de medir cuántas hectáreas cambiaron, ni de guardar ese análisis en un archivo que después pueda reabrir, corregir o comparar con otra fecha. Terminé con la sensación de haber mirado el problema sin resolver nada concreto. Eso fue lo que finalmente me hizo sentarme con R en serio: quería un número y un archivo que pudiera repetir la próxima vez que hubiera humo en el cerro, no solo una impresión visual que se me iba a olvidar en un par de semanas.
Y eso sin contar que hay que lidiar con las nubes tapando la escena justo el día que uno la necesita despejada, algo que preferí resolver aparte en otra entrada sobre cómo quitar nubes de las imágenes de Sentinel-2, porque mezclarlo todo acá habría sido demasiado.
NBR simple versus dNBR, la diferencia que de verdad cuenta
El cálculo en sí es corto: nbr = (nir - swir) / (nir + swir). Con eso solo, usando la imagen de después del incendio, ya se ve harto: un mapa que va del verde al rojo, donde el rojo marca la cicatriz sobre el bosque. Pero la respuesta completa a cómo medir el área quemada no es esa sola imagen, es la diferencia entre el NBR de antes y el NBR de después, lo que se conoce como dNBR. Restar esos dos mapas entre sí muestra el cambio real, no solo el estado final, y esto ya roza comparar series completas de imágenes en el tiempo, un tema que da para una entrada propia y que no voy a abrir acá.
Ahí hay una trampa que vale la pena tener presente: en un lugar como el sur de Chile, donde el bosque cambia harto de una estación a otra, si comparas una imagen de verano seco con una de invierno empapado, el agua en las hojas altera los valores y puedes terminar pensando que el daño fue mayor o menor de lo que fue en realidad. Conviene elegir fechas de clima parecido para el antes y el después, o al menos tener clara esa diferencia estacional al leer el resultado. Para confirmar visualmente que lo que el NBR marca como quemado corresponde a suelo desnudo y no a otra cosa, sigo armando composiciones de color en RStudio para ver el bosque nativo, que es simplemente resaltar en una imagen normal lo que uno quiere confirmar.
¿Sirve el NBR para algo más que medir el área quemada?
Bastante gente me pregunta esto, incluyendo un mensaje reciente de Aldo. La respuesta corta es que sí, aunque cada uso merece su propio detalle en otra parte. Los valores de NBR se pueden usar como una capa más dentro de una clasificación de coberturas vegetales, o como insumo en modelos de distribución de especies, algo que también me han preguntado más de una vez. Si después clasificas coberturas con esa capa, conviene validar el resultado con una matriz de confusión antes de confiar en él, porque un mapa bonito no siempre es un mapa correcto. Y si necesitas comparar el valor de NBR con datos tomados en terreno, en puntos específicos del bosque, existen trucos para extraer esos valores del raster directamente en RStudio que ahorran harto trabajo manual. Llegué a todo esto sin ser programadora ni geógrafa, aprendiendo de a poco los fines de semana con lo que iba encontrando en el camino, y con el tiempo terminé conociendo bastantes paquetes de R para trabajar con imágenes satelitales, no solo terra.
¿Cada cuánto puedo repetir este análisis para ver si el bosque se recupera?
Sentinel-2 vuelve a pasar por el mismo punto cada cinco días más o menos, así que en teoría puedo repetir este script cada vez que haya una imagen nueva y sin nubes, y comparar cómo va evolucionando la cicatriz: si aparece verde nuevo o si el suelo sigue pelado. En la práctica lo hago con menos frecuencia, más que nada porque encontrar una imagen limpia de nubes en esta zona no siempre coincide con esos cinco días, pero la idea de fondo es esa: un script que se puede volver a correr, no un análisis que se hace una sola vez y se guarda en un cajón.
La primera vez que corrí el script completo, sin errores, se lo mostré a Marcia, una compañera de la oficina, que celebra cada mapa que le muestro como si fuera un logro propio. Puede sonar poco, pero después de tanto rato mirando código que no corría, ver a alguien más entusiasmarse con un mapa lleno de rojo y verde ayuda a recordar por qué uno se mete en esto. Si te está costando llegar a tu propio mapa, la recomendación es simple: revisa primero que las dos bandas midan a la misma resolución, después que compartan proyección, y recién ahí preocúpate de interpretar los colores. Lo demás, con calma y un poco de lluvia de fondo, va cayendo solo.