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.
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.
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
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()
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.
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
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.
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")