The raster calculator: scripts, a console and nodes of your own

In this tutorial you write one script that does the work of half a dozen map algebra nodes: from the land use of Sibuyan Island in 1997 and 2011 and a stack of its environmental variables, it finds where the forest stayed and where it was lost, scores the forest of 1997 as habitat for a species that avoids roads and high ground, picks out the best quarter of it, and gives the mean score as a number. You write it in the Calculator pane, which colours the script, completes the names of the layers, shows each result as you type and computes any line you ask it in a console. At the end you save the script as a node of your own, which you can add to any other model like an operator.

If map algebra is new to you, read Map algebra from scratch first: it explains what a raster is, what a local and a focal operation are, and how missing cells behave. This page builds on it.

The idea#

A raster.calc node computes new maps from other maps, cell by cell. Its expression used to be one formula (land == 1 and d.elevation < 1200), and a model that needed six formulas needed six nodes. A script is a short program of several lines: each line gives a name to a value, the later lines use the earlier names, and the lines that begin with out are what leaves the node. A single expression is still a script of one line, so every model written before keeps working.

The language is meant to read like Python and like numpy, the library most ecologists who program already know (Harris et al., 2020): = gives a name a value and == compares; and, or, not; a if c else b; ** for powers; # starts a comment. Like numpy it works on whole arrays at once: x * 2 doubles every cell, and x[x < 0] = 0 changes only the cells where the condition holds.

A script handles three kinds of values, and knowing which one you have explains most of what it does:

So mean(s) and s.mean() are different questions: the first is a map (the mean over the layers, cell by cell), the second a number (the mean of every cell of every layer). Reductions of both kinds leave missing cells out; the standard deviation is the sample one (dividing by \(n - 1\)) and the percentiles interpolate between the two nearest values, as numpy and R do by default.

When to use it#

Use a script when the steps belong together and are small: an index made of a few scores, a mask cleaned in two lines, a set of thresholds, a quick summary. One node keeps the reasoning in one place, and Terraflux computes all its outputs in a single pass over the cells, which is faster than a chain of nodes.

Use separate nodes when a step is worth seeing on the canvas by itself, when another node needs an intermediate result, or when the step is an operator of its own: a reclassification table reads better in raster.reclassify, a distance in distance.euclidean, a moving window larger than a few cells in the focal operators. When no expression can do it (a loop over cells that depends on the previous cell, a flood fill), write a plugin.

What you need#

The data of Map algebra from scratch, from Sibuyan Island in the Philippines at 250 m:

• landuse.tif, the land use in two bands named 1997 and 2011 (1 forest, 2 coconut, 3 grassland, 4 rice, 5 other); the sea has no value.
• drivers.tif, 13 bands named after what they hold: dist_sea, pop_density, diorite, ultramafic, sediments, no_erosion, moderate_erosion, elevation (m), slope, aspect, dist_roads (m), dist_urban and dist_streams.

Both are read as stacks, so the script reaches every band by its name. A band whose name is not a valid name in the language is given one: 1997 becomes _1997 (a name cannot start with a digit), NDVI-2020 becomes NDVI_2020, and a band without a name is b1, b2... The calculator shows both names.

Open the calculator#

The model has two inputs and one node:

format = "terraflux/1"

[inputs.lu]
kind = "stack"
path = "landuse.tif"

[inputs.d]
kind = "stack"
path = "drivers.tif"

[[nodes]]
id = "calc"
op = "raster.calc"
inputs = { lu = "$lu", d = "$d" }
params = { expression = "out years = sum(lu == 1)" }
note = "The forest's years, its loss, and a habitat index of the forest of 1997."

The pane has three parts. At the top, the node's inputs and what they are (lu: stack · 2, d: stack · 13; hover over one to see its layers). In the middle, the script, with line numbers and its words in colour: inputs, names you gave, layers, functions, numbers and comments each have their own. At the right of the script, a preview of every output. At the bottom, the console.

The inputs are known once what feeds them has run (press Apply and run once), or, for rasters read from files, from the files' headers before anything runs.

Write the script#

Type the script below line by line; the sections that follow explain it.

# lu: the land use of 1997 and 2011 (a stack); d: the drivers (a stack)
forest = lu == 1                        # a stack: the forest of each year
out years = sum(forest)                 # years as forest: 0, 1 or 2
out lost = forest._1997 and not forest._2011

