如何自动化计算参考水文站与其他站点的沿河距离？

问题描述 #

需要使用nhdplusTools计算参考水文站（all_gages中的sample_origin）与河流上其他多个水文站的沿河距离。现有教程需手动移除河段且仅计算上下游端点距离，需实现流程自动化并获取参考站点到其他站点的距离。现有代码已完成获取流线、站点 snapping 步骤，需修改最后河段分割与距离计算的代码。

现有前期代码（已完成部分）：

library(sf)
library(nhdplusTools)
library(mapview)

# custom function for snapping points to flowlines
st_snap_points <- function(x, y, namevar, max_dist = 1000) {
  
  # this evaluates the length of the data
  if (inherits(x, "sf")) n = nrow(x)
  if (inherits(x, "sfc")) n = length(x)
  
  # this part: 
  # 1. loops through every piece of data (every point)
  # 2. snaps a point to the nearest line geometries
  # 3. calculates the distance from point to line geometries
  # 4. retains only the shortest distances and generates a point at that intersection
  out = do.call(c,
                lapply(seq(n), function(i) {
                  nrst = st_nearest_points(st_geometry(x)[i], y)
                  nrst_len = st_length(nrst)
                  nrst_mn = which.min(nrst_len)
                  if (as.vector(nrst_len[nrst_mn]) > max_dist) return(st_geometry(x)[i])
                  return(st_cast(nrst[nrst_mn], "POINT")[2])
                })
  )
  # this part converts the data to a dataframe and adds a named column of your choice
  out_xy <- st_coordinates(out) %>% as.data.frame()
  out_xy <- out_xy %>% 
    mutate({{namevar}} := x[[namevar]]) %>% 
    st_as_sf(coords=c("X","Y"), crs=st_crs(x), remove=FALSE)
  
  return(out_xy)
}


# GET FLOWLINES ASSOCIATED WITH COMIDS
  # make a list defining the sourcetype and ID
  feat_list <- list(featureSource = "comid", featureID = "13167914")
  
  # get upstream flowlines
  us_flowlines <- navigate_nldi(nldi_feature = feat_list,
                                mode="UM",
                                data_source = "flowlines")
  
  # get downstream mainstem only
  ds_flowlines <- navigate_nldi(nldi_feature = feat_list,
                                mode="DM",
                                data_source = "flowlines")
  
  # make a list of all the comids we've identified:
  all_comids <- c(us_flowlines$UM_flowlines$nhdplus_comid, ds_flowlines$DM_flowlines$nhdplus_comid)
  
  # download all data and create a geopackage with the comid list
  gpkg <- subset_nhdplus(comids = as.integer(all_comids),
                         simplified = TRUE,
                         overwrite = TRUE,
                         output_file = "13167914_nhdplus.gpkg",
                         nhdplus_data = "download",
                         return_data = FALSE)
  
  # pull the flowlines back in
  streams <- read_sf("13167914_nhdplus.gpkg", "NHDFlowline_Network")

# SNAP GAGE POINTS TO FLOWLINE
  # create sf of all gage locations
  all_gages <- structure(list(location = c("sample_origin", "upstream", "downstream", 
                                           "downstream"), 
                              X = c(-83.609581641841, -83.6125472, -83.54465649, -83.52743409), 
                              Y = c(41.6614141167543, 41.65968056, 41.68865998, 41.70310398)), 
                         row.names = c(NA, 4L), class = "data.frame") %>% 
    sf::st_as_sf(coords = c("X", "Y"), crs = 4326)
  
  
  # project all_gages and streams
  all_gages_proj <- sf::st_transform(all_gages, crs = 26910)
  streams_proj <- sf::st_transform(streams, crs = 26910)
  
  # snap points to the flowlines using a 500 m buffer
  gages_snapped <- st_snap_points(all_gages_proj, streams_proj, namevar = "location", max_dist = 500)
  
  # view map
  mapview(gages_snapped, col.regions="orange", layer.name="Snapped Gages") +
    mapview(streams_proj, color="steelblue", layer.name="Flowlines")

需修改的原始河段处理代码：

library(lwgeom)
  
