Search code examples
rdataframeopenstreetmapgeo

osmdata - how to extract to which city/country/district a point or a way belongs to?

My question is very similar to this one, I just want to learn how to do it in R with osmdata package (if it is possible). Let's say that we have pulled the data on all the primary schools in London.

bb <- getbb("London", featuretype = "city", format_out = "sf_polygon")

schools <- getbb("London", featuretype = "city") %>% 
  opq() %>% 
  add_osm_feature(key = "school", value = "primary") %>% 
  osmdata_sf() %>% 
  trim_osmdata(bb)

Now I want to know which London borough does each school belong to. How can I get to this data?

Solution

  • London boroughs were previously available through OIdata package, but it was removed from CRAN. You can find the shapefile here: https://data.london.gov.uk/dataset/statistical-gis-boundary-files-london

    library(osmdata)
library(sf)
library(tidyverse)

    bb <- getbb("London", featuretype = "city", format_out = "sf_polygon")

schools <- getbb("London", featuretype = "city") %>% 
  opq() %>% 
  add_osm_feature(key = "school", value = "primary") %>% 
  osmdata_sf() %>% 
  trim_osmdata(bb)

    ## get the schools as points
schools_pts <- schools[["osm_points"]] 

## reading the boroughs shapefile
path <- "C:\\Users\\M--\\Downloads\\Unzipped\\London_Borough_Excluding_MHW.shp"
ldn_brg <- st_read(path)

## Check the projection
st_crs(ldn_brg)
#> Coordinate Reference System:
#>   User input: OSGB 1936 / British National Grid 
## ...
st_crs(schools_pts)
#> Coordinate Reference System:
#>   User input: EPSG:4326 
## ...

## set the projections to be the same
ldn_brg <- st_transform(ldn_brg, crs = 4326)

## join the datasets; 
## you can find the borough of each school under NAME column
schools_pts %>% 
  st_join(ldn_brg)-> schools_brg

    ## showing a subset of the data for demonstration
schools_brg %>% 
  select(osm_id, NAME, geometry)

# Simple feature collection with 380 features and 2 fields
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -0.3404579 ymin: 51.39544 xmax: 0.1064406 ymax: 51.57038
# Geodetic CRS:  WGS 84
# First 10 features:
#            osm_id      NAME                    geometry
# 32861089 32861089    Newham  POINT (0.0467075 51.54989)
# 32861229 32861229    Newham  POINT (0.0466267 51.54969)
# 32861230 32861230    Newham  POINT (0.0472586 51.54955)
# 32861231 32861231    Newham  POINT (0.0472972 51.54961)
# 32861232 32861232    Newham   POINT (0.047325 51.54967)
# 32861233 32861233    Newham  POINT (0.0471958 51.54985)
# 32861234 32861234    Newham  POINT (0.0472063 51.54987)
# 32861235 32861235    Newham  POINT (0.0469997 51.54991)
# 32861237 32861237    Newham  POINT (0.0469506 51.54983)
# 49500567 49500567 Islington POINT (-0.0916688 51.53894)

    ## plotting for demonstration purposes only
ggplot() + 
  geom_sf(data = ldn_brg, size = 1, color = "black") + 
  geom_sf(data = schools_pts, size = 2, color = "red") + 
  coord_sf()

    Created on 2022-10-28 by the reprex package (v2.0.1)