Slope-sensitive routing

[Robinlovelace]

Please describe the problem.

It could be useful to be able to create weighting profiles that are sensitive to the slope, gradient, of route segments encoded as linestrings. I've recently developed a way to quickly estimate the slope of 1000s of routes based on a raster dataset with elevations: https://github.com/ITSLeeds/slopes

It could be a fun hackathon challenge. Heads-up @tempospena and @mpadge.

Describe how you currently solve the problem, or how attempted to solve it

I would like to calculate shortest paths as below, but one that responds to steepness:

remotes::install_github("itsleeds/slopes")
remotes::install_github("itsleeds/od")
remotes::install_github("ropensci/stplanr")
library(sfnetworks)
library(sf)

r = slopes::lisbon_road_segments
v = sf::st_coordinates(r)
nrow(v)
set.seed(5)
p = v[sample(nrow(v), size = 3), ]
head(p)
l = od::points_to_odl(p[, 1:2], crs = st_crs(r), interzone_only = TRUE)

net = as_sfnetwork(r)
net_t = net %>%
  activate("edges") %>%
  dplyr::mutate(length = sf::st_length(.))
igraph::shortest_paths(graph = net_t, 1, 200)$vpath

Currently I can do routing on a network as follows:

l_start_points = lwgeom::st_startpoint(l)
l_end_points = lwgeom::st_endpoint(l)
r1 = stplanr::route_local(sln = sln, from = l_start_points[1], to = l_end_points[2])
plot(r1$geometry, col = "red", lwd = 5)

# calculate shortest paths
sp = stplanr::route(
  l = l,
  route_fun = stplanr::route_local,
  sln = sln
)

Describe the desired functionality of sfnetworks A clear and concise description of what functionalities you would like to see in the sfnetworks package such that your problem can be solved more conveniently.

To be able to incorporate slope in the route estimate in a controlled way. Some demonstration code, showing 3 shortest routes calculated with the route() function in stplanr are shown below.

remotes::install_github("itsleeds/slopes")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'slopes' from a github remote, the SHA1 (aaf4d36b) has not changed since last install.
#>   Use `force = TRUE` to force installation
remotes::install_github("itsleeds/od")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'od' from a github remote, the SHA1 (d9e37c4a) has not changed since last install.
#>   Use `force = TRUE` to force installation
remotes::install_github("ropensci/stplanr")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'stplanr' from a github remote, the SHA1 (b2b95040) has not changed since last install.
#>   Use `force = TRUE` to force installation
library(sfnetworks)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0

r = slopes::lisbon_road_segments
r = stplanr::overline(r, "Avg_Slope")
#> 2020-05-02 14:05:12 constructing segments
#> 2020-05-02 14:05:12 building geometry
#> 2020-05-02 14:05:12 simplifying geometry
#> 2020-05-02 14:05:12 aggregating flows
#> 2020-05-02 14:05:12 rejoining segments into linestrings
sln = stplanr::SpatialLinesNetwork(r)
#> Warning in SpatialLinesNetwork.sf(r): Graph composed of multiple subgraphs,
#> consider cleaning it with sln_clean_graph().
sln = stplanr::sln_clean_graph(sln)
#> Input sln composed of 4 graphs. Selecting the largest.
nrow(r)
#> [1] 219
nrow(sln@sl) # simple graph
#> [1] 207
v = sf::st_coordinates(sln@sl)
nrow(v)
#> [1] 2367
set.seed(8)
p = v[sample(nrow(v), size = ), ]
p = st_sample(st_convex_hull(st_union(sln@sl)), size = 3)
l = od::points_to_odl(st_coordinates(p), crs = st_crs(r), interzone_only = TRUE)
l$v = 1
l = od::od_oneway(l)
plot(sln@sl$geometry)
plot(p, add = TRUE)

net = as_sfnetwork(sln@sl)
net_t = net %>%
  activate("edges") %>%
  dplyr::mutate(length = sf::st_length(.))
