Las gotas golpeaban el vidrio con ese ritmo parejo y denso del invierno valdiviano, mientras el ventilador de mi computador empezaba a zumbar con fuerza. En la pantalla, en vez del mapa verde de renovales que esperaba ver surgir capa por capa, apareció esa línea roja que ya me sé de memoria: [rast] extents do not match.
Me quedé mirando el teclado unos segundos, con la taza de té ya fría al lado del mouse pad. Si abría esas mismas dos imágenes en cualquier visor visual de un sistema de información geográfica, se veían una encima de la otra, perfectas. Pero RStudio no funciona con apariencias: funciona con matrices matemáticas puras y duras. Y cuando dos matrices no tienen exactamente las mismas esquinas y la misma cuadrícula, R simplemente se niega a restar un solo píxel.
La ilusión de que dos capas coinciden porque se ven juntas
Todo esto empezó una tarde lluviosa a finales de mayo. Quería hacer algo sencillo —o al menos a mí me parecía sencillo en ese momento—: quería cruzar una escena histórica de Landsat para ver cómo estaba un sector de bosque templado hace algunos años, y compararla con una imagen reciente de Sentinel-2 sobre la misma microcuenca acá en la Región de Los Ríos. Ambas estaban en el mismo sistema de referencia, el clásico EPSG:32718 que usamos para esta franja de Chile, así que en mi cabeza de aprendiz aficionada todo tenía que encajar de inmediato.
Y entonces vino la primera sorpresa. La imagen de Sentinel-2 venía con sus bandas visibles y de infrarrojo cercano en una resolución de 10 metros, mientras que la escena de Landsat traía sus bandas ópticas multiespectrales a 30 metros. En mi ignorancia de los primeros meses, pensé que eso solo significaba que una se veía más nítida que la otra al hacer zoom. Pensé: «bueno, si una tiene píxeles más chicos, R promediará o acomodará las cosas en el fondo». No tenía idea de lo equivocado de esa suposición.
Recuerdo que esa tarde cometí una de mis mayores torpezas de principiante: exporté ambos rasters a GeoTIFF tres veces distintas, cambiando nombres y carpetas, con la secreta esperanza de que al volver a guardar el archivo en el disco el problema del cuadro delimitador se arreglara como por arte de magia. Por supuesto que no cambió nada. El zumbido del computador seguía compitiendo contra la lluvia constante en la ventana, y la consola me seguía arrojando el mismo texto en rojo cada vez que intentaba ejecutar una operación básica como r1 - r2.
Si alguna vez has intentado aprender esto por tu cuenta los fines de semana, sabrás que la tentación en ese punto es cerrar la tapa del notebook e irte a cocinar algo caliente. Yo casi lo hice. Pero después de dar un par de vueltas por el departamento, me senté a revisar con calma qué información me daba la función ext() del paquete terra.
Entendiendo la geometría invisible de un raster
En el paquete terra —que es el estándar que se usa hoy en día en R para mover capas espaciales rápido sin ahogar la memoria del computador—, un objeto raster no es una fotografía decorativa. Es una cuadrícula estricta de números. Para que R pueda restar la banda 8 de Sentinel de una banda de Landsat, o para calcular un índice combinado, el píxel de la fila 1 y columna 1 de la primera imagen debe compartir exactamente las mismas coordenadas físicas en el terreno que el píxel de la fila 1 y columna 1 de la segunda.
Eso significa que no basta con que ambas imágenes cubran Valdivia o la cordillera de la Costa. Tienen que coincidir en cuatro elementos innegociables:
- El sistema de referencia de coordenadas (CRS): Ambas deben hablar el mismo idioma proyectado (por ejemplo, EPSG:32718).
- La extensión espacial (extent): Las coordenadas de los bordes norte, sur, este y oeste deben ser idénticas al milímetro.
- La resolución del píxel: Si una mide 10 metros de lado y la otra 30 metros, una sola celda no puede sumarse con un tercio de otra sin definir antes una regla clara.
- El origen de la malla: Incluso si ambas miden 10 metros y están en el mismo sistema, si una malla parte en la coordenada X=680000 y la otra en X=680003, las líneas de la cuadrícula están desfasadas por tres metros. Para R, eso sigue siendo una incompatibilidad insalvable.
Un domingo por la tarde, a mediados de julio, cuando la humedad en Valdivia no daba tregua y el aire olía a humo de leña en todo el barrio, me puse a experimentar con diferentes comandos para forzar la coincidencia. Casi todos los foros y tutoriales rápidos sugerían lo mismo: «aplica resample() directo con la extensión de la otra imagen y listo». Y sí, al poner esa línea en la consola, el error desaparece al instante y R genera un nuevo mapa sin reclamar. Pero cuando me acerqué a mirar los valores espectrales resultantes, me di cuenta de una trampa silenciosa en la que casi todo el mundo cae al empezar.
El problema de forzar píxeles: por qué prefiero recortar desde el origen
Cuando fuerzas un reescalado directo de una imagen grande para hacerla encajar a la fuerza en la extensión de otra más pequeña o desfasada, el software tiene que inventar o promediar números que no existían. Es lo que hace una interpolación bilineal o un método de vecino más cercano. Y aquí está el meollo del asunto, algo que no entendí hasta bien entrado el invierno: reescalar o ajustar la extensión manualmente a menudo degrada la integridad espectral de tus datos; es preferible realizar un nuevo recorte limpio desde la fuente original usando una máscara común que forzar la interpolación artificial de píxeles sobre bordes artificiales.
Si tomas una imagen de reflectancia y le aplicas transformaciones sucesivas de extensión y resolución solo para callar el mensaje de error de RStudio, estás alterando los valores físicos que midió el sensor en el satélite. Cuando más adelante quieres estudiar la salud del dosel o comparar cambios en la cobertura vegetal, esos pequeños promedios matemáticos acumulados pueden hacerte creer que el bosque perdió vigor cuando en realidad solo fue un artefacto del remuestreo.
Después de casi tres semanas de pruebas durante mis ratos libres, entendí que el camino prolijo consta de tres pasos claros y ordenados en R, sin atajos apurados:
1. Definir una geometría de referencia común
En vez de adivinar cuál imagen debe ceder, elijo cuál de las dos será mi plantilla espacial o defino un polígono vectorial de mi área de estudio (por ejemplo, el límite de la cuenca o la reserva natural). Si ya tienes polígonos dibujados, puedes revisar cómo convertir polígonos a raster en RStudio de forma sencilla para generar esa cuadrícula base perfecta desde cero.
2. Recortar primero, remuestrear después
Si la imagen Sentinel-2 abarca un área gigante y la imagen Landsat cubre otro retazo con límites corridos, lo más sensato es recortar ambas capas a la intersección exacta de sus extensiones usando terra::crop(). Al recortar primero, le ahorras a tu computador la tarea pesada de recalcular miles de píxeles que ni siquiera vas a mirar en el mapa final.
library(terra)
# Supongamos que tenemos r_senti (10 m) y r_landsat (30 m)
# 1. Aseguramos el área de solapamiento común
ext_comun <- intersect(ext(r_senti), ext(r_landsat))
r_senti_recorte <- crop(r_senti, ext_comun)
r_landsat_recorte <- crop(r_landsat, ext_comun)
3. Alinear la malla con el método adecuado
Una vez que ambas capas están circunscritas al mismo cajón geográfico aproximado, recién ahí utilizamos resample() para emparejar la resolución y el origen de las celdas. Y el método que elijas en este punto importa muchísimo:
- Bilineal (
method = "bilinear"): Ideal para variables continuas, como bandas de reflectancia óptica, temperatura de superficie o índices de vegetación. Realiza un promedio ponderado suave entre los cuatro píxeles vecinos. - Vecino más cercano (
method = "near"): Indispensable si estás trabajando con coberturas de suelo categóricas (por ejemplo, 1 = bosque nativo, 2 = plantación, 3 = agua). Si usaras bilineal en una capa temática, terminarías con valores decimales absurdos como 1.4 sobre el mapa.
# Llevamos Landsat a la cuadrícula de Sentinel usando bilineal para bandas continuas
r_landsat_alineado <- resample(r_landsat_recorte, r_senti_recorte, method = "bilinear")
# Ahora ambas comparten extent, crs, resolución y origen exacto
# La resta o apilamiento ya no arrojará ningún error:
diferencia <- r_senti_recorte - r_landsat_alineado
Aprender a mirar la distribución de las bandas después de este paso me salvó de varios disgustos posteriores. Si notas que al alinear se colaron valores extraños o extremos en los bordes por culpa de sombras o esquinas vacías de las escenas satelitales, siempre viene bien darles una pasada con mi método para eliminar valores atípicos en un raster con RStudio antes de seguir calculando índices.
Guardar la receta para cuando regrese la lluvia
No soy programadora profesional ni tengo un título universitario en geomática; todo esto lo he ido desarmando a punta de ensayo, errores tontos y tardes enteras leyendo la documentación oficial de terra mientras la tetera silbaba en la cocina. Durante mucho tiempo le tuve cierto respeto temeroso a la consola de RStudio porque sentía que los errores eran señales de que esto era demasiado difícil para alguien que solo estudia los fines de semana.
Pero con las extensiones de raster me di cuenta de algo liberador: la geometría espacial en un computador no tiene nada de mística ni de caprichosa. R no se equivoca cuando se niega a restar dos imágenes; simplemente te está cuidando la espalda para que no mezcles peras con manzanas sobre el territorio. Si alguna vez te toca comparar bandas como las del infrarrojo y el rojo para ver el verdor de un bosque, tal como cuando vemos cómo comparar índices de vegetación GNDVI y NDVI en RStudio, tener tus cuadrículas calibradas y con su integridad espectral intacta es lo que separa un mapa confiable de un simple dibujo de colores bonitos.
Ahora dejé guardado ese pequeño bloque de recorte e intersección en un script que llamé cariñosamente herramientas_base.R. Afuera ya no llueve con tanta fuerza, solo queda ese goteo suave que baja de los tejados, y en mi pantalla por fin puedo ver cómo el dosel de ulmos y coihues del cerro vecino se dibuja en tonos verdes continuos, píxel por píxel, exactamente en el lugar del mapa donde deben estar.