我有一个多边形列表,我想用它来对terra::rast
砖进行子集化。为此,我使用terra::crop
with lapply
,如下所示,而且速度相当慢。是否有一种矢量化的方式来使用多边形而不是lapply
通过多边形进行子集?
例子
首先,我加载库并创建一个rast
对象。
# Load library
library(terra)
library(geohashTools)
# Create raster
r <- rast(matrix(runif(360 * 2 * 180 *2), ncol = 360 * 2))
# Set extent
ext(r) <- c(-180, 180, -90, 90)
# Examine raster
r
## class : SpatRaster
## dimensions : 360, 720, 1 (nrow, ncol, nlyr)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. :
## source : memory
## name : lyr.1
## min value : 7.853378e-07
## max value : 0.9999981
接下来,我创建了这些对象的砖块。
# Create a brick
b <- c(r, r, r, r, r)
## class : SpatRaster
## dimensions : 360, 720, 5 (nrow, ncol, nlyr)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. :
## sources : memory
## memory
## memory
## ... and 2 more source(s)
## names : lyr.1, lyr.1, lyr.1, lyr.1, lyr.1
## min values : 7.853378e-07, 7.853378e-07, 7.853378e-07, 7.853378e-07, 7.853378e-07
## max values : 0.9999981, 0.9999981, 0.9999981, 0.9999981, 0.9999981
在这里,我只是在创建大量空间多边形。
# All possible coordinates
coords <- expand.grid(-180:180, -90:89)
# Get all unique geohashes for raster
geohashes <- unique(gh_encode(coords$Var2, coords$Var1, precision = 4L))
# Convert to spatial polygons
sp <- geohashTools::gh_to_sp(geohashes)
最后,我遍历每个多边形并用它来裁剪我的砖块。
# Crop raster using geohash polygon
b_cropped <- lapply(seq_along(sp), function(x) terra::crop(b, sp[x]))