Morpher that breaks edge if it shares point with other edge

[idshklein]

Hi,
thanks for the great package. it truly gives new options to use R when analyzing spatial networks.
I downloaded data from osm using R osmdata package and tried to turn it to an sfnetwork object.
It worked, but - the network representation was wrong, due to the fact that some linestrings do share points with other linestrings, and therefore should be connected, but these points are not necessarily the start/endpoint of one of the linestrings.
A simple example:

library(tidyverse)
library(sf)
library(sfnetworks)
horz <- st_linestring(matrix(c(0,0,1,0,2,0),ncol=2,byrow = T))
vert <- st_linestring(matrix(c(1,0,1,1),ncol=2,byrow = T))
bind_rows(horz %>% st_sfc() %>% st_sf(),
           vert %>% st_sfc()%>% st_sf(),
           horz %>% st_sfc() %>% st_cast("POINT")%>% st_sf(),
           vert %>% st_sfc() %>% st_cast("POINT")%>% st_sf()) %>% 
  plot()

As can be seen, there should be a node in the intersection between the lines, that splits the horizontal line into to edges.

when converting to sfnetworks and ploting the network, the same plot is generated:

geom <- st_geometrycollection(list(horz,vert)) %>% 
  st_collection_extract("LINESTRING") %>% 
  st_sf() %>% 
  mutate(rn = row_number())
geom_net <- as_sfnetwork(geom)
plot(geom_net)

However, when I check the number of edges, I see that only two edges exist:

geom_net %>% 
  activate(edges) %>% 
  st_as_sf() %>% 
  select(rn) %>% 
  plot()

I bypassed this problem using sf and lwgeom, doing so:

nodes_to_split_by <- geom %>%
  st_cast("POINT") %>% 
  group_by(rn) %>% 
  filter(row_number() == 1| row_number() == max(row_number())) %>% 
  st_geometry() %>% 
  st_sf() %>% 
  distinct() %>% 
  st_snap(geom ,tolerance = 1e-9)
corrected_geom <- st_split(geom$geometry, nodes_to_split_by$geometry) %>% 
  st_collection_extract("LINESTRING") %>% 
  st_sf() %>% 
  mutate(rn = row_number())
geom_net2 <- as_sfnetwork(corrected_geom)
geom_net2 %>% 
  activate(edges) %>% 
  st_as_sf() %>% 
  select(rn) %>% 
  plot()

However, I suggest that a dedicated spatial morpher will be created for this purpose.

[agila5]

Hi @idshklein, and thanks for your reprex. Unfortunately, that's a known issue (see, for example, the second part of my comment here or these issues in stplanr repo: ropensci/stplanr#412 and ropensci/stplanr#416). A temporary fix (based on the github version of stplanr) is the following:

# upgrade packages
remotes::install_github("ropensci/stplanr")
#> Skipping install of 'stplanr' from a github remote, the SHA1 (3d80f79b) has not changed since last install.
#>   Use `force = TRUE` to force installation

# packages
library(dplyr)
library(sfnetworks)
library(sf)
#> Linking to GEOS 3.7.1, GDAL 2.2.2, PROJ 4.9.2
library(stplanr)

# fake data
horz <- st_linestring(matrix(c(0, 0, 1, 0, 2, 0), ncol = 2, byrow = T))
vert <- st_linestring(matrix(c(1, 0, 1, 1), ncol = 2, byrow = T))
my_geom <- st_geometrycollection(list(horz,vert)) %>% 
  st_collection_extract("LINESTRING") %>% 
  st_sf() %>% 
  mutate(rn = row_number())
# create sfnetwork
(my_geom_snet <- as_sfnetwork(my_geom)) # 2 edges since it reflects the input data
#> # An sfnetwork with 4 nodes and 2 edges
#> #
#> # CRS:  NA 
#> #
#> # A rooted forest with 2 trees with spatially explicit edges
#> #
#> # Node Data:     4 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: 0 ymin: 0 xmax: 2 ymax: 1
#>   geometry
#>    <POINT>
#> 1    (0 0)
#> 2    (2 0)
#> 3    (1 0)
#> 4    (1 1)
#> #
#> # Edge Data:     2 x 4
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: 0 ymin: 0 xmax: 2 ymax: 1
#>    from    to        geometry    rn
#>   <int> <int>    <LINESTRING> <int>
#> 1     1     2 (0 0, 1 0, 2 0)     1
#> 2     3     4      (1 0, 1 1)     2

