Chapter 26: Land Use and Cover Change Modeling¶

Part V — Domain Modeling: Land Use & Coastal Systems

Implemented by the DisSLUCC package (disslucc 0.5.0 on PyPI).

Learning Objectives¶

By the end of this chapter you will be able to:

• Say what land use and land cover are, and what a land use and cover change (LUCC) model does and does not tell you
• Describe the Demand/Potential/Allocation loop that LuccME established and DisSLUCC keeps, and say what each of its components computes, with its formula
• Tell a continuous model (a cell holds a fraction of each class) from a discrete one (a cell holds one class), and choose between them for a research question
• Build and run a continuous LUCC model in a script, and know the executor and TOML route to the same components
• Read the LuccME test labs as a map of the components, and say what the comparison of DisSLUCC with them shows — and what it does not

In [1]:

# Standard imports — add chapter-specific imports below
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Standard imports — add chapter-specific imports below import numpy as np import pandas as pd import matplotlib.pyplot as plt

Every model since Chapter 23 has been a paradigm demonstration — Predator-Prey, Game of Life, a handful of agents. This chapter is the first full domain application: land use and cover change (LUCC), the process that motivated LuccME at INPE/CCST, and the process DisSModel's own migration story (Chapter 34) keeps returning to as its central example.

It is also the first chapter of the book that assumes you know what a land change model is for. You may not, and a book read from the first page to this one has not said. So the chapter starts from the beginning — what land use and land cover are, what a model of their change is and is not, and which kinds exist — and then opens the three components of the models DisSLUCC implements one by one: what each asks of the map, the formula it computes, how it behaves, and which of the LuccME test labs uses it. Readers who already work with LUCC models can go straight to What DisSLUCC Is.

Land Use and Land Cover¶

Two expressions are used almost interchangeably, and the difference matters for modelling.

Land cover is what is physically on the ground: forest, grass, bare soil, water, buildings. It is what a satellite sees. Land use is what people do with that land: grazing cattle, growing soy, logging, living in it. The same cover can serve different uses (a grassland may be a pasture, a nature reserve or an airfield), and one use can leave different covers (a pasture may be green or burned). The maps of Chapter 28 show the difficulty: MapBiomas and PRODES are, strictly, maps of cover, from which a use is inferred — "forest" is a cover; "pasture" is a use.

The field avoids choosing by speaking of land use and cover change, LUCC, and so does this chapter. What matters to a model is not the philosophy of the word but the class: each cell of the map belongs to, or is shared among, a small set of classes, and a land change is a transition of area from one class to another. The classes of the test labs used below are three: forest (f), deforested (d) and others (outros in the Acre lab, o in the Moju lab), the last of which is not allowed to change.

Why model such changes at all? Because they move carbon between vegetation and atmosphere, alter habitats, the water cycle and local climate — and because, unlike the climate, they are the result of decisions: where roads are paved, which areas are protected, what the law allows. A model lets you ask where a change would happen if some assumptions held, and compare one set of assumptions with another.

What a Land Change Model Is, and Is Not¶

A LUCC model of the family in this chapter takes three things:

• a map of classes for the starting year (the share of each class in each cell, or the class of each cell);
• driver layers: maps believed to explain where change happens — distance to roads, to rivers or to a port, protected areas, soil fertility, slope;
• a statement of how much change each class undergoes in each year;

and returns a map of classes for each following year. Notice what is on the input side: how much change is given to the model, not computed by it. This is the most common misreading of these models. Three words help:

• A prediction claims to say what will happen.
• A projection says what happens if stated assumptions hold.
• A scenario is a coherent set of such assumptions — a demand trajectory, a policy, a road plan — that you can run and compare with another.

The models of this chapter produce projections. The amount of deforestation in 2030 comes from the demand you supply (a trend extrapolated from the past, a target, the output of an economic model); what the model adds is where that amount would go, given the drivers. If you ask such a model "how much forest will be lost by 2030?" it can only give you back the number you put in. If you ask "if 12,000 km² are lost, which cells are most likely to be, and what changes if the protected areas are enlarged?", it has something to say.

That is why LUCC models are mostly used to compare scenarios, to test whether a driver matters, and to feed other models downstream (carbon emissions, hydrology). LuccME-BR, the national model behind Chapter 31, builds national land-use scenarios for Brazil in this way (Bezerra et al., 2022). Whether a model is good — whether it reproduces change it was not fitted to — is a separate question, taken up at the end of the chapter.

Families of Land Change Models¶

Models of land change are built on a handful of ideas, and most real models combine them. A short map, with the caveat that the boundaries are blurred:

Family	Core idea	What it is good at	What it leaves out
Statistical (empirical)	Fit the probability or amount of change to the drivers, from past change	Saying which drivers matter and how much, with measurable fit	Dynamics: it describes where change was, and assumes the relation holds
Markov chain	Transition probabilities between classes, estimated from two dates	Simple, fast, gives the quantity of change	Space: it says how much moves between classes, not where
Cellular automata	Each cell changes by rules that read its neighbours (Chapter 24)	Spatial contagion and pattern	Anything that is not local — demand for the whole region, markets
Agent-based	Individual decision-makers (farmers, households) act and the map results (Chapter 25)	Mechanisms and heterogeneity of actors	Data to calibrate individual behaviour; large areas
Economic	Land goes to its most profitable use; prices and rents decide	Responses to prices and policy	Fine spatial detail
Hybrid, top-down + bottom-up	A regional part decides how much of each class is needed; a spatial part decides where	Combining the regional demand with spatial detail	Feedback from the cells back to the demand

The CLUE family (Veldkamp & Fresco, 1996; Verburg et al., 2002) and LuccME (Aguiar et al., 2012; Carneiro et al., 2013) belong to the last row. The top-down part is the demand. The bottom-up part is two components: a potential, which is a statistical model of the drivers, and an allocation, which places the demanded area on the cells in order of potential. DisSLUCC reimplements this family in Python.

Theory: Demand, Potential, Allocation¶

Every LUCC model in this chapter, continuous or discrete, runs the same three-part loop each time step — the structure LuccME took from the CLUE family of models:

1. Demand ("how much?") — the area of each land use in this step. It can come from a historical trend, a constructed scenario, an economic model, or simply a precomputed table read from a CSV.
2. Potential ("where?") — a suitability map for each class, built from driving-factor layers: distance to roads, distance to ports, protected-area status, soil fertility. This is a regression problem — predicting suitability from spatial covariates — with those covariates prepared upstream by DisSCube (Chapter 28).
3. Allocation — spatially distributing the demanded area according to the potential map, by iteration until the allocated area meets the demand within a tolerance (max_difference), or by competition between classes.

The whole idea fits in a few lines. Six cells in a row, a road at the left end, all forest; the demand asks for two cells of new pasture this year. The potential falls with the distance to the road, and the allocation gives the change to the highest potentials first. (This is a toy written with NumPy, not DisSLUCC; the sections below open the real components.)

In [2]:

# Six cells in a row, a road at the left end, all forest. The demand asks for two cells of new pasture.
dist = np.arange(6)                         # distance to the road, in cells
potential = 1.0 - 0.2 * dist                # "where?": the closer to the road, the higher the potential
demand = 2                                  # "how much?": two cells' worth of forest is converted this year

chosen = np.argsort(-potential)[:demand]    # "allocation": give the change to the highest potentials first
cleared = np.zeros(6, dtype=int); cleared[chosen] = 1

pd.DataFrame({"distance to road": dist, "potential": potential.round(2), "cleared this year": cleared}).T

# Six cells in a row, a road at the left end, all forest. The demand asks for two cells of new pasture. dist = np.arange(6) # distance to the road, in cells potential = 1.0 - 0.2 * dist # "where?": the closer to the road, the higher the potential demand = 2 # "how much?": two cells' worth of forest is converted this year chosen = np.argsort(-potential)[:demand] # "allocation": give the change to the highest potentials first cleared = np.zeros(6, dtype=int); cleared[chosen] = 1 pd.DataFrame({"distance to road": dist, "potential": potential.round(2), "cleared this year": cleared}).T

Out[2]:

	0	1	2	3	4	5
distance to road	0.0	1.0	2.0	3.0	4.0	5.0
potential	1.0	0.8	0.6	0.4	0.2	0.0
cleared this year	1.0	1.0	0.0	0.0	0.0	0.0

That is one year. The model repeats it for every step: the map after the allocation is the starting map of the next year, the demand moves on to its next row, and the potential is computed again, because the map it depends on has changed.

Is a land-change model a cellular automaton? By Chapter 24's own six-element definition — grid, neighborhood, finite states, transition rules, initial state, discrete time — yes, formally. The difference from DisSModel CA's models is one of emphasis and disciplinary origin (land-use science, not physics), and one of mechanism: allocation is driven by a demand for the whole region, not by each cell's neighbors alone.

Two ways to hold a class in a cell¶

A cell covers an area, and an area is rarely all of one class. A model has to decide what a cell is, and the two answers give the two styles of model in this chapter.

• In a continuous model, a cell holds a share of each class (a cell that is 40% cleared holds 0.4 of d and 0.6 of f). The shares of a cell add up to 1, and a change is a share that moves.
• In a discrete model, a cell holds one class. A change is a cell that flips from one class to another.

The figure builds both from the same landscape: a map of pixels is grouped into coarse cells, which hold either the share of each class or the class that covers most of the cell.

In [3]:

from matplotlib.colors import ListedColormap

rng = np.random.default_rng(7)
n, k = 60, 10                                   # a 60 x 60 "satellite" map, read as 6 x 6 coarse cells of 10 x 10 pixels
yy, xx = np.mgrid[0:n, 0:n]
dist_road = np.abs(yy - 0.35 * xx - 12) / 10    # distance to a diagonal road, in coarse-cell widths
cleared = rng.random((n, n)) < np.exp(-dist_road / 0.9)   # the chance of clearing a pixel falls with the distance

