diff --git a/R/slopes.R b/R/slopes.R index dfb7ba5..c544302 100644 --- a/R/slopes.R +++ b/R/slopes.R @@ -154,7 +154,7 @@ slope_matrices = function(m_xyz_split, fun = slope_matrix_weighted, ...) { #' cor(routes$Avg_Slope, s) slope_raster = function( routes, - dem = NULL, + dem, lonlat = sf::st_is_longlat(routes), method = "bilinear", fun = slope_matrix_weighted, diff --git a/README.Rmd b/README.Rmd index b68b6ec..ff4fbbb 100644 --- a/README.Rmd +++ b/README.Rmd @@ -21,7 +21,13 @@ knitr::opts_chunk$set( [![Status at rOpenSci Software Peer Review](https://badges.ropensci.org/420_status.svg)](https://github.com/ropensci/software-review/issues/420) -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. +The **slopes** R package calculates the slope (longitudinal steepness, also known as gradient) of linear features such as roads and rivers, based on two main inputs: + +- [vector](https://geocompr.robinlovelace.net/spatial-class.html#vector-data) linestring geometries defined by classes in the [`sf` package](https://r-spatial.github.io/sf/) +- [raster](https://geocompr.robinlovelace.net/spatial-class.html#raster-data) objects with pixel values reporting average height, commonly known as digital elevation model (DEM) datasets, defined by classes in the [`raster`](https://cran.r-project.org/package=raster) or more recent [`terra`](https://rspatial.org/terra) packages + +This README covers installation and basic usage. +For more information about slopes and how to use the package to calculate them, see the [slopes vignette](https://itsleeds.github.io/slopes/). ## Installation @@ -53,7 +59,7 @@ usethis::edit_r_environ() MAPBOX_API_KEY=xxxxx ``` -## Usage +## Basic example Load the package in the usual way: @@ -67,8 +73,87 @@ We will also load the `sf` library: library(sf) ``` -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. +The minimum input data requirement for using the package is an `sf` object containing LINESTRING geometries. +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 object `lisbon_road_segment` that is contained within the package: + +```{r} +class(lisbon_road_segment) +st_geometry_type(lisbon_road_segment) +``` + +Don't worry if you don't yet have your linear features in this class: you can read-in data from a wide range of formats into an `sf` object. +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](https://r-spatial.github.io/sf/articles/sf1.html)): + +```{r, eval=FALSE, echo=FALSE} +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)) +``` + +```{r} +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) +``` + +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: + +```{r linestringmap} +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 imaginary route shown above, starting at the Castelo de São Jorge and descending towards the coast. + +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): + +```{r} +sf_linestring_xyz = elevation_add(sf_linestring) +``` + +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 90 m above sea level to only 24 m in a short distance): + +```{r} +st_coordinates(sf_linestring_xyz) +``` + +Step 2: calculate the average slope of the linestring + +```{r} +slope_xyz(sf_linestring_xyz) +``` + +The result, just over 0.1, tells us that it's quite a steep slope, a 10% gradient *on average*. + +If you already have a DEM, you can calculate the slopes directly as follows, with `slope_raster()`: + +```{r} +class(dem_lisbon_raster) +sf_linestring_proj = st_transform(sf_linestring, st_crs(lisbon_road_segment)) +slope_raster(routes = sf_linestring_proj, dem = dem_lisbon_raster) +``` + +Likewise, if your linestring object already has X, Y and Z coordinates (e.g. from a GPS device), you can use the `slope_` functions directly. +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 the gradient of roads + +Typical use cases for the package are calculating the slopes of geographic objects representing roads or other linear features (not the imaginary line shown in the previous section). These two types of input data are represented in the code output and plot below. ```{r dem-lisbon} @@ -92,7 +177,7 @@ summary(lisbon_road_segments$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 `elevation_add()` in the [3D Analyst extension](https://desktop.arcgis.com/en/arcmap/10.3/tools/3d-analyst-toolbox/slope.htm), 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): +This result is equivalent to that returned by ESRI's `Slope_3d()` in the [3D Analyst extension](https://desktop.arcgis.com/en/arcmap/10.3/tools/3d-analyst-toolbox/slope.htm), 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): ```{r} cor( @@ -123,7 +208,7 @@ mapview::mapview(lisbon_route) ``` -We can convert the `lisbon_route` object into a 3d linestring object as follows: +We can convert the `lisbon_route` object into a 3d linestring object with X, Y and Z coordinates, as follows: ```{r, warning=FALSE, echo=FALSE, eval=FALSE} sln = stplanr::SpatialLinesNetwork(lisbon_road_segments) @@ -131,8 +216,8 @@ points = sf::st_as_sf(crs = 4326, coords = c("X1", "X2"), data.frame(rbind( stplanr::geo_code("Santa Catarina, Lisbon"), stplanr::geo_code("Castelo, Lisbon") ))) -points_projected = sf::st_transform(points, sf::st_crs(lisbon_road_segments)) -coords = sf::st_coordinates(points_projected) +points_proj = sf::st_transform(points, sf::st_crs(lisbon_road_segments)) +coords = sf::st_coordinates(points_proj) nodes = stplanr::find_network_nodes(sln, coords[, 1], coords[, 2]) lisbon_route = stplanr::sum_network_routes(sln = sln, start = nodes[1], end = nodes[2]) mapview::mapview(lisbon_route) + @@ -143,26 +228,32 @@ usethis::use_data(lisbon_route_3d, overwrite = TRUE) ``` ```{r} -lisbon_route_3d = elevation_add(lisbon_route, dem_lisbon_raster) +lisbon_route_xyz = elevation_add(lisbon_route, dem_lisbon_raster) ``` We can now visualise the elevation profile of the route as follows: ```{r plot_slope} -plot_slope(lisbon_route_3d) +plot_slope(lisbon_route_xyz) ``` -If you do not have a raster dataset representing elevations, you can automatically download them as follows. +If you do not have a raster dataset representing elevations, you can automatically download them as follows (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): -```{r, message=FALSE, warning=FALSE} -lisbon_route_3d_auto = elevation_add(lisbon_route) -plot_slope(lisbon_route_3d_auto) +```{r, message=FALSE, warning=FALSE, eval=FALSE} +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: + +```{r plotauto} +lisbon_route_xyz_auto = elevation_add(lisbon_route) +plot_slope(lisbon_route_xyz_auto) +``` ## Code of Conduct Please note that the slopes project is released with a [Contributor Code of Conduct](https://contributor-covenant.org/version/2/0/CODE_OF_CONDUCT.html). By contributing to this project, you agree to abide by its terms. - - diff --git a/README.md b/README.md index 7072cb8..90e6124 100644 --- a/README.md +++ b/README.md @@ -12,9 +12,22 @@ coverage](https://codecov.io/gh/itsleeds/slopes/branch/master/graph/badge.svg)]( Review](https://badges.ropensci.org/420_status.svg)](https://github.com/ropensci/software-review/issues/420) -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. +The **slopes** R package calculates the slope (longitudinal steepness, +also known as gradient) of linear features such as roads and rivers, +based on two main inputs: + +- [vector](https://geocompr.robinlovelace.net/spatial-class.html#vector-data) + linestring geometries defined by classes in the [`sf` + package](https://r-spatial.github.io/sf/) +- [raster](https://geocompr.robinlovelace.net/spatial-class.html#raster-data) + objects with pixel values reporting average height, commonly known + as digital elevation model (DEM) datasets, defined by classes in the + [`raster`](https://cran.r-project.org/package=raster) or more recent + [`terra`](https://rspatial.org/terra) packages + +This README covers installation and basic usage. For more information +about slopes and how to use the package to calculate them, see the +[slopes vignette](https://itsleeds.github.io/slopes/). ## Installation @@ -50,7 +63,7 @@ usethis::edit_r_environ() MAPBOX_API_KEY=xxxxx ``` -## Usage +## Basic example Load the package in the usual way: @@ -62,15 +75,127 @@ We will also load the `sf` library: ``` r library(sf) -#> Warning: package 'sf' was built under R version 4.0.5 #> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1 ``` -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. +The minimum input data requirement for using the package is an `sf` +object containing LINESTRING geometries. 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 +object `lisbon_road_segment` that is contained within the package: + +``` r +class(lisbon_road_segment) +#> [1] "sf" "tbl_df" "tbl" "data.frame" +st_geometry_type(lisbon_road_segment) +#> [1] LINESTRING +#> 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE +``` + +Don’t worry if you don’t yet have your linear features in this class: +you can read-in data from a wide range of formats into an `sf` object. +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](https://r-spatial.github.io/sf/articles/sf1.html)): + +``` r +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) +#> [1] "sf" "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: + +``` r +library(tmap) +tmap_mode("view") +#> tmap mode set to interactive viewing +tm_shape(sf_linestring) + + tm_lines(lwd = 5) + + tm_basemap(leaflet::providers$Esri.WorldTopoMap) +``` + + + +Imagine you want to calculate the gradient of the imaginary route shown +above, starting at the Castelo de São Jorge and descending towards the +coast. + +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): + +``` r +sf_linestring_xyz = elevation_add(sf_linestring) +#> Loading required namespace: ceramic +#> Preparing to download: 16 tiles at zoom = 17 from +#> https://api.mapbox.com/v4/mapbox.terrain-rgb/ +``` + +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 90 m above sea level to only 24 m in a short distance): + +``` r +st_coordinates(sf_linestring_xyz) +#> X Y Z L1 +#> [1,] -9.1333 38.714 91.4 1 +#> [2,] -9.1340 38.712 66.9 1 +#> [3,] -9.1300 38.710 24.0 1 +``` + +Step 2: calculate the average slope of the linestring + +``` r +slope_xyz(sf_linestring_xyz) +#> 1 +#> 0.104923 +``` + +The result, just over 0.1, tells us that it’s quite a steep slope, a 10% +gradient *on average*. + +If you already have a DEM, you can calculate the slopes directly as +follows, with `slope_raster()`: + +``` r +class(dem_lisbon_raster) +#> [1] "RasterLayer" +#> attr(,"package") +#> [1] "raster" +sf_linestring_proj = st_transform(sf_linestring, st_crs(lisbon_road_segment)) +slope_raster(routes = sf_linestring_proj, dem = dem_lisbon_raster) +#> 1 +#> 0.1218691 +``` + +Likewise, if your linestring object already has X, Y and Z coordinates +(e.g. from a GPS device), you can use the `slope_` functions directly. +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 the gradient of roads + +Typical use cases for the package are calculating the slopes of +geographic objects representing roads or other linear features (not the +imaginary line shown in the previous section). These two types of input +data are represented in the code output and plot below. ``` r # A raster dataset included in the package: @@ -105,7 +230,7 @@ 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 `elevation_add()` in the [3D Analyst +to that returned by ESRI’s `Slope_3d()` in the [3D Analyst extension](https://desktop.arcgis.com/en/arcmap/10.3/tools/3d-analyst-toolbox/slope.htm), with a correlation between the ArcMap implementation and our implementation of more than 0.95 on our test dataset (we find higher @@ -138,7 +263,6 @@ the Castelo de Sao Jorge to the West of the map: ``` r library(tmap) -#> Warning: package 'tmap' was built under R version 4.1.0 tmap_mode("view") #> tmap mode set to interactive viewing qtm(lisbon_route) @@ -146,32 +270,44 @@ qtm(lisbon_route) -We can convert the `lisbon_route` object into a 3d linestring object as -follows: +We can convert the `lisbon_route` object into a 3d linestring object +with X, Y and Z coordinates, as follows: ``` r -lisbon_route_3d = elevation_add(lisbon_route, dem_lisbon_raster) +lisbon_route_xyz = elevation_add(lisbon_route, dem_lisbon_raster) ``` We can now visualise the elevation profile of the route as follows: ``` r -plot_slope(lisbon_route_3d) +plot_slope(lisbon_route_xyz) ``` If you do not have a raster dataset representing elevations, you can -automatically download them as follows. +automatically download them as follows (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): + +``` r +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: ``` r -lisbon_route_3d_auto = elevation_add(lisbon_route) +lisbon_route_xyz_auto = elevation_add(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) +plot_slope(lisbon_route_xyz_auto) ``` - + ## Code of Conduct diff --git a/man/figures/README-dem-lisbon-1.png b/man/figures/README-dem-lisbon-1.png index d0d9be1..5e343bc 100644 Binary files a/man/figures/README-dem-lisbon-1.png and b/man/figures/README-dem-lisbon-1.png differ diff --git a/man/figures/README-plot_slope-1.png b/man/figures/README-plot_slope-1.png index 2bd7381..22c3dc4 100644 Binary files a/man/figures/README-plot_slope-1.png and b/man/figures/README-plot_slope-1.png differ diff --git a/man/figures/README-route-1.png b/man/figures/README-route-1.png index 0e06c9a..316f210 100644 Binary files a/man/figures/README-route-1.png and b/man/figures/README-route-1.png differ diff --git a/man/figures/README-slope-vis-1.png b/man/figures/README-slope-vis-1.png index e3c041c..ac4f788 100644 Binary files a/man/figures/README-slope-vis-1.png and b/man/figures/README-slope-vis-1.png differ diff --git a/man/figures/README-unnamed-chunk-10-1.png b/man/figures/README-unnamed-chunk-10-1.png deleted file mode 100644 index 7ee65f7..0000000 Binary files a/man/figures/README-unnamed-chunk-10-1.png and /dev/null differ diff --git a/man/figures/README-unnamed-chunk-12-1.png b/man/figures/README-unnamed-chunk-12-1.png deleted file mode 100644 index 7ee65f7..0000000 Binary files a/man/figures/README-unnamed-chunk-12-1.png and /dev/null differ diff --git a/man/figures/README-unnamed-chunk-9-1.png b/man/figures/README-unnamed-chunk-9-1.png deleted file mode 100644 index 6968c06..0000000 Binary files a/man/figures/README-unnamed-chunk-9-1.png and /dev/null differ