Links to the other two parts of the workshop:
ggplot2 notes: https://verticalmeadows.github.io/ggplot_basics.html
Introduction to R Markdown: https://verticalmeadows.github.io/basics.html

1 Map basics

Maps in R are, in essence, two dimensional plots where the x and y axes correspond to geographic longitude and latitude.

Several packages can be used to produce maps in R.
tmap specialises in thematic maps.
leaflet is the porting of a Javascript library, used for making interactive maps.
And ggmap is the one most useful for people familiar with ggplot2.

1.1 Mapping with ggplot

There is an easy way to make maps also with ggplot2, with the help of the map_data() function. This turns data from the maps package into a data frame suitable for plotting with ggplot2.
You can check if you need to install a package by the command require("maps") and if it is already installed, the command imports the package, acting like library().

Let’s import and plot New Zealand, for instance.

# https://r4ds.had.co.nz/data-visualisation.html#coordinate-systems
 nz <- map_data("nz")
 
ggplot(nz, aes(long, lat, group = group)) +
  geom_polygon(fill = "white", colour = "black")

 ggplot(nz, aes(long, lat, group = group)) +
   geom_polygon(fill = "white", colour = "black") +
   coord_quickmap()

Notice the difference the command coord_quickmap makes. On the second map NZ looks much more usual.

coord_map can take on more parameters, such as the projection, which would play a much more important role if you have a world map. But you might see the difference already in the latitude/longitude grid around new Zealand when zoomed out.
The Mercator projection projects the poles as a line, making longitudinal and latitudinal gridlines straight.
(This projection is used in Google Maps on your phone, if you zoom out all the way)

The aitoff projection is a modified azimuthal map projection.

world <- map_data("world")

 ggplot(nz, aes(long, lat, group = group)) +
   geom_polygon(fill = "white", colour = "black") +
   coord_map(projection="mercator", xlim=c(150,180), ylim=c(-60,-30))

 ggplot(world, aes(long, lat, group = group)) +
   geom_polygon(fill = "white", colour = "black") +
   coord_map(projection="mercator", xlim = c(-181,181))

 ggplot(nz, aes(long, lat, group = group)) +
   geom_polygon(fill = "white", colour = "black") +
   coord_map(projection="aitoff", xlim=c(150,180), ylim=c(-60,-30))

 ggplot(world, aes(long, lat, group = group)) +
   geom_polygon(fill = "white", colour = "black") +
   coord_map(projection="aitoff")

2 Ggmap

What is ggmap?
ggmap is an R package that makes it easy to retrieve raster map tiles from popular online mapping services like Google Maps and Stamen Maps and plot them using the ggplot2 framework.
A non-profit project of David Kahle: https://github.com/dkahle/ggmap

Essentially, it helps us source a background map easily.

Why ggmap?
ggmap is the fastest and most intuitive way to get mapping done for small areas with proper background. Plus as it is part of the gg ecosystem, one could expect it to work similarly to ggplot2.

ggmap used to work nicely with google maps until Google has limited the amount of time one can freely show their maps, one would need to register a credit card with google to get an API. A simple guide included here for that: https://www.littlemissdata.com/blog/maps
So instead we’ll use the provider Stamen for background maps.
Check the options: http://maps.stamen.com/#terrain/10/46.9362/7.7309

2.1 Installation

Its installation has not always been straightforward but now there seems to be a stable version on CRAN: ìnstall.packages("ggmap")

Otherwise, it has to be pulled from Github:

if (!requireNamespace("devtools")) install.packages("devtools") # for an installation from Github, one needs the package 'devtools'
devtools::install_github("dkahle/ggmap")

2.2 Basics

library(tidyverse)
library(ggmap)

myBoundingBox <- c(6.9, 45.9, 10, 47.85) # the coordinates of the encompassing rectangle of your map: LowerLeftLongitude, LowerLeftLatitude, UpperRightLongitude, UpperLeftLatitude
myMap <- get_stamenmap(bbox=myBoundingBox, maptype="terrain-background", crop=TRUE, color = 'color')

# Don't try it with googlemaps before setting it up
#myMap2 <- get_googlemap(location=myBoundingBox, source="google", zoom=10, maptype="roadmap", crop=TRUE)

