Make a bivariate plot using raster data and ggplot2

bivarRasterPlot.R

# The function that produces the colour matrix
colmat <- function(nbreaks = 3, breakstyle = "quantile",
                   upperleft = "#0096EB", upperright = "#820050", 
                   bottomleft = "#BEBEBE", bottomright = "#FFE60F",
                   xlab = "x label", ylab = "y label", plotLeg = TRUE,
                   saveLeg = FALSE) {
   # TODO - replace any tidyr, dplyr etc. functions with data.table #
  library(tidyverse)
  require(ggplot2)
  require(classInt)
  if (breakstyle == "sd") {
    warning("SD breaks style cannot be used.\nWill not always return the correct number of breaks.\nSee classInt::classIntervals() for details.\nResetting to quantile",
            call. = FALSE, immediate. = FALSE)
    breakstyle <- "quantile"}
  # The colours can be changed by changing the HEX codes for:
  # upperleft, upperright, bottomleft, bottomright
  # From http://www.joshuastevens.net/cartography/make-a-bivariate-choropleth-map/
  # upperleft = "#64ACBE", upperright = "#574249", bottomleft = "#E8E8E8", bottomright = "#C85A5A",
  # upperleft = "#BE64AC", upperright = "#3B4994", bottomleft = "#E8E8E8", bottomright = "#5AC8C8",
  # upperleft = "#73AE80", upperright = "#2A5A5B", bottomleft = "#E8E8E8", bottomright = "#6C83B5", 
  # upperleft = "#9972AF", upperright = "#804D36", bottomleft = "#E8E8E8", bottomright = "#C8B35A",
  # upperleft = "#DA8DC8", upperright = "#697AA2", bottomleft = "#E8E8E8", bottomright = "#73BCA0",
  # Similar to Teuling, Stockli, Seneviratnea (2011) [https://doi.org/10.1002/joc.2153]
  # upperleft = "#F7900A", upperright = "#993A65", bottomleft = "#44B360", bottomright = "#3A88B5",
  # Viridis style
  # upperleft = "#FEF287", upperright = "#21908D", bottomleft = "#E8F4F3", bottomright = "#9874A1",
  # Similar to Fjeldsa, Bowie, Rahbek 2012
  # upperleft = "#34C21B", upperright = "#FFFFFF", bottomleft = "#595757",  bottomright = "#A874B8",
  # Default from original source
  # upperleft = "#0096EB", upperright = "#820050", bottomleft= "#BEBEBE", bottomright = "#FFE60F",
  my.data <- seq(0, 1, .01)
  # Default uses terciles (Lucchesi and Wikle [2017] doi: 10.1002/sta4.150)
  my.class <- classInt::classIntervals(my.data,
                                       n = nbreaks,
                                       style = breakstyle,
                                       )
  my.pal.1 <- classInt::findColours(my.class, c(upperleft, bottomleft))
  my.pal.2 <- classInt::findColours(my.class, c(upperright, bottomright))
  col.matrix <- matrix(nrow = 101, ncol = 101, NA)
  for (i in 1:101) {
    my.col <- c(paste(my.pal.1[i]), paste(my.pal.2[i]))
    col.matrix[102 - i, ] <- classInt::findColours(my.class, my.col)
  }
  ## need to convert this to data.table at some stage.
  col.matrix.plot <- col.matrix %>%
    as.data.frame(.) %>% 
    mutate("Y" = row_number()) %>%
    mutate_at(.tbl = ., .vars = vars(starts_with("V")), .funs = list(as.character)) %>% 
    pivot_longer(data = ., cols = -Y, names_to = "X", values_to = "HEXCode") %>% 
    mutate("X" = as.integer(sub("V", "", .$X))) %>%
    distinct(as.factor(HEXCode), .keep_all = TRUE) %>%
    mutate(Y = rev(.$Y)) %>% 
    dplyr::select(-c(4)) %>%
    mutate("Y" = rep(seq(from = 1, to = nbreaks, by = 1), each = nbreaks),
           "X" = rep(seq(from = 1, to = nbreaks, by = 1), times = nbreaks)) %>%
    mutate("UID" = row_number())
  # Use plotLeg if you want a preview of the legend
  if (plotLeg) {
    p <- ggplot(col.matrix.plot, aes(X, Y, fill = HEXCode)) +
      geom_tile() +
      scale_fill_identity() +
      coord_equal(expand = FALSE) +
      theme_void() +
      theme(aspect.ratio = 1,
            axis.title = element_text(size = 12, colour = "black",hjust = 0.5, 
                                      vjust = 1),
            axis.title.y = element_text(angle = 90, hjust = 0.5)) +
      xlab(bquote(.(xlab) ~  symbol("\256"))) +
      ylab(bquote(.(ylab) ~  symbol("\256")))
    print(p)
    assign(
      x = "BivLegend",
      value = p,
      pos = .GlobalEnv
    )
  }
  # Use saveLeg if you want to save a copy of the legend
  if (saveLeg) {
    ggsave(filename = "bivLegend.pdf", plot = p, device = "pdf",
           path = "./", width = 4, height = 4, units = "in",
           dpi = 300)
  }
  seqs <- seq(0, 100, (100 / nbreaks))
  seqs[1] <- 1
  col.matrix <- col.matrix[c(seqs), c(seqs)]
  attr(col.matrix, "breakstyle") <- breakstyle
  attr(col.matrix, "nbreaks") <- nbreaks
  return(col.matrix)
}

