Using Web Services to Work with Geodata in R

Through collaborative mapping, a massive amount of data is accessible. Many individuals contribute information each day. The growing amount of geodata is gathered by volunteers or obtained via crowd-sourcing. One outstanding example of this is the OpenStreetMap (OSM) Project which provides access to big data in geography. Another online mapping service that enables the integration of geodata into the analysis is Google Maps. The expanding content and the availability of geographic information radically changes the perspective on geodata (Chilton (2009)). Recently many application programming interfaces (APIs) have been built on OSM and Google Maps. That leads to a point where it is possible to access sections of geographical information without the usage of a complex database solution, especially if one only requires a small data section for a visualization.
First tools for spatial analysis have been included in the R language very early (Bivand and Gebhardt 2000) and this development will continue to accelerate, underpinning a continual change. Notably, in recent years many tools have been developed to enable the usage of R as a geographic information system (GIS). With a GIS it is possible to process spatial data. QuantumGIS (QGIS) is a free software solution for these tasks, and a user interface is available for this purpose. R is, therefore, an alternative to geographic information systems like QGIS (QGIS Development Team (2009)). Besides, add-ins for QGIS and R-packages (RQGIS) are available, that enables the combination of R and QGIS (Muenchow and Schratz (2017)). It is the target of this article to present some of the most important R-functionalities to download and process geodata from OSM and the Google Maps API. The focus of this paper is on functions that enable the natural usage of these APIs.

Jan-Philipp Kolb (GESIS Leibniz Institute for the Social Sciences)
2019-07-30

1 Introduction and outline

This paper introduces some interesting web services for downloading, processing and visualizing geodata. The focus especially in the second half of the paper is on OpenStreetMap-data, because it is released under the Open Database License (ODbL) 1.0. That allows multiple uses of the data (Schmidt et al. 2013). The study of (Barrington-Leigh and Millard-Ball 2017) shows for example, that the data quality available at OSM is already sufficient in many countries to use it for scientific and analytic purposes. However, (Barron et al. 2014) state that the quality of the OSM-data depends on the individual use case. And (Grippa et al. 2018) mention that it is essential to consider the variations at regional or national scales. One example of a scientific analysis based on OSM-data is the Simulation of Urban MObility (SUMO) project (Behrisch et al. (2011)). (Meijer et al. 2018) for example use OSM-data to analyze global patterns of road infrastructure. (Gervasoni et al. 2018) use OSM-data to generate urban features that help to estimate population density at a higher resolution. (Arsanjani et al. 2015) give an overview of typical and recent examples of studies done with OSM-data. Much more research, carried out in various countries, is listed at (OpenStreetMap Wiki 2017e).

