rasterVis

Table of Contents

Introduction

The raster package defines classes and methods for spatial raster data access and manipulation. The rasterVis package complements raster providing a set of methods for enhanced visualization and interaction.

The stable release of rasterVis can be found at CRAN. The development version is at GitHub.

Let's show some of its functionalities with some examples, using data from the CMSAF project as described here (download data).

library(raster)
library(rasterVis)

##change to your folder...
old <- setwd('~/Datos/CMSAF')
listFich <- dir(pattern='2008')
stackSIS <- stack(listFich)
stackSIS <- stackSIS*24 ##from irradiance (W/m2) to irradiation Wh/m2
setwd(old)

idx <- seq(as.Date('2008-01-15'), as.Date('2008-12-15'), 'month')

SISmm <- setZ(stackSIS, idx)
names(SISmm) <- month.abb

Level plots

The first step is to display the content with a levelplot:

levelplot(SISmm)

levelplot.png

If only one layer is chosen, this method displays a marginal plot of a function across each coordinate:

levelplot(SISmm, layers=1, FUN.margin=median, contour=TRUE)

levelplot_layer1.png

The result of this call is a trellis object. The latticeExtra package provides the layer function to add contents. For example, let's add the administrative borders. This information is available at the GADM service.

library(maptools)
proj <- CRS('+proj=longlat +ellps=WGS84')
##Change to your folder
mapaSHP <- readShapeLines('~/Datos/ESP_adm/ESP_adm2.shp', proj4string=proj)

p <- levelplot(SISmm, layers=1, FUN.margin=median)
p + layer(sp.lines(mapaSHP, lwd=0.8, col='darkgray'))

levelplot_layer_borders.png

Log scale

The zscaleLog argument controls whether the object will be log transformed before being passed to the panel function. Defaults to ‘NULL’, in which case the Raster* is not transformed. Other possible values are any number that works as a base for taking logarithm, ‘TRUE’ (which is equivalent to 10), and ‘"e"’ (for the natural logarithm). As a side effect, the colorkey is labeled differently.

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
levelplot(r^2, zscaleLog=TRUE, contour=TRUE)

Themes

The previous plots used the default theme of rasterVis, rasterTheme. This theme defines a sequential palette with yellow, orange and red. There are three more themes in rasterVis: GrTheme (with a grey palette), BTCTheme (defined with the BTC palette of the hexbin package) and RdBuTheme (with a diverging palette with red and blue).

The irradiation of August is:

Aug <- raster(SISmm, 8)

and its overall mean is calculated with cellStats:

meanAug <- cellStats(Aug, mean)

The diverging palette is specially well suited to this data:

levelplot(Aug-meanAug, par.settings=RdBuTheme)

levelplotAug.png

Besides, it is easy to define a new theme with a different palette. For example, using a sequential palette from colorspace:

library(colorspace)
myTheme=rasterTheme(region=sequential_hcl(10, power=2.2))
levelplot(Aug, par.settings=myTheme, contour=TRUE)

levelplot_colorspace.png

or with the colour-blindness corrections from the dichromat package:

library(dichromat)
myTheme <- rasterTheme(region=dichromat(terrain.colors(15)))
levelplot(Aug, par.settings=myTheme)

levelplot_dichromat.png

Scatterplots and histograms

There are methods to show scatter plots and hexbin plots of the layers and coordinates of a Raster object:

##Relation between the January & February versus July radiation for four
##differents longitude regions.
xyplot(Jan+Feb~Jul|cut(x, 4), data=SISmm, auto.key=list(space='right'))

xyplot_formula.png

##Faster with hexbinplot
hexbinplot(Jan~Jul|cut(x, 6), data=SISmm)

hexbinplot_formula.png

…a method for scatter plot matrices:

splom(SISmm)

splom.png

..and methods for histograms, box-and-whisker and violin plots or density estimates:

histogram(SISmm)

histogram.png

densityplot(SISmm)

density.png

bwplot(SISmm)

bwplot.png

These methods accept a FUN argument to be applied to the z slot of the Raster object. The result of this function is used as the grouping variable of the plot:

histogram(SISmm, FUN=as.yearqtr)

histogram_FUN.png

Space-time plots