igraph::shortest_paths(graph = net_t, 1, 9)$vpath
#> Warning in igraph::shortest_paths(graph = net_t, 1, 9): At
#> structural_properties.c:745 :Couldn't reach some vertices
#> [[1]]
#> + 0/187 vertices, from 05a320d:

# rnet
r_test = stplanr::sum_network_routes(sln = sln, start = 1, end = 9)
plot(r_test)

# calculate shortest paths
# test with route_local
l_start_points = lwgeom::st_startpoint(l)
l_end_points = lwgeom::st_endpoint(l)
r1 = stplanr::route_local(sln = sln, from = l_start_points[1], to = l_end_points[2])
plot(r1$geometry, col = "red", lwd = 5)

# calculate shortest paths
sp = stplanr::route(
  l = l,
  route_fun = stplanr::route_local,
  sln = sln
)
#> Most common output is sf
#> Loading required namespace: data.table
plot(st_geometry(sln@sl))
plot(st_geometry(l), add = TRUE, lwd = 5)
plot(sp["route_number"], add = TRUE, lwd = rev(sp$route_number) * 3)

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

[Robinlovelace]

It would be great to add sfnetworks as a backend that stplanr can use in the route() function, to go from straight lines to routes directly. Here are the results of what I have so far:

remotes::install_github("itsleeds/slopes")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'slopes' from a github remote, the SHA1 (aaf4d36b) has not changed since last install.
#>   Use `force = TRUE` to force installation
library(sfnetworks)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0

r = slopes::lisbon_road_segments
e = slopes::dem_lisbon_raster
r$slope = slopes::slope_raster(r, e)
plot(r["slope"])

v = sf::st_coordinates(r)
nrow(v)
#> [1] 4365
set.seed(5)
p = v[sample(nrow(v), size = 10), ]
head(p)
#>              X         Y  L1
#> [1,] -87754.50 -106090.5 253
#> [2,] -87615.91 -105512.5 154
#> [3,] -87388.56 -106119.0 118
#> [4,] -87673.75 -105314.2  84
#> [5,] -87303.96 -105826.8 249
#> [6,] -88268.24 -106002.7  54
psf = sf::st_as_sf(
  data.frame(p),
  coords = c("X", "Y")
)
l = stplanr::points2line(p = psf)
plot(l)

net = as_sfnetwork(r)
net_t = net %>%
  activate("edges") %>%
  dplyr::mutate(length = sf::st_length(.))
igraph::shortest_paths(graph = net_t, 1, 200)$vpath
#> [[1]]
#> + 13/205 vertices, from ee62bb3:
#>  [1]   1   2   3   4  10  11  12  13  19 157 156 199 200

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

[mpadge]

@luukvdmeer There is also code in osmdata for elevation estimates which you can freely borrow, along with the elevatr package which is primarily US-based, but likely extensible.

[Robinlovelace]

Those links could be useful, thanks for sharing. We found that method = "bilinear" gives much better results when compared with ArcMap 3d Analyst extension, as partially documented here: https://itsleeds.github.io/slopes/reference/slope_raster.html

Also see here re auto-getting elevation data in the common scenario that a DEM data file isn't on your computer - is elevatr still the best in town? Doesn't seem very frequently maintained. https://github.com/ITSLeeds/slopes/issues/6

[luukvdmeer]

Interesting topic!

I am not sure if I understand exactly what you'd like sfnetworks to provide. It can already convert an sf object containing linestrings into a routable graph. What should come on top is a user-friendly wrapper around igraph shortest path function. But would you also like to see the slope calculation functionality to be part of the package (either by importing the slopes package or by moving the source code from that package into sfnetworks)? Or am I getting that wrong?

[Robinlovelace]

Good questions. The first thing I'd like to see is to be able to reproduce the routing code, specifically

r_test = stplanr::sum_network_routes(sln = sln, start = 1, end = 9)

with sfnetworks. That may already be possible but for some reason I couldn't get it to work on my set-up:

igraph::shortest_paths(graph = net_t, 1, 9)$vpath
#> Warning in igraph::shortest_paths(graph = net_t, 1, 9): At
#> structural_properties.c:745 :Couldn't reach some vertices
#> [[1]]
#> + 0/187 vertices, from 05a320d:

