Question
How do I measure distances between multiple coordinates along a stream/flowline?
I am trying to measure the distance between a reference stream gage/location/coordinates and multiple other gages on a stream/river that were discovered using nhdplusTools. I am following a method described here. However, the blogger manually removes stream segments at the end (I need to automate this process) and ultimately only measures the distance between the most upstream and downstream points.
Here is the the code I have used so far:
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")
The blogger then uses the following code to split the stream into segments. But again, they manually remove stream segments (I need to automate this process) and only measure the distance between upstream and downstream points (I want to measure the distance between a reference gage/location - called 'sample_origin' in the all_gages object - and the other stream gages).
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)
How do I alter this final code chunk to find the distances between the reference gage/location - called 'sample_origin' in the all_gages object - and the other stream gages?