myBasemap<-ggmap(myMap, extent = "device", darken = c(0.3, "white")) # if we use `device` as an extent, removes latitude, longitude and their ticks, darken lightens it if we use white as a colour

Loading the map background looks like the image below. Essentially we download raster tiles from the provider, which ggmap puts together as our background map. By cropping it, our background map will correspond to the bounding box supplied, otherwise the entire tile will be shown. Notice 10 in the file names - it means the zoom level. As the tiles are stored as rasters, there is no seamless zooming possible.
(If you often use google maps, you may also notice that when you zoom in, it happens by “jumps” in different levels of details. it means that there are different sets of tiles, depending on the zoom level (1-16).)
One can set this zoom level when downloading with get_stamenmap().
Also, one can set the tiles black and white: color = 'bw'

#display the basemap
myBasemap

3 Different layers on a map

Go through the real life example of mapping perception of dialects.

Let’s bring in some Swiss data, about 550 German-speaking Swiss communities, along with their population in 2000 and 2018

spatialdata <-  read.csv("./Spatial metadata0.csv", header=T, stringsAsFactors = FALSE) # numbers shouldnt be treated as factors.
spatialdata$Comment_AltNR <- NULL
spatialdata[spatialdata=="#N/A"] <- NA # supply NA's that R considers NA's as well
spatialdata$POPS2018 <- as.numeric(spatialdata$POPS2018)
spatialdata$POPS2000 <- as.numeric(spatialdata$POPS2000)

3.1 Points

Let’s map the points in the dataset, and colour them according to their canton, using the palette created in the ggplot part of this workshop.
(This dataset has the half-cantons merged.)

myBasemap +
  geom_point(data=spatialdata, 
             aes(x = LONG, y = LAT,color=SDS_KT),
             size = 4) + 
  scale_color_manual(values=myPalette) +
  labs(x= "Longitude", y = "Latitude")

3.2 Lines

We want to see some boundaries too. Let’s import linear feature from so called shapefiles that are often used in Geographic Information Science. (Another way to import spatial data is through [geo]jsons. Check out )

library(rgdal) # needed for the related commands
schweiz_whole <- readOGR(dsn="shapefiles/", layer="Switzerland_wgs84") # layer corresponds to the filename of the shp file
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/user/Work/Prog/R Visu Workshop/Workshop/shapefiles", layer: "Switzerland_wgs84"
## with 1 features
## It has 4 fields
## Integer64 fields read as strings:  KANTONSNR
schweiz_DE_WGS <- readOGR(dsn="shapefiles/", layer="SADS_DE_wgs84")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/user/Work/Prog/R Visu Workshop/Workshop/shapefiles", layer: "SADS_DE_wgs84"
## with 2 features
## It has 4 fields
## Integer64 fields read as strings:  KANTONSNR
df_schweiz <- fortify(schweiz_whole) # fortify converts the shapefile into a  dataset mappable for ggmap (a dataframe)
df_schweiz_DE <- fortify(schweiz_DE_WGS)

# Mapping them with ggplot, to see how they look
ggplot() +  
  geom_map(data=df_schweiz_DE, map=df_schweiz_DE,
                     aes(map_id=id, x=long, y=lat, group=group),
                     color="#ffffff", fill="#bbbbbb", size=0.25) +
  geom_path(data=df_schweiz,aes(long,lat,group=group),color="black")

myBasemap +
  geom_point(data=spatialdata, 
             aes(x = LONG, y = LAT,color=SDS_KT),
             size = 4) + 
  scale_color_manual(values=myPalette)  +
  geom_path(data=df_schweiz_DE,aes(long,lat,group=group),color="red", linetype=2) +
  geom_path(data=df_schweiz,aes(long,lat,group=group),color="black")

Then we might want to see the communities based on on their population

library(viridis) # for nice colour scales
# first let's remove 0s, distorting our size scale
spatialdata <- spatialdata %>%
  mutate(POPS2018 = na_if(POPS2018, 0))

