rasterVis

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

levelplot(SISmm)

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)

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 here:

library(maptools)
proj <- CRS('+proj=latlon +ellps=WGS84')
##Change to your folder
mapaSHP <- readShapeLines('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'))

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

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)

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)

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

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

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'))
##Faster with hexbinplot
hexbinplot(Jan~Jul|cut(x, 6), data=SISmm)

…a method for scatter plot matrices:

splom(SISmm)

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

histogram(SISmm)
densityplot(SISmm)
bwplot(SISmm)

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)

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:

library(zoo)

url <- "ftp://ftp.wiley.com/public/sci_tech_med/spatio_temporal_data/"
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)
layerNames(SST) <- as.character(idx)
hovmoller(SST, contour=FALSE, panel=panel.levelplot.raster,
          interpolate=TRUE, par.settings=RdBuTheme)

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

horizonplot(SST)

The function slopeAspect 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. This method is not restricted to the RasterLayer class.

df <- expand.grid(x=seq(-2, 2, .1), y=seq(-2, 2, .1))
df$z <- with(df, (3*x^2 + y)*exp(-x^2-y^2))

r <- rasterFromXYZ(df)
projection(r) <- CRS("+proj=longlat +datum=WGS84")

vectorplot(r, par.settings=RdBuTheme)

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

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

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()
summary(reg)