使用R的raster包focal函数计算3x3邻域像元差平方均值异常

你当前的代码错误完全是因为对focal函数的运行机制理解有误：

1. 传入focal的自定义函数，入参是当前滑动到的3×3窗口展开的长度为9的向量，且要求函数返回单个标量值作为当前窗口中心像元的计算结果。你编写的sq_dff完全没有使用传入的窗口参数，反而全局计算了整个栅格的所有邻接像元差平方，返回的是全栅格长度的均值向量，focal只会取固定值作为所有窗口的输出，因此所有有值的像元都得到了相同的0.67结果。
2. focal默认不对栅格边缘做填充，边缘像元无法形成完整3×3窗口时会直接返回NA，因此你只会看到中心9个像元有值。
3. 你在函数内开启的多线程完全没有作用，for循环本身是单线程执行，还额外增加了全局遍历的开销，完全违背了用focal提升效率的初衷。

你只需要针对每个3×3窗口做局部计算即可，逻辑非常简单：

library(raster)
# 你的示例栅格
r <- raster(matrix(1:25,nrow=5))
r[] <-c(2,3,2,3,2,
        3,2,3,2,NA,
        NA,3,2,3,2,
        NA,2,3,2,3,
        2,3,2,3,NA)

# 正确的自定义计算函数
sq_dff <- function(w) {
  # 3×3窗口按行展开后的第5位是中心像元值
  center <- w[5]
  # 排除中心自身，取剩余8个邻域像元计算差平方均值
  mean((w[-5] - center)^2, na.rm = TRUE)
}

# 调用focal计算，pad=TRUE表示对边缘做填充，可根据需求决定是否开启
r1 <- focal(x = r, w = matrix(1, 3, 3), fun = sq_dff, pad = TRUE, padValue = NA)

上述代码的计算效率远高于你原有遍历全栅格邻接关系的方案，13万像元的计算耗时会在秒级。如果你的栅格尺寸后续进一步扩大，可以用focalCluster函数直接调用多核心运行，不需要自己管理多线程进程。

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