Multi Faceted Maps in R: using GGPLOT2 and TMAP

Premise

During my PhD research, I had to create subsets of thematic maps depicting Land Use and Land Cover (LULC) classifications using various machine learning algorithms. Alongside these thematic maps, I aimed to integrate two base maps: Sentinel-2A by ESA (European Space Agency) and NAIP (National Agriculture Imagery Program).

One of the challenges I encountered was the integration of graticule labels on axes (both X and Y) resulted wide gaps and misalignment. I accomplished this first generating several individual maps and then consolidate them using the patchwork package. Before, patchwork I used to rely on cowplot and gridExtra. Avoiding graticule labels, I was close to align these maps (thematic and base maps). however, problem still persisted while presenting figures small polygons were disproportionately highlighted. I will probably write on the latter issue in next blog, but for now. lets focus on combining several maps. My approach about three years ago resulted like this:

Faceted thematic and base maps

Recently, I also came across similar issue of making several maps of cities in the United States. To highlight the issue. I downloaded a city boundaries data of Texas from Texas Department of Transportation (TxDOT)

The Data

library(sf)        # manipulate gepsiatial data
library(magrittr)  # pipe operator
library(dplyr)     # manipulate data
library(ggplot2)   # plot geospatial data geom_sf
library(tmap)      # map spatial data
library(ggspatial) # annotation and arrow
library(patchwork) # plot multiple ggplots

# for simplicity I will be selecting top six cities 

## read data 
dat <- sf::st_read("./TxDOT_City_Boundaries.gpkg") %>% 
        dplyr::filter(CITY_NM %in% c("Corpus Christi", 
        "San Antonio", "Austin", "Dallas", "Houston","Fort Worth"))

Reading layer `Cities' from data source 
  `C:\NRM5404\00_Multiple_maps\TxDOT_City_Boundaries.gpkg' using driver `GPKG'
Simple feature collection with 1223 features and 13 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -11870670 ymin: 2979310 xmax: -10430690 ymax: 4369628
Projected CRS: WGS 84 / Pseudo-Mercator

#-- dat %>% mutate(area = st_area(.)) %>% arrange(desc(area)) %>% head(6)

Examine Data

For simplicity I did not use city council or other features that falls within the city limits

str(dat)

