如何在R中控制栅格数值的图例和颜色?

9

我正在尝试将ESRI Grid绘制为表面的光栅图像。我已经知道如何绘制,但不知道如何控制R的色彩比例。

# open necessary libraries
library("raster")
library("rgdal")
library("ncdf")

# goal: select an ESRI Grid ASCII file and plot it as an image.
infile <- file.choose("Results")
r <- raster(infile)

# read in metadata from ESRI output file, split up into relevant variables
info <- read.table(infile, nrows=6) 
NCOLS <- info[1,2]
NROWS <- info[2,2]
XLLCORNER <- info[3,2]
YLLCORNER <- info[4,2]
CELLSIZE <- info[5,2]
NODATA_VALUE <- info[6,2]
XURCORNER <- XLLCORNER+(NCOLS-1)*CELLSIZE
YURCORNER <- YLLCORNER+(NROWS-1)*CELLSIZE

# plot output data - whole model domain
pal <- colorRampPalette(c("purple","blue","cyan","green","yellow","red"))
par(mar = c(5,5,2,4))  # margins for each plot, with room for colorbars
par(pin=c(5,5))  # set size of plots (inches)
par(xaxs="i", yaxs="i")  # set up axes to fit data plotted
plot(r, xlim=c(XLLCORNER, XURCORNER), ylim=c(YLLCORNER, YURCORNER), ylab='UTM Zone 16 N Northing [m]', xlab='UTM Zone 16 N Easting [m]', col = pal(50))

“infile”的一个例子可能是这样的:”
NCOLS        262  
NROWS        257  
XLLCORNER     304055.000  
YLLCORNER    4792625.000  
CELLSIZE         10.000  
NODATA_VALUE    -9999.000  
42.4 42.6 42.2 0 42.2 42.8 40.5 40.5 42.5 42.5 42.5 42.9 43.0 ...  
42.5 42.5 42.5 0 0 43.3 42.7 43.0 40.5 42.5 42.5 42.4 41.9 ...  
42.2 42.7 41.9 42.9 0 0 43.7 44.0 42.4 42.5 42.5 43.3 43.2 ...  
42.5 42.5 42.5 42.5 0 0 41.9 40.5 42.4 42.5 42.4 42.4 40.5 ...  
41.9 42.9 40.5 43.3 40.5 0 41.9 42.8 42.4 42.4 42.5 42.5 42.5 ...  
...  

问题在于数据中的0值拉伸了颜色轴,超出了我需要的范围。请参见下面:

enter image description here

基本上,我想告诉R将颜色轴从25-45拉伸,而不是0-50。我知道在Matlab中我会使用命令caxis。R有类似的东西吗?


1
你能再解释一下吗?0 值应该映射到哪种颜色? - joran
如果将零替换为NaN会发生什么?我猜0应该映射为紫色(例如25也应该是紫色)?然后,您可以重新调整数据,使任何小于25的值变为0,任何大于45的值变为1。我不知道是否可以自动完成此缩放。 - Tobias
谢谢大家。理想情况下,零值可以是我选择的颜色。但是,如果它们被映射为紫色(最低颜色),那也是可以接受的。重新缩放是个好主意,如果我找不到更快的方法来调整色条,我可能会尝试类似的东西。 - Sam Zipper
1个回答

27

我使用数据集"volcano"生成了一个栅格对象,并使代码具有再现性。

一种选择是将栅格值离散化为类别,然后你可以控制每个类别中的值的范围。以下是两个示例:

1- 绘图.

library(raster)
library(RColorBrewer)

r = raster(volcano) #raster object

cuts=c(100,150,160,170,180,190,200) #set breaks
pal <- colorRampPalette(c("white","black"))

plot(r, breaks=cuts, col = pal(7)) #plot with defined breaks

在此输入图片描述

2- ggplot2

library(ggplot2)
library(raster)
library(RColorBrewer)

r = raster(volcano) #raster object
#preparing raster object to plot with geom_tile in ggplot2
r_points = rasterToPoints(r)
r_df = data.frame(r_points)
head(r_df) #breaks will be set to column "layer"
r_df$cuts=cut(r_df$layer,breaks=c(100,150,160,170,180,190,200)) #set breaks

ggplot(data=r_df) + 
  geom_tile(aes(x=x,y=y,fill=cuts)) + 
  scale_fill_brewer("Legend_title",type = "seq", palette = "Greys") +
  coord_equal() +
  theme_bw() +
  theme(panel.grid.major = element_blank()) +
  xlab("Longitude") + ylab("Latitude")

这里输入图片描述


1
很好,我很喜欢这个解决方案。谢谢。 - Sam Zipper

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接