The focus is on the most important APIs to download geodata. The significant advantage of using these specific APIs is that we can obtain data free of charge. Short examples are used to describe how the data can be imported into R and processed. Some examples show the easiest and fastest way to get the information needed. In other examples I look a bit further behind the scenes. Static maps can be used as background information for geographic visualization and may be used to highlight positions of so-called points of interest (poi). A prerequisite to visualise these points is the availability of their exact spatial location. With the Overpass API (http://wiki.openstreetmap.org/wiki/Overpass_API) for example, we can get the positions for many points of interest. This application programming interface (API) is perfect to download data on very particular topics. For example, if you are looking for special map features.

The used API’s are listed in the individual sections below. I discuss an example where I am interested in public transportation in Amsterdam. In the next section (Background Maps - Download via Map Tile Servers), hints on the download of static maps from so-called map tile servers are presented. In the third section (Geocoding with Application Programming Interfaces (APIs) the functionality of APIs like the Google Maps and OSM Nominatim API is used to realize geocoding. It is shown, how the Nominatim API can be used to search OSM-data by name and address (OpenStreetMap Wiki (2018b)). In the fourth section (Downloading and Importing OSM-data) I show various possibilities to download more general OSM-data. The usage of the main OSM-API is presented as well as some functions of the osmdata package, which also uses the Overpass API, are described in this section. Possibilities to process OSM-data with R are presented in the fifth section (Processing OSM-data). A summary follows at the end.

2 Background maps - download via Map Tile Servers

If a background map is needed, it is possible to use a tile server to download them. Map tiles are quadratic bitmap graphics which are arranged in a grid to show a map. The vector tile is a newer format developed recently which is for example used by Mapbox (https://www.mapbox.com/). Vector tiles have a vector representation (OpenStreetMap Wiki (2018d 1)). The tiles contain vector data instead of the rendered image and provide readable, descriptive, and extensible content (Li et al. (2018)). Vector tiles can be rendered dynamically and allow for an efficient extraction of the relevant data (Gaffuri (2012 94)).

So-called map tile servers offer to download static maps of various types. It is, for example, possible to get maps on such diverse issues as biking, public transportation, or land shading.1 Map Tiles are very suitable for the use as background image. Various R-packages can be used to access map tile servers. One way to get static maps is the package OpenStreetMap. It is a package to access high-resolution raster maps using the OSM protocol (Fellows 2016). A high number of satellite, topographic and road map servers can be accessed directly using the JMapViewer Java component (Stotz 2018). The used map servers are for example CloudMade, Mapnik, Bing, Stamen, and MapQuest. The function openmap can be used to retrieve a map. It is necessary to provide values for the upper left latitude and longitude value as well as for the lower right values. In the example below, this is done for some coordinates in Amsterdam. Also, we have to specify the type of source. That may be the tile server from which to get the map or the uniform resource locator (URL) pattern. However, OSM servers have limited capacity, and heavy use adversely affects the purpose of use. With the package OpenStreetMap, it is also possible to access other web services. Bing Maps, the web mapping service provided by Microsoft is one example. In the following code example, the function openmap is used to get a map based on latitude and longitude coordinates.

We need a geocode to get a map of Amsterdam. In the next section, geocodes will be explained in more detail. We specify the tile server with the argument type. In this example, OSM is chosen, but a Bing map would also be possible. The result of this call is visible in Figure 1.

library("OpenStreetMap")
map <- openmap(c(52.278174, 4.729242),
               c(52.431064, 5.079162),
               type = "osm")
plot(map)
graphic without alt text
Figure 1: Map of Amsterdam - data from OSM

Stamen is an alternative source. Stamen Design publishes maps under a Creative Commons license CC BY-3.0 (Attribution). The maps are based on OSM-data (Lamigueiro (2014 95)). We get a Stamen map when we add further arguments to the openmap call. In the following the source is stamen. The type was specified as toner and watercolor. The resulting Stamen maps are depicted in Figure 2. The downloaded maps are very suitable as background for info graphics. It is possible to add further layers using, for example, the ggplot2 framework (Wickham (2009)). That will be shown later.

map_stt <- OpenStreetMap::openmap(c(52.385914, 4.874383), c(52.35514, 4.92054),
              type = "stamen-toner")
map_st <- OpenStreetMap::openmap(c(52.278174, 4.729242), c(52.431064, 5.079162),
              type = "stamen-watercolor")
plot(map_st)
plot(map_stt)

graphic without alt textgraphic without alt text

Figure 2: Static stamen maps of Amsterdam

Another package to get static maps is ggmap (Kahle and Wickham (2013)). This package provides a collection of functions to visualize spatial data and models on top of static maps from various online sources, like Google Maps, OSM, Stamen Maps and CloudMade Map (Kahle and Wickham 2013). Only a few lines of code are necessary to get a map for a freely selectable location. The default source of ggmap is the Google Static Maps API, and with the download of these images, you agree to the terms of usage (https://developers.google.com/maps/terms - Dorman (2014)). Recently, the Google Maps API terms of use have changed. Now you need an account to use the API for downloading a static map. The development version on Github (https://github.com/dkahle/ggmap) already has the function register_google where you can define your key. Previously you have to register your project at https://cloud.google.com/maps-platform/. The function qmap is a wrapper for ggmap and get_map. It is necessary to specify the place for which the map should be downloaded and a zoom factor, whereas the zoom parameter takes values between three and 21. A whole continent is visible on the map for a zoom factor of three whereas only one building is on the map for a zoom factor of 21. The default value is ten. In this case, a city is visible on the map.

The RgoogleMaps package can be used for querying OSM servers for static maps in the form of portable network graphics (PNGs) (Loecher and Ropkins 2015). Map tiles can be downloaded using, for example, the command GetMap.bbox. In this case, the center and a zoom level have to be provided. The center is determined using so-called geocodes. Geocodes are used to specify a precise location on the map. More information about these codes will be given in the next section.

Also, it is possible to create interactive maps with online mapping services. This can be done for example with the R-package leaflet created by (Cheng and Xie 2016). The package can be used to create interactive maps for websites. These kinds of maps are also known as slippy maps where it is possible to zoom and pan (OpenStreetMap Wiki (2016b)). That means that the map slips around when you drag the mouse. Slippy maps in OSM are based on the AJAX library OpenLayers which is written in JavaScript. Here, the default OSM tiles can be added to the interactive map visualization. It is also possible to use Stamen maps or CartoDB as background for slippy maps (Abernathy (2016 311)). A good starting point for the work with this package is https://rstudio.github.io/leaflet/. In the following example, the pipe-forward operator of package magrittr is used which enables chain operations. The first operation in this chain is the creation of a leaflet map widget (function leaflet). The second operation is to add a layer (function addTiles), and in the last operation a marker is added (function addMarkers). Here we have to specify the position of the marker and the text that pops up.

library(leaflet)
leaflet() %>%
  addTiles() %>%  
  addMarkers(lng=4.891013, lat=52.38054, popup = "Amsterdam")

Map Tiles are a good possibility for geographic data visualization. They can be valuable to get a first impression, but these visualizations may also occlude relevant geodata. Therefore it is useful to know, how to add further information to the map. In the next section, it is shown how to append more information to either the static map or the interactive map. Interactive maps can for example be produced very easily with tmap and mapview (Appelhans et al. (2018)) and many other packages. The functionalities for interactive maps in both mentioned packages are built on top of the R-package leaflet. The interface to Javascript allows a very vital exchange. The R-package mapdeck is for example a very suitable tool, for a browser based visualization of geodata. Like most of the interactive graphics in R this package is also based on a Javascript library. In this case Mapbox GL JS is used (Cooley (2018b)). Like for package deckard (see Hansel (2018)), the package also provides access to the Deck.gl framework of Uber (Lovelace et al. (2018)). A registration is necessary to use the framework. The package lawn provides a client for the geospatial analysis with the javascript library Turfjs (Chamberlain and Hollister (2017)). That allows us to use for example Javascript libraries like geojson-random and geojsonhint, which can be used to randomly create GeoJSON objects and to color them. For example the gr_polygon function can be used to create an example object. Then you can plot the object with the generic function view.

We have seen in the examples above, that we need a set of coordinates to locate a point of interest or pop-up. For a study on the transport system, it is for example good to know which public transport stops (train or bus stops, etc. ) are located in the surrounding of the area under research. The geocoding for an address list of these stops can be done with the Nominatim API of the OSM project. In the next section it is shown how this geocoding process is done in R.

3 Geocoding with Aplication Programming Interfaces (APIs)

Geocoding is the derivation of a structured spatial representation from textual information like postal codes (Aitchison 2009 157). Many possibilities are available to realize geocoding with R. The most popular is perhaps the usage of interfaces like the Google Maps API. The API is described in (Svennerberg 2010). It can be accessed directly from R using the R-package ggmap (Kahle and Wickham 2013 156). The ggmap-package was one of the first R packages to provide an interface for data exchange between R and the Google Maps API. The process can be implemented using the geocode function. Then we just need an adress to get the corresponding coordinates. For the example "Waterlooplein 1, Amsterdam, Zentrum" we get the latitude and longitude coordinates visible in  Table 1 as a result of this call.

Table 1: Latitude and longitude coordinates of the address “Amsterdam, Waterlooplein 1”
lat lon
1 4.901323 52.36896

The Mercator-projection is used in Google Maps (Moore and Drecki (2008 206)). The European Petroleum Survey Group Geodesy (EPSG) published a system of globally unique key numbers of geodetic data records (EPSG codes). The used coordinate reference system and the projection are determined using this EPSG codes (http://epsg.org/). The EPSG code used in this example is 3857 (Harris (2016)). One of the issues when using the Google Maps API is that only 2,500 requests per day can be performed free of charge, which might cause problems with large data sets. The terms of usage of the API can be seen at https://developers.google.com/maps/terms. It is clear that nobody should use such an online service to geocode privacy sensitive data. The googleway package also connects to Google (Cooley (2018a)). To use the package, the registration and an API key is necessary to use most of the functionalities. When we have done this, it is for example possible to query the distance, the elevation or the timezone. Additional information on the visinity of a point is accessible using for example the packages geonames (https://github.com/ropensci/geonames) or RDSTK.

A less known alternative described here in more detail is the usage of the Nominatim API of the OSM project (Warden 2011 25). This tool is an open source system designed on top of OSM-data to search by name and address (Clemens 2015). Nominatim (OpenStreetMap Wiki (2018b)) is the main geocoder maintained by OSM (Abernathy 2016). Detailed information on this geocoder is available at http://wiki.openstreetmap.org/wiki/Nominatim. Similar to the geocoding with the Google Maps API, the service Nominatim can be used to query the name of the reference object to obtain corresponding GPS coordinates. When using the Nominatim API, it is possible to choose between different output formats. We can choose between HTML, XML, JSON and JSONV2. The RJSONIO and jsonlite packages can be used to import JSON files to R (Lang (2014) and Ooms (2014)). The core from the following example is the command fromJSON of package RJSONIO. It converts JSON content to R objects. The code chunk is designated to get the corresponding coordinates for the address “Rozengracht 1” in Amsterdam. With the function url we specify a path to be opened. In this case it is the adress of the Nominatim API http://nominatim.openstreetmap.org/ plus some additional information (format and adress details).

library("RJSONIO")
con <- url("http://nominatim.openstreetmap.org/search?format=json&
addressdetails=1&extratags=1&q=Amsterdam+Niederlande+Rozengracht+1")
geoc <- fromJSON(paste0(readLines(con, warn=F)))
close(con)

The result is an object that contains a lot of information. We can get an overview when we query the names of the first object in geoc.

names(geoc[[1]])
 [1] "place_id"     "licence"      "osm_type"     "osm_id"       "boundingbox" 
 [6] "lat"          "lon"          "display_name" "class"        "type"        
[11] "importance"   "address"      "extratags"   

This object contains for example information on the license which is ODbL 1.0 (http://www.openstreetmap.org/copyright). We get the latitude and longitude coordinates:

geoloc <- c(geoc[[1]][which(names(geoc[[1]]) == "lat")],
            geoc[[1]][which(names(geoc[[1]]) == "lon")])
lat lon
1 52.3737223 4.8826404

We also get the combination of OSM id and the OSM type, which can be useful information, when downloading and processing specific OSM-data. And then we have some key-value pairs, which will be explained below. The package jsonlite can also be used for importing json data:

link <- url("http://nominatim.openstreetmap.org/search?format=json&
           addressdetails=1&extratags=1&q=Amsterdam+Niederlande+Rozengracht+1")
geoc2 <- jsonlite::fromJSON(link)
geoc2df <- with(geoc2, data.frame(osm_id, lat,lon))
geoc2df$house_number <-  geoc2$address$house_number

The dataframe geoc2df contains several addresses in the "Rozengracht" street that all start with a one. Parts of the adresses are visible in the following table.

Table 2: Adresses for house number with one - "Rozengracht"
osm_id lat lon house_number
1 2721815875 52.3737223 4.8826404 1
2 2743624072 52.3719482 4.8755534 237-1
3 2721830930 52.3736673 4.8823914 7-1
4 2721827922 52.3734021 4.8813371 53-1
5 2721824637 52.372232 4.8767542 231-1
6 2721823434 52.3724786 4.8776618 187-1
7 2721820122 52.3727335 4.8786657 137-1
8 2721816644 52.3729874 4.8797588 105E-1
9 2720971311 52.3727658 4.8775263 194-1
10 2720971056 52.3728019 4.8775994 184-1

The package tmaptools offers tool functions to supply the workflow to create thematic maps. It provides the function geocode_OSM which is a wrapper for geocoding using the OSM Nominatim API (Tennekes (2018)). A bounding box can be created with the tmap::bb command. The functions tmaptools::bb_poly, and osmdata::getbb are also worth mentioning here, in particular with regard to extracting bounding polygons rather than mere boxes. In the following example for function geocode_OSM of package tmaptools we get only the coordinates as output because the default value of the argument details is FALSE.

library("tmaptools")
gc_tma <- geocode_OSM("Amsterdam, Buiten Brouwersstraat")
Table 3: Result of the geocoding request for a poi in Amsterdam
Info Nominatim Google Maps
CRS EPSG:4326 EPSG:3857
Projection Mercator Mercator
Longitude 4.891013 4.900478
Latitude 52.380541 52.36859

The result of the request for a postal address in Amsterdam is visible in Table 3 in the second column. A big difference becomes apparent when we compare the result for Nominatim and the result for Google Maps API (third column). The projection and the coordinate reference system (CRS) may be of great importance (Brown 2016). The projection is used to display the three-dimensional earth. Typically, this is a projection onto a two-dimensional map display. Google Maps uses, for example, the Mercator projection, which is good for zoomed-in viewing, but causes distortions when zoomed out (Turner (2006)). For the Nominatim-query we get EPSG:4326 (Maier (2014)).

The projection of the data is often necessary for the work with geodata from different sources, and this is true if we want to combine the gained information with other geodata, for example, static maps or satellite pictures. In the following example a transformation is shown. In a first step we create a data.frame which is called poi. In a second step we use the function coordinates to set spatial coordinates and to create a spatial object. Then we set the projection attributes with the command proj4string. In this case we use the epsg projection 4326. Afterwards it is possible to transform the spatial points using the function spTransform.

library(sp)
poi <- data.frame(lat = gc_tma$coords["x"],
                  lon = gc_tma$coords["y"])
sp::coordinates(poi) <- c("lat", "lon")
sp::proj4string(poi) <- sp::CRS("+init=epsg:4326")
res <- spTransform(poi, CRS("+init=epsg:3035"))

We get the following numbers for the coordinates:

res@coords
      lat     lon
  3973434 3264547

A clearer alternative to realize the coordinate reference system transformation is included in package sf:

library(sf)
poi2 <- st_sfc (st_point(gc_tma$coords), crs = 4326)
res2 <- st_transform (poi2, crs = 3035)

Further functions are available to switch between reference systems. The OpenStreeetMap package has the openproj function to translate from Mercator to another coordinate reference system and can be used to create ggplot2 and base graphics. (Lovelace et al. 2017) present a possibility to transform the reference system. A whole section in the book of (Pebesma and Bivand 2019) is dedicated to this topic. (Brown 2016) also explains how to work with map projections and coordinate reference systems in R. More information is available in (Plant et al. 2012). Further, the mapmisc package provides projection capabilities and utilities for producing maps (Brown (2016)). The function projection supplies information on the coordinate reference system.

The accessed geocodes might be combined with a static map in the next step. As already shown, the openmap function from package OpenStreetMap can be utilized to download a background map. We have to re-project the original OSM map, for example with the function openproj from package OpenStreetMap. The static map can then be combined with the extracted geocode. The result is visible in Figure 3.

poi <- data.frame(lon = gc_tma$coords["x"],
                  lat = gc_tma$coords["y"])
adm_map <- openproj(map_stt,
                    projection = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
library(ggmap)
autoplot(adm_map) + geom_point(aes(x = lon, y = lat), data = poi,size=5, col="red")