Wyciągam obszar i procent pokrycia różnych rodzajów użytkowania gruntów z rastra opartego na kilku tysiącach granic wielokąta. Przekonałem się, że funkcja wyodrębniania działa znacznie szybciej, jeśli iteruję po każdym wielokącie i przycinam, a następnie maskuje raster do rozmiaru określonego wielokąta. Niemniej jednak jest dość powolny i zastanawiam się, czy ktoś ma jakieś sugestie dotyczące poprawy wydajności i szybkości mojego kodu.

Jedyne, co znalazłem w związku z tym, to odpowiedź Rogera Bivanda, który zasugerował użycie, GDAL.open()a GDAL.close()także getRasterTable()i getRasterData(). Przyjrzałem się tym, ale w przeszłości miałem problemy z Gdalem i nie znam go wystarczająco dobrze, aby wiedzieć, jak go wdrożyć.

Powtarzalny przykład:

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) 

Jak dotąd najszybsza metoda

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 

Przetwarzanie równoległe

Przetwarzanie równoległe skróciło czas użytkownika o połowę, ale podważyło korzyści wynikające z podwojenia czasu systemowego. Raster używa tego do funkcji wyodrębniania, ale niestety nie do funkcji przycinania lub maskowania. Niestety pozostawia to nieco więcej całkowitego czasu, który upłynął z powodu „oczekiwania” przez „IO”.

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

uruchom kod na wielu rdzeniach:

  user  system elapsed 
  23.31    0.68   42.01 

następnie zakończ klaster

endCluster()

Powolna metoda: alternatywna metoda wyodrębniania bezpośrednio z funkcji rastra zajmuje dużo więcej czasu i nie jestem pewien, czy zarządzanie danymi doprowadzi ją do żądanej formy:

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

W końcu udało mi się ulepszyć tę funkcję. Przekonałem się, że dla moich celów rasterize()najpierw był najszybszy do wielokąta i używał go getValues()zamiast extract(). Rasteryzacja nie jest dużo szybsza niż oryginalny kod do zestawiania wartości rastrowych w małych wielokątach, ale świeci, gdy dojdzie do dużych obszarów wielokąta, które miały duże rastry do przycięcia i wyodrębnienia wartości. Odkryłem również, że getValues()był znacznie szybszy niż extract()funkcja.

Zorientowałem się również przy użyciu przetwarzania wielordzeniowego foreach().

Mam nadzieję, że jest to przydatne dla innych osób, które chcą rozwiązania R do wyodrębniania wartości rastrowych z nakładki wielokąta. Jest to podobne do „Przecięcia tabelarycznego” ArcGIS, które nie działało dla mnie dobrze, zwracając puste wyniki po godzinach przetwarzania, jak ten użytkownik.

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

Oto funkcja:

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

Aby go użyć, dostosuj, single@data$OWNERaby pasował do nazwy kolumny wielokąta identyfikacyjnego (zgadnij, że mogła zostać wbudowana w funkcję ...) i wstaw:

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

Przyspiesz wyodrębnianie rastra (stosu rastrowego) z punktu, XY lub wielokąta

Świetna odpowiedź Luke. Musisz być czarodziejem R. Oto bardzo drobna poprawka w celu uproszczenia kodu (w niektórych przypadkach może nieco poprawić wydajność). Możesz uniknąć niektórych operacji, używając cellFromPolygon (lub cellFromXY dla punktów), a następnie klip i getValues.

Wyodrębnij dane wielokątów lub punktów ze stosów rastrowych ------------------------

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

system użytkownika upłynął 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()

system użytkownika upłynął 33,038 3,511 3,288

Jeśli utrata precyzji nakładki nie jest niezwykle ważna - zakładając, że jest ona precyzyjna na początku - zazwyczaj można osiągnąć znacznie większe prędkości obliczeniowe stref, najpierw przekształcając wielokąty w raster. rasterPakiet zawiera zonal()funkcję, która powinna działać dobrze dla zamierzonego zadania. Zawsze można jednak podzestawić dwie macierze tego samego wymiaru przy użyciu standardowego indeksowania. Jeśli musisz zachować wielokąty i nie przeszkadza ci oprogramowanie GIS, w rzeczywistości statystyki QGIS muszą być szybsze niż ArcGIS lub ENVI-IDL.

Zmagałem się również z tym przez jakiś czas, próbując obliczyć udział w powierzchni klas pokrycia terenu na mapie o wymiarach ~ 300m x 300m na ​​siatce o długości ~ 1 km x 1 km. Ten ostatni był plikiem wielokąta. Wypróbowałem rozwiązanie wielordzeniowe, ale wciąż było to zbyt wolne dla liczby komórek siatki, które miałem. Zamiast tego:

1. Zrasteryzowano siatkę 1 km x 1 km, nadając wszystkim komórkom siatki unikalny numer
2. Zastosowano allign_rasters (lub bezpośrednio gdalwarp) z pakietu gdalUtils z opcją r = "near", aby zwiększyć rozdzielczość siatki 1 km x 1 km do 300 m x 300 m, ta sama projekcja itp.
3. Ułóż mapę pokrycia terenu o wymiarach 300 x 300 m oraz siatkę o wymiarach 300 x 300 m od kroku 2, korzystając z pakietu rastrowego: plik_pliku <- stos (LC, siatka).
4. Utwórz data.frame, aby połączyć mapy: df <- as.data.frame (rasterToPoints (plik_pliku)), który mapuje numery komórek siatki mapy 1 km x 1 km na mapę pokrycia terenu o wymiarach 300 m x 300 m
5. Użyj dplyr, aby obliczyć udział komórek klasy pokrycia terenu w komórkach 1 km x 1 km.
6. Utwórz nowy raster na podstawie kroku 5, łącząc go z oryginalną siatką 1 km x 1 km.

Ta procedura działa dość szybko i bez problemów z pamięcią na moim komputerze, gdy wypróbowałem ją na mapie pokrycia terenu z> 15 milionami komórek siatki o wymiarach 300m x 300m.

Zakładam, że powyższe podejście zadziała również, jeśli chcemy połączyć plik wielokąta o nieregularnych kształtach z danymi rastrowymi. Być może można połączyć krok 1 i 2, bezpośrednio rasteryzując plik wielokąta do siatki o wymiarach 300 x 300 za pomocą rasteryzacji (raster prawdopodobnie wolny) lub gdal_rasterize. W moim przypadku musiałem również ponownie wykonać rzut, więc użyłem gdalwarp zarówno do zmiany wyglądu, jak i do dezagregacji w tym samym czasie.

Muszę zmierzyć się z tym samym problemem, aby wyodrębnić wartości wewnątrz wielokąta z dużej mozaiki (50k x 50k). Mój wielokąt ma tylko 4 punkty. Najszybszą metodą, jaką znalazłem, jest cropmozaikowanie w ramach wieloboku, trójkątowanie wieloboku w 2 trójkąty, a następnie sprawdzanie, czy punkty w trójkącie (najszybszy algorytm, jaki znalazłem). W porównaniu z extractfunkcją czas pracy skraca się z 20 s do 0,5 s. Jednak funkcja cropnadal wymaga około 2 s dla każdego wielokąta.

Niestety nie mogę podać pełnego odtwarzalnego przykładu. Poniższe kody R nie obejmują plików wejściowych.

Ta metoda działa tylko w przypadku prostych wielokątów.

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