如何使terra::focal()支持偶数边长权重矩阵（含Roberts交叉检测）

问题背景 #

需要用terra::focal()实现Roberts交叉边缘检测，但focal()仅支持奇数边长的窗口。尝试用3×3矩阵填充NA模拟偶数窗口的方法无效，测试代码如下：

library(terra)
r <- rast(matrix(rnorm(10000), ncol=100, nrow=100))
# 尝试用NA填充的3×3矩阵模拟Roberts窗口
robx <- matrix(c(NA,1,0,NA,0,1),nrow=3)
plot(focal(r,w=robx,fun="sum")) # 无法得到预期结果

而不含NA的3×3矩阵（如Sobel算子）可以正常运行：

sobx <- matrix(c(-1,-2,-1,0,0,0,1,2,1) / 4, nrow=3)
plot(focal(r,w=sobx,fun="sum"))

解决方案 #

Roberts交叉本质是计算相邻像素对的差值，不需要依赖偶数窗口的focal()，可以通过以下两种方法实现：

方法1：使用shift()直接计算像素差值 #

这是最直接的方式，通过位移栅格获取目标像素对，再计算差值：

library(terra)

# 创建测试栅格
r <- rast(matrix(rnorm(10000), ncol=100, nrow=100))

# Roberts交叉的两个方向计算
# 方向1：当前像素 - 右下方像素
rob1 <- r - shift(r, dx=1, dy=1)
# 方向2：下方像素 - 右方像素
rob2 <- shift(r, dx=0, dy=1) - shift(r, dx=1, dy=0)

# 合并两个方向的边缘结果（平方和开根号）
roberts_edge <- sqrt(rob1^2 + rob2^2)

# 可视化结果
plot(roberts_edge)

方法2：调整权重矩阵适配focal() #

如果一定要用focal()，可以构造以当前像素为中心的3×3权重矩阵，通过NA和权重值定义有效像素对，并开启na.rm=TRUE忽略边缘NA：

library(terra)

r <- rast(matrix(rnorm(10000), ncol=100, nrow=100))

# 定义Roberts两个方向的3×3权重矩阵
rob_mat1 <- matrix(c(0, 0, 0,
                     0, 1, 0,
                     0, 0, -1), nrow=3)
rob_mat2 <- matrix(c(0, 0, 0,
                     0, 0, 1,
                     0, -1, 0), nrow=3)

# 用focal计算两个方向的差值
rob1 <- focal(r, w=rob_mat1, fun="sum", na.rm=TRUE)
rob2 <- focal(r, w=rob_mat2, fun="sum", na.rm=TRUE)

# 合并结果
roberts_edge <- sqrt(rob1^2 + rob2^2)

plot(roberts_edge)

说明 #

• 原方法失败的原因：focal()的窗口以当前像素为中心，填充NA的3×3矩阵无法正确对应Roberts交叉需要的像素对位置，且未设置na.rm=TRUE导致计算被NA干扰。
• 两种方法中，shift()的方式更高效直观，适合这类基于像素对差值的边缘检测算法。

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