Dinamica EGO on Sibuyan: the whole model

The first two parts of this series built the pieces of Dinamica EGO's classic land change model one by one: the weights of evidence that say where each transition is likely, and the Expander and Patcher that place the change in patches. In this part you put them together into the whole model of Sibuyan Island, from the maps of 1997 and 2011, and run it for fourteen yearly steps. You will write the model with one command (terraflux new dinamica), then read it node by node: what each node does, which output feeds which input and why, and above all how the loop carries the map from one year to the next, recomputing the probabilities on the map of each year. At the end you will have fourteen simulated maps, know how to read them and the tables beside them, how to make a projection beyond the last observed date, and what to check before trusting the result.

The idea#

A model in three pieces, run in a loop#

Dinamica EGO's classic model (Soares-Filho et al., 2002) answers the three questions of every land change model with three pieces fitted on two maps of the past:

These three are fitted once, before the simulation. The simulation is then a loop: at every step, the probabilities of the transitions are computed on the current map (woe.probability) and the change is allocated on it (dinamica.step); the new map is the start of the next step. Fourteen steps take the map of 1997 to a simulated map of 2011, which can be compared with the map that was really observed in 2011.

Why the probabilities are recomputed inside the loop#

The probability of a transition exists only in cells of its origin class: forest to grassland in forest cells. After a step, some forest is grassland, so the map of the probabilities computed on 1997 is wrong for 1998: it gives a probability of becoming grassland to cells that already are. Computing them inside the loop, on the map of the step, keeps them right. The weights themselves (the effect of each range of each driver) do not change: the model assumes that the processes of 1997 to 2011 go on. When a driver itself changes with the map (the distance to the nearest grassland, the share of farmland around a cell), it must be recomputed in the loop too: that is the subject of Dynamic drivers.

Sibuyan Island#

Sibuyan is a mountainous island of the Philippines, mostly forested, with coconut plantations, grassland and rice on its coastal plains. Its land use maps of 1997 and 2011 on cells of 250 m (108 × 128 cells, 7,226 of them on the island) and thirteen drivers are the classic data of CLUE-S (Verburg et al., 2002) and come with Terraflux (see What land change models are for the maps and their change). Between the two dates the forest fell from 65% to 58% of the island, mostly to grassland and coconut.

When to use it#

The whole model suits a landscape where change is a set of transitions with steady rates and patch-shaped outcomes: a deforestation frontier, the expansion of plantations or pastures, regrowth of abandoned land. It is the model of hundreds of studies of tropical deforestation since Soares-Filho et al. (2002), from the Amazon-wide scenarios of Soares-Filho et al. (2006) to studies of single municipalities.

Its assumptions are worth keeping in mind:

• the rates of the calibration period go on (a Markov chain: stationary, each step like the last);
• the drivers' effects go on, and each driver brings evidence the others do not;
• the shape of the change goes on (the same share of expansions, the same patch sizes).

When the amount of each class comes from elsewhere (a scenario of areas, an economic model), a demand-driven model fits better (CLUE-S); when it comes from what people decide, agents. The rates can also be given by hand (markov.matrix) for scenarios of faster or slower change, keeping the rest of the model.

What you need#

• Two maps of land use of the same area on the same grid: the start and the end of the calibration period, with whole class codes and a no-data value outside the area. Here sibuyan_landuse.tif, two bands named 1997 and 2011.
• The drivers: a raster with a band per driver, every band named, on the same grid (sibuyan_drivers.tif, 13 bands: distances to the sea, roads, towns and streams in metres, elevation, slope, aspect, population density, and five soil and erosion maps of 0 and 1).
• The number of time steps between the maps: 14 for yearly steps from 1997 to 2011.

Both files are in Terraflux's test data (tests/data/inputs) and can be downloaded from the first part's data. The application's template Dinamica EGO builds the same model on maps drawn by chance, to try it without data.

Build the model#

With one command#

terraflux new dinamica --t0 sibuyan_landuse.tif:1997 --t1 sibuyan_landuse.tif:2011 \
    --drivers sibuyan_drivers.tif --steps 14 -o sibuyan.toml
terraflux check sibuyan.toml
terraflux run sibuyan.toml --out sibuyan

Dinamica EGO's classic model: Markov rates, weights of evidence, Expander and Patcher

Usage: terraflux new dinamica [OPTIONS] --t0 <T0> --t1 <T1> --drivers <DRIVERS> --steps <STEPS> --output <OUTPUT>

