Introduction

This tutorial introduces geo-spatial data visualization in R. The entire R markdown document for this tutorial can be downloaded here.

This tutorial is based on R. If you have not installed R or are new to it, you will find an introduction to and more information how to use R here. For this tutorials, we need to install certain packages from an R library so that the scripts shown below are executed without errors. Before turning to the code below, please install the packages by running the code below this paragraph. If you have already installed the packages mentioned below, then you can skip ahead ignore this section. To install the necessary packages, simply run the following code - it may take some time (between 1 and 5 minutes to install all of the libraries so you do not need to worry if it takes some time).

# clean current workspace
rm(list=ls(all=T))
# set options
options(stringsAsFactors = F)         # no automatic data transformation
options("scipen" = 100, "digits" = 4) # suppress math annotation
op <- options(gvis.plot.tag='chart')  # set gViz options
# install libraries
install.packages(c("RgoogleMaps", "ggmap", "mapproj", "sf",
                   "dplyr", "OpenStreetMap", "devtools", "DT"))
# install package from github
devtools::install_github("dkahle/ggmap", ref = "tidyup")

Depending on the maps that are used in the visualization, it may also be necessary to access other data bases. One very useful data base for maps is, of course, Google Maps. However, to access Google Maps materials, installation and setting up other pieces of software is necessary. How to get access to Google’s data is discussed below. In the following section, methods that do not require installation of software other than R.

1 Getting started with maps

The most basic way to display geospatial data is to simply download and display a map. In order to do that, we load the libraries necessary for extracting and plotting the map.

The package OpenStreetMap offers a range of maps with different features. To access the OpenStreetMap data base, it is necessary to install the package. Once the package is installed, we can simply extract the map and define the region we want to plot by defining the longitude and latitude of the upper left and lower right corner of the region we want to display. The argument minNumTiles defines the accuracy of the map, the higher the number of tiles, the higher the resolution. The type of map is defined by the type argument. The type argument defines from which server the map is extracted. Once we have extracted a map, we can plot it using the “plot” function.

# load library
library(OpenStreetMap)
# extract map
AustraliaMap <- openmap(c(-8,110),
    c(-45,160),
#   type = "osm",
#   type = "esri",
    type = "nps",
    minNumTiles=7)
# plot map
plot(AustraliaMap)

In order to obtain different map types, we change the type argument. The following options are available for type:

# load package 
library(DT)
opt <- c("osm", "osm-bw","maptoolkit-topo", "waze", "bing", "stamen-toner", "stamen-terrain", "stamen-watercolor", "osm-german", "osm-wanderreitkarte", "mapbox", "esri", "esri-topo", "nps", "apple-iphoto", "skobbler", "hillshade", "opencyclemap", "osm-transport", "osm-public-transport", "osm-bbike", "osm-bbike-german")
opt <- data.frame(opt)
datatable(opt, rownames = FALSE, options = list(pageLength = 5, scrollX=T), filter = "none")

Unfortunately, not all options work. If they do not work, then an error message is shown telling us that the number of tiles is not supported.

We can zoom in or out by either changing the “zoom” or the “minNumTiles” arguments - in both cases, the higher the number, the more fine-grained the dispalyed map. Let’s check out some examples for maps of Queensland.

# extract map
queensland1 <- openmap(c(-8,135),
    c(-30,160),
    type = "osm",
    minNumTiles=6)
queensland2 <- openmap(c(-8,135),
    c(-30,160),
    type = "esri",
    minNumTiles=6)
# plot maps
par(mfrow = c(1, 2)) # display plots in 1 row/2 columns
plot(queensland1); plot(queensland2); par(mfrow = c(1, 1)) # restore original settings

The leaflet function from the leaflet package creates a Leaflet map widget using html-widgets. The widget can be rendered on HTML pages generated from R Markdown, Shiny, or other applications. The advantage in using this function lies in the fact that it offers very detailed maps which enable zooming in to very specific locations.

# load package
library(leaflet)
# load library
m <- leaflet() %>% setView(lng = 153.05, lat = -27.45, zoom = 12)
# display map
m %>% addTiles()

