Automating map-making with tmap

Image credit: Nils Riach

Speed up repetitive map making with the tmap-package by iterating through multiple files with for-loops

Climate Models

Future climatic conditions on the regional scale are commonly modelled with different combinations of global circulation models (GCM) and regional climate model (RCM) in order to account for uncertainties surrounding the projection of future conditions1. A collection of models is called a “model ensemble” and often encompasses 20 or more different GCM-RCM combinations. Naturally, these model ensembles produce a lot of data which makes the (semi-) automated analysis and visualization necessary.

Data Structure

Recently, I found myself analyzing such a model ensemble, basically a bunch of shapefiles all equal in terms of geometric extent and data structure (variables and observations). The only difference between them is the file name and, of course, the values inside the columns.

I wanted to do a “quick-and-dirty” vizualization in order to get a better understanding of my data. Visualizing each file individually would have been too time-consuming. Unfortunately, I only found one example of [automated map creation (with ggplot), so I wrote a foor loop to create maps with the tmap-package.

MWE

Below is a reproducible MWE of how a multitude of shapefiles can be mapped automatically. I boiled it down to the key functions, however, the full functionality of the tmap package (and others) may be explored. I include only 4 layers but the example will also work with more layers. They just need to be included in the folder (and of course the data structure needs to be the same as in the two already existing layers).

load packages

# Dependings
if (!require(sf)) {install.packages("sf"); library(sf)}
if (!require(dplyr)) {install.packages("dplyr"); library(dplyr)}
if (!require(readr)) {install.packages("readr"); library(readr)}
if (!require(tmap)) {install.packages("tmap"); library(tmap)}

get data

First We need the data. They can be downloaded from GitHub. Save them to a folder of your choice. Also, define a directory where your results shall be saved.

# clear global environment
rm(list=ls())

#set wd. This is where the data are stored
setwd("V:/blog")
getwd()

# directory where results are saved
# folder needs to be created 
result <- "V:/blog/results/"

execute the loop

# get list of shp files from working directory
shape_list2=list.files(getwd(), pattern="*shp", full.names=)
#start for loop
for(shapefile in 1:length(shape_list2)){
  assign(unlist(strsplit(shape_list2[shapefile], "[.]"))[1], st_read(shape_list2[shapefile])) 
  # a is a placeholder for the different shapefiles. read them one by one
  a <- st_read(shape_list2[shapefile])
  # remove ".shp of name_string for name
  shape_list2_name <-stringr::str_remove_all(shape_list2, ".shp")
  # remove ending of file name. Use remaining parts for pdf file name and title in map
  shape_list3_name <- stringr::str_remove_all(shape_list2, "_community_mean.shp")
  # map  layer
  # include message to keep track of process
  message("mapping layer", shapefile, ' of ', length(shape_list2))
  # create pdf. file name based on  result + shape_list3_name
  pdf(file = paste(result, unlist(strsplit(shape_list3_name[shapefile], "[.]")),"_ens_50.pdf",
                   sep=""))
  # print is important to tell r to create 
  # from here it is all tmap-functionality
  ## col = "ens_50" refers to column which is used for visualization 
  ## style defines method to process the color scale
  print(tm_shape(a) + tm_fill(col = "ens_50", palette = "Blues", style = "jenks", 
                              title = "Ensemble Percentile 50",
                              legend.hist = TRUE)+
          #add borders 
          tm_borders(alpha= .2, col = "black") + 
          # add compass, scalebar and credits
          tm_compass() +
          tm_scale_bar() +
          tm_credits("Data:  credits", size = 0.7, position=c("left","bottom"))+
          # add title to map based on string in shape_list3_name
          tm_layout(title = paste(unlist(strsplit(shape_list3_name[shapefile], "[.]")),
                                  sep=""),  legend.outside= FALSE,
                    legend.position = c("left", "top"),
                    attr.outside= FALSE,  attr.position= c("RIGHT", "BOTTOM"), 
                    legend.frame = FALSE, legend.width = -0.30, legend.hist.width = 0.24, legend.hist.height = 0.24,
                    inner.margins = c(0.085, 0.2, 0.01, 0.03),
                    frame= TRUE,
                    design.mode= FALSE))
  dev.off()
  # remove a. a will be reused in the next itteration cycle
  rm(a)
  # include message to keep track of process
  message('Processed ', shapefile, ' of ', length(shape_list2), " for ens_50")
  ###
  # if more columns in the shapefile shall be printed to an individual map, 
  # copy the above lines from:
  ## a is a placeholder for the different shapefiles
  ##a <- st_read(shape_list2[shapefile])
  # change the parameters below, especially tm_fill(...)
}

  1. More on climate models, scenarios and projections can be found here.↩︎

Avatar
Nils Riach
Research associate

My research interests are climate change vulnerability, GIS methods, natural hazards and global change research matter.

Related