[FEATURE] Better documentation and maybe support for relations

[Robinlovelace]

It would be nice to get all features associated with a particular relation, e.g. network=ncn.

Reprex below may illustrate this better, one to talk about and test things out @agila5 but, for now, do you see broadly what I mean?

# Aim: get UK's national cycle network

library(osmextract)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright.
#> Check the package website, https://docs.ropensci.org/osmextract/, for more details.

# Test on a small amount of data
query = "SELECT * FROM 'lines' WHERE ncn = 'yes'"
et = c("oneway", "maxspeed", "ref", "ncn")
cycleways_wyca = oe_get(place = "West Yorkshire", query = query, extra_tags = et)
#> The input place was matched with: West Yorkshire
#> The chosen file was already detected in the download directory. Skip downloading.
#> The corresponding gpkg file was already detected. Skip vectortranslate operations.
#> Reading query `SELECT * FROM 'lines' WHERE ncn = 'yes'' from data source `/mnt/57982e2a-2874-4246-a6fe-115c199bc6bd/data/osm/geofabrik_west-yorkshire-latest.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 24 features and 14 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -2.126027 ymin: 53.66803 xmax: -1.345368 ymax: 53.87748
#> Geodetic CRS:  WGS 84
nrow(cycleways_wyca) # 24 rows
#> [1] 24
mapview::mapview(cycleways_wyca)

query = "SELECT * FROM 'multilinestrings' WHERE ncn is not null OR network = 'ncn'"
et = c("oneway", "maxspeed", "ref", "ncn", "network")
cycleways_wyca_relations = oe_get(place = "West Yorkshire", query = query, extra_tags = et)
#> The input place was matched with: West Yorkshire
#> Warning: The query selected a layer which is different from layer argument. We
#> will replace the layer argument.
#> The chosen file was already detected in the download directory. Skip downloading.
#> The corresponding gpkg file was already detected. Skip vectortranslate operations.
#> Reading query `SELECT * FROM 'multilinestrings' WHERE ncn is not null OR network = 'ncn'' from data source `/mnt/57982e2a-2874-4246-a6fe-115c199bc6bd/data/osm/geofabrik_west-yorkshire-latest.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 15 features and 9 fields
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -2.156127 ymin: 53.5308 xmax: -1.313181 ymax: 53.95098
#> Geodetic CRS:  WGS 84
nrow(cycleways_wyca_relations) # 27 rows
#> [1] 15
mapview::mapview(cycleways_wyca_relations)

# How to get the results above as individual lines?
# wyca_relations = oe_get(place = "West Yorkshire", ...)

^{Created on 2022-08-19 by the reprex package (v2.0.1)}

[Robinlovelace]

More specific example: https://www.openstreetmap.org/relation/9960763#map=13/53.8281/-1.5482

Can we get the constituent parts of that relation?

[Robinlovelace]

Being:

Way 145821681
Way Buck Stone Road (55027600)
Way 227675422
...

From the bottom of the page in the link above.

[agila5]

Hi Robin! I think that if you want to extract all elements that compose that relation, this is the easiest way:

# packages
library(osmextract)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright.
#> Check the package website, https://docs.ropensci.org/osmextract/, for more details.
library(sf)
#> Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUE

# data
ncn_leeds <- oe_get(
  "West Yorkshire", 
  force_vectortranslate = TRUE, 
  extra_tags = c("network", "ref"), 
  query = "SELECT * FROM multilinestrings WHERE network = 'ncn' AND ref = 668"
)
#> The input place was matched with: West Yorkshire
#> The chosen file was already detected in the download directory. Skip downloading.
#> The corresponding gpkg file was already detected. Skip vectortranslate operations.
#> Reading query `SELECT * FROM multilinestrings WHERE network = 'ncn' AND ref = 668' from data source `D:\osm-data\geofabrik_west-yorkshire-latest.gpkg' using driver `GPKG'
#> Simple feature collection with 1 feature and 6 fields
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -1.563765 ymin: 53.79892 xmax: -1.532757 ymax: 53.8572
#> Geodetic CRS:  WGS 84

