A complete CLUE-S model of Plum Island

The first three parts of this series built CLUE-S's pieces one at a time: the competition between classes, the suitability of each class from logistic regressions, and the demand. In this part you assemble the whole model of the Plum Island Ecosystems with one command, terraflux new clue, and read it node by node: what each of its seven nodes does, how they connect, and how the loop carries two values from year to year, the map and its history, while R's random number generator keeps drawing. You will run the fourteen years from 1985 to 1999, watch the demand being met year by year, compare the simulated 1999 with the observed one, and move the two settings that change the answer most: the elasticities and the allocation's scale. At the end you will have a CLUE-S model you can adapt to your own maps.

The idea#

The model in one sentence#

Every year, every cell takes the land use whose suitability (plus a bonus for staying as it is, the elasticity, plus a number per class adjusted until the areas match) is highest, so that the map holds the cells the demand asks for that year; the next year starts from that map (Verburg et al., 2002). Three things are fitted or computed once, before the years:

and one thing runs every year, inside a loop: lulcc.step, the allocation of that year's demand.

Why carry the history, and the generator#

The history (lulcc.history) does nothing in this model: it is there because lulcc's conversion rules can depend on it ("forest can become built land only after 5 years as other land": part 5), and because a rule needs the history of the year it is applied in, so the history must travel with the map. The generator gives the small random jitter added to the suitabilities to break ties. It is R's own generator (Mersenne-Twister after set.seed(seed)): Terraflux draws exactly the numbers R would, which is why a run here gives lulcc's maps cell for cell (Moulds et al., 2015). The generator is one object that every step draws from in turn, so the draws of a year continue where the previous year's stopped, as in R.

When to use it#

Use the whole model when you have (or can assume) the amount of each land use year by year, and want its location: scenarios of regional land use, the consequences of a target ("no net loss of forest"), several classes growing at once. The first part discusses when CLUE-S is the right family; Choosing a land change model compares it with the others.

What you need#

• Observed maps of land use, one band per date, named by year (pie_landuse.tif: 1985, 1991, 1999; 434 × 497 cells of about 100 m; 1 forest, 2 built, 3 other; 255 no data).
• Drivers on the same grid (pie_drivers.tif: elevation, slope, dist_built_1985).
• Formulas: which drivers explain each class (from part 2).
• The files are on the first part's page. The application's template CLUE-S builds the same model on maps drawn by chance.

Build the model#

With one command#

terraflux new clue --observed pie_landuse.tif --times 1985,1991,1999 --drivers pie_drivers.tif \
  --formulas '{"1": ["elevation", "slope"], "2": ["elevation", "slope", "dist_built_1985"], "3": ["elevation", "slope"]}' \
  --elasticity 0.2,0.2,0.2 \
  --settings '{"jitter": 0.0002, "scale": 0.000001, "max_iter": 1000, "max_diff": 50, "mean_diff": 50}' \
  -o clue.toml
terraflux run clue.toml --out clue

These are the settings of lulcc's vignette for Plum Island. On Windows' PowerShell, write the JSON with single quotes as here; in cmd, escape the double quotes. Every option:

