Geospatial vector data: Part 1 OSM Query
Apr 26, 2022
Green Space
1 Load libraries
library(dplyr) # data wrangling
library(tidyr) # data wrangling
library(ggplot2) # data visualisation
library(sf) # simple features - geospatial geometries
library(osmdata) # obtaining OpenStreetMap vector data
library(units) # working with units
library(mapview) # interactive geometry viewing
library(ggmap) # downloading raster maps from a variety of sources
library(ggspatial) # map backgrounds and annotations for ggplot
library(tmap) # static/interactive map library with ggplot-like syntax2 OpenStreetMap query
2.1 OSM query
# Get the polygon for Edinburgh
city_polygon <- getbb("City of Edinburgh",
featuretype = "settlement",
format_out = "polygon")
# Get the rectangular bounding box
city_rect <- getbb("City of Edinburgh", featuretype = "settlement")# Get the data from OSM (might take a few seconds)
greensp_osm <-
opq(bbox = city_polygon) %>% # start query, input bounding box
add_osm_feature(key = "leisure",
value = c("park", "nature_reserve", "golf_course")) %>%
# we want to extract "leisure" features that are park, nature reserve OR a golf course
osmdata_sf() %>%
# query OSM and return as simple features (sf)
trim_osmdata(city_polygon)
# limit data to the Edinburgh polygon instead of the rectangular bounding box# Look at query output
greensp_osmThe query returns a list which contains multiple sf object, each for a different geometry type
list_object$polygons
list_object[["polygons"]]glimpse(greensp_osm$osm_polygons)
glimpse(greensp_osm$osm_multipolygons)As we have both POLYGON and MULTIPOLYGON features, it would be easiest to convert the POLYGON features to MULTIPOLYGON and then bind the two sf objects. (Polygons can easily be converted to multipolygons without change, while multipolygons may have to be split into multiple features to become polygons.)
# Convert POLYGON to MULTIPOLYGON and bind into one sf object.
greensp_sf <- bind_rows(st_cast(greensp_osm$osm_polygons, "MULTIPOLYGON"),
greensp_osm$osm_multipolygons) %>%
select(name, osm_id, leisure) # leisure (the type of green space).greensp_sf <- readRDS("greensp_sf.rds")head(greensp_sf)## Simple feature collection with 6 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -3.388398 ymin: 55.90319 xmax: -3.140106 ymax: 55.98513
## Geodetic CRS: WGS 84
## name osm_id leisure
## 4271400 Dundas Park 4271400 park
## 4288119 Harrison Park East 4288119 park
## 4348244 Bruntsfield Links 4348244 park
## 4891768 Baberton Golf Course 4891768 golf_course
## 4891786 Kingsknowe Golf Course 4891786 golf_course
## 4892551 <NA> 4892551 <NA>
## geometry
## 4271400 MULTIPOLYGON (((-3.387083 5...
## 4288119 MULTIPOLYGON (((-3.224193 5...
## 4348244 MULTIPOLYGON (((-3.203498 5...
## 4891768 MULTIPOLYGON (((-3.28719 55...
## 4891786 MULTIPOLYGON (((-3.265052 5...
## 4892551 MULTIPOLYGON (((-3.146338 5...unique(greensp_sf$leisure)## [1] "park" "golf_course" NA "nature_reserve"2.2 Data Cleaning
# Filter out unneeded shapes
greensp_sf <-
greensp_sf %>%
filter(is.na(leisure) == FALSE) %>%
# remove leisure NAs
rename(greensp_type = leisure) %>%
# rename leisure to greensp_type
st_make_valid()
# a good function to use after importing data to make sure shapes are valid2.3 Plot
plot(greensp_sf["greensp_type"])
2.4 Save / Load spatial data from file
sf objects can be saved in the form of spatial vector files easily using the st_write() function.
Here’s how it can be saved as a GeoPackage (.gpkg) format:
st_write(greensp_sf,
dsn = "greenspaces_Edi_OSM.gpkg", # file path
layer="greenspaces", # layer name
layer_options = c(paste0("DESCRIPTION=Contains spatial multipolygons for parks, ",
"nature reserves and golf courses in Edinburgh, Scotland. ",
"Copyright OpenStreetMap constibutors. ODbL ",
"https://www.openstreetmap.org/copyright")),
# add layer description
delete_dsn = TRUE
# to delete the whole file first, because sometimes, we can just
# overwrite or append one layer to an already existing file.
# If the file doesn't exist, this will return a friendly warning.
)## Deleting source `greenspaces_Edi_OSM.gpkg' using driver `GPKG'
## Writing layer `greenspaces' to data source
## `greenspaces_Edi_OSM.gpkg' using driver `GPKG'
## options: DESCRIPTION=Contains spatial multipolygons for parks, nature reserves and golf courses in Edinburgh, Scotland. Copyright OpenStreetMap constibutors. ODbL https://www.openstreetmap.org/copyright
## Writing 274 features with 3 fields and geometry type Multi Polygon.Read the file:
greensp_sf <- st_read(dsn = "greenspaces_Edi_OSM.gpkg", layer = "greenspaces" )## Reading layer `greenspaces' from data source
## `/Users/admin/Desktop/Linh Data Studio/Spatial Data/Spatial Modelling/greenspaces_Edi_OSM.gpkg'
## using driver `GPKG'
## Simple feature collection with 274 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -3.416874 ymin: 55.85796 xmax: -3.079656 ymax: 55.99143
## Geodetic CRS: WGS 843 Coordinate reference systems (CRS)
st_crs(greensp_sf)## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]Let’s transform our data into the British National Grid CRS, using it’s ESPG number (27700):
greensp_sf <- st_transform(greensp_sf, 27700)st_crs(greensp_sf)## Coordinate Reference System:
## User input: EPSG:27700
## wkt:
## PROJCRS["OSGB36 / British National Grid",
## BASEGEOGCRS["OSGB36",
## DATUM["Ordnance Survey of Great Britain 1936",
## ELLIPSOID["Airy 1830",6377563.396,299.3249646,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4277]],
## CONVERSION["British National Grid",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",49,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-2,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996012717,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",400000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",-100000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["United Kingdom (UK) - offshore to boundary of UKCS within 49°45'N to 61°N and 9°W to 2°E; onshore Great Britain (England, Wales and Scotland). Isle of Man onshore."],
## BBOX[49.75,-9,61.01,2.01]],
## ID["EPSG",27700]]