The raster calculator

The raster calculator, raster.calc, computes new maps from other maps, cell by cell, with a short script that reads like Python: forest = lu == 1, out lost = forest._1997 and not forest._2011. This page has a playground where you can write such scripts over real data and see every result at once, then the language part by part, recipes for the tasks that come up in most analyses, and the mistakes worth knowing before they cost you a day.

Try it#

The playground below computes over Sibuyan Island (Philippines) at 250 m, 128 × 108 cells: lu, its land use in 1997 and 2011 (a stack of two layers, lu._1997 and lu._2011: 1 forest, 2 coconut, 3 grassland, 4 rice, 5 other; the sea has no value), and d, seven of its drivers (a stack: d.elevation in m, d.slope in degrees, d.dist_roads, d.dist_urban, d.dist_streams and d.dist_sea in m, d.pop_density). Pick an example, or change the script: it runs as you type.

The idea#

A script is a short program over whole maps. Each line gives a name a value (forest = lu == 1), the later lines read the earlier names, and the lines that begin with out are the node's outputs. Like numpy, the library most ecologists who program know (Harris et al., 2020), it computes on whole arrays at once: x * 2 doubles every cell, and x[x < 0] = 0 changes only the cells where the condition holds. The operations are Tomlin's local and focal ones (Tomlin, 1990): a cell's value from the values of the same cell in other maps, or of the cells around it.

A script handles three kinds of values; knowing which one you have explains most of what it does.

The language#

Values and names#

• Numbers (3, 0.25, 1e6), pi, True and False (1 and 0), and None, a missing value.
• The node's inputs by the names of its ports: maps, stacks and numbers. A stack's layers are s.layer. A layer whose name is not a name of the language gets one: 1997 is _1997 (a name cannot begin with a digit), NDVI-2020 is NDVI_2020, a band without a name b1, b2...
• Names you give: name = value. A name is a value, not a formula: y = x, then x = x * 10, leaves y as x was. A name given in the script hides an input of the same name from there on.
• Lines: one statement a line, or several separated by ;; # begins a comment; a line goes on after an open bracket until it closes.

Arithmetic, comparisons, logic#

Every operation is made cell by cell. A map with a number applies the number to every cell (d.elevation / 1000); two maps pair their cells (lu._1997 * 10 + lu._2011).

Functions#

Changing some cells#

name[condition] = value changes the cells where the condition holds; the others keep their value. It is numpy's rule, and it reads better than a nested a if c else b:

x = d.slope
x[x > 30] = 30                    # slopes above 30 degrees capped
x[lu._2011 == 5] = None           # and "other" left out

Several such lines apply in order, a later one over an earlier one. A cell where the condition has no value (over the sea, lu._2011 == 5 has none) is not changed.

s.layer = value replaces a layer of a stack or adds one: s.changed = lu._1997 != lu._2011.

Missing values#

A cell without a value (the sea, a cloud, a mask) is None, and a value computed from a missing one is missing: None + 1, None > 3 and not None are None, and so is 0 / 0; 1 / 0 is infinite (inf), and sqrt of a negative number is missing. Three tools handle them:

• isnull(x), which is 1 where x has no value: isnull(lu._2011).sum() counts the cells of the sea;
• coalesce(x, y), which fills: coalesce(land, 0);
• x[condition] = None, which takes cells away: hsi[lu._1997 != 1] = None.

The numbers of a whole map (x.mean(), x.sum()) and the reductions of a stack's layers (mean(s)) leave the missing cells out; min(a, b) of several maps is missing where any of them is.

Stacks and their layers#

An operation on a stack is made on each layer, with a number or a map (lu == 1, d * 2), and layer by layer with another stack of as many layers. Its layers keep their names, so (lu == 1)._1997 is the forest of 1997.

At each cell, over the layers (the missing ones left out), a map:

sum(lu == 1) is in how many years a cell was forest. count(lu == 1) is not: it counts the layers with a value, 2 on all of the island.

Numbers of all the cells#

x.mean() x.sd() x.median() x.sum() x.count() x.min() x.max() x.percentile(p), with the brackets, give one number from every cell of x with a value (of every layer of a stack). The standard deviation is the sample one (dividing by \(n - 1\)), and percentiles interpolate between the two nearest values, as numpy and R do by default.

mean(s) and s.mean() answer different questions: the first is a map (the mean of the layers at each cell), the second one number (the mean of every cell of every layer). mean(x) of a single map is refused, and points to x.mean().

A number depends on the cells its map has when it is computed: in the habitat example, hsi.percentile(75) comes after hsi[lu._1997 != 1] = None, so it is the best quarter of the forest of 1997; move it before, and it is of the whole island.

