在R语言中从DEM指定点查找目标海拔最近点的方法问询

核心思路 #

你需要的是基于成本表面的距离变换——先把所有符合目标海拔的栅格单元标记为"目标区域"，然后计算整个区域内每个点到最近目标点的最低成本距离，最后直接提取起点位置的距离值即可。这种方法不需要预先知道具体目标点，比批量计算到每个目标点的路径效率高得多。

具体实现（用gdistance包） #

以下是针对青藏高原DEM场景优化的可运行代码：

1. 加载依赖包

library(gdistance)
library(raster)
library(sf)

2. 读取DEM和起点数据

假设你的DEM是tif格式，起点为sf点对象：

# 读取DEM
dem <- raster("your_dem.tif")

# 定义起点（替换为你的实际点位坐标）
start_point <- st_sfc(st_point(c(90, 30)), crs = st_crs(dem))

3. 标记目标海拔区域

考虑DEM精度误差，筛选接近目标海拔的栅格单元：

target_elev <- 2500
# 允许±5米误差，避免漏选
target_mask <- dem >= (target_elev - 5) & dem <= (target_elev + 5)

4. 创建成本表面

以坡度成本为例（可自定义其他成本函数，比如海拔落差成本）：

# 计算坡度（单位：度）
slope <- terrain(dem, opt = "slope", unit = "degrees")

# 转换为成本系数：坡度越大，移动成本越高
cost <- 1 + (slope / 90) ^ 2

# 创建栅格间的过渡矩阵（处理8方向移动成本）
transition_matrix <- transition(cost, function(x) 1/mean(x), 8)
transition_matrix <- geoCorrection(transition_matrix, type = "c")

5. 计算到最近目标点的最低成本距离

# 提取目标区域的栅格坐标
target_coords <- xyFromCell(dem, which(target_mask[] == TRUE))

# 生成距离变换栅格，每个单元值为到最近目标点的最低成本距离
distance_raster <- distanceFromPoints(transition_matrix, target_coords)

6. 提取起点的距离值

min_distance <- extract(distance_raster, start_point)
cat("最近2500米海拔点的最低成本距离：", min_distance, "米\n")

7. （可选）提取并可视化最短路径

# 找到距离起点最近的目标点坐标
closest_target <- target_coords[which.min(extract(distance_raster, target_coords)), ]

# 生成最短路径
shortest_path <- shortestPath(transition_matrix, 
                              fromCoords = st_coordinates(start_point),
                              toCoords = closest_target,
                              output = "SpatialLines")

# 可视化
plot(dem)
plot(shortest_path, add = TRUE, col = "red", lwd = 2)
plot(start_point, add = TRUE, col = "blue", pch = 19)

效率优化建议 #

针对青藏高原大尺度DEM：

• 先裁剪DEM到起点周围的合理范围（比如以起点为中心，半径500公里的区域）
• 若目标海拔点数量极多，可对目标区域做降采样，减少计算量

替代方案（目标点较少时） #

如果符合目标海拔的点数量有限，可使用leastcostpath批量计算：

library(leastcostpath)

# 转换目标坐标为sf点对象
target_points <- st_as_sf(target_coords, coords = c("x", "y"), crs = st_crs(dem))

# 批量生成从起点到所有目标点的成本路径
lcp_list <- create_lcps(
  dem = dem,
  origin = start_point,
  destination = target_points,
  cost_function = "tobler"
)

# 筛选最短路径
lcp_lengths <- sapply(lcp_list, function(x) st_length(x))
min_lcp <- lcp_list[which.min(lcp_lengths)]

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