寻找山峰

Question

所以我将这 2 个栅格合并为一个包含高程值的 dem 栅格：

dem1 = read_stars("srtm_43_06.tif")
dem2 = read_stars("srtm_44_06.tif")
pol = st_read("israel_borders.shp")
dem = st_mosaic(dem1, dem2)
dem = dem[, 5687:6287, 2348:2948]
names(dem) = "elevation"
dem = st_warp(src = dem, crs = 32636, method = "near", cellsize = 90)

现在我需要通过找到图像中具有最高海拔的像素的质心来计算山峰的点几何，有人知道我可以使用什么函数吗？

Answer 1

让我们使用 terra，但是 raster 包也可以应用类似的方法。出于测试目的，我们将使用 terra 包

提供的光栅

library(terra)
#> terra 1.5.12
f <- system.file("ex/elev.tif", package="terra")
v <- rast(f)
plot(v)

您只需输入栅格对象名称并按回车键即可检查栅格的详细信息，您可以使用 minmax() 函数形式 [=16] 检查 min 和 max 值=]:

minmax(v)
#>      elevation
#> [1,]       141
#> [2,]       547

让我们通过复制原始光栅来创建另一个光栅，但是检查该值是否为 max 高程值：

w <- v == minmax(v)[2]
plot(w)

让我们创建一个替换矩阵，并将所有 FALSE 替换为 NA，将 TRUE 替换为 1:

mx <- matrix(c(FALSE, NA, TRUE, 1), ncol = 2, byrow = TRUE)
w <- classify(w, mx)

plot(v)
plot(as.polygons(w), add=TRUE)

让我们找到这些多边形的质心：

pts <- centroids(as.polygons(w))
plot(pts, add=TRUE)

让我们看看我们的坐标：

as.data.frame(pts, geom = "WKT")
#>   elevation                   geometry
#> 1         1 POINT (6.020833 50.179167)

^{由 reprex package (v2.0.1)}

创建于 2022-01-29

Answer 2

以 Grzegorz Sapijaszko 的示例为基础，这是通往山顶的替代路径。

library(terra)
f <- system.file("ex/elev.tif", package="terra")
x <- rast(f)

如果有单个最大值，可以做

g <- global(x, which.max)
xyFromCell(x, g[,1])
#            x        y
#[1,] 6.020833 50.17917

现在，考虑具有多个最大值的情况。我再添加三个具有最大值的单元格。

x[c(1000, 2500, 5000)] <- 547

我们可以找到四个最高峰：

g <- global(x, which.max)[[1]]
v <- x[g] |> unlist()
y <- ifel(x == v, v, NA)
p <- as.points(y)
crds(p)
#[1,] 6.020833 50.17917
#[2,] 6.154167 50.10417
#[3,] 5.987500 49.97083
#[4,] 6.237500 49.75417

您应该首先不扭曲（project with terra）栅格数据，因为这会改变像元值并可能改变最高峰的位置。你应该找到原始数据的峰值，然后你可以像这样转换结果。

pp <- project(p, "EPSG:32636") 
crds(pp)
#            x       y
#[1,] -1411008 5916157
#[2,] -1404896 5904422
#[3,] -1422145 5894509
#[4,] -1413735 5864236

对于您的文件，您可以从类似

的内容开始

ff <- c("srtm_43_06.tif", "srtm_44_06.tif")
v <- vrt(ff)
g <- global(x, which.max)

然后继续上面的例子。

寻找山峰

Finding the peak of a mountain

r

point

raster