# Example: gradients of a road network for a given city

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 km2. See Other examples below.

This will follow three steps:

1. Download of road network from OpenStreetMap
2. Prepare the network
3. 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)

## Other examples

- São Paulo (Brazil)
- Lisbon (Portugal)
- Oporto (Portugal)
- Leeds (UK)
- Zurich (CH)
- Paris (FR)