# Function to assign colour-codes to raster data
# As before, by default assign tercile breaks
bivariate.map <- function(rasterx, rastery, colourmatrix = col.matrix,
                          export.colour.matrix = TRUE,
                          outname = paste0("colMatrix_rasValues", names(rasterx))) {
  # TO DO - replace raster with terra #
  require(raster)
  require(classInt)
  # export.colour.matrix will export a data.frame of rastervalues and RGB codes 
  # to the global environment outname defines the name of the data.frame
  quanx <- getValues(rasterx)
  tempx <- data.frame(quanx, quantile = rep(NA, length(quanx)))
  brks <- with(tempx, classIntervals(quanx,
                                    n = attr(colourmatrix, "nbreaks"),
                                    style = attr(colourmatrix, "breakstyle"))$brks)
  ## Add (very) small amount of noise to all but the first break
  ## https://stackoverflow.com/a/19846365/1710632
  brks[-1] <- brks[-1] + seq_along(brks[-1]) * .Machine$double.eps
  r1 <- within(tempx, quantile <- cut(quanx,
                                     breaks = brks,
                                     labels = 2:length(brks),
                                     include.lowest = TRUE))
  quantr <- data.frame(r1[, 2])
  quany <- getValues(rastery)
  tempy <- data.frame(quany, quantile = rep(NA, length(quany)))
  brksy <- with(tempy, classIntervals(quany,
                                     n = attr(colourmatrix, "nbreaks"),
                                     style = attr(colourmatrix, "breakstyle"))$brks)
  brksy[-1] <- brksy[-1] + seq_along(brksy[-1]) * .Machine$double.eps
  r2 <- within(tempy, quantile <- cut(quany,
                                     breaks = brksy,
                                     labels = 2:length(brksy),
                                     include.lowest = TRUE
  ))
  quantr2 <- data.frame(r2[, 2])
  as.numeric.factor <- function(x) {
    as.numeric(levels(x))[x]
  }
  col.matrix2 <- colourmatrix
  cn <- unique(colourmatrix)
  for (i in 1:length(col.matrix2)) {
    ifelse(is.na(col.matrix2[i]),
           col.matrix2[i] <- 1, col.matrix2[i] <- which(
             col.matrix2[i] == cn
           )[1]
    )
  }
  # Export the colour.matrix to data.frame() in the global env
  # Can then save with write.table() and use in ArcMap/QGIS
  # Need to save the output raster as integer data-type
  if (export.colour.matrix) {
    # create a dataframe of colours corresponding to raster values
    exportCols <- as.data.frame(cbind(
      as.vector(col.matrix2), as.vector(colourmatrix),
      t(col2rgb(as.vector(colourmatrix)))
    ))
    # rename columns of data.frame()
    colnames(exportCols)[1:2] <- c("rasValue", "HEX")
    # Export to the global environment
    assign(
      x = outname,
      value = exportCols,
      pos = .GlobalEnv
    )
  }
  cols <- numeric(length(quantr[, 1]))
  for (i in 1:length(quantr[, 1])) {
    a <- as.numeric.factor(quantr[i, 1])
    b <- as.numeric.factor(quantr2[i, 1])
    cols[i] <- as.numeric(col.matrix2[b, a])
  }
  r <- rasterx
  r[1:length(r)] <- cols
  return(r)
}

