In ITSLeeds/slopes: Calculate Slopes of Roads, Rivers and Trajectories

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Welcome to the slopes vignette, a type of long-form documentation/article that introduces the core functions and functionality of the slopes package.

Installation

You can install the released version of slopes from CRAN with:

install.packages("slopes")

Install the development version from GitHub with:

# install.packages("remotes")
remotes::install_github("ropensci/slopes")

Installation for DEM downloads

If you do not already have DEM data and want to make use of the package's ability to download them using the ceramic package, install the package with suggested dependencies, as follows:

# install.packages("remotes")
remotes::install_github("ropensci/slopes", dependencies = "Suggests")

Furthermore, you will need to add a MapBox API key to be able to get DEM datasets, by signing up and registering for a key at https://account.mapbox.com/access-tokens/ and then following these steps:

usethis::edit_r_environ()
MAPBOX_API_KEY=xxxxx # replace XXX with your api key

Functions

Elevation

• elevation_add() Take a linestring and add a third dimension (z) to its coordinates
• elevation_get() Get elevation data from hosted maptile services (returns a raster)

• elevation_extract() Extract elevations from coordinates

• z_value() retrieves elevation values for each z (as vector of sequential vertices)

• z_start() retrieves the elevation value of the first linestring vertice
• z_end() retrieves the elevation value of the last linestring vertice
• z_mean() retrieves the elevation mean value
• z_max() retrieves the elevation max value
• z_min() retrieves the elevation min value

Distance

• sequential_dist() Calculate the sequential distances between sequential coordinate pairs

Slope

• slope_vector() calculates the slopes associated with consecutive elements in one dimensional distance and associated elevations.
• slope_distance() calculates the slopes associated with consecutive distances and elevations.
• slope_distance_mean() calculates the mean average slopes associated with consecutive distances and elevations.

• slope_distance_weighted() calculates the slopes associated with consecutive distances and elevations, with the mean value associated with each set of distance/elevation vectors weighted in proportion to the distance between each elevation measurement, so longer sections have proportionally more influence on the resulting gradient estimate.

• slope_raster() Calculate slopes of linestrings based on local raster map

• slope_matrix() Calculate the gradient of line segments from a 3D matrix of coordinates

• slope_matrix_weighted() Calculate the weighted gradient of line segments from a 3D matrix of coordinates

• slope_xyz() Calculates the slope associated with linestrings that have xyz coordinates

Plot

• plot_dz() Plot a digital elevation profile based on xyz data
• plot_slope() Plots the slope profile associated with a linestring with base R graphics

Package datasets

The slopes package comes with some datasets to play with:

Linestrings:

• lisbon_road_segment: a single road segment in Lisbon (XY)
• lisbon_route: a route with some variation in elevation in Lisbon (XY)
• cyclestreets_route: a bike route in Leeds (XY)

Road network:

• lisbon_road_network: a sample of road segments in downtown Lisbon
• magnolia_xy: a sample of road segments in center Seattle, in the Magnolia neighborhood

Digital elevation model (DEM):

• dem_lisbon_raster a DEM of downtown Lisbon (EPSG:3763)

Examples

Load the package in the usual way. We will also load the sf library:

library(slopes)
library(sf)

The minimum input data requirement for using the package is an sf object containing LINESTRING geometries.

