In-Class Exercise 3: Analytical Mapping

Published

January 25, 2023

Modified

April 3, 2023

0.1 Overview

0.2 Installing and loading packages

pacman::p_load(tmap, tidyverse, sf)

The tmap package documentation can be found here.

0.3 Importing data

NGA_wp <- read_rds("data/rds/NGA_wp.rds")

0.4 Visualizing distribution of functional and non-functional water point using Choropleth Maps

p1 <- tm_shape(NGA_wp) + 
  tm_fill("wp_functional",
          n = 10,
          
          #style is used for data classification method
          style = "equal",
          
          #colour palettes are always plural
          palette = "Blues") +  
  
  # line width
  tm_borders(lwd = 0.1,
             
             # opacity
             alpha = 1) + 
  
  # main.title places the title outside the plot
  tm_layout(main.title = "Distribution of functional water point",
            
            # legend.outside puts legends in- or outside of your plot
            legend.outside = FALSE)

# things to take note: tm_fill() and tm_borders() combined is tm_polygon()

p1

p1 <- tm_shape(NGA_wp) + 
  tm_fill("wp_functional",
          n = 10,
          
          #style is used for data classification method
          style = "equal",
          
          #colour palettes are always plural
          palette = "Blues") +  
  
  # line width
  tm_borders(lwd = 0.1,
             
             # opacity
             alpha = 1) + 
  
  # main.title places the title outside the plot
  tm_layout(main.title = "Distribution of functional water points by LGAs",
            
            # legend.outside puts legends in- or outside of your plot
            legend.outside = FALSE)

# things to take note: tm_fill() and tm_borders() combined is tm_polygon()

# p1 is a map object
p1

p2 <- tm_shape(NGA_wp) + 
  tm_fill("total_wp",
          n = 10,
          
          #style is used for data classification method
          style = "equal",
          
          #colour palettes are always plural
          palette = "Blues") +  
  
  # line width
  tm_borders(lwd = 0.1,
             
             # opacity
             alpha = 1) + 
  
  # main.title places the title outside the plot
  tm_layout(main.title = "Distribution of total water points by LGAs",
            
            # legend.outside puts legends in- or outside of your plot
            legend.outside = FALSE)

p2

# show both together on one row
tmap_arrange(p2, p1, nrow=1)

# You can even show both maps as one by using Rates: usually we do this to see if both categories tell the same story
NGA_wp <- NGA_wp %>%
  
  # use mutate() to calc % (i.e., rate) of functional and non-functional water points
  mutate(pct_functional = wp_functional/total_wp) %>%
  mutate(pct_nonfunctional = wp_nonfunctional/total_wp)
tm_shape(NGA_wp) + 
  tm_fill("pct_functional",
          n = 10,
          style = "equal",
          palette = "Blues",
          legend.hist = TRUE) +
  tm_borders(lwd = 0.1,
             alpha = 1) +
  tm_layout(main.title = "Rate map of functional water points by LGAs",
            legend.outside = TRUE)

0.5 Visualization using a Percentile Map

Tells you which areas are the top 10%. The six specific categories are: 0-1%, 1-10%, 10-50%, 50-90%, 90-99%, and 99-100%.

# Data Preparation
  # Exclude NA values
NGA_wp <- NGA_wp %>%
  drop_na()

  # create customised classification and extract values
percent <- c(0, .01, .1, .5, .9, .99, 1)
var <- NGA_wp["pct_functional"] %>%
  
  # 
  st_set_geometry(NULL)
quantile(var[,1], percent)
       0%        1%       10%       50%       90%       99%      100% 
0.0000000 0.0000000 0.2169811 0.4791667 0.8611111 1.0000000 1.0000000 

0.5.1 Writing functions to do the same functionality for specific data sets…

# R function to extract a variable (i.e., wp_nonfunctional as a vector out of an sf data.frame)
get.var <- function(vname, df) {
  v <- df[vname] %>%
    st_set_geometry(NULL)
    v <- unname(v[,1])
    
  # return vector with values (without a col name)
  return(v)
}
# percentile mapping function
percentmap <- function(vnam, df, legtitle=NA, mtitle="Percentile Map"){
  percent <- c(0,.01,.1,.5,.9,.99,1)
  var <- get.var(vnam, df)
  bperc <- quantile(var, percent)
  tm_shape(df) +
  tm_polygons() +
  tm_shape(df) +
     tm_fill(vnam,
             title=legtitle,
             breaks=bperc,
             palette="Blues",
          labels=c("< 1%", "1% - 10%", "10% - 50%", "50% - 90%", "90% - 99%", "> 99%"))  +
  tm_borders() +
  tm_layout(main.title = mtitle, 
            title.position = c("right","bottom"))
}
# test the function
percentmap("total_wp", NGA_wp)

0.6 Visualizing using a Box Plot

# boxplot is as augmented quartile map with an additional lower and upper category.
ggplot(data = NGA_wp,
       aes(x = "",
           y = wp_nonfunctional)) + 
  geom_boxplot()

# creating a boxbreaks function
# arguments - v: vector with observations, mult: multiplier for IQR (default 1.5)
# returns - bb: vector with 7 breakpoints compute quartile and fences

boxbreaks <- function(v,mult=1.5) {
  qv <- unname(quantile(v))
  iqr <- qv[4] - qv[2]
  upfence <- qv[4] + mult * iqr
  lofence <- qv[2] - mult * iqr
  # initialize break points vector
  bb <- vector(mode="numeric",length=7)
  # logic for lower and upper fences
  if (lofence < qv[1]) {  # no lower outliers
    bb[1] <- lofence
    bb[2] <- floor(qv[1])
  } else {
    bb[2] <- lofence
    bb[1] <- qv[1]
  }
  if (upfence > qv[5]) { # no upper outliers
    bb[7] <- upfence
    bb[6] <- ceiling(qv[5])
  } else {
    bb[6] <- upfence
    bb[7] <- qv[5]
  }
  bb[3:5] <- qv[2:4]
  return(bb)
}
# creating the get.var function
# arguments - vname: variable name (in quotes), df: name of sf data.frame
# returns - v: vector with values (without col name)
get.var <- function(vname,df) {
  v <- df[vname] %>% st_set_geometry(NULL)
  v <- unname(v[,1])
  return(v)
}
# test function
var <- get.var("wp_nonfunctional", NGA_wp) 
boxbreaks(var)
[1] -56.5   0.0  14.0  34.0  61.0 131.5 278.0

0.6.1 Boxmap function

boxmap <- function(vnam, df, 
                   legtitle=NA,
                   mtitle="Box Map",
                   mult=1.5){
  var <- get.var(vnam,df)
  bb <- boxbreaks(var)
  tm_shape(df) +
    tm_polygons() +
  tm_shape(df) +
     tm_fill(vnam,title=legtitle,
             breaks=bb,
             palette="Set3",
          labels = c("lower outlier", 
                     "< 25%", 
                     "25% - 50%", 
                     "50% - 75%",
                     "> 75%", 
                     "upper outlier"))  +
  tm_borders() +
  tm_layout(main.title = mtitle, 
            title.position = c("left",
                               "top"))
}
tmap_mode("plot")
boxmap("wp_nonfunctional", NGA_wp)

0.7 Recode to zero

This code chunk is used to recode LGAs with zero total water points into NA.

NGA_wp <- NGA_wp %>%
  mutate(wp_functional = na_if(
    total_wp, total_wp < 0))