Neighbourhoods#

nbsum(x, r), nbmean(x, r), nbmin(x, r), nbmax(x, r) and nbcount(x, r) read the square of radius r cells around each cell, the cell included (r = 1: 3 × 3 cells; r = 2: 5 × 5), the missing cells left out; nbcount counts those with a value. nbmean(lu._2011 == 1, 2) is the share of forest in the 5 × 5 cells around each cell, a moving window as in FRAGSTATS (McGarigal & Marks, 1995); nbmax(farm, 1) == 1 marks the cells next to farmland. r is a whole number or a number input. For larger windows, circles or weights, use the focal operators (raster.focal).

Since the missing cells are left out, a cell without a value next to cells with one gets a value: the sea along the coast has a share of forest. Take it away with share[isnull(lu._2011)] = None.

Outputs#

• out name = value gives a name and makes it an output of the node; out name alone makes an output of a name given before. An output takes the value its name has at the end of the script.
• An output is a map, a stack or a number, an output port of the node named as in the script: calc.years, calc.lost, calc.lost_km2.
• What each output is, is known before the model runs, from what the node's inputs are: a map made of maps is a map, one made of a stack a stack, x.mean() a number, mean(s) of a stack's layers a map. Connecting an output where it does not fit is refused when the model is checked: calc.lost_km2 (a number) where a map is taken, or calc.lost (a map) where a number is. One made of rasters read from files is a map or a stack (by their bands), fit for either, but not for a number.
• A script without out gives its last value as the output map: lu._2011 == 1 and d.elevation > 500 is a script too.

What is refused#

The calculator reads Python's syntax and no other. Each of these is refused before anything runs, with the form to write:

A name that is neither an input nor given before is told with the one you probably meant (d.elevaton: did you mean d.elevation?), and a layer the stack does not have with the layers it has.

Recipes#

In the application, the model file and the assistant#

[[nodes]]
id = "calc"
op = "raster.calc"
inputs = { lu = "landuse.raster", d = "drivers.raster", max_slope = 20 }
params = { expression = """
forest = lu == 1
out years = sum(forest)
out steep = lu._2011 == 1 and d.slope > max_slope
""" }

Common pitfalls#

• = and ==. x = 1 gives x the value 1; x == 1 asks whether it is 1.
• No & or |. Use and and or. In numpy, a > 5 & b < 8 is read as a > (5 & b) < 8, a classic trap.
• count(s) is not a count of true layers. It counts the layers with a value; in how many layers a condition holds is sum(s == 1).
• s.mean() of a stack is of every layer at once. (d - d.mean()) / d.sd() mixes metres, degrees and people; standardize each layer by itself: (d.elevation - d.elevation.mean()) / d.elevation.sd().
• Python's arithmetic, almost. round(2.5) is 3 (a half is rounded away from 0; Python gives 2), and % keeps the sign of the number divided (-7 % 3 is -1; Python gives 2).
• Missing values spread, through and and or too (see above). Look at the missing cells of an output: the maps of the playground draw them checkered.
• Neighbourhoods reach into missing cells. A window around a cell of the sea has land in it, so the sea gets a value; take it away with x[isnull(land)] = None.
• A number of which cells? A percentile or a mean depends on the cells the map has at that line (see Numbers of all the cells).
• Units. Areas of cells are those of the grid: cells of 250 m are 6.25 ha, 0.0625 km². A map in geographic coordinates (degrees) has cells of different areas: reproject it first (coordinate systems).

Going further#

• The vector calculator: the same language over features, with texts and geometries.
• Map algebra from scratch: local and focal operations, missing values, buffers and edges, built step by step.
• The raster calculator, a tutorial: the Calculator pane, its console, and a script saved as a node of your own, with its versions.
• Stacks and layers: stacks from maps and folders of files.
• The operator reference: raster.calc, raster.reclassify, raster.focal.

References#

• Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del Río, J. F., Wiebe, M., Peterson, P., … Oliphant, T. E. (2020). Array programming with NumPy. Nature, 585(7825), 357–362. https://doi.org/10.1038/s41586-020-2649-2
• McGarigal, K., & Marks, B. J. (1995). FRAGSTATS: spatial pattern analysis program for quantifying landscape structure (General Technical Report PNW-GTR-351). USDA Forest Service, Pacific Northwest Research Station. https://doi.org/10.2737/PNW-GTR-351
• Tomlin, C. D. (1990). Geographic Information Systems and Cartographic Modeling. Prentice Hall.