ohsome

Aggregating OSM elements

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 sf object franconia that is included in the mapview package and contains boundary polygons of the 37 districts of the region:

library(mapview)

q <- ohsome_elements_count(franconia, filter = "craft=brewery")

The resulting 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.

ohsome_post(q, strict = FALSE)
#> Warning: The time parameter is not defined and defaults to the latest available timestamp within the underlying OSHDB.
#> You can use set_time() to set the time parameter.
#>             timestamp value
#> 1 2023-08-10 20:00:00   180

As you can see, ohsome_post() issues a warning that the time parameter of the query is not defined. The ohsome API returns the number of elements at the latest available timestamp by default. ¹

Defining the time parameter unlocks the full power of ohsome API by giving access to the OSM history. The time 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:

ohsome_elements_count(franconia, filter = "craft=brewery", time = "2010/2020/P1M")

Alternatively, we can update the existing ohsome_query object q with the set_time() function, pipe ² the modified query directly into ohsome_post() and make a quick visualization with ggplot2:

library(ggplot2)

q |> 
    set_time("2010/2020/P1M") |>
    ohsome_post() |>
    ggplot(aes(x = timestamp, y = value)) +
    geom_line()

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:

The 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 character vector: set_endpoint(q, c("groupBy", "boundary"), append = TRUE) ³.

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

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 named id. You can use mutate() or rename() from the dplyr package to create such a column as in the example below.

By default, ohsome_post() returns an sf 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 return_value to density. This will modify the endpoint path of the query so that ohsome returns the number of breweries per area instead of the absolute value:

franconia |> 
    mutate(id = NAME_ASCI) |>
    ohsome_elements_count(filter = "craft=brewery", return_value = "density") |>
    set_time("2021-06-01") |>
    set_grouping("boundary") |>
    ohsome_post() |>
    mapview(zcol = "value", layer.name = "Breweries per sqkm")

Extracting OSM elements

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 ohsome_elements_centroid() to 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.

Extracting the full history of OSM elements

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 validFrom and validTo timestamps.

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

You may find using clean_names() from the janitor package helpful in order to remove special characters from column names in the parsed ohsome API response – just as in the example above.

Aggregating OSM contributions

With queries to the ohsome API’s contributions aggregation endpoints, you can get counts of the contributions provided by users to OSM. The following code requests the number of deletions of man-made objects at the location of the hypothetical Null Island per year between 2010 and 2020:

ohsome_contributions_count(
    boundary = "0,0,10", 
    filter = "man_made=*", 
    time = "2010/2020/P1Y",
    contributionType = "deletion"
) |>
    ohsome_post()
#>    fromTimestamp toTimestamp value
#> 1     2010-01-01  2011-01-01     0
#> 2     2011-01-01  2012-01-01     1
#> 3     2012-01-01  2013-01-01     0
#> 4     2013-01-01  2014-01-01     0
#> 5     2014-01-01  2015-01-01     0
#> 6     2015-01-01  2016-01-01     3
#> 7     2016-01-01  2017-01-01     1
#> 8     2017-01-01  2018-01-01     8
#> 9     2018-01-01  2019-01-01     4
#> 10    2019-01-01  2020-01-01     4

The contributionType parameter is used to filter for specific types of contributions (in this case: deletions). If it is not set, any contribution is counted. Note that the resulting values apply to time intervals defined by a fromTimestamp and a toTimestamp.

Extracting OSM contributions

The contributions extraction endpoints of the ohsome API can be used to extract feature geometries of contributions.

In the following example, we extract the centroids of all amenities in the Berlin city district of Neukölln that have had contributions in March 2020. Consequently, we filter for features that have had tags changed and visualize their locations:

nominatimlite::geo_lite_sf("Berlin Neukoelln", points_only = FALSE) |>
    ohsome_contributions_centroid() |>
    set_filter("amenity=*") |>
    set_time("2020-03,2020-04") |>
    set_properties("contributionTypes") |> 
    ohsome_post() |>
    filter(`@tagChange`) |>
    mapview(layer.name = "Amenities with Tag Changes")
#> Warning: 60 element(s) with empty geometries omitted.