Classes 'sf' and 'data.frame':  6 obs. of  14 variables:
 $ CITY_NM       : chr  "Austin" "Corpus Christi" "Fort Worth" "San Antonio" ...
 $ TXDOT_CITY_NBR: int  2100 9800 15000 37450 19750 10850
 $ CITY_FIPS     : chr  "4805000" "4817000" "4827000" "4865000" ...
 $ CNTY_SEAT_FLAG: chr  "Y" "Y" "Y" "Y" ...
 $ POP1990       : int  465622 257453 447619 935933 1630553 1006977
 $ POP2000       : int  656562 277454 534694 1144646 1953631 1885580
 $ POP2010       : int  790390 305215 741206 1327407 2099451 1197816
 $ POP2020       : int  995484 327248 927720 1567118 2316120 1343266
 $ POP2022       : int  974447 316239 956709 1472909 2302878 1299544
 $ POP_CD        : chr  "4" "4" "4" "4" ...
 $ MAP_COLOR_CD  : int  1 1 1 1 1 1
 $ COLOR_CD      : int  1 1 1 1 1 1
 $ GID           : num  1210 975 1002 1130 1175 ...
 $ SHAPE         :sfc_MULTIPOLYGON of length 6; first list element: List of 3
  ..$ :List of 34
  .. ..$ : num [1:16044, 1:2] -10877896 -10877868 -10877833 -10877792 -10877735 ...
  .. ..$ : num [1:30, 1:2] -10886198 -10886195 -10886190 -10886193 -10886194 ...
  .. ..$ : num [1:17, 1:2] -10866855 -10866838 -10866837 -10866826 -10866790 ...
  .. ..$ : num [1:32, 1:2] -10881576 -10881563 -10881561 -10881560 -10881559 ...
  .. ..$ : num [1:16, 1:2] -10884485 -10884483 -10884477 -10884476 -10884465 ...
  .. ..$ : num [1:39, 1:2] -10887624 -10887624 -10887621 -10887617 -10887590 ...
  .. ..$ : num [1:30, 1:2] -10885064 -10885078 -10885079 -10885094 -10885101 ...
  .. ..$ : num [1:12, 1:2] -10887679 -10887816 -10887847 -10887846 -10887800 ...
  .. ..$ : num [1:7, 1:2] -10886998 -10886938 -10886844 -10886692 -10886763 ...
  .. ..$ : num [1:30, 1:2] -10866947 -10866955 -10866964 -10866972 -10866981 ...
  .. ..$ : num [1:18, 1:2] -10881734 -10881701 -10881692 -10881659 -10881512 ...
  .. ..$ : num [1:51, 1:2] -10885226 -10885244 -10885294 -10885324 -10885336 ...
  .. ..$ : num [1:34, 1:2] -10884498 -10884507 -10884510 -10884516 -10884522 ...
  .. ..$ : num [1:42, 1:2] -10887845 -10887844 -10887755 -10887728 -10887688 ...
  .. ..$ : num [1:49, 1:2] -10888470 -10888449 -10888455 -10888457 -10888465 ...
  .. ..$ : num [1:56, 1:2] -10887659 -10887672 -10887668 -10887495 -10887318 ...
  .. ..$ : num [1:70, 1:2] -10890682 -10890685 -10890685 -10890639 -10890639 ...
  .. ..$ : num [1:45, 1:2] -10883721 -10883709 -10883701 -10883697 -10883692 ...
  .. ..$ : num [1:161, 1:2] -10872391 -10872389 -10872388 -10872387 -10872385 ...
  .. ..$ : num [1:202, 1:2] -10870313 -10870314 -10870314 -10870315 -10870318 ...
  .. ..$ : num [1:45, 1:2] -10891973 -10891978 -10891981 -10892003 -10892005 ...
  .. ..$ : num [1:79, 1:2] -10865783 -10865779 -10865872 -10865965 -10865961 ...
  .. ..$ : num [1:32, 1:2] -10870530 -10870496 -10870435 -10870315 -10870197 ...
  .. ..$ : num [1:211, 1:2] -10888400 -10888364 -10888401 -10888402 -10888453 ...
  .. ..$ : num [1:141, 1:2] -10889803 -10889814 -10889820 -10889840 -10889852 ...
  .. ..$ : num [1:17, 1:2] -10863511 -10863544 -10863629 -10863686 -10863871 ...
  .. ..$ : num [1:197, 1:2] -10871705 -10871716 -10871723 -10871733 -10871756 ...
  .. ..$ : num [1:332, 1:2] -10889170 -10889168 -10889164 -10889215 -10889234 ...
  .. ..$ : num [1:608, 1:2] -10884623 -10884571 -10884443 -10884216 -10883799 ...
  .. ..$ : num [1:412, 1:2] -10892637 -10892633 -10892660 -10892664 -10892672 ...
  .. ..$ : num [1:587, 1:2] -10867261 -10867261 -10867264 -10867267 -10867537 ...
  .. ..$ : num [1:1079, 1:2] -10885428 -10885423 -10885384 -10885449 -10885502 ...
  .. ..$ : num [1:2733, 1:2] -10873844 -10873721 -10873720 -10873718 -10873716 ...
  .. ..$ : num [1:776, 1:2] -10888335 -10888353 -10888390 -10888397 -10888410 ...
  ..$ :List of 1
  .. ..$ : num [1:28, 1:2] -10867902 -10867859 -10867849 -10867804 -10867598 ...
  ..$ :List of 1
  .. ..$ : num [1:192, 1:2] -10867219 -10867499 -10867499 -10867499 -10867499 ...
  ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
 - attr(*, "sf_column")= chr "SHAPE"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
  ..- attr(*, "names")= chr [1:13] "CITY_NM" "TXDOT_CITY_NBR" "CITY_FIPS" "CNTY_SEAT_FLAG" ...

head(dat[, 1:5])

Simple feature collection with 6 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -10999860 ymin: 3184887 xmax: -10577180 ymax: 3901848
Projected CRS: WGS 84 / Pseudo-Mercator
         CITY_NM TXDOT_CITY_NBR CITY_FIPS CNTY_SEAT_FLAG POP1990
1         Austin           2100   4805000              Y  465622
2 Corpus Christi           9800   4817000              Y  257453
3     Fort Worth          15000   4827000              Y  447619
4    San Antonio          37450   4865000              Y  935933
5        Houston          19750   4835000              Y 1630553
6         Dallas          10850   4819000              N 1006977
                           SHAPE