# MEASURE DISTANCE BETWEEN POINTS  
  # Split stream segments between points
  # Create a 1 m buffer around snapped point
  gages_snapped_buff <- st_buffer(gages_snapped, 1)
  
  # Use lwgeom::st_split to split stream segments
  segs <- st_collection_extract(lwgeom::st_split(streams_proj, gages_snapped_buff), "LINESTRING") %>% 
    tibble::rownames_to_column(var = "rowid") %>% 
    mutate(rowid=as.integer(rowid))
  
  # Calculate river distances
  segs_dist <- segs %>% 
    # drop the "loose ends" on either extent (upstream or downstream) of first/last gage
    # here they manually provided rowids from the example
    # filter(!rowid %in% c(232, 100, 66, 62, 63)) %>% 
    mutate(seg_len_m = units::drop_units(units::set_units(st_length(.), "m")),
           seg_len_km = seg_len_m/1000) %>% 
    arrange(desc(hydroseq)) %>% 
    mutate(total_len_km = cumsum(seg_len_km)) %>% 
    # filter to just cols of interest
    select(rowid, comid, gnis_id:reachcode, streamorde, hydroseq, seg_len_km, total_len_km, geom)

解决方案 #

以下是修改后的自动化代码，无需手动指定rowid，可自动计算sample_origin到其他站点的沿河距离：

library(lwgeom)
library(dplyr)
library(sf)

# MEASURE DISTANCE BETWEEN POINTS  
# 1. 分割河段并保留原始属性
gages_snapped_buff <- st_buffer(gages_snapped, 1)
segs <- st_collection_extract(lwgeom::st_split(streams_proj, gages_snapped_buff), "LINESTRING") %>% 
  tibble::rownames_to_column(var = "rowid") %>% 
  mutate(rowid = as.integer(rowid)) %>%
  # 保留原始流线的hydroseq等关键属性（通过空间关联匹配）
  st_join(streams_proj %>% select(comid, hydroseq, streamorde), left = TRUE) %>%
  # 计算每个分割河段的长度
  mutate(seg_len_km = units::drop_units(units::set_units(st_length(.), "km")))

# 2. 定位参考站点sample_origin的位置与对应的hydroseq
ref_point <- gages_snapped %>% filter(location == "sample_origin")
# 找到参考点所在的原始河段hydroseq
ref_hydroseq <- st_join(ref_point, streams_proj %>% select(hydroseq)) %>% pull(hydroseq) %>% unique()

# 3. 按河流流向排序（hydroseq越大越上游）
segs_sorted <- segs %>% 
  arrange(desc(hydroseq)) %>%
  # 标记每个河段相对于参考点的位置：上游/下游/包含参考点
  mutate(relative_position = case_when(
    hydroseq > ref_hydroseq ~ "upstream",
    hydroseq < ref_hydroseq ~ "downstream",
    hydroseq == ref_hydroseq ~ "reference_segment"
  ))

# 4. 计算从参考点出发的累积距离
# 上游方向：从参考点所在河段开始向上游累加
segs_upstream <- segs_sorted %>% 
  filter(relative_position %in% c("upstream", "reference_segment")) %>%
  # 找到参考点所在河段的起始累积距离为0
  mutate(total_from_ref_km = cumsum(seg_len_km) - sum(seg_len_km[relative_position == "reference_segment"]))

# 下游方向：从参考点所在河段开始向下游累加（注意排序是从上游到下游，所以需要反向累加）
segs_downstream <- segs_sorted %>% 
  filter(relative_position %in% c("downstream", "reference_segment")) %>%
  arrange(hydroseq) %>% # 按hydroseq从小到大（下游方向）排序
  mutate(total_from_ref_km = cumsum(seg_len_km) - sum(seg_len_km[relative_position == "reference_segment"]))

# 合并上下游距离数据
segs_dist <- bind_rows(segs_upstream, segs_downstream) %>%
  arrange(desc(hydroseq)) %>%
  select(rowid, comid, hydroseq, streamorde, seg_len_km, total_from_ref_km, relative_position, geom)

# 5. 匹配每个站点到对应的河段，获取到参考点的距离
gage_distances <- gages_snapped %>%
  st_join(segs_dist %>% select(total_from_ref_km, relative_position), left = TRUE) %>%
  # 去重，每个站点保留唯一的距离值
  distinct(location, .keep_all = TRUE) %>%
  select(location, total_from_ref_km, relative_position, X, Y, geom)

# 查看结果
print(gage_distances)

# 可视化验证
mapview(segs_dist, zcol = "total_from_ref_km", layer.name = "Distance from Reference (km)") +
  mapview(gage_distances, col.regions = "red", layer.name = "Gages & Distances")

代码说明 #

• 自动关联属性：通过空间连接保留原始流线的hydroseq（河流排序字段，值越大越上游），避免手动识别河段。
• 参考点定位：自动找到sample_origin所在河段的hydroseq，以此为基准划分上下游。
• 累积距离计算：分别对上下游河段计算从参考点出发的累积距离，无需手动移除“松散端点”——因为我们只保留了参考点和其他站点之间的河段（通过空间连接过滤）。
• 站点距离匹配：将每个snapped站点与对应的河段连接，直接获取到参考点的沿河距离，输出清晰的结果表格。

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