# a score from 0 to 1 for elevation and for roads, and the index
si_elev = clamp((1200 - d.elevation) / 900, 0, 1)
si_road = clamp(d.dist_roads / 3000, 0, 1)
out hsi = sqrt(si_elev * si_road)
hsi[not forest._1997] = None            # the forest of 1997 only

out best = hsi >= hsi.percentile(75)    # its best quarter
out mean_hsi = hsi.mean()

Enlarge

Stacks: the forest of each year#

forest = lu == 1 compares every layer of lu with 1, so forest is a stack of two layers with the same names as lu's: the forest of 1997 and the forest of 2011, 1 where forest and 0 elsewhere. sum(forest) adds the layers at each cell: 2 where the cell was forest in both years, 1 in one of them, 0 in neither.

forest._1997 and not forest._2011 reads two layers by name. When you type forest. or d., the calculator offers the layers of that stack, with their names in the file beside them; Tab or Enter takes the one chosen, and Ctrl+Space asks for the names at any point. forest is not an input but a stack the script made, and it is completed all the same: a stack changed cell by cell (lu == 1, d * 2, abs(d)) keeps the layers of the one it came from, so forest. offers _1997 and _2011. A reduction over the layers (sum(forest)) is a map, with no layers to offer.

Enlarge

Names and conditions: the habitat index#

The two scores and the index are those of Map algebra from scratch: elevation scores 1 below 300 m and falls to 0 at 1,200 m, distance to roads scores 0 at the road and 1 beyond 3 km, and the index is their geometric mean. Here each score is a name, used on the next line, and nothing but hsi leaves the node.

hsi[not forest._1997] = None is a conditional assignment: only the cells where the condition holds are changed (here, the cells that were not forest in 1997 lose their value); the others keep the one they had. This is numpy's rule, and it replaces the nested a if c else b of a single expression. Several such lines apply in order, a later one over an earlier one.

Numbers: the best quarter and the mean#

hsi.percentile(75) is one number, the value that 75 % of the forest's cells are below (the cells without a value are left out). hsi >= hsi.percentile(75) then marks the best quarter of the forest of 1997, whatever the scores' range. out mean_hsi = hsi.mean() makes a number an output of the node: in the application it shows as a number, and another node can read it, for instance as a threshold in one of its settings.

The outputs#

Every out is an output port of the node, named as in the script: calc.years, calc.lost, calc.hsi, calc.best and calc.mean_hsi. out name = value gives the name and makes it an output; 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, so hsi leaves the node without the cells that line 10 took away.

The preview beside the script is computed about a third of a second after you stop typing, on a sample of up to 160,000 cells, so it stays quick on large maps; the node itself computes every cell when the model runs.

To look at one part closely (a valley, the edge of a reserve), drag a rectangle on any map of the preview: that part alone is computed, every cell of it when it has 160,000 cells or fewer, and a line over the outputs says which columns and rows it shows. Whole map, beside that line, or a double-click on a map goes back to all of it.

Errors, where they are#

A mistake is underlined where it is and told under the script with its line: a name that is neither an input nor given before (with the one you probably meant: hsii → hsi), a layer the stack does not have, a bracket not closed, = where a comparison was meant. Until it is fixed, the preview shows the outputs as they last were; a script given to the node with a mistake in it makes the run stop at the node, with the same message.

Give it to the node and run#

The script typed in the pane is not the node's until you give it: Ctrl+Enter (or Apply and run) gives it and runs the model up to the node; Apply gives it without running, as does leaving the editor. A dot beside the buttons marks a script not given yet.

The console#

Under the script, the console computes one line in the script's scope, on every cell, without changing the node: all the inputs and every name the script gives are there. Type a line and press Enter:

❯ hsi.percentile(75)
= 0.814519
❯ forest._1997.sum() * 0.0625
= 295.1875
❯ lost.sum() / forest._1997.sum()
= 0.119204

