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

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

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)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"
st_geometry_type(sf_linestring)
#> [1] LINESTRING
#> 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

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/

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)
#>               X         Y    Z L1
#>  [1,] -87064.34 -105506.3 88.0  1
#>  [2,] -87065.47 -105514.3 87.7  1
#>  [3,] -87066.60 -105522.3 86.3  1
#>  [4,] -87067.73 -105530.3 86.3  1
#>  [5,] -87068.86 -105538.2 86.2  1
#>  [6,] -87069.99 -105546.2 84.4  1
#>  [7,] -87075.24 -105548.4 83.1  1
#>  [8,] -87080.48 -105550.5 81.5  1
#>  [9,] -87080.06 -105560.1 80.8  1
#> [10,] -87079.65 -105569.6 79.8  1
#> [11,] -87079.23 -105579.2 78.5  1
#> [12,] -87078.81 -105588.8 77.7  1
#> [13,] -87078.39 -105598.3 76.4  1
#> [14,] -87069.73 -105601.7 80.0  1
#> [15,] -87068.93 -105608.4 79.5  1
#> [16,] -87068.14 -105615.1 78.5  1
#> [17,] -87067.34 -105621.7 77.6  1
#> [18,] -87062.16 -105625.7 79.5  1
#> [19,] -87056.99 -105629.6 79.8  1

You can use the z_ functions to extract such values:

z_value(sf_linestring_xyz) # returns all the elevation values between xy coordinates
#>  [1] 88.0 87.7 86.3 86.3 86.2 84.4 83.1 81.5 80.8 79.8 78.5 77.7 76.4 80.0 79.5
#> [16] 78.5 77.6 79.5 79.8

z_mean(sf_linestring_xyz) # elevation mean value
#> [1] 81.66316
z_min(sf_linestring_xyz) # elevation min value 
#> [1] 88
z_max(sf_linestring_xyz) # elevation max value 
#> [1] 88
z_start(sf_linestring_xyz) # first z
#> [1] 88
z_end(sf_linestring_xyz) # last z
#> [1] 79.8

Step 2: calculate the average slope of the linestring

slope_xyz(sf_linestring_xyz)
#>         1 
#> 0.1394946

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)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"
st_geometry_type(sf_route)
#> [1] LINESTRING
#> 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

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/

Step 2: calculate the average slope of the route

slope_xyz(sf_route_xyz)
#>          1 
#> 0.07681812

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)
#> [1] "RasterLayer"
#> attr(,"package")
#> [1] "raster"
slope_raster(routes = sf_route,
             dem = dem_lisbon_raster)
#>          1 
#> 0.07817098

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.

# 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)
#> [1]  0.05  0.00  0.20 -0.10 -0.05

# 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)
#>  [1] -0.047226259 -0.040883072 -0.025032918 -0.061124557 -0.017447060
#>  [6] -0.062426272 -0.123580541  0.033705378  0.004292243 -0.040360003
#> [11] -0.151671893 -0.182367906  0.409246854 -0.034463974 -0.098406640
#> [16] -0.161798173  0.076261379  0.100654228
slope_distance_mean(dist, elevations)
#> [1] 0.09283052
slope_distance_weighted(dist, elevations)
#> [1] 0.09501323

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.

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
#> [1] "RasterLayer"
#> attr(,"package")
#> [1] "raster"
summary(raster::values(dem_lisbon_raster)) # heights range from 0 to ~100m
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>   0.000   8.598  30.233  33.733  55.691  97.906    4241
raster::plot(dem_lisbon_raster)

# A vector dataset included in the package:
class(lisbon_road_network)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"
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)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> 0.00000 0.01246 0.03534 0.05462 0.08251 0.27583

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
)
#> [1] 0.9770436

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)

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
#> 2518.951 [m]
lisbon_route_segments = stplanr::rnet_breakup_vertices(lisbon_route_xyz)
summary(sf::st_length(lisbon_route_segments)) # mean of 50 m
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>   0.000   0.000   7.631  46.647  65.887 409.479

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)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#> 0.01325 0.04628 0.08120 0.09273 0.12339 0.22033      24

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)
#>      Min.   1st Qu.    Median      Mean   3rd Qu.      Max.      NA's 
#> -0.202564 -0.062971  0.009981  0.005019  0.085559  0.202564        24

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)

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

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.