share = cleared.reshape(n // k, k, n // k, k).mean(axis=(1, 3))   # continuous cells: the share of each class
dominant = (share > 0.5).astype(int)                              # discrete cells: the class that covers most

fig, axes = plt.subplots(1, 3, figsize=(11, 3.8))
two = ListedColormap(["#2b7a3d", "#d9a441"])
axes[0].imshow(cleared, cmap=two, interpolation="nearest")
axes[0].set_title("The landscape (60 × 60 pixels)")
im = axes[1].imshow(share, cmap="YlOrBr", vmin=0, vmax=1, interpolation="nearest")
axes[1].set_title("Continuous cells: share cleared")
fig.colorbar(im, ax=axes[1], shrink=0.75)
axes[2].imshow(dominant, cmap=two, interpolation="nearest")
axes[2].set_title("Discrete cells: dominant class")
for ax in axes:
    ax.set_axis_off()
plt.tight_layout()
plt.show()
print(f"Cleared area in the landscape: {cleared.mean():.1%}; "
      f"continuous cells add up to {share.mean():.1%}; discrete cells to {dominant.mean():.1%}")

from matplotlib.colors import ListedColormap rng = np.random.default_rng(7) n, k = 60, 10 # a 60 x 60 "satellite" map, read as 6 x 6 coarse cells of 10 x 10 pixels yy, xx = np.mgrid[0:n, 0:n] dist_road = np.abs(yy - 0.35 * xx - 12) / 10 # distance to a diagonal road, in coarse-cell widths cleared = rng.random((n, n)) < np.exp(-dist_road / 0.9) # the chance of clearing a pixel falls with the distance share = cleared.reshape(n // k, k, n // k, k).mean(axis=(1, 3)) # continuous cells: the share of each class dominant = (share > 0.5).astype(int) # discrete cells: the class that covers most fig, axes = plt.subplots(1, 3, figsize=(11, 3.8)) two = ListedColormap(["#2b7a3d", "#d9a441"]) axes[0].imshow(cleared, cmap=two, interpolation="nearest") axes[0].set_title("The landscape (60 × 60 pixels)") im = axes[1].imshow(share, cmap="YlOrBr", vmin=0, vmax=1, interpolation="nearest") axes[1].set_title("Continuous cells: share cleared") fig.colorbar(im, ax=axes[1], shrink=0.75) axes[2].imshow(dominant, cmap=two, interpolation="nearest") axes[2].set_title("Discrete cells: dominant class") for ax in axes: ax.set_axis_off() plt.tight_layout() plt.show() print(f"Cleared area in the landscape: {cleared.mean():.1%}; " f"continuous cells add up to {share.mean():.1%}; discrete cells to {dominant.mean():.1%}")

Cleared area in the landscape: 28.6%; continuous cells add up to 28.6%; discrete cells to 25.0%

The continuous cells keep information that the discrete ones discard. The cell in the third row, second column is 47% cleared and shows as forest in the discrete map; another cell of the same row, at 51%, shows as pasture — two nearly identical cells with opposite classes. The areas show the same loss: the landscape is 28.6% cleared, the continuous cells add up to exactly that, and the discrete map holds 25.0%. The discrete map, in return, is easier to compare with a classified image, because both have one class per cell. The choice between the two is taken up again in Choosing Between Continuous and Discrete.

A few words are used throughout the chapter:

Word	Meaning
Land-use type (class)	One entry of the map's legend: f, d, outros in the Acre lab
Step	One year (or other fixed period) of the simulation
Driver	A map layer believed to explain change, with one value per cell
Coefficient (betas, const)	The weight of a driver in a regression, and the constant term
Demand	The area each class must hold at a step
Potential	A number per cell and class: how strongly the cell attracts the class
Allocation	The step that turns the demand and the potentials into a new map
Elasticity	A factor that scales a potential into a change, adjusted until the area matches the demand
Region	A subset of cells with its own coefficients; the labs here use one

What DisSLUCC Is¶

DisSLUCC implements LuccME's components — demand, potential, allocation — in Python, on the raster substrate of DisSModel. It is one package, disslucc, in which the continuous and discrete components live side by side, as they do in LuccME itself. (Earlier drafts of this book used two packages, disslucc-continuous and disslucc-discrete; they were merged into this one.) It is meant to be run next to TerraME and LuccME, which remain the reference, not to replace them.

Role	Components
Demand — how much?	DemandInline, DemandPreComputedValues (a table), DemandComputeTwoDates, DemandComputeThreeDates (from land-use layers)
Potential — where?	PotentialLinearRegression (continuous), PotentialSpatialLagRegression (continuous, with a spatial-lag term), PotentialDLogisticRegression (discrete)
Allocation — how is the change placed?	AllocationClueLike (continuous, CLUE-like), AllocationClueLikeSaturation (continuous, with saturation), AllocationDClueSLike (discrete, CLUE-S-like)

Each component is an ordinary DisSModel Model, registered in the same Environment in the order it is built — the lifecycle of Chapter 22, nothing more. There are two entry points to the same components: a script (build the models in Python and call env.run() — the fastest way to experiment) and an executor (disslucc.executors, which wraps the same components in a ModelExecutor with an ExperimentRecord, reads its configuration from a TOML file, and has a command line). Both give the same result; the executor adds the provenance record of Chapter 29.

The Cases: Two Study Areas of the LuccME Test Labs¶

The components below are illustrated with the data of the test labs that ship with LuccME 3.1. There are 21 labs, and they use only two sets of cells, so there are really two cases.

• Acre (csAC, labs 01 to 09): 6,574 square cells of 5 km, in a projected system. Three classes: forest (f), deforested (d) and others (outros), the shares of each in every cell for 2008, and again for 2011 and 2014 (f2011, d2011, f2014, d2014), plus the driver layers. The labs use the continuous components: a cell holds shares. A cell is 25 km².
• Moju (cs_moju, labs 10 to 21): 5,914 square cells of 0.0045° (about 500 m), in western Pará. Three classes, f, d and o, each cell holding one of them, for 1999, and for 2004 and 2007 (f04, d04, f07, d07), plus the drivers. The labs use the discrete components. The lab sets the cell area to 1, so areas are counted in cells.

The data are the ones the DisSLUCC repository vendors; the cell below downloads them from the release the book was checked against (v0.5.0), so the chapter does not copy them.

In [4]:

import tempfile, urllib.request
from pathlib import Path
import geopandas as gpd

# The labs' data are the ones the DisSLUCC repository vendors, pinned to the release the book was checked against.
BASE = "https://raw.githubusercontent.com/DisSModel/disslucc/v0.5.0/data/input/"
CACHE = Path(tempfile.gettempdir()) / "disslucc_v0.5.0_data"; CACHE.mkdir(exist_ok=True)

def fetch(name):
    path = CACHE / name
    if not path.exists():
        urllib.request.urlretrieve(BASE + name, path)
    return path

acre = gpd.read_file(fetch("csAC.zip"))        # lab01 to lab09
moju = gpd.read_file(fetch("cs_moju.zip"))     # lab10 to lab21

for name, gdf, cls in [("Acre (csAC)", acre, ["f", "d", "outros"]), ("Moju (cs_moju)", moju, ["f", "d", "o"])]:
    first = gdf.geometry.iloc[0].bounds
    print(f"{name}: {len(gdf):,} cells of {first[2]-first[0]:g} x {first[3]-first[1]:g} {'metres' if gdf.crs.is_projected else 'degrees'}")
    others = [c for c in gdf.columns if c not in cls + ["geometry"]]
    print("   classes:", cls)
    print("   other columns (drivers, later-year classes, indices):", others)

import tempfile, urllib.request from pathlib import Path import geopandas as gpd # The labs' data are the ones the DisSLUCC repository vendors, pinned to the release the book was checked against. BASE = "https://raw.githubusercontent.com/DisSModel/disslucc/v0.5.0/data/input/" CACHE = Path(tempfile.gettempdir()) / "disslucc_v0.5.0_data"; CACHE.mkdir(exist_ok=True) def fetch(name): path = CACHE / name if not path.exists(): urllib.request.urlretrieve(BASE + name, path) return path acre = gpd.read_file(fetch("csAC.zip")) # lab01 to lab09 moju = gpd.read_file(fetch("cs_moju.zip")) # lab10 to lab21 for name, gdf, cls in [("Acre (csAC)", acre, ["f", "d", "outros"]), ("Moju (cs_moju)", moju, ["f", "d", "o"])]: first = gdf.geometry.iloc[0].bounds print(f"{name}: {len(gdf):,} cells of {first[2]-first[0]:g} x {first[3]-first[1]:g} {'metres' if gdf.crs.is_projected else 'degrees'}") others = [c for c in gdf.columns if c not in cls + ["geometry"]] print(" classes:", cls) print(" other columns (drivers, later-year classes, indices):", others)

Acre (csAC): 6,574 cells of 5000 x 5000 metres
   classes: ['f', 'd', 'outros']
   other columns (drivers, later-year classes, indices): ['object_id0', 'row', 'col', 'f2011', 'd2011', 'f2014', 'd2014', 'rios_todos', 'uc_pi', 'fertilidad', 'rodovias', 'assentamen', 'ti', 'uc_us', 'ap', 'dist_riobr', 'reg_prot', 'd_difReal']
Moju (cs_moju): 5,914 cells of 0.0045 x 0.0045 degrees
   classes: ['f', 'd', 'o']
   other columns (drivers, later-year classes, indices): ['object_id_', 'col', 'lin', 'dist_br', 'dist_area_', 'dist_curua', 'dist_porto', 'dist_mojui', 'dist_rios', 'dist_estra', 'dist_belte', 'dist_santa', 'media_decl', 'd04', 'd07', 'inverse_di', 'dist_rios_', 'f04', 'f07']

In [5]:

fig, axes = plt.subplots(1, 3, figsize=(12, 3.6))
acre.plot(column="d", cmap="YlOrBr", vmin=0, vmax=1, ax=axes[0], linewidth=0)
axes[0].set_title("Acre, 2008: share deforested (d)")
acre.plot(column="d2014", cmap="YlOrBr", vmin=0, vmax=1, ax=axes[1], linewidth=0, legend=True, legend_kwds={"shrink": 0.7})
axes[1].set_title("Acre, 2014: share deforested (d2014)")
from matplotlib.colors import ListedColormap
from matplotlib.patches import Patch
moju.assign(cls=moju[["f", "d", "o"]].values.argmax(axis=1)).plot(
    column="cls", cmap=ListedColormap(["#2b7a3d", "#d9a441", "#888888"]), ax=axes[2], linewidth=0)
axes[2].set_title("Moju, 1999: class of each cell")
axes[2].legend(handles=[Patch(color=c, label=l) for c, l in zip(["#2b7a3d", "#d9a441", "#888888"], ["f", "d", "o"])], loc="lower center", bbox_to_anchor=(0.5, -0.1), ncol=3, frameon=False, fontsize=8)
for ax in axes:
    ax.set_axis_off()
plt.tight_layout()
plt.show()

fig, axes = plt.subplots(1, 3, figsize=(12, 3.6)) acre.plot(column="d", cmap="YlOrBr", vmin=0, vmax=1, ax=axes[0], linewidth=0) axes[0].set_title("Acre, 2008: share deforested (d)") acre.plot(column="d2014", cmap="YlOrBr", vmin=0, vmax=1, ax=axes[1], linewidth=0, legend=True, legend_kwds={"shrink": 0.7}) axes[1].set_title("Acre, 2014: share deforested (d2014)") from matplotlib.colors import ListedColormap from matplotlib.patches import Patch moju.assign(cls=moju[["f", "d", "o"]].values.argmax(axis=1)).plot( column="cls", cmap=ListedColormap(["#2b7a3d", "#d9a441", "#888888"]), ax=axes[2], linewidth=0) axes[2].set_title("Moju, 1999: class of each cell") axes[2].legend(handles=[Patch(color=c, label=l) for c, l in zip(["#2b7a3d", "#d9a441", "#888888"], ["f", "d", "o"])], loc="lower center", bbox_to_anchor=(0.5, -0.1), ncol=3, frameon=False, fontsize=8) for ax in axes: ax.set_axis_off() plt.tight_layout() plt.show()

Demand: How Much?¶

Demand is the one component that does not look at geography. It says, for each class, how much area the region must hold at each step. In DisSLUCC a demand component keeps a table with one row per step and one column per class, and answers three questions for the other components: the demand of a class now, its demand in the previous step, and its direction — increasing, decreasing or static. The direction tells the allocation whether a class may gain cells in this step or lose them.

Component	Where the numbers come from	Labs that use it
DemandPreComputedValues (and DemandInline, the same with the table given in Python instead of a CSV)	A table computed elsewhere: a trend, a target, the output of an economic model, or by hand	01, 02, 03, 06 to 15, 18, 21
DemandComputeTwoDates	The map itself: the areas of the classes at the start and at a final year	04, 16, 19, 20
DemandComputeThreeDates	The map itself, at the start, a middle year and a final year	05, 17

The two last components build the table from land-use layers. For each class $c$ they take its total area $A_c$ at the start year $t_0$ and at the final year $t_1$, and interpolate linearly:

$$\begin{gathered} s_c = \frac{|A_c(t_1) - A_c(t_0)|}{t_1 - t_0}, \\ D_c(t) = \operatorname{round}\big(D_c(t-1) \pm s_c\big) \end{gathered}$$

where the sign is the direction (the sign of $A_c(t_1) - A_c(t_0)$, unless you give it), $D_c(t_0) = A_c(t_0)$, and the demand at the reference year is set to its exact value. With three dates the same is done twice, from the start to the middle year and from the middle to the final year, each segment with its own step. Two details follow LuccME, and matter if you compare with it: every year adds the step to the already rounded demand of the previous year and rounds again, half up; so the demand is a whole number and can differ by a few units from a table computed with exact arithmetic. After the final year, the same step keeps being applied.

The cell computes the demand of d in Acre by both. Lab04 uses two dates (2008 and 2014), lab05 three (2008, 2011 and 2014).

In [6]:

import pandas as pd
from disslucc.components.demand.dates import interpolate_two_dates, interpolate_three_dates

CELL_AREA = 25.0                                    # km² per Acre cell
start = [float(acre[c].sum()) * CELL_AREA for c in ("f", "d", "outros")]       # areas in 2008
mid   = [float(acre[c].sum()) * CELL_AREA for c in ("f2011", "d2011", "outros")]
final = [float(acre[c].sum()) * CELL_AREA for c in ("f2014", "d2014", "outros")]

two   = interpolate_two_dates(start, final, final_offset=6, n_years=7)                       # lab04
three = interpolate_three_dates(start, mid, final, middle_offset=3, final_offset=6, n_years=7)  # lab05

years = list(range(2008, 2015))
table = pd.DataFrame({"d, two dates": [r[1] for r in two], "d, three dates": [r[1] for r in three]}, index=years)
print(table.T.to_string())

# The table that lab01 ships is the three-date demand, rounded the same way
lab1 = pd.read_csv(fetch("examples_demand_lab1.csv"))
print("largest difference between lab01's table and the three-date demand (km²):",
      float(np.abs(lab1.to_numpy() - np.array(three)).max().round(2)))

import pandas as pd from disslucc.components.demand.dates import interpolate_two_dates, interpolate_three_dates CELL_AREA = 25.0 # km² per Acre cell start = [float(acre[c].sum()) * CELL_AREA for c in ("f", "d", "outros")] # areas in 2008 mid = [float(acre[c].sum()) * CELL_AREA for c in ("f2011", "d2011", "outros")] final = [float(acre[c].sum()) * CELL_AREA for c in ("f2014", "d2014", "outros")] two = interpolate_two_dates(start, final, final_offset=6, n_years=7) # lab04 three = interpolate_three_dates(start, mid, final, middle_offset=3, final_offset=6, n_years=7) # lab05 years = list(range(2008, 2015)) table = pd.DataFrame({"d, two dates": [r[1] for r in two], "d, three dates": [r[1] for r in three]}, index=years) print(table.T.to_string()) # The table that lab01 ships is the three-date demand, rounded the same way lab1 = pd.read_csv(fetch("examples_demand_lab1.csv")) print("largest difference between lab01's table and the three-date demand (km²):", float(np.abs(lab1.to_numpy() - np.array(three)).max().round(2)))

                 2008   2009   2010   2011   2012   2013   2014
d, two dates    19983  20254  20525  20796  21067  21338  21607
d, three dates  19983  20239  20495  20750  21036  21322  21607
largest difference between lab01's table and the three-date demand (km²): 1.05

Both reach the same 2008 and 2014 values; they differ in between, by up to 46 km² in 2011, the middle year, because the three-date version bends the line through the area the map had in 2011. Lab01 does not compute its demand: it ships the table as a CSV — and the table is the three-date demand, which the last line confirms within about 1 km² (the rounding of the CSV).

Two points are easy to miss. First, the demand is where the scenario lives: change the table and you get a different future with the same potentials. The interpolation here only extends the observed trend; DisSLUCC has no economic model that would produce a demand from prices or policies, and no component to forecast it. Second, a demand is a total for the region. Nothing in it says where the change goes — that is the potential's job.

Where a Demand Comes From¶

DisSLUCC takes the demand as given, but a reader building a model has to produce it, and the literature on top-down models describes two ways, which mix in practice.

• Trend extrapolation. The rates of change observed in the past are projected forward. It answers "what happens if things continue as they were?", and it is what the interpolation of two or three dates does. It needs no assumption about the future, which is its strength and its limit: it cannot represent a change of policy or of market.
• Scenarios. The demand is built from a set of premises: a rate of change defined for each year by experts, or the output of an economic model that turns population growth, income or prices into areas. It lets one ask "what if?" — a stricter enforcement of the law, a fall in the price of beef — and compare the futures that result.

In the regional model of the Brazilian Amazon on which LuccME's approach builds (Aguiar, 2006), the demand for deforestation is an external variable that stands for the pressure of national and international markets, regionalised by macro-region; roads and protected areas enter only as spatial determinants, in the potential, and do not change the total. Keeping the two apart is what lets the same model separate the effect of how much is cleared from the effect of where: the demand moves with the market scenario, the potential with the paving of a road or the creation of a reserve. The same distinction applies to the potential's drivers: a distance to roads helps to predict where clearing happens, but it is a proximate cause and not the underlying one (Veldkamp & Lambin, 2001), so a model that fits well says little about why the clearing happens.

Potential: Where?¶

The potential is a number for each cell and each class, and it answers: how strongly does this cell attract this class, compared with what it holds now? A positive potential says the cell has room to gain the class; a negative one says it tends to lose it. The allocation reads these maps and nothing else about the drivers.

The drivers enter here. Each potential component turns the driver layers — distances to roads and rivers, protected areas, soil fertility — into one map per class, through coefficients (const and betas) that say how strongly each driver counts. In the labs the coefficients were estimated beforehand, and DisSLUCC takes them as given; it has no routine to fit them. What the three components below differ in is the form of the relation between drivers and potential:

Component	The relation	Style	Labs that use it
PotentialLinearRegression	A linear regression of the share of the class	continuous	01
PotentialSpatialLagRegression	The same, plus a term for the neighbouring cells	continuous	02 to 07
PotentialDLogisticRegression	A logistic regression of the probability of the class, plus a bonus for the cells that already have it	discrete	13, 15, 16, 17, 21

(LuccME's names carry a letter for the style, PotentialCLinearRegression, PotentialDLogisticRegression; DisSLUCC drops it from the continuous ones.)

Linear regression¶

For class $c$ in cell $i$, with drivers $x_{k,i}$ and the share $s_{c,i}$ the cell holds now:

$$\begin{gathered} r_{c,i} = \operatorname{clip}_{[0,1]}\Big(\alpha_c + \sum_k \beta_{c,k}\, x_{k,i}\Big), \\ p_{c,i} = r_{c,i} - s_{c,i} \end{gathered}$$

Here $\alpha_c$ is the constant (const), $\beta_{c,k}$ the coefficient of driver $k$ (betas), $r$ is the share the regression thinks the cell should hold, clipped to a valid share, and the potential $p$ is the difference between that and what the cell holds. A class marked is_log is regressed in log space, $r = 10^{\alpha + \sum \beta x}$, before the clipping. If one class is declared as not transitionable (land_use_no_data, outros in the labs), $r$ is also multiplied by one minus the share the cell already has in that class, so a cell mostly taken by it is asked for little of the others.

At each step after the first, the constant follows the demand, by a small amount: $\alpha_c \leftarrow \alpha_c + 0.01\,\dfrac{D_c(t) - D_c(t-1)}{D_c(t-1)}$. A class whose demand grows by 5% in a year has its constant raised by 0.0005 — a nudge, not a re-fit.

The figure draws one regression, the urban one of the synthetic scenario further below (constant 0.3, coefficient −0.6 on the distance to the road).

In [7]:

dist = np.linspace(0, 1, 101)                       # distance to the road, scaled to 0-1
const, beta = 0.3, -0.6                             # the "urban" regression of the synthetic scenario below, one driver only

reg = np.clip(const + beta * dist, 0, 1)            # what the cell "should" hold, according to the regression
fig, axes = plt.subplots(1, 2, figsize=(10, 3.4))
axes[0].plot(dist, const + beta * dist, "--", color="grey", label="regression, before clipping")
axes[0].plot(dist, reg, color="C1", label="regression, clipped to [0, 1]")
axes[0].set_xlabel("distance to the road"); axes[0].set_ylabel("share the cell should hold"); axes[0].legend()
for past in (0.0, 0.2, 0.5):
    axes[1].plot(dist, reg - past, label=f"cell holds {past:.0%} already")
axes[1].axhline(0, color="k", lw=0.8)
axes[1].set_xlabel("distance to the road"); axes[1].set_ylabel("potential = regression − share now"); axes[1].legend()
plt.tight_layout()
plt.show()

# The constant follows the demand: when the demanded area of a class grows by 5% in a year...
prev, curr = 100.0, 105.0
print("the constant of that class rises by", 0.01 * (curr - prev) / prev)

dist = np.linspace(0, 1, 101) # distance to the road, scaled to 0-1 const, beta = 0.3, -0.6 # the "urban" regression of the synthetic scenario below, one driver only reg = np.clip(const + beta * dist, 0, 1) # what the cell "should" hold, according to the regression fig, axes = plt.subplots(1, 2, figsize=(10, 3.4)) axes[0].plot(dist, const + beta * dist, "--", color="grey", label="regression, before clipping") axes[0].plot(dist, reg, color="C1", label="regression, clipped to [0, 1]") axes[0].set_xlabel("distance to the road"); axes[0].set_ylabel("share the cell should hold"); axes[0].legend() for past in (0.0, 0.2, 0.5): axes[1].plot(dist, reg - past, label=f"cell holds {past:.0%} already") axes[1].axhline(0, color="k", lw=0.8) axes[1].set_xlabel("distance to the road"); axes[1].set_ylabel("potential = regression − share now"); axes[1].legend() plt.tight_layout() plt.show() # The constant follows the demand: when the demanded area of a class grows by 5% in a year... prev, curr = 100.0, 105.0 print("the constant of that class rises by", 0.01 * (curr - prev) / prev)

the constant of that class rises by 0.0005

The left panel is the regression: it asks for 0.3 of the class at the road and nothing from half the scaled distance on. The right panel shows what that means for the allocation. The potential is relative to what the cell holds: the same driver values give a lower potential to a cell that already holds more of the class, and beyond half the distance, where the regression asks for nothing, a cell holding 20% has potential −0.2 — the drivers say that share is not justified, and that cell is where the class is lost if its demand falls.

Lab01 uses this potential for the three classes (the code is in the whole lab section further on). Its coefficients tell a story. For deforestation (d) the coefficient of the roads layer (rodovias, a distance in metres) is negative, so deforestation is lower far from roads; so is the one of the sustainable-use conservation units (uc_us, a share of the cell), while the one of the settlements (assentamen) is positive. The outros class has a constant of 0 and no coefficients, because it does not change.

Linear regression with a spatial lag¶

Land change is rarely scattered at random: it starts where change already happened. The spatial-lag regression adds that to the regression with one more term, the mean of the class in the neighbourhood of the cell:

$$\begin{gathered} Y_{c,i} = \frac{1}{n_i + 1}\Big(\tilde s_{c,i} + \sum_{j \in N(i)} \tilde s_{c,j}\Big), \\ r_{c,i} = \alpha_c + \sum_k \beta_{c,k}\, x_{k,i} + \rho_c\, Y_{c,i}, \\ p_{c,i} = r_{c,i} - s_{c,i} \end{gathered}$$

$N(i)$ are the up to eight cells around $i$ (the Moore neighbourhood), $n_i$ how many of them count, and $\tilde s$ the share rescaled by the part of the cell that is not no-data. $\rho_c$ (ro) is the weight of the neighbourhood. If the value so computed is above max_reg, the cell is asked for 1; below min_reg, for 0. There is no clipping to $[0,1]$ before that, and the result is scaled by the available part of the cell, as in the linear case. Two things follow LuccME on purpose: the clip tests the value built with the original constant but writes 1 or 0 rather than the bounds, and the adaptation of the constant to the demand, the same formula as before, accumulates from year to year.

LuccME-BR uses this potential for every class (Bezerra et al., 2022). In the labs, labs 02 to 07 use it.

The figure shows the effect with no driver at all: only a constant of 0.05 and the neighbours, for a clearing in the middle of a forest.

In [8]:

from disslucc.components.potential.spatial_lag import spatial_lag_regression
from disslucc.schemas import SpatialLagRegressionSpec

past = np.zeros((15, 15)); past[6:9, 6:9] = 1.0                 # a clearing in the middle of a forest
valid = np.ones_like(past, dtype=bool)

fig, axes = plt.subplots(1, 3, figsize=(12, 3.6))
axes[0].imshow(past, cmap="YlOrBr", vmin=0, vmax=1); axes[0].set_title("Share cleared now")
for ax, ro in zip(axes[1:], (0.0, 0.8)):
    spec = SpatialLagRegressionSpec(const=0.05, ro=ro)            # no driver at all: only the constant and the neighbours
    reg, pot = spatial_lag_regression(spec, past, {}, valid)      # reg: the share the cell should hold; pot = reg - past
    im = ax.imshow(reg, cmap="YlOrBr", vmin=0, vmax=1); ax.set_title(f"Share asked of each cell\nro = {ro}")
fig.colorbar(im, ax=axes[1:], shrink=0.75)
for ax in axes: ax.set_axis_off()
plt.show()
print(f"with ro = {ro}: cell touching the clearing along an edge {reg[5, 7]:.2f}, at a corner {reg[5, 5]:.2f}, in its centre {reg[7, 7]:.2f}")

from disslucc.components.potential.spatial_lag import spatial_lag_regression from disslucc.schemas import SpatialLagRegressionSpec past = np.zeros((15, 15)); past[6:9, 6:9] = 1.0 # a clearing in the middle of a forest valid = np.ones_like(past, dtype=bool) fig, axes = plt.subplots(1, 3, figsize=(12, 3.6)) axes[0].imshow(past, cmap="YlOrBr", vmin=0, vmax=1); axes[0].set_title("Share cleared now") for ax, ro in zip(axes[1:], (0.0, 0.8)): spec = SpatialLagRegressionSpec(const=0.05, ro=ro) # no driver at all: only the constant and the neighbours reg, pot = spatial_lag_regression(spec, past, {}, valid) # reg: the share the cell should hold; pot = reg - past im = ax.imshow(reg, cmap="YlOrBr", vmin=0, vmax=1); ax.set_title(f"Share asked of each cell\nro = {ro}") fig.colorbar(im, ax=axes[1:], shrink=0.75) for ax in axes: ax.set_axis_off() plt.show() print(f"with ro = {ro}: cell touching the clearing along an edge {reg[5, 7]:.2f}, at a corner {reg[5, 5]:.2f}, in its centre {reg[7, 7]:.2f}")

with ro = 0.8: cell touching the clearing along an edge 0.32, at a corner 0.14, in its centre 0.85

With $\rho = 0$ every cell is asked for the constant, 0.05, and the clearing is invisible to its neighbours. With $\rho = 0.8$ the cells that touch the clearing are asked for 0.32 along an edge and 0.14 at a corner (the line printed under the figure) — about six and nearly three times the constant — in a halo one cell thick, because the neighbourhood is the 3 × 3 window. The centre of the clearing is asked for 0.85, a little less than the 1.0 it holds: its potential is slightly negative, and only the halo has room to gain. The allocation therefore gives the next change to the halo first, and the clearing grows from its edges. A $\rho$ near 0 makes change scatter where the drivers say; a large $\rho$ makes it spread like a stain.

Logistic regression¶

In the discrete style a cell holds one class, so the question is not "what share?" but "which class, and with what probability?". The logistic regression models exactly that: the probability of the class as an S-shaped function of the drivers, between 0 and 1.

$$\begin{gathered} z_{c,i} = \alpha_c + \sum_k \beta_{c,k}\, x_{k,i}, \\ \pi_{c,i} = \frac{1}{1 + e^{-z_{c,i}}}, \\ p_{c,i} = \pi_{c,i} + \varepsilon_c\,\mathbf{1}[\text{cell } i \text{ is class } c \text{ now}] \end{gathered}$$

The potential is the probability plus the elasticity $\varepsilon_c$ (elasticity in LogisticRegressionSpec) for the cells that already are in the class: how hard a class is to remove from a cell. A high elasticity locks a cell into its class (a forest that is protected, a built-up area that is not abandoned); with 0, only the probability counts. The raw probability is also kept, in a separate layer.

Unlike the continuous potential, this one does not subtract what the cell holds: a cell either is the class or is not, and the allocation compares the potentials of all the classes for the same cell. Labs 13, 15, 16, 17 and 21 use it, with the distance layers of Moju (dist_br, dist_estra, dist_rios, among others) as drivers.

In [9]:

from scipy.special import expit

x = np.linspace(-1, 1, 201)
fig, axes = plt.subplots(1, 2, figsize=(10, 3.4))
for const, beta in [(0.0, 3.0), (-2.0, 3.0), (0.0, 8.0)]:
    axes[0].plot(x, expit(const + beta * x), label=f"const = {const}, beta = {beta}")
axes[0].set_xlabel("driver"); axes[0].set_ylabel("probability of the class"); axes[0].legend(fontsize=8)

prob = expit(-2.0 + 3.0 * x)
axes[1].plot(x, prob, label="probability (elasticity 0)")
axes[1].plot(x, prob + 0.5, "--", label="potential if the cell already has the class (elasticity 0.5)")
axes[1].set_xlabel("driver"); axes[1].set_ylabel("potential"); axes[1].legend(fontsize=8)
plt.tight_layout()
plt.show()

from scipy.special import expit x = np.linspace(-1, 1, 201) fig, axes = plt.subplots(1, 2, figsize=(10, 3.4)) for const, beta in [(0.0, 3.0), (-2.0, 3.0), (0.0, 8.0)]: axes[0].plot(x, expit(const + beta * x), label=f"const = {const}, beta = {beta}") axes[0].set_xlabel("driver"); axes[0].set_ylabel("probability of the class"); axes[0].legend(fontsize=8) prob = expit(-2.0 + 3.0 * x) axes[1].plot(x, prob, label="probability (elasticity 0)") axes[1].plot(x, prob + 0.5, "--", label="potential if the cell already has the class (elasticity 0.5)") axes[1].set_xlabel("driver"); axes[1].set_ylabel("potential"); axes[1].legend(fontsize=8) plt.tight_layout() plt.show()

The left panel shows how the constant moves the curve to the right or left and how the coefficient makes it steeper — a large coefficient makes the class switch on over a narrow range of the driver. The right panel shows the elasticity: for a cell already in the class the potential sits 0.5 above the probability, so the class keeps the cell unless another class beats it by more than that.

LuccME has other potentials that DisSLUCC does not have yet: mixtures of the linear and the lag regression, sample-based potentials, and rules based on the neighbours or the inverse distance, among them. They appear in the table of the 21 labs at the end of the chapter, with the labs they would unlock.

Allocation: Placing the Change¶

The allocation takes the demand and the potential maps and produces the new map. It is the part that decides which cells actually change, and it must do so in a way that the total area of each class meets the demand within the tolerance max_difference. The three components differ in how.

Component	Style	Idea	Labs that use it (in LuccME)
AllocationClueLike	continuous	Scale each potential by an elasticity, adjusted per class until the area meets the demand	01, 02, 04, 05, 07, 08, 09
AllocationClueLikeSaturation	continuous	The same, with the change slowed in cells whose neighbourhood is already mostly converted	03, 06
AllocationDClueSLike	discrete	The classes compete for each cell; a correction per class is adjusted until the areas meet the demand	14, 15, 16, 17, 18, 20

CLUE-like allocation (continuous)¶

Each year after the first (in the starting year nothing is allocated) the allocation runs a loop that has one number per class, the elasticity $\varepsilon_c$, and repeats until the areas fit:

1. Start. $\varepsilon_c$ takes its initial value (0.1 in the labs) for every class.
2. Change. The change of class $c$ in cell $i$ is $\Delta_{c,i} = \varepsilon_c\, p_{c,i}$: a share to add if the potential is positive, to remove if negative. Changes smaller than min_change are ignored and larger ones capped at max_change.
3. Direction. The change is applied only in the direction the demand asks. If the demand of a class increases, only cells with a positive potential gain; if it decreases, only cells with a negative one lose. (This is the static setting of a class: $-1$ follows the demand, $0$ is free, $1$ is fixed and never changes.) The new share is kept within $[0,1]$ and within the class's min_value and max_value.
4. Cell correction. The shares of a cell must add up to 1. Where they do not, a correction spreads the excess in proportion to the changes made (up to 25 rounds). LuccME's own version of this step never runs, because of a typo; DisSLUCC runs it by default and has a switch to reproduce LuccME (see Checked Against LuccME).
5. Compare. The allocated area $a_c$ of each class is the sum of its shares times the cell area. If the largest gap $\max_c |a_c - D_c|$ is within max_difference, the allocation stops.
6. Adapt. Otherwise each elasticity is rescaled: $\varepsilon_c \leftarrow \varepsilon_c \cdot D_c / a_c$ if the demand increases, and $\varepsilon_c \cdot a_c / D_c$ if it decreases. Elasticity is kept between min_elasticity and max_elasticity (0.001 and 1.5); when a limit is hit, the constant of the potential is nudged by 0.1 so that the map itself moves, and the loop repeats. When half of max_iteration has passed the tolerance is doubled, and at max_iteration an error is raised.
7. Complement. At the end the class named complementar_lu takes whatever the others leave: $s^{\mathrm{compl}} = 1 - \sum s^{\mathrm{others}}$. In the labs it is forest.

It is a feedback controller with one knob per class: the potential decides which cells have the highest claim, and the elasticity decides how much of it is honoured, until the area is right. The toy below keeps only that part — one class, 20 cells, no limits on the change, no correction — to make the feedback visible.

In [10]:

# The elasticity loop of CLUE, stripped to its core: 20 cells, one class, no bounds, no correction.
potential = np.array([4, 3, 2, 1, .5, .3, .2, .1, .05, 0, -.1, -.2, -.3, -.4, -.5, -.6, -.7, -.8, -.9, -1.0])   # potential of the class, cell by cell
demand = 5.0                                    # area demanded, in cells
tolerance = 0.05                                # max_difference
share = np.zeros(20)                            # the class holds nothing yet
elasticity = 0.1                                # initial elasticity

rows = []
for it in range(40):
    new = np.where(potential > 0, np.clip(share + potential * elasticity, 0, 1), share)   # change = potential x elasticity, only where potential > 0, share kept in [0, 1]
    area = new.sum()
    rows.append((it, round(float(elasticity), 3), round(float(area), 2)))
    if abs(area - demand) <= tolerance:
        break
    elasticity *= demand / area                  # allocated too little -> more elastic; too much -> less
print(pd.DataFrame(rows, columns=["iteration", "elasticity", "allocated area"]).to_string(index=False))

fig, ax = plt.subplots(figsize=(6.5, 3))
ax.bar(range(20), potential, color="lightgrey", label="potential")
ax.bar(range(20), new, width=0.5, color="C1", label="share allocated")
ax.set_xlabel("cell"); ax.legend()
plt.tight_layout()
plt.show()

# The elasticity loop of CLUE, stripped to its core: 20 cells, one class, no bounds, no correction. potential = np.array([4, 3, 2, 1, .5, .3, .2, .1, .05, 0, -.1, -.2, -.3, -.4, -.5, -.6, -.7, -.8, -.9, -1.0]) # potential of the class, cell by cell demand = 5.0 # area demanded, in cells tolerance = 0.05 # max_difference share = np.zeros(20) # the class holds nothing yet elasticity = 0.1 # initial elasticity rows = [] for it in range(40): new = np.where(potential > 0, np.clip(share + potential * elasticity, 0, 1), share) # change = potential x elasticity, only where potential > 0, share kept in [0, 1] area = new.sum() rows.append((it, round(float(elasticity), 3), round(float(area), 2))) if abs(area - demand) <= tolerance: break elasticity *= demand / area # allocated too little -> more elastic; too much -> less print(pd.DataFrame(rows, columns=["iteration", "elasticity", "allocated area"]).to_string(index=False)) fig, ax = plt.subplots(figsize=(6.5, 3)) ax.bar(range(20), potential, color="lightgrey", label="potential") ax.bar(range(20), new, width=0.5, color="C1", label="share allocated") ax.set_xlabel("cell"); ax.legend() plt.tight_layout() plt.show()

 iteration  elasticity  allocated area
         0       0.100            1.12
         1       0.448            3.86
         2       0.581            4.25
         3       0.683            4.47
         4       0.765            4.64
         5       0.823            4.77
         6       0.863            4.86
         7       0.889            4.91
         8       0.905            4.95
         9       0.915            4.97

Started at 0.1 the class gets 1.12 cells' worth, far from the 5 demanded; the elasticity is multiplied by 5/1.12 and the allocation rises to 3.86. The gain per iteration shrinks as the best cells reach a share of 1 and cannot grow any more, so the last steps are slow. The real component stops as soon as the gap is within max_difference, and that choice matters more than it looks, as the lab below shows.

Saturation¶

AllocationClueLikeSaturation is the same loop with one more brake. For every cell, at the start of each year, it computes a saturation indicator: the share of the available area (not no-data, not protected) in the 3 × 3 window around the cell that is no longer in the complementary class — how converted the neighbourhood already is. Where the indicator exceeds change_limiar_value, a change in the direction of the demand is halved or capped at max_change_above_limiar. The idea is that a neighbourhood that is already largely converted has less left to convert, so the change in its cells is slowed. This variant also runs the cell correction exactly as LuccME writes it, and uses a few different bounds. LuccME-BR uses it with the spatial-lag potential; labs 03 and 06 do.

CLUE-S-like allocation (discrete)¶

In the discrete style a cell holds one class, so the allocation is a competition. For each class the algorithm keeps a correction $\lambda_c$, starting at 0, and each year repeats:

1. Every cell takes the class with the best score among those that its transition matrix allows from the class it had at the start of the year:

$$ \hat c_i = \arg\max_{c \in A(i)} \Big[(1 + \tau_{c,i})\, p_{c,i} + \lambda_c\Big] $$

2. The area $a_c$ of each class is counted and compared with the demand, $\delta_c = D_c - a_c$.
3. If $\max_c |\delta_c|$ is within max_difference, it stops. Otherwise $\lambda_c \leftarrow \lambda_c + \phi\,\delta_c$, with $\phi$ the factor_iteration, and the competition is run again.

$A(i)$ is the set of classes allowed for cell $i$: the transition matrix is indexed [region][from][to], with 1 where a change is allowed, and [1, 1, 0] for forest in the Moju example means that forest may stay forest or become deforested but not become the third class. $\tau$ is an optional attraction or repulsion layer for a class (zero by default). A class that fell short gets a positive $\lambda$ and wins more cells; one that overshot gets a negative one. There is no elasticity to adjust and the potential is never changed: the loop moves only the corrections.

The toy below takes 12 cells, two classes, and a transition matrix that lets forest become pasture but not the reverse.

In [11]:

# The competition of CLUE-S, stripped to its core: 12 cells, 2 classes (forest, pasture), pasture may not go back to forest.
rng = np.random.default_rng(3)
n = 12
pot = np.column_stack([rng.uniform(0.3, 0.7, n),                      # potential of forest in each cell
                       np.linspace(0.95, 0.15, n) + rng.normal(0, 0.05, n)])  # potential of pasture: falls along a gradient
now = np.array([0] * 10 + [1] * 2)                                    # 10 forest cells, 2 pasture cells
allowed = np.array([[True, True], [False, True]])[now]                # row = class now; forest -> {forest, pasture}; pasture -> {pasture}
demand = np.array([6, 6])                                             # the demand: 6 cells of each class
tolerance, factor, iter_vec = 0, 0.005, np.zeros(2)   # tolerance = max_difference, factor = factor_iteration

rows = []
for it in range(200):
    score = np.where(allowed, pot + iter_vec, -np.inf)                # the score of a class in a cell
    best = score.argmax(axis=1)                                       # every cell takes the class with the best score
    area = np.bincount(best, minlength=2)
    diff = demand - area
    rows.append((it, iter_vec.round(3).tolist(), area.tolist()))
    if np.abs(diff).max() <= tolerance:
        break
    iter_vec += diff * factor                                         # push the classes that fell short up, the others down
print(pd.DataFrame(rows, columns=["iteration", "iter_vec", "cells (forest, pasture)"]).to_string(index=False))

# The competition of CLUE-S, stripped to its core: 12 cells, 2 classes (forest, pasture), pasture may not go back to forest. rng = np.random.default_rng(3) n = 12 pot = np.column_stack([rng.uniform(0.3, 0.7, n), # potential of forest in each cell np.linspace(0.95, 0.15, n) + rng.normal(0, 0.05, n)]) # potential of pasture: falls along a gradient now = np.array([0] * 10 + [1] * 2) # 10 forest cells, 2 pasture cells allowed = np.array([[True, True], [False, True]])[now] # row = class now; forest -> {forest, pasture}; pasture -> {pasture} demand = np.array([6, 6]) # the demand: 6 cells of each class tolerance, factor, iter_vec = 0, 0.005, np.zeros(2) # tolerance = max_difference, factor = factor_iteration rows = [] for it in range(200): score = np.where(allowed, pot + iter_vec, -np.inf) # the score of a class in a cell best = score.argmax(axis=1) # every cell takes the class with the best score area = np.bincount(best, minlength=2) diff = demand - area rows.append((it, iter_vec.round(3).tolist(), area.tolist())) if np.abs(diff).max() <= tolerance: break iter_vec += diff * factor # push the classes that fell short up, the others down print(pd.DataFrame(rows, columns=["iteration", "iter_vec", "cells (forest, pasture)"]).to_string(index=False))

 iteration        iter_vec cells (forest, pasture)
         0      [0.0, 0.0]                 [1, 11]
         1 [0.025, -0.025]                 [2, 10]
         2 [0.045, -0.045]                  [4, 8]
         3 [0.055, -0.055]                  [5, 7]
         4   [0.06, -0.06]                  [5, 7]
         5 [0.065, -0.065]                  [5, 7]
         6   [0.07, -0.07]                  [6, 6]

The first pass gives pasture 11 of the 12 cells, because its potential is higher almost everywhere. The correction then moves, 0.025 in each direction, and the pasture cells fall to 10, 8, 7 and, once the correction reaches 0.07, to 6 — the demand. The two pasture cells of the start never leave it, because the matrix does not let them.

Building a Continuous Model¶

With the three components opened, here they are put together in a script. PotentialLinearRegression + AllocationClueLike answer "how much does this cell's land use change" — the right choice whenever a cell can legitimately hold more than one land use at once (a partially deforested cell, a partially urbanized one). Each cell holds a fraction of each class. Potential is one linear regression per land-use type, each with its own intercept and driving-factor coefficients (RegressionSpec); allocation distributes the demand across cells according to that potential, subject to per-class bounds (AllocationSpec).

The example below is the package's own synthetic scenario, run in this notebook: a 30 × 30 grid, mostly forest, with a horizontal road through the middle, a smooth terrain (slope), and one urban seed cell on the road. Two drivers, three classes:

In [12]:

from dissmodel.geo import RasterBackend

LAND_USE_TYPES = ["forest", "agriculture", "urban"]
H = W = 30
rng = np.random.default_rng(42)

road_row = H // 2
dist_road = np.repeat(np.abs(np.arange(H, dtype=np.float32)[:, None] - road_row), W, axis=1)
dist_road /= dist_road.max()                                    # 0 on the road, 1 at the edges

yy, xx = np.mgrid[0:H, 0:W]
slope = (np.sin(xx / 4.0) * 0.5 + np.cos(yy / 5.0) * 0.5 + rng.normal(0, 0.05, (H, W))).astype(np.float32)
slope = (slope - slope.min()) / (slope.max() - slope.min())     # normalised to [0, 1]

backend = RasterBackend(shape=(H, W))
backend.set("dist_road", dist_road.astype(np.float32))
backend.set("slope", slope)
urban = np.zeros((H, W), dtype=np.float32); urban[road_row, W // 2] = 1.0
backend.set("urban", urban)
backend.set("agriculture", np.zeros((H, W), dtype=np.float32))
backend.set("forest", 1.0 - urban)

# Demand: total area per class and step (order of LAND_USE_TYPES); urban and agriculture grow
N_STEPS = 8
demand_matrix = [[float(H * W - (1 + 3 * s) - s), float(s), float(1 + 3 * s)] for s in range(N_STEPS)]
pd.DataFrame(demand_matrix, columns=LAND_USE_TYPES).rename_axis("step").T

from dissmodel.geo import RasterBackend LAND_USE_TYPES = ["forest", "agriculture", "urban"] H = W = 30 rng = np.random.default_rng(42) road_row = H // 2 dist_road = np.repeat(np.abs(np.arange(H, dtype=np.float32)[:, None] - road_row), W, axis=1) dist_road /= dist_road.max() # 0 on the road, 1 at the edges yy, xx = np.mgrid[0:H, 0:W] slope = (np.sin(xx / 4.0) * 0.5 + np.cos(yy / 5.0) * 0.5 + rng.normal(0, 0.05, (H, W))).astype(np.float32) slope = (slope - slope.min()) / (slope.max() - slope.min()) # normalised to [0, 1] backend = RasterBackend(shape=(H, W)) backend.set("dist_road", dist_road.astype(np.float32)) backend.set("slope", slope) urban = np.zeros((H, W), dtype=np.float32); urban[road_row, W // 2] = 1.0 backend.set("urban", urban) backend.set("agriculture", np.zeros((H, W), dtype=np.float32)) backend.set("forest", 1.0 - urban) # Demand: total area per class and step (order of LAND_USE_TYPES); urban and agriculture grow N_STEPS = 8 demand_matrix = [[float(H * W - (1 + 3 * s) - s), float(s), float(1 + 3 * s)] for s in range(N_STEPS)] pd.DataFrame(demand_matrix, columns=LAND_USE_TYPES).rename_axis("step").T

Out[12]:

step	0	1	2	3	4	5	6	7
forest	899.0	895.0	891.0	887.0	883.0	879.0	875.0	871.0
agriculture	0.0	1.0	2.0	3.0	4.0	5.0	6.0	7.0
urban	1.0	4.0	7.0	10.0	13.0	16.0	19.0	22.0

The model is three components in one Environment, built in the order they run each step — demand, then potential, then allocation:

In [13]:

import contextlib, io
from dissmodel.core import Environment
from disslucc import DemandInline, PotentialLinearRegression, AllocationClueLike
from disslucc.schemas import RegressionSpec, AllocationSpec

env = Environment(end_time=N_STEPS - 1)

demand = DemandInline(values=demand_matrix, land_use_types=LAND_USE_TYPES)

potential = PotentialLinearRegression(
    backend=backend, demand=demand, land_use_types=LAND_USE_TYPES,
    potential_data=[[                                                        # one region
        RegressionSpec(const=-0.2, betas={"slope": 0.4}),                    # forest
        RegressionSpec(const=0.4, betas={"dist_road": -0.1, "slope": -0.5}), # agriculture
        RegressionSpec(const=0.3, betas={"dist_road": -0.6, "slope": -0.3}), # urban
    ]],
)

allocation = AllocationClueLike(
    backend=backend, demand=demand, potential=potential, land_use_types=LAND_USE_TYPES,
    static={"forest": 0, "agriculture": -1, "urban": -1},   # forest is free; the others follow demand
    complementar_lu="forest",                                # forest takes whatever the others leave
    cell_area=1.0, max_difference=5.0,                       # stop when within 5 cells of the demand
    allocation_data=[AllocationSpec(static=0), AllocationSpec(static=-1), AllocationSpec(static=-1)],
)

with contextlib.redirect_stdout(io.StringIO()):              # the components log every iteration
    env.run()

for i, lu in enumerate(LAND_USE_TYPES):
    print(f"{lu:>12}: demand {demand_matrix[-1][i]:7.1f}   allocated {float(backend.get(lu).sum()):7.1f}")

import contextlib, io from dissmodel.core import Environment from disslucc import DemandInline, PotentialLinearRegression, AllocationClueLike from disslucc.schemas import RegressionSpec, AllocationSpec env = Environment(end_time=N_STEPS - 1) demand = DemandInline(values=demand_matrix, land_use_types=LAND_USE_TYPES) potential = PotentialLinearRegression( backend=backend, demand=demand, land_use_types=LAND_USE_TYPES, potential_data=[[ # one region RegressionSpec(const=-0.2, betas={"slope": 0.4}), # forest RegressionSpec(const=0.4, betas={"dist_road": -0.1, "slope": -0.5}), # agriculture RegressionSpec(const=0.3, betas={"dist_road": -0.6, "slope": -0.3}), # urban ]], ) allocation = AllocationClueLike( backend=backend, demand=demand, potential=potential, land_use_types=LAND_USE_TYPES, static={"forest": 0, "agriculture": -1, "urban": -1}, # forest is free; the others follow demand complementar_lu="forest", # forest takes whatever the others leave cell_area=1.0, max_difference=5.0, # stop when within 5 cells of the demand allocation_data=[AllocationSpec(static=0), AllocationSpec(static=-1), AllocationSpec(static=-1)], ) with contextlib.redirect_stdout(io.StringIO()): # the components log every iteration env.run() for i, lu in enumerate(LAND_USE_TYPES): print(f"{lu:>12}: demand {demand_matrix[-1][i]:7.1f} allocated {float(backend.get(lu).sum()):7.1f}")

      forest: demand   871.0   allocated   866.3
 agriculture: demand     7.0   allocated    11.3
       urban: demand    22.0   allocated    22.3

In [14]:

fig, axes = plt.subplots(1, 3, figsize=(11, 3.6))
for ax, lu in zip(axes, LAND_USE_TYPES):
    im = ax.imshow(backend.get(lu), cmap="viridis", vmin=0, vmax=1)
    ax.set_title(lu); ax.set_axis_off()
fig.colorbar(im, ax=axes, shrink=0.8, label="share of the cell")
fig.suptitle(f"Land use after {N_STEPS} steps (synthetic scenario)")
plt.show()

fig, axes = plt.subplots(1, 3, figsize=(11, 3.6)) for ax, lu in zip(axes, LAND_USE_TYPES): im = ax.imshow(backend.get(lu), cmap="viridis", vmin=0, vmax=1) ax.set_title(lu); ax.set_axis_off() fig.colorbar(im, ax=axes, shrink=0.8, label="share of the cell") fig.suptitle(f"Land use after {N_STEPS} steps (synthetic scenario)") plt.show()

The allocated areas meet the demand of the last step within the tolerance of 5 cells (forest ends 4.7 cells below its demand and agriculture 4.3 above — a small gap in absolute terms, but more than half of agriculture's own demand of 7, a reminder that a tolerance set for the whole grid says little about a small class). Urban land grows out from the seed along the road, where dist_road is small and its coefficient is most negative; agriculture spreads thinly around the same corridor, and forest gives up what the other two take, as complementar_lu says. Nothing here is calibrated: the coefficients were chosen to make the mechanism visible.

The two continuous variants that matter for real models were opened above: the spatial-lag potential (a term ro times the mean of the class in the neighbouring cells, so that change is drawn towards where the class already is) and the saturation allocation (a limit on how fast a cell can change once its neighbourhood passes a threshold). LuccME-BR, the national model behind Chapter 31, uses both.

A Whole Lab: LuccME's lab01¶

Lab01 is the simplest of the 21 labs, and a good first model: Acre, three classes, a demand table, a linear potential and a CLUE-like allocation. The cell below writes it with DisSLUCC. Every setting comes from the Lua script of the lab: the class names and the no-data class (outros), the demand table of the seven years 2008 to 2014, the coefficients of the three regressions, and the allocation (maxDifference, complementarLU, which class is static). The data are the Acre cells of the cell above; vector_to_raster_backend rasterizes them on a 5 km grid, and the cell area is 25 km².

The allocation is run twice, with the two tolerances discussed above: 5000 km², the value in the lab, and 1643, the value in other labs, which the benchmark uses for the variant it calls lab01_md1643. The switch cell_correction=False runs the allocation the way LuccME does (see Checked Against LuccME).

In [15]:

import contextlib, io
from dissmodel.core import Environment
from dissmodel.io.convert import vector_to_raster_backend
from disslucc import DemandPreComputedValues, PotentialLinearRegression, AllocationClueLike
from disslucc.components.demand import load_demand_csv
from disslucc.schemas import RegressionSpec, AllocationSpec

LU = ["f", "d", "outros"]
DRIVERS = ["assentamen", "uc_us", "uc_pi", "ti", "dist_riobr", "fertilidad", "rodovias"]
demand_table = load_demand_csv(fetch("examples_demand_lab1.csv").read_text(), LU)    # lab01's annual demand, 2008-2014

def run_lab01(max_difference):
    with contextlib.redirect_stdout(io.StringIO()):               # the conversion reports its coverage; not needed here
        backend = vector_to_raster_backend(source=acre, resolution=5000.0, nodata_value=-1,
                                           attrs={c: 0.0 for c in LU + DRIVERS})
    env = Environment(end_time=6)
    demand = DemandPreComputedValues(annual_demand=demand_table, land_use_types=LU)
    potential = PotentialLinearRegression(                      # lab01's PotentialCLinearRegression
        backend=backend, demand=demand, land_use_types=LU, land_use_no_data="outros",
        potential_data=[[
            RegressionSpec(const=0.7392, betas={"assentamen": -0.2193, "uc_us": 0.1754, "uc_pi": 0.09708,
                                                "ti": 0.1207, "dist_riobr": 0.0000002388, "fertilidad": -0.1313}),   # f
            RegressionSpec(const=0.267, betas={"rodovias": -0.0000009922, "assentamen": 0.2294, "uc_us": -0.09867,
                                               "dist_riobr": -0.0000003216, "fertilidad": 0.1281}),                  # d
            RegressionSpec(const=0.0),                                                                               # outros
        ]])
    allocation = AllocationClueLike(                            # lab01's AllocationCClueLike
        backend=backend, demand=demand, potential=potential, land_use_types=LU,
        static={"f": -1, "d": -1, "outros": 1}, complementar_lu="f", cell_area=25.0,
        max_difference=max_difference, cell_correction=False,   # False: the way LuccME runs it (see below)
        allocation_data=[AllocationSpec(static=-1), AllocationSpec(static=-1), AllocationSpec(static=1)])
    start_d = backend.get("d").copy()
    with contextlib.redirect_stdout(io.StringIO()):
        env.run()
    return backend, start_d, allocation, demand

backend, start_d, allocation, demand = run_lab01(max_difference=5000)      # lab01 as shipped
mask = backend.get("mask").astype(bool)
print("iterations per year, maxDifference 5000:", allocation.iterations_per_step)
for i, lu in enumerate(LU):
    print(f"{lu:>7}: demand {demand_table[-1][i]:9.1f} km²   allocated {float(backend.get(lu)[mask].sum()) * 25:9.1f} km²")

gain = (backend.get("d") - start_d)[mask]
print(f"cells that gain more than 1% of their area as deforested: {int((gain > 0.01).sum()):,} of {int(mask.sum()):,}; largest gain: {gain.max():.2f}")

backend_tight, _, allocation_tight, _ = run_lab01(max_difference=1643)    # the variant the benchmark calls lab01_md1643
print("iterations per year, maxDifference 1643:", allocation_tight.iterations_per_step)
print(f"deforested area in 2014: {float(backend_tight.get('d')[mask].sum()) * 25:9.1f} km² (demand {demand_table[-1][1]:.1f})")

change = np.where(mask, backend.get("d") - start_d, np.nan)
fig, axes = plt.subplots(1, 3, figsize=(12, 3.4))
for ax, arr, title, kw in [(axes[0], np.where(mask, start_d, np.nan), "Deforested share, 2008", dict(cmap="YlOrBr", vmin=0, vmax=1)),
                           (axes[1], np.where(mask, backend.get("d"), np.nan), "Deforested share, 2014 (simulated)", dict(cmap="YlOrBr", vmin=0, vmax=1)),
                           (axes[2], change, "Change, 2008 to 2014", dict(cmap="Reds", vmin=0, vmax=0.15))]:
    im = ax.imshow(arr, **kw); ax.set_title(title); ax.set_axis_off(); fig.colorbar(im, ax=ax, shrink=0.7)
plt.tight_layout()
plt.show()

import contextlib, io from dissmodel.core import Environment from dissmodel.io.convert import vector_to_raster_backend from disslucc import DemandPreComputedValues, PotentialLinearRegression, AllocationClueLike from disslucc.components.demand import load_demand_csv from disslucc.schemas import RegressionSpec, AllocationSpec LU = ["f", "d", "outros"] DRIVERS = ["assentamen", "uc_us", "uc_pi", "ti", "dist_riobr", "fertilidad", "rodovias"] demand_table = load_demand_csv(fetch("examples_demand_lab1.csv").read_text(), LU) # lab01's annual demand, 2008-2014 def run_lab01(max_difference): with contextlib.redirect_stdout(io.StringIO()): # the conversion reports its coverage; not needed here backend = vector_to_raster_backend(source=acre, resolution=5000.0, nodata_value=-1, attrs={c: 0.0 for c in LU + DRIVERS}) env = Environment(end_time=6) demand = DemandPreComputedValues(annual_demand=demand_table, land_use_types=LU) potential = PotentialLinearRegression( # lab01's PotentialCLinearRegression backend=backend, demand=demand, land_use_types=LU, land_use_no_data="outros", potential_data=[[ RegressionSpec(const=0.7392, betas={"assentamen": -0.2193, "uc_us": 0.1754, "uc_pi": 0.09708, "ti": 0.1207, "dist_riobr": 0.0000002388, "fertilidad": -0.1313}), # f RegressionSpec(const=0.267, betas={"rodovias": -0.0000009922, "assentamen": 0.2294, "uc_us": -0.09867, "dist_riobr": -0.0000003216, "fertilidad": 0.1281}), # d RegressionSpec(const=0.0), # outros ]]) allocation = AllocationClueLike( # lab01's AllocationCClueLike backend=backend, demand=demand, potential=potential, land_use_types=LU, static={"f": -1, "d": -1, "outros": 1}, complementar_lu="f", cell_area=25.0, max_difference=max_difference, cell_correction=False, # False: the way LuccME runs it (see below) allocation_data=[AllocationSpec(static=-1), AllocationSpec(static=-1), AllocationSpec(static=1)]) start_d = backend.get("d").copy() with contextlib.redirect_stdout(io.StringIO()): env.run() return backend, start_d, allocation, demand backend, start_d, allocation, demand = run_lab01(max_difference=5000) # lab01 as shipped mask = backend.get("mask").astype(bool) print("iterations per year, maxDifference 5000:", allocation.iterations_per_step) for i, lu in enumerate(LU): print(f"{lu:>7}: demand {demand_table[-1][i]:9.1f} km² allocated {float(backend.get(lu)[mask].sum()) * 25:9.1f} km²") gain = (backend.get("d") - start_d)[mask] print(f"cells that gain more than 1% of their area as deforested: {int((gain > 0.01).sum()):,} of {int(mask.sum()):,}; largest gain: {gain.max():.2f}") backend_tight, _, allocation_tight, _ = run_lab01(max_difference=1643) # the variant the benchmark calls lab01_md1643 print("iterations per year, maxDifference 1643:", allocation_tight.iterations_per_step) print(f"deforested area in 2014: {float(backend_tight.get('d')[mask].sum()) * 25:9.1f} km² (demand {demand_table[-1][1]:.1f})") change = np.where(mask, backend.get("d") - start_d, np.nan) fig, axes = plt.subplots(1, 3, figsize=(12, 3.4)) for ax, arr, title, kw in [(axes[0], np.where(mask, start_d, np.nan), "Deforested share, 2008", dict(cmap="YlOrBr", vmin=0, vmax=1)), (axes[1], np.where(mask, backend.get("d"), np.nan), "Deforested share, 2014 (simulated)", dict(cmap="YlOrBr", vmin=0, vmax=1)), (axes[2], change, "Change, 2008 to 2014", dict(cmap="Reds", vmin=0, vmax=0.15))]: im = ax.imshow(arr, **kw); ax.set_title(title); ax.set_axis_off(); fig.colorbar(im, ax=ax, shrink=0.7) plt.tight_layout() plt.show()

iterations per year, maxDifference 5000: [0, 0, 0, 0, 0, 0, 0]
      f: demand  136253.4 km²   allocated  133330.7 km²
      d: demand   21607.4 km²   allocated   24530.1 km²
 outros: demand    6489.2 km²   allocated    6489.2 km²
cells that gain more than 1% of their area as deforested: 3,495 of 6,574; largest gain: 0.26

iterations per year, maxDifference 1643: [0, 0, 8, 26, 18, 17, 17]
deforested area in 2014:   22608.8 km² (demand 21607.4)

What the run shows:

• With the lab's own tolerance the loop never runs. The iteration counts are all 0: the first allocation of every year is accepted. The final deforested area, 24,530 km², is 2,923 km² above the 21,607 km² demanded, and forest is 2,923 km² below its demand — a gap that is within the 5,000 km² tolerance, and about 14% of the area demanded for deforestation. The map is the potential scaled by the initial elasticity of 0.1, and no more: the feedback loop of the previous section does nothing.
• With 1643 it does. The allocation runs 8, 26, 18, 17 and 17 times in the last five years, and the deforested area ends at 22,609 km², 1,001 km² above the demand instead of 2,923. A smaller tolerance is what makes the comparison of iteration counts with TerraME informative, which is why the benchmark adds this variant.
• The change is widespread, and strongest in the east. More than half of the cells (3,495 of 6,574) gain over 1% of their area as deforestation, and the largest gains, up to about a quarter of a cell, are in the east of the state, where the drivers favour it. This is a shape of the model, not an error: the potential is a smooth regression, and an allocation that scales it by a factor spreads the change in proportion to it.

The shapefile also holds the observed 2014 shares (d2014). Comparing them with the simulated map would be a validation against observation, which is a different question from this chapter's comparison with LuccME; it is the subject of the last section.

Building a Discrete Model: the Executor Route¶

PotentialDLogisticRegression + AllocationDClueSLike answer a different question: "which single land use dominates this cell" — a CLUE-S-style allocation, right whenever the data are themselves categorical (a classified land-cover map, not a fractional-cover raster). Potential comes from a logistic regression predicting which class a cell most likely belongs to (LogisticRegressionSpec, which adds an elasticity term — how hard a class is to remove from a cell); allocation runs a competition between classes, governed by a transition matrix ([region][from][to]: which changes are allowed at all).

The executor route reads the same configuration from a TOML file — the format Chapter 33's platform registers — and runs it from the command line. The package's example for the Moju region:

[model]
executor_module   = "disslucc.executors"
name              = "lucc_discrete"
land_use_types    = ["f", "d", "o"]
cell_area         = 1.0
transition_matrix = [[[1, 1, 0], [0, 1, 0], [0, 0, 1]]]   # deforestation is irreversible

[model.parameters]
n_steps          = 6
max_difference   = 100.0
factor_iteration = 0.0001

[[model.potential_data]]
lu         = "f"
const      = -2.34187976925989
elasticity = 0.0
  [model.potential_data.betas]
  dist_br = 3.10319957497883

git clone https://github.com/DisSModel/disslucc.git && cd disslucc
pip install -e ".[examples]"
python -m disslucc.executors.discrete run \
  --toml examples/dissmodel-configs/lucc_discrete.toml \
  --input data/input/cs_moju.zip \
  --param demand_csv=data/input/demand_moju.csv \
  --output outputs/result.tif

--output writes a GeoTIFF with one band per land-use class, georeferenced from the input, and the run leaves an ExperimentRecord with the checksums of the input and the output. One caveat comes with this route: the executor rasterizes the input layer at the resolution given in the TOML, so its grid is close to, but not cell for cell the same as, the original cellular space (5,842 of the 5,914 cells of Moju). That is why the comparison with LuccME below builds each model by hand on the lab's own cells, not through the executor.

Choosing Between Continuous and Discrete¶

Both styles implement the same Demand/Potential/Allocation loop with different algorithms at each step — the decision between them is about your data, not architecture:

	Continuous	Discrete
A cell holds	a fraction of each class	one class
Potential	linear regression (optionally with a spatial lag)	logistic regression
Allocation	CLUE-like, per-class bounds, optional saturation	CLUE-S competition, transition matrix
Natural comparison with a reference	tolerance on fractions	agreement of classes, cell by cell

The fastest way to decide: look at your calibration and validation data first. If they are a class label per cell, discrete allocation is what can be checked against them directly. If they are an area or a share per cell, continuous allocation represents them without a lossy discretization forced on them first.

Checked Against LuccME¶

Writing LuccME's components again in another language is only useful if they compute what the originals compute. The package does not check that itself: a separate repository, disslucc-benchmark, runs DisSLUCC on the test labs that ship with LuccME 3.1 and compares the result, cell by cell and year by year, with the outputs TerraME 2.0.1 produced for the same labs, archived in luccme-goldens. Chapter 30 explains how those reference outputs were produced and checked; here are the results.

LuccME 3.1 ships 21 test labs. DisSLUCC implements every component of ten of them — lab01 to lab07 (continuous, Acre, 6,574 cells) and lab15 to lab17 (discrete, Moju, 5,914 cells) — and those ten match the TerraME outputs with the same number of allocation iterations in every year and a largest cell difference between 5×10⁻¹³ and 3×10⁻⁶, the larger values coming from single-precision arithmetic. The other eleven labs use components DisSLUCC does not have yet: seven lack one component each, four need several.

The 21 labs, read as a map of the components. Each row is a lab; the columns name its three components, without their common prefix (the C or D in LuccME's names — PotentialCSpatialLagRegression, AllocationDClueSLike — is the continuous or discrete style, and DisSLUCC drops it from the continuous ones). The last column says whether DisSLUCC has all its components.

Lab	Cells	Style	Demand	Potential	Allocation	In DisSLUCC
01	Acre, 5 km	continuous	PreComputedValues	LinearRegression	ClueLike	yes
02	Acre	continuous	PreComputedValues	SpatialLagRegression	ClueLike	yes
03	Acre	continuous	PreComputedValues	SpatialLagRegression	ClueLikeSaturation	yes
04	Acre	continuous	ComputeTwoDates	SpatialLagRegression	ClueLike	yes
05	Acre	continuous	ComputeThreeDates	SpatialLagRegression	ClueLike	yes
06	Acre	continuous	PreComputedValues	SpatialLagRegression	ClueLikeSaturation	yes
07	Acre	continuous	PreComputedValues	SpatialLagRegression	ClueLike	yes
08	Acre	continuous	PreComputedValues	SpatialLagLinearRegressionMix	ClueLike	no: the potential
09	Acre	continuous	PreComputedValues	SampleBased	ClueLike	no: the potential
10	Moju, 0.0045°	discrete	PreComputedValues	NeighSimpleRule	SimpleOrdering	no: both
11	Moju	discrete	PreComputedValues	InverseDistanceRule	SimpleOrdering	no: both
12	Moju	discrete	PreComputedValues	NeighInverseDistanceRule	SimpleOrdering	no: both
13	Moju	discrete	PreComputedValues	LogisticRegression	SimpleOrdering	no: the allocation
14	Moju	discrete	PreComputedValues	LogisticRegressionNeighAttract	ClueSLike	no: the potential
15	Moju	discrete	PreComputedValues	LogisticRegression	ClueSLike	yes
16	Moju	discrete	ComputeTwoDates	LogisticRegression	ClueSLike	yes
17	Moju	discrete	ComputeThreeDates	LogisticRegression	ClueSLike	yes
18	Moju	discrete	PreComputedValues	SampleBased	ClueSLike	no: the potential
19	Moju	discrete	ComputeTwoDates	LogisticRegressionNeighAttractRepulsion	ClueSNeighOrdering	no: both
20	Moju	discrete	ComputeTwoDates	LogisticRegressionNeighAttractRepulsion	ClueSLike	no: the potential
21	Moju	discrete	PreComputedValues	LogisticRegression	ClueSNeighOrdering	no: the allocation

Reading across the table is a good way to see what each component is for. Labs 02 and 03 differ only in the allocation, so what separates them is the saturation brake; labs 04 and 05 differ only in the demand, two dates against three; labs 15, 16 and 17 are the discrete version of the same comparison. The table follows what each lab's script instantiates; the comment at the top of some of them (labs 09, 18 and 20) names a different component from the one the code builds.

A continuous and a discrete LuccME lab, DisSLUCC 0.5.0 (top row) against the TerraME 2.0.1 output (bottom row), on the same cells and with the same color limits. Left, lab01: the deforested fraction in 2014 and its change since 2008. Right, lab15_md10: the class of each cell in 2004 and the 237 cells that changed class since 1999. The two rows are not expected to be distinguishable by eye. From disslucc-benchmark, docs/fig_models.png.

Three details of this comparison are worth more than the headline.

The labs barely exercise the convergence loop. In the ten labs as shipped, the allocation is accepted at the first pass every year, so "same number of iterations" says little. The benchmark adds two variants with a smaller maxDifference, in which the loop runs 8 to 67 times per year: lab01_md1643 (continuous) and lab15_md10 (discrete). DisSLUCC reproduces TerraME's iteration count in every year of both. The second answers an old caveat of this chapter: with its original maxDifference of 300, lab15 is nearly non-discriminative — no cell changes class, and a static ranking by the difference of two probabilities already reproduces it — so its match only showed that the coefficients were transcribed correctly. In lab15_md10, 237 cells change class, the competition between classes matters, and every cell still agrees.

One difference is deliberate. LuccME's continuous allocation contains a correction of the change in each cell that never runs: the condition that guards it reads a field named regionregionAloc instead of regionAloc, so it is always false. DisSLUCC runs the correction by default, as the algorithm intends, and has a switch, AllocationClueLike(cell_correction=False), that reproduces LuccME exactly. With the correction on, lab01 with maxDifference 1643 differs from TerraME by a mean of 0.0036 in the forest and deforested fractions in 2014 (largest 0.027), and the iteration counts change. An earlier draft of this chapter reported that 0.0036 as the agreement of DisSLUCC with LuccME and blamed it on the convergence tolerance; it is in fact entirely this correction, and with the switch off the difference disappears.

Parity is not validation of a model. All of this shows that DisSLUCC's components compute what LuccME's compute, on the labs, to the precision measured. It does not show that the labs are good models of Acre or Moju, nor that LuccME itself is right. Those are the questions of the next section.

Calibration, Validation, and Goodness-of-Fit¶

The comparison above is verification — confirming that a Python implementation computes what a reference run computes (Chapter 30). That is a different question from validation: fitting a model to real-world observed data by splitting a timeline into a calibration period (used to fit the parameters) and a held-out validation period, checked only after fitting.

A cell-by-cell comparison between a simulated map and an observed one is often too strict a bar even for scientific validation — two reasonable model runs can disagree pixel for pixel while still being "equally good" at the pattern level. A multiscale comparison — checking agreement in successively larger windows (3×3, then 5×5, then 9×9, and so on) — is the standard alternative: two maps might show weak agreement at the finest resolution but strong agreement at a coarser one, which is real, usable information a strict pixel match would have discarded. For categorical maps, the Pontius & Millones decomposition into quantity and allocation disagreement — available in disslucc.validation.pontius_millones — says how two maps disagree, which a single accuracy number hides.

Two more checks keep a validation honest. The first is a null model: persistence, the map in which nothing changes, is hard to beat at the resolution of the cell, because most cells do not change in a few years, and a model has merit only where it does better than that. In a comparison of thirteen applications, the error was larger than the correctly predicted change in twelve of them at the resolution of the data, and the coarser resolution at which a model starts to beat persistence, the null resolution, was one of the two statistics used to compare them (Pontius et al., 2008). The second is the figure of merit, the cells where both maps show change divided by the cells where either does, which counts only change and is not inflated by the cells that neither map changes. Both are a few lines of NumPy on two class arrays.

Exercises¶

1. Match the component to the question. For each research question, name which allocation style fits and why: (a) "how does the percentage of forest cover in each cell change over the next 20 years," (b) "which cells convert from forest to pasture by 2030."
2. Why two different regressions? PotentialLinearRegression and PotentialDLogisticRegression both take driving-factor coefficients (betas) per land-use type. From the class names alone, explain why one needs a linear regression and the other a logistic one — what is each one actually trying to predict?
3. Move the road. In the synthetic scenario, change the urban coefficient of dist_road from −0.6 to +0.6 and run it again. Where does urban land go now, and does the allocation still meet the demand?
4. Tolerance and iterations. Run the scenario with max_difference=0.5 instead of 5.0. Does the allocated area get closer to the demand? Why does a smaller tolerance make a comparison of iteration counts with a reference more informative?
5. Read a deliberate difference. In your own words, explain why the benchmark keeps lab01/cell_correction as a scenario that is expected to differ from TerraME, instead of removing it.
6. Two dates or three. Using the cell in Demand, compute the demand of forest (f) with both. In which year do they differ most, and why is it the middle year?
7. How far does a neighbour reach. In the spatial-lag cell, try ro equal to 0.3, 0.8 and 1.5. What ro makes the cells that touch the clearing ask for half a share? What goes wrong with a value that is too large?
8. A harder demand. In the CLUE-like toy, raise the demand from 5 to 8. How many iterations does the loop need, and why does each iteration gain less as the loop goes on?
9. Read two labs. Choose two rows of the table of labs that differ in a single component and say what question the difference answers (the saturation brake, the number of dates of the demand, the allocation style).
10. Tolerance in a whole lab. In the lab01 cell, run max_difference=500. How many iterations does the loop take per year, and how far from the demand does the allocation end? Compare with 5000 and 1643.

In [16]:

# Your code here

# Your code here

Summary¶

Key concepts introduced¶

• Land use versus land cover; a land change as a transition of area between classes; LUCC models as projections conditioned on a demand that the user supplies, mostly used to compare scenarios and not to forecast
• The families of land change models, and the hybrid top-down / bottom-up family of CLUE and LuccME that DisSLUCC implements
• The Demand/Potential/Allocation loop, and what each component computes: the interpolated demand of two and three dates; the linear, spatial-lag and logistic potentials; the elasticity feedback of CLUE-like allocation, its saturation variant, and the competition of CLUE-S-like allocation
• Continuous cells (shares of each class) versus discrete cells (one class), and how to choose between them by the type of your data
• The LuccME test labs as two real cases, Acre and Moju, and the table of the 21 labs as a map of the components, with the ten that DisSLUCC covers
• A continuous model built and run in a script and as a whole real lab; the executor and TOML route to the same components, with its command line and its provenance record
• The comparison with ten of the 21 LuccME labs: identical iteration counts, differences at single-precision level, two variants that make the convergence loop and the class competition matter, and one deliberate difference (the cell correction LuccME never runs)
• Verification against a reference implementation versus validation against observations, with a null model and the figure of merit

Chapter 27 stays in domain-application territory but moves from land far inland to the coastline itself — a coupled flood and mangrove-migration model, run on both substrates and validated two different ways at once.

Further Reading¶

• Costa, S. S. (2026). DisSLUCC (Version 0.5.0) [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.23219338
• Costa, S. S. (2026). DisSLUCC benchmark (Version 0.2.2) [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.23224097 — the comparison with the LuccME labs
• Veldkamp, A., & Fresco, L. O. (1996). CLUE: a conceptual model to study the conversion of land use and its effects. Ecological Modelling, 85(2–3), 253–270
• Verburg, P. H. et al. (2002). Modeling the spatial dynamics of regional land use: the CLUE-S model. Environmental Management, 30(3), 391–405
• Costanza, R. (1989). Model goodness of fit: a multiple resolution procedure. Ecological Modelling, 47(3–4), 199–215 — the multiscale comparison
• Pontius Jr., R. G. et al. (2008). Comparing the input, output, and validation maps for several models of land change. The Annals of Regional Science, 42(1), 11–37. https://doi.org/10.1007/s00168-007-0138-2 — persistence as the null model, and the figure of merit
• Pontius Jr., R. G. & Millones, M. (2011). Death to Kappa: birth of quantity disagreement and allocation disagreement for accuracy assessment. International Journal of Remote Sensing, 32(15), 4407–4429
• Aguiar, A. P. D. (2006). Modelagem de mudança do uso da terra na Amazônia: explorando a heterogeneidade intra-regional. Tese (Doutorado em Computação Aplicada), INPE, São José dos Campos — the regional model, with demand and potential kept apart
• Veldkamp, A., & Lambin, E. F. (2001). Predicting land-use change. Agriculture, Ecosystems & Environment, 85(1–3), 1–6 — proximate causes versus underlying driving forces
• Briassoulis, H. (2000). Analysis of land use change: theoretical and modeling approaches. The Web Book of Regional Science, West Virginia University — a review of the theories
• Verburg, P. H., de Koning, G. H. J., Kok, K., Veldkamp, A., & Bouma, J. (1999). A spatial explicit allocation procedure for modelling the pattern of land use change based upon actual land use. Ecological Modelling, 116(1), 45–61 — the allocation procedure behind the continuous components
• Aguiar, A. P. D., Carneiro, T., Andrade, P., & Assis, T. (2012). LuccME-TerraME: An open-source framework for spatially explicit land use change modelling. GLP News, 8, 21–23
• Carneiro, T. G. S., Andrade, P. R., Câmara, G., Monteiro, A. M. V., & Pereira, R. R. (2013). An extensible toolbox for modeling nature–society interactions. Environmental Modelling & Software, 46, 104–117. https://doi.org/10.1016/j.envsoft.2013.03.002 — the TerraME platform under LuccME
• Bezerra, F. G. S., Von Randow, C., Assis, T. O., Bezerra, K. R. A., Tejada, G., Castro, A. A., Gomes, D. M. P., Avancini, R., & Aguiar, A. P. D. (2022). New land-use change scenarios for Brazil: Refining global SSPs with a regional spatially-explicit allocation model. PLOS ONE, 17(4), e0256052. https://doi.org/10.1371/journal.pone.0256052 — LuccME-BR
• LuccME documentation (INPE): http://www.dpi.inpe.br/luccme/