This package calculates longitudinal steepness of linear features such as roads and rivers, based on two main inputs: vector linestring geometries and raster digital elevation model (DEM) datasets.
Install the development version from GitHub with:
# install.packages("remotes") remotes::install_github("itsleeds/slopes")
Load the package in the usual way:
We will also load the
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 #>  "RasterLayer" #> attr(,"package") #>  "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) #>  "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, dem = dem_lisbon_raster) 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 dataset (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 ) #>  0.9770436
We can now visualise the slopes calculated by the
slopes package as follows:
# 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:
library(tmap) tmap_mode("view") #> tmap mode set to interactive viewing qtm(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)
We can now visualise the elevation profile of the route as follows:
If you do not have a raster dataset representing elevations, you can automatically download them as follows.
lisbon_route_3d_auto = slope_3d(lisbon_route) #> Preparing to download: 12 tiles at zoom = 15 from #> https://api.mapbox.com/v4/mapbox.terrain-rgb/ plot_slope(lisbon_route_3d_auto)