¿Aumento de la velocidad de recorte, máscara y extracción de trama por muchos polígonos en R?

26

Estoy extrayendo el área y el porcentaje de cobertura de diferentes tipos de uso de suelo 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 itero a través de cada polígono individual y recorte luego enmascara el ráster hasta el tamaño del polígono en particular. No obstante, 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 por Roger Bivand quien sugirió usar GDAL.open() y GDAL.close() así como getRasterTable() y getRasterData() . Las busqué, pero he tenido problemas con gdal en el pasado y no lo conozco 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) 

El 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([email protected]$NAME,mat,mat2$Freq) #combines into single dataframe
      result<-rbind(final,result)
      })

   user  system elapsed 
  39.39    0.11   39.52 

Procesamiento 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 " esperando "por el" IO "

beginCluster( detectCores() -1) #use all but one core

código de ejecución 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 de realizar un extracto directamente desde la función de ráster toma mucho más tiempo, y no estoy seguro de la administración de datos para incluirlo en el formulario que quiero:

system.time(ext<-extract(c,bound))
   user  system elapsed 
1170.64   14.41 1186.14 
    
pregunta Luke Macaulay 16.01.2015 - 12:38

5 respuestas

21

Finalmente he logrado mejorar esta función. Descubrí que para mis propósitos, era más rápido que el primero fuera del polígono rasterize() y usara getValues() en lugar de extract() . El rasterizado no es mucho más rápido que el código original para tabular los valores ráster en polígonos pequeños, pero brilla cuando se trata de grandes áreas poligonales que tenían rásteres grandes para recortar y los valores extraídos. También encontré que getValues() era mucho más rápido que la función extract() .

También descubrí el procesamiento multi-core 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 tablas" de ArcGIS, que no funcionó bien para mí, devolviendo salidas vacías después de horas de procesamiento, como este usuario.

#initiate multicore cluster and load packages
library(foreach)
library(doParallel)
library(tcltk)
library(sp)
library(raster)

cores<- 7
cl <- makeCluster(cores, output="") #output should make it spit errors
registerDoParallel(cl)

Aquí está la función:

multicore.tabulate.intersect<- function(cores, polygonlist, rasterlayer){ 
  foreach(i=1:cores, .packages= c("raster","tcltk","foreach"), .combine = rbind) %dopar% {

    mypb <- tkProgressBar(title = "R progress bar", label = "", min = 0, max = length(polygonlist[[i]]), initial = 0, width = 300) 

    foreach(j = 1:length(polygonlist[[i]]), .combine = rbind) %do% {
      final<-data.frame()
      tryCatch({ #not sure if this is necessary now that I'm using foreach, but it is useful for loops.

        single <- polygonlist[[i]][j,] #pull out individual polygon to be tabulated

        dir.create (file.path("c:/rtemp",i,j,[email protected]$OWNER), showWarnings = FALSE) #creates unique filepath for temp directory
        rasterOptions(tmpdir=file.path("c:/rtemp",i,j, [email protected]$OWNER))  #sets temp directory - this is important b/c it can fill up a hard drive if you're doing a lot of polygons

        clip1 <- crop(rasterlayer, extent(single)) #crop to extent of polygon
        clip2 <- rasterize(single, clip1, mask=TRUE) #crops to polygon edge & converts to raster
        ext <- getValues(clip2) #much faster than extract
        tab<-table(ext) #tabulates the values of the raster in the polygon

        mat<- as.data.frame(tab)
        final<-cbind([email protected]$OWNER,mat) #combines it with the name of the polygon
        unlink(file.path("c:/rtemp",i,j,[email protected]$OWNER), recursive = TRUE,force = TRUE) #delete temporary files
        setTkProgressBar(mypb, j, title = "number complete", label = j)

      }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) #trycatch error so it doesn't kill the loop

      return(final)
    }  
    #close(mypb) #not sure why but closing the pb while operating causes it to return an empty final dataset... dunno why. 
  }
}

Por lo tanto, para usarlo, ajuste [email protected]$OWNER para que se ajuste al nombre de la columna de su polígono de identificación (supongo que se podría haber incorporado a la función ...) y póngalo en:

myoutput <- multicore.tabulate.intersect(cores, polygonlist, rasterlayer)
    
