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:
vertex
table holding the coordinates of all OSM
nodes;edge
table mapping all edge connections between
vertices;object_link_edge
table linking all edge
entities to the OSM objects of which they are part;object
table holding all ‘key’–‘value’ pairs for
each OSM way object;relation_properties
table holding all ‘key’–‘value’
pairs for each OSM relation object;relation_members
table holding all members of each
OSM relation; andnodes
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.
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
## # A tibble: 345,239 × 3
## x_ y_ vertex_
## <dbl> <dbl> <chr>
## 1 -95.9 41.2 31536366
## 2 -95.9 41.2 31536367
## 3 -95.9 41.2 31536368
## 4 -95.9 41.2 31536370
## 5 -95.9 41.2 31536378
## 6 -95.9 41.2 31536379
## 7 -96.2 41.3 133898322
## 8 -96.2 41.3 133898328
## 9 -96.3 41.3 133898340
## 10 -96.3 41.3 133898342
## # ℹ 345,229 more rows
Incorporating elevation data is then as simple as
dat <- osm_elevation (dat, elev_file = "/path/to/elevation/data/filename.tiff")
## Loading required namespace: raster
## 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_
## # A tibble: 345,239 × 4
## x_ y_ z_ vertex_
## <dbl> <dbl> <dbl> <chr>
## 1 -95.9 41.2 291 31536366
## 2 -95.9 41.2 295 31536367
## 3 -95.9 41.2 297 31536368
## 4 -95.9 41.2 301 31536370
## 5 -95.9 41.2 295 31536378
## 6 -95.9 41.2 300 31536379
## 7 -96.2 41.3 359 133898322
## 8 -96.2 41.3 359 133898328
## 9 -96.3 41.3 358 133898340
## 10 -96.3 41.3 358 133898342
## # ℹ 345,229 more rows
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
## # A tibble: 376,370 × 8
## .vx0 .vx1 edge_ .vx0_x .vx0_y .vx1_x .vx1_y zmn
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1903265686 1903265664 V6kgqvWjtM -96.2 41.3 -96.2 41.3 351
## 2 1903265664 1903265638 mX4HQkykiD -96.2 41.3 -96.2 41.3 352
## 3 1903265638 1903265710 26e5NHT8nI -96.2 41.3 -96.2 41.3 352
## 4 1903265710 1903265636 9TOmVAvGH4 -96.2 41.3 -96.2 41.3 352
## 5 1903265636 1903265685 hYbpf832vX -96.2 41.3 -96.2 41.3 352
## 6 1903265685 1903265678 ctvd1FWGEw -96.2 41.3 -96.2 41.3 352
## 7 1903265678 1903265646 mvaAOdSOKA -96.2 41.3 -96.2 41.3 352
## 8 1903265646 1903265714 dSVFPNDFty -96.2 41.3 -96.2 41.3 352
## 9 1903265714 1903265659 uc8L3jGR87 -96.2 41.3 -96.2 41.3 352
## 10 1903265659 1903265702 MpjXnvIvcF -96.2 41.3 -96.2 41.3 352
## # ℹ 376,360 more rows
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
)
(The result is not shown here, but can be directly inspected by simply running the above lines.)