基于局部最大值绘制边界线——terra包森林砍伐分析技术问询

你已经用terra::distance()计算出森林像素到林缘的距离栅格，接下来可以通过以下步骤实现脊线提取、多边形生成与像素统计：

1. 优化距离栅格计算 #

先简化现有距离计算代码，避免冗余操作：

library(terra)

# 构建示例干扰图层
disturbance <- c(10,250,300,301,302,336,400,500,600)
r <- rast(ncols=36, nrows=18, crs="+proj=utm +zone=1 +datum=WGS84")
r[disturbance[1:2]] <- 1
r[disturbance[3:6]] <- 2
r[disturbance[7:9]] <- 3

# 计算森林像素到干扰区域（林缘）的距离，仅保留森林区域的距离值
dist <- distance(is.na(r)) 
dist <- ifel(is.na(r), NA, dist)

# 可视化原始干扰图层与距离栅格
plot(r, col = c("red", "blue", "purple"), main = "森林干扰图层")
plot(dist, add = TRUE, alpha = 0.5)

2. 提取局部最大值（脊线） #

把距离栅格视为高程数据，用focal()函数识别局部最大值点（即距林缘最远的脊线区域）：

# 定义3x3移动窗口，筛选每个窗口内的中心最大值点
focal_window <- matrix(1, nrow=3, ncol=3)
local_max <- focal(dist, w = focal_window, fun = function(x) {
  if (is.na(x[5])) return(NA)
  if (x[5] == max(x, na.rm = TRUE)) return(x[5]) else return(NA)
})

# 可视化脊线（局部最大值点）
plot(dist, main = "森林距离栅格与脊线")
plot(local_max, add = TRUE, col = "yellow", legend = FALSE)

3. 生成脊线多边形并统计像素 #

提供两种贴合需求的实现方式：

方式一：基于等高线转换多边形 #

直接将距离栅格的等高线转换为多边形，可选择脊线附近的等高线层级：

# 提取等高线并转换为多边形
contour_poly <- as.polygons(as.contour(dist), dissolve = TRUE)

# 可视化等高线多边形
plot(dist, main = "距离栅格与等高线多边形")
plot(contour_poly, add = TRUE, border = "black", lwd = 1)

# 统计每个多边形内的有效像素数量
pixel_count <- extract(dist, contour_poly, fun = function(x) length(na.omit(x)))
contour_poly$pixel_count <- pixel_count[,2]
print(contour_poly[, "pixel_count"])

方式二：基于分水岭算法划分区域（更贴合脊线边界） #

将局部最大值作为种子点，用分水岭算法把距离栅格划分为以脊线为边界的独立区域：

# 将局部最大值转换为矢量点作为分水岭种子
max_points <- as.points(local_max)
max_points$id <- 1:nrow(max_points)

# 执行分水岭算法，生成区域多边形
watershed_poly <- watershed(dist, max_points)
watershed_poly <- as.polygons(watershed_poly, dissolve = TRUE)

# 可视化分水岭区域
plot(dist, main = "距离栅格与脊线划分区域")
plot(watershed_poly, add = TRUE, border = "red", lwd = 1)

# 统计每个区域内的有效像素数量
region_pixel_count <- extract(dist, watershed_poly, fun = function(x) length(na.omit(x)))
watershed_poly$pixel_count <- region_pixel_count[,2]
print(watershed_poly[, "pixel_count"])

两种方式中，分水岭算法更贴合“以脊线为边界划分森林区域”的需求，能精准对应每个脊线辐射的独立森林区块。

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