---
title: "2. Elevation data and OSM: The osmdata_sc function"
author: 
  - "Mark Padgham"
date: "`r Sys.Date()`"
output: 
    html_document:
        toc: true
        toc_float: true
        number_sections: false
        theme: flatly
vignette: >
  %\VignetteIndexEntry{2. osmdata_sc}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

## 1. Introduction


[`silicate`](https://github.com/hypertidy/silicate) is a new
form for representing spatial data in **R**. In contrast to all
other forms (such as [`sp`](https://cran.r-project.org/package=sp) or
[`sf`](https://cran.r-project.org/package=sf)),
[`silicate`](https://github.com/hypertidy/silicate) is *multi-tabular*, and
primarily consists of one table for point entities; one table for binary
relationships between these point entities -- spatial ''edges'' -- and
additional tables for higher-order inter-relationships between these objects.
The new `osmdata` function `osmdata_sc()` returns Open Street Map (OSM) data in
[`silicate`](https://github.com/hypertidy/silicate) form. This form also closely
resembles the data storage scheme of Open Street Map itself, and in this case
consists of the following tables:

1. A `vertex` table holding the coordinates of all OSM nodes;
2. An `edge` table mapping all edge connections between vertices;
3. An `object_link_edge` table linking all `edge` entities to the OSM objects of
   which they are part;
4. An `object` table holding all 'key'--'value' pairs for each OSM *way* object;
5. A `relation_properties` table holding all 'key'--'value' pairs for each OSM
   *relation* object;
6. A `relation_members` table holding all members of each OSM relation; and
7. A `nodes` table holding all 'key'--'value' pairs for each OSM *node* object.

The translation of the underlying OSM data structure -- consisting of nodes,
way, and relations -- into 
[Simple Features (SF)](https://www.ogc.org/standards/sfa/) via the
[`osmdata_sf()`](https://docs.ropensci.org/osmdata/reference/osmdata_sf.html)
function is less than 100% faithful, and results in some representational loss
compared with the original OSM structure (for details, see the [vignette on
translation of OSM into
SF](https://docs.ropensci.org/osmdata/articles/osm-sf-translation.html)). In
contrast, the [`osmdata_sc()`
function](https://docs.ropensci.org/osmdata/reference/osmdata_sc.html) delivers
a representation that is entirely faithful to the underlying OSM representation.

One of the advantages of [`silicate`](https://github.com/hypertidy/silicate)
format offered by the `osmdata` package is enabling elevation data to be
combined with OSM data. The result is
a [`silicate`](https://github.com/hypertidy/silicate)-format object which is
able to be submitted directly to the [`dodgr`
package](https://github.com/UrbanAnalyst/dodgr) to enable routing on street
networks that accounts for elevation changes.


## 2. Elevation Data

Incorporating elevation data with OSM data currently requires local storage of
desired elevation data. These must be downloaded for the desired region from 
[http://srtm.csi.cgiar.org/srtmdata](https://srtm.csi.cgiar.org/srtmdata/) in Geo
TIFF format. Elevation data may then be incorporated with 
[`silicate`](https://github.com/hypertidy/silicate)-format data generated by 
[`x <- osmdata_sc()`](https://docs.ropensci.org/osmdata/reference/osmdata_sc.html) 
through the `osm_elevation()` function. The entire procedure is demonstrated
with the following lines:
```{r omaha, eval = FALSE}
dat <- opq ("omaha nebraska") |>
    add_osm_feature (key = "highway") |>
    osmdata_sc ()
```
This object has a `vertex` table like this:
```{r dat_vertex, eval = FALSE}
dat$vertex
```
```{r dat_vertex_dat, echo = FALSE}
n <- 345239
x_ <- c (-95.9, -95.9, -95.9, -95.9, -95.9, -95.9, -96.2, -96.2, -96.3, -96.3)
y_ <- c (41.2, 41.2, 41.2, 41.2, 41.2, 41.2, 41.3, 41.3, 41.3, 41.3)
z_ <- c (291.0, 295.0, 297.0, 301.0, 295.0, 300.0, 359.0, 359.0, 358.0, 358.0)
vertex_ <- paste0 (c (
    31536366, 31536367, 31536368, 31536370, 31536378, 31536379,
    133898322, 133898328, 133898340, 133898342
))

tibble::tibble (
    x_ = c (x_, rep (NA, n - 10)),
    y_ = c (y_, rep (NA, n - 10)),
    vertex_ = c (vertex_, rep (NA, n - 10))
)
```

Incorporating elevation data is then as simple as
```{r osm_elevation, eval = FALSE}
dat <- osm_elevation (dat, elev_file = "/path/to/elevation/data/filename.tiff")
```
```{r osm_elevation2, echo = FALSE}
message (
    "Loading required namespace: raster\n",
    "Elevation data from Consortium for Spatial Information; see ",
    "https://srtm.csi.cgiar.org/srtmdata/"
)
```
This function then simply appends the elevation values to the `vertex_` table,
so that it now looks like this:
```{r dat_vertex2, eval = FALSE}
dat$vertex_
```
```{r, dat_vertex_dat2, echo = FALSE}
tibble::tibble (
    x_ = c (x_, rep (NA, n - 10)),
    y_ = c (y_, rep (NA, n - 10)),
    z_ = c (z_, rep (NA, n - 10)),
    vertex_ = c (vertex_, rep (NA, n - 10))
)
```

### Example usage of elevation data

The [`silicate`](https://github.com/hypertidy/silicate) format is very easy to
manipulate using standard [`dplyr`](https://dplyr.tidyverse.org) verbs. The
following code uses the [`mapdeck`
package](https://github.com/SymbolixAU/mapdeck) package to colour the street
network and elevation data downloaded and processed in the preceding lines by
the elevation of each network edge. We first join the vertex elevation data on
to the edges, and calculate the mean elevation of each edge.

```{r edges, eval = FALSE}
edges <- dplyr::left_join (dat$edge, dat$vertex, by = c (".vx0" = "vertex_")) |>
    dplyr::rename (".vx0_x" = x_, ".vx0_y" = y_, ".vx0_z" = z_) |>
    dplyr::left_join (dat$vertex, by = c (".vx1" = "vertex_")) |>
    dplyr::rename (".vx1_x" = x_, ".vx1_y" = y_, ".vx1_z" = z_) |>
    dplyr::mutate ("zmn" = (.vx0_z + .vx1_z) / 2) |>
    dplyr::select (-c (.vx0_z, .vx1_z))
edges
```
```{r edges-dat, echo = FALSE}
n <- 376370
x <- paste0 (c (
    1903265686, 1903265664, 1903265638, 1903265710, 1903265636,
    1903265685, 1903265678, 1903265646, 1903265714, 1903265659
))
y <- paste0 (c (
    1903265664, 1903265638, 1903265710, 1903265636, 1903265685,
    1903265678, 1903265646, 1903265714, 1903265659, 1903265702
))
edge <- c (
    "V6kgqvWjtM", "mX4HQkykiD", "26e5NHT8nI", "9TOmVAvGH4", "hYbpf832vX",
    "ctvd1FWGEw", "mvaAOdSOKA", "dSVFPNDFty", "uc8L3jGR87", "MpjXnvIvcF"
)
x0_x <- c (-96.2, -96.2, -96.2, -96.2, -96.2, -96.2, -96.2, -96.2, -96.2, -96.2)
x0_y <- c (41.3, 41.3, 41.3, 41.3, 41.3, 41.3, 41.3, 41.3, 41.3, 41.3)
x1_x <- c (-96.2, -96.2, -96.2, -96.2, -96.2, -96.2, -96.2, -96.2, -96.2, -96.2)
x1_y <- c (41.3, 41.3, 41.3, 41.3, 41.3, 41.3, 41.3, 41.3, 41.3, 41.3)
z <- c (351.0, 352.0, 352.0, 352.0, 352.0, 352.0, 352.0, 352.0, 352.0, 352.0)

tibble::tibble (
    ".vx0" = c (x, rep (NA, n - 10)),
    ".vx1" = c (y, rep (NA, n - 10)),
    "edge_" = c (edge, rep (NA, n - 10)),
    ".vx0_x" = c (x0_x, rep (NA, n - 10)),
    ".vx0_y" = c (x0_y, rep (NA, n - 10)),
    ".vx1_x" = c (x1_x, rep (NA, n - 10)),
    ".vx1_y" = c (x1_y, rep (NA, n - 10)),
    "zmn" = c (z, rep (NA, n - 10))
)
```

Those data can then be submitted directly to
[`mapdeck`](https://github.com/SymbolixAU/mapdeck) to generate an interactive
plot with the following code:
```{r, eval = FALSE}
library (mapdeck)
set_token (Sys.getenv ("MAPBOX_TOKEN")) # load local token for MapBox
mapdeck (style = mapdeck_style ("dark")) |>
    add_line (edges,
        origin = c (".vx0_x", ".vx0_y"),
        destination = c (".vx1_x", ".vx1_y"),
        stroke_colour = "z",
        legend = TRUE
    )
```
(The result is not shown here, but can be directly inspected by simply running
the above lines.)