Pasos para calcular el índice NDVI usando imágenes satelitales en R

Un mapa puede revelar cosas que tus propios ojos todavía no han notado. A mí me pasó con el NDVI, ese índice que en el papel es solo una resta y una división, pero que la primera vez que lo vi pintado sobre un pedazo de bosque me dejó mirando la pantalla más rato del que debería. Ando aprendiendo teledetección para principiantes por mi cuenta, en R, los fines de semana, y aquí junté paso a paso lo que más me costó entender para calcular el NDVI con imágenes Sentinel-2 sobre un trozo de bosque acá en el sur de Chile: qué instalar, qué bandas usar, cómo armar la fórmula y con qué errores rojos te vas a topar en el camino.

Estas son las preguntas que más me hacían cuando recién empezaba, ordenadas más o menos como se me fueron apareciendo a mí. Ninguna la aprendí de un tirón: cada una me costó su propia tarde de ensayo y error, código que no corría y ganas de cerrar el laptop para siempre.

Antes del NDVI: lo que hay que instalar en R paso a paso

R por sí solo no sabe qué hacer con una imagen satelital, así que lo primero es el paquete que la lee: terra. Antes se usaba uno que se llamaba raster, pero ahora se recomienda ir directo a terra porque es más rápido y más liviano. Sumo ggplot2 cuando quiero que el mapa final quede prolijo, aunque para el cálculo del NDVI en sí no hace ninguna falta.

Antes de llegar aquí probé seguir un tutorial en Python, pero daba por sentado que yo ya sabía programar, y a los pocos pasos estaba más perdida que al principio, así que me cambié a R, que tiene una curva bastante más amable para quien parte de cero. Después, cuando la instalación de los paquetes espaciales se me trabó justo en la configuración del PATH, fue mi vecino, Boris, quien me sacó del hoyo: trabaja en infraestructura de servidores y de configurar ambientes de programación sabe muchísimo más que yo.

Si esa parte de leer un archivo raster en R todavía te tiene perdida, tengo aparte cómo leer archivos raster en RStudio para principiantes, que cubre el paso cero de todo esto: el NDVI, en el fondo, no es más que una operación aritmética entre dos capas raster, algo que a mí me costó bastante entender antes de meterme con la fórmula. Terra tampoco anda solo — hay todo un ecosistema de paquetes espaciales en R, lo que mucha gente conoce como r-spatial, pero para lo que necesitas acá con el NDVI, este paquete solo ya te alcanza.

Código en RStudio para calcular el índice NDVI con el paquete terra en R, ejemplo de teledetección para principiantes

Banda roja y banda NIR: las dos que importan en Sentinel-2

Sentinel-2 no saca una foto como la del celular: separa la luz que le llega en varios pedacitos, y para el NDVI solo hacen falta dos. La Banda 4 es la roja, con una longitud de onda central de 665 nm, y es la que nuestros ojos reconocen sin problema. La Banda 8 es el infrarrojo cercano o NIR, centrada en 842 nm, que no vemos pero que las plantas sí aprovechan: absorben buena parte de la luz roja para la fotosíntesis y devuelven casi toda la infrarroja, como si fuera su protector solar.

Tampoco entrega todas sus bandas al mismo tamaño de píxel. Sentinel-2 trae tres resoluciones distintas: 10 metros para las bandas RGB y para esta banda 8 que es el NIR principal, 20 metros para las bandas red-edge y las SWIR, y 60 metros para las de calibración atmosférica, que casi nunca vas a tocar. Como la roja y la NIR principal comparten los 10 metros, calzan perfecto para el NDVI sin tener que remuestrear nada. Ahí está justo la primera trampa que explico más abajo, porque no todas las combinaciones de bandas calzan así de fácil.

¿Por qué R me tira error de extents al restar las bandas?

Porque en algún punto mezclaste una banda de 10 metros con otra de 20, o dos capas que no cubren exactamente el mismo pedazo de terreno. R te va a tirar un mensaje sobre "extents" que no calzan, o te va a pedir que remuestrees, y si no entiendes qué significa eso da bastante susto la primera vez que aparece. En teledetección todo tiene que encajar como un rompecabezas: mismo tamaño de píxel, misma cantidad de filas y columnas, misma posición exacta. Si una capa está corrida un poquito, el cálculo simplemente no corre.

