Как измерить расстояния между несколькими координатами вдоль потока/линии потока?

Вот рабочий процесс, который я использую в этих случаях:

1. Объедините смежные строки, используя sf::st_line_merge(), чтобы создать одну строку.
2. «Привязка» указывает на строку, используя комбинацию lwgeom::st_endpoint() и sf::st_nearest_points().
3. Сместите привязанные точки по обе стороны от линии, используя geosphere::bearing() и spatialEco::bearing.distance(), затем используйте полученные точки для создания «лезвий» для разделения линии.
4. Используйте лезвия, чтобы разделить леску и удалить «свободные концы», которые всегда будут первым и последним рядом.
5. Присвойте значения результату, используя подмножество версий точек привязки, и вычислите совокупную длину.

library(nhdplusTools)
library(dplyr)
library(tidyr)
library(sf)
library(lwgeom)
library(geosphere)
library(spatialEco)
library(ggplot2)

# Load flowlines (previously created from your code example)
streams <- read_sf("13167914_nhdplus.gpkg", "NHDFlowline_Network")

# Create single linestring of streams object
streams1 <- st_transform(streams, crs = 4326) %>%
  st_combine() %>%
  st_line_merge() %>%
  st_as_sf()

# Create sf of 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") %>% 
  st_as_sf(coords = c("X", "Y"), crs = 4326)

# Create sf of nearest points of all_gages on streams1 (replaces st_snap_points())
nearest <- all_gages %>% 
  mutate(rowid1 = 1:n(),
         nearest = st_endpoint(st_nearest_points(., streams1)),
         azimuth1 = bearing(st_coordinates(geometry),
                            st_coordinates(nearest)),
         azimuth2 = (azimuth1 + 180) %% 360) %>%
  st_set_geometry("nearest")

# Filter nearest to remove sample_origin (for joining attritubes to final result)
nearest1 <- filter(nearest, location != "sample_origin") %>%
  mutate(rowid1 = 1:n()) %>%
  select(rowid1, location)

# Offset modified all_gages points by ~100mm using azimuth
offset_p1 <- data.frame(bearing.distance(st_coordinates(nearest)[,1], 
                                         st_coordinates(nearest)[,2], 
                                         0.1/111320, 
                                         nearest$azimuth1),
                        rowid1 = 1:nrow(nearest),
                        azimuth1 = nearest$azimuth1) %>%
  st_as_sf(coords = c("X", "Y"), crs = 4326) %>%
  rename(geometry1 = geometry)

offset_p2 <- data.frame(bearing.distance(st_coordinates(nearest)[,1], 
                                         st_coordinates(nearest)[,2], 
                                         0.1/111320, 
                                         nearest$azimuth2),
                        rowid1 = 1:nrow(nearest),
                        azimuth2 = nearest$azimuth2) %>%
  st_as_sf(coords = c("X", "Y"), crs = 4326) %>%
  rename(geometry2 = geometry)

# Create linestring blades to split stream1
blades <- left_join(offset_p1, data.frame(offset_p2), by = "rowid1") %>%
  select(rowid1, geometry1, geometry2) %>%
  data.frame() %>%
  pivot_longer(starts_with("geo"),
               values_to = "geometry") %>%
  data.frame() %>%
  st_as_sf() %>%
  summarise(geometry = st_combine(geometry), .by = c(rowid1)) %>%
  st_cast("LINESTRING") %>%
  left_join(select(st_drop_geometry(nearest), rowid1, location), by = "rowid1")

# Split streams1 using blades, drop "loose ends",
# assign location values using nearest1,
# calculate individual and cumulative length, project
segs <- st_collection_extract(st_split(streams1, blades), "LINESTRING") %>%
  slice(2:(n()-1)) %>%
  mutate(rowid = 1:n(),
         seg_len_km = as.numeric(st_length(.) / 1000),
         rowid1 = st_nearest_feature(., nearest1)) %>%
  left_join(data.frame(nearest1), by = "rowid1") %>%
  group_by(location) %>%
  arrange(if_else(location == "downstream", rowid, desc(rowid)), .by_group = TRUE) %>%
  mutate(cumu_len_km = cumsum(seg_len_km)) %>%
  ungroup() %>%
  st_transform(26910) %>%
  select(rowid, location, seg_len_km, cumu_len_km)

segs
# Simple feature collection with 3 features and 4 fields
# Geometry type: LINESTRING
# Dimension:     XY
# Bounding box:  xmin: 3790633 ymin: 5432731 xmax: 3795232 ymax: 5441401
# Projected CRS: NAD83 / UTM zone 10N
# # A tibble: 3 × 5
#   rowid location   seg_len_km cumu_len_km                                 x
#   <int> <chr>           <dbl>       <dbl>                  <LINESTRING [m]>
# 1     2 downstream      7.18        7.18  (3790785 5433050, 3790812 543309…
# 2     3 downstream      2.37        9.55  (3794443 5439086, 3794501 543914…
# 3     1 upstream        0.314       0.314 (3790633 5432731, 3790672 543277…


# For illustrative purposes only
nearest2 <- left_join(nearest1, data.frame(segs), by = join_by(rowid1 == rowid))

# Plot
ggplot() +
  geom_sf(data = st_transform(streams1, 26910), colour = "grey", linewidth = .2) +
  geom_sf(data = segs, 
          aes(colour = location),
          linewidth = 1) +
  geom_sf(data = st_transform(nearest, 26910), aes(shape = location), size = 1.3) +
  geom_sf_text(data = nearest2,
               aes(label = paste0(round(cumu_len_km, 2), "km")),
               size = 3,
               nudge_x = -1500) +
  theme(panel.background = element_blank(),
        legend.position = "inside", legend.position.inside =  c(0.25, 0.8),
        legend.title = element_blank(),
        legend.spacing.y = unit(-3, "mm"))