Polygons as lines

[luukvdmeer]

Hi! A small question regarding the extraction of a street network from a pbf file. In omsdata, there is the function osm_poly2line to convert polygons into linestrings. This is needed because circular street segments become polygons. Is that behaviour the same when using osmextract? And if yes, what would be the best way to exract a street network including circular street segments?

[agila5]

Hi Luuk! I'm not 100% sure (since OSM data may be pretty weird and inconsistent), but I think that the data obtained with osmextract do not require the same preprocessing steps like the ones obtained using osmdata. Hence, osm_poly2line is (usually) not required (since, by default, most "closed ways" are saved as lines instead of polygons). I try to summarise an example with the following (quite long) reprex, hope it's clear:

Load packages

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1
library(osmdata)
#> Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
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.

Define a bbox

olginate_bbox <- st_bbox(
  c(xmin = 9.408513, ymin = 45.799540, xmax = 9.410992, ymax = 45.801222), 
  crs = 4326
) %>% 
  st_as_sfc()

Download highway data with osmdata

my_osmdata <- opq(olginate_bbox) %>% 
  add_osm_feature(key = "highway") %>% 
  osmdata_sf()

Check output

my_osmdata
#> Object of class 'osmdata' with:
#>                  $bbox : 45.79954,9.408513,45.801222,9.410992
#>         $overpass_call : The call submitted to the overpass API
#>                  $meta : metadata including timestamp and version numbers
#>            $osm_points : 'sf' Simple Features Collection with 76 points
#>             $osm_lines : 'sf' Simple Features Collection with 15 linestrings
#>          $osm_polygons : 'sf' Simple Features Collection with 1 polygons
#>        $osm_multilines : NULL
#>     $osm_multipolygons : NULL

We can see that there is 1 polygon and I know that it corresponds to a roundabout. Let’s plot it.

par(mar = rep(0, 4))
plot(st_geometry(my_osmdata$osm_lines), reset = FALSE, lwd = 3)
plot(st_geometry(my_osmdata$osm_polygons), border = "red", add = TRUE)

Convert that poly to a line

(my_osmdata <- my_osmdata %>% osm_poly2line())
#> Object of class 'osmdata' with:
#>                  $bbox : 45.79954,9.408513,45.801222,9.410992
#>         $overpass_call : The call submitted to the overpass API
#>                  $meta : metadata including timestamp and version numbers
#>            $osm_points : 'sf' Simple Features Collection with 76 points
#>             $osm_lines : 'sf' Simple Features Collection with 16 linestrings
#>          $osm_polygons : 'sf' Simple Features Collection with 1 polygons
#>        $osm_multilines : NULL
#>     $osm_multipolygons : NULL

Now we can see that the roundabout is defined as a circular line

par(mar = rep(0, 4))
plot(st_geometry(my_osmdata$osm_lines), reset = FALSE)

Same stuff with osmextract (much slower in this case, set quiet = FALSE for progress bar)

my_osmextract <- oe_get(
  place = olginate_bbox, 
  boundary = olginate_bbox, 
  query = "SELECT * FROM lines WHERE highway IS NOT NULL", 
  quiet = TRUE
)

Again, plot

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

We can see that the roundabout (i.e. the circular line previously coded as a polygon) is alreadyrepresented as a line. Hence, I think we don’t need any poly2line function in this case. We can modify the osmconf.ini file and process closed_ways (i.e. the circular lines representing roundabouts) as polygons. For example:

config <- readLines(system.file("osmconf.ini", package = "osmextract"))
config[7] <- paste(config[7], "junction=roundabout", sep = ",")
temp_ini = tempfile(fileext = ".ini")
writeLines(config, temp_ini)

my_osmextract_lines <- oe_get(
  olginate_bbox, 
  boundary = olginate_bbox, 
  query = "SELECT * FROM lines WHERE highway IS NOT NULL", 
  osmconf_ini = temp_ini, 
  quiet = TRUE
)
my_osmextract_poly <- oe_get(
  olginate_bbox, 
  boundary = olginate_bbox, 
  layer = "multipolygons",
  osmconf_ini = temp_ini, 
  query = "
  SELECT *, hstore_get_value(other_tags, 'junction') AS junction 
  FROM multipolygons 
  WHERE osm_way_id IS NOT NULL AND junction = 'roundabout'
  ", 
  quiet = TRUE
)

and now

par(mar = rep(0, 4))
plot(st_geometry(my_osmextract_lines), reset = FALSE)
plot(st_geometry(my_osmextract_poly), border = "red", add = TRUE)

Not sure about the benefits. Also not sure if the same holds for bbbike and openstreetmap_fr providers but I think so.

^{Created on 2021-10-05 by the reprex package (v2.0.1)}