如何有效地从(许多)大栅格中为许多多边形提取值
how to efficiently extract values from (many) large raster(s) for many polygones
我有一组点,我想从几个大栅格中提取值作为这些点周围的缓冲区。 栅格太大,无法保存在内存中(> 1e10 个像元)。我在下面说明了我目前的方法,但如果有更快的方法我会很感兴趣。
library(maps)
library(sf)
library(raster)
library(dplyr)
library(parallel)
# sf object with polygones for which we want values
crs <- "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
map <- sf::st_as_sf( maps::map(regions = "Sweden", plot = FALSE, fill = TRUE))
map <- st_transform(map, crs = crs)
sf_points <- st_sfc(st_sample(map, 100))
sf_points <-
data.frame(A = 1:length(sf_points)) %>%
st_set_geometry(sf_points)
# raster too large to fit in memory
# the raster(s) I am working on has 10m resolution
r <- raster(extent(map), nrow = 15000, ncol = 7000,
crs = crs)
values(r) <- rep(sample(1:10, 77, replace = T), length.out = ncell(r))
#use the parallel package for parallel processing
cluster <- makeCluster(4)
clusterExport(cluster, c("r","sf_points", "as_Spatial"))
List_points <-
sf_points %>%
mutate(split = rep(1:ceiling(n()/4), each=4, length.out=n())) %>% # 4 cores
split(f = .$split) %>%
parLapply(cl = cluster, X = ., function(x) raster::extract(r, y = as_Spatial(x), buffer = 5000)) %>%
unlist(recursive = F)
我对每个栅格重复提取。由于这些值是有序的,因此我可以汇总列表中的像素值。我不能(轻松地)创建栅格堆栈,因为栅格具有不同的范围。
使用 velox
包已经 sugested 似乎不起作用,因为它试图将光栅加载到失败的内存中。我可以尝试按块加载它,但随后我需要弄清楚哪些点在哪个块上...
也许你可以使用多边形来加快速度而不聚合(溶解)缓冲区
library(raster)
swe <- getData("GADM", country="SWE", level=0)
set.seed(0)
pts <- spsample(swe, 100, "regular")
r <- raster(swe, nrow = 15000, ncol=7000)
values(r) <- rep(sample(1:10, 77, replace = T), length.out = ncell(r))
b1 <- buffer(pts, 5000, dissolve=FALSE)
b2 <- buffer(pts, 5000, dissolve=TRUE)
system.time(e1 <- extract(r, pts, buffer=5000))
# user system elapsed
# 1.39 0.02 1.40
system.time(e2 <- extract(r, b1))
# user system elapsed
# 0.88 0.00 0.88
system.time(e3 <- extract(r, b2))
# user system elapsed
# 26.34 25.02 51.52
显然 b1
比 b2
表现得更好;但并不比第一种方法快多少。
你说你不能制作RasterStack,因为栅格有不同的范围。但是,如果(且仅当!)它们具有相同的原点和分辨率,您可以先将所有区域转换为 xy 坐标,然后使用它们。
像这样:
z <- rasterize(b, r)
pts <- rasterToPoints(z, xy=TRUE)
以上需要时间,但之后
system.time(a <- extract(r, zz[,1:2]))
user system elapsed
0.04 0.00 0.04
每个点平行执行此操作可能会更快,并在光栅化之前使用 crop(raster(r), polygon)
。
我有一组点,我想从几个大栅格中提取值作为这些点周围的缓冲区。 栅格太大,无法保存在内存中(> 1e10 个像元)。我在下面说明了我目前的方法,但如果有更快的方法我会很感兴趣。
library(maps)
library(sf)
library(raster)
library(dplyr)
library(parallel)
# sf object with polygones for which we want values
crs <- "+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
map <- sf::st_as_sf( maps::map(regions = "Sweden", plot = FALSE, fill = TRUE))
map <- st_transform(map, crs = crs)
sf_points <- st_sfc(st_sample(map, 100))
sf_points <-
data.frame(A = 1:length(sf_points)) %>%
st_set_geometry(sf_points)
# raster too large to fit in memory
# the raster(s) I am working on has 10m resolution
r <- raster(extent(map), nrow = 15000, ncol = 7000,
crs = crs)
values(r) <- rep(sample(1:10, 77, replace = T), length.out = ncell(r))
#use the parallel package for parallel processing
cluster <- makeCluster(4)
clusterExport(cluster, c("r","sf_points", "as_Spatial"))
List_points <-
sf_points %>%
mutate(split = rep(1:ceiling(n()/4), each=4, length.out=n())) %>% # 4 cores
split(f = .$split) %>%
parLapply(cl = cluster, X = ., function(x) raster::extract(r, y = as_Spatial(x), buffer = 5000)) %>%
unlist(recursive = F)
我对每个栅格重复提取。由于这些值是有序的,因此我可以汇总列表中的像素值。我不能(轻松地)创建栅格堆栈,因为栅格具有不同的范围。
使用 velox
包已经 sugested
也许你可以使用多边形来加快速度而不聚合(溶解)缓冲区
library(raster)
swe <- getData("GADM", country="SWE", level=0)
set.seed(0)
pts <- spsample(swe, 100, "regular")
r <- raster(swe, nrow = 15000, ncol=7000)
values(r) <- rep(sample(1:10, 77, replace = T), length.out = ncell(r))
b1 <- buffer(pts, 5000, dissolve=FALSE)
b2 <- buffer(pts, 5000, dissolve=TRUE)
system.time(e1 <- extract(r, pts, buffer=5000))
# user system elapsed
# 1.39 0.02 1.40
system.time(e2 <- extract(r, b1))
# user system elapsed
# 0.88 0.00 0.88
system.time(e3 <- extract(r, b2))
# user system elapsed
# 26.34 25.02 51.52
显然 b1
比 b2
表现得更好;但并不比第一种方法快多少。
你说你不能制作RasterStack,因为栅格有不同的范围。但是,如果(且仅当!)它们具有相同的原点和分辨率,您可以先将所有区域转换为 xy 坐标,然后使用它们。
像这样:
z <- rasterize(b, r)
pts <- rasterToPoints(z, xy=TRUE)
以上需要时间,但之后
system.time(a <- extract(r, zz[,1:2]))
user system elapsed
0.04 0.00 0.04
每个点平行执行此操作可能会更快,并在光栅化之前使用 crop(raster(r), polygon)
。