1 MULTIPOLYGON (((-10877896 3...
2 MULTIPOLYGON (((-10869713 3...
3 MULTIPOLYGON (((-10860084 3...
4 MULTIPOLYGON (((-10941677 3...
5 MULTIPOLYGON (((-10602764 3...
6 MULTIPOLYGON (((-10757313 3...

Take I: Using GGplot

Let’s plot using ggplot2

dat %>% ggplot()+geom_sf(aes(fill = POP2022))

Take I - faceting

dat %>% ggplot()+geom_sf(aes(fill = POP2022))+
  facet_wrap(vars(CITY_NM))+
  scale_fill_continuous(name = "Population\nCensus (2022)")

with geom_sf facet_wrap can’t be use with scales = 'free'. In such case we can try tmap package. Note that tmap3.x is retiring. Version 4 may already be able to solve such problem we can see below.

TakeI: using tmap

  dat %>% 
  tm_shape() +
  tm_polygons(col = "POP2022",
              title = "Population\nCensus (2022)",
              style = "cont") +                  
  tm_facets("CITY_NM", ncol = 3) +
  tm_layout(legend.outside.position =  "right")        

The map is much better already. Now lets add graticule, and control labelling along with scale and orientation.

  dat %>% 
  tm_shape() +
  tm_graticules(n.x = 3, n.y = 3)+
    tm_xlab("Longitude")+
  tm_ylab("Latitude")+
  tm_polygons(col = "POP2022",
              title = "Population\nCensus (2022)",
              style = "cont") +   
  tm_facets("CITY_NM", ncol = 3) +
  # tm_xlab("Longitude")+
  # tm_ylab("Latitude")+
  tm_scale_bar(width = .15, text.size = 1, text.color = "orangered4", lwd = 2)

  tm_layout(legend.outside.position =  "right")   

$tm_layout
$tm_layout$legend.outside.position
[1] "right"

$tm_layout$style
[1] NA


attr(,"class")
[1] "tm"

Take I - tmap Final

For publication quality maps lets modify several stuffs

  dat %>% 
  tm_shape() +
  tm_graticules(n.x = 3, n.y = 3)+
  tm_polygons(col = "POP2022",
              title = "Population\nCensus (2022)",
              style = "cont") +   
  tm_facets("CITY_NM", ncol = 2) +
  tm_borders(alpha = 0.5) + # tm_scale_bar width or breaks
  tm_scale_bar(position = c("left", "bottom"), 
               breaks = c(0, 5, 10, 15), text.size = 1) +
  tm_compass(position = c("left", "top"), size = 2)+
  tm_layout(fontface = "bold",
            title.size = 12,
            panel.label.size = 2,
            inner.margins = 0,
            legend.text.size = 1,
            legend.outside = TRUE,
            legend.outside.position = "right",
            #legend.outside.size = 2,
            panel.label.height = 1.5,
            legend.width = 2.5,
            legend.height = 5)

Warning in pre_process_gt(x, interactive = interactive, orig_crs =
gm$shape.orig_crs): legend.width controls the width of the legend within a map.
Please use legend.outside.size to control the width of the outside legend

Warning: One tm layer group has duplicated layer types, which are omitted. To
draw multiple layers of the same type, use multiple layer groups (i.e. specify
tm_shape prior to each of them).

Here we are so close but legend outside of layout area is difficult to control as we don’t have many options. For final map we can use same code to plot publication quality map; however, I would like to increase legend height.

Take II: Using GGPLot

Two things I want to improve are: 1. x, and y axis label to each sub-figure 2. ability to control legend (combine legend)

For this for better naming, I want to rename CITY_NM into city, and pop2022 to population. Since, scale can’t be free, I want to write a function such that bounding box for each city is same. In my actual example, I tried to put only two labels in x (longitude), and y (latitude) axis. This labeling control is way easy in tmap. But in ggplot, it is slighlty difficult. I am creating one function for this.

# Calculate the bounding boxes for each city

all_data <- dat %>% rename(city = CITY_NM, 
                           population = POP2022,
                           geometry = SHAPE) %>% st_transform(.,crs = "epsg:4326")
# Create function for labeling
custom_label_transparent <- function(values, is_x_axis = TRUE) {
  labels <- rep("", length(values))  
  if (is_x_axis) {
    labels[c(1, length(values))] <- paste0(sprintf("%.2f", values[c(1, length(values))]), "° W")  
  } else {
    labels[c(1, length(values))] <- paste0(sprintf("%.2f", values[c(1, length(values))]), "° S") 
  }
  return(labels)
}
## create equivalent bounding box

bounding_boxes <- list()
for (city_name in unique(all_data$city)) {
  city_data <- all_data %>% filter(city == city_name)
  city_bbox <- sf::st_bbox(city_data$geometry)
  bounding_boxes[[city_name]] <- city_bbox
}

# calculate the maximum bounding box size among all cities
max_bbox_size <- c(
  max(sapply(bounding_boxes, function(bb) bb[3] - bb[1])),
  max(sapply(bounding_boxes, function(bb) bb[4] - bb[2]))
)

# calculate the required expansion for each city and adjust the bounding boxes
for (city_name in names(bounding_boxes)) {
  city_bbox <- bounding_boxes[[city_name]]
  x_expansion <- max_bbox_size[1] - (city_bbox[3] - city_bbox[1])
  y_expansion <- max_bbox_size[2] - (city_bbox[4] - city_bbox[2])
  city_bbox[1] <- city_bbox[1] - x_expansion / 2
  city_bbox[3] <- city_bbox[3] + x_expansion / 2
  city_bbox[2] <- city_bbox[2] - y_expansion / 2
  city_bbox[4] <- city_bbox[4] + y_expansion / 2
  bounding_boxes[[city_name]] <- city_bbox
}


# Create the plots and store them in a list
plot_list <- list()
for (city_name in names(bounding_boxes)) {
  city_data <- all_data %>% filter(city == city_name)
  
  city_plot <- ggplot(data = city_data) +
    geom_sf(aes(fill = population), show.legend = TRUE) +
    coord_sf(xlim = c(bounding_boxes[[city_name]][1], bounding_boxes[[city_name]][3]),
             ylim = c(bounding_boxes[[city_name]][2], bounding_boxes[[city_name]][4])) +
    scale_x_continuous(
      breaks = round(seq(bounding_boxes[[city_name]][1], bounding_boxes[[city_name]][3], 0.2), 2),
      labels = custom_label_transparent(round(seq(bounding_boxes[[city_name]][1], bounding_boxes[[city_name]][3], 0.2), 2), is_x_axis = TRUE)
    ) +
    scale_y_continuous(
      breaks = round(seq(bounding_boxes[[city_name]][2], bounding_boxes[[city_name]][4], 0.2), 2),
      labels = custom_label_transparent(round(seq(bounding_boxes[[city_name]][2], bounding_boxes[[city_name]][4], 0.2), 2), is_x_axis = FALSE)
    ) +
    scale_fill_gradientn(name = "population\n", limits = c(316239, 2302878 ), 
                         colors = viridis::viridis(n = 10, option = "D")) +
    annotation_scale(
      location = "bl", width_hint = 0.5,
      height = unit(0.25, "cm"),
      pad_x = unit(0.25, "cm"),
      pad_y = unit(0.25, "cm"),
      text_pad = unit(0.25, "cm"),
      text_cex = 0.7,
      tick_height = 0.6
    ) +
    annotation_north_arrow(
      location = "bl", which_north = "true", 
      pad_x = unit(0.75, "cm"), pad_y = unit(0.5, "cm"),
      style = north_arrow_fancy_orienteering
    ) +
    xlab(" ") +
    ylab(" ") +
    facet_wrap(vars(city))+
   # ggtitle(city_name) +
    theme_bw() +
    theme(
      plot.margin = unit(c(0.5, 0.5, 0.5, 0.5), "cm"),
      panel.grid.major = element_line(color = gray(.5), linetype = "dotted", linewidth = 0.4),
      panel.background = element_rect(fill = "grey80"),
      strip.text.x = element_text(size = 16, color = "black", face = "bold"),
      strip.background = element_rect(color = "black", fill = "aliceblue", size = 1, linetype = "solid"),
      legend.key.height = unit(1, "cm"),
      legend.key.width = unit(1, "cm"),
      axis.title.x = element_text(size = 14, face = "bold"),
      axis.title.y = element_text(size = 14, face = "bold"),
      axis.text.x = element_text(face = "bold", size = 12),
      axis.text.y = element_text(face = "bold", size = 12),
      legend.text = element_text(size = 12, face = "bold"),
      axis.ticks.x = element_line(linewidth = 1),
      axis.ticks.y = element_line(linewidth = 1),
      legend.position = c(0.88, 0.190),
      legend.background = element_rect(fill = "grey80"),
      plot.title = element_text(size = 18, face = "bold", hjust = 0.5)
    )
  
  plot_list[[city_name]] <- city_plot
}

Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
ℹ Please use the `linewidth` argument instead.

Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.

# Combine the plots using patchwork with fixed aspect ratio and larger titles
# Combine the plots using patchwork with fixed aspect ratio and larger titles
combined_plots <- wrap_plots(plot_list, ncol = 2) +
  plot_layout(guides = 'collect') &
  theme(
    aspect.ratio = 1,
    plot.margin = unit(c(1, 1, 1, 1), "cm"),
    plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
    legend.position = "right",
    legend.background = element_rect(fill = "white", color = NA),
    legend.key.size = unit(3, "cm"),
    legend.key.width = unit(3, "cm"),
    legend.key.height = unit(8, "cm"),
    legend.title = element_text(size = 14),
    legend.text = element_text(size = 14)
  )

# Save the combined plot with consistent dimensions
width <- 20  
height <-24
ggsave("combined_plots_aligned.png",
combined_plots, 
width = width, height = height,units = "in",dpi = 300)

Final map

Although I haven’t discussed in detail about shared legends and map sizes, the examples presented in previous sections (implicitly) demonstrate the use of accurate shared legends alongside aligned maps. If you have any further questions or would like more information on these topics, please feel free to reach out to me via email.