La lluvia golpeaba el vidrio con esa insistencia suave y constante que tenemos en Valdivia casi todo el invierno, y yo estaba ahí, con la taza de té negro enfriándose en la esquina de la mesa, mirando una capa de polígonos que simplemente no quería calzar con mi grilla satelital en RStudio.
No tengo estudios en cartografía ni vengo del mundo de la informática; trabajo medio tiempo en una oficina y los fines de semana me siento a intentar entender qué le pasa a los bosques que veo desde el camino. Y esa tarde, lo único que quería era algo conceptualmente muy simple: tomar unas figuras vectoriales que dibujaban retazos de bosque nativo y transformarlas en una matriz de píxeles para poder cruzarlas con una escena satelital. Parecía un trámite de cinco minutos, pero en R, cuando uno empieza a ciegas, cinco minutos se convierten fácilmente en un sábado entero lleno de frustración.
¿Por qué querríamos convertir un vector en una grilla de píxeles?
Al principio me costaba entender por qué había que dar este paso. Si ya tenía los polígonos delimitados de manera impecable en un archivo vectorial, ¿para qué degradar esos bordes tan limpios a una cuadrícula de cuadritos rígidos? La respuesta me llegó de golpe cuando quise cruzar esos límites con datos satelitales.
Las imágenes satelitales ópticas más comunes, como las de Sentinel-2, trabajan con una cuadrícula fija; de hecho, las bandas visibles tienen una resolución espacial de 10 meters por celda. Cuando intentas hacer operaciones matemáticas directas —como pedirle a R que extraiga valores de reflectancia o calcule estadísticas sólo dentro de tus zonas de interés— mezclar vectores continuos con matrices discretas suele poner a prueba la paciencia de la máquina (y la mía). Si ambos objetos están en el mismo formato matricial, todo fluye en una fracción de segundo.
Básicamente, rasterizar un polígono consiste en 'quemar' o proyectar la silueta de nuestros vectores sobre un molde cuadriculado. Cada celda que caiga bajo el manto del polígono recibe un valor determinado (por ejemplo, 1 si es bosque nativo y 0 si es pastizal o suelo desnudo). Pero para que esa transferencia funcione sin descalces absurdos, ambos elementos tienen que compartir exactamente las mismas reglas de ubicación espacial en el mapa.
Mis tropiezos con paquetes antiguos y matrices en blanco
Durante un fin de semana lluvioso en julio pasado, caí de lleno en la trampa de los tutoriales desactualizados. En internet abundan foros que todavía recomiendan paquetes antiguos como raster, sp y rgdal. Esos paquetes fueron fundamentales en su momento, pero hoy en día varias de esas librerías están discontinuadas o son innecesariamente lentas para quien recién intenta comprender los conceptos básicos de un sistema de información geográfica.
Esa tarde pasé horas pegada frente a un error en rojo que decía algo sobre discrepancias en la extensión espacial. Intenté forzar la conversión exportando el resultado a un archivo GeoTIFF; cuando por fin abrió el archivo para revisar el mapa, el monitor mostró un cuadrado completamente blanco, repleto de valores nulos (NA), sin una sola celda que reflejara los árboles. El olor del té frío en la mesa y ese silencio frustrante de la pantalla en blanco son de esas cosas que te hacen cerrar la tapa del computador por el resto de la jornada.
El problema no era la computadora ni la complejidad de la tarea: era que yo no le estaba entregando al programa un 'molde' o plantilla adecuada sobre la cual dibujar los polígonos. R no puede adivinar qué tamaño de píxel quieres, ni qué coordenadas de origen usar, a menos que se lo digas con precisión milimétrica.
La llegada de terra y la lógica del molde
Después de un par de semanas probando distintos scripts y leyendo la documentación con calma (algo que siempre recomiendo hacer en vez de copiar líneas al azar en foros), me cambié por completo al paquete terra. En consejos para aprender R para geografía si no sabes programar mencionaba lo intimidante que resulta la sintaxis al comienzo, pero terra simplificó todo el flujo de trabajo porque agrupa las funciones espaciales bajo nombres intuitivos y corre sobre un motor mucho más veloz.
La clave para convertir cualquier polígono a raster sin perder los estribos radica en entender una sola idea: la función necesita dos cosas en la mano, el vector que contiene las formas geométricas y un raster vacío que sirva de plantilla o rejilla receptora. Si esa plantilla tiene el sistema de coordenadas correcto y la extensión adecuada, el proceso nunca falla.
Para la zona sur de Chile donde vivo, casi siempre trabajamos bajo la proyección UTM correspondiente; en mi caso, el código EPSG 32718 (que corresponde a WGS 84 / UTM Zona 18S). Si tanto el polígono como la plantilla raster descansan sobre esa misma proyección, la quema de los datos sobre la cuadrícula es limpia y directa.
El problema oculto de los bordes: no todos los píxeles caen al centro
Hacia finales de septiembre, cuando ya dominaba la función básica, me topé con un detalle que casi ningún manual elemental advierte. Estaba analizando unos parches de renovales muy delgados que crecen junto a una quebrada cerca de Niebla, y al rasterizar mis polígonos de forma directa, noté que la superficie total resultante en el mapa rasterizado era bastante menor que la superficie calculada sobre los vectores originales. Habían desaparecido franjas enteras de vegetación ribereña.
¿Qué había ocurrido? Por defecto, la mayoría de los algoritmos de rasterización evalúan si el centro del píxel cae dentro del polígono. Si tu polígono es angosto o tiene bordes muy sinuosos, muchas celdas quedan cruzadas por la mitad por la línea del límite, pero como su punto central queda justo afuera, la función las descarta como si no existieran. Esto introduce un sesgo silencioso pero relevante en áreas pequeñas o fragmentadas.
Convertir polígonos a raster directamente en R a menudo genera estos errores de descarte en el perímetro; por eso descubrí que es mucho más seguro trabajar con cuidado los bordes. En terra, el truco está en utilizar el argumento touches = TRUE dentro de la función rasterize() cuando trabajamos con coberturas lineales o bordes irregulares, o bien rasterizar las líneas del contorno por separado para asegurar que ninguna celda perimetral quede excluida antes de consolidar la máscara final.
El paso a paso en cuatro instrucciones claras
Para no perderse en explicaciones teóricas extensas, el flujo elemental en RStudio se puede condensar en apenas cuatro pasos que no requieren instalaciones complejas ni paquetes en conflicto. Solo necesitamos tener instalado y cargado el paquete principal:
library(terra)
# 1. Cargar el polígono vectorial (por ejemplo, un archivo geopackage o shapefile)
mis_bosques <- vect("bosques_nativos.gpkg")
# 2. Crear una plantilla raster con la extensión del vector y resolución definida
# En este caso, celdas de 10 metros y EPSG:32718 para calzar con Sentinel-2
plantilla <- rast(mis_bosques, resolution = 10, crs = "EPSG:32718")
# 3. Rasterizar el polígono sobre la plantilla
# El argumento 'field' define el valor que tendrán los píxeles; 'touches' rescata bordes
bosque_raster <- rasterize(mis_bosques, plantilla, field = 1, background = 0, touches = TRUE)
# 4. Graficar para comprobar visualmente que todo calce
plot(bosque_raster, col = c("white", "forestgreen"))
Si ya cuentas con una imagen satelital cargada de la misma zona, el paso 2 es todavía más sencillo: en vez de inventar un raster vacío desde cero, simplemente pasas esa imagen como plantilla en el segundo argumento de rasterize(). De esa forma te aseguras de que cada píxel de tu polígono transformado coincida exactamente con la celda correspondiente de la escena satelital, sin desfases de medio milímetro ni errores de remuestreo.
Mirar la pantalla cuando todo por fin encaja
Ver aparecer esa silueta verde sobre el fondo blanco en la pestaña de gráficos de RStudio produce una satisfacción silenciosa, de esas que solo entendemos quienes nos hemos quedado pegados con un script un fin de semana entero. Ver que las pequeñas quebradas de renovales no se perdieron en la cuadrícula y que los bordes respetan la topografía del terreno recompensa con creces las horas de tanteo.
A veces, luego de rasterizar y calcular valores sobre capas derivadas (como índices de verdor o humedad), aparecen valores extraños en los extremos de la escena satelital que pueden arruinar los promedios de la vegetación; para esos casos suelo recurrir a mi método para eliminar valores atípicos en un raster con RStudio antes de ponerme a extraer números definitivos.
También me sirvió mucho aplicar lo que compartí en lo que aprendí organizando mis scripts de R para no perder tiempo, guardando esta pequeña rutina de cuatro líneas en un archivo limpio que puedo volver a consultar sin tener que reinventar la rueda cada vez que descargo un nuevo límite predial. Al final, aprender a manejar datos espaciales sin una formación técnica previa es un proceso lento, lleno de días nublados y pruebas fallidas, pero cuando logras que las piezas conversen entre sí en la pantalla, el mapa deja de ser un problema informático y vuelve a ser simplemente una ventana para mirar el bosque.