Skip to content

Instantly share code, notes, and snippets.

@agila5
Last active September 28, 2021 15:35
Show Gist options
  • Save agila5/219445de4fb5c26b832dbb9cbb7414cd to your computer and use it in GitHub Desktop.
Save agila5/219445de4fb5c26b832dbb9cbb7414cd to your computer and use it in GitHub Desktop.
# 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.

# bounding poly
bb <- st_bbox(c(xmin = 9.20, ymin = 45.51, xmax = 9.21, ymax = 45.52), crs = 4326)

# get all road segments in that bb using osmdata
my_osmdata <- opq(bb) %>% 
  add_osm_feature(key = "highway") %>% 
  osmdata_sf()

# get all road segments in that bb using osmextract (warning: it might take some
# time, especially if you need to download the pbf file which is quite huge.
# osmdata is typically much better when you need to work with small areas)
my_osmextract <- oe_get(
  place = bb, 
  boundary = st_as_sfc(bb), 
  query = "SELECT * FROM lines WHERE highway IS NOT NULL", 
  quiet = TRUE
)

# Compare the number of rows
nrow(my_osmdata$osm_lines)
#> [1] 483
nrow(my_osmextract)
#> [1] 486

# Plot the result obtained with osmdata
par(mfrow = c(1, 2), mar = rep(0, 4))
plot(st_geometry(my_osmdata$osm_lines))

# Plot the result obtained with osmextract
plot(st_geometry(my_osmextract))

Created on 2021-09-28 by the reprex package (v2.0.1)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment