如何使用移动 window 在 R 中迭代裁剪栅格？

Question

我想使用 bbox 或已知范围裁剪栅格，即行和列各 10 个像素。

您可以在下面看到一个可重现的示例：

library(terra)
r <- terra::set.ext(rast(volcano), terra::ext(0, nrow(volcano), 0, ncol(volcano)))
plot(r)

xmin <- xmin(r)
ymax <- ymax(r)
rs <- res(r)

n <- 10 # i.e. of rows and cols size
xmax <- xmin + n * rs[1] 
ymin <- ymax - n * rs[2]

x <- crop(r, c(xmin, xmax, ymin, ymax))
plot(x)

预期的循环是：

遍历所有光栅（r）长度裁剪并将每个光栅片段临时保存到data.frame（或data.table，光栅，spatraster，列表）

Answer 1

使用 terra 数据，按照你的方法，我喜欢：

library(terra)
r <- terra::set.ext(rast(volcano), terra::ext(0, nrow(volcano), 0, ncol(volcano)))

然后我们需要具有 4 个值 (xmin,xmax,ymin,ymax) 的范围，[对于索引 {i,j,k,l}} 然后我们想要从左上角到右上角遍历栅格，然后再次下降，直到我们到达右下角 [for indices {o,p}]，然后进一步裁剪 [for index {q}]。这是 for 循环的相当深的嵌套，为了保护精神...最好做一些助手。

x_mtx <- matrix(c(seq(0,70,10), seq(10,80,10)), ncol = 2)
y_mtx <- matrix(c(seq(51, 1, -10), seq(61, 11, -10)), ncol = 2)

这些消除索引i:l，然后我们使我们的范围：

crop_exts <- vector(mode = 'list', length= 48L)
for(p in 1:6){
  for(o in 1:8){
    crop_exts[[p]][[o]] <- ext(x_mtx[o, 1], x_mtx[o, 2], y_mtx[p, 1], y_mtx[p, 2])
  }
}
# check our work
crop_exts[[1]][[1]] |> as.vector()
#xmin xmax ymin ymax 
#0   10   51   61 
crop_exts[[6]][[8]] |> as.vector()
#xmin xmax ymin ymax 
#  70   80    1   11

看起来是从左上到右下。并注意到 [[1]][[1]]、[[6]][[8]]，我们有第二个嵌套的 {p, o} 循环，因为我们在此处进行索引以提取作物：

cropped_r <- vector(mode='list', length = 48L)# edit, establish outside `for` 
for(p in 1:6){
 for(o in 1:8){
   cropped_r[[p]][[o]] <- crop(r, crop_exts[[p]][[o]])
 }
}
plot(cropped_r[[1]][[1]])
plot(cropped_r[[6]][[8]])

总的来说，可能更容易的是生成完整的范围矩阵，在这种情况下为 dim(48,4)，但我无法绕过它，记住索引，并保持我的精神.. .

Answer 2

了解更多关于为什么要这样做的背景信息会很有用。如果您希望裁剪区域重叠，您的问题也不清楚；我假设你不想要那个。

你可以这样：

library(terra)
r <- rast(volcano)
n <- 10

获取感兴趣的起始细胞

rows <- seq(1, nrow(r), by=n)
cols <- seq(1, ncol(r), by=n)    
cells <- cellFromRowColCombine(r, rows, cols)

获取坐标

# upper-left coordinates of the starting cells 
xy <- xyFromCell(r, cells)
rs <- res(r)
xy[,1] <- xy[,1] - rs[1]/2
xy[,2] <- xy[,2] + rs[2]/2

# add the lower-right coordinates of the end cell
xy <- cbind(xy[,1], xy[,1] + n*rs[1], xy[,2] - n*rs[2], xy[,2])

然后循环

x <- lapply(1:nrow(xy), function(i) {
         crop(r, xy[i,])
       })

验证

e <- lapply(x, \(i) ext(i) |> as.polygons()) |> vect()
plot(r)
lines(e, col="blue", lwd=2)

sapply(x, dim) |> t() |> head()
#     [,1] [,2] [,3]
#[1,]   10   10    1
#[2,]   10   10    1
#[3,]   10   10    1
#[4,]   10   10    1
#[5,]   10   10    1
#[6,]   10   10    1

或者使用基于起始编号和 end-cell 编号的替代方法（为此，您需要 terra 1.5-25，目前是您可以使用 install.packages('terra', repos='https://rspatial.r-universe.dev') 安装的开发版本）

srows <- seq(1, nrow(r), by=n)
scols <- seq(1, ncol(r), by=n)
erows <- pmin(nrow(r), srows+n-1)
ecols <- pmin(ncol(r), scols+n-1)
scell <- cellFromRowColCombine(r, srows, scols)
ecell <- cellFromRowColCombine(r, erows, ecols)
cells <- cbind(scell, ecell)

x <- lapply(1:nrow(cells), function(i) {
        e <- ext(r, cells[i,])
        crop(r, e)
    })

如何使用移动 window 在 R 中迭代裁剪栅格？

How to iterate crop raster in R using a moving window?

r

for-loop

raster

r-raster

terra