Я извлекаю площадь и процент покрытия различных типов землепользования из растра, основанного на нескольких тысячах границ полигонов. Я обнаружил, что функция извлечения работает намного быстрее, если я перебираю каждый отдельный многоугольник и обрезаю, а затем маскирую растр до размера определенного многоугольника. Тем не менее, это довольно медленно, и мне интересно, есть ли у кого-нибудь предложения по повышению эффективности и скорости моего кода.

Единственное, что я нашел в связи с этим, - это ответ Роджера Биванда, который предложил использовать GDAL.open()и GDAL.close()так же, как getRasterTable()и getRasterData(). Я изучал их, но в прошлом у меня были проблемы с gdal, и я не знаю этого достаточно хорошо, чтобы знать, как это реализовать.

Воспроизводимый пример:

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) 

Самый быстрый метод на данный момент

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 

Параллельная обработка

Параллельная обработка сократила время пользователя вдвое, но свела на нет преимущество, удвоив системное время. Растр использует это для функции извлечения, но, к сожалению, не для функции обрезки или маски. К сожалению, из-за «ожидания» с помощью «IO» остается немного больше общего истекшего времени .

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

запустить код на нескольких ядрах:

  user  system elapsed 
  23.31    0.68   42.01 

затем завершить кластер

endCluster()

Медленный метод: альтернативный метод извлечения непосредственно из растровой функции занимает намного больше времени, и я не уверен насчет управления данными, чтобы получить его в нужной форме:

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

Я наконец-то дошел до улучшения этой функции. Я обнаружил, что для моих целей он был быстрее всего rasterize()на полигоне и использовал getValues()вместо extract(). Растеризация не намного быстрее, чем исходный код для табулирования растровых значений в маленьких полигонах, но она сияет, когда дело доходит до больших областей многоугольника, в которых были обрезаны большие растры и извлекались значения. Я также обнаружил, что getValues()намного быстрее, чем extract()функция.

Я также разобрался с использованием многоядерной обработки foreach().

Я надеюсь, что это полезно для других людей, которые хотят R-решение для извлечения растровых значений из наложения полигонов. Это похоже на «пересечение таблиц» в ArcGIS, которое у меня не получалось, возвращая пустые выходные данные после нескольких часов обработки, как этот пользователь.

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

Вот функция:

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$OWNERчтобы оно соответствовало имени столбца вашего идентифицирующего многоугольника (угадайте, что могло быть встроено в функцию ...) и вставьте:

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

Ускорение извлечения растра (растрового стека) из точки, XY или Polygon

Отличный ответ Люк. Вы должны быть волшебником R! Вот очень незначительный трюк для упрощения вашего кода (в некоторых случаях может немного повысить производительность). Вы можете избежать некоторых операций, используя cellFromPolygon (или cellFromXY для точек), а затем обрезать и получить значения.

Извлечение данных многоугольника или точек из растровых стеков ------------------------

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

Пользовательская система прошла 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()

Пользовательская система прошла 33.038 3.511 3.288

Если потеря точности наложения не очень важна - при условии, что она точна с самого начала - обычно можно достичь гораздо более высоких скоростей зонального расчета, сначала преобразовав полигоны в растр. В rasterпакете содержится zonal()функция, которая должна хорошо работать для поставленной задачи. Однако вы всегда можете поднастроить две матрицы одного и того же измерения, используя стандартную индексацию. Если вам нужно поддерживать полигоны и вы не против программного обеспечения ГИС, QGIS на самом деле должен работать быстрее в зональной статистике, чем ArcGIS или ENVI-IDL.

Я также некоторое время боролся с этим, пытаясь вычислить долю площади классов земного покрова на карте сетки ~ 300 х 300 м в сетке ~ 1 км х 1 км. Последний был полигональным файлом. Я попробовал многоядерное решение, но оно все еще было слишком медленным для количества ячеек сетки. Вместо этого я:

1. Растеризовал сетку 1 км x 1 км, дав всем ячейкам сетки уникальный номер
2. Использовал allign_rasters (или непосредственно gdalwarp) из пакета gdalUtils с опцией r = "near", чтобы увеличить разрешение сетки 1 км x 1 км до 300 x 300 м, та же проекция и т. Д.
3. Сложите карту земного покрова размером 300 х 300 м и сетку размером 300 х 300 м, начиная с шага 2, используя растровый пакет: stack_file <- stack (lc, grid).
4. Создайте файл data.frame для объединения карт: df <- as.data.frame (rasterToPoints (stack_file)), который отображает номера ячеек сетки карты 1 км x 1 км на карту земного покрова размером 300 х 300 м.
5. Используйте dplyr для расчета доли ячеек класса растительного покрова в ячейках 1 км х 1 км.
6. Создайте новый растр на основе шага 5, связав его с исходной сеткой 1 км х 1 км.

Эта процедура выполняется довольно быстро и без проблем с памятью на моем компьютере, когда я попробовал ее на карте земного покрова с> 15-миллиметровыми ячейками сетки при 300 х 300 м.

Я предполагаю, что описанный выше подход также будет работать, если вы захотите объединить файл многоугольника с неправильными формами с растровыми данными. Возможно, можно объединить шаги 1 и 2, напрямую растеризовав файл полигона в сетку 300xx300, используя растеризацию (растр, вероятно, медленный) или gdal_rasterize. В моем случае мне также нужно было перепроектировать, поэтому я использовал gdalwarp для перепроектирования и дезагрегации одновременно.

Мне приходится сталкиваться с этой же проблемой, чтобы извлечь значения внутри многоугольника из большой мозаики (50k x 50k). У моего многоугольника всего 4 очка. Самый быстрый метод, который я нашел, состоит в том, чтобы cropналожить мозаику на границу многоугольника, триангулировать многоугольник на 2 треугольника, а затем проверить наличие точек в треугольнике (самый быстрый алгоритм, который я нашел). Сравните с extractфункцией, время работы уменьшено с 20 с до 0,5 с. Однако функция cropвсе еще требует около 2 с для каждого многоугольника.

Извините, я не могу предоставить полный воспроизводимый пример. Коды R ниже не включают входные файлы.

Этот метод работает только для простых полигонов.

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