Example: gradients of a road network for a given city • slopes

An example of the demand for data provided by the slopes package is a map showing gradients across São Paulo (Brazil, see image below), with a simplistic classification for cycling difficulty.

using slopes() to create a road network gradient for cycling for São Paulo (Brazil)

This vignette will guide through the production of an interactive slope map for a road network, using slopes, osmextract, sf, stplanr and tmap.

For the convenience of sample, we will use Isle of Wight case, with 384 km². See Other examples below.

This will follow three steps:

Download of road network from OpenStreetMap
Prepare the network
Compute slopes and export the map in html

Extract the OSM network from geofabrik

For this step you may use osmextract package which downloads the most recent information available at OSM (https://download.geofabrik.de/index.html) and converts to GeoPackage (.gpkg), the equivalent to shapefile.

library(dplyr)
library(sf)
# remotes::install_github("ITSLeeds/osmextract")
library(osmextract)

# get the network
iow_osm = oe_get("Isle of Wight", provider = "geofabrik", stringsAsFactors = FALSE, 
                 quiet = FALSE, force_download = TRUE, force_vectortranslate = TRUE) # 7 MB

# filter the major roads
iow_network = iow_osm %>% 
  dplyr::filter(highway %in% c('primary', "primary_link", 'secondary',"secondary_link", 
                               'tertiary', "tertiary_link", "trunk", "trunk_link", 
                               "residential", "cycleway", "living_street", "unclassified", 
                               "motorway", "motorway_link", "pedestrian", "steps", "track")) #remove: "service"

Clean the road network

These are optional steps that give better results, although they may slow down the process since they increase the number of segments present in the network.

Filter the unconnected segments

The rnet_group() function from stplanar package assesses the connectivity of each segment assigns a group number (similar to a clustering process). Then we may filter the main group, the one with more connected segments.

# remotes::install_github("ropensci/stplanr")
library(stplanr)

# filter unconnected roads
iow_network$group = rnet_group(iow_network)
iow_network_clean = iow_network %>% filter(group == 1) # the network with more connected segments

Break the segments on vertices

A very long segment will have an assigned average slope, but a very long segment can be broken into its nodes and have its own slope in each part of the segment. On one hand, we want the segments to break at their nodes. On the other hand, we don’t want artificial nodes to be created where two lines cross, in particular where they have different z levels (eg. brunels: bridges and tunnels).

The rnet_breakup_vertices from stplanr breaks the segments at their inner vertices, preserving the brunels.

iow_network_segments = rnet_breakup_vertices(iow_network_clean)

In this case, there are around 1.6 x the number of segments than in the original network.

Get slope values for each segment

For this case we will use slope_raster() function, to retrieve the z values from a digital elevation model. This raster was obtained from STRM NASA mission.

The SRTM (Shuttle Radar Topography Mission) NASA’s mission provides freely available worldwide DEM, with a resolution of 25 to 30m and with a vertical accuracy of 16m - more. The resolution for USA might be better.

Alternatively, COPERNICUS ESA’s mission also provides freely available DEM for all Europe, with a 25m resolution and a vertical accuracy of 7m - more.

Depending of how large is your road network, you can use elevation_add() function - this will require a valid Mapbox api key.

# Import and plot DEM
u = "https://github.com/U-Shift/Declives-RedeViaria/releases/download/0.2/IsleOfWightNASA_clip.tif"
f = basename(u)
download.file(url = u, destfile = f, mode = "wb")
dem = raster::raster(f)
# res(dem) #27m of resolution
network = iow_network_segments

library(raster)
plot(dem)
plot(sf::st_geometry(network), add = TRUE) #check if they overlay

All the required data is prepared to estimate the road segments’ gradient.

# Get the slope value for each segment (abs), using slopes package
library(slopes)
library(geodist)
network$slope = slope_raster(network, dem)
network$slope = network$slope*100 #percentage
summary(network$slope) #check the values

Half of the road segments in Isle of Wight have a gradient below 3.1%.

We will adopt a simplistic qualitative classification for cycling effort uphill, and compare the number of segments in each class.

# Classify slopes
network$slope_class = network$slope %>%
  cut(
    breaks = c(0, 3, 5, 8, 10, 20, Inf),
    labels = c("0-3: flat", "3-5: mild", "5-8: medium", "8-10: hard", 
               "10-20: extreme", ">20: impossible"),
    right = F
  )
round(prop.table(table(network$slope_class))*100,1)

It means that 49% of the roads are flat or almost flat (0-3%) and about 75% of the roads are easily cyclable (0-5%).

Now let us put this information on a map (see here for interactive version).

# more useful information
network$length = st_length(network)

# make an interactive map
library(tmap)
palredgreen = c("#267300", "#70A800", "#FFAA00", "#E60000", "#A80000", "#730000") #color palette
# tmap_mode("view")
tmap_options(basemaps = leaflet::providers$CartoDB.Positron) #basemap

slopemap =
  tm_shape(network) +
  tm_lines(
    col = "slope_class",
    palette = palredgreen,
    lwd = 2, #line width
    title.col = "Slope [%]",
    popup.vars = c("Highway" = "highway",
                   "Length" = "length",
                  "Slope: " = "slope",
                  "Class: " = "slope_class"),
    popup.format = list(digits = 1),
    # id = "slope"
    id = "name" #if it gets too memory consuming, delete this line
  )

slopemap

#export to html
tmap_save(slopemap, "html/SlopesIoW.html") 

# export information as geopackage
st_write(network, "shapefiles/SlopesIoW.gpkg", append=F)

Result:

Isle of Wight (UK)