Tengo una conocida de un foro de R, la Denise, que comparte sus pantallazos de error sin ningún filtro ni vergüenza en el grupo. Así que cuando a mí me apareció ese error de extents por primera vez, mandé la captura ahí mismo, taza de té ya helada junto al teclado mientras la consola seguía devolviéndome esa línea roja que no entendía. En un rato alguien me explicó que el problema casi siempre es la proyección de las capas, no la fórmula.

La proyección es justamente el tema que más dolores de cabeza me dio, y por eso escribí aparte sobre cómo corregir proyecciones cartográficas en R de forma sencilla: si tus dos bandas tienen sistemas de proyección distintos, el NDVI te va a salir mal aunque la fórmula esté perfecta, y eso casi no se nota a simple vista en el mapa final.

Cuaderno con la fórmula del NDVI anotada a mano, NIR menos rojo dividido por NIR más rojo, paso a paso

La fórmula del NDVI y qué significan sus valores

Una vez que las dos bandas están alineadas, calcular el NDVI es aritmética pura. La fórmula es (NIR − Rojo) / (NIR + Rojo), y en R con el paquete terra se escribe más o menos así: ndvi <- (banda_nir - banda_roja) / (banda_nir + banda_roja). El resultado es una imagen nueva donde cada píxel tiene un valor entre −1 y +1.

La vegetación sana suele quedar entre 0,3 y 0,8: mientras más alto el número dentro de ese rango, más densa o más vigorosa está la planta. Los valores cercanos a 0 casi siempre son suelo desnudo, pasto seco o superficie construida. Los negativos, sobre todo los que se acercan a −1, suelen ser agua, nieve o nubes, cosas que reflejan la luz de forma completamente distinta a una hoja.

El momento en que corres plot(ndvi) y aparece esa mancha verde marcando justo dónde termina el pasto y empieza el bosque, sin que tú se lo hayas dicho, es la parte que a mí más me convenció de seguir metida en esto.

Comparar NDVI de dos fechas es la trampa que casi todos pisamos

Un curso o un tutorial en algún momento te va a advertir de esto, y aun así casi todo el mundo cae: calcular el NDVI directo sobre dos imágenes de fechas distintas y compararlas sin más, pensando que si el número bajó es porque la planta está peor. La atmósfera (humedad, humo, polvo) ensucia la luz antes de que llegue al satélite, así que un día más húmedo que otro te puede cambiar el NDVI aunque la vegetación esté exactamente igual de sana.

Para comparar fechas en serio hay que usar imágenes que ya vengan corregidas, lo que en Sentinel-2 se conoce como Nivel 2A, o hacer tú misma esa corrección atmosférica antes de restar nada. Yo caí en esta trampa igual, con dos escenas de otoños distintos que no calzaban entre sí, y solo después de leer harto sobre el tema entendí que el problema no era el bosque sino la atmósfera del día en que se tomó cada imagen.

Si quieres el contexto completo de cómo fui llegando a todo esto sin venir de una carrera de programación, ahí quedó mi experiencia aprendiendo teledetección con R sin ser programadora.

¿Me conviene calcular el NDVI sobre toda la escena de Sentinel-2?

Casi nunca. Cada escena de Sentinel-2 cubre una franja enorme de terreno, y si lo que te interesa es un bosque puntual, calcular el NDVI sobre todo eso es puro gasto de memoria y de paciencia esperando a que termine de procesar. Lo que conviene es recortar el raster al área que en verdad te importa antes de hacer la resta, a mí calcularlo sobre la escena completa se me hacía un desperdicio cuando en realidad solo me interesaba el bosque cerca de mi casa, tema que dejo para otro texto porque da para bastante más.

A mí me pasó calculando el NDVI de un sector del Jardín Botánico de la Universidad Austral: apareció un parche rojizo justo donde se abre el sendero que la gente usa en verano, ahí donde tanto paso apisona el pasto hasta dejarlo sin vida. El mapa lo mostraba clarito, sin que yo tuviera que ir a caminar hasta allá para comprobarlo.

Si más adelante te metes en modelado de distribución de especies, el NDVI suele aparecer ahí también como una de las variables de entrada, pero eso ya es otro nivel de complejidad que no toca en este texto.

De todo esto, si hay una sola idea que vale la pena llevarse, es que el NDVI vale lo que valen tus bandas: misma resolución, mismo sistema de proyección, mismo recorte de área, y recién ahí la fórmula hace su trabajo. Todo lo demás, el bosque que quieras mirar y la fecha que quieras comparar, se arma sobre esa base.

Artículos relacionados