Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Controlling legend and colors for raster values in R?

I'm trying to plot an ESRI Grid as a raster image of a surface. I've figured out how to make the plot, but not how to control R's color scale.

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

An example of the 'infile' would be something like this:

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 ...  
...  

The problem is that the 0 values in the data stretch the color axis beyond what's useful to me. See below:

enter image description here

Basically, I would like to tell R to stretch the color axis from 25-45, rather than 0-50. I know in Matlab I would use the command caxis. Does R have something similar?

like image 411
Sam Zipper Avatar asked Apr 04 '13 13:04

Sam Zipper


1 Answers

I used the dataset "volcano" to generate a raster object, and to make the code reproducible.

Ope option is to discretize raster values in classes, then, you can control the range of values in each class. Two examples below:

1- plot.

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

enter image description here

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

enter image description here

like image 88
Andre Silva Avatar answered Nov 14 '22 13:11

Andre Silva