This tutorial is about making multivariate maps. All data can be downloaded in R and most parameters are custom. These plots hold an abundance of information and well suited to explain complex environmental patterns. In fact, we use them as a starting point for things like natural risk delineation. When things like soil are added, it strengthens the meaning of such maps.

In this tutorial, we use the country of Iran to map precipitation, temperature and elevation in ggplot. Since the colors need to make sense and I am color blind, the colors maybe not the best combination but this can be changed by the user.

library(geodata) #download data
library(sf) #hndle vector data
library(terra) #handle raster data
library(tidyterra) #plot raster data
library(biscale) #for biplots
library(cowplot) #to add to ggplot
library(ggnewscale) #create multiple color scales
library(ggspatial) #mapping
library(tidyverse) #tidying and many essentials

First step is to load all the data, you can use your own or packages like elevatr, geodata, tigris, etcs. Here we only need geodata. We need DEM, Climate and if you want continuous soil data or non ranked ordinal soil data (e.g., soil texture).

We will make a multivariate map of climate (MAP and MAT) and elevation. Large areas will take some time as they need to be downloaded but in general it is a fairly rapid code.

#Region polygon
iran = gadm("Iran", level = 0, path = tempdir())%>%
  st_as_sf()%>%
  select(COUNTRY)

#DEM
dem = elevation_30s("Iran", path = tempdir())%>%
  project(., "EPSG:4326")%>%
  setNames("DEM")

#Climate - first temp (c)
mat = worldclim_country("Iran", var = "tavg", path = tempdir())%>%
  median()%>%
  crop(., iran)%>%
  mask(., iran)%>%
  setNames("MAT")

#add precipitation (mm) together
clim = worldclim_country("Iran", var = "prec", path = tempdir())%>%
  median()%>%
  crop(., iran)%>%
  mask(., iran)%>%
  setNames("MAP")%>%
  c(mat, .)%>%
  resample(., dem, method = 'lanczos')

#See how climate looks
plot(clim)

Now we need to make hillshade from slope and aspect of the DEM. This is to apply terrain to the map to add more information in a simple way. Here we use degrees and not radians as in larger regions it enhances the topography so it is more clear. At the regional and local scale radians are suggested because degrees will give a speckled appearance.

First we make the multi angle hillshade from the DEM.

#hillshade function
hillmulti = function(x, angles = c(270, 15, 60, 330), units = "radians"){
  
  top = terrain(x, v = c('slope',"aspect"), unit = units)
  
  hst = sapply(angles, function(dir){
    shade(top[[1]], top[[2]],
          angle = 45,
          direction = dir,
          normalize = TRUE)
  })
  sum(rast(hst))
}

#Run function with degrees
hillSlope = hillmulti(dem, units = "degrees")

#Plot to see if it will work
plot(hillSlope, col = grey(1:100/100), legend = F, main = "Hillshade")

There are a couple ways to make a multivariate map. For this tutorial we will have precipitation on the y axis, temp on the x axis and elevation as the z (seperate to the right). Another way of doing this is make a color triangle; however, that works best when the characteristics have peaks in different regions. Since these factors are often correlated in the same direction we use this method. The first step is to make a biclass of the the quantiles (4) of MAT and MAP but keep the DEM continuous. We than make a figure and legend.

#create a data frame of the raster images
df = as.data.frame(c(clim, dem), xy = T)%>%
  na.omit()

#create classes
cdat = bi_class(df, MAT, MAP, style = "quantile", dim = 4)

#Create a ggplot
img = ggplot()+
  geom_raster(data = cdat, aes(x=x, y = y, fill = bi_class), 
              show.legend = F, inherit.aes = F)+
  bi_scale_fill(pal = "BlueYl", dim = 4, guide = "none")+
  new_scale_fill()+
  geom_spatraster(data = dem, alpha = 0.4, show.legend = T, inherit.aes = F)+
  scale_fill_gradientn(name = "Elevation (m)",colours = rev(rainbow(25)), na.value = NA)+
  new_scale_fill()+
  geom_spatraster(data = hillSlope, alpha = 0.4, guide = 'none', 
                  show.legend = F, inherit.aes = F)+
  scale_fill_distiller(palette = "Greys", na.value = NA, guide = 'none')+
  geom_sf(data = iran, fill = NA, color = "black", inherit.aes = F)+
  annotation_north_arrow(location = 'tr')+
  annotation_scale(location = "bl")+
  coord_sf(crs=4326)+
  labs(x = "", y = "", title = "Multivariate climate map")+
  theme_bw()

#plot to see if you like it
img

Now we need to add a legend for the MAT and MAP which is bivariate, so a square gradient. We add it to the map with cowplot and ggdraw.

#Make legend with same dimensions and color palette. 
legend <- bi_legend(pal = "BlueYl",
                    dim = 4,
                    ylab = "MAP",
                    xlab = "MAT",
                    size = 10)

#now we draw it onto the map
fig = ggdraw() +
  draw_plot(img, 0, 0, 1, 1) +
  draw_plot(legend, 0.17, .17, 0.2, 0.2)

#Plot the final map
fig

And we are finished, pretty quick way to learn a lot from one map. Time is the next dimension that can be added, which is a little trickier.