绘制色带在零附近发散的栅格
Plotting a raster with the color ramp diverging around zero
我正在尝试绘制具有正值和负值的地图。
所有正值都应该是红色,而负值应该是蓝色,零应该是白色,就像这个样本图中的离散颜色一样
下面是我使用的代码:
library (rasterVis)
ras1 <- raster(nrow=10,ncol=10)
set.seed(1)
ras1[] <- rchisq(df=10,n=10*10)
ras2=ras1*(-1)/2
s <- stack(ras1,ras2)
levelplot(s,par.settings=RdBuTheme())
非常感谢您提供了一个通用的解决方案,该解决方案也可以应用于其他映射练习。
我写 a gist 来做到这一点。它采用由 rasterVis::levelplot
生成的 trellis
对象和色带,并绘制颜色围绕零发散的对象。
使用你的s
,你可以这样使用它:
devtools::source_gist('306e4b7e69c87b1826db')
p <- levelplot(s)
diverge0(p, ramp='RdBu')
ramp
应该是 RColorBrewer
调色板的名称、要插值的颜色矢量或 colorRampPalette
.
来源如下:
diverge0 <- function(p, ramp) {
# p: a trellis object resulting from rasterVis::levelplot
# ramp: the name of an RColorBrewer palette (as character), a character
# vector of colour names to interpolate, or a colorRampPalette.
require(RColorBrewer)
require(rasterVis)
if(length(ramp)==1 && is.character(ramp) && ramp %in%
row.names(brewer.pal.info)) {
ramp <- suppressWarnings(colorRampPalette(brewer.pal(11, ramp)))
} else if(length(ramp) > 1 && is.character(ramp) && all(ramp %in% colors())) {
ramp <- colorRampPalette(ramp)
} else if(!is.function(ramp))
stop('ramp should be either the name of a RColorBrewer palette, ',
'a vector of colours to be interpolated, or a colorRampPalette.')
rng <- range(p$legend[[1]]$args$key$at)
s <- seq(-max(abs(rng)), max(abs(rng)), len=1001)
i <- findInterval(rng[which.min(abs(rng))], s)
zlim <- switch(which.min(abs(rng)), `1`=i:(1000+1), `2`=1:(i+1))
p$legend[[1]]$args$key$at <- s[zlim]
p$par.settings$regions$col <- ramp(1000)[zlim[-length(zlim)]]
p
}
请注意,正如@LucasFortini 的 中所建议的,如果您乐于让色键在零上方和下方延伸相同的距离,则该过程会简单得多,例如:levelplot(s,par.settings=RdBuTheme(), at=seq(-max(abs(cellStats(s, range))), max(abs(cellStats(s, range))), len=100))
。
这是我经常使用以下脚本执行的操作:
library(colorRamps)
col5 <- colorRampPalette(c('blue', 'gray96', 'red')) #create color ramp starting from blue to red
color_levels=20 #the number of colors to use
max_absolute_value=0.4 #what is the maximum absolute value of raster?
plot(img, col=col5(n=color_levels), breaks=seq(-max_absolute_value,max_absolute_value,length.out=color_levels+1) , axes=FALSE)
使用来自 here 的数据,这是一个示例输出和实际脚本:
library(raster)
library(colorRamps)
mask_data=shapefile("D:/temp/so/Main_Hawaiian_Islands_simple3.shp")
img=raster("D:/temp/so/PPT_wet_minus_dry.tif")
col5 <- colorRampPalette(c('blue', 'gray96', 'red')) #create color ramp starting from blue to red
color_levels=10 #the number of colors to use
max_absolute_value=max(abs(c(cellStats(img, min), cellStats(img, max)))) #what is the maximum absolute value of raster?
color_sequence=seq(-max_absolute_value,max_absolute_value,length.out=color_levels+1)
plot(img, col=col5(n=color_levels), breaks=color_sequence, axes=FALSE)
plot(mask_data, add=T)
这可能会打扰一些人,因为在未使用的负范围内有很多颜色箱(如您提供的示例)。下面的修改允许从地图图例中排除空颜色:
n_in_class=hist(img, breaks=color_sequence, plot=F)$counts>0
col_to_include=min(which(n_in_class==T)):max(which(n_in_class==T))
breaks_to_include=min(which(n_in_class==T)):(max(which(n_in_class==T))+1)
plot(img, col=col5(n=color_levels)[col_to_include], breaks=color_sequence[breaks_to_include] , axes=FALSE)
plot(mask_data, add=T)
我正在尝试绘制具有正值和负值的地图。
所有正值都应该是红色,而负值应该是蓝色,零应该是白色,就像这个样本图中的离散颜色一样
下面是我使用的代码:
library (rasterVis)
ras1 <- raster(nrow=10,ncol=10)
set.seed(1)
ras1[] <- rchisq(df=10,n=10*10)
ras2=ras1*(-1)/2
s <- stack(ras1,ras2)
levelplot(s,par.settings=RdBuTheme())
非常感谢您提供了一个通用的解决方案,该解决方案也可以应用于其他映射练习。
我写 a gist 来做到这一点。它采用由 rasterVis::levelplot
生成的 trellis
对象和色带,并绘制颜色围绕零发散的对象。
使用你的s
,你可以这样使用它:
devtools::source_gist('306e4b7e69c87b1826db')
p <- levelplot(s)
diverge0(p, ramp='RdBu')
ramp
应该是 RColorBrewer
调色板的名称、要插值的颜色矢量或 colorRampPalette
.
来源如下:
diverge0 <- function(p, ramp) {
# p: a trellis object resulting from rasterVis::levelplot
# ramp: the name of an RColorBrewer palette (as character), a character
# vector of colour names to interpolate, or a colorRampPalette.
require(RColorBrewer)
require(rasterVis)
if(length(ramp)==1 && is.character(ramp) && ramp %in%
row.names(brewer.pal.info)) {
ramp <- suppressWarnings(colorRampPalette(brewer.pal(11, ramp)))
} else if(length(ramp) > 1 && is.character(ramp) && all(ramp %in% colors())) {
ramp <- colorRampPalette(ramp)
} else if(!is.function(ramp))
stop('ramp should be either the name of a RColorBrewer palette, ',
'a vector of colours to be interpolated, or a colorRampPalette.')
rng <- range(p$legend[[1]]$args$key$at)
s <- seq(-max(abs(rng)), max(abs(rng)), len=1001)
i <- findInterval(rng[which.min(abs(rng))], s)
zlim <- switch(which.min(abs(rng)), `1`=i:(1000+1), `2`=1:(i+1))
p$legend[[1]]$args$key$at <- s[zlim]
p$par.settings$regions$col <- ramp(1000)[zlim[-length(zlim)]]
p
}
请注意,正如@LucasFortini 的 levelplot(s,par.settings=RdBuTheme(), at=seq(-max(abs(cellStats(s, range))), max(abs(cellStats(s, range))), len=100))
。
这是我经常使用以下脚本执行的操作:
library(colorRamps)
col5 <- colorRampPalette(c('blue', 'gray96', 'red')) #create color ramp starting from blue to red
color_levels=20 #the number of colors to use
max_absolute_value=0.4 #what is the maximum absolute value of raster?
plot(img, col=col5(n=color_levels), breaks=seq(-max_absolute_value,max_absolute_value,length.out=color_levels+1) , axes=FALSE)
使用来自 here 的数据,这是一个示例输出和实际脚本:
library(raster)
library(colorRamps)
mask_data=shapefile("D:/temp/so/Main_Hawaiian_Islands_simple3.shp")
img=raster("D:/temp/so/PPT_wet_minus_dry.tif")
col5 <- colorRampPalette(c('blue', 'gray96', 'red')) #create color ramp starting from blue to red
color_levels=10 #the number of colors to use
max_absolute_value=max(abs(c(cellStats(img, min), cellStats(img, max)))) #what is the maximum absolute value of raster?
color_sequence=seq(-max_absolute_value,max_absolute_value,length.out=color_levels+1)
plot(img, col=col5(n=color_levels), breaks=color_sequence, axes=FALSE)
plot(mask_data, add=T)
n_in_class=hist(img, breaks=color_sequence, plot=F)$counts>0
col_to_include=min(which(n_in_class==T)):max(which(n_in_class==T))
breaks_to_include=min(which(n_in_class==T)):(max(which(n_in_class==T))+1)
plot(img, col=col5(n=color_levels)[col_to_include], breaks=color_sequence[breaks_to_include] , axes=FALSE)
plot(mask_data, add=T)