#### TESTS ####
example <- FALSE
if (example) {
  library(data.table)
  library(tidyverse)
  library(raster)
  library(classInt)
  library(patchwork)
  
  # Retrieve some raster data from worldclim
  # and clip to SE Asia and Northern Australia
  clipExt <- extent(52, 180, -22, 39)
  r <- getData("worldclim", var = "bio", res = 10)
  r <- crop(r[[c(1, 12)]], clipExt)
  r[[1]] <- r[[1]]/10
  names(r) <- c("Temp","Prec")
  
  # Define the number of breaks
  nBreaks <- 5
  
  # Create the colour matrix
  col.matrixQ <- colmat(nbreaks = nBreaks, breakstyle = "quantile",
                       xlab = "Temperature", ylab = "Precipiation", 
                       bottomright = "#F7900A", upperright = "#993A65",
                       bottomleft = "#44B360", upperleft = "#3A88B5",
                       saveLeg = FALSE, plotLeg = TRUE)
  
  # create the bivariate raster
  bivmapQ <- bivariate.map(rasterx = r[["Temp"]], rastery = r[["Prec"]],
                          export.colour.matrix = FALSE,
                          colourmatrix = col.matrixQ)
  
  # Convert to dataframe for plotting with ggplot
  bivMapDFQ <- setDT(as.data.frame(bivmapQ, xy = TRUE))
  colnames(bivMapDFQ)[3] <- "BivValue"
  bivMapDFQ <- melt(bivMapDFQ, id.vars = c("x", "y"),
                   measure.vars = "BivValue",
                   value.name = "bivVal",
                   variable.name = "Variable")
  
  # Make the map using ggplot
  map_q <- ggplot(bivMapDFQ, aes(x = x, y = y)) +
    geom_raster(aes(fill = bivVal)) +
    scale_y_continuous(breaks = seq(-20, 60, by = 10), 
                       labels = paste0(seq(-20, 60, 10), "°")) +
    scale_x_continuous(breaks = seq(50,175,25), 
                       labels = paste0(seq(50,175,25), "°")) +
    scale_fill_gradientn(colours = col.matrixQ, na.value = "transparent") + 
    theme_bw() +
    theme(text = element_text(size = 10, colour = "black")) +
    borders(colour = "black", size = 0.5) +
    coord_quickmap(expand = FALSE, xlim = clipExt[1:2], ylim = clipExt[3:4]) +
    theme(legend.position = "none",
          plot.background = element_blank(),
          strip.text = element_text(size = 12, colour = "black"),
          axis.text.y = element_text(angle = 90, hjust = 0.5),
          axis.text = element_text(size = 12, colour = "black"),
          axis.title = element_text(size = 12, colour = "black")) +
    labs(x = "Longitude", y = "Latitude")
  
  # Using different breaks algorithm
  # Create the colour matrix
  col.matrixF <- colmat(nbreaks = nBreaks, breakstyle = "fisher",
                       xlab = "Temperature", ylab = "Precipiation", 
                       bottomright = "#F7900A", upperright = "#993A65",
                       bottomleft = "#44B360", upperleft = "#3A88B5",
                       saveLeg = FALSE, plotLeg = TRUE)
  
  # create the bivariate raster
  bivmapF <- bivariate.map(rasterx = r[["Temp"]], rastery = r[["Prec"]],
                          export.colour.matrix = FALSE,
                          colourmatrix = col.matrixF)
  bivMapDFF <- setDT(as.data.frame(bivmapF, xy = TRUE))
  colnames(bivMapDFF)[3] <- "BivValue"
  bivMapDFF <- melt(bivMapDFF, id.vars = c("x", "y"),
                   measure.vars = "BivValue",
                   value.name = "bivVal",
                   variable.name = "Variable")
  
  # Make the map using ggplot
  map_F <- ggplot(bivMapDFF, aes(x = x, y = y)) +
    geom_raster(aes(fill = bivVal)) +
    scale_y_continuous(breaks = seq(-20, 60, by = 10), 
                       labels = paste0(seq(-20, 60, 10), "°")) +
    scale_x_continuous(breaks = seq(50,175,25), 
                       labels = paste0(seq(50,175,25), "°")) +
    scale_fill_gradientn(colours = col.matrixF, na.value = "transparent") + 
    theme_bw() +
    theme(text = element_text(size = 10, colour = "black")) +
    borders(colour = "black", size = 0.5) +
    coord_quickmap(expand = FALSE, xlim = clipExt[1:2], ylim = clipExt[3:4]) +
    theme(legend.position = "none",
          plot.background = element_blank(),
          strip.text = element_text(size = 12, colour = "black"),
          axis.text.y = element_text(angle = 90, hjust = 0.5),
          axis.text = element_text(size = 12, colour = "black"),
          axis.title = element_text(size = 12, colour = "black")) +
    labs(x = "Longitude", y = "Latitude")
  
  fig <- {{map_q + ggtitle("Quantile breaks")} / {map_F + ggtitle("Fisher breaks")} } + 
    inset_element(BivLegend + theme(plot.background = element_rect(fill = "white",
                                                                   colour = NA)), 
                  left = 0.5, bottom = 0.5, right = 1.2, top = 0.90,
                  align_to = "full") +
    plot_annotation(caption = "Both plots are made on the same data, but have breaks defined differently")
  fig
  # Save
  ggsave(plot = fig,
         filename = "BivariatePlot_ggsave.pdf",
         device = "pdf", path = "./",
         width = 6, height = 7, units = "in",
         dpi = 320)
}

