Update README for #32

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,

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)
 <!-- badges: end -->
 
-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,16 +208,16 @@ 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)
 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.
 
-
-