i.e. we filter using the fields named “ncn” and “ref”. The plot is exactly the same plot shown here.

par(mar = rep(0, 4))
plot(st_geometry(ncn_leeds))

We can also plot it piece by piece, but unfortunately we cannot recover the OSM tags for each element that compose the relationship (since, AFAIK, the relevant links are stripped by Geofabrik and the other providers. In these cases, I suggest the use of osmdata).

ncn_leeds_lines <- st_cast(ncn_leeds, "LINESTRING")
#> Warning in st_cast.sf(ncn_leeds, "LINESTRING"): repeating attributes for all
#> sub-geometries for which they may not be constant
st_geometry(ncn_leeds_lines)
#> Geometry set for 60 features 
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -1.563765 ymin: 53.79892 xmax: -1.532757 ymax: 53.8572
#> Geodetic CRS:  WGS 84
#> First 5 geometries:
#> LINESTRING (-1.563703 53.8572, -1.563708 53.857...
#> LINESTRING (-1.557441 53.85267, -1.557564 53.85...
#> LINESTRING (-1.557441 53.85267, -1.557449 53.85...
#> LINESTRING (-1.555827 53.85138, -1.55572 53.851...
#> LINESTRING (-1.553986 53.84935, -1.553984 53.84...

We can clearly see that there are 60 elements as in the original relation, but the corresponding fields are meaningless. If you want, we can also plot it

plot(st_geometry(ncn_leeds_lines), col = sf.colors(60, categorical = TRUE), lwd = 2)

^{Created on 2022-09-03 with reprex v2.0.2}

[Robinlovelace]

This is really useful Andrea and, as I suspected, there seems to be no way to do it directly. Heads-up @husiejames re. OpenInfra.
I think a function that

• Gets the relation as a multi object
• Gets the objects in the bounding polygon around it
• Finds the sub-objects that compose the relation using spatial relations

Could be a possible work-around. I could have a go at writing such a function if you think it's a good idea for this package.

[agila5]

Just for future reference, the following code uses osmdata for the same purposes:

# packages
suppressPackageStartupMessages({
  library(osmdata)
  library(sf)
})

# select id and download data
ncn <- opq_osm_id(id = 9960763, type = "relation") |>
  osmdata_sf()

The object osm_lines containes each way and the relevant fields

ncn$osm_lines[, 1:2]
#> Simple feature collection with 60 features and 2 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -1.563765 ymin: 53.79892 xmax: -1.532757 ymax: 53.8572
#> Geodetic CRS:  WGS 84
#> First 10 features:
#>            osm_id                         name                       geometry
#> 10204524 10204524 Definitive Footpath LEEDS 94 LINESTRING (-1.555489 53.82...
#> 10204527 10204527            Scott Hall Street LINESTRING (-1.544284 53.81...
#> 12488960 12488960              Sugar Well Road LINESTRING (-1.55844 53.822...
#> 13267772 13267772                Farm Hill Way LINESTRING (-1.558057 53.82...
#> 14411253 14411253                         <NA> LINESTRING (-1.551293 53.83...
#> 22957171 22957171                Oatland Close LINESTRING (-1.54151 53.810...
#> 24866172 24866172                Oatland Place LINESTRING (-1.542262 53.81...
#> 25573765 25573765          King Alfred's Drive LINESTRING (-1.556258 53.83...
#> 54905405 54905405              Deanswood Drive LINESTRING (-1.555827 53.85...
#> 54905450 54905450                   Saxon Road LINESTRING (-1.557487 53.84...

while the complete relation is summarised in

ncn$osm_multilines[, 1:2]
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -1.563765 ymin: 53.79892 xmax: -1.532757 ymax: 53.8572
#> Geodetic CRS:  WGS 84
#>                    osm_id                   name                       geometry
#> 9960763-(no role) 9960763 NCN National Route 668 MULTILINESTRING ((-1.563703...

^{Created on 2022-09-04 with reprex v2.0.2}