根据局部最高值绘制边界线

Question

我正在使用 terra 来探索森林砍伐地区。我有两个栅格，第一个栅格包含按树木覆盖年份划分的森林砍伐面积。第二个栅格是剩余的树木覆盖率。我目前正在使用函数

terra::distance()

来计算森林的任何像素距边缘的距离。

这里是一个例子，r代表森林扰动层，dist代表距离边缘的森林冠层距离。

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
plot(r, col =c("red", "blue", "purple"))

dist <- rast(ncols=36, nrows=18, crs="+proj=utm +zone=1 +datum=WGS84")
dist[disturbance] <- 1
dist <- distance(r) 
dist <- ifel(dist == 0, NA, dist)
plot(dist, add=T)

我认为下一步是将距离函数的结果视为高程，并将最大值（即距离森林边缘最远）作为山脊线。我希望能够创建沿着局部最高像素值运行的多边形，类似于

plot(as.contour(dist), add = TRUE)

。这样我就能够量化落在每个边界内的所有光栅像素。

Answer 1

一种方法如下：

首先使用距离函数计算距森林边缘的距离。然后使用

as.contour

从距离栅格创建等高线，并使用

as.polygons

将等高线转换为多边形。最后将栅格 r 转换为点以方便空间操作，并使用 intersect 函数查找每个多边形内的点并聚合这些点以计算每个多边形内有多少个点。

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
plot(r, col=c("red", "blue", "purple"), main="Forest Disturbance")

dist <- distance(r) 
dist <- ifel(dist == 0, NA, dist)
plot(dist, main="Distance from Forest Edge")

contour_lines <- as.contour(dist, maxpixels=100000)
plot(contour_lines, add=TRUE, col="black")
contour_polygons <- as.polygons(contour_lines)
plot(contour_polygons, add=TRUE, col=rainbow(length(unique(contour_polygons$id))), border="black")


r_points <- as.points(r)

intersected_points <- intersect(r_points, contour_polygons)

polygon_summary <- aggregate(intersected_points, by=contour_polygons, FUN=length)

print(polygon_summary)

这将返回以下图：

和