CRS为美制英尺、SpatRaster分辨率10m时创建圆形公顷焦点矩阵

需求：通过focal()函数计算大片林区中圆形公顷范围内有林非零二进制像素的总和，所用SpatRaster分辨率为32.8084英尺（即10米），且包含NA值。

数据准备代码 #

library(terra)

# 栅格尺寸与CRS设置
w <- round(2210 * 32.8084)# 还原尺寸与分辨率
h <- round(1829 * 32.8084)
ext <- c(0, w, 0, h)
crs <- "EPSG:2261" # 英尺单位UTM坐标系
res <- 32.8084 # 10米对应的英尺数
num_rows <- ceiling(h / res)
num_cols <- ceiling(w / res)

# 创建栅格
forest <- rast(nrows = num_rows, ncols = num_cols, ext = ext, crs = crs, res = res)

# 生成0-41的随机高度值
set.seed(123)
forest[] <- runif(ncell(forest), min = 0, max = 41)

# 手动设置森林/非森林区块
forest[100:500, 100:500] <- 0
forest[200:750, 500:900] <- 0
forest[800:900, 1200:1700] <- 0
forest[1400:1500, 1700:1800] <- 0
forest[400:900, 1000:1200] <- 41 # 高植被区域
forest[1200:1300, 1700:1800] <- 41
forest[1600:1700, 2000:2100] <- 41
forest[300:870, 1600:1800] <- 41

# 添加NA值模拟 vernal pools
num_cells <- ncell(forest)
sample_indices <- sample(num_cells, 2500, replace = TRUE)
forest[sample_indices] <- NA

创建森林/非森林二进制栅格 #

threshold <- 4 # 阈值：低于4米的灌木视为非森林
forest <- forest > threshold 
forest <- as.numeric(forest)
forest[forest == 1] <- 100 # 1个像素对应100平方米（10m×10m）
plot(forest)

核心计算代码（当前存在问题） #

# 创建临时栅格用于生成权重矩阵
r <- rast(ncols =12, nrow =12, xmin = 0, crs = crs(forest)) 
# 计算公顷圆形的半径（转换为英尺）
radius <- sqrt(10000 / pi) * 3.28084 # 约56.41897米，对应英尺数
# 生成圆形权重矩阵
fw <- focalMat(r, radius, type = "circle")
# 执行focal计算
forest_ha <- terra::focal(x = forest, w = fw, fun = sum, na.policy="omit", na.rm = TRUE, expand =TRUE)
plot(forest_ha)

当前问题与疑问 #

• 输出范围仅为0到100，并非预期的圆形范围内的像素总和，如何验证计算正确性？
• 应该选择12×12的窗口矩阵，还是将栅器重采样至11.28米后使用10×10矩阵，以获得更准确的每公顷计算结果？
• 结合数据的CRS和分辨率，是否遗漏了关键的初始步骤？
• 当前输出刻度不合理，能否直接乘以100将其缩放到10000平方米？

问题分析与解决方案 #

1. 输出范围异常的原因 #

你将森林像素赋值为100（代表100㎡），但focalMat生成的权重矩阵是0-1的加权掩码：完全在圆内的像素权重为1，边缘像素会根据圆形覆盖比例生成0-1之间的权重值。由于你选用的12×12窗口半宽仅为5.5个像素（对应55米），刚好略小于56.42米的半径，导致大部分圆形范围未被窗口覆盖，最终加权求和后最大值仅为100（单个全权重像素的贡献）。

2. 验证计算正确性的方法 #

• 单点验证：找一个完全被森林覆盖的区块（比如你设置的全41的区域），手动计算该点圆形范围内的加权像素数：
  1. 确认圆形公顷半径约为56.42米，对应5.64个10米像素；
  2. 统计窗口内所有像素的权重之和，乘以100㎡，对比forest_ha的输出值；
  3. 如果该点周围全是森林，理论值应接近10000㎡，实际输出会因边缘像素的加权略低。
• 对比空白区域：找完全非森林的区域，输出应接近0，以此验证权重矩阵的有效性。

3. 窗口尺寸与重采样的选择 #

• 优先调整窗口尺寸：12×12窗口的覆盖范围不足，建议改用13×13窗口（半宽6个像素，对应60米），完全覆盖56.42米的半径，确保圆形范围被完整包含。修改代码中临时栅格的尺寸：

  r <- rast(ncols =13, nrow =13, xmin = 0, crs = crs(forest)) 
  

• 不建议盲目重采样：将10米分辨率重采样到11.28米会损失原始数据精度，除非你有更高分辨率的原始数据。调整窗口尺寸是更高效准确的方案。

4. 初始步骤的遗漏点 #

• CRS单位处理正确：你已经将半径转换为英尺（匹配EPSG:2261的单位），这部分没问题；
• 二进制栅格赋值可优化：建议用1代表森林像素，最后再乘以100㎡，避免计算混淆：

  forest <- as.numeric(forest > threshold)
  # 计算后再转换为平方米
  forest_ha <- focal(forest, w=fw, fun=sum, na.rm=TRUE) * 100
  

5. 输出刻度的调整 #

不能直接乘以100——你当前的输出已经是加权像素数×100，但因为窗口尺寸不足导致最大值仅为100。调整窗口尺寸后，全森林覆盖区域的输出会接近10000㎡，无需额外缩放。

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