CLUE-S (lulcc's CluesModel): logistic regressions, a linear demand, CLUE allocation

Usage: terraflux new clue [OPTIONS] --observed <OBSERVED> --drivers <DRIVERS> --output <OUTPUT>

Options:
      --observed <OBSERVED>            The observed maps, one band per date (the demand is interpolated from them)
      --threads <THREADS>              Threads to use (default and most: a quarter of the cores left free; never every core)
      --times <TIMES>                  The date of each band of `observed`, e.g. `0,6,14`
      --drivers <DRIVERS>              The drivers: a raster with one band per variable
      --start <START>                  The map the simulation starts from and the regressions are fitted on (default: the first observed band)
      --formulas <FORMULAS>            The drivers of each class's regression as JSON, e.g. `{"1": ["slope"], "2": ["slope", "dist"]}`
      --elasticity <ELASTICITY>        CLUE: the elasticity of each class (one value, or a list in class order), e.g. `0.2,0.2,0.2`
      --order <ORDER>                  Ordered: the order the classes are served in, e.g. `2,1,3`
      --deterministic                  Ordered: take the most suitable cells instead of drawing them
      --settings <SETTINGS>            The allocation's settings as JSON (CLUE: jitter, scale, max_iter, max_diff, mean_diff)
      --allowed <ALLOWED>              The conversion matrix as JSON, in class order, e.g. `[[1,1,0],[0,1,0],[1,1,1]]`
      --neighbourhood <NEIGHBOURHOOD>  Neighbourhood rules as JSON, e.g. `{"rules": {"1": 0.3, "2": 0.2}}`
      --restrictions <RESTRICTIONS>    A map of restricted cells (non-zero keep their class)
      --history <HISTORY>              A map of the steps each cell has had its class
      --seed <SEED>                    [default: 1]
  -o, --output <OUTPUT>                
  -h, --help                           Print help

The run takes about 15 seconds on Plum Island (113,563 cells, 14 years; the slowest years iterate a few hundred times) and writes clue/maps.tif (a band per year, 1986 to 1999) and clue/steps.csv.

Step 1: suitability, demand, generator and history#

format = "terraflux/1"

[inputs.observed]
kind = "stack"
path = "pie_landuse.tif"

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

[inputs.start]
kind = "map"
path = "pie_landuse.tif"
band = "1985"

[inputs.seed]
kind = "value"
value = 1

[[nodes]]
id = "regression"
op = "glm.fit"
inputs = { landscape = "$start", drivers = "$drivers" }
params = { formulas = { 1 = ["elevation", "slope"], 2 = ["elevation", "slope", "dist_built_1985"], 3 = ["elevation", "slope"] } }
note = "Logistic regression of the presence of each class on its drivers (lulcc's glmModels)."

[[nodes]]
id = "suitability"
op = "glm.predict"
inputs = { model = "regression.model", drivers = "$drivers", landscape = "$start" }
note = "The suitability of each class at every cell; the drivers are static, so it is computed once."

[[nodes]]
id = "demand"
op = "demand.from_maps"
inputs = { maps = "$observed" }
params = { times = [1985, 1991, 1999] }
note = "The cells of each class at each step interpolated linearly from the observed maps."

[[nodes]]
id = "rng"
op = "rng.r"
inputs = { seed = "$seed" }
note = "R's generator after set.seed(seed), as lulcc draws."

[[nodes]]
id = "history"
op = "lulcc.history"
inputs = { landscape = "$start" }
note = "How long each cell has had its class (1 where not given)."

The model's four inputs are the observed maps (a stack: every band), the drivers (a stack), the start (one band of the observed maps, 1985) and the seed.

1. regression (glm.fit) fits, on the start map, one logistic regression per class: presence of the class (1) or not (0) against the drivers listed in formulas (a table of class to drivers). Its model output holds the coefficients (part 2 reads them).
2. suitability (glm.predict) turns the regressions into a stack with one layer per class (named 1, 2, 3), the probability of each class in each cell. It needs the drivers again, and the start map for its grid and its cells. Because Plum Island's drivers do not change during the simulation, it is computed once, outside the loop; a dynamic driver would move it inside (see Dynamic drivers).
3. demand (demand.from_maps) counts the cells of each class in each observed map and interpolates them linearly year by year (times: the date of each band). It has two outputs: demand, the table of cells per class and year, and steps, its number of years after the first (14), which the loop uses as its times. So the loop runs exactly as many years as the demand has.
4. rng (rng.r) makes R's generator from the seed.
5. history (lulcc.history) gives every cell a history of 1 (or reads a map of histories from its optional input).

Step 2: the loop#

format = "terraflux/1"

[inputs.observed]
kind = "stack"
path = "pie_landuse.tif"

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

[inputs.start]
kind = "map"
path = "pie_landuse.tif"
band = "1985"

[inputs.seed]
kind = "value"
value = 1

[[nodes]]
id = "regression"
op = "glm.fit"
inputs = { landscape = "$start", drivers = "$drivers" }
params = { formulas = { 1 = ["elevation", "slope"], 2 = ["elevation", "slope", "dist_built_1985"], 3 = ["elevation", "slope"] } }

[[nodes]]
id = "suitability"
op = "glm.predict"
inputs = { model = "regression.model", drivers = "$drivers", landscape = "$start" }

[[nodes]]
id = "demand"
op = "demand.from_maps"
inputs = { maps = "$observed" }
params = { times = [1985, 1991, 1999] }

[[nodes]]
id = "rng"
op = "rng.r"
inputs = { seed = "$seed" }

[[nodes]]
id = "history"
op = "lulcc.history"
inputs = { landscape = "$start" }

[[nodes]]
id = "simulation"
op = "control.repeat"
inputs = { times = "demand.steps", landscape = "$start", history = "history.history" }

[nodes.body]
carry = { landscape = "allocate.landscape", history = "allocate.history" }
collect = { maps = "allocate.landscape", info = "allocate.info" }

[[nodes.body.nodes]]
id = "allocate"
op = "lulcc.step"
inputs = { landscape = "@landscape", history = "@history", suitability = "suitability.suitability", demand = "demand.demand", step = "@step", rng = "rng.rng" }
params = { allocation = { method = "clue", jitter = 0.0002, scale = 0.000001, max_iter = 1000, max_diff = 50, mean_diff = 50 }, elasticity = [0.2, 0.2, 0.2] }

[outputs.maps]
from = "simulation.maps"

[outputs.steps]
from = "simulation.info"

The loop simulation (control.repeat) has three inputs: times, connected to demand.steps (a number that comes from a node, not typed), and two carried values, landscape (the start map) and history (node 5). Its body:

• carry = { landscape = "allocate.landscape", history = "allocate.history" }: after each year, both carried values take what the allocation gives, so next year's @landscape and @history are this year's results;
• collect = { maps = "allocate.landscape", info = "allocate.info" }: the map and the table of every year are kept (simulation.maps, simulation.info).

allocate (lulcc.step) is the year. Its inputs, and why each comes from where it does:

Its settings: allocation (method = "clue" and the five numbers of part 1), elasticity (one per class, in the order of the demand's classes, or { 2 = 0.9 } by class), and the rules of part 5 (allowed, neighbourhood; the input restrictions).

Read the results#

Fourteen years, and the demand met#

Enlarge

PlayYear

Every year the allocation met the demand within its tolerance of 50 cells (steps.csv: converged = 1, from 52 iterations in 1986 to a few hundred later). The map, though, tells a second story: the first year changed 13,023 cells (11% of the landscape) for a net change of a few hundred, and the following years about 450 to 600 each. In the first year every cell whose suitability favours another class by more than the elasticity of 0.2 switches: the map is "redrawn" to agree with the regressions. After that the map is consistent with them and only the demand's yearly increments move cells. This first-year reshuffle is the strongest sign that the elasticities are too low (below).

Against the observed map of 1999#

Observed 1999

Simulated 1999

Enlarge

The comparison of Validating a simulation (validate.agreement and validate.three_map, from 1985): 1,609 hits, 5,948 misses, 1,021 wrong hits and 16,892 false alarms, a figure of merit of 0.063. The simulated map agrees with 79.0% of the observed one, while persistence (the map of 1985 unchanged) agrees with 92.4%: the false alarms are the reshuffle of the first year. Measured at coarser blocks, the figure of merit rises:

This is lulcc's published example, run exactly; it is a demonstration of the method, not a calibrated model. The next section shows the first thing to change.

Replicates#

The jitter is tiny (0.0002 against suitabilities that differ by hundredths), so the seed matters little here: runs with seeds 2 to 5 differ from seed 1 in 27 to 33 cells of 1999, and 48 cells are built in some of the five runs but not all. CLUE-S's allocation is nearly deterministic; its uncertainty is in its inputs (the demand, the regressions, the elasticities), which scenarios explore better than seeds.

What to report#

The demand (its source, its table); the regressions (formulas, coefficients, their fit, part 2); the elasticities and any rules; the allocation's settings; whether every year converged (steps.csv); the validation on an observed map at several resolutions, against persistence.

Settings that change the answer#

The elasticities#

Enlarge

PlayElasticity (forest, built, other)

Enlarge

The observed change from 1985 to 1999 was 8,578 cells. With no elasticity, a third of the landscape (33,430 cells) ends in another class: the model is free to redraw everything its regressions disagree with. Raising the elasticity to 0.4 brings the change down to 8,267 cells, the observed amount, and the agreement with the map of 1999 from 79% to 87%. Higher elasticities change little more but make the allocation struggle: with 0.8, the iteration variables must climb far above the elasticity before a cell will move, and with scale = 0.000001 the 1,000 iterations are not enough in seven years of fourteen (converged = 0; Terraflux warns). Raise max_iter, or the scale, with high elasticities. Giving built land alone a high elasticity (0.9, last frame) keeps the towns from turning into forest, a realistic rule, and cuts the reshuffle by a quarter.

The figure of merit barely moves (0.058 to 0.065): the elasticities decide how much of the landscape is reshuffled, not whether the new built land lands on the right cells, which is the regressions' job. Choose them from what is known of each class (in CLUE-S's logic, near 1 for land uses that are costly to convert, such as built land or orchards, lower for land that changes easily, such as shifting farmland; Verburg et al., 2002) and check the gross change of the first year against the observed yearly change. Part 5 adds conversion rules for finer control.

The scale#

Enlarge

PlayScale of the iteration variables

scale × (the difference in cells) is how far an iteration variable moves in one pass. With 0.00001, a difference of 5,000 cells moves it by 0.05, about the gaps between suitabilities: the allocation converges in 5 passes. With 0.0001 the same difference moves it by 0.5, more than the whole range of the suitabilities, and the variables swing from one extreme to the other: in the last pass every cell of the landscape is forest, and the map of that year is useless. lulcc's default (0.0005) is far too large here; on maps of 100,000 cells and more, start from 0.000001 to 0.00001 and check the first year. Always read converged in steps.csv.

Common pitfalls#

• Ignoring the first year. A gross change far larger than the net change in the first year means the elasticities let the regressions redraw the map. Compare it with the observed yearly change.
• A start map the demand does not begin from. The first row of the demand must be the start map's counts: demand.from_maps guarantees it when the start is the first observed map.
• Regressions fitted on another map than the start. The regression and suitability nodes read $start; when projecting from a later map (--start of terraflux new clue), the regressions are fitted on it.
• converged = 0. The map of that year does not have the demanded areas; fix scale, max_iter or the tolerances.
• Treating lulcc's vignette as calibrated. Its settings show the method; calibrate yours (elasticities, formulas) against an observed map (Calibration and many runs).
• Classes of the map missing from the demand, or a stray code: the step stops (part 1).

Going further#

• Part 5: Rules: conversion matrices, neighbourhoods and restrictions.
• Validating a simulation and Calibration and many runs.
• A demand per region: CLUE-S with a demand for each part of the map.
• Dynamic drivers: the suitability recomputed in the loop.
• The operators: lulcc.step, glm.fit, glm.predict, demand.from_maps, rng.r, lulcc.history.

References#

• Moulds, S., Buytaert, W., & Mijic, A. (2015). An open and extensible framework for spatially explicit land use change modelling: the lulcc R package. Geoscientific Model Development, 8(10), 3215–3229. https://doi.org/10.5194/gmd-8-3215-2015
• Verburg, P. H., Soepboer, W., Veldkamp, A., Limpiada, R., Espaldon, V., & Mastura, S. S. A. (2002). Modeling the spatial dynamics of regional land use: the CLUE-S model. Environmental Management, 30(3), 391–405. https://doi.org/10.1007/s00267-002-2630-x