3 min read

Create custom maps from openstreetmap

by

Openstreetmaps is a gigantic collection of geographical information. Most often, it is used to create base maps for plots. This small tutorial shows how you can generate your own custom openstreetmap base map for ggplot2 and how to extract data from OSM.

Our goal is to create a map of Freetown with just the roads. In addition we will search openstreetmap for hospitals in Freetown and display them on a map as well.

All this is based on the excellent packages osmplotr and osmdata by Mark Padgham et al., that makes it incredibly easy to interact with openstreetmap.

Creating the basemap

First we define the bounding box around Freetown, then we query objects tagged with highway and add them to our map. The resulting object is a ggplot2 object and so we can further process it.

library(osmplotr)

# A bounding box of Freetown
bbox <- osmdata::getbb("Freetown")

# get highways by type
dat_highway_prim <- extract_osm_objects(key = "highway", value = "primary", bbox = bbox)
Sys.sleep(10) # be nice to overpass API
dat_highway_sec <- extract_osm_objects(key = "highway", value = "secondary", bbox = bbox)
Sys.sleep(10) # be nice to overpass API
dat_highway_tert <- extract_osm_objects(key = "highway", value = "tertiary", bbox = bbox)
Sys.sleep(10) # be nice to overpass API
dat_highway_res <- extract_osm_objects(key = "highway", value = "residential", bbox = bbox)
    
# now create the map
# more important roads are darker
map <- osm_basemap(bbox = bbox, bg = "white")
map <- add_osm_objects(map, dat_highway_res, col = "gray90")
map <- add_osm_objects(map, dat_highway_tert, col = "gray80")
map <- add_osm_objects(map, dat_highway_sec, col = "gray60")
map <- add_osm_objects(map, dat_highway_prim, col = "gray40")

map

Adding hospital data as points

We can query openstreetmap data with the package osmdata. In this example we look for all points taged with amenity = hospital. In this case, the data is a simple feature sf object that can be processed with familiar verbs, like filter or select.

library(magrittr)
library(osmdata)
# query the overpass API for hospitals
# and return a sf object
hospitals <- bbox %>%
  opq () %>% 
  add_osm_feature("amenity", "hospital") %>% 
  osmdata_sf()

library(dplyr)
library(sf)
# only use the points with type == "hospital"
hospitals <- filter(hospitals$osm_points, !is.na(name), type == "hospital")
select(hospitals, name, amenity, type)
## Simple feature collection with 9 features and 3 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -13.28041 ymin: 8.385998 xmax: -13.13049 ymax: 8.489861
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
##                                    name  amenity     type
## 1                  Blue Shield Hospital hospital hospital
## 2     China-SL Friendship Hospital, Jui hospital hospital
## 3                        DOCTOR QUARTER hospital hospital
## 4                          Goderich ETC hospital hospital
## 5                    Lakka Hospital ETU hospital hospital
## 6            Lumley Government Hospital hospital hospital
## 7    Police Training Sch-Hastings 1 ETC hospital hospital
## 8 Princess Christian Maternity Hospital hospital hospital
## 9          SL SACRIFICE HOSPITAL MENTAL hospital hospital
##                         geometry
## 1 POINT (-13.2466716766357 8....
## 2 POINT (-13.1427555084229 8....
## 3 POINT (-13.1937084197998 8....
## 4 POINT (-13.2804069519043 8....
## 5 POINT (-13.2639999389648 8....
## 6 POINT (-13.2624807357788 8....
## 7 POINT (-13.1304883956909 8....
## 8 POINT (-13.2193479537964 8....
## 9 POINT (-13.1938190460205 8....

The only thing now left to do is add the hospitals to the map and display the labels. Here we use ggrepel to prevent labels from overlapping.

library(ggrepel)
map <- add_osm_objects(map, hospitals, col = "red")

plot_hospitals <- as.data.frame(hospitals) %>% 
  rowwise() %>% 
  mutate(lat = geometry[1], lon = geometry[2]) 
map + 
  geom_label_repel(data = plot_hospitals, aes(lat, lon, label = name), size = 2)