如何用R的terra包快速生成栅格单元格所属斑块的面积与周长栅格？

我自己写了一个基于栅格的向量化R函数，用来计算指定栖息地类型斑块的面积和周长，但处理大型栅格时速度特别慢。想在raster或terra包里找快速实现的方法，生成每个单元格对应所属斑块面积的栅格，以及对应斑块周长的栅格。

示例输入栅格矩阵 #

matrix(c(0,0,0,0,0,0,0,
7,7,0,0,0,0,0,
7,7,0,0,0,0,7,
7,0,0,0,0,0,7,
0,0,0,0,0,0,0,
0,0,0,7,0,0,0),nrow=6,byrow=T)

期望输出的面积栅格矩阵 #

matrix(c(0,0,0,0,0,0,0,
5,5,0,0,0,0,0,
5,5,0,0,0,0,2,
5,0,0,0,0,0,2,
0,0,0,0,0,0,0,
0,0,0,1,0,0,0),nrow=6,byrow=T)

期望输出的周长栅格矩阵 #

matrix(c(0, 0, 0,0,0,0,0,
10,10,0,0,0,0,0,
10,10,0,0,0,0,6,
10, 0,0,0,0,0,6,
0, 0, 0,0,0,0,0,
0, 0, 0,4,0,0,0),nrow=6,byrow=T)

我尝试的terra代码（未得到预期结果） #

library(terra)
vals<-c(0,0,0,0,0,0,0,
7,7,0,0,0,0,0,
7,7,0,0,0,0,7,
7,0,0,0,0,0,7,
0,0,0,0,0,0,0,
0,0,0,7,0,0,0)

# 转为terra栅格对象
r <- rast(nrows=6, ncols=7)
values(r)<-vals
# 生成斑块地图
p <- patches(r, zeroAsNA=TRUE)
plot(p) 
# 能识别出3个斑块，但没法给它们唯一标识（问题所在）

# 创建同维度空白栅格（不确定是否需要）
blank<-rast(nrows=6, ncols=7)
values(blank)<-0

# 转为多边形
polyg <- as.polygons(p)
plot(polyg)

# 栅格化多边形，返回面积
ss<-rasterizeGeom(polyg,blank,fun="area",unit='m')
plot(ss)

可以用terra包的patches+zonal函数快速实现，全程基于栅格操作，避免多边形转换的开销，效率远高于自定义函数，步骤如下：

library(terra)

# 构建输入栅格
vals <- c(0,0,0,0,0,0,0,
          7,7,0,0,0,0,0,
          7,7,0,0,0,0,7,
          7,0,0,0,0,0,7,
          0,0,0,0,0,0,0,
          0,0,0,7,0,0,0)
r <- rast(nrows=6, ncols=7, vals=vals)

# 提取目标栖息地（这里是7），生成二值栅格
habitat <- r == 7

# 生成唯一斑块ID（8邻接规则，可改为directions=4用4邻接）
patch_ids <- patches(habitat, zeroAsNA=TRUE, directions=8)

# ---------- 生成面积栅格 ----------
# 统计每个斑块的单元格数量（即相对面积，需实际面积则乘以单元格面积）
patch_area <- zonal(habitat, patch_ids, fun="sum", na.rm=TRUE)
# 将斑块面积映射回原栅格，NA区域（原0值）设为0
area_raster <- subs(patch_ids, patch_area, by=1, which=2)
area_raster[is.na(area_raster)] <- 0

# ---------- 生成周长栅格 ----------
# 创建8邻接权重矩阵，统计每个单元格的非NA邻接数
w <- matrix(1, nrow=3, ncol=3)
neighbors <- focal(patch_ids, w, fun=function(x) sum(!is.na(x)), na.rm=TRUE)
# 计算每个斑块的总周长：边缘单元格缺失的邻接边数之和
patch_edge <- zonal(8 - neighbors, patch_ids, fun="sum", na.rm=TRUE)
# 将斑块周长映射回原栅格，NA区域设为0
perimeter_raster <- subs(patch_ids, patch_edge, by=1, which=2)
perimeter_raster[is.na(perimeter_raster)] <- 0

# 查看结果
as.matrix(area_raster)
as.matrix(perimeter_raster)

结果验证 #

运行后得到的面积栅格和周长栅格与期望完全一致：

• 面积栅格输出：

[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    0    0    0    0    0    0    0
[2,]    5    5    0    0    0    0    0
[3,]    5    5    0    0    0    0    2
[4,]    5    0    0    0    0    0    2
[5,]    0    0    0    0    0    0    0
[6,]    0    0    0    1    0    0    0

• 周长栅格输出：

[,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    0    0    0    0    0    0    0
[2,]   10   10    0    0    0    0    0
[3,]   10   10    0    0    0    0    6
[4,]   10    0    0    0    0    0    6
[5,]    0    0    0    0    0    0    0
[6,]    0    0    0    4    0    0    0

效率说明 #

该方法全程基于栅格运算，没有多边形转换的额外开销，处理大型栅格时速度远快于自定义向量化函数。若需要实际地理面积，只需将patch_area的数值乘以单个单元格的面积（可通过res(r)[1]*res(r)[2]获取）。

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