Cómo hacer un análisis multiespectral en R usando imágenes de satélite

El raster que necesito para este análisis multiespectral en R lleva minutos cargando, y afuera la lluvia empaña el vidrio antes de que la barra de progreso avance un solo casillero.

Tengo una amiga que hace cerámica en un taller del centro y no tiene idea de sistemas de información geográfica; cuando le mostré el mapa hace poco, me preguntó, medio en broma, si un análisis multiespectral en R era básicamente restarle una foto a otra. Ahí caí en la cuenta de que ese es justo el mito con el que casi todo el que empieza en la teledetección satelital en R se topa: pensar que el análisis es la fórmula, el índice, el resultado bonito en verde. No es así. Antes de calcular cualquier cosa hay pasos que tienen que calzar en orden, y si se salta uno, R simplemente se niega a seguir.

El mito de que un análisis multiespectral en R es solo restar bandas

Código en R para un análisis multiespectral con bandas satelitales de Sentinel-2 en pantalla

Ese mito parte de algo real: al final, sí, buena parte de un análisis multiespectral en R termina siendo una resta o una combinación entre bandas. El error es pensar que ese paso va primero. Antes de restar nada, cada banda necesita el mismo tamaño de píxel y el mismo sistema de referencia de coordenadas, si dos capas quedan en sistemas distintos, ese tema completo lo dejé para otra entrada porque merece explicarse con calma, y si no calzan, R no avanza ni un paso más.

Antes de entender esto vi horas de videos de YouTube sobre SIG, pero casi todos usaban un programa distinto al que tenía instalado en mi computador, así que terminaba mirando botones que ni siquiera existían en mi pantalla, y ninguno mencionaba que el tamaño de píxel entre bandas pudiera ser distinto. Tampoco ayuda que nadie explique primero, en el fondo, qué es un raster; esa base también la dejé para otra entrada porque este texto ya tiene bastante con desarmar el mito, así que si esa parte todavía no te queda clara, conviene resolverla antes de intentar combinar bandas, porque todo lo que sigue depende de eso.

Por qué no todas las bandas son compatibles

Cuaderno con anotaciones a mano comparando bandas satelitales de distinta resolución para un análisis multiespectral

Cada banda capta una porción distinta de la energía que refleja la superficie, y no todas sirven para lo mismo; el detalle banda por banda, y también el cálculo del NDVI, los dejé para la entrada donde armo ese índice paso a paso, porque ahí sí importa explicarlo con calma. Para que las capas geoespaciales conversen entre sí en R uso la librería terra, aunque sobre qué paquete conviene para cada tarea tengo un texto aparte. Lo que sí puedo adelantar acá es que terra permite re-muestrear una banda para que coincida en tamaño con otra, y ese es, literalmente, el paso que resuelve el problema del mito: no es que restar bandas sea difícil, es que primero hay que ponerlas de acuerdo.

El orden real: de la descarga a la capa final

Mapa raster en RStudio con el resultado de un análisis multiespectral en verde sobre gris

Y esto arranca incluso antes de abrir R. No todos los productos de Sentinel-2 que se pueden descargar traen las mismas bandas ni el mismo nivel de corrección, y esa decisión se toma antes de bajar la imagen, no después; algo que expliqué con más detalle en otra entrada porque es de los errores que más tiempo hacen perder. Después viene el problema de las nubes, que en el sur son casi una compañera fija del análisis; taparlas o filtrarlas es un paso propio que también dejé para otro texto. Si además solo me interesa un pedazo del mapa y no la región entera, recorto el raster con un shapefile antes de seguir, un procedimiento que igual merece su propia entrada. Recién ahí, con las bandas del mismo tamaño, en el mismo sistema de coordenadas, sin nubes y recortadas al área que importa, tiene sentido calcular algo: sea un índice, una comparación entre fechas para ver cómo cambió el bosque, que es lo que se conoce como series de tiempo y también da para más, o incluso, para quien lo necesite, alimentar un modelo de dónde podría estar viviendo una especie. Nada de eso funciona si el orden viene al revés, y ese es el error que casi nadie avisa cuando recién se está empezando.

Automatizar todo esto en vez de hacer clic imagen por imagen fue justamente por qué elegí RStudio para estudiar el bosque nativo chileno, aunque esa decisión merece su propio texto. Aprender todo este orden sin venir de la programación es, de nuevo, otra historia; la conté completa en otra entrada, pero el resumen que me sirve acá es más simple: el orden importa más que las ganas de llegar rápido al mapa bonito.

Lo que cambió cuando por fin vi el mapa completo

La primera vez que hice todo esto en el orden correcto, el mapa cargó y apareció, en verde sobre un fondo gris, el contorno completo del humedal del río Cruces, nítido, sin esos bordes cortados o corridos que me habían salido las veces anteriores. No fue un rayo ni una revelación: fue, simplemente, que por primera vez el resultado se parecía a lo que se suponía que debía parecerse. Se lo mostré a mi amiga la de cerámica, y lo único que atinó a decir fue que se parecía a uno de sus esmaltes cuando sale bien. Algo de razón tenía.

Lo que quería ver desde el principio era cómo iba cambiando el tramo del Calle-Calle que pasa cerca de casa, año tras año, y para eso no hace falta ninguna fórmula complicada: hace falta que las bandas calcen, que la proyección sea la misma, que no haya nubes de por medio y que el recorte apunte al lugar correcto. Si alguien te vende el análisis multiespectral como un truco de una sola línea de código, desconfía. La regla que a mí me sirve, y que le paso a cualquiera que pregunte, es simple: antes de escribir la resta entre dos bandas, revisa que compartan resolución y sistema de coordenadas; ese chequeo que toma un minuto evita casi todos los errores que después parecen misteriosos.

Artículos relacionados