So in some senses this hack idea isn't so much about slope-sensitive routing as it is about routing and modifying weighting profiles with sfnetworks. It's more of an application of generic capabilities, a use case, and opportunity to demonstrate the flexibility of the approach. Do you think I should change the title accordingly, e.g. to

Use case: routing with modifiable network weighting profiles (e.g. to discourage routing on steep sections)

[luukvdmeer]

Do you think I should change the title accordingly

No, I think the slope-sensitive routing is a very interesting hackathon topic on its own! So that is useful as a hackathon labeled issue.

Regarding the example with igraph::shortest_paths. This should indeed work already. I think it is a bug in the sf method of as_sfnetwork(), that seems to create a not well-connected network. I have already an idea of what might cause this bug. I'll take a look at that this week and create an issue for that separately.

[Robinlovelace]

Great, thanks @luukvdmeer. It's certainly an interesting use case that made me delve into stplanr's route() function, I was sure I had found a bug there also. I'm working with others including @temospena and @joeytalbot on projects that use gradient in routing, although we use routing services such as CycleStreet.net rather than igraph on a local network. Will be interesting to see if it's possible to do routing on a larger scale using the sf/igraph approach if sfnetworks/stplanr. If anyone is interested in this as a hackathon topic, ideas and other datasets / use cases are very welcome in comments below.

[luukvdmeer]

@Robinlovelace

That may already be possible but for some reason I couldn't get it to work on my set-up

In the end I found out the issue was not a bug, but a design difference between stplanr and sfnetworks. When creating an sfnetwork object, by default directed = TRUE. Hence, it creates a directed graph by default. From what I've seen, stplanr by default creates an undirected graph.

When changing the net object to an undirected graph, the results are the same for sfnetworks and stplanr.

Also, when calculating shortest paths with sfnetworks, we have to define which weights to use, since there is no weight attribute like in stplanr.

library(sfnetworks)

# Create an sf object with slopes.
r = slopes::lisbon_road_segments
r = stplanr::overline(r, "Avg_Slope")

# Create and clean SpatialLinesNetwork
sln = stplanr::SpatialLinesNetwork(r)
sln = stplanr::sln_clean_graph(sln)

# Convert to undirected sfnetwork
net = as_sfnetwork(sln@sl, directed = FALSE)
  activate("edges") %>%
  dplyr::mutate(length = sf::st_length(.))

# Shortest path with SpatialLinesNetwork
igraph::shortest_paths(graph = sln@g, 1, 9)$vpath

#> [[1]]
#> + 14/187 vertices, from ed56407:
#>  [1]   1   2  44 108  87  60  20 119  52  51 106   8   7   9

# Shortest path with sfnetwork
weights = net %>% activate('edges') %>% dplyr::pull(length)
igraph::shortest_paths(graph = net, 1, 9, weights = weights)$vpath

#> [[1]]
#> + 14/187 vertices, from 8284e6c:
#>  [1]   1   2  44 108  87  60  20 119  52  51 106   8   7   9

Findings:

• There is a need for a as_sfnetwork.SpatialLinesNetwork method, to directly convert a SpatialLinesNetwork into a sfnetwork, without having to worry about directedness etc.
• There is a need for a st_shortest_paths function in sfnetworks, being a user-friendly wrapper around the igraph function, and without the need to first pull the weights column.

[Robinlovelace]

Great discoveries @luukvdmeer, agree about the conversion functions.

[agila5]

In the end I found out the issue was not a bug, but a design difference between stplanr and sfnetworks. When creating an sfnetwork object, by default directed = TRUE. Hence, it creates a directed graph by default. From what I've seen, stplanr by default creates an undirected graph.

Would you reconsider setting the default value of directed to FALSE? IMO it's not a proper default since "directed" routing should be managed with more planned and careful options. The problem is that at the moment we can't event estimate something so simple:

# packages
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(sfnetworks)
library(igraph)
#> 
#> Attaching package: 'igraph'
#> The following objects are masked from 'package:stats':
#> 
#>     decompose, spectrum
#> The following object is masked from 'package:base':
#> 
#>     union

