Raster algebra

Many generic functions that allow for simple and elegant raster algebra have been implemented for Raster* objects, including the normal algebraic operators such as +, -, *, /, logical operators such as >, >=, <, ==, !} and functions such as abs, round, ceiling, floor, trunc, sqrt, log, log10, exp, cos, sin, max, min, range, prod, sum, any, all. In these functions you can mix raster objects with numbers, as long as the first argument is a raster object.

# create an empty RasterLayer
library(raster)
r <- raster(ncol=10, nrow=10)
# assign values to cells
values(r) <- 1:ncell(r)
s <- r + 10
s <- sqrt(s)
s <- s * r + 5
r[] <- runif(ncell(r))
r <- round(r)
r <- r == 1

You can also use replacement functions:

s[r] <- -0.5
s[!r] <- 5
s[s == 5] <- 15

If you use multiple Raster* objects (in functions where this is relevant, such as range), these must have the same resolution and origin. The origin of a Raster* object is the point closest to (0, 0) that you could get if you moved from a corners of a Raster* object towards that point in steps of the x and y resolution. Normally these objects would also have the same extent, but if they do not, the returned object covers the spatial intersection of the objects used.

When you use multiple multi-layer objects with different numbers or layers, the ‘shorter’ objects are ‘recycled’. For example, if you multuply a 4-layer object (a1, a2, a3, a4) with a 2-layer object (b1, b2), the result is a four-layer object (a1b1, a2b2, a3b1, a3b2).

r <- raster(ncol=5, nrow=5)
r[] <- 1
s <- stack(r, r+1)
q <- stack(r, r+2, r+4, r+6)
x <- r + s + q
x
## class      : RasterBrick
## dimensions : 5, 5, 25, 4  (nrow, ncol, ncell, nlayers)
## resolution : 72, 36  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs
## source     : memory
## names      : layer.1, layer.2, layer.3, layer.4
## min values :       3,       6,       7,      10
## max values :       3,       6,       7,      10

Summary functions (min, max, mean, prod, sum, Median, cv, range, any, all) always return a RasterLayer object. Perhaps this is not obvious when using functions like min, sum or mean.

a <- mean(r,s,10)
b <- sum(r,s)
st <- stack(r, s, a, b)
sst <- sum(st)
sst
## class      : RasterLayer
## dimensions : 5, 5, 25  (nrow, ncol, ncell)
## resolution : 72, 36  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs
## source     : memory
## names      : layer
## values     : 11.5, 11.5  (min, max)

Use cellStats if instead of a RasterLayer you want a single number summarizing the cell values of each layer.

cellStats(st, 'sum')
## layer.1.1 layer.1.2 layer.2.1 layer.2.2   layer.3
##      25.0      25.0      50.0      87.5     100.0
cellStats(sst, 'sum')
## [1] 287.5