Extraer valores raster por polígonos, líneas y puntos
¿ArcGIS/QGIS para todo?
Muchos de los procesos básicos disponibles en los programas típicos (ArcMap/QGIS) tienen sus homólogos en el entorno R. Con un mínimo esfuerzo, se pueden desarrollar herramientas propias para resultados muy satisfactorios. Si eres de los que dicen que desde que se inventó el modelbuilder escribir código es tontería te diré lo que le dijo El Gran Lewosky a Jesús Quintana: Well…that’s your opinion man.
Este pretende ser el primero, de hecho lo es, de una serie de artículos sobre herramientas básicas de R para mezclar y hacer chisporrotear procesos geográficos. La idea es lanzar un artículo semana, pero quién sabe, porque sobran ideas pero falta tiempo. En cualquier caso, si te pasas por aquí y crees que podría ser interesante algún artículo sobre alguna herramienta geográfica concreta, pasate por el formulario del contacto al final.
La previa…
Primero cargaremos las librerías y capas que utilizaremos para los
consiguientes procesos. Símplemente utiliza raster para las capas
raster y st_read para las capas de información vectorial. En caso de
que fuesen información .gpkg pues haces lo mismo.
# Carga librerías
library(raster)
library(sf)
library(mapview)
library(dplyr)
library(htmlwidgets)
# Carga las capas (corrige la ruta donte tengas las tuyas)
r <- raster('../../../../00_BASE/MDT25/MDT25.tif')
municipios <- st_read('../../../../00_BASE/Concellos/Concellos.shp')
La extracción…
¿Qué tal si calculamos la cota media sobre este MDT de resolución 25 m
para todos los municipios? Podrías hacerlo para todos, pero Galicia
tiene 313 concellos y no se descarta que pueda haber más (o menos)
debido al espíritu troceador de esta preciosa tierra. Por eso, aunque
totos ellos son igual de importantes, lo mejor será extraer los datos
para los 5 primeros registros. Para eso haces un subset ([1:5,]) sobre
la capa, y cuando veas que todo sale a pedir de (Millhouse te tiras a la
piscina)1.
msel <- municipios[1:5,] # subconjunto con los 5 primeros registros
cotamedia <- extract(r, msel, fun=mean)
La función extract devuelve para cada polígono el valor de la media de los registros dentro del mismo. El resultado por lo tanto será una lista de valores (726, 398, 616, 401, 678) que después habrá que agregar como atributo.
La herramienta extract tiene parámetros interesantes como por ejemplo:
na.rmsi queremos que en el cálculo se dejen los valores nulos fuera.method, por defecto está en ‘simple’, lo cual básicamente hace un celda a punto (centro) y si el punto cae dentro del polígono, se cuenta.
Para añadir estos datos como atributo a la capa original habría varías formas. Esta es una (tened en cuenta que solo estarán bien los registros para las tres primeras filas por el subset anterior):
msel$cotamedia <- round(cotamedia)
results <- dplyr::select(data.frame(msel), Concello, cotamedia)
Mmmm, unos datos espectaculares:
| Concello | cotamedia |
|---|---|
| Sarreaus | 726 |
| Taboadela | 398 |
| Bola, A | 616 |
| Maside | 401 |
| Trasmiras | 678 |
En cualquier caso, las posibilidad de esta herramienta son muchas, ya que se pueden extraer datos a partir de polígonos, líneas y puntos, por lo que habría aplicaciones para:
- Extraer la cota de todas las cacas de perro (información vectorial dependiente del dueño), para después hacer un bonito histograma.
- Extraer la precipitación media en los tramos del Camino de Santiago a partir de un raster de precipitaciones. El resultado ya sabemos que es cuanto más cerca de la catedral, más lluvia.
- Extraer extraer y contar un determinado valor de raster por polígono. Por ejemplo si tenemos un raster de color (bandeado) extraer los valores del rojo que están entre 100 y 150. Sentinelfreaks va por ustedes.
La exportación…
Exportar la capa con el nuevo atributo es tan sencillo como:
filename <- 'ConcellosCota.gpkg'
st_write(msel, paste0(cachepath, filename), 'ConcellosCota', append=FALSE)
Vista rápida de los resultados:
-
En el futuro haré algún artículo sobre cómo ejecutar este proceso a toda velocidad, en paralelo, quemando el pc. ↩