The z slot of this Raster object stores a time index. This 3D space-time Raster object can be displayed with a hovmoller diagram.

The hovmoller method uses the function xyLayer, which creates a RasterLayer from a function of the coordinates.

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
dirXY <-xyLayer(r, sqrt(x^2 + y^2))
dirXY

For example, the next code builds a hovmoller diagram showing the time evolution of the mean value along the latitude (data available at ftp://ftp.wiley.com/public/sci_tech_med/spatio_temporal_data/):

library(zoo)

url <- "~/Datos/Cressie/"
sst.dat = read.table(paste(url, "SST011970_032003.dat", sep=''), header = FALSE) 
sst.ll = read.table(paste(url, "SSTlonlat.dat", sep=''), header = FALSE)

spSST <- SpatialPointsDataFrame(sst.ll, sst.dat)
gridded(spSST) <- TRUE
proj4string(spSST) = "+proj=longlat +datum=WGS84"
SST <- brick(spSST)

idx <- seq(as.Date('1970-01-01'), as.Date('2003-03-01'), by='month')
idx <- as.yearmon(idx)
SST <- setZ(SST, idx)
names(SST) <- as.character(idx)
hovmoller(SST, contour=FALSE, panel=panel.levelplot.raster,
          yscale.components=yscale.raster.subticks,
          interpolate=TRUE, par.settings=RdBuTheme)

hovmoller.png

The horizonplot and xyplot methods also are useful for the space-time Raster objects:

horizonplot(SST)

horizon.png

Vector field plots

The function terrain from raster provides the vector field (gradient) from a scalar field stored in a RasterLayer object. The magnitude (slope) and direction (aspect) of the vector field is usually displayed with a set of arrows (e.g. quiver in Matlab).

rasterVis includes a method, vectorplot, to calculate and display this vector field.

proj <- CRS('+proj=longlat +datum=WGS84')
df <- expand.grid(x=seq(-2, 2, .01), y=seq(-2, 2, .01))

df$z <- with(df, (3*x^2 + y)*exp(-x^2-y^2))
r <- rasterFromXYZ(df, crs=proj)
vectorplot(r, par.settings=RdBuTheme())

vectorplot.png

If the Raster* object passed to vectorplot is a vector field (isField=TRUE), the terrain calculation is skipped.

An alternative method to display a vector field plots streamlines along the field lines. Streamlines, a family of curves that are tangent to the vector field, show the direction an element (droplet) will follow under the effect of the field. streamplot displays streamlines with a procedure inspired by the FROLIC algorithm: for each point (droplet) of a jittered regular grid, a short streamline portion (streamlet) is calculated by integrating the underlying vector field at that point. The main color of each streamlet indicates local vector magnitude (slope). Besides, streamlets are composed of points whose sizes, positions and color degradation encode the local vector direction (aspect).

streamplot(r)

streamplot.png

streamplot accepts two arguments (droplets and streamlets) to control the number of droplets, the length of the streamlets and the streamlet calculation step. The streamlet colour palette and the panel background color are defined with an specific theme for streamplot, streamTheme. The default options can be changed easily:

df$z <- with(df, sqrt(x^2 + y^2))
df$phi <- with(df, atan2(-y, x))
r2 <- rasterFromXYZ(df, crs=proj)

streamplot(r2, isField=TRUE, streamlet=list(L=30), droplet=list(pc=.3),
           par.settings=streamTheme(symbol=brewer.pal(n=5, name='Reds')))

streamplotReds.png

Interaction

This package includes two functions to interact with the trellis objects.

The identifyRaster method labels and returns points of a trellis graphic according to mouse clicks. It is commonly used after levelplot, although it can be also used after xyplot, hexbinplot or even splom:

levelplot(SISmm)

## Do not close the last graphical window.  Use the left button of the
## mouse to identify points and the right button to finish

chosen <- identifyRaster(SISmm, layer=3, values=TRUE)

The chooseRegion function provides a set of points (in the form of a SpatialPoints object) inside a region defined by several mouse clicks. Use the left button of the mouse to build a border with points, and the right button to finish. The points enclosed by the border will be highlighted and returned as a SpatialPoints object.

reg <- chooseRegion()

FAQs

View the Project on GitHub

Maintained by Oscar Perpiñán.