Ich extrahiere die Fläche und die prozentuale Abdeckung verschiedener Landnutzungstypen aus einem Raster, das auf mehreren tausend Polygongrenzen basiert. Ich habe festgestellt, dass die Extrahierungsfunktion viel schneller funktioniert, wenn ich durch jedes einzelne Polygon iteriere und es beschneide und dann das Raster auf die Größe des jeweiligen Polygons verkleinere. Trotzdem ist es ziemlich langsam und ich frage mich, ob jemand Vorschläge zur Verbesserung der Effizienz und Geschwindigkeit meines Codes hat.

Das einzige, was ich im Zusammenhang damit gefunden habe, ist die Antwort von Roger Bivand, der die Verwendung von GDAL.open()und GDAL.close()sowie getRasterTable()und vorschlug getRasterData(). Ich habe mir diese angeschaut, hatte aber in der Vergangenheit Probleme mit gdal und kenne sie nicht gut genug, um zu wissen, wie man sie umsetzt.

Reproduzierbares Beispiel:

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) 

Schnellste Methode bisher

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 

Parallelverarbeitung

Durch die parallele Verarbeitung wurde die Benutzerzeit halbiert, der Vorteil jedoch durch die Verdoppelung der Systemzeit zunichte gemacht. Raster verwendet dies für die Extraktionsfunktion, aber leider nicht für die Zuschneide- oder Maskenfunktion. Leider verbleibt dadurch etwas mehr Gesamtzeit, da das "IO" "herumwartet".

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

Code auf mehreren Kernen ausführen:

  user  system elapsed 
  23.31    0.68   42.01 

Beenden Sie dann den Cluster

endCluster()

Langsame Methode: Die alternative Methode zum Erstellen eines Extrakts direkt aus der Rasterfunktion dauert viel länger, und ich bin nicht sicher, ob die Datenverwaltung das gewünschte Format annimmt:

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

Ich habe es endlich geschafft, diese Funktion zu verbessern. Ich fand, dass es für meine Zwecke am schnellsten war, zuerst rasterize()das Polygon zu verwenden und getValues()stattdessen extract(). Das Rastern ist nicht viel schneller als der ursprüngliche Code zum Tabellieren von Rasterwerten in kleinen Polygonen, aber es scheint, wenn es um große Polygonbereiche geht, in denen große Raster zugeschnitten und die Werte extrahiert werden müssen. Ich fand auch getValues()war viel schneller als die extract()Funktion.

Ich habe auch die Multi-Core-Verarbeitung mit herausgefunden foreach().

Ich hoffe, dies ist nützlich für andere Personen, die eine R-Lösung zum Extrahieren von Rasterwerten aus einer Polygon-Überlagerung wünschen. Dies ähnelt dem "Tabulate Intersection" von ArcGIS, das für mich nicht gut funktioniert hat und nach stundenlanger Verarbeitung leere Ausgaben zurückgibt , wie dieser Benutzer.

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

Hier ist die Funktion:

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. 
  }
}

single@data$OWNERUm es zu verwenden, passen Sie es an den Spaltennamen Ihres identifizierenden Polygons an (raten Sie, dass dies in die Funktion integriert sein könnte ...) und geben Sie Folgendes ein:

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

Beschleunigen Sie das Extrahieren des Rasters (Raster-Stack) von Punkt, XY oder Polygon

Tolle Antwort, Luke. Sie müssen ein R-Assistent sein! Hier ist eine sehr geringfügige Änderung, um Ihren Code zu vereinfachen (kann in einigen Fällen die Leistung geringfügig verbessern). Sie können einige Operationen vermeiden, indem Sie cellFromPolygon (oder cellFromXY für Punkte) verwenden und dann clip und getValues.

Polygon- oder Punktedaten aus Raster-Stapeln extrahieren ------------------------

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

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

Benutzersystem verstrichen 33.038 3.511 3.288

