---
title: "Get started"
output: bookdown::html_vignette2
vignette: >
  %\VignetteIndexEntry{Get started}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

Welcome to the slopes vignette, a type of long-form documentation/article that introduces the core functions and functionality of the `slopes` package.

# Installation

You can install the released version of slopes from [CRAN](https://CRAN.R-project.org) with:

```r
install.packages("slopes")
```

Install the development version from [GitHub](https://github.com/) with:

```{r, eval=FALSE}
# install.packages("remotes")
remotes::install_github("ropensci/slopes")
```

### 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:

```{r, eval=FALSE}
# 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:

```{r, eval=FALSE}
usethis::edit_r_environ()
MAPBOX_API_KEY=xxxxx # replace XXX with your api key
```


# Functions 

## 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

# Package datasets

The `slopes` package comes with some datasets to play with:

**Linestrings:**

- `lisbon_road_segment`: a single road segment in Lisbon (XY)
- `lisbon_route`: a route with some variation in elevation in Lisbon (XY)
- `cyclestreets_route`: a bike route in Leeds (XY)

**Road network:**

- `lisbon_road_network`: a sample of road segments in downtown Lisbon
- `magnolia_xy`: a sample of road segments in center Seattle, in the Magnolia neighborhood

**Digital elevation model (DEM):**

- `dem_lisbon_raster` a DEM of downtown Lisbon (EPSG:3763)


# Examples

Load the package in the usual way. We will also load the `sf` library:

```{r message=FALSE, warning=FALSE}
library(slopes)
library(sf)
```

The minimum input data requirement for using the package is an `sf` object containing LINESTRING geometries.  

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)
```

> maybe remove this? or add step 1 and step 2 again.

## 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:

```{r}
sf_linestring = lisbon_road_segment
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, message=FALSE, warning=FALSE}
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](https://account.mapbox.com/access-tokens/) key):

```{r, eval=FALSE}
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/
```

```{r, echo=FALSE}
# note: the following should be TRUE
# identical(sf_linestring_xyz, lisbon_road_segment_xyz_mapbox)
sf_linestring_xyz = lisbon_road_segment_xyz_mapbox
```

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).

```{r}
st_coordinates(sf_linestring_xyz)
```

You can use the `z_` functions to extract such values:

```{r}
z_value(sf_linestring_xyz) # returns all the elevation values between xy coordinates

z_mean(sf_linestring_xyz) # elevation mean value
z_min(sf_linestring_xyz) # elevation min value 
z_max(sf_linestring_xyz) # elevation max value 
z_start(sf_linestring_xyz) # first z
z_end(sf_linestring_xyz) # last z
```

**Step 2**: calculate the average slope of the linestring

```{r}
slope_xyz(sf_linestring_xyz)
```

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](https://ropensci.github.io/intro-to-slopes/) vignette).
  
Again, let us use the same function with a entire route, `lisbon_route`, also available in the package:

```{r message=FALSE, warning=FALSE}
sf_route = lisbon_route
class(sf_route)
st_geometry_type(sf_route)

tm_shape(sf_route) +
  tm_lines(lwd = 3) +
  tm_basemap(leaflet::providers$Esri.WorldTopoMap)
```

**Step 1**: add elevations to each coordinate in the route:

```{r, eval=FALSE}
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/
```

```{r, echo=FALSE}
# note: the following should be TRUE
# identical(sf_route_xyz, lisbon_road_segment_xyz_mapbox)
sf_route_xyz = lisbon_route_xyz_mapbox
```

**Step 2**: calculate the average slope of the route

```{r}
slope_xyz(sf_route_xyz)
```

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()`:

```{r}
class(dem_lisbon_raster)
slope_raster(routes = sf_route,
             dem = dem_lisbon_raster)
```

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](https://ropensci.github.io/slopes/articles/verification.html).) 

## 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. 

```{r eval=FALSE, include=FALSE}
#not to use like this... it would ge good to have a gps example to demonstrate

slope_vector(sf_route_xyz)
slope_distance(sf_route_xyz)
slope_distance_mean(sf_route_xyz)
slope_distance_weighted(sf_route_xyz)

slope_vector(sf_linestring_xyz)
slope_distance(sf_linestring_xyz)
slope_distance_mean(sf_linestring_xyz)
slope_distance_weighted(sf_linestring_xyz)
```

```{r}
# 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)

# 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)
slope_distance_mean(dist, elevations)
slope_distance_weighted(dist, elevations)
```

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 and plotting gradients

## 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.

```{r dem-lisbon}
# A raster dataset included in the package:
class(dem_lisbon_raster) # digital elevation model
summary(raster::values(dem_lisbon_raster)) # heights range from 0 to ~100m
raster::plot(dem_lisbon_raster)

# A vector dataset included in the package:
class(lisbon_road_network)
plot(sf::st_geometry(lisbon_road_network), add = TRUE)
```

Calculate the average gradient of **each road segment** as follows:

```{r}
lisbon_road_network$slope = slope_raster(lisbon_road_network, dem = dem_lisbon_raster)
summary(lisbon_road_network$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 `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 - see the [verification of slopes](https://ropensci.github.io/slopes/articles/verification.html article):

```{r}
cor(
  lisbon_road_network$slope,    # slopes calculates by the slopes package
  lisbon_road_network$Avg_Slope # slopes calculated by ArcMap's 3D Analyst extension
)
```

We can now visualise the average slopes of each route calculated by the `slopes` package as follows:

```{r slope-vis}
raster::plot(dem_lisbon_raster)
plot(lisbon_road_network["slope"], add = TRUE, lwd = 5)
```

## Elevation profile

Taking the [first route example](#route--available-dem), 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:

```{r route}
# library(tmap)
# tmap_mode("view")
qtm(lisbon_route)
```

```{r, echo=FALSE, eval=FALSE}
# Removed because it's not rendering in RMarkdown
mapview::mapview(lisbon_road_network["slope"], map.types = "Esri.WorldStreetMap")
mapview::mapview(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:

```{r, eval=TRUE}
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_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](https://github.com/paleolimbot/qgisprocess/issues/26), but here we'll use the `stplanr` function `rnet_breakup_vertices()` (see  [`vignette("roadnetworkcycling")`](https://ropensci.github.io/slopes/articles/roadnetworkcycling.html) for an example of this function working on a large road network):

```{r}
sf::st_length(lisbon_route_xyz) # check route length: 2.5 km
lisbon_route_segments = stplanr::rnet_breakup_vertices(lisbon_route_xyz)
summary(sf::st_length(lisbon_route_segments)) # mean of 50 m
```

We can now calculate the slope for each of these segments.

```{r}
lisbon_route_segments$slope = slope_xyz(lisbon_route_segments)
summary(lisbon_route_segments$slope)
```

## 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.

```{r}
lisbon_route_segments$slope_directed = slope_xyz(lisbon_route_segments, directed = TRUE)
summary(lisbon_route_segments$slope_directed)
```

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.

```{r, fig.show='hold', out.width="50%"}
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)
```

```{r, eval=FALSE, echo=FALSE}
# test code
z = sf::st_make_grid(lisbon_route_xyz, cellsize = 100)
sampled_points = sf::st_line_sample(lisbon_route_xyz, n = 30)
points_sf = sf::st_sf(geometry = sf::st_cast(sampled_points, "POINT"))
plot(points_sf)
lisbon_route_segments = stplanr::route_split(lisbon_route_xyz, p = points_sf[2:3, ])
lisbon_route_segments = stplanr::rnet_breakup_vertices(lisbon_route_xyz)
library(tmap)
tmap_mode("view")
qtm(lisbon_route_segments, lines.lwd = 9, lines.col = 1:nrow(lisbon_route_segments))
plot(lisbon_route_segments, col = 1:nrow(lisbon_route_segments))
```


```{r, eval=FALSE, echo=FALSE}
# Test: try using QGIS
remotes::install_github("paleolimbot/qgisprocess")
library(qgisprocess)
qgis_configure()
algorithms = qgis_algorithms()
View(algorithms)
result = qgis_run_algorithm(
  algorithm = "grass7:v.split",
  INPUT = lisbon_route_xyz,
  LENGTH = 500
  )
route_segments = sf::st_read(result$OUTPUT)
route_segments
plot(lisbon_route_xyz$geometry)
plot(route_segments$geom, add = T, lwd = 3)
mapview::mapview(route_segments)
```


## 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):

```{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, eval=FALSE}
lisbon_route_xyz_auto = elevation_add(lisbon_route) #dem = NULL
```
```{r, echo=FALSE}
lisbon_route_xyz_auto = lisbon_route_xyz_mapbox
```
```{r}
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.

# Commulative elevation change

The following example calculate the elevations of a route in Leeds, and plots its commutative sum along the route (not evaluated).

```{r, eval=FALSE}
cyclestreets_xyz = elevation_add(cyclestreets_route) 
plot_slope(cyclestreets_xyz)
plot(cumsum(cyclestreets_xyz$distances), cumsum(cyclestreets_xyz$elevation_change))
```


# See more in vignettes

-   [slopes package](https://ropensci.github.io/slopes/index.html)
-   [An introduction to slopes](https://ropensci.github.io/slopes/articles/intro-to-slopes.html)
-   [Reproducible example: gradients of a road network for a given city](https://ropensci.github.io/slopes/articles/roadnetworkcycling.html)
-   [Verification of slopes](https://ropensci.github.io/slopes/articles/verification.html)
-   [Benchmarking slopes calculation](https://ropensci.github.io/slopes/articles/benchmark.html)