Cómo calcular el índice de vegetación MSAVI para suelos con poca densidad

En las dunas camino a Corral, frente a la bahía de Valdivia, hay manchones de matorral costero separados por metros de arena pelada: ni bosque cerrado ni desierto, un término medio que vuelve loco a cualquier índice de vegetación calculado a la rápida. Ese fue el tipo de suelo que me hizo entender que no, el NDVI no sirve para cualquier terreno, y que el MSAVI existe justamente para los casos donde la vegetación es poca y el suelo se ve de sobra. Llevo un tiempo metida en esto de la teledetección por mi cuenta, con RStudio abierto casi todos los sábados, y hace poco me tocó aprender, otra vez a la fuerza, que confiar ciegamente en un solo índice es el error más común que uno comete cuando recién empieza a trabajar con Landsat 8.

Fue Raimundo, un conocido del grupo de senderismo al que a veces me sumo los fines de semana, quien me habló por primera vez de algo llamado NDVI — anda siempre con la cámara tras algún hongo nuevo entre los árboles y quería saber si se podía notar el cambio de la vegetación de un año a otro. En su momento el NDVI me pareció la respuesta a todo. Con el tiempo aprendí que no: funciona de maravilla en un bosque tupido, pero se enreda apenas el terreno se abre.

Por qué el NDVI no funciona bien donde la vegetación es rala

La creencia instalada es que basta con sacar el NDVI y listo, sea cual sea el terreno. Funciona de maravilla en un bosque cerrado, donde la cámara no alcanza a ver ni un pedazo de tierra. El problema aparece en terrenos como el de Corral, o en cualquier ladera con matorral disperso: ahí el suelo desnudo entra en la ecuación y arruina la lectura. El NDVI no distingue bien entre "poca planta" y "tierra que brilla mucho", así que terminas con píxeles que deberían mostrar vida y en cambio aparecen casi muertos en el mapa.

Un suelo claro, arenoso como el de las dunas costeras, refleja mucha luz por sí solo, sin que haya ninguna hoja de por medio. Uno muy oscuro hace lo contrario. En ninguno de los dos casos el NDVI logra aislar bien la señal de la planta, porque está pensado para coberturas densas, no para esos parches donde se ve más tierra que arbusto. Ahí es donde entra la familia de índices ajustados al suelo, que intentan restarle ese ruido de fondo al cálculo.

La familia de índices ajustados al suelo: del SAVI al MSAVI2

Mi primer intento de corregir el problema fue con el SAVI, que pide un factor "L" fijo para compensar el suelo. Si quieres ver cómo me fue exactamente con ese cálculo, lo dejé todo anotado en cómo calcular el índice de vigor de la vegetación SAVI usando RStudio, que fue el paso anterior a este. El problema con el SAVI es que ese número L hay que más o menos adivinarlo según cuánta vegetación crees que hay, y adivinar es exactamente lo que uno no quiere hacer cuando recién está aprendiendo a confiar en sus propios mapas.

Comparación en pantalla entre los valores de NDVI y MSAVI2 sobre vegetación rala captados con Landsat 8

Calculando el MSAVI2 en RStudio, banda por banda

El MSAVI2 resuelve esto solo, sin pedirte ese factor L adivinado: calcula la corrección banda por banda, a partir de las mismas dos bandas que ya usas para el NDVI — el rojo y el infrarrojo cercano, banda 4 y banda 5 si trabajas con Landsat 8. La fórmula asusta un poco la primera vez que la escribes en el editor de RStudio:
MSAVI2 = (2 * NIR + 1 - sqrt((2 * NIR + 1)^2 - 8 * (NIR - Red))) / 2

Lo que más me costó entender es que el 8 que aparece ahí adentro no es un número cualquiera: viene de una simplificación matemática del factor L que antes había que estimar a mano, y es justamente lo que permite que la fórmula se autoajuste sin que una tenga que medir nada en terreno.

Primer plano de manos escribiendo en RStudio el script en R para calcular el índice MSAVI2

Tengo el escritorio de mi pieza armado con el notebook conectado a un monitor prestado, frente a una ventana que en invierno se empaña entera por el lado del patio, y al lado del trackpad siempre hay algo que ya se enfrió. Ahora mismo tengo abiertas en RStudio las bandas de una escena Sentinel-2 que llevo horas sin cerrar ni reiniciar, cosa que hace un tiempo me habría dado pánico — antes cualquier ventana que se quedara pegada me hacía cerrar todo y partir de nuevo por si acaso.

Antes de llegar a esto intenté meterme a artículos científicos sobre teledetección para entender bien de dónde salía toda esta matemática, y la verdad es que el vocabulario técnico me dejó completamente bloqueada — ni iba a entender media página a la primera ni fingir que sí. Terminé aprendiendo la fórmula tal cual, probándola sobre mis propios datos, y recién ahí le fui encontrando el sentido a cada parte.

Correr esto sobre una escena completa no es instantáneo, sobre todo si el raster pesa harto y R se pone lento con la memoria — si te ha pasado lo mismo, dejé anotado aparte cómo evitar errores de memoria al procesar rasters grandes en RStudio, porque ese problema solo me tomó una tarde entera resolverlo la primera vez.

El factor L no es un número fijo, y eso cambia todo

La costumbre más repetida en los tutoriales es usar el SAVI con un L fijo de 0.5, porque así lo enseñan en casi todos lados, como si fuera una regla universal. Lo que vi comparando mis propios mapas es que ese 0.5 funciona razonable en un suelo promedio, pero se desarma en terrenos extremos — en las dunas cerca de Corral, por ejemplo, donde el suelo es arena clara casi todo el año, un L fijo sobreestima la vegetación en algunos sectores y la subestima en otros, según cómo le pegue la luz ese día. El MSAVI2, al recalcular su propio ajuste píxel por píxel, no carga con ese sesgo: se adapta a lo que hay realmente en cada punto del mapa, no a un promedio inventado.

Cuaderno con apuntes a mano sobre las bandas de Landsat 8 y la fórmula del MSAVI2

Cambiar de NDVI a MSAVI2

Conviene mirar la imagen antes de calcular nada: si entre los arbustos se alcanza a ver tierra pelada, o si tu zona de estudio es de pastoreo, de matorral costero o de un área recién quemada, el NDVI probablemente te va a mentir un poco. En cambio, si trabajas sobre un bosque cerrado donde no se ve ni un pedazo de suelo, seguir con el NDVI de toda la vida sigue siendo lo más simple, y no hay necesidad de complicarse con una fórmula más larga.

Isidora, una compañera de este mismo camino autodidacta que conocí comentando los mismos tutoriales en inglés que ninguna de las dos entendía del todo, me escribió hace poco preguntando justo esto para su propio estudio, y fue bueno darme cuenta de que la respuesta no es "usa siempre MSAVI" ni "usa siempre NDVI", sino mirar primero qué tan pelado se ve tu suelo.

Guardando el script para la próxima escena

Dejé el script guardado para la próxima vez, porque mi plan ahora es aplicar esto mismo sobre varias escenas seguidas y ver cómo cambia la cobertura con los años — ando metida en cómo procesar varias imágenes satelitales en R usando bucles sencillos para no tener que repetir cada paso a mano cada vez que baje una escena nueva. Por ahora me quedo con esto: el índice que uses no depende de cuál suene más sofisticado, sino de cuánta tierra se alcanza a ver en tu terreno.

Artículos relacionados