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. ropensci/slopes#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.

[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 ***@***.***> 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 <#46 (comment)>, 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!