Options:
      --t0 <T0>                      The map at the start of the calibration period: `file.tif` or `file.tif:band` (a name or a 0-based index)
      --threads <THREADS>            Threads to use (default and most: a quarter of the cores left free; never every core)
      --t1 <T1>                      The map at its end
      --drivers <DRIVERS>            The drivers: a raster with one band per variable
      --steps <STEPS>                Time steps in the calibration period
      --simulate <SIMULATE>          Steps to simulate (default: `steps`)
      --start <START>                The map the simulation starts from (default: t0)
      --increment <INCREMENT>        The weights' buffer increment (default: a hundredth of each variable's range)
      --expansion <EXPANSION>        The Expander's share of every transition (default: measured on the observed change)
      --isometry <ISOMETRY>          [default: 1.5]
      --prune-factor <PRUNE_FACTOR>  [default: 10]
      --seed <SEED>                  [default: 1]
  -o, --output <OUTPUT>              The model file to write (.toml or .json)
  -h, --help                         Print help

--t0 and --t1 name a file and a band (by name, or by its 0-based index). The model written is a plain text file in which every choice is a node with a note; the run takes a third of a second on Sibuyan and writes sibuyan/maps.tif (a band per step), sibuyan/steps.csv (what dinamica.step did at every step) and sibuyan/weights.csv (the weights report of part 1).

The inputs#

The model starts with four inputs, the values a user is expected to change (in the application they are the column at the left of the canvas; from the command line, --set name=value replaces one):

Inside the model they are referred to as $t0, $t1, $drivers and $seed.

Step 1: the three fitted pieces#

format = "terraflux/1"

[inputs.t0]
kind = "map"
path = "sibuyan_landuse.tif"
band = "1997"

[inputs.t1]
kind = "map"
path = "sibuyan_landuse.tif"
band = "2011"

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

[[nodes]]
id = "markov"
op = "markov.fit"
inputs = { t0 = "$t0", t1 = "$t1", steps = 14 }
note = "Transition rates per step: the 14th root of the matrix of the calibration period."

[[nodes]]
id = "woe"
op = "woe.fit"
inputs = { t0 = "$t0", t1 = "$t1", drivers = "$drivers", markov = "markov.markov" }
params = { ranges = { dist_sea = 250, dist_roads = 250, dist_urban = 250, dist_streams = 250, elevation = 50, slope = 2, aspect = 30, pop_density = 50 }, categorical = ["diorite", "ultramafic", "sediments", "no_erosion", "moderate_erosion"] }
note = "Weights of evidence of every transition with a rate."

[[nodes]]
id = "patches"
op = "patches.observe"
inputs = { t0 = "$t0", t1 = "$t1", markov = "markov.markov" }
params = { isometry = 1.5 }
note = "The Expander's share and the patch sizes measured on the observed change."

1. markov (markov.fit) reads t0 and t1 and steps = 14. It cross-tabulates the two maps (how many cells of each class in 1997 hold each class in 2011), turns the counts into the probabilities of the period, and takes the 14^th root of that matrix: the rates of one year whose fourteen repetitions give the period's change (see Markov chains). Its output markov carries both. On Sibuyan it has ten transitions with a positive rate.
2. woe (woe.fit) reads the two maps, the drivers and the chain. The connection from markov.markov to its markov input is what restricts the weights to the transitions the chain will ask for; without it every observed transition would be fitted. Its settings, ranges and categorical, are those chosen in part 1. terraflux new dinamica writes the node without them (every continuous driver cut into buffers of a hundredth of its range, --increment another width); add them in the Node pane or the file, as here. Its outputs: model (the weights, used by the loop) and table (the report).
3. patches (patches.observe) measures, for every transition of the chain, the Expander's share and the mean and variance of the patches' sizes on the observed change (part 2), and gives them with the isometry set (1.5). Its markov input, again from the chain, says which transitions to measure.

Why connect the chain to both? Because the three pieces must agree on the list of transitions: the allocation asks, for every transition with a rate, for its probability layer and its patch settings, and stops if one is missing. Feeding the same chain to all three guarantees they match.

Step 2: the loop#

format = "terraflux/1"

[inputs.t0]
kind = "map"
path = "sibuyan_landuse.tif"
band = "1997"

[inputs.t1]
kind = "map"
path = "sibuyan_landuse.tif"
band = "2011"

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

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

[[nodes]]
id = "simulation"
op = "control.repeat"
inputs = { times = 14, landscape = "$t0" }
note = "One iteration per time step, the map carried from step to step."

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

[[nodes.body.nodes]]
id = "probability"
op = "woe.probability"
inputs = { model = "woe.model", landscape = "@landscape", drivers = "$drivers" }
note = "The transition probabilities of the current map."

