如何选取多边形边界向内指定距离的栅格单元并计算其平均值？

需求说明 #

需要从多边形边界向内选取指定距离（如1公里、5公里）范围内的栅格单元，最终计算这些单元的平均值，而非对整个多边形覆盖的栅格求平均。

实现思路 #

• 创建原始多边形的向内负缓冲区（缓冲区距离设为负的指定值）
• 用原始多边形减去该负缓冲区，得到仅包含边界指定距离范围的新多边形
• 用新多边形对目标栅格进行掩膜和裁剪，筛选出目标区域的栅格单元
• 计算裁剪后栅格单元的平均值

完整代码示例 #

# 加载所需包
library(raster)
library(rgeos)

# 创建示例栅格并赋予随机值（替代全0值，便于验证计算结果）
r <- raster(xmn=1035792, xmx=1116792, ymn=825303.6, ymx=937803.6, resolution=12.5, crs="+init=epsg:3174")
r <- setValues(r, runif(ncell(r), min=0, max=100)) # 赋予0-100的随机数值

# 创建示例多边形
x <- c(1199999, 1080000, 1093067, 1090190, 1087977, 1070419, 1180419)
y <- c(957803.6,937803.6, 894366.9, 872153.9, 853703.0, 825353.6, 805353.6)
poly.lake <- SpatialPolygons(list(Polygons(list(Polygon(data.frame(x,y))), ID = 1)))
crs(poly.lake) <- crs(r) # 确保多边形与栅格坐标系一致

# 初始栅格掩膜与裁剪（保留多边形内的栅格）
r_masked <- mask(r, poly.lake)
r_cropped <- crop(r_masked, poly.lake)

# --------------------------核心步骤--------------------------
# 1. 设置向内缓冲距离（单位：米，此处以1公里为例）
buffer_distance <- -1000 # 负数值表示向内缓冲

# 2. 生成向内负缓冲区
poly_buffer <- gBuffer(poly.lake, width = buffer_distance, byid = TRUE)

# 3. 原始多边形减去缓冲区，得到仅含边界范围的多边形
poly_boundary <- gDifference(poly.lake, poly_buffer, byid = TRUE)

# 4. 用边界多边形筛选栅格
r_boundary <- mask(r_cropped, poly_boundary)
r_boundary_cropped <- crop(r_boundary, poly_boundary)

# 5. 计算边界范围内栅格的平均值
mean_value <- cellStats(r_boundary_cropped, stat='mean', na.rm=TRUE)
cat("边界向内1公里范围内栅格的平均值：", mean_value, "\n")
# -----------------------------------------------------------

# 可视化结果（可选）
par(mfrow=c(1,2))
plot(poly.lake, main="原始多边形与边界范围")
plot(poly_buffer, add=TRUE, col="lightblue")
plot(poly_boundary, add=TRUE, col="orange")
legend("topleft", legend=c("原始多边形", "向内缓冲区", "边界范围"), fill=c("transparent", "lightblue", "orange"))

plot(r_boundary_cropped, main="边界范围内的栅格")
plot(poly_boundary, add=TRUE, border="red")

示例配图说明 #

示例配图1：展示了原始多边形与被其掩膜后的栅格的叠加效果
示例配图2：展示了仅保留多边形边界向内指定距离范围的栅格区域示意

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