Generating maps with rworldmap and base R

Another data base that is very useful when certain maps is the rworldmap package. The rworldmap package contains the shape files for countries but also more fine grained-shape files that display the states of selected countries. The most basic data, however, simply represents the shapes of the countries in the world.

Using the worldmap package has the advantage that one is not dependent on third parties and their servers but can operate within R without being denied access due to e.g. copy right issues or server maintenance.

# load library
library(rworldmap)
# get map
worldmap <- getMap(resolution = "coarse")
# plot world map
plot(worldmap, col = "lightgrey", 
     fill = T, border = "darkgray",
     xlim = c(-180, 180), ylim = c(-90, 90),
     bg = "aliceblue",
     asp = 1, wrap=c(-180,180))

The basic map shown above can then be modified and enriched with color coding to convey geospatial data. The following shows how to customize the world map.

# load library
library(maps)
# plot maps
par(mfrow = c(1, 2)) # display plots in 1 row/3 columns
# show map with Latitude 200 as center
map('world', xlim = c(100, 300))
# add axes
map.axes()
# show filled map with Latitude 200 as center
ww2 <- map('world', wrap=c(0,360), plot=FALSE, fill=TRUE)
map(ww2, xlim = c(100, 300), fill=TRUE)

par(mfrow = c(1, 1)) # restore original settings

Generating maps with rnaturalearth and ggplot2

We can also use the data provided by the rnaturalearth and the rnaturalearthdata and use ggplot function from the ggplot2 package as well as the sf package to create very nice visualizations of geospatial data. The advantage over using rworldmap and base R lies in the fact that the code is easier to interpret and the visualizations are more appealing.

# load packages
library(ggplot2)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
# load data
world <- ne_countries(scale = "medium", returnclass = "sf")
# gene world map
ggplot(data = world) +
  geom_sf() +
  labs( x = "Longitude", y = "Latitude") +
  ggtitle("World map", subtitle = paste0("(", length(unique(world$admin)), " countries)"))

We can also easily zoom in on certain areas in the map using the coord_sf function and also prettify the map by adding some custom features such as a compass.

library(rgeos)
library(ggspatial)
# gene world map
ggplot(data = world) +
  geom_sf() +
  labs( x = "Longitude", y = "Latitude") +
  coord_sf(xlim = c(100.00, 160.00), ylim = c(-45.00, -10.00), expand = FALSE) +
  annotation_scale(location = "bl", width_hint = 0.5) +
  annotation_north_arrow(location = "bl", which_north = "true", 
                         pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
                         style = north_arrow_fancy_orienteering) +
  theme_bw()

We will now customize these basic maps and add information to them.

2 Customizing Maps

Displaying basic maps is usually less interesting because, typically, we want to add different layers to a map. In order to add layers to a map, we need to have a shape file, i.e. a file which contains information about borders or locations that can then be displayed in different colors. In other words, we need to have a shape object to add information to the map.

# extract locations
world_points<- st_centroid(world)
# extract labels
world_points <- cbind(world, st_coordinates(st_centroid(world$geometry)))
# generate annotated world map
ggplot(data = world) +
  geom_sf(fill= "gray90") +
  labs( x = "Longitude", y = "Latitude") +
  coord_sf(xlim = c(100.00, 180.00), ylim = c(-45.00, -10.00), expand = FALSE) +
  annotation_scale(location = "bl", width_hint = 0.5) +
  annotation_north_arrow(location = "bl", which_north = "true", 
                         pad_x = unit(0.75, "in"), pad_y = unit(0.5, "in"),
                         style = north_arrow_fancy_orienteering) +
  coord_sf(xlim = c(100.00, 180.00), ylim = c(-45.00, -10.00)) +
   theme(panel.grid.major = element_line(color = "gray60", linetype = "dashed", size = 0.25), 
         panel.background = element_rect(fill = "aliceblue")) +
  geom_text(data= world_points,aes(x=X, y=Y, label=name),
            color = "gray20", fontface = "italic", check_overlap = T, size = 3)

