An introductory tutorial by Agustin.Lobo@ictja.csic.es
20140406
updated: 20140512
file: /media/alobo/Iomega_HDD/RrasterGuide/rasterVisTutorial/rasterVisTut1_log.R
rasterVis is an R package for the visualization of raster layers authored by Oscar Perpinan
This introductory material that does not substitute the standard documentation provided with the package and the documents provided by Oscar Perpinan:
require(utils)
require(rgdal)
require(raster)
require(rasterVis)
Example 1. Probability field
Probability that a given pixel has got > 1kg/m2 of ash at a given volcanic eruption during the period of study
The probability field is a raster layer in GeoTif format
Download the file:
download.file("https://dl.dropboxusercontent.com/u/3180464/rprob520.tif", "rprob520.tif",
method = "curl")
Import the raster file to R:
The rgdal package (Bindings for GDAL) let us read and write raster files in different (virtually any) formats.
We first check the header information:
GDALinfo("rprob520.tif")
## rows 161
## columns 200
## bands 1
## lower left origin.x 418200
## lower left origin.y 6231355
## res.x 10000
## res.y 10000
## ysign -1
## oblique.x 0
## oblique.y 0
## driver GTiff
## projection +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs
## file rprob520.tif
## apparent band summary:
## GDType hasNoDataValue NoDataValue blockSize1 blockSize2
## 1 Float64 TRUE -1.7e+308 5 200
## apparent band statistics:
## Bmin Bmax Bmean Bsd
## 1 0.007692 1 0.3396 0.2335
## Metadata:
## AREA_OR_POINT=Area
And now create the raster layer object in R, using package raster (which uses rgdal for input/output)
rprob <- raster("rprob520.tif")
# CAUTION: by some reason, using a remote file delme <-
# raster('https://dl.dropboxusercontent.com/u/3180464/rprob520.tif') CRASHES
# RStudio
Exceptionally in R, entering the name of the raster layer object does not print the contents of the object but just some information on it:
rprob
## class : RasterLayer
## dimensions : 161, 200, 32200 (nrow, ncol, ncell)
## resolution : 10000, 10000 (x, y)
## extent : 418200, 2418200, 6231355, 7841355 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## data source : /media/alobo/Iomega_HDD/RrasterGuide/rasterVisTutorial/rprob520.tif
## names : rprob520
## values : 0.007692, 1 (min, max)
Beware though that you should not do this for large R objects in general.
Also exceptional in R context, raster layer objects do not hold their values in memeory. This is precisley
a unique feature of package raster aiming to do not collapse R with potentially huge objects.
We discus more on this and other features of package raster elsewhere.
The “standard” R raster plot would just be:
plot(rprob)
Alternative color tables:
par(mfrow = c(2, 2))
plot(rprob, col = rev(heat.colors(32)), zlim = c(0, 1))
plot(rprob, col = rev(terrain.colors(32)), zlim = c(0, 1))
plot(rprob, col = grey.colors(32, start = 0, end = 1, gamma = 1))
plot(rprob, col = rev(grey.colors(32, start = 0, end = 1, gamma = 1)))
See Wikipedia Gamma correction for background on the use of the gamma function for graphic display.
The simplest call to the function results into:
levelplot(rprob)
This plot implies no logarithmic transformation of the raster values (zscaleLog=NULL) and no level contours added (contour=FALSE).
We can modify this by:
levelplot(rprob, zscaleLog = TRUE, contour = TRUE)
and/or change the function for the margin distributions
levelplot(rprob, zscaleLog = NULL, contour = TRUE, FUN.margin = median)
In most cases, we prefer not to plot the margin distributions to gain space for the plot and add a title in the usual way in R graphics:
levelplot(rprob, contour = TRUE, margin = FALSE, main = "p(dep > 1kg/m2 per eruption event)")
The most minimalist raster plot with rasterVis would eliminate margins, labels and color key:
levelplot(rprob, margin = FALSE, scales = list(draw = FALSE), colorkey = FALSE)
We often want to control where the colors change:
miat = c(0, 0.25, 0.5, 0.75, 1)
levelplot(rprob, contour = TRUE, margin = FALSE, at = miat)
Note that the colors have changed at the values defined by miat, but values in the color key are independent.
We need the colorkey argument to control the color key. We start by a simple colorkey object in which we define where the colors should change in the color key:
myColorkey <- list(at = c(0, 0.5, 1))
levelplot(rprob, zscaleLog = NULL, contour = TRUE, margin = FALSE, at = miat,
colorkey = myColorkey)
which is kind of useless, as we most often want coincident breaks in the plot and in the color key,
but serves us to show that graphic and key can be stylized independently.
In a more useful example, ee define where the colors should change, which labels to write and where to write them:
myColorkey <- list(at=miat, ## where the colors change
labels=list(labels=miat, ##what to print
at=miat)) ##where to print
levelplot(rprob,zscaleLog=NULL,contour=TRUE, margin=FALSE,at=miat,colorkey=myColorkey)
This complexity provides flexibility, i.e.
myColorkey <- list(at = miat, labels = list(labels = c("Low", "Medium-low",
"Medium-high", "High"), at = miat + 0.125))
levelplot(rprob, zscaleLog = NULL, contour = TRUE, margin = FALSE, at = miat,
colorkey = myColorkey)
Color tables are controlled by argument par.settings, which requires a Theme. The default theme we have just used is named rasterTheme.
levelplot(rprob, contour = TRUE, margin = FALSE, at = (0:10)/10, par.settings = GrTheme)
levelplot(rprob, contour = TRUE, margin = FALSE, at = (0:10)/10, par.settings = BTCTheme)
levelplot(rprob, contour = TRUE, margin = FALSE, at = (0:10)/10, par.settings = RdBuTheme)
Themes require having selected a palette. In this example, we will use a standard 10-level rainbow palette in reversed order
Define function pal() to visualize the palette before actually making the Theme.
pal() is reproduced from HCL-Based Color Palettes in R
pal <- function(col, border = "light gray", ...) {
n <- length(col)
plot(0, 0, type = "n", xlim = c(0, 1), ylim = c(0, 1), axes = FALSE, xlab = "",
ylab = "", ...)
rect(0:(n - 1)/n, 0, 1:n/n, 1, col = col, border = border)
}
Visualize the palette
pal(rev(rainbow(n = 10)))
Make the Theme
rainbTheme10 <- rasterTheme(region = rev(rainbow(n = 10)))
Display with levelplot
levelplot(rprob, margin = FALSE, contour = TRUE, par.settings = rainbTheme10,
at = (0:10)/10, main = "p(dep > 1kg/m2 per eruption event)")
The best is saving the rasterVis graphic object for ulterior use and then either print() to the R graphic window or to a graphic “device” using the standard R method. For example:
rprobrv1 <- levelplot(rprob, margin = FALSE, contour = TRUE, par.settings = rainbTheme10,
at = (0:10)/10, main = "p(dep > 1kg/m2 per eruption event)")
Note saving to a graphic object implies that the graphic is not displayed on the graphic window. To do so:
print(rprobrv1)
Actually, the print() is not strictly required, just the name of the object would make the display. By default, entering the name of any R object shows its contents, which in the case of graphic objects implies displaying it.
In order to save the graphic to a graphic file, we first “open” the graphics device (in this case, a bmp file), then display it and finally close the device:
bmp("rprobrv1.bmp", height = 1024 * 0.707, width = 1024) #0.707 is a convenient aspect.ratio
rprobrv1
dev.off()
Important: as with the rest of standard graphics in R, the file is not actually written until dev.off()