-
-
Save scbrown86/2779137a9378df7b60afd23e0c45c188 to your computer and use it in GitHub Desktop.
# 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! 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.
Great! Thank you very much!
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
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
Would you be able to send me a link to an example file and the code you're using when this happens?