J'extrais la superficie et le pourcentage de couverture de différents types d'utilisation des terres à partir d'un raster basé sur plusieurs milliers de limites de polygones. J'ai trouvé que la fonction d'extraction fonctionne beaucoup plus rapidement si j'itère à travers chaque polygone individuel et recadre puis masque le raster jusqu'à la taille du polygone particulier. Néanmoins, c'est assez lent, et je me demande si quelqu'un a des suggestions pour améliorer l'efficacité et la vitesse de mon code.

La seule chose que j'ai trouvée à ce sujet est cette réponse de Roger Bivand qui a suggéré d'utiliser GDAL.open()et GDAL.close()ainsi que getRasterTable()et getRasterData(). J'ai étudié ces problèmes, mais j'ai eu des problèmes avec gdal dans le passé et je ne le connais pas assez bien pour savoir comment le mettre en œuvre.

Exemple reproductible:

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éthode la plus rapide à ce jour

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 

Traitement parallèle

Le traitement parallèle a réduit de moitié le temps de l'utilisateur, mais a annulé l'avantage en doublant le temps système. Raster l'utilise pour la fonction d'extrait, mais malheureusement pas pour la fonction de recadrage ou de masque. Malheureusement, cela laisse un peu plus de temps total écoulé en raison de «l'attente» par le «IO».

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

exécuter du code sur plusieurs cœurs:

  user  system elapsed 
  23.31    0.68   42.01 

puis terminez le cluster

endCluster()

Méthode lente: la méthode alternative de faire un extrait directement à partir de la fonction raster prend beaucoup plus de temps, et je ne suis pas sûr de la gestion des données pour les mettre sous la forme que je veux:

system.time(ext<-extract(c,bound))
   user  system elapsed 
1170.64   14.41 1186.14 

J'ai enfin réussi à améliorer cette fonction. J'ai trouvé que pour mes besoins, c'était le plus rapide au rasterize()polygone en premier et à utiliser à la getValues()place de extract(). La pixellisation n'est pas beaucoup plus rapide que le code d'origine pour tabuler les valeurs de raster dans de petits polygones, mais elle brille quand il s'agit de grandes zones de polygone qui avaient de grands rasters à recadrer et les valeurs extraites. J'ai également trouvé que getValues()c'était beaucoup plus rapide que la extract()fonction.

J'ai également compris le traitement multicœur à l'aide de foreach().

J'espère que cela est utile pour d'autres personnes qui souhaitent une solution R pour extraire des valeurs raster d'une superposition de polygones. Ceci est similaire à la "Tabulation Intersection" d'ArcGIS, qui ne fonctionnait pas bien pour moi, renvoyant des sorties vides après des heures de traitement, comme cet utilisateur.

#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)

Voici la fonction:

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,single@data$OWNER), showWarnings = FALSE) #creates unique filepath for temp directory
        rasterOptions(tmpdir=file.path("c:/rtemp",i,j, single@data$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(single@data$OWNER,mat) #combines it with the name of the polygon
        unlink(file.path("c:/rtemp",i,j,single@data$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. 
  }
}

Donc, pour l'utiliser, ajustez le single@data$OWNERpour qu'il corresponde au nom de la colonne de votre polygone d'identification (devinez qui aurait pu être intégré à la fonction ...) et insérez:

myoutput <- multicore.tabulate.intersect(cores, polygonlist, rasterlayer)

Accélérez l'extraction du raster (pile de raster) à partir d'un point, d'un XY ou d'un polygone

Excellente réponse Luke. Vous devez être un assistant R! Voici un ajustement très mineur pour simplifier votre code (peut améliorer légèrement les performances dans certains cas). Vous pouvez éviter certaines opérations en utilisant cellFromPolygon (ou cellFromXY pour les points), puis clip et getValues.

Extraire des données de polygones ou de points à partir de piles de rasters ------------------------

 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()

système utilisateur écoulé 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()

système utilisateur écoulé 33,038 3,511 3,288

Si une perte de précision de la superposition n'est pas très importante - en supposant qu'elle soit précise pour commencer - on peut généralement atteindre des vitesses de calcul zonales beaucoup plus importantes en convertissant d'abord les polygones en raster. Le rasterpackage contient la zonal()fonction, qui devrait bien fonctionner pour la tâche prévue. Cependant, vous pouvez toujours sous-définir deux matrices de la même dimension à l'aide d'une indexation standard. Si vous devez conserver des polygones et que le logiciel SIG ne vous dérange pas, QGIS est en fait plus rapide aux statistiques zonales qu'ArcGIS ou ENVI-IDL.

J'ai également eu du mal avec cela pendant un certain temps, en essayant de calculer la part de la superficie des classes de couverture terrestre d'une carte de grille de ~ 300mx300m dans une grille de ~ 1kmx1km. Ce dernier était un fichier polygonal. J'ai essayé la solution multicœur mais c'était encore trop lent pour le nombre de cellules de grille que j'avais. Au lieu de cela, je:

1. Pixellisation de la grille de 1 km x 1 km donnant à toutes les cellules de la grille un numéro unique
2. Utiliser les allign_rasters (ou gdalwarp directement) du package gdalUtils avec l'option r = "near" pour augmenter la résolution de la grille 1kmx1km à 300mx300m, même projection etc.
3. Empilez la carte de couverture terrestre de 300mx300m et la grille de 300mx300m de l'étape 2, en utilisant le package raster: stack_file <- stack (lc, grid).
4. Créez un data.frame pour combiner les cartes: df <- as.data.frame (rasterToPoints (stack_file)), qui mappe les numéros de cellules de la grille de la carte de 1 km x 1 km à la carte de couverture terrestre de 300 m x 300 m
5. Utilisez dplyr pour calculer la part des cellules de classe d'occupation du sol dans les cellules de 1 km x 1 km.
6. Créez un nouveau raster sur la base de l'étape 5 en le liant à la grille d'origine de 1 km x 1 km.

Cette procédure s'exécute assez rapidement et sans problème de mémoire sur mon PC, lorsque je l'ai essayée sur une carte de couverture terrestre avec> 15 cellules de grille de moulin à 300 mx 300 m.

Je suppose que l'approche ci-dessus fonctionnera également si l'on veut combiner un fichier polygonal avec des formes irrégulières avec des données raster. Peut-être, on pourrait combiner les étapes 1 et 2 en tramant directement le fichier polygonal sur une grille de 300mx300 en utilisant rasterize (raster probablement lent) ou gdal_rasterize. Dans mon cas, j'avais également besoin de reprojeter, alors j'ai utilisé gdalwarp pour reprojeter et désagréger en même temps.

Je dois faire face à ce même problème pour extraire des valeurs à l'intérieur du polygone d'une grande mosaïque (50k x 50k). Mon polygone n'a que 4 points. La méthode la plus rapide que j'ai trouvée est de cropmosaïquer en limite de polygone, de trianguler le polygone en 2 triangles, puis de vérifier si les points dans le triangle (l'algorithme le plus rapide que j'ai trouvé). Comparé à la extractfonction, le temps d'exécution passe de 20 s à 0,5 s. Cependant, la fonction cropnécessite encore environ 2 s pour chaque polygone.

Désolé, je ne peux pas fournir l'exemple reproductible complet. Les codes R ci-dessous n'incluent pas les fichiers d'entrées.

Cette méthode ne fonctionne que pour les polygones simples.

par_dsm <- function(i, image_tif_name, field_plys2) {
    library(raster)
    image_tif <- raster(image_tif_name)
    coor <- field_plys2@polygons[[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()