You can also create sf objects from a matrix of coordinates, as illustrated below (don't worry about the details for now, you can read up on how all this works in the sf package documentation):

m = st_coordinates(sf::st_transform(lisbon_road_segment, 4326))
s = seq(from = 1, to = nrow(m), length.out = 4)
round(m[s, 1:2], 5)
dput(round(m[s, 1], 4))
dput(round(m[s, 2], 4))

m = cbind(
  c(-9.1333, -9.134, -9.13),
  c(38.714, 38.712, 38.710)
)
sf_linestring = sf::st_sf(
  data.frame(id = 1),
  geometry = st_sfc(st_linestring(m)),
  crs = 4326
)
class(sf_linestring)
st_geometry_type(sf_linestring)

maybe remove this? or add step 1 and step 2 again.

Single road segment + no DEM

You can check your input dataset is suitable with the functions class() from base R and st_geometry_type() from the sf package, as demonstrated below on the example object lisbon_road_segment that is contained within the package:

sf_linestring = lisbon_road_segment
class(sf_linestring)
st_geometry_type(sf_linestring)

A quick way of testing if your object can have slopes calculated for it is to plot it in an interactive map and to check that underneath the object there is indeed terrain that will give the linestrings gradient:

library(tmap)
tmap_mode("view")
tm_shape(sf_linestring) +
  tm_lines(lwd = 5) +
  tm_basemap(leaflet::providers$Esri.WorldTopoMap)

Imagine you want to calculate the gradient of the route shown above. You can do this as a two step process as follows.

Step 1: add elevations to each coordinate in the linestring (requires a MapBox API key):

sf_linestring_xyz = elevation_add(sf_linestring) # dem = NULL
#> Loading required namespace: ceramic
#> Preparing to download: 9 tiles at zoom = 18 from 
#> https://api.mapbox.com/v4/mapbox.terrain-rgb/

# note: the following should be TRUE
# identical(sf_linestring_xyz, lisbon_road_segment_xyz_mapbox)
sf_linestring_xyz = lisbon_road_segment_xyz_mapbox

With the argument dem = NULL, the function downloads the necessary elevation information from Mapbox. You can use this argument with a local digital elevation model (dem = ...).

You can check the elevations added to the new sf_linestring_xyz object by printing its coordinates, as follows (note the new Z column that goes from above 87 m above sea level to only 79 m in a short distance).

st_coordinates(sf_linestring_xyz)

You can use the z_ functions to extract such values:

z_value(sf_linestring_xyz) # returns all the elevation values between xy coordinates

z_mean(sf_linestring_xyz) # elevation mean value
z_min(sf_linestring_xyz) # elevation min value 
z_max(sf_linestring_xyz) # elevation max value 
z_start(sf_linestring_xyz) # first z
z_end(sf_linestring_xyz) # last z

Step 2: calculate the average slope of the linestring

slope_xyz(sf_linestring_xyz)

The result, just over 0.2, tells us that it's quite a steep slope: a 21% gradient on average.

Route + available DEM

Using the slopes package we can estimate the gradient of individual road segments. When these segments are combined into routes, we then need a means of assessing the hilliness of the entire route. A range of indices can be used to represent route hilliness. The choice of which index is most appropriate may be context dependent (see the introducion to slopes vignette).

Again, let us use the same function with a entire route, lisbon_route, also available in the package:

sf_route = lisbon_route
class(sf_route)
st_geometry_type(sf_route)

tm_shape(sf_route) +
  tm_lines(lwd = 3) +
  tm_basemap(leaflet::providers$Esri.WorldTopoMap)

Step 1: add elevations to each coordinate in the route:

sf_route_xyz = elevation_add(sf_route)
#> Loading required namespace: ceramic
#> Preparing to download: 12 tiles at zoom = 15 from 
#> https://api.mapbox.com/v4/mapbox.terrain-rgb/

# note: the following should be TRUE
# identical(sf_route_xyz, lisbon_road_segment_xyz_mapbox)
sf_route_xyz = lisbon_route_xyz_mapbox

Step 2: calculate the average slope of the route

slope_xyz(sf_route_xyz)

The result shows a 7.7% gradient on average.

Now, if you already have a DEM, you can calculate the slopes directly as follows, with slope_raster():

class(dem_lisbon_raster)
slope_raster(routes = sf_route,
             dem = dem_lisbon_raster)

The result shows a 7.8% gradient on average. As you can see, the retrieved result from elevation information available in Mapbox and in this Digital Elevation Model, is quite similar. (See more about these differences in Verification of slopes.)

Route with xyz coordinates

If your linestring object already has X, Y and Z coordinates (e.g. from a GPS device), you can use the slope_ functions directly.

#not to use like this... it would ge good to have a gps example to demonstrate

slope_vector(sf_route_xyz)
slope_distance(sf_route_xyz)
slope_distance_mean(sf_route_xyz)
slope_distance_weighted(sf_route_xyz)

slope_vector(sf_linestring_xyz)
slope_distance(sf_linestring_xyz)
slope_distance_mean(sf_linestring_xyz)
slope_distance_weighted(sf_linestring_xyz)

# for a line xz
x = c(0, 2, 3, 4, 5, 9)
elevations = c(1, 2, 2, 4, 3, 1) / 10
slope_vector(x, elevations)

# for a path xyz
xy = st_coordinates(sf_linestring)
dist = sequential_dist(xy, lonlat = FALSE)
elevations = elevation_extract(xy, dem_lisbon_raster)

slope_distance(dist, elevations)
slope_distance_mean(dist, elevations)
slope_distance_weighted(dist, elevations)

In any case, to use the slopes package you need elevation points, either as a vector, a matrix or as a digital elevation model (DEM) encoded as a raster dataset.

Calculating and plotting gradients

Road network

Typical use cases for the package are calculating the slopes of geographic objects representing roads or other linear features. These two types of input data are represented in the code output and plot below.

# A raster dataset included in the package:
class(dem_lisbon_raster) # digital elevation model
summary(raster::values(dem_lisbon_raster)) # heights range from 0 to ~100m
raster::plot(dem_lisbon_raster)

# A vector dataset included in the package:
class(lisbon_road_network)
plot(sf::st_geometry(lisbon_road_network), add = TRUE)

Calculate the average gradient of each road segment as follows:

lisbon_road_network$slope = slope_raster(lisbon_road_network, dem = dem_lisbon_raster)
summary(lisbon_road_network$slope)

This created a new column, slope that represents the average, distance weighted slope associated with each road segment. The units represent the percentage incline, that is the change in elevation divided by distance. The summary of the result tells us that the average gradient of slopes in the example data is just over 5%.

This result is equivalent to that returned by ESRI's Slope_3d() in the 3D Analyst extension, with a correlation between the ArcMap implementation and our implementation of more than 0.95 on our test dataset (we find higher correlations on larger datasets - see the verification of slopes:

cor(
  lisbon_road_network$slope,    # slopes calculates by the slopes package
  lisbon_road_network$Avg_Slope # slopes calculated by ArcMap's 3D Analyst extension
)

We can now visualise the average slopes of each route calculated by the slopes package as follows:

raster::plot(dem_lisbon_raster)
plot(lisbon_road_network["slope"], add = TRUE, lwd = 5)

Elevation profile

Taking the first route example, imagine that we want to go from from the Santa Catarina area in the East of the map to the Castelo de São Jorge in the West. This route goes down a valley and up the other side:

# library(tmap)
# tmap_mode("view")
qtm(lisbon_route)

# Removed because it's not rendering in RMarkdown
mapview::mapview(lisbon_road_network["slope"], map.types = "Esri.WorldStreetMap")
mapview::mapview(lisbon_route)

We can convert the lisbon_route object into a 3d linestring object with X, Y and Z coordinates, using the elevation values stored in the DEM, as follows:

lisbon_route_xyz = elevation_add(lisbon_route, dem_lisbon_raster) 

We can now visualise the elevation profile of the route as follows:

plot_slope(lisbon_route_xyz)

Splitting the network

The lisbon_route_xyz example is useful but often you will want to calculate the slopes not of an entire route (in this case one that is 2.5 km long) but of segments. There are various ways to split segements, including using algorithms from other packages or GIS programs, but here we'll use the stplanr function rnet_breakup_vertices() (see vignette("roadnetworkcycling") for an example of this function working on a large road network):

sf::st_length(lisbon_route_xyz) # check route length: 2.5 km
lisbon_route_segments = stplanr::rnet_breakup_vertices(lisbon_route_xyz)
summary(sf::st_length(lisbon_route_segments)) # mean of 50 m

We can now calculate the slope for each of these segments.

lisbon_route_segments$slope = slope_xyz(lisbon_route_segments)
summary(lisbon_route_segments$slope)

Directed slopes

The route has a direction that is implicit in the order of the vertices and segments. From the perspective of someone travelling along the route, the slopes have a direction which is important: it's easier to go uphill than downhill. To calculate the slopes with direction, add the directed argument as follows.

lisbon_route_segments$slope_directed = slope_xyz(lisbon_route_segments, directed = TRUE)
summary(lisbon_route_segments$slope_directed)

Plotting the directed and undirected slopes side-by-side shows the importance of considering slope direction for route planning, which may want to avoid steep hills going uphill but not downhill for certain types of travel, for example.

breaks = c(0, 3, 5, 8, 10, 20, 50)
breaks_proportion = breaks / 100
breaks_directed = c(-rev(breaks_proportion), (breaks_proportion[-1]))
plot(lisbon_route_segments["slope"], breaks = breaks_proportion)
plot(lisbon_route_segments["slope_directed"], breaks = breaks_directed)

# test code
z = sf::st_make_grid(lisbon_route_xyz, cellsize = 100)
sampled_points = sf::st_line_sample(lisbon_route_xyz, n = 30)
points_sf = sf::st_sf(geometry = sf::st_cast(sampled_points, "POINT"))
plot(points_sf)
lisbon_route_segments = stplanr::route_split(lisbon_route_xyz, p = points_sf[2:3, ])
lisbon_route_segments = stplanr::rnet_breakup_vertices(lisbon_route_xyz)
library(tmap)
tmap_mode("view")
qtm(lisbon_route_segments, lines.lwd = 9, lines.col = 1:nrow(lisbon_route_segments))
plot(lisbon_route_segments, col = 1:nrow(lisbon_route_segments))

# Test: try using QGIS
remotes::install_github("paleolimbot/qgisprocess")
library(qgisprocess)
qgis_configure()
algorithms = qgis_algorithms()
View(algorithms)
result = qgis_run_algorithm(
  algorithm = "grass7:v.split",
  INPUT = lisbon_route_xyz,
  LENGTH = 500
  )
route_segments = sf::st_read(result$OUTPUT)
route_segments
plot(lisbon_route_xyz$geometry)
plot(route_segments$geom, add = T, lwd = 3)
mapview::mapview(route_segments)

Using elevation_add() with and without a dem = argument

If you do not have a raster dataset representing elevations, you can automatically download them by omitting the argument dem = NULL (a step that is automatically done in the function elevation_add() shown in the basic example above, results of the subsequent code chunk not shown):

dem_mapbox = elevation_get(lisbon_route)
lisbon_road_proj = st_transform(lisbon_route, raster::crs(dem_mapbox))
lisbon_route_xyz_mapbox = elevation_add(lisbon_road_proj, dem = dem_mapbox)
plot_slope(lisbon_route_xyz_mapbox)

As outlined in the basic example above this can be done more concisely, as:

lisbon_route_xyz_auto = elevation_add(lisbon_route) #dem = NULL

lisbon_route_xyz_auto = lisbon_route_xyz_mapbox

plot_slope(lisbon_route_xyz_auto)

Note that the elevations shown in both plots differ, since the first is based on DEM elevation available, and the second is based in Mapbox elevation.

Commulative elevation change

The following example calculate the elevations of a route in Leeds, and plots its commutative sum along the route (not evaluated).

cyclestreets_xyz = elevation_add(cyclestreets_route) 
plot_slope(cyclestreets_xyz)
plot(cumsum(cyclestreets_xyz$distances), cumsum(cyclestreets_xyz$elevation_change))

See more in vignettes

• slopes package
• An introduction to slopes
• Reproducible example: gradients of a road network for a given city
• Verification of slopes
• Benchmarking slopes calculation

ITSLeeds/slopes documentation built on Oct. 13, 2024, 3:54 a.m.