In this note I focus on how to use the packages OSMDATA and SF. My aim is not to provide in depth descriptions of the packages and their functions, but instead to proceed as quickly as possible to cover the basics needed for wrangling spatial data in R.

The packages and their excellent vignettes are available here:

https://cran.r-project.org/web/packages/sf/index.html

https://cran.r-project.org/web/packages/osmdata/

The OSMDATA is used for getting data from OpenStreetMaps (OSM).

The SF package is used for data wrangling on datasets with geometric features.

The geometric features are primarily

  1. points
  2. lines
  3. polygons

The first step is to use the OSM-PACKAGE to get a spatial dataset by downloading data from OSM.

# OSMDATA EXAMPLE 1
library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
# The following code is based on this vignette
# https://cran.r-project.org/web/packages/osmdata/vignettes/osmdata.html

# In this code I use the package "osmdata" to download data from OSM
# And I investigate this data ... 

# Step (1): Find an area of interest
# The area is defined as a bounding box so get coordinates from OSM
# https://www.openstreetmap.org/export#map=15/55.6763/12.5705

# I choose Copenhagen.
# The bbox - square boundary box is defined in this way
# c(xmin,ymin,xmax,ymax) - where x is longitude and y is latitude

# Step (2): Build up the call to Overpass API
# First ingredient of the call is the box
call <- opq(bbox = c(12.4,55.5,12.8,55.9)) 

# It is possible to build call using name
# call <- opq(bbox="Copenhagen")
# And it is possible to check which boundary box is returned for the name
# getbb("Copenhagen")

# Second is the feature of interest
call <- add_osm_feature(call, key = "highway",value=c("motorway",
"primary","secondary"))

# Step (3): Make the query/call to Overpass API
#mydata <- osmdata_sf(call)
#mydata

When working with OSMDATA you can download the data from OSM within a bounding box. So the first step is to define the bounding box, which amount to defining the area you are interested in. There are two limitations in this

  1. The bounding box is really a box (not a circle, not a city border etc.)
  2. The bounding box has to be relatively small (not including an entire country)

For larger downloads of data from OSM you can use the package OSMEXTRACT which is described here

https://github.com/ITSLeeds/osmextract

If you want to have other shapes than a box, then there are techniques available but I do not go into these here. Here I simply choose to be satisified with working on an area defined by a box.

The easiest way to get a bounding box is simply to go to OpenStreetMaps and zoom in on the area of interest and read off latituede and longitude.

The second step is to add the osm-feature of interest by specifying a key and a value. To get an idea of the data that is available at OSM you can look at the Wiki-list

https://wiki.openstreetmap.org/wiki/Map_Features

and secondly you can use Overpass-turbo to see the results of a query to the Overpass API, which is how OSMDATA-package retrieves it information from OSM. The Overpass-turbo is available here

https://overpass-turbo.eu/.

Let us return to the example for which the code was

call <- opq(bbox = c(12.4,55.5,12.8,55.9)) 
call <- add_osm_feature(call, key = "highway",value=c("motorway",
"primary","secondary"))
mydata <- osmdata_sf(call)
mydata
## Object of class 'osmdata' with:
##                  $bbox : 55.5,12.4,55.9,12.8
##         $overpass_call : The call submitted to the overpass API
##                  $meta : metadata including timestamp and version numbers
##            $osm_points : 'sf' Simple Features Collection with 7592 points
##             $osm_lines : 'sf' Simple Features Collection with 1582 linestrings
##          $osm_polygons : 'sf' Simple Features Collection with 6 polygons
##        $osm_multilines : NULL
##     $osm_multipolygons : NULL

The result is an object of class osmdata which include geometric features in form of points,lines and polygons.

The result can be plotted using *

# The returned object stores some points, lines and polygons.
# Plot this,
# First get background map
library(ggmap)
mad_map <- get_map(getbb("Copenhagen"), maptype = "roadmap",source="osm",color="bw")

#Then plot:
ggmap(mad_map)+
  geom_sf(data = mydata$osm_lines,
          inherit.aes = FALSE)

Let us try to get subway stations. The key is station and the value is subway.

call <- opq(bbox = c(12.4,55.5,12.8,55.9)) 
call <- add_osm_feature(call, key = "station",value=c("subway"))
mydata <- osmdata_sf(call)
mad_map <- get_map(getbb("Copenhagen"), maptype="roadmap",source="osm",color="bw")

#Then plot:
ggmap(mad_map)+
  geom_sf(data = mydata$osm_points,
          inherit.aes = FALSE)

How do I save the data? One option is to save the dataset mydata as an R-object after which it can be loaded using. This is done as usual with the following commands

save(mydata,file="cph_metros")
load("cph_metros")

However another option is to use osmdata_xml and save it as an xml-object. This has the advantage that it can be loaded with programs like QGIS. The file has to be stored using the extension .osm. This can be done using the code

osmdata_xml(call,file="cph.osm")
readLines('cph.osm')[1:2]
## [1] "<?xml version=\"1.0\" encoding=\"UTF-8\"?>"                        
## [2] "<osm version=\"0.6\" generator=\"Overpass API 0.7.56.7 b85c4387\">"

which also read the first two lines of the file.

When the data has been stored in the xml-format it can also be read using the SF-package. But because osmdata object can include points, lines and polygons simultaneously the layer has to be selected

cph_data <- sf::st_read('cph.osm', layer = 'points')
## Reading layer `points' from data source `C:\Users\np83zg\OneDrive - Aalborg Universitet\Skrivebord\cars-master\rspace_osm\getdata\cph.osm' using driver `OSM'
## Simple feature collection with 41 features and 10 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 12.4925 ymin: 55.61937 xmax: 12.64922 ymax: 55.71184
## geographic CRS: WGS 84

We can the start to wrangle a little by simply dropping the columns that we are not interested in

cph_data <- cph_data[,c("osm_id","name","geometry")]
plot(cph_data)

Summary

Using the package OSMDATA relatively small datasets containing information from OSM in its current state can be generated. They can be generated in the format xml which can be used in for example QGIS and the layers can be read into R as sf-objects using the SF package.

A basic limitation is that this does not give the user access to large datasets generated from OSM. One solution for this limitation is the OSMEXTRACT package, however this package only delivers sf-objects. Another solution is to run a local version of OSM as is described in this tutorial https://www.r-spatial.org/2017/07/14/large_scale_osm_in_r. Another possibility is to use Python and the package PYROSM as described here https://pyrosm.readthedocs.io/en/latest/.