Finding nearest altitude point in a DEM in R

Question

I'm looking for a function that would allow me to locate the closest distance point of a given altitude from a specified point on a DEM in R. For instance I have a site on the interior Tibetan plateau located at 4500 meters above sea level and I would like to know at what distance the closest 2500 meter above sea level point is located. I'm essentially asking to find the least cost path but I dont know where the point is that I would be asking to go to.

I've checked out the topodistance package, however, this presumes that I already know where the point is and only then allows me to calculate the distance. I've also checked out the leastcostpath package, however, there is an assumption that I must already know the location of the point to be travelled to. Is there a function that will allow me to do an automated search for the nearest low altitude point?

Answer 1

I built a relatively simple function ('NearAlt') that can perform this task. It is possible to increase its complexity by considering the slope, path cost, as well as the presence of obstacles. However, if I understood correctly, you already have another solution to take these factors into consideration.

# 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)