MOLUSCE: transition potentials and a Monte Carlo cellular automaton

MOLUSCE is the land change plugin of QGIS, and for many ecologists the first land change model they meet: it takes two maps and a few drivers and, in a handful of dialogue boxes, gives a map of the future. In this part you rebuild its simulation in Terraflux on Sibuyan Island: the cross table of the observed change, the transition potential as MOLUSCE's weights of evidence (or a logistic regression), and its simulator, run from 2011 for two periods (2025 and 2039). You will check it on the calibration period, see where its logic differs from Dinamica EGO's and CLUE-S's, and learn to read the one behaviour that surprises its users most: the simulated change can be much smaller than the change it was asked for.

The idea#

The plugin#

MOLUSCE (Modules for Land Use Change Simulations) was written by NextGIS for Asia Air Survey and released as a QGIS plugin in 2013. It is organised as its modules: the inputs (the maps of two dates and the drivers, which it calls factors), the area change (the cross table of the two maps), the transition potential modelling (an artificial neural network, a logistic regression, weights of evidence, or multi-criteria evaluation), the cellular automaton simulation (a Monte Carlo simulator) and the validation (kappa statistics and an error budget). Its ease of use made it common in studies of land use change in Asia and Africa; Muhammad et al. (2022), for example, used it to project the land use of Linyi, China. Terraflux reproduces its potential (weights of evidence) and its simulator through landschaft, and offers any of its own potentials (a logistic regression here) in their place.

The quantities: the cross table#

MOLUSCE's quantity is the observed cross table of the period: the cells of each class at the first date that hold each class at the second. On Sibuyan from 1997 to 2011:

(made by change.crosstab). Each iteration of the simulator applies these counts once: an iteration is the whole period (14 years), not a year. From 2011, one iteration gives 2025, two give 2039.

The potential: total weights per transition#

MOLUSCE's weights of evidence cut each driver into categories by bounds you give (elevation at 50, 150, 300 and 600 m, slope at 4, 10, 16 and 22°, the distance to roads at 250 m to 4 km) and give each category a weight per transition, persistence included (forest to forest has weights too):

The weights of the distance to roads tell a familiar story: forest to grassland is most likely 500 m to 2 km from a road (weights of 0.85 and 0.72) and ruled out beyond 4 km (−33), while forest to forest becomes certain far from roads (4.9 beyond 4 km):

The simulator#

That last line is the essential difference from Dinamica EGO (whose rates always change their cells, placed by probability) and from CLUE-S (whose allocation iterates until the areas are met). In MOLUSCE the counts are a ceiling: the potential decides how much of it is reached.

MOLUSCE's other potentials#

The plugin offers four ways to make the potential. Its weights of evidence are reproduced exactly by molusce.total. A logistic regression has its equivalent in glm.fit and glm.predict (below). Its artificial neural network (a multi-layer perceptron) and multi-criteria evaluation (weighted overlays of standardised factors) have no node of their own: a neural network's potential computed elsewhere can be read as a stack of a layer per class, and a weighted overlay is a raster.calc expression per class put together by stack.of. Whatever its origin, a potential enters molusce.step the same way, and the simulator then behaves as described above.

When to use it#

Use this model to reproduce or extend a MOLUSCE study, or as a quick first model of a landscape with two maps and a few categorical drivers. Its potentials are easy to explain (a weight per category), and its quantities are the observed transitions. Its limits: one iteration is a whole period (no yearly steps), the quantities are a ceiling that the potential may not reach, the allocation is cell by cell (no patches), and its validation in the plugin rests on kappa. For yearly steps, patches and the quantities always met, use Dinamica EGO; for scenarios of areas, CLUE-S.

What you need#

• The two maps of the calibration period, and the map to start from (Sibuyan's of 1997 and 2011).
• Drivers, and the bounds of their categories (as in the plugin's Weights of evidence dialogue).
• The weights: fitted by MOLUSCE (the plugin's table, or landschaft's MolusceWoE in Python, which gives the same values) and written into molusce.total's settings; or another potential (a logistic regression).

Build the model#

format = "terraflux/1"

[[nodes]]
id = "start"
op = "raster.read"
params = { path = "sibuyan_landuse.tif", band = "2011" }

[[nodes]]
id = "drivers"
op = "raster.read"
params = { path = "sibuyan_drivers.tif" }

[[nodes]]
id = "weights"
op = "molusce.total"
inputs = { drivers = "drivers.raster" }
params = { classes = [1, 2, 3, 4, 5], drivers = ["elevation", "slope", "dist_roads"], bins = { elevation = [50, 150, 300, 600], slope = [4, 10, 16, 22], dist_roads = [250, 500, 1000, 2000, 4000] }, weights = [{ from = 1, to = 1, driver = "dist_roads", categories = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0], weights = [-3.203498, -2.363637, -1.334145, -0.114809, 1.482993, 4.906252] }, { from = 1, to = 3, driver = "dist_roads", categories = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0], weights = [-0.155531, 0.439672, 0.845118, 0.715742, -1.307703, -33.191493] }] }