• a number is answered at once: hsi.mean(), forest._1997.sum() * 0.0625 (the forest's area in km², with cells of 6.25 ha), d.elevation.max();
• a map or a stack is answered with its least, mean and largest values, its missing cells and its image: d.dist_roads < 500, forest * 2;
• help lists the language, clear empties the console, and Up brings back the lines typed before.

Use it to check a value before you write it into the script: the share of forest that is lost, the range of a driver, a threshold.

Ask the assistant#

Under the inputs, a field marked ✦ passes a request to the application's own assistant. Type what to compute ("the forest lost between the two years, and its share of the forest of 1997") and press Write or Enter: the assistant writes the node's script, checked and tried on a sample of the cells before the node is given it, and says what each output gives. Explain asks it to explain the script there, line by line: what it reads, what each line does and what each name holds. Both give the node the script typed in the pane first, and the answer comes in the Assistant pane. A script the assistant wrote is a change like its others, which you review and keep or undo; read it as you would a colleague's before you trust its numbers.

Read the results#

In this run, 260 km² were forest in both years, and 35.2 km² of the 295.2 km² of forest of 1997 were gone by 2011: 11.9 %. The index of the forest of 1997 has a mean of 0.627 (mean_hsi, the number among the outputs); it is low on the high slopes of the centre, above 1200 m, and along the roads, and its best quarter, the cells at 0.81 or more, covers 73.8 km² in patches on the east and the north.

As in Map algebra from scratch, the index is a relative score made for this tutorial, not a published model of a species, and the "best quarter" is a choice: report the rule (the script) with the area it gives.

Save it as a node of your own#

The script now does one job well: from a land use stack and the drivers, it gives a habitat index and its best quarter. To use it again in another model, save it as a node of your own.

The dialog shows what the node takes and gives, as Terraflux found it, for you to change before saving:

Enlarge

• Name and What it does: how the node appears in your library (letters, digits and _).
• Inputs: each one's type (a map, a stack, a number...), the layers it reads, found in the script (lu must have _1997 and _2011, d must have elevation and dist_roads), and its value when none is given: an input that was a number in the model (a threshold, say) keeps it as its own value, and then need not be connected.
• Outputs: each one's type, from what it gave in the last run.

Save writes the node into your library as forest_hsi@1.toml, its first version, and the node of the model becomes that node: a calculator becomes a group of one node, its inputs the group's and each out an output. On the canvas its box names it forest_hsi@1, and its ports have the colours of their types.

Enlarge

Using it in another model#

Your nodes are listed in the Operators pane under My nodes (the search finds them by name and by what they do) and in the commands (Add forest_hsi@1). A click adds one after the node selected, as an operator is added. Hover over it to see what it takes and gives. The quick search finds them too: press / anywhere, type hsi, and your nodes come first, named my:forest_hsi; Enter adds the one chosen.

When you connect a stack to it, Terraflux checks the layers it reads: a land use stack whose bands are named otherwise (lu_1997, say) is told at once, in the node's pane and when the model runs, with the layers missing and those the stack has.

A node of yours is set like an operator. Its inputs that take a value (a number, a whole number, yes or no, a text) and that no other node feeds are its settings: in the Node pane they are under Settings, each field showing the node's own value (the one given in the dialog) until you type another, and emptying the field gives that value back. On the canvas they are folded away with the optional inputs, so the box shows only the data it takes. Connect another node's output to one, and it is an input again. A threshold saved with the node, say min_hsi = 0.6, is then a setting the next user can change without opening the node.

Versions#

Each save is a new version, and a version once saved never changes. A model keeps a copy of the version it uses, in the model file:

[[nodes]]
id = "calc"
op = "control.group"
inputs = { lu = "$lu", d = "$d" }
[nodes.uses]
name = "forest_hsi"
version = 1
description = "The forest's years, its loss, and a habitat index of the forest of 1997."
[nodes.uses.inputs.lu]
type = "stack"
bands = ["_1997", "_2011"]
[nodes.uses.inputs.d]
type = "stack"
bands = ["elevation", "dist_roads"]
[nodes.uses.outputs.hsi]
type = "map"
# ... the other outputs, then the body: the raster.calc node and its script

This has two consequences that matter for science. A model sent to a colleague runs without your library. And a model of last year gives the same result today, even if you have improved the node since: a newer version reaches a model only when you press Update to forest_hsi@2 in the node's pane (what the node is given is kept where the new version still takes it). The pane also tells when the node was changed inside the model since its version, and offers to save the change as a new version, or to make it a plain group again.

Libraries shared with your group#

Your library is a folder of files, one per version (forest_hsi@1.toml, forest_hsi@2.toml...). A lab that uses the same habitat index, the same cleaning of a land cover map or the same set of thresholds can keep such nodes in a shared library: a folder on a shared drive, or a repository that everyone clones.

A few things to know:

• The same name in two libraries. Your own library comes first, then the shared ones in their order in the settings; a node of a name is taken from the first library that has it.
• Taking a node out. A right click on it in the Operators pane offers to take it out of its library, every version. The models that use it, yours and others', keep their copies and still run.
• A colleague does not need the library to run your model, since the model keeps its copy of the version; a shared library is for building new models with the same nodes.
• The assistant uses them too: its tools (and those of terraflux mcp) list the nodes of your libraries, describe one by its name (my:forest_hsi) and add it to a model with its settings.

• Who saves in it. A shared library's folder may hold a library.toml: its name (the section's name in the Operators pane) and its writers, the user names of those who save in it and take nodes out of it:

  Copy

  name = "Habitat lab"
  writers = ["ana", "bruno"]
  

