Двусторонняя карта хороплета в R - PullRequest
Новая
       31

Двусторонняя карта хороплета в R

0 голосов
/ 11 января 2019

Я ищу общее решение для создания двумерных картограмм в R с использованием растровых файлов.

Я нашел следующий код здесь , который почти делает то, что мне нужно, но он ограничен: он может обрабатывать только данные, которые находятся между 0 и 1 по обеим осям. В моем конкретном случае использования одна ось охватывает 0-1, а другая - от -1 до 1. Независимо от моего конкретного случая использования, я думаю, что более общая функция, которая может обрабатывать различные диапазоны данных, была бы полезна для многих людей.

Я уже пытался обновить код в функции colmat для обработки отрицательных данных, но я не могу заставить его работать. В целях ясности я избегал публикации всех моих неудачных попыток и перечитал скопированный код под кодом, который я нашел по ссылке выше, в надежде, что кто-то сможет предложить решение.

Текущий код сначала создает цветовую матрицу, используя 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)) ````

В идеале, более общая функция должна принимать два растровых файла, определять диапазоны данных для обоих, а затем создавать двумерную карту последовательности, основываясь на количестве бинов / квантилей, указанных пользователем.

Ответы [ 2 ]

0 голосов
/ 23 мая 2019

Вот несколько идей, основанных на вашем коде

Три функции 

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)

plots

0 голосов
/ 21 мая 2019

Возможно, вы уже разработали решение этой проблемы самостоятельно, и в этом случае я был бы очень рад получить какое-то обновление (сейчас я подхожу к той же проблеме). В противном случае, я просто нашел несколько ссылок, которые, надеюсь, помогут вам разобраться в этом

http://lenkiefer.com/2017/04/24/bivariate-map/

http://rpubs.com/apsteinmetz/prek

Лучший

...