Wenn ein Verlust in der Genauigkeit der Überlagerung nicht besonders wichtig ist - vorausgesetzt, es ist von vornherein genau -, können Sie in der Regel viel höhere Berechnungsgeschwindigkeiten für Zonen erzielen, indem Sie zuerst die Polygone in ein Raster konvertieren. Das rasterPaket enthält die zonal()Funktion, die für die beabsichtigte Aufgabe gut funktionieren sollte. Sie können jedoch immer zwei Matrizen derselben Dimension mit der Standardindizierung unterteilen. Wenn Sie Polygone pflegen müssen und die GIS-Software nichts ausmacht, muss QGIS in Bezug auf Zonenstatistiken tatsächlich schneller sein als ArcGIS oder ENVI-IDL.

Ich hatte auch einige Zeit damit zu kämpfen und versuchte, den Flächenanteil der Landbedeckungsklassen einer ~ 300mx300m-Gitterkarte in einem ~ 1kmx1km-Gitter zu berechnen. Letzteres war eine Polygondatei. Ich habe die Multicore-Lösung ausprobiert, aber diese war immer noch zu langsam für die Anzahl der Gitterzellen, die ich hatte. Stattdessen ich:

1. Rasterisierte das 1 km x 1 km-Raster und gab allen Rasterzellen eine eindeutige Nummer
2. Verwenden Sie allign_rasters (oder gdalwarp direkt) aus dem gdalUtils-Paket mit der Option r = "near", um die Auflösung des Rasters von 1 km x 1 km auf 300 m x 300 m zu erhöhen, dieselbe Projektion usw.
3. Stapeln Sie die 300 x 300 m große Landbedeckungskarte und das 300 x 300 m große Raster aus Schritt 2 mit dem Rasterpaket: stack_file <- stack (lc, grid).
4. Erstellen Sie einen data.frame, um die Karten zu kombinieren: df <- as.data.frame (rasterToPoints (stack_file)), der die Rasterzellennummern der 1 km x 1 km-Karte auf die 300 m x 300 m große Landbedeckungskarte abbildet
5. Verwenden Sie dplyr, um den Anteil der Zellen der Landbedeckungsklasse an den Zellen von 1 km x 1 km zu berechnen.
6. Erstellen Sie auf der Grundlage von Schritt 5 ein neues Raster, indem Sie es mit dem ursprünglichen Raster von 1 km x 1 km verknüpfen.

Dieser Vorgang läuft ziemlich schnell und ohne Speicherprobleme auf meinem PC ab, als ich ihn auf einer Landbedeckungskarte mit> 15 Mill-Grid-Zellen bei 300 x 300 m ausprobierte.

Ich gehe davon aus, dass der oben beschriebene Ansatz auch funktioniert, wenn eine Polygondatei mit unregelmäßigen Formen und Rasterdaten kombiniert werden soll. Möglicherweise können Sie die Schritte 1 und 2 kombinieren, indem Sie die Polygondatei mit rasterize (wahrscheinlich langsam) oder gdal_rasterize direkt in ein 300 x 300-Raster rastern. In meinem Fall musste ich ebenfalls neu projizieren, also benutzte ich gdalwarp, um gleichzeitig neu zu projizieren und zu disaggregieren.

Ich muss mich dem gleichen Problem stellen, um Werte innerhalb eines Polygons aus einem großen Mosaik (50k x 50k) zu extrahieren. Mein Polygon hat nur 4 Punkte. Die schnellste Methode, die ich gefunden habe, besteht darin, ein cropMosaik in einen Polygonrand zu setzen, ein Polygon in zwei Dreiecke zu triangulieren und dann zu prüfen, ob Punkte im Dreieck vorhanden sind (der schnellste Algorithmus, den ich gefunden habe). Im Vergleich zur extractFunktion reduziert sich die Laufzeit von 20 s auf 0,5 s. Die Funktion cropbenötigt jedoch noch ca. 2 s für jedes Polygon.

Leider kann ich nicht das vollständige reproduzierbare Beispiel bereitstellen. Die folgenden R-Codes enthalten keine Eingabedateien.

Diese Methode funktioniert nur für einfache Polygone.

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