However, it is often the case that we want to add information that is not already available. Therefore, we load the airports data set which contains the longitude and latitude of airports across the world. We will then use this data to show the locations of airports across the globe.

# load data
airports <- read.delim("https://slcladal.github.io/data/airports.txt", 
                       sep = "\t", header = T)
# inspect data
datatable(airports, rownames = FALSE, options = list(pageLength = 5, scrollX=T), filter = "none")

To display the locations of airports on a map, we first plot the map and then add a layer of points to indicate the location of airports. In addition, the “plot” functions offers various arguments for customizing the display, e.g. by changing the backgroundcolor (bg), defining the color of borders (borders), defining the color of the shapes (fill and col).

# plot data on world map
plot(worldmap, xlim = c(-80, 160), ylim = c(-50, 100), 
     asp = 1, bg = "lightblue", col = "black", fill = T)
# add points
points(airports$Longitude, airports$Latitude, 
       col = "red", cex = .01)

It is, of course, also possible to highlight individual countries.

# create data frame with iso3 country codes and number of visits
countriesvisited <- data.frame(country = c("AUS", "JPN", "FIN", 
                                           "CZE", "POL", "AUT", 
                                           "USA", "GBR", "IRL", 
                                           "DEU", "DNK", "FRA", 
                                           "NDL", "BEL", "ESP",
                                           "HRV", "SVN", "NOR", 
                                           "ITA", "HUN", "ROU", 
                                           "BGR", "GRC", "TUR", 
                                           "CHE", "ARE"),
  visited = c(5, 1, 2, 1, 1, 3, 4, 4, 5, 11, 1, 1, 2, 2, 4, 4, 
              1, 1, 3, 1, 1, 2, 1, 1, 3, 2))
# inspect data
datatable(countriesvisited, rownames = FALSE, options = list(pageLength = 5, scrollX=T), filter = "none")
# combine data frame with map
visitedMap <- joinCountryData2Map(countriesvisited, 
                                  joinCode = "ISO3",
                                  nameJoinColumn = "country")
## 25 codes from your data successfully matched countries in the map
## 1 codes from your data failed to match with a country code in the map
## 218 codes from the map weren't represented in your data
# def. map parameters, e.g. def. colors
mapParams <- mapCountryData(visitedMap, 
                            nameColumnToPlot="visited",
                            oceanCol = "azure2",
                            catMethod = "categorical",
                            missingCountryCol = gray(.8),
                            colourPalette = c("coral",
                                              "coral2",
                                              "coral3", "orangered", 
                                              "orangered3", "orangered4"),
                            addLegend = F,
                            mapTitle = "",
                            border = NA)
# add legend and display map
do.call(addMapLegendBoxes, c(mapParams,
                             x = 'bottom',
                             title = "No. of visits",
                             horiz = TRUE,
                             bg = "transparent",
                             bty = "n"))

It is, of course also possible to show only a part of the map by defining the x- and y-axes limits of the plot window.

# get map
newmap <- getMap(resolution = "low")
# plot map
plot(newmap, xlim = c(-20, 59), ylim = c(35, 71), 
     asp = 1, fill = T, border = "darkgray", 
     col = "wheat2", bg = "gray95")
# add points
points(airports$Longitude, airports$Latitude, col = "red", cex = .5, pch = 20)

This is of course also possible to show Australian airports.

# plot data on world map
plot(worldmap, xlim = c(110, 160), ylim = c(-45, -10), 
     asp = 1, bg = "azure2", border = "lightgrey", col = "wheat1", 
     fill = T, wrap=c(-180,180))
points(airports$Longitude, airports$Latitude, 
       col = "darkblue", cex = .5, pch = 20)

In addition to the location of airports, it is also possible to show how many flights arrive at an airport. As this information is not provided in the airport data, we load the routes data which contains information about the routes that airlines fly.

# read in routes data
routes <- read.delim("https://slcladal.github.io/data/routes.txt", 
                     sep = "\t", header=T)
# inspect data
datatable(routes, rownames = FALSE, options = list(pageLength = 5, scrollX=T), filter = "none")