Mis trucos para extraer valores de un raster a puntos en RStudio

Un extract() simple en R devuelve el valor del píxel que queda justo bajo cada punto; un extract() con buffer devuelve el promedio de varios píxeles vecinos alrededor de ese mismo punto, y esa diferencia, que en el papel parece mínima, es la que decide si tu tabla de valores tiene sentido o es puro ruido. Llevo un tiempo extrayendo valores de raster a puntos en RStudio, cruzando coordenadas GPS de unos arrayanes y coigües que tengo marcados cerca de casa con imágenes Sentinel-2, y los paquetes terra y sf terminaron siendo mis herramientas de cabecera para esto. La teledetección con R no perdona los atajos, y este es uno de esos casos donde el método que elijas para pasar de imagen satelital a tabla de datos cambia el resultado más de lo que uno esperaría al principio.

Extraer valores de raster a puntos: dos caminos posibles

El primer camino es el más directo. R busca el píxel bajo cada punto y entrega ese valor tal cual, sin tocarlo, rápido, simple, fácil de explicar si alguien te pregunta qué hiciste. Doy por hecho acá que ya sabes leer un raster en RStudio (esa base la dejé ordenada en otra entrada) y que tus puntos y tu imagen están en el mismo sistema de referencia de coordenadas, porque si no la tabla te sale vacía y ni te avisa por qué. Con eso resuelto, la extracción simple funciona bien cuando el píxel es chico comparado con lo que estás midiendo y el parche de bosque alrededor es bastante parejo.

Ahí es donde entra el paquete terra, que reemplazó al antiguo raster y es harto más rápido para este tipo de operación (no me voy a meter en la lista completa de paquetes espaciales de R, ese repaso ya lo hice en otro lado). Uso sf para manejar los puntos como objetos espaciales y terra para leer la imagen y cruzarla, combinación que para quien recién parte es de las más simples de aprender sin perderse en el camino.

Si todavía estás en la etapa de definir bien tus puntos o tus polígonos antes de llegar a esto, vale la pena revisar cómo manejar datos vectoriales en R para delimitar áreas forestales, porque tener esa base ordenada evita buena parte de los problemas que describo más abajo.

Consola de RStudio mostrando un error de sistema de referencia de coordenadas al extraer valores de raster

Cuándo conviene promediar los píxeles vecinos

La otra opción promedia de forma ponderada los píxeles que rodean al punto, en vez de quedarse con uno solo. Sirve mucho más cuando trabajas con datos continuos, como un índice armado a partir de distintas bandas espectrales de la imagen, el NDVI es el ejemplo típico, y ya lo expliqué con calma en otra entrada, porque un solo píxel puede caer justo en el borde de dos coberturas distintas y darte un valor que no representa nada real.

Noté esto la primera vez que comparé mis dos tablas de valores lado a lado: con el método simple, algunos puntos que yo sabía que estaban en pleno bosque cerrado me daban números raros, como si estuvieran en un claro. Con el promedio de vecinos, esos mismos puntos se acomodaban a lo que yo veía cuando caminaba por ahí. La diferencia no es cosmética, cambia literalmente el número que después usas para cualquier análisis posterior.

Poner los dos métodos a prueba con los mismos puntos

Antes de llegar a esta comparación pasé harto tiempo dando vueltas en falso. Vi horas de videos de YouTube sobre SIG que parecían justo lo que necesitaba, hasta que me daba cuenta de que usaban un programa distinto al que yo tenía instalado en el computador, y ahí quedaba de nuevo en cero, tratando de traducir menús que no calzaban con los míos.

Un domingo cualquiera, de esos grises típicos acá, corrí el script completo de principio a fin y no apareció ni una sola línea roja en la consola. No fue gran cosa como para contarlo en detalle, pero después de tantas tardes cortadas a la mitad por un error, se sintió distinto tener el script corriendo solo, sin que yo tuviera que rescatarlo cada dos minutos.

Una amiga que hace cerámica en un taller del centro me preguntó una vez por qué me demoraba tanto en algo que sonaba tan simple como sacarle un número a una imagen. Le expliqué que es parecido a cuando ella prueba un esmalte en una baldosa de prueba antes de aplicarlo en la pieza completa: si no calibras bien el primer paso, todo lo que sigue sale torcido. Se lo dije caminando entre los puestos del Mercado Fluvial de Valdivia un sábado, más para matar el tiempo que porque tuviera una razón clara para estar ahí.

Puntos GPS superpuestos a una imagen Sentinel-2 antes de aplicar extract() en R

Lo que falla antes de llegar a comparar nada

La mayoría de los tropiezos antes de llegar a esta comparación no tenían que ver con qué método de extracción usar, sino con revisar lo básico antes de confiar en cualquier número. Ahora, antes de mirar un solo valor de la tabla, cuento que el número de filas que salió coincida con el número de puntos que metí; si no calzan, prefiero encontrar el problema ahí antes de seguir adelante con cualquier análisis.

No vengo de la programación ni de la geografía, así que aprendí este flujo por mi cuenta, revisando foros y guardando cada solución que encontraba para no repetir el mismo error (esa parte de aprender sola desde cero la conté completa en otra entrada, acá me quedo con la parte técnica).

El paisaje decide qué método usar

Tabla de valores NDVI extraídos por punto junto a un mapa impreso para comparar métodos

Si tus puntos caen en parches grandes y bastante uniformes (una plantación pareja, un cuerpo de agua grande) el método simple alcanza y de paso corre más rápido, algo que se agradece cuando estás procesando una imagen completa con varias bandas cargadas. Si en cambio tus puntos están cerca de un borde, en un mosaico de coberturas chicas o mezcladas, o si vas a usar esos valores para algo más sensible, por ejemplo alimentar un modelo de distribución de especies donde un solo píxel raro puede torcer todo el resultado, el promedio de vecinos vale la pena aunque tome un poco más de tiempo correr.

Todo esto asume que ya tienes la imagen lista para trabajar, descargada y con las nubes más molestas ya descartadas, porque esos pasos previos son otro tema aparte que también dejé documentado en su momento. Una vez que la tabla de valores está limpia, se puede usar para hartas cosas: yo la ocupé después para reclasificar el raster y simplificar el mapa de vegetación, siguiendo los pasos para reclasificar un raster en R y simplificar mapas de vegetación, pero también sirve como insumo para comparar el mismo punto a través de distintas fechas y ver cómo cambia en el tiempo, o para recortar el mismo raster a la forma de un polígono si lo que necesitas es otra cosa distinta a extraer valores puntuales.

Por ahora me quedo con esta comparación resuelta; el próximo paso que tengo pendiente es meterme a crear composiciones de color en RStudio para ver el bosque nativo con más contraste entre especies, pero eso es cuento para otro sábado.

Artículos relacionados