Un polígono no es un raster, aunque en la pantalla de RStudio ocupen exactamente el mismo pedazo de mapa. Uno es una forma cerrada, la línea que dibuja el borde de un bosque nativo, mientras que el otro es una grilla de píxeles, cada uno con su propio valor. Calcular estadísticas zonales es resolver esa diferencia: pedirle a R que tome todos los píxeles que caen dentro del polígono y devuelva un número que resuma lo que hay ahí adentro, ya sea el promedio de verdor, la variación entre partes del bosque o la suma total. En este tutorial de r-spatial para teledetección en R explico cómo hacer ese cálculo sin que los píxeles del borde inflen o desinflen el resultado sin que te des cuenta.
Antes de meterme en el cálculo mismo doy por sentadas un par de cosas. Que ya sabes abrir un archivo raster en RStudio sin que la sesión se caiga. Que tu imagen y tu polígono de bosque nativo quedaron en el mismo sistema de referencia — ese es un problema aparte, que ya conté en otra entrada del blog, y que en mi caso terminé de confirmar no por un mensaje de error sino porque el borde de un polígono coincidió, en el mapa, con una carretera que me sé de memoria de tanto recorrerla cerca del Parque Saval, en Isla Teja.
Doy por sentado, además, que bajaste tu imagen de Sentinel-2, que le quitaste las nubes de encima y que reconoces sus bandas espectrales lo suficiente como para armar un índice tipo NDVI. Que ese raster ya está recortado al tamaño de tu área de estudio, y que no pesa tanto como para que R se quede pegado al abrirlo. Los paquetes espaciales de R que uso — `sf`, `terra` y `exactextractr` — llegaron después de una prueba fallida: antes había instalado una versión de evaluación de ArcGIS, que expiró antes de que lograra sacar mi primer mapa.
Por qué `extract` se queda corto con un polígono de bosque
La función `extract` de toda la vida hace el trabajo, pero lento, y deja una pregunta sin responder. Cuando un polígono corta un píxel por la mitad, ¿el píxel cuenta o no cuenta? La respuesta no siempre es la misma según la función que uses, y ahí es donde empiezan los problemas: si tu bosque nativo tiene una forma irregular — y casi siempre la tiene — vas a tener decenas de píxeles justo en esa situación.
El paquete `exactextractr` resuelve esa pregunta de otra manera: en vez de decidir si un píxel entra o no entra, calcula qué fracción de cada píxel cae dentro del polígono. Es más rápido que `raster::extract` en polígonos complejos, pero la velocidad no es lo que importa de verdad. Lo que importa es que, si ignoras esa fracción, tus resultados quedan sesgados sin que te des cuenta — un polígono con bordes irregulares puede terminar con un promedio de verdor más alto o más bajo del real, dependiendo de qué haya justo afuera.
El peso exacto de cada píxel del borde
Imagina un píxel que queda solo parcialmente dentro del polígono, con el resto cayendo sobre un camino de tierra vecino. Si R lo cuenta como si estuviera completo adentro, el promedio de verdor de todo el polígono baja artificialmente. Si en cambio solo se cuentan los píxeles que están completamente adentro, se deja fuera justo el borde del bosque, que es donde más cambios ocurren de un año a otro.
`exact_extract()` calcula qué fracción de cada píxel cae dentro del polígono y arma un promedio ponderado con eso, así que ni exagera ni esconde lo que pasa en el borde. Es la diferencia entre un número que se ve bien en la tabla y uno que de verdad representa el estado del bosque que estás midiendo.
Cómo elegir qué estadística pedir por polígono
`exact_extract()` no te da un solo número: te deja elegir qué estadística calcular sobre los píxeles de cada polígono, y ahí conviene pensar antes de tipear. La media es la que más uso, porque resume el estado general de verdor de todo el polígono en un solo valor y es fácil de comparar entre bosques o entre años. La desviación estándar cuenta una historia distinta: si sale alta, es señal de que dentro de ese mismo polígono hay partes sanas y partes que no lo están tanto — un claro, una zona más seca, un sector recién intervenido — mientras que un polígono con desviación baja es más parejo por dentro. La suma, en cambio, casi no la uso para un índice como el NDVI, porque no tiene una lectura intuitiva; sirve más para variables de conteo que para promedios de un rango acotado.
Antes de pedir cualquiera de esas estadísticas conviene revisar que todos los rasters que vas a combinar tengan la misma resolución. Si mezclas capas de distintas fuentes, escribí por separado sobre cómo cambiar la resolución de un raster en r con remuestreo espacial, que es justo el paso que evita que tus estadísticas zonales salgan mezclando peras con manzanas.
Cuando la tabla apareció por fin en el visor de RStudio, sin ninguna línea roja, le mandé captura a Marcia, mi compañera de oficina — no importa qué mapa o qué tabla le mande, siempre contesta con el mismo entusiasmo genuino, y esta vez no fue la excepción.
Después de la tabla: hacia dónde seguir
Esto no es lo mismo que sacar el valor de un raster en un punto exacto, que es otro truco distinto y más simple. Tampoco es clasificar coberturas vegetales completas, ni validar esa clasificación con una matriz de confusión, ni modelar dónde podría vivir una especie con Maxent — para eso primero hay que limpiar los puntos de presencia que vas a usar. Si lo que buscas es la fórmula exacta detrás de índices como el NDVI, el NDWI o el EVI, ese es otro tema aparte, igual que armar un mapa prolijo para un informe técnico.
Aldo, un lector que escribe desde Temuco y que suele desaparecer meses del hilo de comentarios para reaparecer con una duda nueva como si no hubiera pasado el tiempo, me preguntó hace poco algo parecido a esto: cómo comparar el mismo polígono de bosque de un año a otro. Ese es justo el siguiente paso después de tener tu tabla de estadísticas zonales, y sobre eso escribí en cómo analizar la fenología de la vegetación en r con series temporales. Cómo aprendí a armar todo este flujo sin venir de la programación es otra historia, que ya conté aparte.
Si tu polígono es chico o tiene una forma irregular — como casi todos los parches de bosque nativo — conviene usar `exact_extract()` con la media ponderada por fracción de píxel, porque ahí es donde más se nota el sesgo del borde. Si en cambio trabajas con polígonos grandes y compactos, la diferencia entre un extractor y otro es mucho menor, y no vale la pena complicarse. La estadística zonal no es el paso más vistoso de la teledetección, pero es el que convierte una imagen bonita en una tabla que de verdad se puede comparar.