tidyversegisr-sf
sf methods to delineate a centerline from a channel
I am hoping to use Voronoi method to delineate a river channel centerline. I came across multiple sources but failed to replicate the work due to retired pacakges (rgeos).
i.e. https://gist.github.com/FloodHydrology/829d3846f0722610ff69d896ad65af3e Identify a linear feature on a raster map and return a linear shape object using R
library(sf)
library(tibble)
library(sf)
channel <- st_as_sf(tribble(
~x, ~y,
656873.7, 5553711,
656993.1, 5553738,
656993.6, 5553738,
656994.1, 5553738,
656994.6, 5553738,
656995.1, 5553738,
656995.7, 5553738,
656996.2, 5553738,
656996.7, 5553738,
656997.2, 5553738,
656997.7, 5553738,
656998.2, 5553738,
656998.7, 5553738,
656999.2, 5553738,
656999.7, 5553737,
657000.2, 5553737,
657000.6, 5553737,
657001.0, 5553737,
657001.5, 5553736,
657001.9, 5553736,
657002.3, 5553736,
657002.6, 5553735,
657003.0, 5553735,
657003.3, 5553734,
657003.6, 5553734,
657126.8, 5553550,
657127.0, 5553550,
657127.3, 5553549,
657127.5, 5553549,
657127.7, 5553548,
657127.9, 5553548,
657128.1, 5553548,
657128.2, 5553547,
657128.3, 5553546,
657128.4, 5553546,
657128.4, 5553545,
657128.4, 5553545,
657129.8, 5553466,
657129.8, 5553465,
657129.7, 5553465,
657129.7, 5553464,
657129.6, 5553464,
657129.5, 5553463,
657129.4, 5553463,
657129.2, 5553462,
657129.0, 5553462,
657128.8, 5553461,
657128.6, 5553461,
657128.3, 5553460,
657128.1, 5553460,
657127.8, 5553459,
657127.4, 5553459,
657127.1, 5553459,
657126.7, 5553458,
657126.4, 5553458,
657126.0, 5553458,
656995.2, 5553355,
656971.7, 5553303,
656975.2, 5553281,
657061.9, 5553229,
657233.8, 5553179,
657234.3, 5553179,
657234.8, 5553179,
657235.3, 5553179,
657235.9, 5553178,
657236.3, 5553178,
657312.0, 5553130,
657312.5, 5553130,
657312.9, 5553130,
657313.3, 5553129,
657313.7, 5553129,
657314.1, 5553128,
657314.5, 5553128,
657314.8, 5553128,
657315.1, 5553127,
657315.4, 5553127,
657315.6, 5553126,
657315.9, 5553126,
657316.1, 5553125,
657316.2, 5553125,
657406.2, 5552837,
657467.5, 5552813,
657534.9, 5552855,
657535.4, 5552855,
657535.9, 5552855,
657536.3, 5552855,
657536.8, 5552855,
657537.3, 5552856,
657537.8, 5552856,
657538.3, 5552856,
657538.8, 5552856,
657539.4, 5552856,
657539.9, 5552856,
657540.4, 5552856,
657540.9, 5552856,
657541.4, 5552856,
657542.0, 5552856,
657542.5, 5552856,
657543.0, 5552856),
coords = c("x", "y")) %>%
st_combine() %>%
st_cast("LINESTRING")
channel <- st_buffer(channel, 100)
st_voronoi(channel) %>%
st_collection_extract()
I got as far as :
But not sure how to use the edges of the Voronoi object to obtain centerpoints, then centerlines.
Solution
Here's a simplified approach just to illustrate general workflow, though for the actual task, the centerline package suggested in comments might save you a ton of time and you'd probably get better results, too.
Suitability of Voronoi tesselation for centerline detection depends highly on vertex density of the shape and vertex distribution, so here I'm just adding more points, in the first example with st_line_sample() (uniform distribution, nice for visualization, might miss some critical points at bends), and in the 2nd example with st_segmentize() (probably introduces some artefacts at centerline ends).
As we are after Vornoi edges and do not really care about polygons, we can go wiht st_voronoi(... , bOnlyEdges = TRUE)
From the resulting edges, we are only interested in those that are within the channel polygon, so we can use a spatial filter to remove most (but likely not all) of the clutter. From there, we can think of it as a graph problem -- centerline is the longest geodesic, the graph diameter.
tidygraph / igarph and sf are neatly brought together in sfnetworks package, which lets us conveniently convert our sf object of remaining Voronoi edges to a graph. And after graph processing, we can convert it back to sf.
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(dplyr)
library(sfnetworks)
library(tidygraph)
library(ggplot2)
# uniform point distribution along polygon line
# for simplification & visualization
edge_points <-
channel %>%
st_cast("LINESTRING") %>%
st_line_sample(density = 0.02) %>%
st_union()
edge_points
#> Geometry set for 1 feature
#> Geometry type: MULTIPOINT
#> Dimension: XY
#> Bounding box: xmin: 656774.6 ymin: 5552713 xmax: 657642.6 ymax: 5553838
#> CRS: NA
#> MULTIPOINT ((656930.3 5553621), (656880.9 55536...
# we can set bOnlyEdges to get edges instead of polygons
vor_edges <-
st_voronoi(edge_points, bOnlyEdges = TRUE) %>%
# from single MULTILINESTRING to a set of LINESTRINGs
st_cast("LINESTRING")
vor_edges
#> Geometry set for 183 features
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 655649.7 ymin: 5551588 xmax: 658767.4 ymax: 5554963
#> CRS: NA
#> First 5 geometries:
#> LINESTRING (655649.7 5554589, 656873.7 5553711)
#> LINESTRING (656873.7 5553711, 656873.7 5553711)
#> LINESTRING (656873.7 5553711, 656873.7 5553711)
#> LINESTRING (656873.7 5553711, 655649.7 5553857)
#> LINESTRING (656873.7 5553711, 656873.6 5553711)
# keep only edges within channel
edge_within <-
vor_edges %>%
# for convenience of st_filter we need to convert sfc to sf
st_sf() %>%
st_filter(channel, .predicate = st_within)
edge_within
#> Simple feature collection with 68 features and 0 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 656873.6 ymin: 5552816 xmax: 657543.1 ymax: 5553731
#> CRS: NA
#> First 10 features:
#> geometry
#> 1 LINESTRING (656873.7 555371...
#> 2 LINESTRING (656873.7 555371...
#> 3 LINESTRING (656873.7 555371...
#> 4 LINESTRING (656996.3 555330...
#> 5 LINESTRING (656997.3 555329...
#> 6 LINESTRING (656873.8 555371...
#> 7 LINESTRING (656873.8 555371...
#> 8 LINESTRING (656996.6 555330...
#> 9 LINESTRING (656900.8 555371...
#> 10 LINESTRING (656900.8 555371...
# to remove edges that are not part of the centre line,
# yet are too short to intersect with channel polygon outline,
# we can convert the remaining list of edges to a graph,
# keep only nodes (and edges) that are part of a graph diameter (longest geodesic),
# remove pseudo nodes (nodes with 2 incident edges) with to_spatial_smooth morpher and
# convert resulting graph to a single LINESTRING
channel_cline <-
edge_within %>%
as_sfnetwork(directed = FALSE) %>%
activate(nodes) %>%
# get_diameter returns a list node tbl indices, so we can use it with slice
slice(igraph::get_diameter(.) %>% as.vector()) %>%
convert(to_spatial_smooth) %>%
st_as_sf(active = "edges")
channel_cline
#> Simple feature collection with 1 feature and 3 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 656873.6 ymin: 5552831 xmax: 657543.1 ymax: 5553725
#> CRS: NA
#> # A tibble: 1 × 4
#> from to .tidygraph_edge_index geometry
#> <int> <int> <list> <LINESTRING>
#> 1 1 2 <int [63]> (656873.6 5553711, 656873.7 5553711, 656873…
# visualize previous steps ---------------------------------------------------
# cropped variant for plot
vor_crop <- st_crop(vor_edges, st_buffer(channel, 100))
p1 <- ggplot() +
geom_sf(data = channel, fill = "lightblue") +
geom_sf(data = edge_points) +
# Calssify edges (within / not within channel polygon)
geom_sf(data = vor_crop, aes(color = st_within(vor_crop, channel, sparse = FALSE)), show.legend = FALSE) +
ggtitle("Vor. edges from \nunif. sampled points") +
theme_void()
p2 <- ggplot() +
geom_sf(data = channel, fill = "lightblue") +
geom_sf(data = edge_within) +
ggtitle("Vor. edges within \nchannel poly") +
theme_void() +
theme(legend.position = "bottom", legend.title.position = "top")
p3 <- ggplot() +
geom_sf(data = channel, fill = "lightblue") +
geom_sf(data = channel_cline) +
ggtitle("Longest geodesic\n") +
theme_void() +
theme(legend.position = "bottom", legend.title.position = "top")
patchwork::wrap_plots(p1, p2, p3)
As st_voronoi works with a polygon, we can just add some vertex density to that channel edge with st_segmentize(), so in a single pipeline it might look something like this:
cline_sfc <-
channel %>%
# we can add more points to the polygon outline and pass it directly to st_voronoi
st_segmentize(dfMaxLength = 10) %>%
st_voronoi(bOnlyEdges = TRUE) %>%
st_cast("LINESTRING") %>%
st_sf() %>%
st_filter(channel, .predicate = st_within) %>%
as_sfnetwork(directed = FALSE) %>%
activate(nodes) %>%
slice(igraph::get_diameter(.) %>% as.vector()) %>%
# there's a good chance that there are few artefacts at the beginning and end,
# naive but perhaps good enough approach would be just to remove few nodes
# from both ends of the diameter list:
# slice(igraph::get_diameter(.) %>% as.vector() %>% head(-2) %>% tail(-2)) %>%
convert(to_spatial_smooth) %>%
st_as_sf(active = "edges") %>%
st_geometry()
cline_sfc
#> Geometry set for 1 feature
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 656873.7 ymin: 5552826 xmax: 657553.8 ymax: 5553729
#> CRS: NA
#> LINESTRING (656873.7 5553711, 656873.7 5553711,...
ggplot() +
geom_sf(data = channel, fill = "lightblue") +
geom_sf(data = cline_sfc) +
ggtitle("cline_sfc") +
theme_void()
Input shape:
channel <- st_as_sf(tribble(
~x, ~y,
656873.7, 5553711,
656993.1, 5553738,
656993.6, 5553738,
656994.1, 5553738,
656994.6, 5553738,
656995.1, 5553738,
656995.7, 5553738,
656996.2, 5553738,
656996.7, 5553738,
656997.2, 5553738,
656997.7, 5553738,
656998.2, 5553738,
656998.7, 5553738,
656999.2, 5553738,
656999.7, 5553737,
657000.2, 5553737,
657000.6, 5553737,
657001.0, 5553737,
657001.5, 5553736,
657001.9, 5553736,
657002.3, 5553736,
657002.6, 5553735,
657003.0, 5553735,
657003.3, 5553734,
657003.6, 5553734,
657126.8, 5553550,
657127.0, 5553550,
657127.3, 5553549,
657127.5, 5553549,
657127.7, 5553548,
657127.9, 5553548,
657128.1, 5553548,
657128.2, 5553547,
657128.3, 5553546,
657128.4, 5553546,
657128.4, 5553545,
657128.4, 5553545,
657129.8, 5553466,
657129.8, 5553465,
657129.7, 5553465,
657129.7, 5553464,
657129.6, 5553464,
657129.5, 5553463,
657129.4, 5553463,
657129.2, 5553462,
657129.0, 5553462,
657128.8, 5553461,
657128.6, 5553461,
657128.3, 5553460,
657128.1, 5553460,
657127.8, 5553459,
657127.4, 5553459,
657127.1, 5553459,
657126.7, 5553458,
657126.4, 5553458,
657126.0, 5553458,
656995.2, 5553355,
656971.7, 5553303,
656975.2, 5553281,
657061.9, 5553229,
657233.8, 5553179,
657234.3, 5553179,
657234.8, 5553179,
657235.3, 5553179,
657235.9, 5553178,
657236.3, 5553178,
657312.0, 5553130,
657312.5, 5553130,
657312.9, 5553130,
657313.3, 5553129,
657313.7, 5553129,
657314.1, 5553128,
657314.5, 5553128,
657314.8, 5553128,
657315.1, 5553127,
657315.4, 5553127,
657315.6, 5553126,
657315.9, 5553126,
657316.1, 5553125,
657316.2, 5553125,
657406.2, 5552837,
657467.5, 5552813,
657534.9, 5552855,
657535.4, 5552855,
657535.9, 5552855,
657536.3, 5552855,
657536.8, 5552855,
657537.3, 5552856,
657537.8, 5552856,
657538.3, 5552856,
657538.8, 5552856,
657539.4, 5552856,
657539.9, 5552856,
657540.4, 5552856,
657540.9, 5552856,
657541.4, 5552856,
657542.0, 5552856,
657542.5, 5552856,
657543.0, 5552856),
coords = c("x", "y")) %>%
st_combine() %>%
st_cast("LINESTRING")
channel <- st_buffer(channel, 100)
- Avoiding NAs in the output of dplyr
- Create Unevaluated Calls with Tidyverse Defusing Operators
- what happened when `separate_wider_regex` in multiple match
- How do you code the creation of fixture lists for a sports league?
- Create new variable by multiplying corresponding variables and summing
- How to print the number of observations dropped by tidyverse's functions like filter or drop_na?
- How to I remove networks with less than 10 nodes using the egor package?
- tidyverse: Cross tables of one variable with all other variables in data.frame
- Variable grouping in gtsummary with dichotomous variable
- readr has no reasons to say there are less columns than expected
- Error out if filtering for a non-existent level
- R: Assign column value based on range of values found in columns in another dataframe
- Problems with mutate(), across (), and Winsorize ()
- Find total length of time in multiple values of a dataset that overlaps with a given period of time
- Mutate multiple variable to create multiple new variables
- Using vector with strings in mutate with paste
- Can you use dplyr across() to iterate across pairs of columns?
- Filter data based on an external column using dplyr
- How to calculate rolling standard deviations in dplyr
- round_any equivalent for dplyr?
- Multiple stacked bar chart with ggplot
- Use R function separate() on all dataframes in a list
- Re-write tidyverse function for clarity and brevity
- 'INDEX' is of length zero error in gtsummary::tbl_svysummary
- Apply a custom function to any number of columns using dplyr
- Pass a named list or character value to dplyr’s mutate to create multiple columns in a dataframe
- Check if element of one dataframe is in another dataframe, within group
- How to regroup lines and values when there is no gap between dates?
- Split a list into separate data frame in R
- retrieve .RData file of nameless object created with `magrittr::%>%`