Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Plot city streets
library(sf)
library(osmdata)
location = "Tel-Aviv, Israel"
# Get data
q = opq(bbox = location)
dat = add_osm_feature(q, key = "highway")
dat = osmdata_sf(dat)
lines = dat$osm_lines
pol = dat$osm_polygons
pol = st_cast(pol, "MULTILINESTRING")
pol = st_cast(pol, "LINESTRING")
lines = rbind(lines, pol)
# Crop
bb = q$bbox
bb = strsplit(bb, ",")[[1]]
bb = as.numeric(bb)
bb = bb[c(2, 4, 1, 3)]
names(bb) = c("xmin", "xmax", "ymin", "ymax")
bb = st_bbox(bb, crs = 4326)
bb = st_as_sfc(bb)
lines = st_intersection(lines, bb)
lines = st_transform(lines, 32636)
# Plot
plot(lines[, "highway"], key.pos = NULL)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment