A veces un mapa de temperatura superficial hecho con Landsat 8 y R dice que el suelo está a mil grados Celsius. Me hacen esa pregunta seguido, o variaciones de ella, desde mi escritorio en un departamento del centro de Valdivia, con los vidrios empañados por el frío del patio interior y el café que se me olvida al lado del trackpad, frío hace rato. No es la única duda que llega: hay quien pregunta por qué su mapa sale completamente rojo, por qué el número no cuadra con lo que siente al caminar afuera, o de plano por dónde empezar con esto de la teledetección y R. Junté las preguntas que más se repiten y las voy respondiendo acá, con lo que a mí me costó entender primero.
¿Por qué el mapa sale con valores de mil grados?
La respuesta corta es que los datos de Landsat no vienen en grados de ningún tipo. Vienen como Niveles Digitales, números enteros sin sentido físico propio, y si los tratas como si fueran temperatura terminas con un mapa que parece el cráter de un volcán en vez de una tarde nublada en el sur. Para que esos números se conviertan en algo real hay que pasar primero por la radiancia espectral, usando los factores de reescalado que vienen en el archivo de la escena, y de ahí recién se llega a Kelvin. Restas 273.15 y ahí sí, por fin, tienes Celsius.
Antes de este paso había instalado, toda entusiasmada, una versión de prueba de ArcGIS pensando que sería más fácil que aprender código, pero la licencia se venció antes de que lograra armar mi primer mapa completo. Terminé en R casi por descarte, no porque me pareciera lo más simple, y ahora no me imagino haciendo esto de otra forma: al menos el código no caduca a mitad de camino. Todo esto además asume que ya tienes la imagen descargada desde el sitio correcto — conseguir la escena indicada también tiene su maña, pero es cuento aparte.
El archivo MTL, la pieza que nadie explica
Lo que a mí me salvó fue toparme con el archivo MTL, un archivo de texto plano que la mayoría de los tutoriales mencionan de pasada y que en realidad es donde está todo lo importante. Ahí vienen los coeficientes de reescalado y las constantes K1 y K2, que cambian de una escena a otra porque dependen de las condiciones exactas en que el satélite tomó esa imagen puntual. Si usas las constantes de una imagen distinta, por más que se parezcan las fechas, los números salen mal sin que el código marque ningún error — esa es la parte traicionera.
Uso hartos paquetes de R para todo este proceso, sobre todo para leer y manipular las bandas, aunque cuáles me han funcionado mejor da para su propia lista en otro momento. Lo mismo con abrir un archivo raster por primera vez en RStudio: si nunca lo has hecho, ese paso solo ya merece su espacio aparte, porque marea antes de que entiendas bien qué estás mirando. Una vez que la Banda 10 está cargada y los metadatos leídos, la consola deja de escupir números de cinco cifras y empiezan a aparecer valores que se parecen a temperaturas de verdad.
¿Banda 10 o banda 11 de Landsat 8?
Landsat 8 trae dos bandas térmicas, la 10 y la 11, y la respuesta corta es: usa la 10. En foros de gente que entiende más que yo de esto explican que la banda 11 arrastra un problema de calibración por luz dispersa que puede desviar bastante el resultado, así que salvo que estés haciendo algo más avanzado combinando las dos bandas, quédate con la 10 no más. Antes de esto solo sabía crear composiciones de color en RStudio para ver el bosque nativo, y pasar de eso a trabajar con una banda térmica sola fue otro aprendizaje completo.
Ojo con algo que a mí se me olvidó la primera vez: si tu imagen tiene nubes encima, ahí donde tapan, la temperatura sale mal no más, sin ningún aviso. Sacar las nubes de una imagen antes de meterte en estos cálculos es un paso que ahora hago siempre, aunque explicarlo completo es tema para otra entrada. La imagen que uses tiene que venir razonablemente despejada, algo que aquí en el sur, la verdad, cuesta encontrar.
Temperatura de brillo versus temperatura real
Básicamente, lo primero que sale del cálculo es la temperatura de brillo: lo que el satélite "ve" si asumiera que la Tierra entera es un cuerpo negro perfecto, algo que por supuesto no es cierto en ningún lugar del planeta. Ahí entra la longitud de onda central de la Banda 10, que anda por los 10.89 micrómetros — un número minúsculo pero que hace toda la diferencia en la física del cálculo. Sin corregir por esto, la temperatura de brillo sola tiende a subestimar lo complicada que es realmente la superficie que estás mirando.
Corregir esto pasa por calcular la emisividad de la superficie, que en simple es cuánto calor emite de verdad cada tipo de cobertura comparado con ese cuerpo negro ideal. El sensor TIRS de Landsat capta esto con una resolución nativa de 100 metros, aunque el archivo que descargas ya viene remuestreado a 30 para que calce con las demás bandas — y en esos metros de por medio puede haber de todo mezclado: techo, árbol, calle.
¿Alcanza con el NDVI para estimar la emisividad?
Casi todo tutorial enseña lo mismo para estimar la emisividad: sacar el NDVI y usarlo para calcular qué proporción de vegetación hay en cada píxel. El paso a paso completo para calcular el NDVI con bandas espectrales de Sentinel-2 lo dejé para otra entrada; acá simplemente lo uso como un ingrediente más dentro de la fórmula. El problema es que esta simplificación asume que todo lo que no es verde se comporta igual, y esa idea se cae rápido apenas sales de un bosque homogéneo.
Fue un artículo que me pasó Marcia — compañera de la oficina que anda siempre mandándome papers que después ella misma no alcanza a leer — el que me hizo dudar en serio de esta simplificación. Ahí entendí que corregir la emisividad solo con el NDVI sirve para partir, pero se queda corto si de verdad quieres algo preciso, sobre todo apenas se mete cemento o asfalto en la escena.
Zonas urbanas: donde el método se cae
En zonas urbanas densas — pienso en el centro de Valdivia, con sus techos de zinc apretados uno junto al otro — este método simplificado se equivoca bastante, porque el cemento, el asfalto y el metal emiten calor de formas muy distintas entre sí y muy distintas a un suelo desnudo de tierra. Al final, mi primer mapa completo decía que ciertas cuadras estaban más frías de lo que yo sentía caminando por ahí un día de sol, que es la señal más clara de que la corrección se quedó corta.
Si quieres limpiar un poco tus datos antes de llegar a este punto, ayuda revisar los pasos para reclasificar un raster en R y simplificar mapas de vegetación, para separar mejor lo que es bosque de lo que es construcción antes de sacar cualquier promedio. Y si además tu mapa sale corrido de lugar, sin calzar con la realidad, probablemente el problema esté en la proyección — eso también es tema aparte, uno bien común entre quienes recién empiezan.
¿Qué pasa si el archivo pesa más de 2 GB?
Un lector de Temuco, Aldo, me escribió hace poco con una pregunta bien puntual: qué hacer cuando el archivo de Landsat pesa más de 2 GB y RStudio empieza a ir lento con cada operación. La respuesta práctica es que casi nunca necesitas la escena completa: si tu zona de interés es un barrio o una cuenca chica, recortar el raster con un shapefile antes de calcular nada ahorra memoria y tiempo de espera, aunque ese proceso de recorte también merece su propia explicación en detalle más adelante.
Nada de esto tiene que ver, ojo, con modelar dónde podría vivir una especie usando algo como Maxent, aunque alguna vez me preguntaron si la temperatura superficial sirve como variable para ese tipo de modelo — puede ser un dato más entre varios, pero es un tema aparte que no voy a abrir acá.
Ver el cemento y el bosque desde la pantalla
Recuerdo la primera vez que el raster completo cargó sin errores y apareció en el visor, verde sobre gris, la curva reconocible del río Cruces: ahí supe que los números por fin correspondían a algo real y no a un experimento fallido más. Ese sábado había ido a caminar un rato por la costanera del Calle-Calle antes de sentarme de nuevo frente al computador, y volver con la cabeza despejada ayudó más que cualquier tutorial.
No soy geógrafa ni programadora, y cómo llegué a meterme sola en todo esto es una historia que ya conté en otra entrada. Cada vez que salgo y siento el calor que rebota del cemento en pleno centro, pienso en esos 10.89 micrómetros y en todo lo que hay que corregir antes de confiar en un número. Con el tiempo he ido guardando estos mapas para comparar uno con otro, aunque armar una serie completa que muestre cómo cambia el bosque año tras año es un proyecto que dejo pendiente. Si te preguntas por qué estudiar teledetección satelital con R para ver cambios ambientales vale la pena a pesar de lo tedioso del proceso, para mí la respuesta está en ese momento exacto en que el mapa por fin se parece a lo que uno ve por la ventana.