# Create a sf object with LINESTRING geometry
es <- st_sf(
  data.frame(ID = c(1, 2)), 
  geometry = st_sfc(
    st_linestring(rbind(c(0, 0), c(0.5, 0), c(1, 0))), 
    st_linestring(rbind(c(1, 0), c(1, 0.5), c(1, 1)))
  )
)
plot(es, col = sf.colors(2, categorical = TRUE), lwd = 2)

# transform into sfnetwork
es_sfnetworks <- as_sfnetwork(es)
es_sfnetworks
#> # A sfnetwork with 3 nodes and 2 edges
#> #
#> # CRS:  NA 
#> #
#> # A rooted tree
#> # and spatially explicit edges
#> #
#> # Node Data:     3 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: 0 ymin: 0 xmax: 1 ymax: 1
#>   geometry
#>    <POINT>
#> 1    (0 0)
#> 2    (1 0)
#> 3    (1 1)
#> #
#> # Edge Data:     2 x 4
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: 0 ymin: 0 xmax: 1 ymax: 1
#>    from    to    ID          geometry
#>   <int> <int> <dbl>      <LINESTRING>
#> 1     1     2     1 (0 0, 0.5 0, 1 0)
#> 2     2     3     2 (1 0, 1 0.5, 1 1)

# The first node is at (0, 0) and the last point is at (1, 1)

# Estimate the shortest path between (0, 0) and (2, 1)
shortest_paths(es_sfnetworks, 1, 3, output = "both")
#> $vpath
#> $vpath[[1]]
#> + 3/3 vertices, from c32f054:
#> [1] 1 2 3
#> 
#> 
#> $epath
#> $epath[[1]]
#> + 2/2 edges from c32f054:
#> [1] 1->2 2->3
#> 
#> 
#> $predecessors
#> NULL
#> 
#> $inbound_edges
#> NULL
shortest_paths(es_sfnetworks, 3, 1, output = "both")
#> Warning in shortest_paths(es_sfnetworks, 3, 1, output = "both"): At
#> structural_properties.c:745 :Couldn't reach some vertices
#> $vpath
#> $vpath[[1]]
#> + 0/3 vertices, from c32f054:
#> 
#> 
#> $epath
#> $epath[[1]]
#> + 0/2 edges from c32f054:
#> 
#> 
#> $predecessors
#> NULL
#> 
#> $inbound_edges
#> NULL

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

I know that the example is trivial (and maybe factious), but I think it can represent why I would modify the default argument to FALSE.

[mpadge]

I agree with @agila5 - a default of directed = TRUE requires an awful lot of work, and a completely different approach to routing. Defaults should almost always be directed = FALSE.

[luukvdmeer]

Thanks for the input! Having directed = TRUE was never really a well-considered choice that I made myself, but rather an inheritance from tidygraph and igraph. An sfnetwork subclasses a tbl_graph, which has directed = TRUE as default, probably because a tbl_graph subclasses again igraph, which also has directed = TRUE as default.

Since sfnetworks is meant as a general purpose spatial network analysis package, and not specifically focused on routing, I would argue that it makes most sense to stick with the defaults of tidygraph and igraph, on whose shoulders the package is largely standing. Changing that default might cause confusion, I think. But of course I am very open to change if the majority of users agree that the current default is not intuitive. Let's poll this during the hackathon!

[mpadge]

All very fair points @luukvdmeer, and definitely something that would be very worthwhile raising as a general discussion point. Just to have this for the record here, I note that the Simple Features standard is in no way intended to, nor really capable of, holding or representing concepts of "directedness". A LINESTRING is nothing more than a direction-agnostic ordering a points along a line, with no geometric operations (for example via GDAL) paying any attention to the direction of ordering. So sf is and will always remain very strictly hard-coded for undirected networks (when built from LINESTRING objects). That ought in turn logically suggest that any graph built from Simple Features Standards-compliant objects should therefore always use directed = FALSE as a default.

[Robinlovelace]