hi @lukasbaumbach - sorry for the slow reply. I'm not sure how I missed your message.

Long story short, I don't think it does matter. Here is a comparison of my method, and the equivalent method from the biscale package. The biscale method uses the underlying data to define the breaks and then produce the classification. You can see the plots are effectively the same - I think the slight difference is due to the way I interpolate the colours between the high and low values along each of the axes.

require(classInt)
require(data.table)
require(ggplot2)
require(raster)
library(biscale)
library(cowplot)
library(patchwork)

# get some data
clipExt <- extent(52, 180, -22, 39)
r <- getData("worldclim", var = "bio", res = 10)
r <- crop(r[[c(1, 12)]], clipExt)
r[[1]] <- r[[1]]/10
names(r) <- c("Temp","Prec")
dat <- setDT(as.data.frame(r, xy = TRUE, na.rm = TRUE))

# Define the number of breaks
nBreaks <- 3

# my method
# Create the colour matrix
col.matrixQ <- colmat(nbreaks = nBreaks, breakstyle = "fisher",
                      xlab = "Temperature", ylab = "Precipiation", 
                      bottomright = "#AE3A4E", upperright = "#3F2949",
                      bottomleft = "#CABED0", upperleft = "#4885C1",
                      saveLeg = FALSE, plotLeg = TRUE)

# create the bivariate raster
bivmapQ <- bivariate.map(rasterx = r[["Temp"]], rastery = r[["Prec"]],
                         export.colour.matrix = FALSE,
                         colourmatrix = col.matrixQ)

# Convert to dataframe for plotting with ggplot
bivMapDFQ <- setDT(as.data.frame(bivmapQ, xy = TRUE))
colnames(bivMapDFQ)[3] <- "BivValue"
bivMapDFQ <- melt(bivMapDFQ, id.vars = c("x", "y"),
                  measure.vars = "BivValue",
                  value.name = "bivVal",
                  variable.name = "Variable")