[[nodes.body.nodes]]
id = "allocate"
op = "dinamica.step"
inputs = { landscape = "@landscape", probability = "probability.probability", rates = "markov.markov", patches = "patches.patches", step = "@step", seed = "$seed" }
params = { prune_factor = 10.0, window = [3, 3] }
note = "Expander and Patcher: each transition changes its rate times the cells of its origin class."

The whole model, as written and run for the figures, is dinamica.toml. The loop node, simulation (control.repeat), has two inputs:

• times = 14: how many times the body runs;
• landscape = "$t0": a carried value, the map of 1997 at the first step. Any input of a loop other than times is carried: inside the body it is read as @landscape.

And two tables that define what moves between steps and what is kept:

• carry = { landscape = "allocate.landscape" }: at the end of every step, the carried value landscape takes the map the allocation gives; the next step's @landscape is that map. This is the line that makes the years follow each other.
• collect = { maps = "allocate.landscape", info = "allocate.info" }: the values to keep from every step. After the loop, simulation.maps is the series of the fourteen maps and simulation.info the fourteen tables (one CSV of fourteen rows when written). simulation.landscape, the carried value after the last step, is the map of 2011.

The body has two nodes:

1. probability (woe.probability) reads the weights from outside the loop (woe.model: fitted once), the map from inside (@landscape: this step's map) and the drivers ($drivers). Its output is the stack of a probability layer per transition, on this step's map.
2. allocate (dinamica.step) reads the same @landscape, the probabilities just computed, the rates (markov.markov) and the patch settings (patches.patches) from outside, the step number (@step, 1 to 14, a value every loop gives its body) and the seed ($seed). It gives the map after the step (landscape) and the table of what it did (info).

Why does step matter? The allocation's random numbers are drawn from the seed and the step, so every year gets its own draws while the whole run stays reproducible: the same seed gives the same fourteen maps on any computer.

Step 3: the outputs#

[outputs.maps]
from = "simulation.maps"      # maps.tif, a band per step

[outputs.steps]
from = "simulation.info"      # steps.csv, a row per step

[outputs.weights]
from = "woe.table"            # weights.csv, the weights report

The model's outputs say what terraflux run writes into the output folder. Any port can be one: add [outputs.patches] from = "patches.patches" to keep the patch settings as JSON, or a raster.write node for a map at a path of your choice. In the application, the maps of a loop show in the Map pane as they are computed, and the Charts pane draws the area of each class along the steps.

Read the results#

Fourteen years#

Enlarge

PlayYear

Enlarge

Move the slider and watch the north-west, where grassland spreads from its old patches into the forest, and the coasts, where coconut grows around older coconut and rice spreads along the south-western shore. Every year about 57 cells change (356 ha, 0.8% of the island); after fourteen years the shares are those of the observed map of 2011 to the tenth of a percent (forest 57.6%, coconut 19.4%, grassland 15.1%, rice 4.0%). That agreement is built in: the rates come from the same two maps. It says the chain was applied correctly, not that the model is good.

What the steps table says#

steps.csv has a row per step and, for every transition, four columns: requested_1_3 (the cells the rate asks for), expanded_1_3 and patched_1_3 (made by each function) and patches_1_3 (patches formed). Check that every row has expanded + patched = requested; on Sibuyan they do at every step. The requests fall slowly, as the origin classes shrink:

The simulated and the observed change#

Observed 2011

Simulated 2011

The simulated change sits in the same regions as the observed one: the grassland of the north-west, the coconut of the south-east, the rice of the south-western coast. It also invents some (grassland patches in the north-east that did not appear) and misses some (the scattered small changes of the east). Measured cell by cell, 461 of the 743 cells that really changed were simulated with the right class (a figure of merit of 0.47); within a window of 5 × 5 cells, the fuzzy similarity of the two changes is 0.74. What these numbers mean, and why the 93% of cells that agree is not one of them, is the subject of the next series, Validating a simulation. Add its two nodes to the model to compute them at every run:

[[nodes]]
id = "agreement"
op = "validate.agreement"
inputs = { start = "$t0", observed = "$t1", simulated = "simulation.landscape" }

[[nodes]]
id = "similarity"
op = "validate.similarity"
inputs = { start = "$t0", observed = "$t1", simulated = "simulation.landscape" }

Enlarge

What to report#

The calibration period and the time step; the classes; the transitions and their rates per step (the chain); the drivers, their ranges and the significant weights (part 1); the patch settings (part 2); the seed and the number of replicates; the validation against an observed map with its metrics at several resolutions (the next series). Publish the model file: it holds every one of these choices, readable.

Settings that change the answer#

The seed#

The allocation is stochastic. With another seed, the same model gives another landscape, the same amounts of each class in other places:

Enlarge

A single run is one draw among many possible landscapes. The cells that change in every run are consequences of the rates, the weights and the shapes; those that change in a few runs are where the allocation happened to fall. Run replicates and report what is common to them (Calibration and many runs): terraflux submit sibuyan.toml --replicates 50 runs fifty seeds on the cores of the computer.

The ranges of the weights#

The model as terraflux new dinamica writes it cuts every continuous driver into a hundred buffers (1,435 ranges in all on Sibuyan); with the increments chosen in part 1 there are 677. On the same seed the chosen ranges give a figure of merit of 0.47 against 0.41 for the automatic ones, and a similarity at the cell level of 0.60 against 0.56: noisy weights of thin ranges place the change less well. Choose the increments (or use --increment on the command line) rather than keep the default.

Projecting beyond the observed period#

The point of fitting a model on 1997 to 2011 is usually to ask what comes after. --start sets the map the simulation starts from and --simulate the steps to run, while the pieces are still fitted on --t0 and --t1:

terraflux new dinamica --t0 sibuyan_landuse.tif:1997 --t1 sibuyan_landuse.tif:2011 \
    --drivers sibuyan_drivers.tif --steps 14 --start sibuyan_landuse.tif:2011 --simulate 14 -o projection.toml

The model gets a fifth input, start (the map of 2011), which the loop carries instead of $t0, and times = 14. The rates, the weights and the patches are those of 1997 to 2011:

A projection is a scenario, "if the processes of 1997 to 2011 go on", not a forecast. Report it as such, with its replicates, and compare it with other scenarios (rates given by markov.matrix, a protected area as a weight set by hand, a road as a changed driver).

Other settings#

• prune_factor and window of dinamica.step (part 2); --prune-factor on the command line.
• The isometry and the Expander's share of patches.observe (--isometry, --expansion): calibrate them against the observed map, as Soares-Filho et al. (2013) do for Dinamica EGO's models (Calibration and many runs).
• The number of steps: yearly steps from a 14-year period are the 14^th root of the period's matrix; two steps of seven years are its square root. Choose the step that matches the questions (and the drivers that change).

Common pitfalls#

• Validating on the calibration period. A model fitted on 1997-2011 and compared with 2011 has only been checked for consistency. Its real test is a period it did not see: with three dates, fit on the first two and validate on the third (Validating a simulation).
• Classes the chain cannot change. A class with no observed transitions has no rates; it stays as it is in every simulation (Sibuyan's Other stays at 3.9%). That is the model's assumption, not a fact.
• Drivers in another order or other names. The weights are tied to the drivers' names. A drivers file for 2025 must have bands of the same names, or the weights no longer find their driver.
• Maps not on the same grid. The two dates and the drivers must share the grid (extent, cell size, rows and columns). terraflux check and the Wizard warn; terraflux resample drivers.tif --like landuse.tif puts a raster on another's grid.
• The loop reading the start map. Inside the body, read @landscape (see the warning above).
• One run reported as the result. Report replicates, and what is common to them.

Going further#

• Validating a simulation: the figure of merit, fuzzy similarity and the honest test on another period.
• Calibration and many runs: replicates, scenarios, and settings chosen against an observed map.
• Dynamic drivers: drivers recomputed in the loop.
• A demand per region: rates of their own for parts of the map.
• Maps larger than memory: the same model on a continent.
• The operators: markov.fit, woe.fit, woe.probability, patches.observe, dinamica.step, control.repeat.

References#

• Bonham-Carter, G. F. (1994). Geographic Information Systems for Geoscientists: Modelling with GIS. Pergamon.
• Soares-Filho, B. S., Coutinho Cerqueira, G., & Lopes Pennachin, C. (2002). DINAMICA: a stochastic cellular automata model designed to simulate the landscape dynamics in an Amazonian colonization frontier. Ecological Modelling, 154(3), 217–235. https://doi.org/10.1016/S0304-3800(02)00059-5
• Soares-Filho, B. S., Nepstad, D. C., Curran, L. M., Cerqueira, G. C., Garcia, R. A., Ramos, C. A., et al. (2006). Modelling conservation in the Amazon basin. Nature, 440(7083), 520–523. https://doi.org/10.1038/nature04389
• Soares-Filho, B., Rodrigues, H., & Follador, M. (2013). A hybrid analytical-heuristic method for calibrating land-use change models. Environmental Modelling & Software, 43, 80–87. https://doi.org/10.1016/j.envsoft.2013.01.010
• 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