El mensaje aparece apenas aprieto enter: 'dimensions do not match'. Ahí me quedo, mirando la consola de RStudio, con el café ya frío al lado del trackpad y dos capas —una de Sentinel-2, otra de Landsat— que se niegan a superponerse en el mapa. Llevo un rato dándole vueltas al mismo problema de teledetección que persigue a media comunidad de R (el #rstats de Twitter está lleno de gente peleándose con esto): los píxeles de una imagen no calzan con los de la otra, y antes de comparar nada tengo que hacer remuestreo espacial con el paquete terra.
Antes de meterme en R, cuando recién estaba curioseando cómo ver el bosque cambiar año tras año, probé comparar imágenes en Google Earth. Ahí no se puede medir nada en serio, y mucho menos guardar algo que después me sirviera para otra cosa, así que terminé por rendirme y buscar algo más serio. Mis imágenes de Sentinel-2 vienen con una resolución de 10 metros en las bandas visibles, mientras que las de Landsat 8 traen 30 metros, y para que ambas hablen el mismo idioma espacial hay que hacer un remuestreo que iguale las grillas.
El paquete terra y la función resample()
Empecé usando el paquete raster porque era el que aparecía en casi todos los tutoriales que encontraba, pero se demoraba una eternidad con cualquier archivo grande y terminé cambiándome a terra, que es más rápido y, por suerte, usa una sintaxis parecida. No voy a comparar acá todos los paquetes de R para trabajar con datos espaciales —hay hartos y cada quien tiene su preferido— pero para el remuestreo, terra es el que me resultó.
La función resample() toma el raster con los píxeles 'que sobran' y proyecta sus valores sobre la grilla que tú le das como modelo: no es estirar la imagen como si fuera un chicle, sino que R recalcula qué valor le corresponde a cada cuadradito nuevo. La primera vez que la corrí con un archivo pesado me tocó esperar mirando la ventana empañada que da al patio interior mientras el procesador hacía lo suyo, sin saber si iba a salir un mapa decente o un error rojo en la consola.
Esto importa más de lo que parece: cuando estaba revisando cómo detectar la pérdida de bosque con imágenes Landsat en RStudio, me di cuenta de que si el remuestreo queda mal hecho, todo el cálculo de cambio de cobertura después no sirve para nada: un píxel que dice 'bosque' en una capa y 'claro' en la otra, y no hay forma de saber a cuál creerle.
Bilinear o vecino más cercano: cuál elegir
Hace un tiempo Marcia, una compañera de la oficina que estudia agronomía a distancia, me mandó un artículo sobre métodos de interpolación espacial que todavía tengo pendiente de leer (me pasa seguido con lo que ella me manda). Pero lo que sí aprendí a la mala es que el método de remuestreo que uses cambia por completo el resultado, y no da lo mismo elegir cualquiera.
El método bilinear funciona como si pasaras un pincel mojado sobre los datos: suaviza los valores y promedia entre píxeles vecinos. Es perfecto para variables continuas —temperatura, elevación, ese tipo de cosas— porque ahí un valor intermedio tiene sentido.
Ahí está el problema real: cuando trabajo con categorías, como cobertura vegetal, donde el valor 1 es bosque y el 2 es suelo desnudo. Si uso bilinear ahí, R me puede inventar un 1.5, y un 1.5 no significa nada —¿medio árbol?, ¿un arbusto confundido?—. Por eso, para mapas de cobertura, uso siempre el método 'near' (vecino más cercano), que toma el valor del píxel original más pegado y lo deja tal cual, sin inventar promedios raros.
Esa diferencia me sirvió después para los pasos para reclasificar un raster en R y simplificar mapas de vegetación, porque si el remuestreo ya viene 'sucio', la reclasificación sale mal de entrada.
Revisa el sistema de referencia antes de remuestrear
Antes de correr resample(), ambos rasters tienen que estar en el mismo sistema de referencia —eso es tema aparte y ya hay harto escrito sobre proyecciones en otra parte, pero acá basta con saber que si no calzan, R igual intenta el remuestreo, y lo que sale en pantalla parece un tablero de ajedrez roto en los bordes. Lo aprendí viendo mi propio mapa salir así de torcido, y ahora lo primero que hago, antes de cualquier otra cosa, es revisar que las dos capas coincidan.
El archivo pesa más de 2 GB: ahí cambia el problema
Aldo Pinto, que escribe seguido desde Temuco, me preguntó hace poco qué pasa cuando el archivo pesa más de 2 GB al momento de remuestrear —es guardaparques y anda metido en monitoreo forestal, así que se topa con esto más que yo. La respuesta corta es que ahí el problema deja de ser el método y pasa a ser la memoria del computador, que es un tema aparte y da para una entrada propia.
Subir la resolución: cuándo no conviene
En los foros muchas veces dicen que lo ideal es llevar todo a la resolución más fina posible —pasar de 30 a 10 metros, por ejemplo— pero eso tiene truco: al subir la resolución (upsampling) no estás agregando información nueva, solo estás repitiendo lo que ya decía el píxel grande en cuadraditos más chicos. Se ve más fino, pero no sabe más que antes.
Para análisis climáticos, o cuando lo que quiero es conservar la señal original sin ruido inventado, prefiero ir al revés: bajar todo a la resolución más gruesa (los 30 metros de Landsat, en mi caso). Suena menos vistoso tener un mapa con píxeles más grandes, pero es la forma más honesta de trabajar con los datos que realmente tengo, no con los que me gustaría tener.
Lo que me quedó claro después de tanto pixel movido
Hace poco le mandé una captura del mapa ya alineado a mi hermana, y le expliqué qué significaba cada color sin mirar mis apuntes ni una vez —antes eso me habría tomado horas de repasar.
Si tuviera que resumirlo en una frase para alguien que recién parte: usa 'near' cuando trabajas con categorías (cobertura, uso de suelo, cualquier cosa donde inventar un valor intermedio no tiene sentido), usa 'bilinear' cuando la variable es continua, revisa que ambos rasters compartan el mismo sistema de referencia antes de tocar resample(), y no asumas que más resolución es siempre mejor: a veces es solo más pixelado sin ser más verdadero.