Anyone else uses its nodes but does not save in it: the dialog's Library shows it greyed (its reason on hover), and a node of it edited and saved again goes to your own library. A user name is the login of the computer (TERRAFLUX_USER when set). Without writers, anyone who can write the folder saves in it.

A remote library#

A library can also be fetched from elsewhere, and is then only read: a git repository (on GitHub, GitLab, Codeberg or a server of your lab), or a folder on the web. Terraflux keeps a copy in its settings folder and fetches it again each time it starts.

• A git repository is cloned, then pulled (only forward): Terraflux runs git, which must be installed. An address ending in .git, starting with git@ or ssh://, or a repository's page on GitHub, GitLab, Codeberg or Bitbucket is one; a private repository needs your git to have its credentials already, since Terraflux asks for none.
• A folder on the web is read by its index.txt, a node's file a line (forest_hsi@2.toml). Terraflux writes that index in each library it saves in, so a shared folder copied as it is to a web server (GitHub Pages, your lab's site) is a remote library. A version already fetched is not fetched again: versions never change.
• To add a node to a remote library, save it where the library is kept: in your clone of the repository (added as a shared library), then commit and push. The others get it the next time Terraflux starts, or with Update.

The folder of your own library is nodes in Terraflux's settings folder, unless the environment variable TERRAFLUX_NODES names another (see Environment and files).

The language in full#

Stacks and maps mix as you would expect: a stack with a number or a map acts on each layer; two stacks act layer by layer and must have the same number of layers.

Common pitfalls#

• = and ==. x = 1 gives x the value 1; x == 1 asks whether it is 1. In a condition, always ==.
• Python's syntax, and no other. ^, &&, ||, !, null, ln, if c then a else b and c ? a : b are refused, each with the form to write (**, and, or, not, None, log, a if c else b). A model written by an older Terraflux in those forms is rewritten in Python's when it is opened.
• No & or |. In numpy, a > 5 & b < 8 is read as a > (5 & b) < 8, a classic trap; the calculator has and and or, which mean what they say.
• mean(s) is not s.mean(). The first is a map, the mean of the layers at each cell; the second one number. mean(x) of a single map is refused and points to x.mean().
• min(a, b) and min(s). With several maps, min is the least of them at each cell, and a missing value in any gives a missing cell; with one stack, it is the least of its layers, the missing ones left out.
• Layer names. A band named 1997 is _1997 in a script. If the calculator marks a layer as unknown, hover over the input at the top of the pane to see the names it has.
• The preview is a sample. It is computed on at most 160,000 cells spread over the map; a neighbourhood's radius is scaled to the sample. The run computes every cell: look at the result, not only at the preview.
• A number of which cells? hsi.percentile(75) depends on the cells hsi has when it is computed: here only the forest of 1997, because line 10 took the others away before. Move that line after best, and the quarter is of every cell of the island, the farms and the grassland included.

Going further#

• The raster calculator: every part of the language with a playground that runs Terraflux's engine in your browser, recipes and pitfalls; and the vector calculator for the attributes of features.
• Map algebra from scratch: the ideas behind the operations, missing values, edges and buffers.
• Stacks and layers: build stacks from maps and folders of files.
• Groups and submodels: groups, and other model files used as one node.
• Plugins in WebAssembly: when a script is not enough.
• The operator reference: raster.calc, control.group.

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
• Tomlin, C. D. (1990). Geographic Information Systems and Cartographic Modeling. Prentice Hall.