# Make the map using ggplot
map_q <- ggplot(bivMapDFQ, aes(x = x, y = y)) +
  geom_raster(aes(fill = bivVal)) +
  scale_y_continuous(breaks = seq(-20, 60, by = 10), 
                     labels = paste0(seq(-20, 60, 10), "°")) +
  scale_x_continuous(breaks = seq(50,175,25), 
                     labels = paste0(seq(50,175,25), "°")) +
  scale_fill_gradientn(colours = col.matrixQ, na.value = "transparent") + 
  theme_bw() +
  theme(text = element_text(size = 10, colour = "black")) +
  borders(colour = "black", size = 0.5) +
  coord_quickmap(expand = FALSE, xlim = clipExt[1:2], ylim = clipExt[3:4]) +
  theme(legend.position = "none",
        plot.background = element_blank(),
        strip.text = element_text(size = 12, colour = "black"),
        axis.text.y = element_text(angle = 90, hjust = 0.5),
        axis.text = element_text(size = 12, colour = "black"),
        axis.title = element_text(size = 12, colour = "black")) +
  labs(x = "Longitude", y = "Latitude")
map_q

myPlot <- cowplot::ggdraw() +
  cowplot::draw_plot(map_q, 0, 0, 1, 1) +
  cowplot::draw_plot(BivLegend, 0.075, .3, 0.2, 0.2)
myPlot

# bi_scale method
bi <- biscale::bi_class(dat, x = "Temp", y = "Prec",
                        style = "fisher", dim = nBreaks)
map_bi <- ggplot(bi) +
  geom_tile(aes(x, y, fill = bi_class), show.legend = FALSE) +
  bi_scale_fill(pal = "DkViolet") +
  theme_bw() +
  borders(colour = "black", size = 0.5) +
  coord_quickmap(expand = FALSE, xlim = clipExt[1:2], ylim = clipExt[3:4]) +
  scale_y_continuous(breaks = seq(-20, 60, by = 10), 
                     labels = paste0(seq(-20, 60, 10), "°")) +
  scale_x_continuous(breaks = seq(50,175,25), 
                     labels = paste0(seq(50,175,25), "°")) +
  theme(legend.position = "none",
        plot.background = element_blank(),
        strip.text = element_text(size = 12, colour = "black"),
        axis.text.y = element_text(angle = 90, hjust = 0.5),
        axis.text = element_text(size = 12, colour = "black"),
        axis.title = element_text(size = 12, colour = "black")) +
  labs(x = "Longitude", y = "Latitude")
map_bi

bi_leg <- bi_legend(pal = "DkViolet",
                    dim = 3,
                    xlab = "Temperature ",
                    ylab = "Precipitation ",
                    size = 8)
biPlot <- cowplot::ggdraw() +
  cowplot::draw_plot(map_bi, 0, 0, 1, 1) +
  cowplot::draw_plot(bi_leg, 0.075, .3, 0.2, 0.2)
biPlot

patchwork::wrap_plots(map_q, map_bi, ncol = 1)

Hello, thank you very much for the function. Overall, it works fine I just don't understand why bivValue of the generated raster contains integer values with a range of 1 to 29 if I use a 5x5 color matrix with equal breaks. Shouldn't it be 25 because one for every colour within the matrix or 26 if we would include another one for NA? Best and thank you in advance esoxx

This gist seems to be derived from this post from 2015: https://rfunctions.blogspot.com/2015/03/bivariate-maps-bivariatemap-function.html. Just a heads up that the author of that made it into a package that might be of interest: https://cran.r-project.org/package=bivariatemaps

Dear Stuart, Thank you for such a amazing function in the map-plot! It is so usful for my map. I used it for gird map, it works well. This time, I want to plot for hexagon plots, but it doesnor work. Is there something wrong in my code?
[bivmap<-bivariate.map(ABS.snake.richness.5,ABS.snake.abs.5,colormatrix=col.matrix, nquantiles=10)
sorry， i failed to attach thefile， but i attached a picture of thedetails of the two hexagon files.
Can you help me have a look?

Any help woule be appreciate
Best
Tao

Hello, thanks for the function!
after I fixed the breaks problem, had this Error when trying to run the bivmap function using my rasters.
Did any one have the same error? know how to fix this?

Error in if (n < 2 & needn) stop("n less than 2") :
argument is of length zero

Cheers
J