In R raster calculations, involving one or more raster layers, can be conducted in 3 different ways:

I. Using ‘Raster Algebra’

II. Using Higher Level Functions in the package raster

III. Using specialised functions from particular package

A brief review of these methods is presented below. Further details and examples can be found in TERN’s ‘Using Raster Data in R’ Tutorial.


I. 'Raster Algebra’ (also known as ‘Raster Maths’)

Used for small rasters and simple calculations.

Directly combines rasters using:



II. Raster Calculations using Higher Level Functions in the package ‘raster’

Directly adding, subtracting, etc. layers of big raster objects is not recommended. For big raster objects High Level Functions should be used instead.

On the other hand, in many cases, what can be achieved with Higher Level Functions, can also be accomplished with a more intuitive ‘raster-algebra’ notation. For example, the log10 of the cells in a raster can be computed as log10(r), which is simpler than using calc with this function. However, calcshould be faster when using large datasets and/or complex formulas.


II.A. For a Single Raster Object

The following ‘raster’ package functions can be used:


calc’: Uses a single Raster* object and a provided formula to calculate values for a new Raster* object.

The function type will depend of the input object type:

The object type returned depends on the input object type and function type:

The function can even be a statistical function, such as lm. See calc help for details.


stackApply’: Computes summary type layers for subsets of a RasterStack or RasterBrick.

To do so it uses a summary function applied on subsets of layers of a RasterStack or RasterBrick object. These subsets are defined by vector indices. If the number of indices in the vector is smaller than the number of layers the indices are recycled.

The summary function must return a single value; thus, the number of layers in the output Raster* equals then number of unique values in the indices vector. If the output raster contains a single layer it would be a RasterLayer object. If it contains multiple layers it would be a RasterBrick object.


overlay’: It is typically used with multiple Raster* objects, but it can also be used with a single Raster* object. However, it is often simpler to use ‘calc’ (or “Raster Algebra”) for single Raster* objects. See next section for its application to multi-layered raster objects.


II.B. For a Multiple Raster Objects

Rasters in operations including multiple rasters must have the same Dimensions, Projection, Extent, and Resolution.

The ‘raster’ package function to perform calculations involving multiple Raster* objects is ‘overlay’. It is more efficient that ‘Raster Maths’ when rasters are large and/or calculations complex and/or rasters stored in separate individual files. Do not confuse with the deprecated ‘overlay’ function in the ‘sp’ package. They are different.

Syntax: output.ras = overlay(input1.ras, input2.ras, fun=FunctionCalcName)

The raster objects used as inputs can be single-layered or multi-layered rasters (with the same number of layers). Both case produce similar results. In the multi-layered rasters case the function is applied to each pair of layers.



III. Raster Calculations using functions from other packages

R raster packages other than ‘raster’ contain functions that perform useful raster calculations. For example, the functions ‘normImage’ and ‘spectralIndices’ in the package ‘RStoolbox’ can standardise the raster and compute multi-spectral indices (e.g. NDVI, EVI, SAVI,…) respectively.





EXAMPLES


Below examples of raster calculations are shown below. These examples are taken from the “Effects of Cyclone Yasi on Green Cover at Mission Beach” tutorial. Putting the code snippets in context, by looking at a broader section of the R script, could be benefitial. Boxes with grey background contain code snippets, and boxes with white background containt code (text) outputs.



Single-Layer Raster Calculations

Re-scale raster cell values from a 0-255 scale to a 0-100 scale (to more appropriate represent the Green Cover Fraction in Percentages). To do so we divide by 255 and multiply by 100. This is a particular case of the general formula for Normalisation: ((X - Xmin) * 100) / (Xmax - Xmin). In our particular case, Xmin = 0 and Xmax = 255.


Example 1: Re-scale SPGC for the winter of year 2010

.

#==================================================================================
# Re-scale raster cell values
#==================================================================================
# The rasters cell values are in scale 0-255, so we re-scale them a 0-100 scale to represent the Green Cover Fraction in %.
# Generalized version: ((X - Xmin) * 100) / (Xmax - Xmin)
SPGC.StudyArea.2010q3.rl = SPGC.StudyArea.2010q3.rl * 100 / 255
SPGC.StudyArea.2010q3.rl

.

## class       : RasterLayer 
## dimensions  : 295, 206, 60770  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1484025, 1490205, -1996985, -1988135  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : SPGC_2010q3 
## values      : 39.21569, 77.2549  (min, max)

.

.

Example 2: Re-scale SPGC for the winter of year 2011

.

SPGC.StudyArea.2011q3.rl = SPGC.StudyArea.2011q3.rl * 100 / 255
SPGC.StudyArea.2011q3.rl

