Querying single element groups from OpenStreetMap using Overpass and osmdata in R

Question

I'm trying to extract administrative boundary data from OpenStreetMap using osmdata in R. The Overpass API lets you determine the element group (i.e. nodes, ways, relations) to return. For example, in the following code, only relations are returned:

area(3600062761)->.searchArea;
(
    relation["admin_level"="6"]["boundary"="administrative"](area.searchArea);
); 
out geom;

However, I couldn't replicate this with osmdata in R. The following R code:

osm_bounds <- getbb(place_name = 'NRW', format_out = 'polygon')[[1]] %>%
    opq(timeout = 1000) %>%
    add_osm_feature(key = 'admin_level', value = '6') %>%
    add_osm_feature(key = 'boundary', value = 'administrative') %>%
    osmdata_sf()

... generates an Overpass query that returns nodes, ways and relations:

[out:xml][timeout:1000];
    (
        node  ["admin_level"="6"] ["boundary"="administrative"] (50.3226897,5.8663153,52.5314923,9.4617417);
        way  ["admin_level"="6"] ["boundary"="administrative"] (50.3226897,5.8663153,52.5314923,9.4617417);
        relation  ["admin_level"="6"] ["boundary"="administrative"] (50.3226897,5.8663153,52.5314923,9.4617417);
    );
(._;>;);
out body;

Thus, I have to download all elements and select the ones of interest afterwards. The first Overpass query downloads approximately 120 MB while the second query downloads 700 MB and takes an eternity.

I also tried formatting the Overpass code as a string and then passing it to osmdata::osmdata_sf():

osm_id <- tmaptools::geocode_OSM('NRW', details = TRUE)$osm_id
overpass_string <- paste0('area(%s)->.searchArea;',
                          '(relation["admin_level"="%s"]["boundary"="administrative"](area.searchArea);',
                          '); out geom;') %>%
  sprintf(3600000000 + as.integer(osm_id), 6)
osmdata_sf(overpass_string)

However, this returns an empty osmdata object or an error (seemingly randomly one of these two outcomes):

Object of class 'osmdata' with:
                 $bbox : 
        $overpass_call : The call submitted to the overpass API
                 $meta : metadata including timestamp and version numbers
           $osm_points : 'sf' Simple Features Collection with 0 points
            $osm_lines : NULL
         $osm_polygons : 'sf' Simple Features Collection with 0 polygons
       $osm_multilines : NULL
    $osm_multipolygons : NULL

Error in rcpp_osmdata_sf(doc): way can not be found

My question is if there are other ways of determining the element groups that are to be selected using osmdata or if there are workarounds to accomplish this.

Answer 1

I figured that there seems to be no alternative to osmdata which is why I decided to craft my own solution.

First, I formatted the query similar to the code I posted in my question:

query.osm.boundaries <- function(admin_level, region) {
  osm_id <- tmaptools::geocode_OSM(region, details = TRUE)$osm_id
  query <- paste0(
    '[out:json][timeout:100];',
    'area(%s)->.searchArea;',
    '(relation["admin_level"="%s"]["boundary"="administrative"](area.searchArea););',
    'out geom;'
  ) %>%
    sprintf(3600000000 + as.integer(osm_id), admin_level)
  response <- query.overpass(query)
  osm_bounds <- relations.to.sf(response)
  return(osm_bounds)
}

The query is then passed to a function that queries Overpass and returns a parsed json using the packages httr and jsonlite:

query.overpass <- function(query) {
  # Get Overpass-URL from osmdata and parse it
  url <- osmdata::get_overpass_url()
  parsed_url <- httr::parse_url(url)
  # Create a url from the query
  parsed_url$query <- list(data = query)
  request <- httr::build_url(parsed_url)
  # Send GET request and parse the response json
  response <- httr::GET(request) %>%
    httr::content(as = 'text', encoding = 'UTF-8', type = 'application/json') %>%
    jsonlite::fromJSON()
  return(response)
}

And finally, the coordinates are extracted from the parsed JSON and converted to sf data frames. This is done by extracting coordinates iteratively, turning them to multipoints, linestrings, multilinestrings and finally to polygons/multipolygons.

relations.to.sf <- function(osm_json) {
  get.multilines <- function(way) {
    way %>%
      # Convert coordinates to a matrix,  
      as.matrix() %>%
      # (small custom function to swap lon and lat)
      swap.xy() %>%
      # ... then to multipoints, 
      sf::st_multipoint() %>%
      # ... and then to a linestring
      sf::st_linestring()
  }
  get.multipolygons <- function(relations_index) {
    # Finally, put linestrings together to form a multilinestring
    multilines <- map_if(relations$members[[relations_index]]$geometry, ~!is.null(..1), get.multilines) %>%
      dplyr::discard(sapply(., is.null)) %>%
      sf::st_multilinestring()
    # Extract column names from the json
    cols <- relations$tags[relations_index, ]
    # Convert multilinestrings to a single boundary polygon
    polygon <- sf::st_polygonize(multilines) %>%
      sf::st_collection_extract() %>%
      # Pass geometries and tags to sf
      sf::st_sfc() %>%
      sf::st_sf() %>%
      cbind(cols)
    if (nrow(polygon) > 1) {
      # If a boundary consists of more than one polygon, glue them together to form a multipolygon
      polygon <- polygon %>%
        dyplr::summarise(geometry = sf::st_union(geometry)) %>%
        cbind(cols)
    }
    return(polygon)
  }
  relations <- osm_json$elements
  output_polygons <- map(seq_len(length(relations$members)), get.multipolygons) %>%
    # Summarise single boundary polygons to a big sf dataset
    do.call(rbind, .) %>%
    # ... and set its CRS
    sf::st_set_crs(4326)
  return(output_polygons)
}

Using these functions, I managed to significantly improve the processing time compared to downloading everything and then selecting the elements of interest:

start <- proc.time()
query.osm.boundaries(6, 'NRW') # new function
proc.time() - start
#      user      system     elapsed 
#      2.64        0.00       18.72 

start <- proc.time()
query.osm.boundaries2(6, 'NRW') # old function
proc.time() - start
#      user      system     elapsed 
#     12.00        0.86       81.18 

This whole process seems unnecessarily cumbersome. It seems that, coming from Python, which offers a multitude of different packages to query Overpass and convert Overpass JSONs to GeoJSON, I was falsely expecting R to have a similar flexibility. If there is anything I overlooked, I would appreciate a heads up!