This ohsome R package grants access to the power of the ohsome API from R. ohsome lets you analyze the rich data source of the OpenStreetMap (OSM) history. It aims to leverage the tools of the OpenStreetMap History Database (OSHDB).
With ohsome, you can …
Get aggregated statistics on the evolution of OpenStreetMap elements and specify your own temporal, spatial and/or thematic filters. The data aggregation endpoint allows you to access functions, e.g., to calculate the area of buildings or the length of streets at any given timestamp.
Retrieve the geometry of the historical OpenStreetMap data, e.g., to visualize the evolution of certain OpenStreetMap elements over time. You can get the geometries for specific points in time or all changes within a timespan (full-history).
Upon attaching the ohsome package, a metadata request is sent to the ohsome API. The package message provides some essential metadata information, such as the current temporal extent of the underlying OSHDB:
The metadata is stored in
.ohsome_metadata. You can
print it to the console to get more details.
You can create any ohsome API query using the generic
ohsome_query() function. It takes the endpoint path and any
query parameters as inputs. For information on all available endpoints
with their parameters, consult the ohsome API documentation or print
ohsome_endpoints to the console.
However, this ohsome R package provides specific wrapper functions for queries to all available endpoints.
The elements aggregation endpoints of the ohsome API allow querying for the aggregated amount, length, area or perimeter of OpenStreetMap elements with given properties, within given boundaries and at given points in time.
Let us create a query for the total amount of breweries on OSM in the
region of Franconia. The first argument to
ohsome_elements_count() is the
franconia that is included in the mapview
package and contains boundary polygons of the 37 districts of the
ohsome_query object can be sent to the
ohsome API with
ohsome_post(). By default,
ohsome_post() returns the parsed API response. In this
case, this is a simple
data.frame of only one row.
As you can see,
ohsome_post() issues a warning that the
time parameter of the query is not defined. The
returns the number of elements at the latest available timestamp by
time parameter unlocks the full power of
ohsome API by giving access to the OSM history. The
parameter requires one or more ISO-8601 conform timestring(s). Here is how to create
a query for the number of breweries at the first of each month between
2010 and 2020:
Alternatively, we can update the existing
q with the
set_time() function, pipe 2 the
modified query directly into
ohsome_post() and make a quick
This is how to query the total number of breweries in all of Franconia. But what if we want to aggregate the amount per district? The ohsome API provides specific endpoints for different grouped calculations, such as aggregation grouped by bounding geometry.
There are several ways to define a query for an aggregation grouped by boundary:
set_endpointfunction is used to change or append the
endpoint path of an API request. In this case, we could append
groupBy/boundary to the existing query to the
elements/count endpoint. The endpoint path can either be
given as a single string (
/groupBy/boundary) or as a
set_endpoint(q, c("groupBy", "boundary"), append = TRUE) 3.
More comfortable, however, is the use of either the grouping argument
with an elements aggregation function (e.g.
ohsome_elements_count(grouping = "boundary)) or of the
set_grouping() function to modify an existing query
library(dplyr) franconia |> mutate(id = NAME_ASCI) |> ohsome_elements_count(filter = "craft=brewery", time = "2021-06-01") |> set_grouping("boundary") |> ohsome_post() #> Simple feature collection with 37 features and 3 fields #> Geometry type: MULTIPOLYGON #> Dimension: XY #> Bounding box: xmin: 8.975926 ymin: 48.8625 xmax: 12.27535 ymax: 50.56422 #> Geodetic CRS: WGS 84 #> First 10 features: #> groupByBoundaryId value timestamp geometry #> 1 Bamberg, Kreisfreie Stadt 6 2021-06-01 MULTIPOLYGON (((10.92581 49... #> 2 Bayreuth, Kreisfreie Stadt 6 2021-06-01 MULTIPOLYGON (((11.58157 49... #> 3 Coburg, Kreisfreie Stadt 0 2021-06-01 MULTIPOLYGON (((10.95355 50... #> 4 Hof, Kreisfreie Stadt 1 2021-06-01 MULTIPOLYGON (((11.93067 50... #> 5 Bamberg, Landkreis 13 2021-06-01 MULTIPOLYGON (((10.87615 50... #> 6 Bayreuth, Landkreis 13 2021-06-01 MULTIPOLYGON (((11.70656 50... #> 7 Coburg, Landkreis 6 2021-06-01 MULTIPOLYGON (((10.88654 50... #> 8 Forchheim 8 2021-06-01 MULTIPOLYGON (((11.26376 49... #> 9 Hof, Landkreis 4 2021-06-01 MULTIPOLYGON (((11.91989 50... #> 10 Kronach 1 2021-06-01 MULTIPOLYGON (((11.36979 50...
If you want your own identifiers for the geometries returned by
ohsome, your input
sf object needs a column explicitly
id. You can use
rename() from the dplyr package to create such a column as in the
ohsome_post() returns an
object whenever the ohsome API is capable of delivering GeoJSON data.
This is the case for elements extraction queries as well as for
aggregations grouped by boundaries.
Thus, you can easily create a choropleth map from the query results.
In addition, you can set the argument
density. This will modify the endpoint path of the query so
that ohsome returns the number of breweries per area instead of the
The elements extraction endpoints of the ohsome API allow obtaining geometries, bounding boxes or centroids of OSM elements with given properties, within given boundaries and at given points in time. Together with the elements, you can choose to query for their tags and/or their metadata such as the changeset ID, the time of the last edit or the version number.
The following query extracts the geometries of buildings within 500 m
of Heidelberg main station with their tags. The response is used to
visualize the buildings and the values of their
building:levels tag (if available):
hd_station_500m <- ohsome_boundary("8.67542,49.40347,500") ohsome_elements_geometry( boundary = hd_station_500m, filter = "building=* and type:way", time = "2021-12-01", properties = "tags", clipGeometry = FALSE ) |> ohsome_post() |> transmute(level = factor(`building:levels`)) |> mapview(zcol = "level", lwd = 0, layer.name = "Building level")
Similarly, you can use
extract centroids of OSM elements and
ohsome_elements_bbox() for their bounding boxes. Note that
OSM node elements (with point geometries) are omitted from the results
if querying for bounding boxes.
While the elements extraction endpoints provide geometries and
properties of OSM elements at specific timestamps, results of queries to
the full history endpoints will include all changes to
matching OSM features with corresponding
Here, we request the full history of OSM buildings within 500 m of Heidelberg main station, filter for features that still exist and visualize all building features with their year of creation:
hd_station_1km <- ohsome_boundary("8.67542,49.40347,1000") ohsome_elements_geometry( boundary = hd_station_1km, filter = "building=* and type:way", time = "2021-12-01", properties = "tags", clipGeometry = FALSE ) |> ohsome_post() |> transmute(level = factor(`building:levels`)) |> mapview(zcol = "level", lwd = 0, layer.name = "Building level")