The goal of slopes is to enable fast, accurate and user friendly calculation longitudinal steepness of linear features such as roads and rivers, based on commonly available input datasets such as road geometries and digital elevation model (DEM) datasets.

Installation

Install the development version from GitHub with:

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

Usage

Load the package in the usual way:

We will also load the sf library:

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0

The minimum data requirements for using the package are elevation points, either as a vector, a matrix or as a digital elevation model (DEM) encoded as a raster dataset. Typically you will also have a geographic object representing the roads or similar 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_segments)
#> [1] "sf"         "tbl_df"     "tbl"        "data.frame"
plot(sf::st_geometry(lisbon_road_segments), add = TRUE)

Calculate the average gradient of each road segment as follows:

lisbon_road_segments$slope = slope_raster(lisbon_road_segments, e = dem_lisbon_raster)
#> [1] TRUE
summary(lisbon_road_segments$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 datast (we find higher correlations on larger datasets):

cor(
  lisbon_road_segments$slope,    # slopes calculates by the slopes package
  lisbon_road_segments$Avg_Slope # slopes calculated by ArcMap's 3D Analyst extension
)
#> [1] 0.9770436

We can now visualise the slopes calculated by the slopes package as follows:

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

# mapview::mapview(lisbon_road_segments["slope"], map.types = "Esri.WorldStreetMap")

Imagine that we want to go from Santa Catarina to the East of the map to the Castelo de Sao Jorge to the West of the map:

mapview::mapview(lisbon_route)

We can convert the lisbon_route object into a 3d linestring object as follows:

lisbon_route_3d = slope_3d(lisbon_route, dem_lisbon_raster)
#> [1] TRUE

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

plot_slope(lisbon_route_3d)

If you do not have a raster dataset representing elevations, you can automatically download them as follows.

lisbon_route_3d_auto = slope_3d(r = lisbon_route)
#> Preparing to download: 12 tiles at zoom = 15 from 
#> https://api.mapbox.com/v4/mapbox.terrain-rgb/
#> [1] TRUE
plot_slope(lisbon_route_3d_auto)

Performance

For this benchmark we will download the following small (< 100 kB) .tif file:

u = "https://github.com/ITSLeeds/slopes/releases/download/0.0.0/dem_lisbon.tif"
if(!file.exists("dem_lisbon.tif")) download.file(u, "dem_lisbon.tif")

A benchmark can reveal how many route gradients can be calculated per second:

e = dem_lisbon_raster
r = lisbon_road_segments
et = terra::rast("dem_lisbon.tif")
res = bench::mark(check = FALSE,
  slope_raster = slope_raster(r, e),
  slope_terra = slope_raster(r, et)
)
res
#> # A tibble: 2 x 6
#>   expression        min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr>   <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 slope_raster   39.3ms   44.5ms      22.6    32.7MB     9.02
#> 2 slope_terra    60.5ms   61.6ms      14.4    29.2MB     4.80

That is approximately

round(res$`itr/sec` * nrow(r))
#> [1] 6113 3900

routes per second using the raster and terra (the default if installed, using RasterLayer and native SpatRaster objects) packages to extract elevation estimates from the raster datasets, respectively.

The message: use the terra package to read-in DEM data for slope extraction if speed is important.

To go faster, you can chose the simple method to gain some speed at the expense of accuracy:

e = dem_lisbon_raster
r = lisbon_road_segments
res = bench::mark(check = FALSE,
  bilinear1 = slope_raster(r, e),
  bilinear2 = slope_raster(r, et),
  simple1 = slope_raster(r, e, method = "simple"),
  simple2 = slope_raster(r, et, method = "simple")
)
res
#> # A tibble: 4 x 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 bilinear1    38.5ms   40.6ms      24.2    32.7MB     2.20
#> 2 bilinear2    61.7ms   63.4ms      15.8    29.2MB     5.25
#> 3 simple1      31.2ms   32.4ms      30.6      29MB     5.09
#> 4 simple2      61.9ms   63.3ms      15.8    29.2MB     5.26