.

## class       : RasterLayer 
## dimensions  : 295, 206, 60770  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1484025, 1490205, -1996985, -1988135  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : SPGC_2011q3 
## values      : 39.21569, 76.86275  (min, max)

.

.

.

Multiple-Layer Raster Calculations

.

Example 1: Compute the change in Raw (in %) Seasonal Persistent Ground Cover (SPGC) fraction between the winters of 2011 and 2010

.

#----------------------------------------------------------------------------------
# Raw Change (in %)
#----------------------------------------------------------------------------------
# We use Raster Algegra to calculate Raw Change (i.e. change in %) in SPGC

SPGC.StudyArea.Diffq3.rl = SPGC.StudyArea.2011q3.rl - SPGC.StudyArea.2010q3.rl
SPGC.StudyArea.Diffq3.rl

.

## class       : RasterLayer 
## dimensions  : 295, 206, 60770  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1484025, 1490205, -1996985, -1988135  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -7.843137, 7.058824  (min, max)

.

.

Example 2: Compute the change in Standardised (in Standard Deviations) Seasonal Persistent Ground Cover (SPGC) fraction between the winters of 2011 and 2010

.

Method 1: Compute using Raster Algebra

.

# Method 1: Compute using Raster Algebra
# --------------------------------------
SPGC.StudyArea.StdDiffq3.rl = ((SPGC.StudyArea.2011q3.rl - cellStats(SPGC.StudyArea.2011q3.rl,mean)) / cellStats(SPGC.StudyArea.2011q3.rl,sd)) - 
                               ((SPGC.StudyArea.2010q3.rl - cellStats(SPGC.StudyArea.2010q3.rl,mean)) / cellStats(SPGC.StudyArea.2010q3.rl,sd))
SPGC.StudyArea.StdDiffq3.rl

.

## class       : RasterLayer 
## dimensions  : 295, 206, 60770  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1484025, 1490205, -1996985, -1988135  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -3.245601, 2.273073  (min, max)

.

.

Method 2: Compute using Higher Level Functions from the ‘raster’ package

.

# Method 2: Compute using Higher Level Functions from the Package `raster`
# ------------------------------------------------------------------------

calc.DiffStdRaster.f = function(rla, r1b)
{
    diff.std.rl = ((r1b - cellStats(r1b, mean)) / cellStats(r1b, sd)) - ((rla - cellStats(rla, mean)) / cellStats(rla, sd))
    return(diff.std.rl)
} # calc.DiffStdRaster.f = function(rla, r1b)
SPGC.StudyArea.StdDiffq3.rl.m2 = calc.DiffStdRaster.f(SPGC.StudyArea.2010q3.rl, SPGC.StudyArea.2011q3.rl)
SPGC.StudyArea.StdDiffq3.rl.m2

.

## class       : RasterLayer 
## dimensions  : 295, 206, 60770  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1484025, 1490205, -1996985, -1988135  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -3.245601, 2.273073  (min, max)

.

.

Method 3: Compute using ‘normImage’ function in ‘RStoolbox’ package

.

# Method 3: Compute using function `normImage` in package 'RStoolbox'
# -------------------------------------------------------------------
SPGC.StudyArea.StdDiffq3.rl.m3 = normImage(SPGC.StudyArea.2011q3.rl) - normImage(SPGC.StudyArea.2010q3.rl)
SPGC.StudyArea.StdDiffq3.rl.m3

.

## class       : RasterLayer 
## dimensions  : 295, 206, 60770  (nrow, ncol, ncell)
## resolution  : 30, 30  (x, y)
## extent      : 1484025, 1490205, -1996985, -1988135  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=aea +lat_1=-18 +lat_2=-36 +lat_0=0 +lon_0=132 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs 
## data source : in memory
## names       : layer 
## values      : -3.245601, 2.273073  (min, max)

.

.

Compare Raster resulting from the different Standarisation Methods

.

all.equal(SPGC.StudyArea.StdDiffq3.rl, SPGC.StudyArea.StdDiffq3.rl.m2)

.

## [1] TRUE

.

.

all.equal(SPGC.StudyArea.StdDiffq3.rl, SPGC.StudyArea.StdDiffq3.rl.m3)

.

## [1] TRUE

.

.

Remove unwanted objects

.

# Remove unnecessary (i.e. repeated) objects
#ls(pattern="SPGC.StudyArea.StdDiffq3")
rm(SPGC.StudyArea.StdDiffq3.rl.m2, SPGC.StudyArea.StdDiffq3.rl.m3)
#ls(pattern="SPGC.StudyArea.StdDiffq3")

.

.

.