Group ID for all edges between points using sfnetwork in r
I have a rooted tree with spatially explicit edges (ln_sfnetwork) with additional edges created by adding a point layer (pt).
I would like to give all of the edges between each point on the network the same ID so I can calculate the total length of the network between points. I have a manual solution but this has to be done with large datasets > 20,000 points.
install.packages("remotes"); library("remotes")
install_github("luukvdmeer/sfnetworks")
library(sf)
library(sfnetworks)
library(dplyr)
ln <- structure(list(River_ID = c(159, 160, 161, 186, 196), geometry = structure(list(
structure(c(
289924.625, 289924.5313, 289922.9688, 289920.0625,
289915.7499, 289912.7188, 289907.4375, 289905.3438, 289901.1251,
289889, 289888.5, 289887.5938, 289886.5, 289886.4063, 289885.3124,
289884.0938, 289884.0001, 289882.8125, 289881.625, 289878.6875,
289877.9688, 289876.25, 289874.5625, 289874.25, 289872.7188,
289871.2813, 289871.1875, 289870.0313, 289869, 289868.5939,
289867.8436, 289865.8438, 289864.0625, 289862.5939, 289862.375,
289861.5, 289860.7812, 289860.5625, 289859.5313, 289858.375,
289857.7813, 289855.4063, 289854.25, 289850.8749, 289846.4376,
289841.9064, 289836.0625, 289828.1562, 289822.8438, 289816.625,
289812.4376, 289807.9064, 289798.75, 289793.125, 289786.2188,
289781.375, 289777.3124, 289770.0313, 289765.4375, 289762.2188,
289759.25, 289755.5938, 289753.0625, 289747.9687, 289743.7499,
289741.5938, 289739.5, 289736.1874, 289732.75, 289727, 289723.7499,
289719.625, 289715.5626, 289713.7499, 202817.531300001, 202817.2031,
202815.1094, 202812.468699999, 202809.3906, 202806.7656,
202799.7969, 202797.906300001, 202794.093800001, 202783.515699999,
202783.125, 202782.4844, 202781.906300001, 202781.8125, 202781.3594,
202781.093800001, 202780.9999, 202780.5469, 202780, 202777.625,
202777.0469, 202775.718800001, 202774.1875, 202773.906300001,
202772.1875, 202770.4531, 202770.25, 202768.5156, 202766.6719,
202766, 202764.0469, 202759.6719, 202755.8749, 202752.781300001,
202752.1875, 202749.953199999, 202748.297, 202747.906300001,
202746.0625, 202744.2344, 202743.5625, 202740.4375, 202738.8125,
202734.5, 202727.9844, 202723.5625, 202719.1875, 202714.9845,
202713.031300001, 202710.6875, 202710.0469, 202711.406300001,
202714.5626, 202716.9845, 202718.718900001, 202719.5469,
202718.734300001, 202716.4531, 202715.125, 202713.7344, 202712.093800001,
202709.8749, 202708.875, 202709.2655, 202710.7031, 202712.375,
202712.375, 202712.2344, 202711.0469, 202707.906300001, 202705.406300001,
202703.0469, 202701.468800001, 202700.7656
), .Dim = c(
74L,
2L
), class = c("XY", "LINESTRING", "sfg")), structure(c(
289954.375,
289953.5, 289950.6562, 289949.7499, 289949, 289948.125, 289946.0625,
289945.9688, 289944.5313, 289943.4063, 289941.3438, 289939.4375,
289937.4375, 289935.1875, 289932.75, 289930.625, 289928.8125,
289928.25, 289926.7188, 289925.5313, 289925.7813, 289925.625,
289925.4063, 289925.1251, 289924.625, 202872.75, 202872.031400001,
202868.7031, 202867.343699999, 202864.906199999, 202861.515699999,
202858.297, 202854.406300001, 202851.9375, 202849.468800001,
202847.703, 202846.75, 202845.4531, 202843.6719, 202843.0625,
202841.593900001, 202839.7344, 202839.2344, 202838, 202835.9375,
202832.875, 202825.7344, 202822.9531, 202819.4531, 202817.531300001
), .Dim = c(25L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(
290042.6563, 290042.3437, 290041.5313, 290038.4376,
290037.625, 290036.5313, 290035.5313, 290034.8438, 290034.5313,
290033.7188, 290032.9375, 290032.125, 290030.3437, 290030.0313,
290028.625, 290027.5626, 290027.3438, 290026.7188, 290024.5313,
290023.625, 290020.625, 290018.0001, 290014.9375, 290012.0938,
290008.5625, 290004.375, 290000.0001, 289999.875, 289997.625,
289993.7188, 289990.5, 289987.1562, 289985.4063, 289980.375,
289973.3124, 289966.375, 289961.8438, 289959, 289954.375,
202884.0625, 202884.25, 202884.843800001, 202888.4531, 202889.75,
202891.0469, 202892.0469, 202892.656300001, 202892.843800001,
202893.2501, 202893.5469, 202893.656300001, 202893.4531,
202893.4531, 202893.343699999, 202893.093800001, 202893.0469,
202892.843800001, 202891.953199999, 202891.5469, 202889.843800001,
202888.218800001, 202885.1094, 202880.9219, 202877.5625,
202873.968800001, 202872.5469, 202872.5156, 202872.625, 202874.5469,
202876.734300001, 202878.1719, 202877.953199999, 202876.3125,
202873.468800001, 202872.031400001, 202872.906199999, 202873.0781,
202872.75
), .Dim = c(39L, 2L), class = c(
"XY", "LINESTRING",
"sfg"
)), structure(c(
290054.125, 290053.4375, 290052.5313,
290051.625, 290050.0313, 290048.125, 290044.125, 290040.4376,
290039.4375, 290036.9688, 290031.4375, 290027.5312, 290024.8125,
290021.7499, 290020.9688, 290018.3437, 290015, 290010.25,
290006.0313, 290002.4376, 290000.0001, 289999.2187, 289996.6875,
289995.3438, 289994.125, 289991.1875, 289989.2187, 289987.9688,
289986.125, 289980.5313, 289975.0314, 289970.9063, 289968.5625,
289961.0312, 289948.0001, 289939.625, 289933.1563, 289928.3125,
289926.5313, 289924.625, 202835.953199999, 202835.656300001,
202835.4531, 202835.343699999, 202835.5469, 202835.7656,
202836.25, 202836.4531, 202836.5469, 202836.5469, 202835.953199999,
202836.031400001, 202836.625, 202837.7969, 202838.4844, 202839.343699999,
202836.25, 202832.7656, 202832.3125, 202833.4844, 202834.4844,
202834.8125, 202834.2344, 202832.625, 202830.625, 202828.593800001,
202828.968800001, 202831.0625, 202833.2655, 202835.5781,
202838, 202838.906199999, 202839.125, 202836.4531, 202830.781300001,
202827.093800001, 202823.625, 202818.5, 202817.5625, 202817.531300001
), .Dim = c(40L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(
290042.625, 290042.0313, 290041.2187, 290040.3125,
290038.4063, 290037.7188, 290035.8125, 290033.7188, 290030.9063,
290028.2187, 290021.5313, 290021.2187, 290014.2188, 290013.4063,
290012.3125, 290010.0625, 290007.9375, 290005.9688, 290004.125,
290000.0001, 289999.4063, 289998.3125, 289997.5312, 289996.8438,
289993.625, 289993.0314, 289989.7188, 289989.3438, 289987.625,
289987.2187, 289984.0313, 289978.125, 289977.9375, 289974.3437,
289972.7188, 289970.9375, 289967.9375, 289965.2187, 289965.1563,
289962.3437, 289960.5313, 289959.1251, 289959.0314, 289959.3438,
289959.4375, 289959.4375, 289959.3438, 289959.2187, 289958.9375,
289958.5313, 289956.125, 289954.375, 202953.781300001, 202952.4844,
202951.281300001, 202950.0781, 202948.1875, 202947.5781,
202945.8749, 202944.281300001, 202941.781300001, 202940.1875,
202936.375, 202936.1875, 202931.968800001, 202931.4844, 202930.875,
202929.093800001, 202927.1094, 202925.031300001, 202922.734300001,
202917.2031, 202916.4375, 202915.2031, 202914.5469, 202914.4531,
202911.4531, 202910.843800001, 202908.0469, 202907.75, 202906.75,
202906.5469, 202904.843800001, 202901.843800001, 202901.75,
202900.0469, 202899.156400001, 202898.0469, 202894.656300001,
202892.0469, 202891.9844, 202889.343699999, 202887.656300001,
202885.75, 202884.5469, 202883.343699999, 202882.5469, 202881.343699999,
202880.0469, 202879.343699999, 202877.656300001, 202876.25,
202874.25, 202872.75
), .Dim = c(52L, 2L), class = c(
"XY",
"LINESTRING", "sfg"
))
), n_empty = 0L, crs = structure(list(
input = "OSGB 1936 / British National Grid", wkt = "PROJCRS[\"OSGB 1936 / British National Grid\",\n BASEGEOGCRS[\"OSGB 1936\",\n DATUM[\"OSGB 1936\",\n ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4277]],\n CONVERSION[\"British National Grid\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",49,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-2,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",400000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",-100000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"unknown\"],\n AREA[\"UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E\"],\n BBOX[49.75,-9.2,61.14,2.88]],\n ID[\"EPSG\",27700]]"
), class = "crs"), class = c(
"sfc_LINESTRING",
"sfc"
), precision = 0, bbox = structure(c(
xmin = 289713.7499,
ymin = 202700.7656, xmax = 290054.125, ymax = 202953.781300001
), class = "bbox"))), row.names = c(NA, -5L), class = c(
"sf",
"data.frame"
), sf_column = "geometry", agr = structure(c(River_ID = NA_integer_), .Label = c(
"constant",
"aggregate", "identity"
), class = "factor"))
pt <- structure(list(lat = c(
202805.8942, 202836.136, 202872.9487,
202905.3284
), lng = c(
289912.0584, 290014.8446, 290001.2364,
289984.9382
), id = 1:4, geometry = structure(list(structure(c(
289912.058400425,
202805.894199679
), class = c("XY", "POINT", "sfg")), structure(c(
290014.844597566,
202836.136003318
), class = c("XY", "POINT", "sfg")), structure(c(
290001.236395958,
202872.948712436
), class = c("XY", "POINT", "sfg")), structure(c(
289984.938209474,
202905.32838227
), class = c("XY", "POINT", "sfg"))), n_empty = 0L, crs = structure(list(
input = "OSGB 1936 / British National Grid", wkt = "PROJCRS[\"OSGB 1936 / British National Grid\",\n BASEGEOGCRS[\"OSGB 1936\",\n DATUM[\"OSGB 1936\",\n ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4277]],\n CONVERSION[\"British National Grid\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",49,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-2,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",400000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",-100000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"unknown\"],\n AREA[\"UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E\"],\n BBOX[49.75,-9.2,61.14,2.88]],\n ID[\"EPSG\",27700]]"
), class = "crs"), class = c(
"sfc_POINT",
"sfc"
), precision = 0, bbox = structure(c(
xmin = 289912.058400425,
ymin = 202805.894199679, xmax = 290014.844597566, ymax = 202905.32838227
), class = "bbox"))), row.names = c(NA, 4L), class = c(
"sf",
"data.frame"
), sf_column = "geometry", agr = structure(c(
lat = NA_integer_,
lng = NA_integer_,
id = NA_integer_
), .Label = c(
"constant",
"aggregate", "identity"
), class = "factor"))
Make it an sfnetwork
ln_sfnetwork <- as_sfnetwork(ln)
Add the point geometry
ln_sfnetwork <- st_network_blend(ln_sfnetwork, st_geometry(pt))
Use the point geom to create new edges by essentially splitting the sfnetwork
ln_sf <- ln_sfnetwork %>%
activate("edges") %>%
mutate(new_river_id = as.character(1:n())) %>%
st_as_sf()
The manual bit to group all of the edges between the nodes from pt. This groups all of the edges occuring between the points.
ln_sf[["new_river_id"]][2] <- "1"
ln_sf[["new_river_id"]][c(1, 3, 7, 5, 9)] <- "2"
ln_sf[["new_river_id"]][6] <- "3"
ln_sf[["new_river_id"]][4] <- "4"
ln_sf[["new_river_id"]][8] <- "5"
Desired output
Then group_by to get the combined length
ln_sf_merged <- ln_sf %>%
group_by(new_river_id) %>%
summarise(fraglen = sum(seglen)) %>%
as.data.frame() %>%
select(-geometry)
Finally
out_sf <- ln_sf %>%
left_join(ln_sf_merged, by = "new_river_id")
Solution
Interesting problem and nice to see sfnetworks being used on other type of networks than road networks! I gave this some thought. Beware, the answer is long and may be complicated sometimes.
It is clear what you want the output of your example to be, but the example does not cover many different possibilities of how the network structure can look like after blending in the given points. Therefore it is hard to generalize the rules that should be used for merging. But I made some assumptions regarding that:
As far as I understand you want to group edges based on what the first point is that you pass when travelling downstream on the river network. Hence, if the paths downstreams (towards the root of the tree) from edge 1 and edge 2 both pass through point 1 before passing any of the other points, they should be in the same group. All edges from which downstream travelling reaches the root without passing through any of the points should be together in a single group as well.
So the ultimate goal is to assign a group index to each of the edges. In tidygraph (the library on which sfnetworks is build), there are several grouping functions that can be applied to either the nodes or edges of the network. Such grouping functions are meant to be used inside tidyverse verbs like mutate() and filter(), where the network that is being worked on is known and not needed as an argument to the grouping function. Hence, you can run network %>% mutate(group = group_components()), without the need to explicitly forward the network as argument to group_components().
All nice and well, but as far as I know tidygraph does not have a grouping function implemented that addresses your problem. That means, we should create our own! I gave it a try. It may not be exactly what you need, but at least it should give you a good starting point to build further on.
The idea is that we first calculate the distances from the nodes that correspond with the points you added to the network, to all nodes in the network. Then, we select for each point a set of nodes such that the travel time from that point to these nodes is lower than from any other point. The edges that connect the nodes in such a set form the group of edges belonging to that point.
I created another example network which covers more possible cases. The first part of the code below is mainly constructing that network. Then we get to writing the custom group function, and finally applying it to my example network, and also to yours.
library(sfnetworks)
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.0, PROJ 7.2.0
library(dplyr)
#>
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#>
#> filter, lag
#> The following objects are masked from 'package:base':
#>
#> intersect, setdiff, setequal, union
library(tidygraph)
#>
#> Attaching package: 'tidygraph'
#> The following object is masked from 'package:stats':
#>
#> filter
# Create an example network.
n01 = st_sfc(st_point(c(0, 0)))
n02 = st_sfc(st_point(c(1, 2)))
n03 = st_sfc(st_point(c(1, 3)))
n04 = st_sfc(st_point(c(1, 4)))
n05 = st_sfc(st_point(c(2, 1)))
n06 = st_sfc(st_point(c(2, 3)))
n07 = st_sfc(st_point(c(2, 4)))
n08 = st_sfc(st_point(c(3, 2)))
n09 = st_sfc(st_point(c(3, 3)))
n10 = st_sfc(st_point(c(3, 4)))
n11 = st_sfc(st_point(c(4, 2)))
n12 = st_sfc(st_point(c(4, 4)))
from = c(1, 2, 2, 3, 3, 5, 5, 8, 8, 9, 9)
to = c(5, 3, 6, 4, 7, 2, 8, 9, 11, 10, 12)
nodes = st_as_sf(c(n01, n02, n03, n04, n05, n06, n07, n08, n09, n10, n11, n12))
edges = data.frame(from = from, to = to)
G_1 = sfnetwork(nodes, edges)
#> Checking if spatial network structure is valid...
#> Spatial network structure is valid
n01 = st_sfc(st_point(c(0, 0)))
n02 = st_sfc(st_point(c(-1, 2)))
n03 = st_sfc(st_point(c(-1, 3)))
n04 = st_sfc(st_point(c(-1, 4)))
n05 = st_sfc(st_point(c(-2, 1)))
n06 = st_sfc(st_point(c(-2, 3)))
n07 = st_sfc(st_point(c(-2, 4)))
n08 = st_sfc(st_point(c(-3, 2)))
n09 = st_sfc(st_point(c(-3, 3)))
n10 = st_sfc(st_point(c(-3, 4)))
n11 = st_sfc(st_point(c(-4, 2)))
n12 = st_sfc(st_point(c(-4, 4)))
from = c(1, 2, 2, 3, 3, 5, 5, 8, 8, 9, 9)
to = c(5, 3, 6, 4, 7, 2, 8, 9, 11, 10, 12)
nodes = st_as_sf(c(n01, n02, n03, n04, n05, n06, n07, n08, n09, n10, n11, n12))
edges = data.frame(from = from, to = to)
G_2 = sfnetwork(nodes, edges)
#> Checking if spatial network structure is valid...
#> Spatial network structure is valid
G = st_network_join(G_1, G_2) %>%
convert(to_spatial_explicit, .clean = TRUE)
# Create a set of points.
p1 = st_sfc(st_point(c(1, 0.5)))
p2 = st_sfc(st_point(c(2.5, 1.5)))
p3 = st_sfc(st_point(c(1, 2.5)))
p4 = st_sfc(st_point(c(1, 3.5)))
p5 = st_sfc(st_point(c(1.5, 3.5)))
p6 = st_sfc(st_point(c(-1.5, 1.5)))
p7 = st_sfc(st_point(c(-1.5, 2.5)))
p8 = st_sfc(st_point(c(-1.5, 3.5)))
p9 = st_sfc(st_point(c(-1, 0.5)))
# Important!
# Add a column to the points with TRUE values.
# Such that we know which nodes correspond to the points after we inlcude them in the network.
points = st_as_sf(c(p1, p2, p3, p4, p5, p6, p7, p8, p9))
points$is_point = TRUE
plot(G)
plot(st_geometry(points), col = "red", pch = 20, add = TRUE)
# Blend the points into the network.
# Update the is_point column so that all other nodes get a value of FALSE.
G = st_network_blend(G, points) %>%
morph(to_subgraph, is.na(is_point)) %>%
mutate(is_point = FALSE) %>%
unmorph()
#> Warning: st_network_blend assumes attributes are constant over geometries
#> Subsetting by nodes
G
#> # A sfnetwork with 32 nodes and 31 edges
#> #
#> # CRS: NA
#> #
#> # A rooted tree with spatially explicit edges
#> #
#> # Node Data: 32 x 2 (active)
#> # Geometry type: POINT
#> # Dimension: XY
#> # Bounding box: xmin: -4 ymin: 0 xmax: 4 ymax: 4
#> geometry is_point
#> <POINT> <lgl>
#> 1 (1 0.5) TRUE
#> 2 (2 1) FALSE
#> 3 (0 0) FALSE
#> 4 (1 2.5) TRUE
#> 5 (1 3) FALSE
#> 6 (1 2) FALSE
#> # … with 26 more rows
#> #
#> # Edge Data: 31 x 3
#> # Geometry type: LINESTRING
#> # Dimension: XY
#> # Bounding box: xmin: -4 ymin: 0 xmax: 4 ymax: 4
#> from to geometry
#> <int> <int> <LINESTRING>
#> 1 1 2 (1 0.5, 2 1)
#> 2 3 1 (0 0, 1 0.5)
#> 3 4 5 (1 2.5, 1 3)
#> # … with 28 more rows
plot(G)
# Define our own custom grouping function.
# Remember that in the tidygraph style these functions are meant to be
# used inside tidyverse verbs like mutate() and filter(), in which the
# graph that is being worked on is known and not needed as an input to
# the function.
# The graph being worked on can be accessed with .G().
# Its nodes and edges respectively with .N() and .E().
# The mode argument specifies how to route upstreams on the river network.
# If your network is directed upstreams (i.e. starting at the root of the tree):
# --> Use only outbound edges.
# If your network is directed downstreams (i.e. ending at the root of the tree):
# --> Use only inbound edges.
group_custom = function(mode = "out") {
# First we get the node indices of the nodes we want to route from.
# These are:
# --> All points that were blended into the network.
# --> The root node at the start of the network tree.
# Including the root will group all edges that dont have a blended point downstreams.
origins = which(.N()$is_point | with_graph(.G(), node_is_root()))
# Calculate the cost matrix from the origins, to all nodes in the network.
costs = st_network_cost(.G(), from = origins, Inf_as_NaN = TRUE, mode = mode)
# For each node in the network:
# --> Define which of the origins is the first to be reached when travelling downstreams.
# Remember that in the cost matrix:
# --> The origins (the blended points + the root node) are the rows.
# --> The destinations (all nodes in the network) are the columns.
# Hence, we loop over the columns and keep only the minimum cost value per column.
# We should first remove the zeros, which are the cost values from and to the same node.
keep_minimum_cost = function(i) {
i[i == 0] = NaN
if (any(!is.na(i))) i[i != min(i, na.rm = TRUE)] = NaN
i
}
costs = apply(costs, 2, keep_minimum_cost)
# For each origin we know now which nodes are in its group.
# However, we want to know which edges are in the group.
# The cost matrix does not provide that information.
# Finding the paths from the origins to the nodes in its group will give us this.
# Hence, for each origin:
# --> We compute the paths to all nodes in its group.
# --> We extract the edge indices that are part of these paths.
get_edge_indices = function(i) {
orig = origins[i]
dest = which(!is.na(costs[i, ]))
if (length(dest) > 0) {
paths = st_network_paths(.G(), from = orig, to = dest, mode = mode)
edge_idxs = do.call(c, paths$edge_paths)
unique(edge_idxs)
}
}
groups = lapply(seq_len(nrow(costs)), get_edge_indices)
# In tidygraph the largest group always gets group index number 1.
# We can achieve the same by ordering the groups by number of edges.
groups = groups[order(lengths(groups), decreasing = TRUE)]
# Now we can assign a group index to each edge.
edge_idxs = do.call(c, groups)
group_idxs = rep(seq_along(groups), lengths(groups))
# The last thing left to do is to return the group indices in the correct order.
# That is: the order of the edges in the edge table of the network.
group_idxs[order(edge_idxs)]
}
# Group the edges with our custom grouping function.
G = G %>%
activate("edges") %>%
mutate(group = group_custom())
G
#> # A sfnetwork with 32 nodes and 31 edges
#> #
#> # CRS: NA
#> #
#> # A rooted tree with spatially explicit edges
#> #
#> # Edge Data: 31 x 4 (active)
#> # Geometry type: LINESTRING
#> # Dimension: XY
#> # Bounding box: xmin: -4 ymin: 0 xmax: 4 ymax: 4
#> from to geometry group
#> <int> <int> <LINESTRING> <int>
#> 1 1 2 (1 0.5, 2 1) 2
#> 2 3 1 (0 0, 1 0.5) 6
#> 3 4 5 (1 2.5, 1 3) 5
#> 4 6 4 (1 2, 1 2.5) 2
#> 5 6 7 (1 2, 2 3) 2
#> 6 8 9 (1 3.5, 1 4) 7
#> # … with 25 more rows
#> #
#> # Node Data: 32 x 2
#> # Geometry type: POINT
#> # Dimension: XY
#> # Bounding box: xmin: -4 ymin: 0 xmax: 4 ymax: 4
#> geometry is_point
#> <POINT> <lgl>
#> 1 (1 0.5) TRUE
#> 2 (2 1) FALSE
#> 3 (0 0) FALSE
#> # … with 29 more rows
# Plot the results.
nodes = st_as_sf(G, "nodes")
edges = st_as_sf(G, "edges")
edges$group = as.factor(edges$group) # Such that sf uses categorical colors.
plot(st_geometry(edges))
plot(edges["group"], lwd = 4, key.pos = NULL, add = TRUE)
plot(nodes[nodes$is_point, ], pch = 8, add = TRUE)
plot(nodes[!nodes$is_point, ], pch = 20, add = TRUE)
# Reproduce your example.
# Remember to add the is_point column to the created points.
ln <- structure(list(River_ID = c(159, 160, 161, 186, 196), geometry = structure(list(
structure(c(
289924.625, 289924.5313, 289922.9688, 289920.0625,
289915.7499, 289912.7188, 289907.4375, 289905.3438, 289901.1251,
289889, 289888.5, 289887.5938, 289886.5, 289886.4063, 289885.3124,
289884.0938, 289884.0001, 289882.8125, 289881.625, 289878.6875,
289877.9688, 289876.25, 289874.5625, 289874.25, 289872.7188,
289871.2813, 289871.1875, 289870.0313, 289869, 289868.5939,
289867.8436, 289865.8438, 289864.0625, 289862.5939, 289862.375,
289861.5, 289860.7812, 289860.5625, 289859.5313, 289858.375,
289857.7813, 289855.4063, 289854.25, 289850.8749, 289846.4376,
289841.9064, 289836.0625, 289828.1562, 289822.8438, 289816.625,
289812.4376, 289807.9064, 289798.75, 289793.125, 289786.2188,
289781.375, 289777.3124, 289770.0313, 289765.4375, 289762.2188,
289759.25, 289755.5938, 289753.0625, 289747.9687, 289743.7499,
289741.5938, 289739.5, 289736.1874, 289732.75, 289727, 289723.7499,
289719.625, 289715.5626, 289713.7499, 202817.531300001, 202817.2031,
202815.1094, 202812.468699999, 202809.3906, 202806.7656,
202799.7969, 202797.906300001, 202794.093800001, 202783.515699999,
202783.125, 202782.4844, 202781.906300001, 202781.8125, 202781.3594,
202781.093800001, 202780.9999, 202780.5469, 202780, 202777.625,
202777.0469, 202775.718800001, 202774.1875, 202773.906300001,
202772.1875, 202770.4531, 202770.25, 202768.5156, 202766.6719,
202766, 202764.0469, 202759.6719, 202755.8749, 202752.781300001,
202752.1875, 202749.953199999, 202748.297, 202747.906300001,
202746.0625, 202744.2344, 202743.5625, 202740.4375, 202738.8125,
202734.5, 202727.9844, 202723.5625, 202719.1875, 202714.9845,
202713.031300001, 202710.6875, 202710.0469, 202711.406300001,
202714.5626, 202716.9845, 202718.718900001, 202719.5469,
202718.734300001, 202716.4531, 202715.125, 202713.7344, 202712.093800001,
202709.8749, 202708.875, 202709.2655, 202710.7031, 202712.375,
202712.375, 202712.2344, 202711.0469, 202707.906300001, 202705.406300001,
202703.0469, 202701.468800001, 202700.7656
), .Dim = c(
74L,
2L
), class = c("XY", "LINESTRING", "sfg")), structure(c(
289954.375,
289953.5, 289950.6562, 289949.7499, 289949, 289948.125, 289946.0625,
289945.9688, 289944.5313, 289943.4063, 289941.3438, 289939.4375,
289937.4375, 289935.1875, 289932.75, 289930.625, 289928.8125,
289928.25, 289926.7188, 289925.5313, 289925.7813, 289925.625,
289925.4063, 289925.1251, 289924.625, 202872.75, 202872.031400001,
202868.7031, 202867.343699999, 202864.906199999, 202861.515699999,
202858.297, 202854.406300001, 202851.9375, 202849.468800001,
202847.703, 202846.75, 202845.4531, 202843.6719, 202843.0625,
202841.593900001, 202839.7344, 202839.2344, 202838, 202835.9375,
202832.875, 202825.7344, 202822.9531, 202819.4531, 202817.531300001
), .Dim = c(25L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(
290042.6563, 290042.3437, 290041.5313, 290038.4376,
290037.625, 290036.5313, 290035.5313, 290034.8438, 290034.5313,
290033.7188, 290032.9375, 290032.125, 290030.3437, 290030.0313,
290028.625, 290027.5626, 290027.3438, 290026.7188, 290024.5313,
290023.625, 290020.625, 290018.0001, 290014.9375, 290012.0938,
290008.5625, 290004.375, 290000.0001, 289999.875, 289997.625,
289993.7188, 289990.5, 289987.1562, 289985.4063, 289980.375,
289973.3124, 289966.375, 289961.8438, 289959, 289954.375,
202884.0625, 202884.25, 202884.843800001, 202888.4531, 202889.75,
202891.0469, 202892.0469, 202892.656300001, 202892.843800001,
202893.2501, 202893.5469, 202893.656300001, 202893.4531,
202893.4531, 202893.343699999, 202893.093800001, 202893.0469,
202892.843800001, 202891.953199999, 202891.5469, 202889.843800001,
202888.218800001, 202885.1094, 202880.9219, 202877.5625,
202873.968800001, 202872.5469, 202872.5156, 202872.625, 202874.5469,
202876.734300001, 202878.1719, 202877.953199999, 202876.3125,
202873.468800001, 202872.031400001, 202872.906199999, 202873.0781,
202872.75
), .Dim = c(39L, 2L), class = c(
"XY", "LINESTRING",
"sfg"
)), structure(c(
290054.125, 290053.4375, 290052.5313,
290051.625, 290050.0313, 290048.125, 290044.125, 290040.4376,
290039.4375, 290036.9688, 290031.4375, 290027.5312, 290024.8125,
290021.7499, 290020.9688, 290018.3437, 290015, 290010.25,
290006.0313, 290002.4376, 290000.0001, 289999.2187, 289996.6875,
289995.3438, 289994.125, 289991.1875, 289989.2187, 289987.9688,
289986.125, 289980.5313, 289975.0314, 289970.9063, 289968.5625,
289961.0312, 289948.0001, 289939.625, 289933.1563, 289928.3125,
289926.5313, 289924.625, 202835.953199999, 202835.656300001,
202835.4531, 202835.343699999, 202835.5469, 202835.7656,
202836.25, 202836.4531, 202836.5469, 202836.5469, 202835.953199999,
202836.031400001, 202836.625, 202837.7969, 202838.4844, 202839.343699999,
202836.25, 202832.7656, 202832.3125, 202833.4844, 202834.4844,
202834.8125, 202834.2344, 202832.625, 202830.625, 202828.593800001,
202828.968800001, 202831.0625, 202833.2655, 202835.5781,
202838, 202838.906199999, 202839.125, 202836.4531, 202830.781300001,
202827.093800001, 202823.625, 202818.5, 202817.5625, 202817.531300001
), .Dim = c(40L, 2L), class = c("XY", "LINESTRING", "sfg")), structure(c(
290042.625, 290042.0313, 290041.2187, 290040.3125,
290038.4063, 290037.7188, 290035.8125, 290033.7188, 290030.9063,
290028.2187, 290021.5313, 290021.2187, 290014.2188, 290013.4063,
290012.3125, 290010.0625, 290007.9375, 290005.9688, 290004.125,
290000.0001, 289999.4063, 289998.3125, 289997.5312, 289996.8438,
289993.625, 289993.0314, 289989.7188, 289989.3438, 289987.625,
289987.2187, 289984.0313, 289978.125, 289977.9375, 289974.3437,
289972.7188, 289970.9375, 289967.9375, 289965.2187, 289965.1563,
289962.3437, 289960.5313, 289959.1251, 289959.0314, 289959.3438,
289959.4375, 289959.4375, 289959.3438, 289959.2187, 289958.9375,
289958.5313, 289956.125, 289954.375, 202953.781300001, 202952.4844,
202951.281300001, 202950.0781, 202948.1875, 202947.5781,
202945.8749, 202944.281300001, 202941.781300001, 202940.1875,
202936.375, 202936.1875, 202931.968800001, 202931.4844, 202930.875,
202929.093800001, 202927.1094, 202925.031300001, 202922.734300001,
202917.2031, 202916.4375, 202915.2031, 202914.5469, 202914.4531,
202911.4531, 202910.843800001, 202908.0469, 202907.75, 202906.75,
202906.5469, 202904.843800001, 202901.843800001, 202901.75,
202900.0469, 202899.156400001, 202898.0469, 202894.656300001,
202892.0469, 202891.9844, 202889.343699999, 202887.656300001,
202885.75, 202884.5469, 202883.343699999, 202882.5469, 202881.343699999,
202880.0469, 202879.343699999, 202877.656300001, 202876.25,
202874.25, 202872.75
), .Dim = c(52L, 2L), class = c(
"XY",
"LINESTRING", "sfg"
))
), n_empty = 0L, crs = structure(list(
input = "OSGB 1936 / British National Grid", wkt = "PROJCRS[\"OSGB 1936 / British National Grid\",\n BASEGEOGCRS[\"OSGB 1936\",\n DATUM[\"OSGB 1936\",\n ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4277]],\n CONVERSION[\"British National Grid\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",49,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-2,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",400000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",-100000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"unknown\"],\n AREA[\"UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E\"],\n BBOX[49.75,-9.2,61.14,2.88]],\n ID[\"EPSG\",27700]]"
), class = "crs"), class = c(
"sfc_LINESTRING",
"sfc"
), precision = 0, bbox = structure(c(
xmin = 289713.7499,
ymin = 202700.7656, xmax = 290054.125, ymax = 202953.781300001
), class = "bbox"))), row.names = c(NA, -5L), class = c(
"sf",
"data.frame"
), sf_column = "geometry", agr = structure(c(River_ID = NA_integer_), .Label = c(
"constant",
"aggregate", "identity"
), class = "factor"))
pt <- structure(list(lat = c(
202805.8942, 202836.136, 202872.9487,
202905.3284
), lng = c(
289912.0584, 290014.8446, 290001.2364,
289984.9382
), id = 1:4, geometry = structure(list(structure(c(
289912.058400425,
202805.894199679
), class = c("XY", "POINT", "sfg")), structure(c(
290014.844597566,
202836.136003318
), class = c("XY", "POINT", "sfg")), structure(c(
290001.236395958,
202872.948712436
), class = c("XY", "POINT", "sfg")), structure(c(
289984.938209474,
202905.32838227
), class = c("XY", "POINT", "sfg"))), n_empty = 0L, crs = structure(list(
input = "OSGB 1936 / British National Grid", wkt = "PROJCRS[\"OSGB 1936 / British National Grid\",\n BASEGEOGCRS[\"OSGB 1936\",\n DATUM[\"OSGB 1936\",\n ELLIPSOID[\"Airy 1830\",6377563.396,299.3249646,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4277]],\n CONVERSION[\"British National Grid\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",49,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-2,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996012717,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",400000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",-100000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"unknown\"],\n AREA[\"UK - Britain and UKCS 49°46'N to 61°01'N, 7°33'W to 3°33'E\"],\n BBOX[49.75,-9.2,61.14,2.88]],\n ID[\"EPSG\",27700]]"
), class = "crs"), class = c(
"sfc_POINT",
"sfc"
), precision = 0, bbox = structure(c(
xmin = 289912.058400425,
ymin = 202805.894199679, xmax = 290014.844597566, ymax = 202905.32838227
), class = "bbox"))), row.names = c(NA, 4L), class = c(
"sf",
"data.frame"
), sf_column = "geometry", agr = structure(c(
lat = NA_integer_,
lng = NA_integer_,
id = NA_integer_
), .Label = c(
"constant",
"aggregate", "identity"
), class = "factor"))
ln_sfnetwork <- as_sfnetwork(ln)
pt$is_point <- TRUE
ln_sfnetwork <- st_network_blend(ln_sfnetwork, pt["is_point"]) %>%
morph(to_subgraph, is.na(is_point)) %>%
mutate(is_point = FALSE) %>%
unmorph()
#> Warning: st_network_blend assumes attributes are constant over geometries
#> Subsetting by nodes
# NOTE! Your network is directed towards the root of the tree!
# Therefore we route over inbound edges instead of outbound edges to get correct results.
ln_sfnetwork <- ln_sfnetwork %>%
activate("edges") %>%
mutate(group = group_custom(mode = "in"))
# Plot the results.
nodes = st_as_sf(ln_sfnetwork, "nodes")
edges = st_as_sf(ln_sfnetwork, "edges")
edges$group = as.factor(edges$group) # Such that sf uses categorical colors.
plot(st_geometry(edges))
plot(edges["group"], lwd = 4, key.pos = NULL, add = TRUE)
plot(nodes[nodes$is_point, ], pch = 8, add = TRUE)
plot(nodes[!nodes$is_point, ], pch = 20, add = TRUE)
Created on 2021-02-03 by the reprex package (v0.3.0)
- grid.arrange's relative height argument didn't works well in RMarkdown
- How apply a function to every element of a matrix?
- Polar plot in R using ggplot2
- Find total length of time in multiple values of a dataset that overlaps with a given period of time
- Faster than R's `simplify2array`?
- What is the most simple way to label a data.table frequency chart in R?
- How to show dates in a secondary x-axis while the primary x-axis is numerical?
- online compiler for r-markdown document
- Unecessarily large column widths in HTML tables published by R Quarto when the table header contains special characters
- Bootstrapping a vector of results, by group in R
- Speeding Up Monte Carlo Simulations for Weather-Dependent Claims
- Is there an R class which allows to list 2 elements but shows only the first of them by default when called?
- Convert day of week number to weekday name in R
- How do you use "<<-" (scoping assignment) in R?
- Using reactivePoll to monitor and process file changes
- How to find out which graph is most similar to another graph?
- ggplot2: rectangular background behind axis text in circular bar plot
- How to simulate button click using rvest
- CMIP6, automating by month
- Covert dummy variables to single categorical in R?
- Animated reducing bar chart with gganimate
- Create an object class matrix/array from data frame
- How to add a vector of names to `sort_by()` in base R
- Among many subfolders, find CSVs starting with a given string and within ZIPs or not, and merge them while adding their names in a new column with R
- How to add a greek letter to an axis label to a 3D plot in plotly?
- bigrquery + dplyr error when using functions and sourcing from separate R script
- Shiny App Error: /v1/applications/ 400 - Validation Error Execution halted
- Rupee Unicode replacement not working in R
- Mutate multiple variable to create multiple new variables
- How to use the Stata equivalent of a foreach loop in R with mutate()