R中的双变量等值线图

Question

我正在寻找使用光栅文件在 R 中创建 bivariate choropleth maps 的通用解决方案。

我找到了以下代码here，它几乎可以满足我的需要，但它是有限的：它只能处理两个轴上介于 0 和 1 之间的数据。在我的具体用例中，一个轴跨越 0-1，而另一个轴跨越 -1 和 1。不管我的具体用例如何，我认为可以处理不同数据范围的更通用的函数对许多人有用。

我已经尝试更新函数colmat 中的代码来处理负面数据，但我终其一生都无法让它发挥作用。为了清楚起见，我避免发布所有失败的尝试，并在我在上面链接中找到的代码下方复制了 insread，希望有人能够提供解决方案。

当前代码首先使用colmat 创建一个颜色矩阵。然后在bivariate.map 中使用生成的颜色矩阵以及包含数据的两个光栅文件。我认为理想的解决方案是首先基于两个栅格创建颜色矩阵（以便它可以根据您的实际数据正确地对数据进行分类，而不是当前介于 0 和 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)) ````

理想情况下，更通用的函数将采用两个栅格文件，确定两者的数据范围，然后根据用户指定的 bin/quantiles 的数量创建一个二元 chorpleth 地图。

Answer 1

这里有一些基于您的代码的想法

三个功能

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

示例数据

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)

现在使用函数

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)