respondido por el Luke Macaulay 17.01.2015 - 04:16
3

Acelera la extracción de ráster (pila de ráster) desde el punto, XY o Polígono

Gran respuesta Luke. ¡Debes ser un mago R! Aquí hay una modificación muy pequeña para simplificar su código (en algunos casos puede mejorar el rendimiento). Puedes evitar algunas operaciones usando cellFromPolygon (o cellFromXY para puntos) y luego recortar y obtenerValores.

Extrae datos de puntos o polígonos de pilas de ráster ------------------------

 library(raster)  
 library(sp)   

  # create polygon for extraction
  xys= c(76.27797,28.39791,
        76.30543,28.39761,
        76.30548,28.40236,
        76.27668,28.40489)
  pt <- matrix(xys, ncol=2, byrow=TRUE)
  pt <- SpatialPolygons(list(Polygons(list(Polygon(pt)), ID="a")));
  proj4string(pt) <-"+proj=longlat +datum=WGS84 +ellps=WGS84"
  pt <- spTransform(pt, CRS("+proj=sinu +a=6371007.181 +b=6371007.181 +units=m"))
  ## Create a matrix with random data & use image()
  xy <- matrix(rnorm(4448*4448),4448,4448)
  plot(xy)

  # Turn the matrix into a raster
  NDVI_stack_h24v06 <- raster(xy)
  # Give it lat/lon coords for 36-37°E, 3-2°S
  extent(NDVI_stack_h24v06) <- c(6671703,7783703,2223852,3335852)
  # ... and assign a projection
  projection(NDVI_stack_h24v06) <- CRS("+proj=sinu +a=6371007.181 +b=6371007.181 +units=m")
  plot(NDVI_stack_h24v06)
  # create a stack of the same raster
  NDVI_stack_h24v06 = stack( mget( rep( "NDVI_stack_h24v06" , 500 ) ) )


  # Run functions on list of points
  registerDoParallel(16)
  ptm <- proc.time()
  # grab cell number
  cell = cellFromPolygon(NDVI_stack_h24v06, pt, weights=FALSE)
  # create a raster with only those cells
  r = rasterFromCells(NDVI_stack_h24v06, cell[[1]],values=F)
  result = foreach(i = 1:dim(NDVI_stack_h24v06)[3],.packages='raster',.combine=rbind,.inorder=T) %dopar% {
     #get value and store
     getValues(crop(NDVI_stack_h24v06[[i]],r))
  }
  proc.time() - ptm
  endCluster()
  

sistema de usuario transcurrido    16.682 2.610 2.530

  registerDoParallel(16)
  ptm <- proc.time()
  result = foreach(i = 1:dim(NDVI_stack_h24v06)[3],.packages='raster',.inorder=T,.combine=rbind) %dopar% {
        clip1 <- crop(NDVI_stack_h24v06[[i]], extent(pt)) #crop to extent of polygon
        clip2 <- rasterize(pt, clip1, mask=TRUE) #crops to polygon edge & converts to raster
         getValues(clip2) #much faster than extract
  }
  proc.time() - ptm
  endCluster()
  

sistema de usuario transcurrido    33.038 3.511 3.288

    
respondido por el mmann1123 11.04.2016 - 05:55
2

Si una pérdida en la precisión de la superposición no es terriblemente importante, asumiendo que es preciso para comenzar, uno puede lograr velocidades de cálculo zonales mucho mayores al convertir primero los polígonos en un ráster. El paquete raster contiene la función zonal() , que debería funcionar bien para la tarea prevista. Sin embargo, siempre puede subcontratar dos matrices de la misma dimensión utilizando la indexación estándar. Si debe mantener polígonos y no le importa el software GIS, QGIS es en realidad más rápido en las estadísticas zonales que ArcGIS o ENVI-IDL.

    
respondido por el Adam Erickson 16.01.2016 - 08:28
2

