Bivariate Choropleth Map in R

Question

I am looking for a general solution to create bivariate choropleth maps in R using raster files.

I have found the following code here which nearly does what I need but it is limited: it can only handle data which are between 0 and 1 on both axes. In my specific use-case one axis spans 0-1 while another spans between -1 and 1. Regardless as to my specific use-case, I think a more general function which can handle different data ranges would be useful to many people.

I have already tried updating the code within the function colmat to handle negative data but for the life of me cannot get it to work. In the interests of clarity I have avoided posting all of my failed attempts and have insread copied below the code I found at the link above in the hope that someone may be able to offer a solution.

The current code first creates a colour matrix using colmat. The colour matrix generated is then used in bivariate.map along with your two raster files containing the data. I think the ideal solution would be to create the colour matrix based on the two rasters first (so that it can correctly bin the data based on your actual data, not the current solution which is between 0 and 1).

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

# Creates dummy rasters
rasterx<- raster(matrix(rnorm(400),5,5))
rasterx[rasterx <=0]<-1
rastery<- raster(matrix(rnorm(400),5,5))


# This function creates a colour matrix
# At present it cannot handle negative values i.e. the matrix spans from 0 to 1 along both axes
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,xlim=c(0,1),ylim=c(0,1),frame.plot=F, 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)]

}

# Creates colour matrix
col.matrix<-colmat(nquantiles=2, upperleft="blue", upperright="yellow", bottomleft="green", bottomright="red", xlab="Species Richness", ylab="Change in activity hours")

# Function to create bivariate map, given the colour ramp created previously
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)

}

# Creates map
bivmap<-bivariate.map(rasterx,rastery, colormatrix=col.matrix, nquantiles=2)

# Plots a map
plot(bivmap,frame.plot=F,axes=F,box=F,add=F,legend=F,col=as.vector(col.matrix)) ````

Ideally,a more general function would take two raster files, determine the data ranges of both and then create a bivariate chorpleth map based on the number of bins/quantiles specified by the user.

Robert Hijmans Robert Hijmans · Accepted Answer · 2019-05-23T06:24:32

Here are some ideas based on your code

Three functions

makeCM <- function(breaks=10, upperleft, upperright, lowerleft, lowerright) { 
   m <- matrix(ncol=breaks, nrow=breaks)
   b <- breaks-1
   b <- (0:b)/b
   col1 <- rgb(colorRamp(c(upperleft, lowerleft))(b), max=255)
   col2 <- rgb(colorRamp(c(upperright, lowerright))(b), max=255)
   cm <- apply(cbind(col1, col2), 1, function(i) rgb(colorRamp(i)(b), max=255))
   cm[, ncol(cm):1 ]

}

plotCM <- function(cm, xlab="", ylab="", main="") {
    n <- cm
    n <- matrix(1:length(cm), nrow=nrow(cm), byrow=TRUE)
    r <- raster(n)
    cm <- cm[, ncol(cm):1 ]
    image(r, col=cm, axes=FALSE, xlab=xlab, ylab=ylab, main=main)
}


rasterCM <- function(x, y, n) {
    q1 <- quantile(x, seq(0,1,1/(n)))
    q2 <- quantile(y, seq(0,1,1/(n)))
    r1 <- cut(x, q1, include.lowest=TRUE)
    r2 <- cut(y, q2, include.lowest=TRUE)
    overlay(r1, r2, fun=function(i, j) {
        (j-1) * n + i
    })
}

Example data

library(raster)
set.seed(42)
r <- raster(ncol=50, nrow=50, xmn=0, xmx=10, ymn=0,ymx=10, crs="+proj=utm +zone=1")
x <- init(r, "x") * runif(ncell(r), .5, 1)
y <- init(r, "y") * runif(ncell(r), .5, 1)

And now used the functions

breaks <- 5
cmat <- makeCM(breaks, "blue", "yellow", "green", "red")
xy <- rasterCM(x, y, breaks)

par(mfrow=c(2,2), mai=c(.5,.5,.5,.5), las=1)
plot(x)
plot(y)
par(mai=c(1,1,1,1))
plotCM(cmat, "var1", "var2", "legend")
par(mai=c(.5,.5,.5,.5))
image(xy, col=cmat, las=1)

Bivariate Choropleth Map in R

1 Answers