Raster Calculations: Raster Algebra and Calculations using Functions

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:

Operators (for more information on R Operators see this link):

Arithmetic operators (+,’-,*,/,^,%%,%/%)
Relational operators (>, <, >=, <=, ==, !=)
Logical operators (!,&,&&,|,||)

Functions (algebraic (but ‘abs’), summary, and logical functions need > 1 layer):

Algebraic functions (e.g. abs, sum, prod)
Rounding functions (round, ceiling, floor, trunc)
Re-scaling functions (e.g. log, log10, sqrt, exp)
Trigonometric functions (e.g. sin, cos, tan, atan)
Summary functions (e.g. min, max, range)
Logical functions(any, all)

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:

For a uni-layer (RasterLayer) input, the function is typically a function that can take a single vector as input and return a vector of values of the same length (e.g. sqrt).
For a multi-layer (RasterStack or RasterBrick) input, the function should operate on a vector of values, with one vector for each cell.

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

When input is a RasterLayer, it returns a RasterLayer.
When input is a RasterStack or RasterBrick and the Function a Summary type Function (i.e. returns a single value), it returns a RasterLayer.
When input is a RasterStack or RasterBrick and the Function returns more than one number (e.g. fun=quantile), it returns a RasterBrick.

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. It can be beneficial to put the code snippets in context by looking at a broader section of the R script. Code snippets have a grey background, outputs a white background.

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

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

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

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 2: Compute using Higher Level Functions from the ‘raster’ package

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

Compare Raster resulting from the different Standarisation Methods

Remove unwanted objects