rspatial
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?
Solution
This is the workflow I use in these cases:
- Combine contiguous linestrings using
sf::st_line_merge()to create a single linestring - 'Snap' points to linestring using a combination of
lwgeom::st_endpoint()andsf::st_nearest_points() - Offset snapped points either side of linestring using
geosphere::bearing()andspatialEco::bearing.distance(), then use resulting points to create 'blades' to split linestring - Use blades to split linestring and remove "loose ends", which will always be the first and last row
- Assign values to result using subsetted version of snapped points and calculate cumulative length
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"))
- Convert hex colour to HCL
- dumbell plot with more than one start and end point
- Sequence color blocks inside ggplot2 plot area at specific coordinates
- What does the symbol ::: mean in R
- Reproducing lattice dendrogram graph with ggplot2
- Make a common x-axis for multiple plots in a cowplot
- Is it possible to suppress just one tick mark label?
- R: Color branches of dendrogram while preserving the color legend
- Likert Graphs disappear when moving section labels to top instead of left?
- Adding normal density plot to histogram in ggplot
- Plotting 2 histogram on top and bottom part of a single plot
- How can I create a single paged landscape document in officedown?
- Does ctable in R not work with native pipes?
- How to find the mean number of consecutive days using dates in R?
- How to fix the value of a specific color of Legend?
- ggpplot2 join all points in a graph with a line
- Downloading Shiny Plot as PNG
- How strongly can system software interfere with R package functionality?
- R Markdown doesn't plot leaflet map in Word document
- How to switch from geom_point to geom_bin2d() or geom_tile?
- Extract non null elements from a list
- Error in passing a variable to ggplot in an r custom function
- Plotting stm model expected topic probability using ggplot2
- How to show the number of dropped observations in R?
- R: Replicating samples for a function then cumulatively adding the result
- Optional dependency in R package
- Unnest list column: In row 2, can't recycle input of size 2 to size 4
- Extract gbif occurrences within a geographic administrative area with rgbif
- PDF convert big vector file to small pixel file within R with magick
- Replace variables in a formula with their definitions