如何在R语言terra包中通过行/列索引提取栅格数据子集？

你之前遇到的两个问题原因如下：

1. 使用crop报错是因为terra::ext()默认接收的是地理坐标范围（xmin/xmax/ymin/ymax），你传入的行列号和栅格实际的地理范围不匹配，所以触发「范围不重叠」报错。
2. 方括号索引仅返回数值是因为默认开启了drop=TRUE参数，会自动将栅格子集转为矩阵/向量，关闭该参数即可保留栅格对象格式。

最简实现方案 #

直接在方括号索引中添加drop=FALSE参数，一行代码即可完成操作：

library(terra)
EXT <- c( -108, -105, 39, 42 )
R <- rast( extent=EXT, ncol=14, nrow=14, crs="epsg:4326" )
R[] <- 1:ncell(R)

# 提取第3-12行、第3-12列，保留栅格格式
R_sub <- R[3:12, 3:12, drop = FALSE]

验证结果 #

# 查看输出栅格的行列数
nrow(R_sub)
# [1] 10
ncol(R_sub)
# [1] 10

# 查看输出栅格的地理范围和CRS，和原栅格逻辑一致
crs(R_sub) == crs(R)
# [1] TRUE

可选方案：转地理范围后crop #

如果更习惯使用crop函数，可以先将目标行列号转为对应的地理坐标范围再裁剪：

target_r <- 3:12
target_c <- 3:12

# 构建目标范围对应的地理extent
sub_ext <- ext(
  xyFromCell(R, cellFromRowCol(R, min(target_r), min(target_c))),
  xyFromCell(R, cellFromRowCol(R, max(target_r), max(target_c)))
)

R_sub <- crop(R, sub_ext)

内容的提问来源于stack exchange，提问作者Jim Worrall