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)
}

Hello! I have been using your code in some maps, and sometimes the bivariate.map function returns an error saying that "Breaks' are not unique". I have interpreted, based on this thread (https://stackoverflow.com/questions/16184947/cut-error-breaks-are-not-unique) that it happens when some classes do not have any pixels, and then the classification does not work. The same thread recommended to use the function "unique()" to wrap the "quantile()" function, however I am afraid that then the colors in the map will not match the colors of the bivariate legend made with the "colmat()" function. For a while, I am working around it decreasing the number of breaks until the error stops.

Would you be able to send me a link to an example file and the code you're using when this happens?

Hi! Yes, here are the codes and files in the attachments.
https://app.box.com/s/illt27n052u9fvzdecu112zacw51vw5n
Best regards!
Vitor

Hi Vitor. The bivariate.map() function won't throw the error anymore. I followed another piece of advice from your link, and it fixed the problem. Your data plots up a little strange though, there are a lot of NA classes.

I tested the function on my example dataset included in the script and it worked fine.

Hi! I explored my data and I discovered that most of the white areas are the ones that have 0 (zero) values in at least one of the input rasters. Maybe the code could be changed to include the 0 (zero) values in the color ramp. Just some few areas such as the interior of Greenland have "no value" pixels, so these ones should really not be colored.

Fixed. The issue was the code to fix the breaks was altering the first break value which meant anything less than the adjusted break value was being excluded.

Great! Thank you very much!

Updated 2020-06-26: The bottom right and upper left colours were being inadvertently swapped during creating of the colour matrix. Issue has now been fixed. Maps were displaying correctly - the issue only affected the legends.

Hi scbrown86, thank you very much for this great function. I would like to add labeled ticks with values to the legend axes (e.g. as here https://stackoverflow.com/questions/48572744/plot-a-bivariate-map-in-r). Is there a way to implement this? I didn't see a direct way, since the legend does not contain the actual data anymore, only the respective color codes.

Hi Stuart! Thanks for this, it's a really useful piece of code! Unfortunately I am having the same issues as influence-of-canopy-on-temperature regarding 0 (zero) values in at least one of the input rasters leading to NAs. Where/what does the code need to be changed to to fix this? Forgive me, I'm fairly new to GIS! Thanks in advance!

@rhys-pa, can you attach an image/example dataset, and the code you've used that's giving the errors? Using the code below, where I force a large extent to all have 0 values across both layers works fine.

clipExt <- extent(52, 180, -22, 39)
  r <- getData("worldclim", var = "bio", res = 10)
  r <- r[[c(1, 12)]] %>%
    crop(., clipExt)
  r[[1]] <- r[[1]]/10
  names(r) <- c("Temp","Prec")
  # Make all values 0 for the following cells
  idx <- cellsFromExtent(r, extent(70,87,22,36))
  r[idx] <- 0
  plot(r)
  # Define the number of breaks
  nBreaks <- 5
  col.matrix <- colmat(nquantiles = nBreaks, xlab = "Temp", ylab = "Precip", 
                       ## non default colours
                       saveLeg = FALSE, plotLeg = TRUE)
  bivmap <- bivariate.map(rasterx = r[["Temp"]], rastery = r[["Prec"]],
                          export.colour.matrix = TRUE, outname = "bivMapCols",
                          colormatrix = col.matrix, nquantiles = nBreaks)
  bivMapDF <- as.data.frame(bivmap, xy = TRUE) %>%
    tbl_df() %>%
    dplyr::rename("BivValue" = 3) %>%
    pivot_longer(., names_to = "Variable", values_to = "bivVal", cols = BivValue)
  map <- ggplot(bivMapDF, 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.matrix, 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 <- ggdraw(map) + 
    draw_plot(BivLegend +
                theme(plot.background = element_rect(fill = "white", colour = NA)),
              width = 0.25, height = 0.25, x = 0.75, y = 0.55)
  fig

Updated 2021-07-09: Can now specify a number of different types of breaks. Previously only quantiles could be used. The default method is still terciles, but with the exception of sd any of the methods available in classInt::classIntervals can now be used. To simplify things when using bivariate.map you no longer need to specify the number and style of breaks, these are automatically imported from the colourmatrix argument.

Hi Stuart,

This is a great job! I am using your code for a wildlife study. Everything is going well. I have a couple of questions about your code - has the X and Y axis data been normalized? Also, about your recent update, what break styles can be done other than the ones you are showing in the tutorial?

Cheers from Costa Rica!

Hi @Randall-HYLA, the data is not normalised in any way before plotting. If you want to normalise the data I would do so when creating the bivariate.map object; e.g. something like

bivariate.map(rasterx = r[["Temp"]], rastery = sqrt(r[["Prec"]]),
                          export.colour.matrix = FALSE,
                          colourmatrix = col.matrixF)

Note though, that the breaks on the legend will always appear equally spaced because of the way I make the legend, but if you want to get an idea of what they look like with regards to you data, you could do something like this:

nClass <- 5
tempCols <- c("#44B360", "#F7900A")
rainCols <- c("#44B360", "#3A88B5")
tempDensity <- stats::density(values(r[["Temp"]]), na.rm = TRUE)
rainDensity <- stats::density(values(r[["Prec"]]), na.rm = TRUE)
par(mfrow = c(2, 1), mai = c(0.8,0.5,0.5,0.1), mgp = c(1, 1, 0))
plot(tempDensity, 
     ylim = range(tempDensity$y),
     xlim = range(tempDensity$x),
     ylab = "", main = "", axes = FALSE,
     xlab = "Temperature")
rug(x = values(r[["Temp"]]), side = 3, col = "#A3A3A368")
clBreaks <- {set.seed(54);classInt::classIntervals(values(r[["Temp"]]), 
                                                   n = nClass,
                                                   style = "fisher",
                                                   samp_prop = 1)}
abline(v = clBreaks$brks[-1],
       col = attr(classInt::findColours(clBreaks, pal = tempCols), "palette"))
axis(2); axis(3); box()
par(mai = c(0.8, 0.5, 0.5, 0.1), mgp = c(1, 1, 0))
plot(rainDensity, 
     ylim = range(rainDensity$y),
     xlim = range(rainDensity$x),
     ylab = "", main = "", axes = FALSE,
     xlab = "Precipitation")
rug(x = values(r[["Prec"]]), side = 3, col = "#A3A3A368")
clBreaks <- {set.seed(54);classInt::classIntervals(values(r[["Prec"]]), 
                                                   n = nClass,
                                                   style = "fisher",
                                                   samp_prop = 1)}
abline(v = clBreaks$brks[-1],
       col = attr(classInt::findColours(clBreaks, pal = rainCols), "palette"))
axis(2); axis(3); box()

You can read the documentation for classInt::classIntervals() to see which break styles are supported. The only ones that won't work (in my limited testing) are fixed and sd.

Hi @scbrown86, thank you very much for your help above. I have been playing around with this code for my analysis and it has been going very well. I have one more question. Do you think it is possible to add an extra color in the middle of the color matrix to create a divergent matrix. This would significantly improve the patterns I am looking for in my analysis. So far I have not been able to do this. Thanks!

Hi @Randall-HYLA, it can be done. You just need to make sure you've got some diverging colours specified. It might take a bit of trial and error, but here is one I've used previously that has a diverging x-axis and an increasing y axis

col.matrix <- colmat(nbreaks = 4, 
                     breakstyle = "fisher",
                     xlab = "x", 
                     ylab = "y",
                     bottomleft = "#c3b3d8", upperleft = "#240d5e",
                     bottomright = "#ffcc80", upperright = "#b30000",
                     saveLeg = FALSE, plotLeg = TRUE,
                     ## sampling argument to classInt::classIntervals
                     samp_prop = 1)

Hi @scbrown86, thanks a lot for this code! I would like to ask, is it impossible to produce a bivariate map for rasters that have very little variance? I have tried all changes suggested in your post and in the link (https://stackoverflow.com/questions/16184947/cut-error-breaks-are-not-unique/19846365#19846365) without success. With some of the suggested changes I still get the "Cut() error - 'breaks' are not unique", and when using .bincode I don't get any error, the bivariate.map function seems to run but the final plot is empty. I added two of the rasters I would like to plot to this link in case you have a chance to check them: https://www.dropbox.com/sh/ucxajtyf7o0wqyf/AAADMXU8XGQyHLJVkPRH5muQa?dl=0

Thanks!!

Hi @ilolivaresa, I can run your rasters through the function without getting any errors/warnings. I tried 3 and 5 breaks using both quantile and fisher breaks and the map plots up fine. How many breaks are you trying to use?

In my testing the function works fine on the original resolution using quantile breaks and at an aggregated 0.5° resolution using quantile and fisher breaks. Given enough time, I suspect the function would work fine using fisher breaks at the original resolution but it takes a long time to cut the vector of values into their classes using the fisher algorithm.

The plot below is the quantile and fisher breaks at the aggregated resolution.

Without any more detail on what you've tried I really can't provide any more help.

@scbrown86 I have been using your code with a soil dataset and I would say it works perfectly. Thank you so much for the contribution and time.

Hi @scbrown86 Thanks a lot for your reply! Yes, you were right, the maps I need are those produced with fisher breaks, it simply took a long time. Thanks so much for taking the time to check my rasters.

Hello, I have this error. Could you help me solve it? I don't change anything.

Warning messages: 1: In classIntervals(quanx, n = attr(colourmatrix, "nbreaks"), style = attr(colourmatrix, : var has missing values, omitted in finding classes 2: In classIntervals(quany, n = attr(colourmatrix, "nbreaks"), style = attr(colourmatrix, : var has missing values, omitted in finding classes

@HaoLi2025 thats a warning not an error, just letting you know that there is missing data in your rasters.

@scbrown86 First off, thanks a lot for the effort you put into preparing this function. Maybe I'm mistaken, but it seems to me, that the legend produced only matches up with the data for the case when breakstyle="equal". Since the legend breaks will depend on the result of classInt::classIntervals for raster X and raster Y, the legend should only be drawn after identifying those breaks, right? Best, Lukas

I replace the NA with the '0', and it is successful.

…

________________________________ 发件人: Lukas Baumbach ***@***.***> 发送时间: 2021年11月30日 14:25 收件人: scbrown86 ***@***.***> 抄送: Hao Li ***@***.***>; Mention ***@***.***> 主题: Re: scbrown86/bivarRasterPlot.R @lukasbaumbach commented on this gist.

________________________________ @scbrown86<https://github.com/scbrown86> First off, thanks a lot for the effort you put into preparing this function. Maybe I'm mistaken, but it seems to me, that the legend produced only matches up with the data for the case when breakstyle="equal". Since the legend breaks will depend on the result of classInt::classIntervals for raster X and raster Y, the legend should only be drawn after identifying those breaks, right? Best, Lukas ― You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub<https://gist.github.com/2779137a9378df7b60afd23e0c45c188#gistcomment-3978804>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/AR25WEBYEIP7FIWQLWQD5I3UOTGC5ANCNFSM4MNOHKLQ>. Triage notifications on the go with GitHub Mobile for iOS<https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android<https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

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