2. Elevation data and OSM: The osmdata_sc function
In osmdata: Import 'OpenStreetMap' Data as Simple Features or Spatial Objects

1. Introduction

silicate is a new form for representing spatial data in R. In contrast to all other forms (such as sp or sf), 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 form. This form also closely resembles the data storage scheme of Open Street Map itself, and in this case consists of the following tables:

A vertex table holding the coordinates of all OSM nodes;
An edge table mapping all edge connections between vertices;
An object_link_edge table linking all edge entities to the OSM objects of which they are part;
An object table holding all 'key'--'value' pairs for each OSM way object;
A relation_properties table holding all 'key'--'value' pairs for each OSM relation object;
A relation_members table holding all members of each OSM relation; and
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) via the osmdata_sf() 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). In contrast, the osmdata_sc() function delivers a representation that is entirely faithful to the underlying OSM representation.

One of the advantages of silicate format offered by the osmdata package is enabling elevation data to be combined with OSM data. The result is a silicate-format object which is able to be submitted directly to the dodgr package 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 in Geo TIFF format. Elevation data may then be incorporated with silicate-format data generated by x <- osmdata_sc() through the osm_elevation() function. The entire procedure is demonstrated with the following lines:

dat <- opq ("omaha nebraska") %>%
    add_osm_feature (key = "highway") %>%
    osmdata_sc ()

This object has a vertex table like this:

dat$vertex

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

dat <- osm_elevation (dat, elev_file = "/path/to/elevation/data/filename.tiff")

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:

dat$vertex_

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 format is very easy to manipulate using standard dplyr verbs. The following code uses the mapdeck package 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.

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

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 to generate an interactive plot with the following code:

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
    )