También luché con esto durante algún tiempo, tratando de calcular la proporció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 de polígonos. Probé la solución multinúcleo, pero esto seguía siendo demasiado lento para el número de celdas de la cuadrícula que tenía. En su lugar yo:

  1. Rasterizó la cuadrícula de 1kmx1km dando a todas las celdas de la cuadrícula un número único
  2. Usó allign_rasters (o gdalwarp directamente) del paquete gdalUtils con la opción r="near" para aumentar la resolución de la cuadrícula de 1kmx1km a 300mx300m, la misma proyección, etc.
  3. Apile el mapa de cobertura terrestre de 300mx300m y la cuadrícula de 300mx300m del paso 2, usando el paquete ráster: stack_file < - stack (lc, grid).
  4. Cree un data.frame para combinar los mapas: df < - as.data.frame (rasterToPoints (stack_file)), que asigna los números de celdas de la cuadrícula del mapa de 1kmx1km al mapa de cobertura terrestre de 300mx300m
  5. Use dplyr para calcular la proporción de celdas de clase de cobertura terrestre en las celdas de 1kmx1km.
  6. Cree un nuevo ráster sobre la base del paso 5 vinculándolo a la cuadrícula original de 1kmx1km.

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 > Celdas de 15 molinos a 300mx300m.

Supongo que el enfoque anterior también funcionará si uno quiere combinar un archivo de polígono con formas irregulares con datos ráster. Tal vez, se podría combinar el paso 1 y el 2 rasterizando directamente el archivo del polígono a una cuadrícula de 300mx300 usando rasterize (raster probablemente lento) o gdal_rasterize. En mi caso, también necesitaba reproyectar, así que usé gdalwarp para reproyectar y desagregar al mismo tiempo.

    
respondido por el Michiel van Dijk 19.01.2018 - 11:24
0

Tengo que enfrentar este mismo problema para extraer los 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 mosaico de crop en el límite del polígono, triangular el polígono en 2 triángulos, luego verificar si hay puntos en el triángulo (el algoritmo más rápido que encontré). En comparación con la función extract , el tiempo de ejecución se reduce de 20 s a 0.5 s. Sin embargo, la función crop todavía 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.

par_dsm <- function(i, image_tif_name, field_plys2) {
    library(raster)
    image_tif <- raster(image_tif_name)
    coor <- [email protected][[i]]@Polygons[[1]]@coords
    ext <- extent(c(min(coor[,1]), max(coor[,1]), min(coor[,2]), max(coor[,2])))

    extract2 <- function(u, v, us, vs) {
        u1 <- us[2]  - us[1]
        u2 <- us[3]  - us[2]
        u3 <- us[1]  - us[3]
        v1 <- vs[1]  - vs[2]
        v2 <- vs[2]  - vs[3]
        v3 <- vs[3]  - vs[1]
        uv1 <- vs[2] * us[1] - vs[1] * us[2]
        uv2 <- vs[3] * us[2] - vs[2] * us[3]
        uv3 <- vs[1] * us[3] - vs[3] * us[1]

        s1 <- v * u1 + u * v1 + uv1
        s2 <- v * u2 + u * v2 + uv2
        s3 <- v * u3 + u * v3 + uv3
        pos <- s1 * s2 > 0 & s2 * s3 > 0
        pos 
    }

    system.time({
        plot_rect <- crop(image_tif, ext, snap ='out')
        system.time({
        cell_idx <- cellFromXY(plot_rect, coor[seq(1,4),])
        row_idx <- rowFromCell(plot_rect, cell_idx)
        col_idx <- colFromCell(plot_rect, cell_idx)

        rect_idx <- expand.grid(lapply(rev(dim(plot_rect)[1:2]), function(x) seq(length.out = x)))

        pixel_idx1 <- extract2(
            rect_idx[,2], rect_idx[,1], 
            row_idx[c(1,2,3)], col_idx[c(1,2,3)])
        pixel_idx2 <- extract2(
            rect_idx[,2], rect_idx[,1], 
            row_idx[c(1,4,3)], col_idx[c(1,4,3)])
        pixel_idx <- pixel_idx1 | pixel_idx2
        })
    })
    mean(values(plot_rect)[pixel_idx])
}

# field_plys2: An object of polygons
# image_taf_name: file name of mosaic file
library(snowfall)
sfInit(cpus = 14, parallel = TRUE)
system.time(plot_dsm <- sfLapply(
    seq(along = field_plys2), par_dsm, image_tif_name, field_plys2))
sfStop()
    
respondido por el Bangyou 27.06.2016 - 00:15

Lea otras preguntas en las etiquetas