myBasemap +
  geom_path(data=df_schweiz_DE,aes(long,lat,group=group),color="red", linetype=2) +
  geom_path(data=df_schweiz,aes(long,lat,group=group),color="black") +
  geom_point(data=spatialdata, 
             aes(x = LONG, y = LAT,size=POPS2018, color=POPS2018), alpha = 0.5) + 
  scale_radius(range = c(1,20), trans = "sqrt") + # behold this scaling function for ggplot. let's transform the population values such that the big cities' rule seems not so overwhelming, by putting it on an exponential scale
  scale_color_viridis(option="B")

We can also label them and have the labels as non-overlapping as possible using ggrepel.
Label the subset of towns larger than 20’000 inhabitants.

library(ggrepel)
# first let's remove 0s, distorting our size scale
spatialdata <- spatialdata %>%
  mutate(POPS2018 = na_if(POPS2018, 0))

myBasemap +
  geom_path(data=df_schweiz_DE,aes(long,lat,group=group),color="red", linetype=2) +
  geom_path(data=df_schweiz,aes(long,lat,group=group),color="black") +
  geom_point(data=spatialdata, 
             aes(x = LONG, y = LAT,size=POPS2018, color=POPS2018), alpha = 0.5) + 
  scale_radius(range = c(1,20), trans = "sqrt") + 
  scale_color_viridis(option="B") +
  geom_label_repel(data=spatialdata[which(spatialdata$POPS2018>20000),], aes(label=SDS_NAME, x=LONG, y=LAT))

3.3 Creating Voronoi-polygons

(a.k.a. Thiessen-polygons a.k.a. Dirichlet-tessellation)
A method often used in cartography for getting an extra visual variable is to create abutting polygons around points. In dialectology it is often used for showing an assumed area of influence of a certain survey site.

library(ggvoronoi)

inhab <- spatialdata %>%
  filter(!(is.na(POPS2018) | is.na(POPS2000))) %>% # only keep the rows where data is available for both years
  mutate(proportion = POPS2018/POPS2000) %>% # create a new coloumn which we map
  filter(proportion<10) # remove the outliers

myBasemap +
  geom_path(data=df_schweiz_DE,aes(long,lat,group=group),color="red", linetype=2) +
  geom_voronoi(data=inhab,aes(x = LONG, y = LAT), outline = df_schweiz_DE, color="white") + 
  geom_path(data=df_schweiz,aes(long,lat,group=group),color="black") +
  geom_point(data=spatialdata, 
             aes(x = LONG, y = LAT,size=POPS2018),color = "green", alpha = 0.5) + 
  scale_radius(range = c(1,20), trans = "sqrt") + 
  geom_label_repel(data=spatialdata[which(spatialdata$POPS2018>20000),], aes(label=SDS_NAME, x=LONG, y=LAT))

Note that we only draw a polygon around those sites where the population for both times is known.

Then, let’s use the polygons to plot the proportion by which population increased in between 2000 and 2018.

myBasemap +
  geom_path(data=df_schweiz_DE,aes(long,lat,group=group),color="red", linetype=2) +
  geom_voronoi(data=inhab,aes(x = LONG, y = LAT,fill=POPS2018/POPS2000), outline = df_schweiz_DE) + 
  scale_fill_viridis(option = "inferno") + # this could be mapped as categories, same goes with the inhabitants
  geom_path(data=df_schweiz,aes(long,lat,group=group),color="black") +
  geom_point(data=spatialdata, 
             aes(x = LONG, y = LAT,size=POPS2018),color = "green", alpha = 0.5) + 
  scale_radius(range = c(1,20), trans = "sqrt") + 
  geom_label_repel(data=spatialdata[which(spatialdata$POPS2018>20000),], aes(label=SDS_NAME, x=LONG, y=LAT))

Density plotting is also possible.

myBasemap +
   geom_point(data=spatialdata, 
             aes(x = LONG, y = LAT,color=SDS_KT),
              size = 4)+ 
  scale_color_manual(values=myPalette) +
  geom_density2d(data=spatialdata, 
              aes(x = LONG, y = LAT), size = 0.3, geom = "polygon")

4 Further reading, by topics

4.1 Maps in R

From Robin Lovelace, Jakub Nowosad & Jannes Muenchow, ‘Geocomputation with R’: https://geocompr.robinlovelace.net/adv-map.html#interactive-maps