[[nodes]]
id = "rng"
op = "rng.numpy"
inputs = { seed = 1 }

[[nodes]]
id = "simulation"
op = "control.repeat"
inputs = { times = 2, landscape = "start.raster" }

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

[[nodes.body.nodes]]
id = "allocate"
op = "molusce.step"
inputs = { landscape = "@landscape", total = "weights.total", rng = "rng.rng", step = "@step" }
params = { classes = [1, 2, 3, 4, 5], counts = [{ from = 1, to = 1, cells = 4160 }, { from = 1, to = 2, cells = 126 }, { from = 1, to = 3, cells = 320 }, { from = 1, to = 4, cells = 116 }, { from = 1, to = 5, cells = 1 }, { from = 2, to = 2, cells = 1140 }, { from = 2, to = 3, cells = 18 }, { from = 3, to = 2, cells = 77 }, { from = 3, to = 3, cells = 745 }, { from = 3, to = 4, cells = 16 }, { from = 3, to = 5, cells = 1 }, { from = 4, to = 2, cells = 61 }, { from = 4, to = 3, cells = 7 }, { from = 4, to = 4, cells = 156 }, { from = 5, to = 5, cells = 282 }] }

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

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

The graph shows two transitions' weights of one driver only, to stay readable; the model run for the figures has every transition and driver (model.toml).

The counts come from change.crosstab of the two observed maps; copy its table into counts (the persistence entries, such as { from = 1, to = 1 }, are ignored by the allocation but harmless).

A logistic regression instead#

The potential input takes a stack with a layer per class, such as glm.predict's suitability: each cell is then predicted the class of largest potential, and the rest is the same. Replace node 3 by a glm.fit of every class on the three drivers and a glm.predict, and connect its suitability to potential.

Read the results#

On the calibration period#

Run one iteration from 1997 and compare it with 2011 (validate.three_map, blocks of 1 to 8 cells):

With the weights of evidence the figure of merit is 0.14 at 250 m and 0.29 at 2 km; with the logistic regression half that. Dinamica EGO's model on the same period reaches 0.47 at 250 m (Validating a simulation), with thirteen drivers instead of three and patches instead of single cells; the comparison is about the whole model, not only the simulator.

Observed 2011

Simulated 2011

The quantities that were not reached#

The info table of the projection from 2011 says how many cells each transition changed. With the weights of evidence, the first iteration (2011 to 2025) changed forest to coconut 126 times, forest to grassland 320 and grassland to coconut 77, as asked, but forest to rice only 61 times of 116: only 61 forest cells were predicted to become rice. In the second iteration forest to rice changed no cell, and forest to coconut 70 of 126. With the logistic regression it is starker: forest to grassland changed 9 cells of 320 in the first iteration, none in the second, because the class of largest potential of almost every forest cell is forest.

This is MOLUSCE's design, not an error, but it changes the meaning of the result: the simulated map of 2039 is not "the trend of 1997-2011 continued", it is the part of that trend the potential agrees with. Always read info:

What to report#

The drivers and their categories; the weights (a table per transition, as the plugin shows them); the cross table used as counts; the number of iterations and the period each stands for; the cells each transition actually changed (info), against the counts; and the validation, with the figure of merit and the quantity and allocation disagreement rather than kappa alone (validate.kappa gives MOLUSCE's kappas; validate.error_budget its error budget).

Settings that change the answer#

• The categories' bounds. Coarse categories smooth the weights; fine ones fit noise, and categories without transitions get the weight of −33, which forbids the transition outright. Choose bounds with enough cells in each.
• The potential. As above, it decides how much of the counts is reached.
• The counts. They can be edited for a scenario (half the forest loss), as a ceiling.
• The number of iterations. Each is a whole period; there is no yearly step.

Common pitfalls#

• Reading the counts as the result. Compare info with the counts.
• Iterations taken for years. One iteration of a 14-year cross table is 14 years.
• Drivers with missing values. A cell missing a driver gets no total weight and never changes.
• Validating with kappa only. Use the figure of merit and the disagreement components (Validating a simulation).

Going further#

• Weights of evidence: Dinamica EGO's weights, with generalised ranges and significance.
• A demand per region, the next part.
• FUTURES and PLUS.
• The operators: molusce.total, molusce.step, rng.numpy, change.crosstab.

References#

• Muhammad, R., Zhang, W., Abbas, Z., Guo, F., & Gwiazdzinski, L. (2022). Spatiotemporal change analysis and prediction of future land use and land cover changes using QGIS MOLUSCE plugin and remote sensing big data: a case study of Linyi, China. Land, 11(3), 419. https://doi.org/10.3390/land11030419