如何逐块栅格化 sf 几何？答案

【问题标题】：How to rasterize sf geometries tile by tile?如何逐块栅格化 sf 几何？
【发布时间】：2020-08-28 15:43:42
【问题描述】：

我正在寻找一种方法来逐步光栅化（或加快）一定范围内每个图块内的几何图形，并将光栅化几何图形的各个部分连接到整个光栅对象。

## create sample data

# create a frame
library(sf)
ob = st_sf(st_sfc(st_polygon(list(rbind(c(0,0), c(0,9), c(6,9), c(6,0), c(0,0)))))
# create tiles
library(GSIF)
tl <- getSpatialTiles(as(ob, 'Spatial'), block.x=3, overlap.percent=2)
plot(tl)
# create sample polygons
g <- st_sfc(st_point(c(1,2)), st_point(c(5,6)), st_point(c(2,4)), st_point(c(3,3)), st_point(c(3,4)), st_point(c(4,5)))
g.b <- st_buffer(g,0.6)
p <- st_sf(value = ceiling(10*runif(6)),
           geometry = st_sfc(g.b))
plot(p, add=TRUE)

我猜，裁剪每个图块内的多边形（或多边形的一部分）效果很好。我的实际目标是逐块处理数据块。看起来下面的循环执行此操作，并将多边形的部分连接到整个多边形。好吧，输出（cr）与输入（p）相同......我认为这一步可能是必要的，以便在下一步中栅格化生成的（部分）多边形。我没想到多边形会被连接起来。我试图根据中间结果构建一个循环，因此添加了以下几行：

# crop polygons (or parts) for each tile
result <- p
for(i in 1:length(tl)) {cr <- rbind(result, st_crop(p, tl[i]))}
cr <- cr[-7,]
library(scales)
plot(st_geometry(cr), col=alpha("white", 0.5), add=TRUE)

但是，我很难栅格化（中间）结果。

【问题讨论】：

标签： r loops spatial sf rasterize

【解决方案1】：

试试这个：

myF <- function(i){return(st_crop(p, tl[i]))}
cr <- do.call(rbind, lapply(1:length(tl), myF))

【讨论】：