有识别山脉山脊线的R方法吗？

Question

我正在做一些地形分析，我需要找到从某些点到山脊的距离。我已经使用

geodata

来访问 SRTM 数据，所以我有一个 DEM。

terra

有很棒的内置方法来计算斜率和 TRI。

aspect <- terra::terrain(dem, v = "aspect")
tri <- terra::terrain(dem, v = "TRI")
slope <- terra::terrain(dem, v = "slope")

但是，我没有找到提取脊线的方法。我认为一种方法是找到斜率的变化，但是 1）这里峰和谷都有 0 斜率（即没有负斜率），2）

focal

方法返回局部最大值，而我需要找到最高的山脊。

stacks上的其他答案都集中在寻找最高点，但我需要线条，而不是点值。此外，我使用的地理范围很广，因此为山脊设置阈值海拔是不可行的。

有更多地理经验的人可以提出解决方案吗？

谢谢！

Answer 1

感谢this我能够创建以下代码。请参阅他的帖子以获得进一步的解释。增大

ridge_threshold

可使脊线更粗且颗粒度更小。

simplifyTol=800

表示将合并为一个的多边形或线之间的距离[m]。

library(terra)
library(geodata)
library(tidyverse)
library(sf)
library(tmap)
library(elevatr)
library(whitebox)

get_ridges_from_dem = function(
    dem=0, 
    ridge_threshold=450, # please adjust depending on your dem!
    ridge_color = "red", # rdiges are red
    ridge_lines_width=2,
    write_ridge_lines = False,
    simplifyTol = 200, # simplify tolerance in meters
    usePolygons = T
    ) 
  {
  # set the current script's location as working directory
  setwd(dirname(rstudioapi::getSourceEditorContext()$path))
  
  # Check current working directory
  wd <- getwd()
  
  output_dir <- wd # Directory to save outputs
  corrected_dem_path <- file.path(output_dir, "corrected_dem.tif")
  
  dem_path <- file.path(output_dir, "dem.tif")
  flow_direction_path <- file.path(output_dir, "flow_direction.tif")
  flow_accumulation_path <- file.path(output_dir, "flow_accumulation.tif")
  ridge_output_path <- file.path(output_dir, "ridge_lines.tif")
  
  
  if(dem != 0){
    print("using provided dem!")
  }else{
    print("Fetching a predefined dem!")
    locations <- st_point(c(-105, 40)) %>%
      st_sfc(crs = 4326)  # WGS 84 coordinate system
    
    # Fetch elevation raster using `elevatr`
    dem <- get_elev_raster(locations = locations, z = 10, src = "aws")  # Ensure you have internet access for this step
  }
    
  
  writeRaster(dem, dem_path, format = "GTiff", overwrite = TRUE)
  
  
  # Step 2: Hydrological Correction of DEM (Depression Breaching)
  # Use whitebox's breach_depressions method to avoid disrupting local flow paths.
  whitebox::wbt_breach_depressions(dem, corrected_dem_path)
  
  # Step 3: Calculate Flow Direction
  # Flow direction grid calculation
  whitebox::wbt_d8_pointer(corrected_dem_path, flow_direction_path)
  
  # Step 4: Calculate Flow Accumulation
  # Flow accumulation raster calculation
  whitebox::wbt_d8_flow_accumulation(corrected_dem_path, flow_accumulation_path)
  
  # Step 5: Identify Ridges
  # Ridges can be identified as areas with low flow accumulation.
  # Load the flow accumulation raster
  flow_accum <- rast(flow_accumulation_path)
  
  # Apply a threshold to define ridges
  # Ridges will have very low flow accumulation; you can adjust the threshold value.
  ridges <- flow_accum < ridge_threshold
 
  
  ridge_lines_sf <- ""
  if (usePolygons) {
    # Convert raster ridges to polygons and simplify geometry
    ridge_polygons <- terra::as.polygons(ridges, dissolve = TRUE)  # Convert raster to polygons
    ridge_polygons <- terra::simplifyGeom(ridge_polygons, tolerance = simplifyTol)  # Simplify geometry
    
    # Convert polygons to contour lines
    ridge_lines <- terra::as.contour(ridges, levels = 1)  # Create contours at level 1
    ridge_lines_sf <- sf::st_as_sf(ridge_lines)  # Convert to sf object for better visualization
    
  } else {
    # Convert logical raster ridges to contour lines
    ridge_lines <- terra::as.contour(ridges, levels = 1)  # Create contours at level 1
    ridge_lines_sf <- sf::st_as_sf(ridge_lines)  # Convert to sf object for better visualization
    
    # Simplify and smooth ridge lines
    ridge_lines_sf <- sf::st_simplify(ridge_lines_sf, dTolerance = simplifyTol)  # Adjust tolerance for smoothness
    ridge_lines_sf <- sf::st_union(ridge_lines_sf)  # Merge connected lines
  }
  
  if(write_ridge_lines){
    # save ridge lines
    vect_ridge_output_path <- file.path(output_dir, "ridge_lines.shp")
    writeVector(ridge_lines, vect_ridge_output_path, overwrite = TRUE)
    
    # Save the ridge raster
    writeRaster(ridges, ridge_output_path, overwrite = TRUE)
  }
  
  
  # Plot the DEM and ridge lines using tmap
  tmap_mode("view")  # Set tmap mode to interactive view
  
  tm <- tm_shape(dem) +
    tm_raster(style = "cont", palette = gray.colors(10), title = "Hillshade") +
    tm_shape(ridge_lines_sf) +  # Add ridge lines to the map
    tm_lines(col = ridge_color, lwd = ridge_lines_width, title.col = "Ridges") +
    tm_layout(
      main.title = "Refined Mountain Ridges",
      legend.outside = TRUE
    )
  
  print(tm)
  
}

get_ridges_from_dem(ridge_threshold = 1800, ridge_color="blue", ridge_lines_width=2,  write_ridge_lines=F, simplifyTol=800,usePolygons = T)