在 R 中的 DEM 中查找最近的高度点

Question

我正在寻找一个函数，允许我从 R 中 DEM 上的指定点定位给定高度的最近距离点。例如，我在青藏高原内陆有一个地点，海拔 4500 米我想知道最近的海拔2500米点的距离是多少。我本质上是要求找到成本最低的路径，但我不知道我要求去的重点在哪里。

我已经检查了 topodistance 包，但是，这假设我已经知道该点在哪里，然后才允许我计算距离。我还检查了leastcostpath包，但是，假设我必须已经知道要前往的点的位置。有没有一个功能可以让我自动搜索最近的低海拔点？

Answer 1

我构建了一个相对简单的函数（'NearAlt'）可以执行此任务。通过考虑坡度、路径成本以及障碍物的存在，可以增加其复杂性。但是，如果我理解正确的话，您已经有了另一种解决方案来考虑这些因素。

# install.packages("raster")
# install.packages("rgdal")
# install.packages("sf")

library(raster)
library(rgdal)
library(sf)


dem <- raster("n28_e082_1arc_v3.tif") # Load your DEM

NearAlt <- function(x, CoordStart, targetAlt) {
  RastDiff <- abs(x - targetAlt) # Difference to target altitude
  
  # Extract the value at the start coordinate
  startDiff <- extract(RastDiff, matrix(c(CoordStart[1], CoordStart[2]), ncol=2))
  
  # Check if is already at the desired altitude
  if (startDiff == 0) {
    return(list(distance = 0, coord = CoordStart))
  }
  
  # Convert the raster to a data frame
  diffMatrix <- as.data.frame(RastDiff, xy=TRUE)
  
  # Find the minimum difference
  cells <- diffMatrix[diffMatrix[3] == min(diffMatrix[3], na.rm = TRUE),]

  # Calculate distances and find the closest
  point_orign <- st_as_sf(data.frame(x = c(CoordStart[1]), y = c(CoordStart[2])), coords = c("x", "y"), crs = crs(x))
  points_destiny <- st_as_sf(data.frame(x = cells[,1], y = cells[,2]), coords = c("x", "y"), crs = crs(x))
  distances <- st_distance(point_orign, points_destiny)
  minDistanceIndex <- which.min(as.integer(distances))
  
  return(list(distance = distances[minDistanceIndex], coord = c(cells[minDistanceIndex,1], cells[minDistanceIndex,2])))
}


CoordStart <- c(82.5367721, 28.8306240) # Replace with your coordinates
targetAlt <- 2500 # Replace with your target altitude

result <- NearAlt(x = dem, CoordStart, targetAlt)
print(paste("Closest distance:", result$distance))
print(paste("Coordinates:", result$coord[1], ",", result$coord[2]))

# Check the reslts
extract_start <- extract(dem, SpatialPoints(matrix(CoordStart, ncol = 2)), sp = T)
extract_destiny <- extract(dem, SpatialPoints(matrix(result$coord, ncol = 2)), sp = T)
alt_start <- extract_start@data$n28_e082_1arc_v3 # Change to match your DEM
alt_destiny <- extract_destiny@data$n28_e082_1arc_v3 # Change to match your DEM

# Tabulate the results
results_df <- data.frame("X" = c(CoordStart[1], result$coord[1]), "Y" = c(CoordStart[2], result$coord[2]), "Altitud"=c(alt_start, alt_destiny), "Position"= c("start", "destiny"))
results_df_points <- SpatialPointsDataFrame(results_df[,1:2], data = results_df, proj4string = crs(dem))

# Write as shapefile
writeOGR(obj=results_df_points, dsn="./results_df_points.shp", layer="torn", driver="ESRI Shapefile")

# Visualize the results
plot(dem)
plot(results_df_points, col="blue", add=T)