Mapping with ggplot2

We will use what we've already learned about ggplot2 to create maps in R using data from Del Ponte and Belachew (2020) and from http://GADM.org, a source for administrative map data.

First, load the libraries needed for this work.

library(dplyr) # Manipulate data in R
library(ggplot2) # Create great plots in R!
library(ggsn) # Add scale bars and north arrows to sf maps in R
library(here) # Simplify file paths in R
library(readr) # Import tabular data in R
library(sf) # sf stands for "Simple Features", spatial data in R

GADM is a free resource for research and non-commercial purposes. Read more at https://gadm.org/about. Here we'll set up our R session to download the file. Just for the purposes of the presentation, I've broken up the URL here so that we can fit everything on the presentation slide. Ordinarily there's no need to use paste0() to create a URL!

# Remote file information
u_remote <- "https://biogeo.ucdavis.edu/"
p_remote <- "data/gadm3.6/Rsf/"
f_name <- "gadm36_ETH_3_sf.rds"

GADM distributes native R files in .rds format that we then import.

ethiopia_sf <- readRDS(ethiopia_rds)

We can see that there are several zones, shown as the smaller polygons within the larger country outline, in this data set. We can find out which columns are available to use in our maps.

names(ethiopia_sf)

 [1] "GID_0"     "NAME_0"    "GID_1"     "NAME_1"    "NL_NAME_1" "GID_2"    
 [7] "NAME_2"    "NL_NAME_2" "GID_3"     "NAME_3"    "VARNAME_3" "NL_NAME_3"
[13] "TYPE_3"    "ENGTYPE_3" "CC_3"      "HASC_3"    "geometry" 

We can group the GADM data using dplyr::group_by() to group the polygons into larger regions and add labels to them. In this data set, there are three levels and three columns we might be interested in.

• NAME_0 refers to the country level and gives a whole country outline,
• NAME_1 refers to the second level of administrative boundaries, e.g. states within the USA,
• NAME_2 refers to the third level of administrative boundaries, e.g. counties or parishes within a US state, and
• NAME_3 refers to the fourth level of administrative boundaries, e.g. districts within counties in a US state.

ethiopia_regions_ggplot <- 
  ethiopia_regions_ggplot +
  labs(x = "Longitude",
       y = "Latitude") +
  theme_bw()

Now we will build on the maps we created by adding the survey location points that Del Ponte and Belachew (2020) used with latitude and longitude values from their CSV file.

coffee_survey <-
  read_csv(here("data", "survey_clean.csv"))

While we could plot everything using geom_point() directly from the .CSV file data.frame, but if we convert the data.frame to an sf object it is much more flexible for mapping with other sf objects in ggplot2.

coffee_survey_sf <-
  st_as_sf(coffee_survey,
           coords = c("lon", "lat"),
           crs = 4326)

Using the already existing ggplot2 object, ethiopia_regions_ggplot, we can add the new coffee_survey_sf object to our map!

ethiopia_regions_ggplot +
  geom_sf(data = coffee_survey_sf,
          size = 2,
          color = alpha("black", 0.35))

Warning message:
In st_point_on_surface.sfc(sf::st_zm(x)) :
  st_point_on_surface may not give correct results for longitude/latitude data

Ethiopia is very near the equator, the distortion is minimal and the approximation of a flat two-dimensional surface is reasonable, so the points will be accurately displayed and we can ignore the warning.

Because sf is a data.frame, it works well with dplyr so we can filter() the regions that we want to display.

We will create two layers. The first will be the outline of the country. The second will be the regions we wish to highlight.

# First layer, country outline
ethiopia_simple_sf <-
  ethiopia_sf %>%
  group_by(NAME_0) %>%
  summarise() %>%
  ungroup() %>%
  st_as_sf()

snnp_oromia_ggplot2 <- ggplot() +
  geom_sf(data = ethiopia_simple_sf,
          col = NA) +
  geom_sf(data = oromia_snnp_sf,
          aes(linetype = NAME_1,
              fill = NAME_1)) +
  scale_fill_brewer(palette = "Dark2") +
  labs(
    caption = "Data from GADM (gadm.org).",
    fill = "Region",
    linetype = "Region",
    x = "Longitude",
    y = "Latitude"
  ) +
  theme_bw()

Del Ponte and Belechew (2020) display a map that highlights the two regions, Oromia and SNNP with the zones within them where surveys were conducted also coloured.

final_map <- 
  final_map +
  scale_fill_brewer(palette = "Set3") # this colours the zones, https://colorbrewer2.org/

final_map <- 
  final_map +
  labs(
    x = "Longitude",
    y = "Latitude",
    title = "Ethiopian Coffee Leaf Rust Survey Sites",
    caption = "Data from GADM, https://gadm.org, and Del Ponte & Belachew, 2020,
    https://doi.org/10.17605/OSF.IO/XEJAZ",
    fill = "Zone",
    linetype = "Region"
  )

Lastly, add standard map items like a scale bar and north arrow using functions from ggsn, scalebar() and north().

final_map <-
  final_map +
  north(ethiopia_simple_sf) + # from ggsn
  scalebar( # from ggsn
    ethiopia_simple_sf,
    dist = 250,
    st.size = 2.5,
    dist_unit = "km",
    transform = TRUE,
    model = "WGS84"
  ) +
theme_bw()