El cielo sobre Isla Teja tenía esa luz plana y gris que no deja sombra en nada, ni en el pasto del Parque Saval ni en los árboles que bordean el sendero de al lado. Llevaba un rato mirando una imagen de teledetección de los cerros cerca de mi casa, y la pregunta que me daba vueltas era simple: de toda la luz que le llega al bosque, ¿cuánta se queda ahí adentro y cuánta rebota de vuelta hacia el espacio? Esa fracción que rebota tiene nombre, se llama albedo, y la receta corta para calcularla con imágenes Landsat en R es más o menos esta: tomar varias bandas ya corregidas, multiplicarlas cada una por su propio coeficiente, sumarlas, y restarles al final un ajuste chiquito. Lo largo, como casi siempre en este análisis espacial de aprendiz, es todo lo que hay que entender antes de llegar a esa línea de código.
Lo que el NDVI no alcanza a contarte
Durante bastante tiempo creí que con el NDVI ya tenía resuelto todo lo que necesitaba para seguir mis bosques desde la pantalla. Ese índice de vegetación dice qué tan verde se ve algo, nada más; no dice cuánta energía devuelve esa superficie hacia arriba. Antes de llegar a R probé meterme a esto por otro lado: seguí un tutorial en Python que daba por hecho que uno ya sabía programar, y como esa base yo no la tenía, lo abandoné a las pocas páginas sin entender casi nada de lo que estaba escribiendo. Fue más tarde, leyendo hilos de foros donde gente que sabe mucho más que yo discutía el balance de energía en la superficie, que me topé con el concepto del albedo: una propiedad física sin unidades, que va de 0 a 1, y que representa la fracción de radiación que una superficie refleja. La nieve fresca puede llegar a 0.9, mientras que un bosque nativo denso como el que tengo cerca se queda casi siempre bajo 0.15, porque absorbe la mayor parte de esa energía para sus propios procesos.
Las bandas de Landsat 8 que tuve que aprender a nombrar
Calcular esto no es soplar y hacer botellas, como decimos aquí. Al principio pensé que iba a existir alguna función en R que se llamara calcular_albedo() y ya, lista la pega. Qué ingenua. Lo primero con que uno choca es la corrección atmosférica: no se pueden usar los números digitales crudos que trae la imagen, hay que llegar antes a la reflectancia de superficie. Menos mal que ya me había dado varios cabezazos con lo que aprendí al quitar nubes de imágenes Landsat en RStudio, porque al menos entendía que las imágenes se limpian antes de pedirles peras al olmo.
Para este experimento usé Landsat 8, que trae un sensor llamado OLI bastante conocido entre la gente de teledetección, pero para el albedo de banda ancha, el que sirve para mirar el balance de energía completo, no basta con una sola banda. Se necesitan varias: la B2 (azul), la B3 (verde), la B4 (roja), la B5 (infrarroja cercana), y las B6 y B7, que son las de infrarrojo de onda corta.
Alinear esas seis bandas me tomó varias tardes de puro ensayo y error. Usé el paquete terra, que es el sucesor del paquete raster, y mi laptop sonaba como que iba a pedir licencia médica: el ventilador a full mientras afuera la lluvia caía parejo sobre el techo de zinc. Hubo momentos en que RStudio simplemente se cerraba sin aviso, sin ningún mensaje que uno pudiera buscar en algún lado. Ahí me acordé de unos apuntes que había escrito antes sobre cómo evitar errores de memoria al procesar rasters grandes en RStudio, y tuve que aplicarme mis propios consejos: procesar por pedazos, limpiar la memoria temporal, tener paciencia con una máquina que no está hecha para esto.
La importancia de la reflectancia de superficie
Una trampa de principiante en la que caí fue meter los valores originales de la banda directo en la fórmula. Las imágenes de Landsat 8 vienen escaladas, y si no se aplica ese factor, que para los datos de Nivel 2 suele ser 0.0000275, el resultado no tiene ningún sentido físico. Me acuerdo perfecto de la sensación de vacío cuando mi primer mapa de albedo me tiró valores de 1500 en algunos píxeles. El albedo máximo posible es 1, así que ver un 1500 en pantalla me hizo sentir que nunca iba a entender esto. Un café, volver al código línea por línea, y ahí apareció el error.
Armando la fórmula de Liang
Después de buscar bastante, encontré que para Landsat 8 una de las formas más aceptadas es usar los coeficientes de Liang, publicados originalmente en 2001: una fórmula de integración que combina las bandas para darte ese valor de banda ancha, y que no hay que tenerle miedo aunque de primera vista los números intimiden un poco. Para este sensor, el coeficiente de la banda B2 es 0.356, el de la B3 es 0.130, el de la B4 es 0.373, el de la B5 es 0.085 y el de la banda que cierra la suma es 0.072. Es, en el fondo, una suma ponderada: cada banda multiplicada por su coeficiente, y todas sumadas entre sí.
Cuando corrí esa combinación por primera vez y el mapa cargó sin manchones rojos ni error de sintaxis, apareció un gradiente de reflectividad bastante parejo. El agua absorbía casi toda la energía, con un albedo bajísimo, mientras que techos de bodegas industriales en las zonas más urbanas de alrededor brillaban con valores mucho más altos.
La cobertura del suelo también pasa la cuenta
Mientras comparaba mis resultados con algunos mapas globales, noté algo que me dejó pensando: calcular el albedo de banda ancha con una fórmula genérica sobreestima el error en zonas como esta, porque la humedad y lo complejo que es la estructura del bosque acá cambian bastante cómo se comporta la señal. Los coeficientes de Liang son un punto de partida, no una verdad absoluta para cualquier tipo de cobertura.
Por ejemplo, si hay mucha sombra proyectada por los cerros, el albedo sale bajísimo, casi como si fuera agua, y eso no es necesariamente cierto: es solo que el sensor no está recibiendo luz de vuelta desde ese píxel. Me sirvió mucho volver a mis propios apuntes sobre cómo organicé mis scripts de R para no perder tiempo buscando entre carpetas, porque pude cruzar el mapa de albedo con un modelo de elevación y revisar dónde la sombra me estaba engañando. Fue por esas semanas, cerca de la séptima, cuando además se me arregló un error de proyección que me tenía dando vueltas: cuando quedó bien alineada, el borde del bosque en la pantalla calzaba justo con la curva de la carretera que conozco de memoria de tanto recorrerla en auto, y ahí el mapa dejó de sentirse como un dibujo abstracto. Por esos mismos días Denise, una chica que conocí en un grupo de R, me mandó una captura de su propio mapa con la proyección completamente torcida, sin ningún drama en mostrar el error, y eso me quitó algo de la vergüenza de equivocarme tanto.
El script, al final
Al final, lo que hago es cargar las bandas con rast(), multiplicarlas por su factor de escala para que queden entre 0 y 1, y después aplicar la ecuación de Liang. Queda, más o menos, así:
albedo <- (0.356 * B2) + (0.130 * B3) + (0.373 * B4) + (0.085 * B5) + (0.072 * B7) - 0.0018
Ese número negativo tan chico al final es una constante de ajuste que venía en el paper original, no algo que uno se saque de la manga. Ver que los valores finales quedaban mayoritariamente entre 0.05 y 0.4 para la zona me dio una paz mental increíble: estaba dentro del rango físico que corresponde.
Boris, el vecino del piso de abajo que trabaja en TI, se asomó justo cuando tenía este mapa abierto en la pantalla y se quedó mirando el contraste de colores sin entender mucho, preguntando nomás por qué un lado se veía tan oscuro.
Una caminata por Parque Saval y la lección que me llevo
Esa tarde cerré el script y me quedé mirando el mapa final con una calma que entiende cualquiera que haya pasado semanas aprendiendo a leer píxeles desde cero. Al día siguiente salí a caminar por el Parque Saval, en Isla Teja, y mirando los árboles de cerca pensé en lo raro que es que ahora yo sepa, aunque sea un poco, cuánta energía de esa que cae sobre ellos se queda adentro y cuánta rebota de vuelta.
Landsat 8 tiene una resolución temporal de dieciséis días que en esta zona nos hace sufrir, porque a veces pasan meses antes de conseguir una imagen sin nubes, pero cuando se logra capturar el momento justo, calcular su albedo se siente como ver algo que a simple vista no se puede. Si hay algo que me llevo de estas semanas es que un valor que sale fuera de rango no es un fracaso: es información, te está diciendo que falta revisar un paso atrás en la cadena, una escala, una corrección, una proyección, antes de seguir para adelante. Ya estoy pensando en cómo usar todo esto para mirar el estrés hídrico en verano, pero ese será tema para otro rato.