Extraigo el área y el porcentaje de cobertura de diferentes tipos de uso de la tierra de un ráster basado en varios miles de límites de polígonos. Descubrí que la función de extracción funciona mucho más rápido si recorro cada polígono individual y recorto y luego enmascaro el ráster hasta el tamaño del polígono en particular. Sin embargo, es bastante lento, y me pregunto si alguien tiene alguna sugerencia para mejorar la eficiencia y la velocidad de mi código.
Lo único que he encontrado relacionado con esto es esta respuesta de Roger Bivand, quien sugirió usar GDAL.open()y GDAL.close()así como getRasterTable()y getRasterData(). Los investigué, pero he tenido problemas con gdal en el pasado y no lo sé lo suficiente como para saber cómo implementarlo.
Ejemplo reproducible:
library(maptools) ## For wrld_simpl
library(raster)
## Example SpatialPolygonsDataFrame
data(wrld_simpl) #polygon of world countries
bound <- wrld_simpl[1:25,] #name it this to subset to 25 countries and because my loop is set up with that variable
## Example RasterLayer
c <- raster(nrow=2e3, ncol=2e3, crs=proj4string(wrld_simpl), xmn=-180, xmx=180, ymn=-90, ymx=90)
c[] <- 1:length(c)
#plot, so you can see it
plot(c)
plot(bound, add=TRUE)
Método más rápido hasta ahora
result <- data.frame() #empty result dataframe
system.time(
for (i in 1:nrow(bound)) { #this is the number of polygons to iterate through
single <- bound[i,] #selects a single polygon
clip1 <- crop(c, extent(single)) #crops the raster to the extent of the polygon, I do this first because it speeds the mask up
clip2 <- mask(clip1,single) #crops the raster to the polygon boundary
ext<-extract(clip2,single) #extracts data from the raster based on the polygon bound
tab<-lapply(ext,table) #makes a table of the extract output
s<-sum(tab[[1]]) #sums the table for percentage calculation
mat<- as.data.frame(tab)
mat2<- as.data.frame(tab[[1]]/s) #calculates percent
final<-cbind(single@data$NAME,mat,mat2$Freq) #combines into single dataframe
result<-rbind(final,result)
})
user system elapsed
39.39 0.11 39.52
Procesamiento en paralelo
El procesamiento paralelo reduce el tiempo del usuario a la mitad, pero niega el beneficio al duplicar el tiempo del sistema. Raster usa esto para la función de extracción, pero desafortunadamente no para la función de recorte o máscara. Desafortunadamente, esto deja una cantidad ligeramente mayor de tiempo total transcurrido debido a la "espera" del "IO".
beginCluster( detectCores() -1) #use all but one coreejecutar código en múltiples núcleos:
user system elapsed
23.31 0.68 42.01
luego termina el clúster
endCluster()Método lento: el método alternativo para hacer un extracto directamente desde la función de trama tarda mucho más, y no estoy seguro acerca de la gestión de datos para obtener el formato que quiero:
system.time(ext<-extract(c,bound))
user system elapsed
1170.64 14.41 1186.14
Respuestas:
Finalmente pude mejorar esta función. He descubierto que para mis propósitos, era más rápida para
rasterize()el polígono primera y utilizargetValues()en lugar deextract(). La rasterización no es mucho más rápida que el código original para tabular valores de trama en polígonos pequeños, pero brilla cuando se trata de áreas de polígonos grandes que tenían grandes rásteres para recortar y extraer los valores. También encontré quegetValues()era mucho más rápido que laextract()función.También descubrí el procesamiento multinúcleo usando
foreach().Espero que esto sea útil para otras personas que desean una solución R para extraer valores ráster de una superposición de polígonos. Esto es similar a la "Intersección de tabulación" de ArcGIS, que no funcionó bien para mí, devolviendo salidas vacías después de horas de procesamiento, como este usuario.
Aquí está la función:
Por lo tanto, para usarlo, ajústelo
single@data$OWNERpara que se ajuste al nombre de la columna de su polígono de identificación (supongo que podría haberse integrado en la función ...) y agregue:fuente
getValuesfue mucho más rápida que laextractno parece válida porque si la usasextractno tienes que hacercropyrasterize(omask). El código en la pregunta original hace ambas cosas, y eso debería ser el doble de tiempo de procesamiento.Acelere la extracción del ráster (pila de ráster) del punto, XY o polígono
Gran respuesta Luke. ¡Debes ser un mago R! Aquí hay un pequeño ajuste para simplificar su código (puede mejorar ligeramente el rendimiento en algunos casos). Puede evitar algunas operaciones utilizando cellFromPolygon (o cellFromXY para los puntos) y luego recortar y obtener valores.
Extraiga datos de polígonos o puntos de pilas de ráster ------------------------
fuente
Si una pérdida en la precisión de la superposición no es terriblemente importante, suponiendo que sea precisa para comenzar, normalmente se pueden lograr velocidades de cálculo zonales mucho mayores al convertir primero los polígonos en una trama. El
rasterpaquete contiene lazonal()función, que debería funcionar bien para la tarea prevista. Sin embargo, siempre puede subconjugar dos matrices de la misma dimensión utilizando la indexación estándar. Si debe mantener los polígonos y no le importa el software GIS, QGIS es realmente más rápido en estadísticas zonales que ArcGIS o ENVI-IDL.fuente
También luché con esto durante algún tiempo, tratando de calcular la porción de área de las clases de cobertura terrestre de un mapa de cuadrícula de ~ 300mx300m en una cuadrícula de ~ 1kmx1km. Este último era un archivo poligonal. Probé la solución multinúcleo, pero todavía era demasiado lenta para la cantidad de celdas de cuadrícula que tenía. En cambio yo:
Este procedimiento se ejecuta bastante rápido y sin problemas de memoria en mi PC, cuando lo probé en un mapa de cobertura terrestre con> 15 celdas de cuadrícula de molino a 300mx300m.
Supongo que el enfoque anterior también funcionará si se desea combinar un archivo de polígono con formas irregulares con datos ráster. Quizás, uno podría combinar los pasos 1 y 2 rasterizando directamente el archivo de polígono a una cuadrícula de 300mx300 usando rasterizar (raster probablemente lento) o gdal_rasterize. En mi caso, también necesitaba volver a proyectar, así que utilicé gdalwarp para reproyectar y desagregar al mismo tiempo.
fuente
Tengo que enfrentar este mismo problema para extraer valores dentro del polígono de un gran mosaico (50k x 50k). Mi polígono solo tiene 4 puntos. El método más rápido que encontré es hacer un
cropmosaico en el límite del polígono, triangular el polígono en 2 triángulos, luego verificar si los puntos en el triángulo (El algoritmo más rápido que encontré). Compare con laextractfunción, el tiempo de ejecución se reduce de 20 s a 0.5 s. Sin embargo, la funcióncropaún requiere aproximadamente 2 s para cada polígono.Lo siento, no puedo proporcionar el ejemplo reproducible completo. Los códigos R a continuación no incluyen los archivos de entrada.
Este método solo funciona para polígonos simples.
fuente