# pre-process input data
my_geom <- rnet_breakup_vertices(my_geom)
#> Splitting rnet object at the shared boundary points.
# create sfnetwork
(my_geom_snet <- as_sfnetwork(my_geom)) # 3 edges now
#> # An sfnetwork with 4 nodes and 3 edges
#> #
#> # CRS:  NA 
#> #
#> # A rooted tree with spatially explicit edges
#> #
#> # Node Data:     4 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: 0 ymin: 0 xmax: 2 ymax: 1
#>   geometry
#>    <POINT>
#> 1    (0 0)
#> 2    (1 0)
#> 3    (2 0)
#> 4    (1 1)
#> #
#> # Edge Data:     3 x 4
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: 0 ymin: 0 xmax: 2 ymax: 1
#>    from    to    rn     geometry
#>   <int> <int> <int> <LINESTRING>
#> 1     1     2     1   (0 0, 1 0)
#> 2     2     3     1   (1 0, 2 0)
#> 3     2     4     2   (1 0, 1 1)

^{Created on 2020-10-05 by the reprex package (v0.3.0)}

The important part here is my_geom <- rnet_breakup_vertices(my_geom) that "breaks" the input LINESTRING at the shared points. Check ?stplanr::rnet_breakup_vertices and this preprint for more details. Unfortunately, as I said in the first paragraph, this solution may not work for very large network objects, and I'm not sure how to change rnet_breakup_vertices so any suggestion/contribution/issue/PR is welcome.

EDIT: I think that this is also one of the steps that could be outlined in a "pre-preprocess" vignette. I think that, for the moment, the code behind that spatial morpher could be taken from stplanr::rnet_breakup_vertices (or if you think about a better approach that's even better. Robin suggested some solutions in ropensci/stplanr#416 (comment) but I didn't check anything after that discussion 😅)

[luukvdmeer]

I agree this is definitely something that should be implemented. In general, the action of splitting an edge by a point located on that edge (some function like st_split_edge(point), see also #27), and then the specific case explained is this issue indeed as a morpher. I will add it to the next milestone (hope to finish a new version this fall).

[luukvdmeer]

Anyone with a good idea for a name for this morpher? I am out of inspiration ;)

[agila5]

Hi @idshklein, if you are still interested in this problem please check ropensci/stplanr#416 (comment). That's the first draft of a new function that could fix the interactions between stplanr/sfnetworks and OSM data (even if I still need to test several other cases and situations). I could try working on a morpher that wraps that function during the next weekends.

[luukvdmeer]

@agila5 I updated the to_spatial_subdivision morpher with a first implementation of the what resulted from the discussion in ropensci/stplanr#416

For sfnetworks object it is more complicated since we also need to join original node attributes back in after splitting the edges. Let me know if you have feedback/improvements

[agila5]

Hi @luukvdmeer! Thank you very much again for helping with this function.

For sfnetworks object it is more complicated since we also need to join original node attributes back in after splitting the edges. Let me know if you have feedback/improvements

I think that the implementation is clear and I cannot think of any particular improvement. I would just add more comments here and there since the algorithm is not trivial (especially for sfnetworks objects), and I think it will be very difficult to modify in a few months. I will work on a PR!

Anyway, if you agree, I would also keep working on the PR on stplanr repo to fix that issue (i.e. reimplement rnet_breakup) using some of the ideas that you improved.

[luukvdmeer]

In the new release this is implemented in the to_spatial_subdivision() morpher. See this vignette for explanation. I think the code is already quite neat and fast, thanks to stplanr.

One point of discussion: edges are only subdivided if an interior point of their geometry is shared with another interior point or endpoint of another edge geometry. That is, no subdivision takes place when edges intersect with another edge at locations where there is no interior point in at least one of them. @agila5 I was wondering if implementing it like that in stplanr was purely to make it easier to implement or also because of conceptual reasons.

[agila5]

One point of discussion: edges are only subdivided if an interior point of their geometry is shared with another interior point or endpoint of another edge geometry. That is, no subdivision takes place when edges intersect with another edge at locations where there is no interior point in at least one of them. @agila5 I was wondering if implementing it like that in stplanr was purely to make it easier to implement or also because of conceptual reasons.

Hi. I decided to "design" the breakup/subdivide algorithm in that way since I usually work with OSM data and I didn't want to create artificial intersections between streets located at different levels (such as an overpass or a bridge). Indeed, the implicit assumption of that algorithm is that all "existing" intersections share one interior-point, while roads located at different levels share no point. I excluded the "split at all intersecting locations" approach since it could create non-existing links between streets at different levels. Does it make sense? Does it answer your question?

[luukvdmeer]

Makes sense! I am just thinking there might be other applications where splitting at all intersections does make sense. Maybe could be implemented with a all = TRUE argument or similar. Needs some more thought, so I will keep this issue open.

[luukvdmeer]

I think our current implementation of to_spatial_subdivision is now mature enough to close this issue. I opened a new issue regarding the subdivision of edges at any crossing, see #134