La consola tira el mismo error por segunda vez y freno el cursor justo en esa línea de código. Tengo dos escenas Landsat abiertas en RStudio, una de hace un año y otra de esta temporada, y todavía no se quieren restar entre sí. Detectar pérdida de bosque nativo con teledetección no es apretar un botón: es una cadena de decisiones chicas de análisis espacial (qué escena de Landsat sirve, cómo se alinean las imágenes, qué mancha roja es deforestación de verdad y cuál es solo el paso del invierno) y cualquier eslabón se puede reventar solo. Antes de llegar a este punto probé seguir un tutorial de Python pensado para gente que ya sabía programar, y terminé más perdida que al principio; R, con RStudio abierto al lado y la posibilidad de ir probando línea por línea, me hizo más sentido.
¿Qué imagen de Landsat sirve para este trabajo?
Antes de programar nada, lo primero es elegir bien la escena. Landsat 8 tiene la ventaja de ser gratuita en EarthExplorer y de tener una resolución espacial de 30 metros, así que cada pixel representa un cuadrado de 30 por 30 en el terreno — suficiente para notar si desapareció un parche de bosque, no tanto para distinguir árbol por árbol. El sensor pasa por el mismo punto cada 16 días, y en el sur de Chile eso importa: si justo ese día hay nubes, toca esperar la vuelta siguiente. Filtro las escenas buscando menos de un 10% de cobertura de nubes y, cuando puedo, elijo la misma época del año que la imagen con la que voy a comparar. Bajar y organizar bien las imágenes de origen es, en sí, todo un tema aparte que no alcanzo a cubrir acá.
Apilar las bandas antes de comparar nada
Cada descarga de Landsat 8 trae varios archivos .tif sueltos, uno por banda, y el sensor capta once bandas espectrales distintas, cada una sensible a una parte distinta de la luz — cuál es cuál y para qué sirve cada una da para su propio texto. En R, lo que se hace es apilarlas como si fueran capas de un mismo mapa para poder trabajar con ellas juntas. No me voy a detener en qué es exactamente un raster o cómo funciona esa estructura de grilla — eso también tiene su propio texto —, pero para lo que sigue basta con saber que cada banda es una matriz de números que se puede combinar con las demás. Uso el paquete terra para apilar y operar estas capas porque es más rápido que las alternativas más viejas, aunque hay varios paquetes de R pensados para este tipo de datos y cada uno tiene sus propias mañas.
Dos fechas, un mismo idioma
Restar dos imágenes solo funciona si las dos hablan el mismo idioma de coordenadas. Si una escena quedó en una proyección y la otra en otra, R se niega a operarlas — no tira un error bonito, simplemente no cuadra nada, y esa fue una tarde entera que preferiría no repetir. Corregir el sistema de referencia de coordenadas es un ajuste con su propio paso a paso que no desarrollo acá; lo que importa para este texto es que antes de comparar cualquier par de fechas, hay que confirmar que ambas compartan la misma proyección, o el resultado sale en blanco.
El índice que muestra el pulso del bosque
Con las bandas ya alineadas, calculo un índice de vegetación para cada fecha — el NDVI, que combina dos de esas bandas espectrales que mencioné antes. Cómo se arma ese cálculo paso a paso y qué significa cada valor es material para otro texto; acá lo uso simplemente como una forma de resumir, en un solo número por pixel, qué tan verde y sana se ve la vegetación en ese punto. Antes de meterme en ese cálculo, me gusta armar una composición de falso color con las bandas 5, 4 y 3: el bosque queda en un rojo intenso que ayuda a distinguir a ojo dónde termina un predio y empieza otro, mucho antes de correr cualquier número. Si además necesito delimitar bien esas zonas con un límite exacto, ahí entra manejar datos vectoriales en R para delimitar áreas forestales — recortar un raster con un shapefile es otro proceso aparte que tampoco desarrollo hoy.
La trampa estacional: cuando el invierno parece deforestación
Acá está el verdadero riesgo de este método. Cuando comparé una escena de fines de invierno con otra de pleno verano, aparecieron manchas rojas por todos lados que parecían gritar pérdida de vegetación, y por un momento pensé que habían talado medio cerro. La regla que uso ahora es simple: nunca clasifico pérdida de bosque con un umbral fijo de NDVI, tipo "si bajó de tal número es deforestación", porque en el sur las variaciones estacionales — árboles que pierden hojas, pasto que se seca, luz que cambia — se mezclan demasiado con la señal real. Comparo siempre fechas de la misma época del año, o en su defecto reviso con calma cada mancha roja antes de darla por pérdida real, cruzándola con lo que sé del lugar.
Uno de esos chequeos lo hice hace poco con una escena que incluía el humedal del río Cruces: el mapa terminó de cargar y ahí quedó, quieto, el contorno del río dibujado en verde sobre el gris del resto de la imagen, con la taza de té ya fría al lado del teclado mientras la consola seguía mostrando el mismo mensaje en rojo de un intento anterior. Para confirmar que una mancha realmente correspondía a un cambio en el suelo y no a una sombra o a una nube mal filtrada, reviso algunos puntos que conozco caminando cerca del Río Calle-Calle, en Valdivia, y comparo el valor del pixel exacto ahí con lo que recuerdo haber visto. Eso es, básicamente, extraer valores de un raster a puntos en RStudio, algo que uso más como control de calidad que como método principal.
Hay otros caminos posibles que no tomé para este texto. Si en vez de comparar dos fechas quiero ver la tendencia completa a lo largo de varios años, eso ya es trabajo de series temporales, con su propia forma de armar el análisis. Quien prefiera clasificar coberturas de forma más formal — bosque, suelo desnudo, agua, cada pixel con su etiqueta — puede hacer una clasificación supervisada y después medir qué tan bien salió con una matriz de confusión, algo que yo dejo para cuando necesito más precisión que la que da mirar manchas rojas. También existen formas de enmascarar nubes de manera automática antes de meter las imágenes al análisis, en vez de filtrarlas a mano como hago yo. Y si la pregunta no es dónde se perdió bosque sino dónde podría vivir cierta especie según el hábitat que le queda, eso entra al terreno del modelado de distribución de especies, que es otro mundo aparte.
Guardar el script para la próxima pasada
Lo que sí guardo siempre es el script completo, con los nombres de archivo y los pasos en orden, porque el próximo año, cuando quiera revisar cómo sigue el bosque, solo voy a tener que cambiar las fechas y volver a correrlo. Nadie me enseñó esto en un aula: lo fui armando sola, revisando foros y anotando cada error para no repetirlo — esa forma de aprender por cuenta propia da para su propio texto, así que la dejo ahí. Se lo comenté una vez a Denise Alarcón, una conocida de un foro de R, y me mandó un meme a las dos de la mañana cuando por fin le funcionó algo parecido en su propio proyecto. Y cuando necesito ver dos capas al mismo tiempo y la pantalla del portátil se me queda chica, le pido prestado el monitor extra a Boris, el vecino que trabaja en TI del edificio.
Si hay una regla que me llevo de todo esto es esa: la fecha que elijas para comparar pesa tanto como el cálculo mismo. Un NDVI perfecto entre dos épocas distintas del año no va a decir nada confiable sobre deforestación — va a decir que existe el invierno. Antes de confiar en cualquier mancha roja, primero reviso el calendario de las dos escenas, después el terreno que conozco, y solo al final el número.