Example: gradients of a road network for a given citySource:
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.
This vignette will guide through the production of an interactive slope map for a road network, using
For the convenience of sample, we will use Isle of Wight case, with 384 km2. 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.
# 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
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.
# 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).
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)