bivariate_map.R

# Originally from http://rfunctions.blogspot.fi/2015/03/bivariate-maps-bivariatemap-function.html

# Install:
#install.packages("classInt")
#install.packages("raster")
#install.packages("rgdal")
#install.packages("dismo")
#install.packages("XML")
#install.packages("maps")
#install.packages("sp")

# Load:

library(classInt)
library(dismo)
library(maps)
library(raster)
library(rgdal)
library(sp)
library(XML)

colmat <- function(nquantiles = 10,
                   upperleft = rgb(0,150,235, maxColorValue = 255),
                   upperright = rgb(130,0,80, maxColorValue = 255),
                   bottomleft = "grey",
                   bottomright = rgb(255,230,15, maxColorValue = 255),
                   xlab = "x label", ylab = "y label") {

  my.data <- seq(0,1,.01)
  my.class <- classIntervals(my.data, n = nquantiles, style = "quantile")
  my.pal.1 <- findColours(my.class, c(upperleft, bottomleft))
  my.pal.2 <- 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,] <- findColours(my.class, my.col)
  }

  plot(c(1, 1), pch = 19, col = my.pal.1, cex = 0.5, xli = c(0,1),
       ylim = c(0,1), frame.plot = FALSE, xlab = xlab, ylab = ylab,
       cex.lab = 1.3)

  for (i in 1:101) {
    col.temp <- col.matrix[i - 1,]
    points(my.data, rep((i - 1) / 100, 101), pch = 15, col = col.temp, cex = 1)}
  seqs <- seq(0, 100, (100 / nquantiles))
  seqs[1] <- 1
  col.matrix <- col.matrix[c(seqs), c(seqs)]
}

# You can specify the number of quantiles, colors and labels of your color
# matrix. Example
col.matrix <- colmat(nquantiles = 10, upperleft = "blue", upperright = "yellow",
                     bottomleft = "green", bottomright = "red",
                     xlab = "My x label", ylab = "My y label")

# But let's use this simple code. You can change "nquantiles" to generate color
# matrices with different color schemes. For example, change it to 4 to produce
# a 4x4 color scheme.
col.matrix <- colmat(nquantiles = 10)

# Plot bivariate map

bivariate.map <- function(rasterx, rastery, colormatrix = col.matrix,
                          nquantiles = 10) {

  quanmean <- getValues(rasterx)
  temp <- data.frame(quanmean, quantile = rep(NA, length(quanmean)))
  brks <- with(temp, quantile(temp, na.rm = TRUE,
                              probs = c(seq(0, 1, 1 / nquantiles))))
  r1 <- within(temp, quantile <- cut(quanmean, breaks = brks,
                                     labels = 2:length(brks),
                                     include.lowest = TRUE))
  quantr <- data.frame(r1[,2])

  quanvar <- getValues(rastery)
  temp <- data.frame(quanvar, quantile = rep(NA, length(quanvar)))
  brks <- with(temp, quantile(temp,na.rm = TRUE,
                              probs = c(seq(0, 1, 1 / nquantiles))))
  r2 <- within(temp, quantile <- cut(quanvar, breaks = brks,
                                     labels = 2:length(brks),
                                     include.lowest = TRUE))
  quantr2 <- data.frame(r2[, 2])

  as.numeric.factor <- function(x) {as.numeric(levels(x))[x]}
  col.matrix2 <- colormatrix
  cn <- unique(colormatrix)

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

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

# Get the data
dir.create("data")
download.file("http://download1324.mediafire.com/uo8sq368sb7g/4dpxqa38vue7ne6/rasters.zip",
              "data/rasters.zip")
unzip("data/rasters.zip", exdir = "data")


# Load the first raster ("amphibians.grd"):
raster.amphibians <- raster("data/amphibians.grd")

# Load the second raster ("reptiles.grd")
raster.reptiles <- raster("data/reptiles.grd")

# Plot the first raster. You can see that most amphibians are in northwest
# (Amazon =]) and southeast (Atlantic forest =]) Brazil.
my.colors = colorRampPalette(c("white", "lightblue", "yellow", "orangered",
                               "red"))
plot(raster.amphibians, frame.plot = FALSE, axes = FALSE, box = FALSE,
     add = FALSE, legend.width = 1, legend.shrink = 1, col = my.colors(255))
map(interior = TRUE, add = TRUE)

# Plot the second raster. You can see that most reptiles are in northwest
# (Amazon) Brazil.
plot(raster.reptiles, frame.plot = FALSE, axes = FALSE, box = FALSE,
     add = FALSE, legend.width = 1, legend.shrink = 1, col = my.colors(255))
map(interior = TRUE, add = TRUE)

# Use the bivariate.map function:
bivmap <- bivariate.map(raster.amphibians, raster.reptiles,
                        colormatrix = col.matrix, nquantiles = 10)

# Plot the bivariate map:
p1 <- plot(bivmap, frame.plot = FALSE, axes = FALSE, box = FALSE, add = FALSE,
           legend = FALSE, col = as.vector(col.matrix))
m <- map(interior = TRUE, add = TRUE)