Earlier in the year there was a function added to change the direction of LINESTRINGs to sf: https://github.com/r-spatial/sf/issues/1246

Not sure if that can be used to infer the direction of roads on a complex street network, but it seems to work:

remotes::install_github("luukvdmeer/sfnetworks")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'sfnetworks' from a github remote, the SHA1 (a5e345af) has not changed since last install.
#>   Use `force = TRUE` to force installation
r = sfnetworks::roxel[3, ]
class(r)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"
sf::st_coordinates(r)
#>             X        Y L1
#> [1,] 7.532709 51.95209  1
#> [2,] 7.532823 51.95239  1
#> [3,] 7.532869 51.95257  1
r_rev = sf::st_reverse(r)
sf::st_coordinates(r_rev)
#>             X        Y L1
#> [1,] 7.532869 51.95257  1
#> [2,] 7.532823 51.95239  1
#> [3,] 7.532709 51.95209  1

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

[luukvdmeer]

I have hidden the "directed" discussion (which moved to #50), to keep this clean.

@Robinlovelace I was thinking it might also be an option to use the Z-coordinate of a geometry for slope sensitive routing? I have never worked with that, but might be interesting.

[Robinlovelace]

Yes, currently the slopes package calculates slopes based on the Z dimension, using the function slopes_3d(). An interesting element of this topic is that slopes can be calculated between each coordinate, but slope dependent routing depends on attributes at the linestring level. @temospena and I have talked about using different measures of slope, including the mean, median, max and 90 percentile slope to take slope into account in a more sophisticated way than just the 'mean gradient' metric used by some approaches such as the Propensity to Cycle Tool.

[Robinlovelace]

Place where we can put code/data/tests for this: https://github.com/sfnetworks/sloperouting

If you want to get involved in this hack let me know and I can add you to the repo!

[temospena]

cool!

On Tue, Jun 16, 2020 at 11:52 AM Robin notifications@github.com wrote:

Place where we can put code/data/tests for this: https://github.com/sfnetworks/sloperouting

If you want to get involved in this hack let me know and I can add you to the repo!

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/luukvdmeer/sfnetworks/issues/46#issuecomment-644689389, or unsubscribe https://github.com/notifications/unsubscribe-auth/AJKLUXXFTQVSL55I6FWULS3RW5FHBANCNFSM4MXFHS5Q .

[Robinlovelace]

Great just added you @temospena if anyone else wants to get involved let us know!

[Robinlovelace]

Starting point for the hack topic:

library(sfnetworks)
library(tidygraph)
#> 
#> Attaching package: 'tidygraph'
#> The following object is masked from 'package:stats':
#> 
#>     filter
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(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
r = slopes::lisbon_road_segments
sf::st_is_longlat(r)
#> [1] FALSE
class(r)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"
names(r)
#>  [1] "OBJECTID"   "fid_1"      "gradient_s" "Shape_Leng" "Z_Min"     
#>  [6] "Z_Max"      "Z_Mean"     "SLength"    "Min_Slope"  "Max_Slope" 
#> [11] "Avg_Slope"  "z0"         "z1"         "gradverifi" "query"     
#> [16] "lat"        "lon"        "lat_min"    "lat_max"    "lon_min"   
#> [21] "lon_max"    "bbox"       "geom"
plot(r["Avg_Slope"])

net = as_sfnetwork(r)
p1 = net %>%  
  activate(nodes) %>%  
  st_as_sf() %>%  
  slice(1)  
p2 = net %>%  
  activate(nodes) %>%  
  st_as_sf() %>%  
  slice(9)
mapview::mapview(p1) + mapview::mapview(p2)

path1 = net %>%  
  activate("edges") %>%  
  mutate(weight = edge_length()) %>%  
  convert(to_spatial_shortest_paths, p1, p2)
plot(path1)

mapview::mapview(st_as_sf(path1)) +
  mapview::mapview(p1) + mapview::mapview(p2)

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

[luukvdmeer]

Since the hackathon is now behind us, I will close the hackathon issues. But I look forward to new specific feature request issues that resulted from the hackathon!