Мне нужно рассчитать долю каждой растровой ячейки в сетке высокого разрешения (растровый стек с 8 слоями), покрытой многоугольником, используя R.
Мой стандартный подход заключается в использовании raster :: rasterize (..., getCover = TRUE), однако этот подход занимает очень много времени, особенно когда размер многоугольника увеличивается.
В качестве альтернативы я попытался обрезать растровый стек до размера многоугольникапреобразование набора растров в многоугольники и вычисление пропорции от пересечения результирующих фигур с исходным многоугольником.Это хорошо работает для небольших полигонов, но ломается при увеличении полигона, потому что R не хватает памяти (я ограничен 16 ГБ) или потому, что вычисление пересечения занимает слишком много времени.
Здесь воспроизводимый пример, использующий мойтекущее решение с очень маленьким файлом формы.
library(raster)
library(spex)
library(dplyr)
library(sf)
library(data.table)
# setup a dummy example
r <- raster(nrow = 21600, ncol = 43200)
r[] <- 1:933120000
r_stack <- stack(r,r,r,r,r,r,r,r)
# get a small dummy shapefile
shp_small <- raster::getData(name = "GADM", country = "CHE", level = 2, download = TRUE)
shp_small <- st_as_sf(shp_small)[1, ]
# for comparison, use a big dummy shapefile
# shp_big <- raster::getData(name = "GADM", country = "BRA", level = 0, download = TRUE)
## Approach for a small shape file
stack_small <- raster::crop(r_stack, shp_small, snap = "out")
## transform to polygon
stack_small_poly <- spex::polygonize(stack_small)
stack_small_poly$ID <- 1:nrow(stack_small_poly)
## I can then perform the necessary calculations on the polygons to obtain
## the proportional overlay
# transform to mollweide for area calculation
mollw <- "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
stack_small_crs <- st_crs(stack_small_poly)
stack_small <- st_transform(stack_small_poly, mollw)
stack_small_poly <- st_transform(stack_small_poly, mollw)
# calculate area for each cell
stack_small_poly$area_org <- st_area(stack_small_poly)
# transform to world equidistant cylindrical
stack_small_poly<- st_transform(stack_small_poly, 4087)
shp_small <- st_transform(shp_small, 4087)
# get call the cells that intersect with the shape (this might take a while)
i <- st_intersects(stack_small_poly, shp_small, sparse = FALSE)
stack_small_poly <- dplyr::filter(stack_small_poly, i)
# now calculate the extact intersection (this might take a while)
st_agr(stack_small_poly) <- "constant"
stack_small_poly <- st_intersection(stack_small_poly, st_geometry(shp_small))
# calculate the new areas and backtransform
stack_small_poly <- st_transform(stack_small_poly, mollw)
stack_small_poly$new_area <- st_area(stack_small_poly)
stack_small_poly <- st_transform(stack_small_poly, stack_small_crs)
# calculate proportion
stack_small_poly$proportion <- as.numeric(stack_small_poly$new_area/stack_small_poly$area_org)
# finally transform to data.table for subsequent analysis
st_geometry(stack_small_poly) <- NULL
setDT(stack_small_poly)
Я ищу решение на R, способное выполнить задачу за 10-15 минут (желательно быстрее) с ограничением памяти 16 ГБ ОЗУдля шейп-файла, представляющего Бразилию (см. shp_big в коде выше).
Мне хорошо известно, что этот оптимум может быть недостижим, и каждое предложение, приводящее к сокращению времени выполнения и / или использования памяти, более чем выгодно.