Saya mengekstraksi area dan persen tutupan berbagai jenis penggunaan lahan dari raster berdasarkan beberapa ribu batas poligon. Saya telah menemukan bahwa fungsi ekstrak bekerja lebih cepat jika saya mengulangi setiap poligon dan memotong masing-masing lalu menutupi raster ke ukuran poligon tertentu. Meskipun demikian, ini sangat lambat, dan saya bertanya-tanya apakah ada yang punya saran untuk meningkatkan efisiensi dan kecepatan kode saya.

Satu-satunya hal yang saya temukan terkait dengan ini adalah tanggapan oleh Roger Bivand yang menyarankan menggunakan GDAL.open()dan GDAL.close()juga getRasterTable()dan getRasterData(). Saya melihat itu, tetapi memiliki masalah dengan gdal di masa lalu dan tidak cukup tahu untuk menerapkannya.

Contoh yang Dapat Diproduksi:

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) 

Metode tercepat sejauh ini

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 

Proses paralel

Pemrosesan paralel memotong waktu pengguna hingga setengahnya, tetapi meniadakan manfaatnya dengan menggandakan waktu sistem. Raster menggunakan ini untuk fungsi ekstrak, tetapi sayangnya tidak untuk fungsi crop atau mask. Sayangnya, ini meninggalkan jumlah total waktu berlalu yang sedikit lebih besar karena "menunggu" oleh "IO."

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

jalankan kode pada banyak core:

  user  system elapsed 
  23.31    0.68   42.01 

lalu akhiri cluster

endCluster()

Metode Lambat: Metode alternatif melakukan ekstrak langsung dari fungsi raster membutuhkan banyak waktu lebih lama, dan saya tidak yakin tentang manajemen data untuk memasukkannya ke dalam formulir yang saya inginkan:

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

Saya akhirnya berhasil memperbaiki fungsi ini. Saya menemukan bahwa untuk tujuan saya, itu yang tercepat ke rasterize()poligon pertama dan menggunakannya getValues()sebagai gantinya extract(). Rasterisasi tidak jauh lebih cepat daripada kode asli untuk mentabulasi nilai raster dalam poligon kecil, tetapi bersinar ketika datang ke area poligon besar yang memiliki raster besar untuk dipangkas dan nilai diekstraksi. Saya juga menemukan getValues()jauh lebih cepat daripada extract()fungsinya.

Saya juga menemukan proses multi-core menggunakan foreach().

Saya harap ini berguna untuk orang lain yang menginginkan solusi R untuk mengekstraksi nilai raster dari hamparan poligon. Ini mirip dengan "Tabulasi Intersection" dari ArcGIS, yang tidak bekerja dengan baik untuk saya, mengembalikan output kosong setelah jam pemrosesan, seperti pengguna ini.

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

Inilah fungsinya:

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

Jadi untuk menggunakannya, sesuaikan single@data$OWNERagar sesuai dengan nama kolom poligon pengidentifikasi Anda (tebak yang dapat dibangun ke dalam fungsi ...) dan masukkan:

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

Mempercepat penggalian raster (raster stack) dari titik, XY atau Polygon

Jawaban yang bagus Luke. Anda harus menjadi penyihir R! Berikut ini adalah tweak yang sangat kecil untuk menyederhanakan kode Anda (mungkin sedikit meningkatkan kinerja dalam beberapa kasus). Anda dapat menghindari beberapa operasi dengan menggunakan cellFromPolygon (atau cellFromXY untuk poin) dan kemudian klip dan getValues.

Ekstrak data poligon atau poin dari tumpukan raster ------------------------

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

sistem pengguna berlalu 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()

sistem pengguna berlalu 33.038 3.511 3.288

Jika kehilangan dalam presisi overlay tidak terlalu penting - dengan asumsi tepat untuk memulainya - seseorang biasanya dapat mencapai kecepatan perhitungan zona yang jauh lebih besar dengan terlebih dahulu mengubah poligon menjadi raster. The rasterpaket berisi zonal()fungsi, yang harus bekerja dengan baik untuk tugas dimaksud. Namun, Anda selalu dapat mengelompokkan dua matriks dari dimensi yang sama menggunakan pengindeksan standar. Jika Anda harus memelihara poligon dan tidak keberatan dengan perangkat lunak SIG, QGIS sebenarnya harus lebih cepat pada statistik zona daripada ArcGIS atau ENVI-IDL.

Saya juga berjuang dengan ini selama beberapa waktu, mencoba menghitung pangsa area kelas tutupan lahan dari peta grid ~ 300mx300m dalam grid ~ 1kmx1km. Yang terakhir adalah file poligon. Saya mencoba solusi multicore tetapi ini masih terlalu lambat untuk jumlah sel jaringan yang saya miliki. Sebaliknya saya:

1. Rasterisasi kisi 1kmx1km memberi semua sel kisi nomor unik
2. Menggunakan allign_rasters (atau gdalwarp secara langsung) dari paket gdalUtils dengan opsi r = "near" untuk meningkatkan resolusi grid 1kmx1km menjadi 300mx300m, proyeksi yang sama, dll.
3. Tumpuk peta penutup lahan 300mx300m dan kisi 300mx300m dari langkah 2, menggunakan paket raster: stack_file <- stack (lc, grid).
4. Buat data.frame untuk menggabungkan peta: df <- as.data.frame (rasterToPoints (stack_file)), yang memetakan nomor sel jaringan dari peta 1kmx1km ke peta penutup lahan 300mx300m
5. Gunakan dplyr untuk menghitung bagian sel kelas tutupan lahan dalam sel 1kmx1km.
6. Buat raster baru berdasarkan langkah 5 dengan menautkannya ke kisi 1kmx1km asli.

Prosedur ini berjalan sangat cepat dan tanpa masalah memori pada pc saya, ketika saya mencobanya pada peta tutupan lahan dengan sel jaringan> 15 mill pada 300mx300m.

Saya menganggap pendekatan di atas juga akan berfungsi jika seseorang ingin menggabungkan file poligon dengan bentuk tidak teratur dengan data raster. Mungkin, seseorang dapat menggabungkan langkah 1 & 2 dengan langsung meraster file poligon ke grid 300mx300 menggunakan rasterize (raster mungkin lambat) atau gdal_rasterize. Dalam kasus saya, saya perlu memproyeksi ulang juga jadi saya menggunakan gdalwarp untuk memproyeksi dan memilah pada saat yang sama.

Saya harus menghadapi masalah yang sama untuk mengekstrak nilai di dalam poligon dari mosaik besar (50k x 50k). Poligon saya hanya memiliki 4 poin. Metode tercepat yang saya temukan adalah dengan cropmosaic menjadi terikat poligon, triangulasi poligon menjadi 2 segitiga, kemudian periksa apakah titik dalam segitiga (Algoritma tercepat yang saya temukan). Bandingkan dengan extractfungsi, run time dikurangi dari 20 detik menjadi 0,5 detik. Namun, fungsi ini cropmasih membutuhkan sekitar 2 detik untuk setiap poligon.

Maaf saya tidak bisa memberikan contoh lengkap yang bisa direproduksi. Kode R di bawah ini tidak termasuk file input.

Metode ini hanya berfungsi untuk poligon sederhana.

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