Land change simulation

Land use and land cover change models: analysing the changes between maps, estimating where change is likely, allocating future change and validating the simulated maps: the ls.simulate module. It brings the tools of the field together as interchangeable pieces of one model (CLUE-S, CLUMondo, Dinamica EGO, PLUS, FUTURES), each reproducing its reference. The module is being built phase by phase; this page grows with it.

import landschaft as ls

lu = ls.load_raster(ls.simulate.example_data("pie_landuse"))   # 1985, 1991, 1999
ch = ls.simulate.changes(lu)                                   # transitions, areas, change counts
ia = ch.intensity(gain=2, loss=1)                              # intensity analysis
ls.simulate.difference_components(ch.matrix())                 # quantity, exchange, shift
mk = ch.markov(steps=14)                                       # transition rates per year

ef = ls.load_raster(ls.simulate.example_data("pie_drivers"))
woe = ls.simulate.weights_of_evidence(lu[0], lu[2], ef, {"dist_built_1985": 100})
p = woe.probability(lu[2], ef)                                 # Dinamica EGO's transition probabilities

forms = {1: ["elevation", "slope"], 2: ["elevation", "slope", "dist_built_1985"], 3: ["elevation", "slope"]}
model = ls.simulate.Model(                                     # CLUE-S, as lulcc
    potential=ls.simulate.LogisticRegression(forms).fit(lu[0], ef),
    demand=ls.simulate.AreaDemand.from_maps(lu, times=[0, 6, 14]),
    allocation=ls.simulate.CLUE(jitter=0.0002, scale=1e-6, max_diff=50, mean_diff=50),
    rules=ls.simulate.Rules(elasticity=[0.2, 0.2, 0.2]))
sim = model.run(lu[0], ef, seed=1)                             # 14 yearly maps

Land systems that supply goods and services (CLUMondo), given as overlapping binary layers:

S = ls.simulate
st = S.LandSystemStack({"forest": f, "extensive": pe, "intensive": pi, "soy": soy})
st.overlaps()                                                  # the combinations (mosaics) found
roles = S.Roles(natural=["forest"], families={"cattle": ["extensive", "intensive"], "crops": ["soy"]})
model = S.Model(
    potential=st.suitability({"forest": sf, "extensive": se, "intensive": si, "soy": ss}),
    demand=S.ServiceDemand(demand_df),                         # heads of cattle, tons of soy per year
    allocation=S.CLUMondo(st.supply({"extensive": {"cattle": 0.6}, "intensive": {"cattle": 2.5},
                                     "soy": {"soy": 3.4}})),
    rules=S.Rules(elasticity=0.5, roles=roles, conversion_order=["intensification", "expansion"],
                  caps=[S.Cap("expansion", cells=500)], restrictions=parks))
sim = model.run(st, seed=1)
sim.layers["forest"]; sim.shortfall; sim.trajectory()

Urban growth by patches (FUTURES):

demand = S.DevelopmentDemand.from_population(observed_maps, counties, observed_pop, projected_pop,
                                             range(2020, 2051))
model = S.Model(
    potential=S.DevelopmentPotential(coefficients_df),         # intercept, devpressure, predictors per county
    demand=demand,
    allocation=S.PatchGrowing(patch_library, subregions=counties, neighborhood=10,
                              compactness_mean=0.4, compactness_range=0.05, incentive_power=2))
sim = model.run(developed_2020, {"slope": slope, "road_dist": road_dist}, seed=1)
sim.development                                                # step in which each cell was developed

Change between maps

changes

changes(
    maps,
    times=None,
    *,
    cell_area: float | None = None,
    nodata=None,
    labels: dict | None = None,
    n_threads: int | None = None,
) -> Changes

Land change along a series of categorical maps.

The transitions between consecutive maps and from the first to the last, the area of each class, the intensity analysis of Aldwaik & Pontius (2012) and the number of changes per cell: the analyses of OpenLand (contingencyTable, intensityAnalysis, acc_changes), lulcc (crossTabulate, total) and diffeR (crosstabm, diffTablej), with the same numbers.

Parameters:

Name	Type	Description	Default
maps		A 3-D xarray.DataArray with a band per date (as ls.load_raster reads a multi-band file) or a list of 2-D maps on the same grid.	required
times		The date of each map (default: the band names when they are numbers, else 0, 1, 2...).	None
cell_area	float | None	Area of a cell in square map units (default: the product of the resolutions; km2 columns are cells times this over 10⁶, so give it in m²).	None
nodata		Value marking missing cells (default: the rasters' no-data value).	None
labels	dict | None	Names of the classes, {value: name}, for matrices and tables.	None
n_threads	int | None	Threads (default: all cores but one).	None

Returns:

Type	Description
Changes	With table, one_step, total_area, areas, matrix(), intensity(), components(), change_count() and markov().

Examples:

>>> import landschaft as ls
>>> lu = ls.load_raster(ls.simulate.example_data("pie_landuse"))
>>> ch = ls.simulate.changes(lu)
>>> ch.times
[1985, 1991, 1999]
>>> ia = ch.intensity(gain=2, loss=1)       # gains of built land, losses of forest
>>> ia.interval.columns.tolist()
['Period', 'PercentChange', 'St', 'U']

Changes dataclass

Changes(
    maps: list[Map],
    times: list,
    classes: ndarray,
    cell_area: float,
    labels: dict | None = None,
    _pairs: dict = dict(),
)

Land change along a series of categorical maps (see changes()).

table property

table: DataFrame

Transitions between consecutive maps, one row per pair of classes (OpenLand's lulc_Multistep).

Columns Period ("t0-t1"), From, To, km2 (cells times the cell area, over 10⁶), cells, Interval, yearFrom, yearTo; pairs without cells are left out.

one_step property

one_step: DataFrame

Transitions from the first map to the last (OpenLand's lulc_Onestep).

total_area property

total_area: dict

Area of the first period's cells: {"km2": ..., "cells": ...} (OpenLand's totalArea).

total_interval property

total_interval

From the first time to the last.

areas property

areas: DataFrame

Cells of each class (columns) in each map (rows, by time), as lulcc's total.

matrix

matrix(t0=None, t1=None) -> DataFrame

Transition matrix (cells) from the map at t0 to the map at t1 (default: first to last).

intensity

intensity(gain, loss, *, area: bool = True) -> IntensityAnalysis

Intensity analysis at the interval, category and transition levels (OpenLand's intensityAnalysis).

Parameters:

Name	Type	Description	Default
gain		The class whose gains the transition level examines (OpenLand's category_n): which classes it takes over more intensively than uniformly.	required
loss		The class whose losses the transition level examines (category_m): which classes take it over.	required
area	bool	Areas in km² (default) or in cells.	True

Returns:

Type	Description
IntensityAnalysis

change_count

change_count() -> tuple[DataArray | ndarray, DataFrame]

How many times each cell changed along the series (OpenLand's acc_changes).

Returns:

Type	Description
(map, table)	The count per cell (int32, -1 where any map is missing; a DataArray on the maps' grid when they came as rasters) and a table of PxValue, Qt (cells) and Percent.

components

components(t0=None, t1=None, *, digits: int | None = None)

Gain, persistence, loss and the quantity, exchange and shift components of the change.

difference_components() of matrix() with analysis="change".

markov

markov(t0=None, t1=None, *, steps: int = 1)

The Markov chain of the change from t0 to t1 (see Markov).

IntensityAnalysis dataclass

IntensityAnalysis(
    interval: DataFrame,
    category_gain: DataFrame,
    category_loss: DataFrame,
    transition_gain: DataFrame,
    transition_loss: DataFrame,
    category_gain_stationarity: DataFrame,
    category_loss_stationarity: DataFrame,
    transition_gain_stationarity: DataFrame,
    transition_loss_stationarity: DataFrame,
    gain: object = None,
    loss: object = None,
)

Intensity analysis of land change (Aldwaik & Pontius 2012), as OpenLand computes it.

Every rate is a percentage per time unit (per year when the times are years). Areas are in km² (area=True) or cells.

Attributes:

Name	Type	Description
interval	DataFrame	Per period: PercentChange (changed area, % of the study area), St (annual change intensity), U (uniform intensity over the whole series).
category_gain, category_loss	DataFrame	Per period and class: gross gain (GG) or loss (GL), its intensity Gtj / Lti and the period's St.
transition_gain	DataFrame	Transitions to the class gain from each other class: area T_i2n, intensity Rtin and the uniform intensity Wtn.
transition_loss	DataFrame	Transitions from the class loss to each other class: area T_m2j, intensity Qtmj and the uniform intensity Vtm.
category_gain_stationarity, category_loss_stationarity, transition_gain_stationarity,
transition_loss_stationarity	DataFrame	In how many periods each class is active / dormant (category level) or targeted / avoided (transition level), and Test = "Y" when that holds in every period.

crosstab

crosstab(
    from_map,
    to_map,
    *,
    nodata=None,
    labels: dict | None = None,
    n_threads: int | None = None,
) -> DataFrame

Cells of each class of one map that hold each class of another (a transition matrix of counts).

Rows are the classes of from_map (the earlier map, or the comparison map), columns those of to_map (the later map, or the reference); only cells with a value in both maps count. The classes are those found in such cells, in increasing order (the set diffeR::crosstabm and raster::crosstab use).

Parameters:

Name	Type	Description	Default
from_map		Categorical maps (xarray.DataArray or 2-D arrays) on the same grid.	required
to_map		Categorical maps (xarray.DataArray or 2-D arrays) on the same grid.	required
nodata		Value marking missing cells (default: the rasters' no-data value; NaN is always missing).	None
labels	dict | None	Names of the classes ({value: name}) for the rows and columns.	None
n_threads	int | None	Threads (default: all cores but one).	None

Returns:

Type	Description
DataFrame	Counts (int64), index from, columns to.

Examples:

>>> import landschaft as ls
>>> lu = ls.load_raster(ls.simulate.example_data("sibuyan_landuse"))
>>> m = ls.simulate.crosstab(lu[0], lu[1])
>>> int(m.values.sum()), m.shape
(7226, (5, 5))

Components of difference

difference_components

difference_components(
    ctmatrix, *, analysis: str = "error", digits: int | None = None
) -> DataFrame

Agreement and the components of difference of each category and overall (diffeR's diffTablej).

Parameters:

Name	Type	Description	Default
ctmatrix		A square cross-tabulation (a DataFrame from crosstab(), or an array): rows the comparison map (or the earlier map), columns the reference map (or the later map); counts, areas or percentages.	required
analysis	('error', 'change')	Names of the first three columns: Omission, Agreement, Comission (diffeR's spelling), or Gain, Persistence, Loss for change between two dates.	"error"
digits	int | None	Round to this many decimals as R does (diffeR rounds to 0 by default; here nothing is rounded unless asked).	None

Returns:

Type	Description
DataFrame	One row per category and an Overall row (the column sums, with quantity, exchange and shift halved, as each difference counts in two categories): Category, omission (gain), agreement (persistence), commission (loss), Quantity, Exchange, Shift.

Examples:

>>> import landschaft as ls
>>> lu = ls.load_raster(ls.simulate.example_data("pie_landuse"))
>>> d = ls.simulate.difference_components(ls.simulate.crosstab(lu[0], lu[2]), analysis="change")
>>> d.columns.tolist()
['Category', 'Gain', 'Persistence', 'Loss', 'Quantity', 'Exchange', 'Shift']

overall_difference

overall_difference(ctmatrix) -> dict

The overall difference and its components (diffeR's overallDiff, overallQtyD...).

Returns:

Type	Description
dict	quantity, exchange, shift, allocation (exchange plus shift) and overall (all the off-diagonal entries), in the units of the matrix.

Examples:

>>> import landschaft as ls
>>> lu = ls.load_raster(ls.simulate.example_data("sibuyan_landuse"))
>>> d = ls.simulate.overall_difference(ls.simulate.crosstab(lu[0], lu[1]))
>>> d["overall"] == d["quantity"] + d["exchange"] + d["shift"]
True

difference_resolutions

difference_resolutions(
    comp, ref, *, factor: int = 2, percent: bool = True, nodata=None, population=None
) -> DataFrame

Quantity, exchange and shift differences at increasingly coarse resolutions (diffeR's differenceMR).

At the original resolution the cross-tabulation of the two maps; at factor, factor²... cells per side (up to the larger side of the map) the composite cross-tabulation (composite_matrix()). Shift shrinks as the resolution coarsens; quantity does not change.

Parameters:

Name	Type	Description	Default
comp		The comparison and the reference map (same grid). diffeR gives NaN at the coarser resolutions when a map has missing cells (see composite_matrix()); here they are left out of the blocks.	required
ref		The comparison and the reference map (same grid). diffeR gives NaN at the coarser resolutions when a map has missing cells (see composite_matrix()); here they are left out of the blocks.	required
factor	int	Growth of the block side from one resolution to the next.	2
percent	bool	Differences in % of the domain (default) or as proportions.	True
population		Population of each category (see sample_to_population()) to scale the original-resolution matrix.	None

Returns:

Type	Description
DataFrame	Multiples, Resolution (block side in map units), Quantity, Exchange, Shift, Overall.

composite_matrix

composite_matrix(comp, ref, factor: int, *, nodata=None) -> DataFrame

Cross-tabulation at a coarser resolution by the composite method (diffeR's composite).

Each block of factor x factor cells is a mix of categories; the agreement of a category in a block is the smaller of its two shares, and the rest is distributed among the pairs of categories as Pontius & Cheuk (2006) propose. Blocks with a missing cell or cut by the edge are left out.

Parameters:

Name	Type	Description	Default
comp		The comparison and the reference map. The categories are those of either map, matched by value (diffeR counts a missing value as a category, which leaves a row and a column of NA, and names the rows in the order the values first appear).	required
ref		The comparison and the reference map. The categories are those of either map, matched by value (diffeR counts a missing value as a category, which leaves a row and a column of NA, and names the rows in the order the values first appear).	required
factor	int	Cells per side of a block.	required

Returns:

Type	Description
DataFrame	Proportions of the domain (rows comp, columns ref).

memberships

memberships(map, factor: int = 2, *, nodata=None) -> dict

The share of each category in each block of factor x factor cells (diffeR's memberships).

terra::segregate then aggregate(fun = mean): blocks with a missing cell, and the partial blocks at the right and bottom edges, are NaN.

Returns:

Type	Description
dict	{category: 2-D array}, the categories in increasing order.

mad

mad(grid1, grid2, *, strata=None, multiple: bool = True, nodata=None) -> DataFrame

Mean absolute deviation between two continuous maps and its components (diffeR's MAD).

The deviation grid1 - grid2 is split into quantity (the mean deviation), strata (deviations that cancel within each stratum) and element (the rest, which coarser resolutions remove); multiple repeats it at blocks of 2, 4, 8... cells (sums of the deviations).

Parameters:

Name	Type	Description	Default
grid1		Continuous maps on the same grid (NaN = missing).	required
grid2		Continuous maps on the same grid (NaN = missing).	required
strata		A map of strata (integers; 0 and missing are outside every stratum); default: one stratum over the cells of grid1.	None
multiple	bool	All resolutions (default) or only the original one.	True

Returns:

Type	Description
DataFrame	Multiples (with multiple), Resolution, Quantity, Strata, Element, Total.

sample_to_population

sample_to_population(ctmatrix, population) -> DataFrame

A sample's cross-tabulation scaled to the population of each row category (diffeR's sample2pop).

Each row j becomes population[j] * row / sum(row). diffeR takes the population of the row's position instead of the row's category (a wrong result when the population table is in another order); here the category decides.

Parameters:

Name	Type	Description	Default
ctmatrix	DataFrame	Sample counts, rows and columns named by category.	required
population		The population of each row category: a two-column table (category, size), a dict or a Series.	required

Returns:

Type	Description
DataFrame

Markov chains

Markov

Markov(matrix, classes, steps: int = 1)

A Markov chain of land change: the transition probabilities of a period and of each time step.

single holds the probability that a cell of class i at the start of the period is of class j at its end (every cell counted once), multi the probabilities per time step when the period spans steps steps: the principal steps-th root of single (Dinamica EGO's DetermineTransitionMatrix, single-step and multi-step matrices).

Create one with from_maps() (or Changes.markov), or from a matrix of probabilities.

Parameters:

Name	Type	Description	Default
matrix		Transition probabilities of the period (rows sum to 1).	required
classes		The class of each row and column.	required
steps	int	Time steps in the period.	1

single property

single: DataFrame

Transition probabilities over the whole period.

multi property

multi: DataFrame

Transition probabilities per time step (the principal root of single).

from_counts classmethod

from_counts(counts, classes, steps: int = 1) -> Markov

From a matrix of cell counts (rows: classes at the start); classes without cells at the start are dropped.

from_maps classmethod

from_maps(
    t0, t1, *, steps: int = 1, nodata=None, n_threads: int | None = None
) -> Markov

The chain of the change from map t0 to map t1 over steps time steps.

Examples:

>>> import landschaft as ls
>>> lu = ls.load_raster(ls.simulate.example_data("pie_landuse"))
>>> mk = ls.simulate.Markov.from_maps(lu[0], lu[2], steps=14)
>>> mk.single.shape
(3, 3)

rates

rates(*, multi: bool = True) -> DataFrame

The transition rates as Dinamica EGO lists them: off-diagonal, positive, From, To, Rate.

Negative rates (roots that are not stochastic matrices) are left out, as Dinamica does; a cell then keeps its class with the remaining probability.

project

project(areas, steps: int | None = None, *, multi: bool = True) -> DataFrame

Areas of each class after each step: areas @ P repeated.

Parameters:

Name	Type	Description	Default
areas		The area (or cells) of each class now: a sequence in the order of classes, a dict or a Series indexed by class.	required
steps	int | None	Steps to project (default: the chain's steps with multi, else 1).	None
multi	bool	Use the per-step matrix (default) or the period's.	True

Returns:

Type	Description
DataFrame	One row per step (0 = the given areas), one column per class.

Transition potential

WeightsOfEvidence

WeightsOfEvidence(
    ranges: dict | None = None,
    *,
    transitions=None,
    increment="auto",
    min_delta: float = 1.0,
    max_delta: float = 100000.0,
    tolerance: float = 5.0,
    categorical=(),
    fix_abnormal: bool = False,
    possible: str = "any",
)

Weights of evidence of land use transitions, calibrated on two maps (Dinamica EGO's method).

Parameters:

Name	Type	Description	Default
ranges	dict | None	How each explanatory variable is cut into ranges, by name: "categorical" (a range per category); a number (the increment of a continuous variable, the other settings by default); a dict with increment, min_delta, max_delta, tolerance; or a list of range limits. Variables not named are continuous with the defaults below; categorical names more categorical ones.	None
transitions		The transitions (from, to) to calibrate (default: every one observed between the maps).	None
increment		Default buffer increment of continuous variables, in the variable's units ("auto": a hundredth of its range).	'auto'
min_delta	float	Smallest and largest range, in cells (Dinamica's defaults: 1 and 100 000).	1.0
max_delta	float	Smallest and largest range, in cells (Dinamica's defaults: 1 and 100 000).	1.0
tolerance	float	Tolerance angle, degrees (default 5): a smaller angle gives more, narrower ranges.	5.0
categorical		Names of categorical variables. ranges may also be keyed by transition, {(1, 2): {"dist": 30, "soil": "categorical"}, (1, 3): {...}}: each transition with its own variables and settings (Dinamica's skeleton).	()
fix_abnormal	bool	Give ranges with no transition, or only transitions, finite weights (-ln(B N + 1) and ln(B (A - N) + 1), B the range's cells, N the transitions, A the possible cells) instead of 0.	False
possible	('any', 'stay')	Cells where a transition was possible but did not happen: every other cell of the initial class ("any", default; Dinamica's "potential transition criteria") or only those that kept their class.	"any"

Examples:

>>> import landschaft as ls
>>> lu = ls.load_raster(ls.simulate.example_data("sibuyan_landuse"))
>>> drivers = ls.load_raster(ls.simulate.example_data("sibuyan_drivers"))
>>> woe = ls.simulate.WeightsOfEvidence({"dist_roads": 250, "dist_sea": 250}, transitions=[(1, 3)])
>>> woe = woe.fit(lu[0], lu[1], drivers[[0, 10]])
>>> woe.table.columns[:5].tolist()
['Transition_From', 'Transition_To', 'Variable', 'Range_Lower_Limit', 'Range_Upper_Limit']

table property

table: DataFrame

Ranges and weights, as Dinamica EGO's report: transition, variable, range limits, possible and executed transitions, weight (W+), contrast and its significance (1 = 95 %).

ranges property

ranges: DataFrame

The range limits of each variable and transition.

fit

fit(t0, t1, drivers, *, nodata=None, n_threads: int | None = None) -> WeightsOfEvidence

Calibrate the ranges and weights from the map at the start (t0) and at the end (t1).

Parameters:

Name	Type	Description	Default
t0		Categorical maps of the two dates.	required
t1		Categorical maps of the two dates.	required
drivers		The explanatory variables on the same grid: an xarray.Dataset, a dict of maps, or a DataArray with a band per variable (the band names name them).	required
nodata		Value marking missing cells in the land use maps.	None
n_threads	int | None	Threads (default: all cores but one).	None

Returns:

Type	Description
WeightsOfEvidence	self, fitted.

correlation

correlation(*, n_threads: int | None = None) -> DataFrame

Association between the ranges of each pair of variables (Dinamica EGO's correlation report).

For each transition and pair of variables, the cross-tabulation of their ranges over every cell where both have a value gives chi-square, Cramér's V, the contingency coefficient, the joint entropy (natural logarithms) and the joint information uncertainty 2 (H(A) + H(B) - H(A, B)) / (H(A) + H(B)) (Bonham-Carter 1994). Values of Cramér's V or the joint uncertainty above ~0.5 suggest dropping or merging one of the two variables.

Returns:

Type	Description
DataFrame

probability

probability(
    landscape,
    drivers,
    *,
    transitions=None,
    dtype="float64",
    nodata=None,
    n_threads: int | None = None,
)

Transition probabilities from the weights (Dinamica EGO's CalcWOfEProbabilityMap).

With S the sum of the weights of the cell's range in every variable of a transition, the probability of the transition i -> j at a cell of class i is exp(S_ij) / (1 + sum_k exp(S_ik)) over the transitions from i asked for together (1 / (1 + exp(-S)) for a single one); a variable missing at a cell adds nothing to its sum (as in Dinamica).

Parameters:

Name	Type	Description	Default
landscape		The land use map the probabilities are for.	required
drivers		The explanatory variables (as in fit(); they may have changed since).	required
transitions		The transitions (default: all calibrated).	None
dtype	('float64', 'float32', 'uint8', 'uint16')	float: probabilities 0-1, NaN elsewhere; integer: trunc(p * max), 0 elsewhere (Dinamica's default uint8).	"float64"
nodata		Value marking missing cells in the landscape.	None
n_threads	int | None	Threads (default: all cores but one).	None

Returns:

Type	Description
DataArray	(transition, y, x), the transitions named "probability_<from>_to_<to>".

weights_of_evidence

weights_of_evidence(
    t0, t1, drivers, ranges: dict | None = None, **kwargs
) -> WeightsOfEvidence

Calibrate weights of evidence in one call: WeightsOfEvidence(ranges, **kwargs).fit(t0, t1, drivers).

See WeightsOfEvidence for the settings.

LogisticRegression

LogisticRegression(
    formulas: dict | None = None, *, maxit: int = 25, epsilon: float = 1e-08
)

Logistic regression of the presence of each land use class on the drivers (lulcc's glmModels).

Parameters:

Name	Type	Description	Default
formulas	dict | None	The drivers of each class, {class: [driver, ...]} (default: every class of the map on every driver), as the right-hand sides of lulcc's formulas class ~ d1 + d2; the intercept is always included.	None
maxit	int	glm.fit's iteration limit and convergence tolerance (R's defaults).	25
epsilon	int	glm.fit's iteration limit and convergence tolerance (R's defaults).	25

Examples:

>>> import landschaft as ls
>>> lu = ls.load_raster(ls.simulate.example_data("pie_landuse"))
>>> ef = ls.load_raster(ls.simulate.example_data("pie_drivers"))
>>> lr = ls.simulate.LogisticRegression({2: ["elevation", "slope", "dist_built_1985"]}).fit(lu[0], ef)
>>> lr.coefficients.columns.tolist()
['(Intercept)', 'elevation', 'slope', 'dist_built_1985']

coefficients property

coefficients: DataFrame

The coefficients of each class's model (NaN for aliased drivers, R's NA).

summary property

summary: DataFrame

Per class: cells used, iterations, convergence, deviance, null deviance.

fit

fit(
    landscape, drivers, *, cells=None, nodata=None, n_threads: int | None = None
) -> LogisticRegression

Fit a model per class: presence (1) or absence (0) of the class in the cells.

Parameters:

Name	Type	Description	Default
landscape		The categorical map (lulcc fits on the map of the first date).	required
drivers		The explanatory variables: a Dataset, a dict of maps or a DataArray with a band per variable.	required
cells		The cells to fit on: a boolean map or flat indices (row by row); default every cell. Cells missing in the landscape or in a driver of the model are left out (glm's na.omit).	None
nodata		Value marking missing cells in the landscape.	None

Returns:

Type	Description
LogisticRegression	self, fitted.

from_coefficients classmethod

from_coefficients(coefficients) -> LogisticRegression

A model with given coefficients (CLUMondo's and CLUE-S's alloc1.reg).

Parameters:

Name	Type	Description	Default
coefficients		{class: {"(Intercept)": b0, driver: b, ...}} or a DataFrame with a row per class and a column per term (NaN: not in the model), as coefficients gives it.	required

Examples:

>>> import landschaft as ls
>>> lr = ls.simulate.LogisticRegression.from_coefficients({1: {"(Intercept)": -1.5, "slope": 0.2}})
>>> lr.coefficients.loc[1, "slope"]
np.float64(0.2)

predict

predict(
    drivers, *, classes=None, landscape=None, nodata=None, n_threads: int | None = None
) -> DataArray

Probability of each class in every cell (predict(type = "response") of each model).

Parameters:

Name	Type	Description	Default
drivers		The explanatory variables (they may have changed since the fit).	required
classes		The classes to predict (default: all fitted).	None
landscape		Optional map whose missing cells stay NaN.	None

Returns:

Type	Description
DataArray	(class, y, x), NaN where a driver of the model is missing.

Estimator

Estimator(
    estimator,
    *,
    target: str = "presence",
    classes=None,
    sample=None,
    stratify: bool = False,
    balance: bool = False,
    seed=None,
)

Potential of each class from a scikit-learn classifier, one model per class.

Parameters:

Name	Type	Description	Default
estimator		A scikit-learn classifier (or any object with fit and predict_proba); a clone is fitted for each class.	required
target	('presence', 'expansion')	What the models learn: "presence", where each class is in the map (as lulcc's models); "expansion", where it appeared between two maps among the cells that did not have it (PLUS's land expansion analysis, LEAS), which needs the later map in fit().	"presence"
classes		The classes to model (default: all).	None
sample		Cells to train each model on: None (all), an int (at most this many, drawn at random) or a fraction in (0, 1).	None
stratify	bool	Draw sample from the presences and the absences separately, in their proportions (at least one of each): rare changes keep their cells.	False
balance	bool	Draw as many absences as presences (after sample).	False
seed		Seed of the draws (numpy's generator).	None

Examples:

>>> import landschaft as ls
>>> from sklearn.ensemble import RandomForestClassifier
>>> lu = ls.load_raster(ls.simulate.example_data("pie_landuse"))
>>> ef = ls.load_raster(ls.simulate.example_data("pie_drivers"))
>>> est = ls.simulate.Estimator(RandomForestClassifier(n_estimators=20, random_state=0),
...                             target="expansion", classes=[2], sample=5000, seed=1)
>>> p = est.fit(lu[0], ef, t1=lu[2]).predict(ef)
>>> p.dims
('class', 'y', 'x')

fit

fit(
    landscape,
    drivers,
    *,
    t1=None,
    cells=None,
    nodata=None,
    n_threads: int | None = None,
) -> Estimator

Fit a model per class.

Parameters:

Name	Type	Description	Default
landscape		The map at the start.	required
drivers		Explanatory variables (Dataset, dict or DataArray with a band per variable).	required
t1		The later map (target="expansion").	None
cells		Optional boolean map of the cells to train on.	None
nodata		Value marking missing cells in the maps.	None

predict

predict(
    drivers, *, classes=None, landscape=None, nodata=None, n_threads: int | None = None
) -> DataArray

Probability of each class (predict_proba of the class's model), NaN where a driver is missing.

importance

importance(method: str = 'impurity', *, n_repeats: int = 5, seed=None) -> DataFrame

Contribution of each driver to each class's model (PLUS's driver contributions).

Parameters:

Name	Type	Description	Default
method	('impurity', 'permutation')	"impurity": the estimator's feature_importances_ (a random forest's mean decrease in impurity); "permutation": the drop in accuracy on the training cells when a driver's values are shuffled (scikit-learn's permutation_importance, any estimator).	"impurity"
n_repeats	int	Shuffles per driver and their seed ("permutation").	5
seed	int	Shuffles per driver and their seed ("permutation").	5

Returns:

Type	Description
DataFrame	A row per driver, a column per class.

MultiCriteria

MultiCriteria(comparisons, *, drivers=None, transition=None)

Potential from weighted factors, the weights from pairwise comparisons (AHP; MOLUSCE's MCE).

Parameters:

Name	Type	Description	Default
comparisons		The pairwise comparison matrix: a square DataFrame indexed and named by the drivers (a_ij = how much more important driver i is than j, 1 to 9, a_ji = 1 / a_ij, a_ii = 1), or an array with drivers.	required
drivers		The driver names, in the order of an array's rows.	None
transition		(from, to): the potential is given only to the cells of the initial class (default: every cell).	None

Examples:

>>> import pandas as pd
>>> import landschaft as ls
>>> a = pd.DataFrame([[1, 3, 5], [1 / 3, 1, 2], [1 / 5, 1 / 2, 1]],
...                  index=["roads", "slope", "soil"], columns=["roads", "slope", "soil"])
>>> mce = ls.simulate.MultiCriteria(a)
>>> round(float(mce.weights["roads"]), 3), mce.consistency < 0.1
(0.648, True)

weights property

weights: Series

The weight of each driver (sum 1).

consistency property

consistency: float

Saaty's consistency ratio (acceptable below 0.1; 0 for two drivers, -1 beyond the table).

fit

fit(*args, **kwargs) -> MultiCriteria

Nothing to fit: the weights come from the comparisons.

predict

predict(
    drivers, *, landscape=None, nodata=None, n_threads: int | None = None
) -> DataArray

The potential: sum of weight x driver rescaled to 0-1 (by its minimum and maximum).

Parameters:

Name	Type	Description	Default
drivers		The factors (Dataset, dict or DataArray with a band per factor).	required
landscape		With transition: only the cells of its initial class get a potential (the others 0, missing cells NaN).	None

Returns:

Type	Description
DataArray	(transition, y, x) with one layer.

RasterPotential

RasterPotential(maps)

A potential given as maps: one layer per class (or transition), as computed elsewhere.

Parameters:

Name	Type	Description	Default
maps		A DataArray with a first dimension naming the classes (class, band or transition values), or a dict {class: map}.	required

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> pot = ls.simulate.RasterPotential({1: np.full((3, 4), 0.2), 2: np.full((3, 4), 0.7)})
>>> pot.predict().sel({"class": 2}).values.max()
np.float64(0.7)

fit

fit(*args, **kwargs) -> RasterPotential

Nothing to fit: the potential is given.

predict

predict(*args, classes=None, **kwargs) -> DataArray

The given maps (classes selects some).

Allocation

Model

Model(
    potential,
    demand,
    allocation,
    rules: Rules | None = None,
    dynamic: dict | None = None,
)

A land change model: potential, demand, allocation and rules.

Parameters:

Name	Type	Description	Default
potential		The suitability of each class: a fitted LogisticRegression, Estimator, RasterPotential... (anything with predict(drivers, landscape=...) giving a layer per class).	required
demand		The cells of each class at each step (AreaDemand).	required
allocation		CLUE or Ordered.	required
rules	Rules | None	Rules (allowed conversions, elasticities, restrictions, neighbourhood, history).	None
dynamic	dict | None	{name: piece}: drivers computed again from the map at the start of every step (distance_to(), share_of() or any function of the map), replacing the drivers of the same name; the potential is predicted again at every step.	None

Examples:

>>> import landschaft as ls
>>> lu = ls.load_raster(ls.simulate.example_data("pie_landuse"))
>>> ef = ls.load_raster(ls.simulate.example_data("pie_drivers"))
>>> forms = {1: ["elevation", "slope"], 2: ["elevation", "slope", "dist_built_1985"], 3: ["elevation", "slope"]}
>>> model = ls.simulate.Model(
...     potential=ls.simulate.LogisticRegression(forms).fit(lu[0], ef),
...     demand=ls.simulate.AreaDemand.from_maps(lu, times=[0, 6, 14]),
...     allocation=ls.simulate.CLUE(jitter=0.0002, scale=1e-6, max_diff=50, mean_diff=50),
...     rules=ls.simulate.Rules(elasticity=[0.2, 0.2, 0.2]))
>>> sim = model.run(lu[0], ef, seed=1)
>>> sim.maps.shape[0]
15

run

run(
    landscape,
    drivers=None,
    *,
    seed: int | None = None,
    nodata=None,
    n_threads: int | None = None,
) -> Simulation

Simulate every step of the demand from landscape.

Parameters:

Name	Type	Description	Default
landscape		The map at the start (step 0 of the demand).	required
drivers		The explanatory variables the potential needs.	None
seed	int | None	Seed of R's generator (set.seed): the same seed gives the same simulation, and the same as lulcc's with set.seed(seed). Default: drawn at random (kept in the result).	None
nodata		Value marking missing cells in the landscape.	None

Returns:

Type	Description
Simulation

save

save(
    path,
    landscape,
    drivers=None,
    *,
    seed=None,
    nodata=None,
    name=None,
    description=None,
)

Write the fitted model as a Terraflux <https://github.com/sacridini/terraflux>_ graph (JSON).

The fitted pieces are written as they are (weights of evidence, regressions, Markov chains, demands, patch settings, rules), the data the model runs on as GeoTIFF files in a folder next to it (<name>_data), so that terraflux run gives the same maps as run() with the same seed. Terraflux runs the CLUE-S, ordered, Expander/Patcher, CARS, MOLUSCE, CLUMondo and FUTURES allocations with weights of evidence, logistic regressions or given maps as potential, area demands or Markov chains (or a Regional of them), lulcc's rules and the dynamic variables distance_to() (Euclidean or the accumulated cost over a transition) and share_of(). Other potentials (an Estimator, MultiCriteria...) are written as their predictions on the starting map, which needs a model without dynamic variables.

Parameters:

Name	Type	Description	Default
path		The model file (.json).	required
landscape		The map the simulation starts from (as for run()), or the path of a raster (relative to the model file).	required
drivers		The explanatory variables, or the path of a raster with a band per variable (relative to the model file).	None
seed		The seed (default: drawn at random and written in the model).	None
nodata		Value marking missing cells in the landscape.	None

Returns:

Type	Description
Path

Examples:

>>> import landschaft as ls
>>> S = ls.simulate
>>> lu = ls.load_raster(S.example_data("sibuyan_landuse"))
>>> dr = ls.load_raster(S.example_data("sibuyan_drivers"))
>>> model = S.presets.dinamica(lu[0], lu[1], dr, steps=14, woe=250)
>>> path = model.save(tmp_path / "dinamica.json", lu[0], dr, seed=1)
>>> again, inputs = S.Model.load(path)

load classmethod

load(path)

Read a Terraflux graph written by save(): (model, inputs), the inputs by name.

Run it as before with model.run(inputs["start"], inputs.get("drivers"), seed=inputs["seed"]); graphs that are not made of ls.simulate's pieces are run with terraflux().

Simulation dataclass

Simulation(
    maps: DataArray,
    steps: DataFrame,
    seed: int,
    history: ndarray | None = None,
    ages: DataArray | None = None,
    changed: ndarray | None = None,
    supply_table: object = None,
    demand_table: object = None,
    stack: object = None,
    roles: object = None,
)

The maps of a simulation and what happened at each step.

Attributes:

Name	Type	Description
maps	DataArray	(step, y, x): the starting map and one map per step.
steps	DataFrame	Per step: passes of the allocation, whether it met the demand, the largest and mean difference (cells), the cells of each class.
seed	int	The seed of R's generator the simulation used.

layers property

layers: dict

With a LandSystemStack: {layer: (step, y, x)} binary maps, as the input.

supply property

supply: DataFrame

The goods and services supplied at each step (CLUMondo).

shortfall property

shortfall: DataFrame

Demand minus supply at each step (positive: short).

cover

cover() -> DataFrame

Cells of each layer (a stack) or class at each step.

trajectory

trajectory(roles=None) -> DataFrame

Cells changing at each step, by the role of their transition (Roles).

outputs

outputs(table, *, initial=None) -> DataArray

Per step, the value in table of each cell's class (CLUMondo's output maps).

Parameters:

Name	Type	Description	Default
table		A DataFrame with a row per class and a column per output (carbon, jobs...), or a Series for one output.	required
initial		Optional maps ({output: map}) of the real values at the start: cells that never changed class keep them (CLUMondo's initoutputmap).	None

Returns:

Type	Description
DataArray	(output, step, y, x).

AreaDemand

AreaDemand(table, classes=None)

Cells each class must cover at each time step (CLUE-S's demand).

Parameters:

Name	Type	Description	Default
table		A DataFrame with a row per time step (the first is the starting map) and a column per class, in cells; or an array with classes.	required
classes		The classes of an array's columns.	None

Examples:

>>> import landschaft as ls
>>> d = ls.load_table(ls.simulate.example_data("sibuyan_demand"))
>>> demand = ls.simulate.AreaDemand(d[d.scenario == "demand1"].drop(columns="scenario").set_index("time")
...                                 .rename(columns={"forest": 1, "coconut": 2, "grassland": 3, "rice": 4,
...                                                  "other": 5}))
>>> demand.steps
14

from_maps classmethod

from_maps(maps, times=None, out_times=None, *, nodata=None) -> AreaDemand

The areas of observed maps interpolated and extrapolated linearly (lulcc's approxExtrapDemand).

Each class's area is interpolated (or extrapolated from the first or last two maps) at out_times, negative values set to 0, and each time's areas rounded to integers that add up to the number of cells (the largest fractions rounded up).

Parameters:

Name	Type	Description	Default
maps		Observed maps (a band per date, or a list).	required
times		Their times (as in changes()).	None
out_times		Times of the demand (default: every integer time from the first to the last map).	None

at

at(step: int, classes: list) -> ndarray

The demand of each class at a step (0 = the start).

Rules dataclass

Rules(
    allowed: object = None,
    elasticity: object = None,
    restrictions: object = None,
    history: object = None,
    neighbourhood: Neighbourhood | None = None,
    allow_maps: dict | None = None,
    roles: object = None,
    transition_resistance: dict | None = None,
    conversion_order: list | None = None,
    caps: list | None = None,
    succession: list | None = None,
    _matrix: ndarray | None = None,
)

Decision rules of the allocation (lulcc's rules, elas, mask, hist, neighb).

Parameters:

Name	Type	Description	Default
allowed	object	Conversion matrix, rows the current class, columns the next (a DataFrame named by class, or an array in the model's class order): 1 allowed, 0 not; -1 the cell may leave its class only while the demand for that class decreases; -2 the next class may be taken only while its demand increases; 100 + t: allowed only while the cell has had its class for less than t steps; 1000 + t: allowed only once it has had it for t steps or more.	None
elasticity	object	Added to the suitability of each cell's own class (CLUE-S's conversion elasticity, 0 to 1; one value for all, {class: value} or a list in the class order). Used by the CLUE and CLUMondo allocations (not by the ordered one, as in lulcc).	None
restrictions	object	Map of restricted cells (True or non-zero: the cell keeps its class), lulcc's mask with the opposite convention.	None
history	object	Map of the steps each cell has had its class (default 1 everywhere, as lulcc), updated at each step: +1 where the class stays, 0 where it changes.	None
neighbourhood	Neighbourhood | None	Neighbourhood thresholds (lulcc's allocations).	None
allow_maps	dict | None	CLUMondo: {(from, to): map}, the conversion code of each cell for that pair (0 not allowed, 1 allowed, or an age code below).	None

Notes

With the CLUMondo allocation the conversion matrix takes CLUMondo's codes instead: 1 allowed, 0 not; 100 + t: only once the cell has had its system for t years; -100 - t: only until then (the system must leave after t years when the code is on the diagonal); 1000 + t: autonomous change (succession) after t years, whatever the demand. elasticity is CLUMondo's conversion resistance, history the starting ages (a map, or a number: random ages from 1 to it, as CLUMondo draws them).

With CLUMondo, rules may also be given by the role of the transitions (roles, a Roles):

• allowed as {role: code} ({"expansion": 0}; other pairs 1);
• transition_resistance: {role or (from, to): value} subtracted from the score of making that transition;
• conversion_order: roles from the first tried to the last, as a resistance of 0.1 per position (["intensification", "expansion"]: expansion 0.1); transition_resistance overrides;
• caps: Cap limits per year (Cap("expansion", area=...));
• succession: Succession (autonomous changes).

mask_values

mask_values(cells: ndarray, size: int) -> ndarray

lulcc's mask on the cells: 0 restricted, 1 free.

Neighbourhood dataclass

Neighbourhood(rules: dict, weights: object = None)

Neighbourhood thresholds (lulcc's NeighbRasterStack and nb.rules).

A class may be allocated to a cell only where the mean of weights x presence of the class over the window around the cell reaches its threshold (cells outside the map and missing cells are left out of the mean).

Parameters:

Name	Type	Description	Default
rules	dict	{class: threshold} (0-1 for a window of ones: the share of the window in the class).	required
weights	object	The window (odd sides), the same for every class, or {class: window} (default 3 x 3 ones).	None

allow

allow(landscape: ndarray, classes: list) -> ndarray

1 where each class may go, NaN where not: (cells of the map, classes).

CLUE dataclass

CLUE(
    jitter: float = 0.0001,
    scale: float = 0.0005,
    max_iter: int = 1000,
    max_diff: float = 5.0,
    mean_diff: float = 5.0,
)

CLUE-S allocation (Verburg et al. 2002), as lulcc's CluesModel.

Every cell takes the class of highest suitability plus an iteration variable per class; the variables move by scale x (demand - allocated cells) until the allocation meets the demand. Each pass jitters the suitabilities by U(-jitter, jitter).

Parameters:

Name	Type	Description	Default
jitter	float	Half-width of the random perturbation of the suitabilities.	0.0001
scale	float	Step of the iteration variables per cell of difference.	0.0005
max_iter	int	Passes before giving up (the allocation of the last pass is kept).	1000
max_diff	float	Largest and mean absolute difference between demand and allocation (cells) accepted. lulcc averages the signed differences, which add up to 0, so its mean criterion never acts; here it is the mean of the absolute differences, as CLUE-S defines it.	5.0
mean_diff	float	Largest and mean absolute difference between demand and allocation (cells) accepted. lulcc averages the signed differences, which add up to 0, so its mean criterion never acts; here it is the mean of the absolute differences, as CLUE-S defines it.	5.0

Ordered dataclass

Ordered(order: list, stochastic: bool = True)

The ordered allocation of Fuchs et al. (2013) with stochastic choice, as lulcc's OrderedModel.

The classes are served one at a time in order: a class whose demand grows takes the cells of other classes in decreasing suitability, a class whose demand shrinks releases its cells of lowest suitability (to the classes served after it). With stochastic, cells are drawn by comparing their suitability with uniform numbers (R's generator), in suitability order.

Parameters:

Name	Type	Description	Default
order	list	The classes, in the order they are served (all of them).	required
stochastic	bool	Draw the cells (default) or take the most suitable ones.	True

Dynamic variables and regions

distance_to

distance_to(classes, *, cost=None, unit: str = 'map')

Distance from each cell to the nearest cell of some classes.

Parameters:

Name	Type	Description	Default
classes		A class or a list of classes.	required
cost		None: the exact Euclidean distance (Felzenszwalb & Huttenlocher's transform), in map units (cells without a georeference); or a ls.cost Transition of conductances: the accumulated cost from the nearest such cell (gdistance's accCost), a distance along roads, rivers or any network.	None
unit	('map', 'cells')	Euclidean distance in the map's units (square cells) or in cells.	"map"

Notes

Calling it with a map gives the distance map (NaN where the map has no value). With no cell of the classes, every cell gets the length of the map's diagonal (Euclidean) or inf (cost).

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> land = np.array([[1, 1, 1, 3], [1, 1, 1, 1]])
>>> ls.simulate.distance_to(3)(land)
array([[3.        , 2.        , 1.        , 0.        ],
       [3.16227766, 2.23606798, 1.41421356, 1.        ]])

share_of

share_of(classes, window=3)

Share of the cells around each cell that hold some classes.

Parameters:

Name	Type	Description	Default
classes		A class or a list of classes.	required
window		Side of a square window of ones (odd), or the window's weights.	3

Notes

The mean of weights x presence over the window, the cells without a value left out (the centre included), as Neighbourhood.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> land = np.array([[1, 2, 2], [1, 1, 2], [1, 1, 1]])
>>> float(ls.simulate.share_of(2)(land)[1, 1])
0.3333333333333333

gravity_of

gravity_of(classes, exponent: float = 2.0, *, mass=None, unit: str = 'cells')

Gravitational interaction potential of the cells of some classes.

The sum over the cells k of the classes of m_k / d^exponent (gravity(), Dinamica EGO's CalcInteractionPotentialMap): the pull of towns or cleared land that grows as they grow.

Parameters:

Name	Type	Description	Default
classes		A class or a list of classes.	required
exponent	float	The distance decay exponent.	2
mass	DataArray or array_like	The mass of each cell (population, production...) where it holds one of the classes; default 1.	None
unit	('cells', 'map')	Distances in cells (Dinamica) or in map units.	"cells"

Notes

Calling it with a map gives the potential map (NaN where the map has no value; the cells without a value are no source).

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> land = np.array([[1, 1, 1, 3], [1, 1, 1, 1]])
>>> ls.simulate.gravity_of(3)(land)
array([[0.11111111, 0.25      , 1.        , 1.        ],
       [0.1       , 0.2       , 0.5       , 1.        ]])

Regional

Regional(regions, demands: dict)

A demand per region (municipalities, biomes, basins...).

Parameters:

Name	Type	Description	Default
regions		Map of region ids (cells without a value, or with an id not in demands, belong to no region and keep their class).	required
demands	dict	{region id: demand}: AreaDemand, Markov or a rate matrix / {(from, to): rate} as the allocation takes them, all with the same number of steps.	required

Notes

The allocations meet each region's demand on its cells: CLUE and Ordered allocate region by region (in the order of the ids, with the one generator), CARS with a demand and an adaptive coefficient per region and class, ExpanderPatcher with each region's rates and its probabilities set to 0 outside the region (Dinamica's regionalization).

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> regions = np.repeat([[1, 1, 2, 2]], 2, axis=0)
>>> t0 = np.array([[1, 1, 1, 2], [1, 2, 1, 1]])
>>> t1 = np.array([[1, 2, 1, 2], [2, 2, 2, 2]])
>>> d = ls.simulate.Regional.from_maps(regions, [t0, t1])
>>> d.at(1, [1, 2], region=1), d.at(1, [1, 2], region=2)
(array([1., 3.]), array([1., 3.]))

mask

mask(region) -> ndarray

Boolean map of a region's cells.

index

index(shape) -> ndarray

Region index (0 .. n - 1, in the order of the ids) of each cell, -1 outside the regions.

at

at(step: int, classes: list, region=None) -> ndarray

A region's area demand at a step (region=None: all regions added up).

from_maps classmethod

from_maps(regions, maps, times=None, out_times=None, *, nodata=None) -> Regional

Each region's areas in observed maps interpolated linearly (from_maps()).

markov classmethod

markov(regions, t0, t1, *, steps: int = 1, nodata=None) -> Regional

A Markov chain per region from two maps (from_maps()).

Land systems (CLUMondo)

CLUMondo dataclass

CLUMondo(
    supply: object,
    order: object = None,
    weights: object = None,
    deviation: str = "percent",
    mean_deviation: float = 0.5,
    max_deviation: float = 1.0,
    speed: float = 0.05,
    max_loops: int = 20000,
    kinds: dict | None = None,
    neighbourhood: NeighbourhoodEffect | None = None,
    preferences: dict | None = None,
    patches: dict | None = None,
    _cache: dict = dict(),
)

CLUMondo's demand-driven allocation of land systems.

Parameters:

Name	Type	Description	Default
supply	object	Supply (or its table): what each land system supplies per cell; its systems are the model's.	required
order	object	CLUMondo's lusconv: per service, the rank of each system (higher supplies more; -1: the system does not take part in that service's changes). Default: the rank of the systems' supply of the service.	None
weights	object	Weight of each service in the iteration variables ({service: w} or a list; default 1).	None
deviation	str	"percent": the deviations are percentages of the demand; "absolute": in the demand's units.	'percent'
mean_deviation	float	Largest mean deviation over the services and largest deviation of one service accepted. max_deviation may also be {service: percent}: each service within its own (the 2017 executable's iteration mode 2; mean_deviation is then not used).	0.5
max_deviation	float	Largest mean deviation over the services and largest deviation of one service accepted. max_deviation may also be {service: percent}: each service within its own (the 2017 executable's iteration mode 2; mean_deviation is then not used).	0.5
kinds	dict | None	{service: "exact" | "minimum" | "maximum"} (default exact): a minimum is met by any supply above it, a maximum by any below (the 2017 executable's demand types -1 and 101).	None
speed	float	Starting step of the iteration (0.001-0.9; CLUMondo's default 0.05).	0.05
max_loops	int	Iterations per year before giving up (the simulation stops there, as CLUMondo).	20000
neighbourhood	NeighbourhoodEffect | None	NeighbourhoodEffect.	None
preferences	dict | None	Location-specific additions to the suitability: {system: (map, weight)} (CLUMondo's locspec files; a map per step may be given as {step: map}).	None
patches	dict | None	{system: (radius, threshold, min_cells)}: the 2017 executable's patches. Each cell newly allocated to the system in an iteration brings along the cells of the square window of that radius that may convert to it and whose suitability reaches the threshold (whatever their resistance and the demand); a patch smaller than min_cells is dropped, seed included. Written from the manual (its source is not published).	None

Notes

The conversion resistance of each system is the model's Rules(elasticity=...); the conversion matrix uses CLUMondo's codes (see Rules); Rules(history=n) draws the starting ages as CLUMondo does (rand() % n + 1, the runtime's generator with the run's seed; CLUMondo's own runs use seed 1).

lusconv

lusconv() -> DataFrame

The conversion order (systems x services).

ServiceDemand

ServiceDemand(table, services=None)

The demand for goods and services at each time step (CLUMondo's demand.in).

Parameters:

Name	Type	Description	Default
table		A DataFrame with a row per time step (the first is the start) and a column per good or service, in their own units (tons, heads, hectares...); or an array with services.	required
services		The names of an array's columns.	None

Examples:

>>> import landschaft as ls
>>> d = ls.simulate.ServiceDemand({"crops": [100, 104, 108], "cattle": [50, 50, 55]})
>>> d.steps, d.services
(2, ['crops', 'cattle'])

Supply

Supply(table)

What one cell of each land system supplies of each good or service (CLUMondo's lusmatrix).

Parameters:

Name	Type	Description	Default
table		A DataFrame with a row per land system (its value in the map) and a column per good or service; or {step: table} for a supply that changes over time (CLUMondo's dynamic lusmatrix: each table holds from its step until the next one given; the first must be step 1).	required

Examples:

>>> import landschaft as ls
>>> s = ls.simulate.Supply({0: {"crops": 0.0}, 1: {"crops": 1.0}, 2: {"crops": 3.0}})
>>> s.systems, s.services
([0, 1, 2], ['crops'])

at

at(step: int) -> DataFrame

The table in force at a step (1 = the first simulated year).

of

of(landscape, *, step: int = 1, nodata=None) -> Series

The goods and services a map of land systems supplies.

NeighbourhoodEffect dataclass

NeighbourhoodEffect(systems: dict)

CLUMondo's neighbourhood function (neighmat.txt and alloc2.reg).

For each system i it applies to, the suitability of cells whose current system is in from_ gets weight x (constant + sum_k b_k x share_k), where share_k is the weighted share of system k in the kernel around the cell (over the cells with data), recomputed every year from the current map.

Parameters:

Name	Type	Description	Default
systems	dict	{system: {"weight": w, "kernel": radius or square array, "from": [current systems], "constant": c, "coefficients": {k: b}}}. The kernel is a radius (a square of ones) or the weights (odd side); default radius 1. from defaults to every system. As in CLUMondo, the share of a system is only measured for systems listed here with a positive weight (0 for the others).	required

LandSystemStack

LandSystemStack(
    layers: dict,
    *,
    overlap: dict | None = None,
    default_overlap: str = "mosaic",
    priority: list | None = None,
    min_cells: int = 0,
    allow_new_mosaics=False,
)

Binary layers, possibly overlapping, as land systems (bit masks).

Parameters:

Name	Type	Description	Default
layers	dict	{name: map}: 1 (or True) where the layer is present (at most 31 layers). Cells where every layer is missing (NaN) are no data.	required
overlap	dict | None	{(a, b): policy} for pairs of layers: "mosaic" (both stay: a system of its own), the name of the layer that wins, or "error".	None
default_overlap	str	The policy of the other pairs: "mosaic", "priority" (the layer earlier in priority wins) or "error".	'mosaic'
priority	list | None	The layers from the strongest (default: the order of layers).	None
min_cells	int	Combinations with fewer cells lose their weakest layers until they reach a combination with at least that many (or a single layer).	0
allow_new_mosaics		Systems absent from the map that the simulation may create: True (every pair of layers whose overlap is "mosaic") or a list of tuples of layer names.	False

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> forest = np.array([[1, 1, 0], [1, 0, 0]])
>>> pasture = np.array([[0, 1, 1], [0, 1, 1]])
>>> st = ls.simulate.LandSystemStack({"forest": forest, "pasture": pasture})
>>> st.overlaps()["cells"].tolist()
[2, 3, 1]
>>> st.names_of(3)
('forest', 'pasture')

map property

map: DataArray

The systems (bit masks) as a map, NaN where no data: the landscape for Model.

code_of

code_of(layers) -> int

The bit mask of a set of layer names.

names_of

names_of(code: int) -> tuple

The layers of a system (bit mask).

overlaps

overlaps() -> DataFrame

The systems of the map (after the overlap policies): layers and cells, by bit mask.

layers

layers(maps) -> dict

Binary layers of maps of systems ((..., y, x)): {name: same shape, 1/0, NaN no data}.

supply

supply(per_layer: dict, *, mosaic='split')

The supply of each system from the supply of its layers (a Supply).

Parameters:

Name	Type	Description	Default
per_layer	dict	{layer: {service: per cell}}.	required
mosaic		"split": a mosaic of k layers supplies 1/k of each layer's (each covers a share of the cell); "full": the sum of its layers'; {system or tuple of layers: {service: value}}: given (the others split).	'split'

suitability

suitability(per_layer: dict, *, combine: str = 'mean')

A potential for the systems from maps per layer (mosaics: the mean or max of their layers).

from_categorical classmethod

from_categorical(
    landscape, names: dict | None = None, *, nodata=None
) -> LandSystemStack

A stack without overlap from a categorical map (each class a layer; names: {class: name}).

Roles dataclass

Roles(natural: list, families: dict | None = None)

What each transition between land systems does to their layers.

Parameters:

Name	Type	Description	Default
natural	list	The natural layers (forest, native grassland...).	required
families	dict | None	{family: [layers from the least to the most intensive]}: the intensity ladders ({"cattle": ["extensive", "intensive"]}).	None

Notes

A transition is an expansion when it loses a natural layer, a regeneration when it gains one (without losing any), an intensification when, keeping the natural layers, it moves up a family's ladder or gains a productive layer, a deintensification when it moves down or loses one; stable for no change, other otherwise. set corrects single pairs.

Examples:

>>> import landschaft as ls
>>> r = ls.simulate.Roles(natural=["forest"], families={"cattle": ["extensive", "intensive"]})
>>> r.classify({"forest"}, {"extensive"}), r.classify({"extensive"}, {"intensive"})
('expansion', 'intensification')
>>> r.classify({"forest", "extensive"}, {"forest", "intensive"})
'intensification'

set

set(pair, role: str) -> Roles

Give a pair of systems ((from, to), codes or tuples of layers) its role.

table

table(systems, layers_of=None) -> DataFrame

The role of every pair of systems (layers_of: system -> its layers; default: itself).

matrix

matrix(systems, layers_of=None) -> ndarray

Roles as a (from, to) array of strings.

Cap dataclass

Cap(
    role: str | None = "expansion",
    cells: int | None = None,
    area: float | None = None,
    pairs: list | None = None,
    cell_area: float | None = None,
)

At most so many cells (or so much area) a year may make a role's transitions (or given pairs).

Parameters:

Name	Type	Description	Default
role	str | None	A role ("expansion": deforestation) or None with pairs.	'expansion'
cells	int | None	The limit per step, in cells, or in area units (the cells' area in the map's units, as cell_area: area=50_000, cell_area=1e4 / 1e4 for hectares).	None
area	int | None	The limit per step, in cells, or in area units (the cells' area in the map's units, as cell_area: area=50_000, cell_area=1e4 / 1e4 for hectares).	None
pairs	list | None	[(from, to), ...] systems instead of a role.	None
cell_area	float | None	Area of a cell for area (default: from the map's resolution).	None

Notes

When more cells would make them, those whose change scores highest over staying keep it; the others stay, and the demand iteration finds the supply elsewhere.

Succession dataclass

Succession(
    after_years: int, from_: object = None, to: object = None, role: str | None = None
)

Autonomous change after some years, whatever the demand (CLUMondo's code 1000 + years).

Parameters:

Name	Type	Description	Default
after_years	int	Years a cell must have had its system.	required
from_	object	The systems (values, or tuples of layer names with a stack); or	None
to	object	The systems (values, or tuples of layer names with a stack); or	None
role	str | None	"regeneration": each system goes to its regeneration target with the fewest layers (an error when that is ambiguous).	None

Expander and Patcher (Dinamica)

ExpanderPatcher dataclass

ExpanderPatcher(
    expansion: object = 0.5,
    expander: object = (1.0, 0.0, 1.0),
    patcher: object = (1.0, 0.0, 1.0),
    size_unit: str = "cells",
    prune_factor: float = 10.0,
    window: tuple = (3, 3),
    steps: int | None = None,
)

Dinamica's Expander and Patcher allocation of transitions (AllocateTransitions).

Each step, the changes of a transition are round(rate x cells of the origin class) (the demand's rates, e.g. a Markov chain); a share expansion goes to the Expander, the rest (and what the Expander could not place) to the Patcher.

• Expander: a cell's probability becomes P where at least 5 cells of its 3 x 3 window hold the target class, else P n_j / 5 (Dinamica EGO 8, measured; Soares-Filho et al. 2002 give n_j / 4 from 4). The cells are ranked, the first prune_factor x changes kept, drawn in proportion to that probability and each tried once: a seed with that probability.
• Patcher: seeds drawn in proportion to the probability among the first prune_factor x changes cells, until its changes are made.
• Each seed grows a patch of lognormal size (mean, variance) cell by cell, each new cell drawn among the origin-class cells around the patch in proportion to P x isometry^k, k the patch cells in its window: isometry 0 leaves single cells, 1 follows the probabilities, 2 makes compact patches.

Parameters:

Name	Type	Description	Default
expansion	object	Share of the changes made by the Expander (one value or {(from, to): share}).	0.5
expander	object	(mean, variance, isometry) of the patches (one tuple or per transition), the sizes in size_unit.	(1.0, 0.0, 1.0)
patcher	object	(mean, variance, isometry) of the patches (one tuple or per transition), the sizes in size_unit.	(1.0, 0.0, 1.0)
size_unit	str	"cells", "hectares" (Dinamica's: from the map's cell area in square metres) or "area" (the map's area units).	'cells'
prune_factor	float	How many cells per change are ranked for the draw (Dinamica's 10).	10.0
window	tuple	The neighbour window (rows, columns), odd.	(3, 3)
steps	int | None	Steps to simulate (default: the demand's steps, e.g. a Markov chain's).	None

Notes

The functors of Dinamica EGO 8.13 crash for patches of more than one cell, so only their seed selection could be measured (Expander: acceptance n_j / 5 per cell and the draw by probability, from 12 x 5000 trials; Patcher: the draw by probability among the pruned cells); the patch growth follows the published descriptions.

Stock taken by patches (Dinamica's continuous Patcher)

ContinuousPatcher dataclass

ContinuousPatcher(
    mean: float = 1.0,
    variance: float = 0.0,
    isometry: float = 1.0,
    unit: str = "hectares",
    prune_factor: float = 10.0,
    window: tuple = (3, 3),
)

New patches that take a stock (biomass, carbon, timber) up to a quota.

Dinamica EGO's PatcherContinuousState: the Patcher of a continuous state. Each patch seed is drawn among the cells with stock and probability, the patch grows to a lognormal size, every cell of it gives all its stock, and patches are made until the stock taken reaches the quota.

• Candidates: the cells with stock > 0 and probability > 0 ranked from the highest probability; the first ones kept while their summed stock is at most prune_factor x quota.
• Seeds: drawn from those without replacement in proportion to the probability (measured in Dinamica EGO 8.13). The list is not ranked again when it runs out: the quota is then left unmet.
• Patches: lognormal size of the given mean and variance; each new cell drawn among the cells with stock around the patch in proportion to P x isometry^k, k the patch cells in its window (as ExpanderPatcher); a patch stops as soon as the quota is met, so the last cell may exceed it.

Parameters:

Name	Type	Description	Default
mean	float	Mean and variance of the patch sizes, in unit.	1.0
variance	float	Mean and variance of the patch sizes, in unit.	1.0
isometry	float	0 leaves single cells, 1 follows the probabilities, 2 makes compact patches (0-2).	1.0
unit	('hectares', 'cells', 'area')	Units of the patch sizes and of the stock's area: Dinamica's stock per hectare and sizes in hectares (from the map's cell area in square metres), or per cell and in cells, or per map area unit and in map area units.	"hectares"
prune_factor	float	How much more stock than the quota the ranked cells hold (Dinamica's 10).	10.0
window	tuple	The neighbour window of patch growth (rows, columns), odd.	(3, 3)

Notes

Dinamica's functor crashes for patches of more than one cell (as its Patcher does), so the candidates, the pruning, the seed draw and the quota were measured with one-cell patches; the patch growth is the published Patcher's. Its random numbers depend on more than its seed (the same run in another folder differs), so the comparison is statistical.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> rng = np.random.default_rng(1)
>>> biomass = np.full((20, 20), 9.0)           # per cell
>>> p = rng.random((20, 20))
>>> cp = ls.simulate.ContinuousPatcher(mean=4, variance=1, isometry=1.5, unit="cells")
>>> a = cp.allocate(biomass, p, quota=100, seed=1)
>>> a.n_cells, a.remaining
(12, -8.0)

allocate

allocate(
    stock, probability, quota: float, *, seed: int | None = None, nodata=None
) -> StockAllocation

Take quota units of stock as new patches.

Parameters:

Name	Type	Description	Default
stock	DataArray or array_like	The stock of each cell (per hectare, per cell or per area unit: see unit; NaN, nodata or <= 0: none).	required
probability	DataArray or array_like	The probability of each cell (NaN or <= 0: never taken).	required
quota	float	The stock to take, in the stock's units times the area unit.	required
seed	int	Seed of the generator (xoshiro256**); default: drawn at random.	None

Returns:

Type	Description
StockAllocation

run

run(
    stock, probability, quotas, *, regrowth=None, seed: int | None = None, nodata=None
) -> StockSimulation

Take a quota per step, the probabilities depleted from step to step.

Parameters:

Name	Type	Description	Default
stock		As in allocate().	required
probability		As in allocate().	required
quotas	sequence of float	The quota of each step.	required
regrowth	callable	regrowth(stock, step) -> stock applied to the stock left after each step (numpy arrays, in the stock's units), e.g. growth of logged forest.	None
seed	int		None

Returns:

Type	Description
StockSimulation

StockAllocation dataclass

StockAllocation(
    stock: object,
    probability: object,
    taken: object,
    patches: object,
    remaining: float,
    n_patches: int,
    n_cells: int,
    seed: int,
)

What allocate() took.

Attributes:

Name	Type	Description
stock	DataArray or ndarray	The stock left (0 in the cells taken), in the stock's units.
probability	DataArray or ndarray	The probabilities with the cells taken at 0 (Dinamica's "depleted" probabilities, for the next step).
taken	DataArray or ndarray	The stock taken from each cell in the quota's units (the stock times the cell's area; 0 elsewhere).
patches	DataArray or ndarray	The number of the patch of each cell taken (1, 2, ...; 0 elsewhere).
remaining	float	The quota not met: negative when the last cell exceeded it, positive when the ranked cells ran out.
n_patches, n_cells	int	Patches and cells taken.
seed	int

StockSimulation dataclass

StockSimulation(stock: DataArray, taken: DataArray, steps: DataFrame, seed: int)

What run() did, step by step.

Attributes:

Name	Type	Description
stock	DataArray	The stock at the start and after each step (dimension step).
taken	DataArray	The stock taken in each step (quota units per cell).
steps	DataFrame	Per step: quota, taken, remaining, patches, cells.
seed	int

Roads built by the model (Dinamica's road constructor)

RoadBuilder dataclass

RoadBuilder(
    friction: object,
    attraction: object,
    quadrants: tuple = (2, 10),
    distance: tuple = (1500.0, 3000.0),
    spacing: float = 1000.0,
    per_quadrant: int = 1,
    max_attempts: int = 15,
    lottery: bool = True,
    diagonal: bool = True,
    relative: bool = False,
    mask: object = None,
)

New roads each step, from the network towards attractive places.

Dinamica's road constructor (ClassicRoadConstructor): each step, the attraction is kept where the distance to the roads is within distance (and mask allows); destinations are drawn in quadrants (destinations()); each is joined to the network down the cost of reaching it over the friction (friction_cost(), pathways(), ties at random with lottery), and the new roads are merged into the network.

Parameters:

Name	Type	Description	Default
friction	DataArray or array_like	The effort to build across each cell (> 0; NaN: never).	required
attraction	(DataArray, array_like or callable)	How attractive each cell is as a road's end, or a function of the simulated map giving it (e.g. a potential of deforestation).	required
quadrants	(int, int)	Horizontal and vertical divisions of the map (Dinamica's default).	(2, 10)
distance	(float, float)	Range of distances from the roads (map units) where new roads may end.	(1500, 3000)
spacing	float	Minimum distance between destinations of a step.	1000
per_quadrant	int	As in destinations().	1
max_attempts	int	As in destinations().	1
lottery	bool		True
diagonal	bool	As in friction_cost().	True
relative	bool	As in friction_cost().	True
mask	(DataArray, array_like or callable)	Where roads may end (true / non-zero), or a function of the simulated map giving it.	None

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> rng = np.random.default_rng(0)
>>> net = np.zeros((40, 40)); net[0, :] = 1
>>> rb = ls.simulate.RoadBuilder(np.ones((40, 40)), rng.random((40, 40)), quadrants=(2, 2),
...                              distance=(5, 20), spacing=5)
>>> new = rb.step(net, seed=1)
>>> int(np.nansum(new) > np.nansum(net))
1

step

step(network, *, landscape=None, seed: int | None = None)

One step of building: the network with the new roads.

landscape feeds attraction and mask when they are functions of the map.

run

run(network, steps: int, *, landscape=None, seed: int | None = None) -> DataArray

steps steps of building: the network at the start and after each step.

dynamic

dynamic(network, *, seed: int | None = None) -> RoadDistance

The distance to the roads as a dynamic variable of a Model that builds roads each step.

RoadDistance

RoadDistance(builder: RoadBuilder, network, seed: int)

Distance to a growing road network: a dynamic variable (see dynamic()).

Called with a map it gives the distance (map units) to the roads; called with a different map than the last time (the next step of a simulation), it first builds a step of roads, with the map feeding the builder's attraction and mask. networks keeps the network of each step; reset() goes back to the first.

friction_cost

friction_cost(
    sources,
    friction,
    *,
    diagonal: bool = True,
    relative: bool = False,
    directions: bool = False,
    nodata=None,
)

Cost of reaching every cell from the nearest source over a friction surface.

A step into a cell costs its friction, times sqrt(2) on a diagonal (diagonal), 8 neighbours; every cell gets the cost of its cheapest path from a source (Dinamica EGO's CalcCostMap; Dijkstra's algorithm, the same sums as Dinamica's push-broom run to convergence).

Parameters:

Name	Type	Description	Default
sources	DataArray or array_like	Cells with a value other than 0 are sources (cost 0).	required
friction	DataArray or array_like	The effort to cross each cell (> 0; NaN: impassable).	required
diagonal	bool	Diagonal steps cost sqrt(2) times more.	True
relative	bool	Frictions per map unit (metres): multiplied by the cell size (needs square cells).	False
directions	bool	Also give the direction of each cell's predecessor on its cheapest path (Dinamica's codes 1 NW, 2 N, 3 NE, 4 E, 5 SE, 6 S, 7 SW, 8 W; 0 at the sources).	False
nodata		Missing value of sources.	None

Returns:

Type	Description
xarray.DataArray or numpy.ndarray, or a tuple (cost, directions)	NaN where the friction is missing or no source is reached (directions: 255).

Notes

Ties between predecessors go to the first neighbour in row order (NW, N, NE, W, E, SW, S, SE); Dinamica's follow the order of its sweeps.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> roads = np.zeros((3, 4)); roads[0, 0] = 1
>>> ls.simulate.friction_cost(roads, np.ones((3, 4)), diagonal=False)
array([[0., 1., 2., 3.],
       [1., 1., 2., 3.],
       [2., 2., 2., 3.]])

destinations

destinations(
    attraction,
    quadrants=(1, 1),
    *,
    per_quadrant: int = 1,
    min_distance: float = 0.0,
    max_attempts: int = 15,
    seed: int | None = None,
    nodata=None,
)

Cells drawn by attraction in quadrants of a map.

The map is cut into quadrants = (horizontal, vertical) parts (the columns into horizontal, the rows into vertical); in each, up to per_quadrant cells are drawn in proportion to their attraction (cells of attraction <= 0 or missing never); a draw closer than min_distance (map units) to an earlier destination is drawn again, up to max_attempts times. Dinamica EGO's CalcDestinationMap.

Parameters:

Name	Type	Description	Default
attraction	DataArray or array_like		required
quadrants	(int, int)	Horizontal (column) and vertical (row) divisions.	(1, 1)
per_quadrant	int		1
min_distance	float	In map units (cells without a georeference).	0.0
max_attempts	int		15
seed	int		None

Returns:

Type	Description
DataArray or ndarray	The attraction of the destinations, NaN elsewhere.

Notes

Dinamica EGO 8.13 gives each quadrant n // q rows or columns (the last n % q never get a destination) and draws an integer of the summed attraction with rounding (the first cell of a quadrant gets half its share, and a draw past the end takes the quadrant's last cell even with no attraction, about 1 / (2 x summed attraction) of the time). Here the quadrants cover the map (floor(k n / q) boundaries) and the draw is in exact proportion to the attraction.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> a = np.zeros((4, 6)); a[1, 1] = 5; a[2, 4] = 1
>>> d = ls.simulate.destinations(a, (2, 1), seed=1)
>>> np.argwhere(~np.isnan(d)).tolist()
[[1, 1], [2, 4]]

pathways

pathways(
    network,
    sources,
    cost=None,
    *,
    directions=None,
    lottery: bool = False,
    merge: bool = True,
    value: float = 1,
    seed: int | None = None,
    nodata=None,
)

Paths from source cells to a network down a cost surface.

From each source (in row order) a path goes to the cheapest of its 8 neighbours until a cell of cost <= 0 (the network, e.g. the sources of friction_cost()); or, given directions, follows them until a cell of the network. Dinamica EGO's CalcPathwayMap.

Parameters:

Name	Type	Description	Default
network	DataArray or array_like	The network (cells with a value other than 0 or missing).	required
sources	DataArray or array_like	Where paths start (cells with a value other than 0), e.g. destinations().	required
cost	DataArray or array_like	The cost of reaching the network (<= 0 on it), e.g. friction_cost() from the network.	None
directions	DataArray or array_like	Instead of cost: the direction codes of friction_cost().	None
lottery	bool	Ties between cheapest neighbours drawn at random (else the first in row order: NW, N, NE, W, E, SW, S, SE).	False
merge	bool	The network with the new paths; else only the new paths.	True
value	float	The value of the new path cells.	1
seed	int	For the lottery.	None

Returns:

Type	Description
DataArray or ndarray	NaN off the network and the paths.

Notes

A source of missing cost is skipped; a local minimum of the costs (no cheaper neighbour before the network) is an error, as in Dinamica.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> net = np.zeros((3, 5)); net[:, 0] = 1
>>> src = np.zeros((3, 5)); src[1, 4] = 1
>>> cost = ls.simulate.friction_cost(net, np.ones((3, 5)))
>>> np.nan_to_num(ls.simulate.pathways(net, src, cost, value=2))
array([[1., 2., 2., 2., 0.],
       [1., 0., 0., 0., 2.],
       [1., 0., 0., 0., 0.]])

Multi-type random patch seeds (PLUS)

CARS dataclass

CARS(
    weights: object = 1.0,
    seeds: object = 0.1,
    decay: float = 0.9,
    step_size: float = 500.0,
    window: int = 3,
    max_iter: int = 5000,
    steps: int | None = None,
)

PLUS's CA based on multi-type random patch seeds (Liang et al. 2021).

Each step, CARS iterates from the current map until every class has the cells the demand asks for. In an iteration, each free cell i gives each class k the overall probability

• OP = P(i,k) x Omega(i,k) x D(k), Omega = w(k) x the cells of k in the window x window neighbourhood (the centre left out) over window^2 - 1;
• where Omega = 0, a random seed: OP = P x (r x mu(k)) x D when a uniform r < P, else 0;

a roulette wheel over OP picks a class c, and the cell turns to c when P(i,c) > tau and the conversion is allowed; tau = decay^l x r, r normal of mean 1 in [0, 2], l raised by one after each iteration that cut the total gap to the demand by less than step_size cells. D(k) starts at 1 and grows while the class's shortage grows (shrinks while its excess grows).

The growth probabilities P come from the model's potential (one layer per class, as an Estimator with target="expansion" gives: PLUS's LEAS); the demand is an AreaDemand or a Markov chain projected from the starting map (PLUS's default); the rules give the conversions allowed (allowed, 0/1, rows the current class) and the restricted cells (restrictions).

Parameters:

Name	Type	Description	Default
weights	object	Neighbourhood weight w(k) of each class, 0 to 1 (a number for all or {class: weight}; default 1).	1.0
seeds	object	Patch generation threshold mu(k) of each class, 0 to 1 (the paper's 0.1): how strongly random seeds compete with the growth of existing patches.	0.1
decay	float	Decay factor delta of the threshold, in (0, 1] (the paper's 0.9).	0.9
step_size	float	Gap reduction (cells per iteration) below which the threshold decays (PLUS's "Step"; the paper's 500).	500.0
window	int	Side of the neighbourhood, odd (3 = the Moore neighbourhood).	3
max_iter	int	Iterations per step before giving up (the map is kept as it is).	5000
steps	int | None	Steps to simulate (default: the demand's).	None

Notes

From the paper alone (PLUS has no source code): the standard deviation of the threshold's normal r is not given (here 1/3, values outside [0, 2] drawn again); the changes proposed in an iteration are made in a random order, a class taking no cell once it has its demand and losing none while it has no more than its demand (the paper does not say how a class is held at its demand). The random numbers of a cell come from a hash of (seed, step, iteration, cell): the same seed gives the same simulation with any number of threads.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> S = ls.simulate
>>> land = np.ones((40, 40)); land[15:25, 15:25] = 2
>>> p = S.RasterPotential({1: np.full((40, 40), 0.1), 2: np.full((40, 40), 0.8)})
>>> demand = S.AreaDemand(np.array([[1500, 100], [1400, 200]]), classes=[1, 2])
>>> sim = S.Model(p, demand, S.CARS()).run(land, seed=1)
>>> int((sim.maps[-1] == 2).sum())
200

MOLUSCE

MolusceWoE

MolusceWoE(bins: dict | None = None, unit_cell: float = 1)

MOLUSCE's weights of evidence (WoeManager): one model per transition, persistence included.

Each factor is cut into categories (bins: (-inf, b0) -> 1, [b0, b1) -> 2... [bn, inf) -> n + 2; or its integer values themselves); for each transition code observed between the two maps (the cells that went from i to j, i == j included), every category gets W+ + sum(W-) - W- over the factor's categories (halved with two categories), and a cell's total weight is the sum over the factors. A cell of class i is predicted to take the transition i -> j of largest total weight; its confidence is 100 x (sigmoid(largest) - sigmoid(second largest)), truncated.

Parameters:

Name	Type	Description	Default
bins	dict | None	{driver: [bounds]} (sorted); drivers not named are used as integer categories.	None
unit_cell	float	MOLUSCE's unit cell (areas in unit cells).	1

Notes

MOLUSCE keeps the largest weight seen before the largest as the second, which depends on the order of the classes: here the second largest weight (tested against MOLUSCE with tools/molusce_fixed.patch). MOLUSCE also stores the factors as 8-bit integers, wrapping integer categories outside 0-255 and truncating fractions; here integer categories keep their values and fractional values need bins.

fit

fit(t0, t1, drivers, *, nodata=None) -> MolusceWoE

Fit the weights on the change between two maps.

Parameters:

Name	Type	Description	Default
t0		The maps before and after (the same classes).	required
t1		The maps before and after (the same classes).	required
drivers		The factors (a Dataset, dict or DataArray with a band per variable), in the order MOLUSCE gets them.	required

weights_table

weights_table() -> DataFrame

The weight of each category of each driver, per transition.

total

total(drivers=None) -> dict

Total weight map of each transition code (NaN where a factor or the maps had none).

Without drivers, the maps of the fit (MOLUSCE's); with them, the weights applied to these drivers (dynamic variables).

potentials

potentials() -> dict

MOLUSCE's transition potentials: sigmoid(total weight) x 100 as 8-bit integers, per (from, to).

prediction

prediction(state, drivers=None, *, nodata=None) -> tuple[ndarray, ndarray]

The transition each cell is predicted to take and its confidence (0-100).

Returns:

Type	Description
(ndarray, ndarray)	The predicted class (NaN where the cell has none) and the confidence (NaN outside the cells with a class of the model).

predict

predict(drivers=None, *, landscape=None, nodata=None)

Total weight per transition as probability_i_to_j layers (sigmoid), for the other allocations.

MolusceCA dataclass

MolusceCA(iterations: int = 1)

MOLUSCE's simulator: each iteration makes the observed transitions where the model is most confident.

At each iteration the potential predicts every cell's transition and its confidence (MolusceWoE, or any potential with a layer per class: the class of largest potential, confidence 100 x the gap to the second); then, for each pair of classes in order, the n cells (the demand: the cross table of the observed change, the same every iteration) predicted to take that transition with the highest confidence plus U(0, 1) / 1000 change.

Parameters:

Name	Type	Description	Default
iterations	int	Iterations (steps) to simulate.	1

Notes

The random numbers are numpy's legacy generator (RandomState(seed)), as MOLUSCE draws them after numpy.random.seed(seed). Among cells of equal confidence and draw (which the 1/1000 rarely leaves), the later cells change first (MOLUSCE leaves the order to numpy's quicksort).

Urban growth (FUTURES)

PatchGrowing dataclass

PatchGrowing(
    patch_sizes: object,
    subregions: object = None,
    potential_subregions: object = None,
    pressure: object = None,
    neighborhood: int = 10,
    approach: str = "gravity",
    gamma: float = 1.5,
    scaling_factor: float = 1.0,
    discount_factor: float = 1.0,
    compactness_mean: float = 0.4,
    compactness_range: float = 0.05,
    num_neighbors: int = 4,
    seed_search: str = "probability",
    incentive_power: float = 1.0,
    weight: object = None,
    zoning: object = None,
    zoning_effects: object = None,
    climate: object = None,
)

FUTURES' patch growing allocation (r.futures.simulation).

Parameters:

Name	Type	Description	Default
patch_sizes	object	The library of patch sizes (cells): a list for every subregion, or {subregion: list} (a DataFrame with a column per subregion, NaN padding, as r.futures.calib writes it).	required
subregions	object	The subregions map (default: one subregion, id 1).	None
potential_subregions	object	Subregions of the potential coefficients, when they differ.	None
pressure	object	The starting development pressure (default: computed from the landscape with development_pressure(..., pad=True) and the parameters below).	None
neighborhood	int	The development pressure added around every new cell (the module's n_dev_neighbourhood, development_pressure_approach, gamma, scaling_factor).	10
approach	int	The development pressure added around every new cell (the module's n_dev_neighbourhood, development_pressure_approach, gamma, scaling_factor).	10
gamma	int	The development pressure added around every new cell (the module's n_dev_neighbourhood, development_pressure_approach, gamma, scaling_factor).	10
scaling_factor	int	The development pressure added around every new cell (the module's n_dev_neighbourhood, development_pressure_approach, gamma, scaling_factor).	10
discount_factor	float	Multiplies the patch sizes of the library.	1.0
compactness_mean	float	The distance exponent of the candidate ranking, drawn uniformly in mean +- range / 2 for every candidate.	0.4
compactness_range	float	The distance exponent of the candidate ranking, drawn uniformly in mean +- range / 2 for every candidate.	0.4
num_neighbors	int	4 or 8: the cells a patch grows to.	4
seed_search	str	"probability" (seeds drawn on the cumulative probability) or "random".	'probability'
incentive_power	float	Infill (> 1) or sprawl (< 1): the probability raised to this power (on a lookup of 1001 values, as the module).	1.0
weight	object	Map of -1 to 1 changing the probability (potential_weight).	None
zoning	object	Map of zone ids and their effects: None (the module's predefined zones) or a DataFrame indexed by subregion with a stringency column (0-2) and one column per zone id.	None
zoning_effects	object	Map of zone ids and their effects: None (the module's predefined zones) or a DataFrame indexed by subregion with a stringency column (0-2) and one column per zone id.	None
climate	object	Flooding: floods after each step's development, damage, retreat (abandonment), adaptation and migration between subregions.	None

DevelopmentDemand

DevelopmentDemand(table, subregions=None, population=None)

Cells each subregion develops at each step (r.futures.simulation's demand file).

Parameters:

Name	Type	Description	Default
table		A DataFrame with a row per simulated step (indexed by year) and a column per subregion (its value in the subregions map); or an array with subregions. Values are rounded to whole cells (lroundf).	required
subregions		The subregion of each column of an array.	None
population		The population each subregion adds per step, as table (used by the redistribution of Flooding).	None

Examples:

>>> import landschaft as ls
>>> d = ls.simulate.DevelopmentDemand({1: [10, 12], 2: [5, 0]})
>>> d.steps, d.subregions
(2, [1, 2])

from_population classmethod

from_population(
    developed,
    subregions,
    observed_population,
    projected_population,
    simulation_times,
    *,
    methods=("linear", "logarithmic"),
    nodata=None,
) -> DevelopmentDemand

The demand from the relation of developed cells to population (r.futures.demand).

For each subregion, the developed cells of the observed maps are regressed on the observed population ("linear", "logarithmic", "exponential", "logarithmic2", "exp_approach"; the smallest RMSE wins; a decreasing fit falls back to the latest cells per person), the population is interpolated linearly at simulation_times and the demand is the yearly difference of the predicted cells, truncated to whole cells.

Parameters:

Name	Type	Description	Default
developed		The observed maps of developed cells (1 developed), one per observed time (a list or a DataArray with a band per time).	required
subregions		The subregions map.	required
observed_population		DataFrames with the time as first column (or index name) and a column per subregion.	required
projected_population		DataFrames with the time as first column (or index name) and a column per subregion.	required
simulation_times		The times of the simulation, the first being the start.	required
methods		The candidate regressions.	('linear', 'logarithmic')

Returns:

Type	Description
DevelopmentDemand	With .population (the population added per step) and .curves (method, coefficients and RMSE per subregion).

DevelopmentPotential

DevelopmentPotential(table)

The development potential's coefficients per subregion (r.futures.potential's table).

p = 1 / (1 + exp(-(intercept + b_pressure * pressure + sum(b_k * x_k)))) with the coefficients of the cell's (potential) subregion.

Parameters:

Name	Type	Description	Default
table		A DataFrame indexed by subregion with the columns intercept, devpressure (or FUTURES' Intercept, development_pressure) and one per predictor (named as the drivers); or one row (a dict) for every subregion.	required

Examples:

>>> import landschaft as ls
>>> pot = ls.simulate.DevelopmentPotential({"intercept": -2.0, "devpressure": 0.5, "slope": -0.1})
>>> pot.predictors
['slope']

fit classmethod

fit(
    samples: DataFrame,
    *,
    response: str = "developed",
    pressure: str = "devpressure",
    predictors=None,
    subregions: str | None = "subregion",
    random: str | None = None,
    dredge: bool = False,
    min_variables: int = 1,
    max_variables: int | None = None,
    fixed=None,
) -> DevelopmentPotential

Fit the potential on sampled cells (r.futures.potential).

A binomial mixed model of new development: response ~ predictors + (1 | subregion) (R's glmer, Laplace), or with a random slope (random | subregion); with one subregion, R's glm (bit for bit). The coefficients per subregion are the fixed effects plus the subregion's random effects (lme4's coef).

Parameters:

Name	Type	Description	Default
samples	DataFrame	One row per sampled cell (rows with missing values are dropped).	required
response	str	The 0/1 column of new development.	'developed'
pressure	str	The development pressure column (the devpressure coefficient).	'devpressure'
predictors		The predictor columns (default: all other columns).	None
subregions	str | None	The subregion column (None: one subregion).	'subregion'
random	str | None	Column with a random slope per subregion (usually the pressure, which then has no fixed effect, as r.futures.potential).	None
dredge	bool	Choose the predictors by AIC over every subset (MuMIn's dredge) with min_variables to max_variables terms and the fixed ones always in.	False

Returns:

Type	Description
DevelopmentPotential	With attrs: aic, deviance (-2 log likelihood), and dredge (the ranked models) when used.

surface

surface(drivers, pressure, subregions=None) -> DataArray

The development probability (r.futures.potsurface): 1 / (1 + exp(-(intercept + b x ...))) in double.

Parameters:

Name	Type	Description	Default
drivers		The predictors (named as the table's columns).	required
pressure		The development pressure map.	required
subregions		The subregions map (cells of subregions without coefficients are NaN); default: the only row.	None

development_pressure

development_pressure(
    developed,
    size: int,
    *,
    method: str = "gravity",
    gamma: float | None = 1.5,
    scaling_factor: float | None = 1.0,
    pad: bool = True,
    nodata=None,
) -> DataArray

The development pressure of each cell (r.futures.devpressure).

The sum over the developed cells within size cells (circular window, the cell itself left out) of 1 ("occurrence"), scaling_factor / d^gamma ("gravity") or scaling_factor * exp(-2 d / gamma) ("kernel"), as GRASS's r.mfilter computes it.

Parameters:

Name	Type	Description	Default
developed		1 developed, 0 not, NaN excluded.	required
size	int	Radius of the window, in cells.	required
method	str	The decay of the pressure with distance.	'gravity'
gamma	str	The decay of the pressure with distance.	'gravity'
scaling_factor	str	The decay of the pressure with distance.	'gravity'
pad	bool	The module's -n flag (default): the map is padded with undeveloped cells so the cells near its edges get their pressure, and excluded cells count as undeveloped. Without it, as r.mfilter does, cells within size of the edges keep their input value and a window touching an excluded cell gives NaN.	True

Returns:

Type	Description
DataArray	float32, NaN where excluded.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> dev = np.zeros((7, 7))
>>> dev[3, 3] = 1
>>> p = ls.simulate.development_pressure(dev, 2, method="occurrence")
>>> int(p.values.sum())
12

Flooding dataclass

Flooding(
    catchments: object,
    adaptive_capacity: object,
    damage: object,
    response: tuple,
    response_stddev: float = 0.0,
    depths: dict | None = None,
    probabilities: dict | None = None,
    hand: object = None,
    hand_percentile: float | None = None,
    damage_regions: object = None,
    adaptation: object = None,
    redistribution: object = None,
)

FUTURES' flood submodel (Sanchez et al. 2023): floods, damage, response, abandonment, migration.

After the development of each step, a flood may strike each catchment: with probability max (the largest flood probability of the catchment's cells, or of the depth maps), at a probability p drawn uniformly. Its depth at each cell comes from the depth maps of the probabilities (the map of the largest probability below p) or, from flood probability maps, as the HAND percentile of the cells flooded at p minus the cell's HAND. The damage is the depth-damage curve's (%) unless the cell was adapted to such floods; a damaged developed cell then retreats (abandoned, its demand moved to a subregion drawn from the redistribution matrix, scaled by the population densities), adapts (protected up to the return period of the flood) or stays, from its adaptive capacity and the response lines, perturbed by a Gaussian.

Parameters:

Name	Type	Description	Default
catchments	object	Map of catchments (HUCs): floods are drawn per catchment.	required
adaptive_capacity	object	Map from -1 (most vulnerable) to 1 (most resilient).	required
damage	object	Depth-damage curves: a DataFrame with the depth levels as columns and the damage (%) per damage region in rows (one row: the same curve everywhere).	required
response	tuple	(vulnerability_a, vulnerability_b, resilience_a, resilience_b): the lines damage = a * capacity + b above which a cell retreats (vulnerable cells otherwise stay, resilient ones adapt).	required
response_stddev	float	Standard deviation of the Gaussian added to capacity and damage.	0.0
depths	dict | None	{step: {probability: map}}: flood depths for floods of each annual probability, from that step on (steps count from 1).	None
probabilities	dict | None	Instead of depths: {step: map} of flood probability, the HAND (height above nearest drainage) map and the percentile giving the flood level.	None
hand	dict | None	Instead of depths: {step: map} of flood probability, the HAND (height above nearest drainage) map and the percentile giving the flood level.	None
hand_percentile	dict | None	Instead of depths: {step: map} of flood probability, the HAND (height above nearest drainage) map and the percentile giving the flood level.	None
damage_regions	object	Which curve each cell uses: None (the only row), "subregions", "potential" (potential subregions) or a map of regions.	None
adaptation	object	Map of the return period each cell is already protected against.	None
redistribution	object	DataFrame with a row per subregion people leave and a column per destination: the probabilities (%). Without it, retreating cells are abandoned and their demand leaves the area.	None

patch_library

patch_library(
    start,
    end,
    *,
    subregions=None,
    per_subregion: bool = False,
    threshold: float = 0.0,
    nodata=None,
)

The sizes of the patches developed between two maps (r.futures.calib's patch library).

Patches are the 4-connected groups of cells developed in end and not in start (inside a subregion); their sizes, in cells, are what PatchGrowing draws from.

Parameters:

Name	Type	Description	Default
start		Developed (1) maps at two dates.	required
end		Developed (1) maps at two dates.	required
subregions		Map of subregions (cells outside, or in subregion 0, are ignored).	None
per_subregion	bool	A library per subregion ({subregion: sizes}), each patch cut to its part in the subregion.	False
threshold	float	Smallest patch area kept (in the map's area units).	0.0

Returns:

Type	Description
list or dict	The sizes, sorted.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> a = np.zeros((5, 6))
>>> b = a.copy()
>>> b[0, :3] = 1
>>> b[3:5, 4:6] = 1
>>> ls.simulate.patch_library(a, b)
[3, 4]

calibrate_patches

calibrate_patches(
    model,
    start,
    end,
    drivers=None,
    *,
    compactness_means=(0.4,),
    compactness_ranges=(0.05,),
    discount_factors=(1.0,),
    repeat: int = 10,
    seed: int = 1,
    threshold: float = 0.0,
    per_subregion: bool = False,
) -> DataFrame

Calibrate FUTURES' patch compactness and discount factor (r.futures.calib).

The patches observed between start and end give the patch library and the histograms of patch area and compactness (perimeter / (2 sqrt(pi area))); each combination of the parameters is simulated repeat times from start (the model's demand and potential) and its patches compared with the observed ones by the histograms' chi-square distance.

Parameters:

Name	Type	Description	Default
model		A Model with a PatchGrowing allocation (its patch library is replaced by the observed one).	required
start		The observed maps of development.	required
end		The observed maps of development.	required
drivers		The predictors of the potential.	None
compactness_means		The values to try (every combination).	(0.4,)
compactness_ranges		The values to try (every combination).	(0.4,)
discount_factors		The values to try (every combination).	(0.4,)
repeat	int	Simulations per combination (seeds (seed + k) * 10000 + i).	10
threshold	float	Smallest patch area kept.	0.0
per_subregion	bool	A patch library per subregion.	False

Returns:

Type	Description
DataFrame	Per combination: discount factor, compactness mean and range, the area and compactness errors (normalized by their largest) and their mean, from the best; attrs["library"] holds the patch library.

Validation

three_map_comparison

three_map_comparison(
    t0,
    t1,
    simulated,
    factors=(1,),
    *,
    categories=None,
    nodata=None,
    n_threads: int | None = None,
) -> ThreeMapComparison

Compare the reference maps at the start and the end with the simulated end at several resolutions.

Parameters:

Name	Type	Description	Default
t0		The reference maps at the start and at the end.	required
t1		The reference maps at the start and at the end.	required
simulated		The simulated map of the end.	required
factors		Block sides (cells) of the resolutions (1: the cells themselves).	(1,)
categories		The classes (default: those found in the three maps).	None

Returns:

Type	Description
ThreeMapComparison	With table, agreement_budget and figure_of_merit.

Notes

lulcc takes the cells with a value in t0 and, where t1 or the simulation lack one there, misaligns their blocks; here the cells without a value in any of the three maps are left out of all three.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> t0 = np.array([[1, 1, 1, 1], [1, 1, 2, 2]])
>>> t1 = np.array([[2, 2, 1, 1], [1, 1, 2, 2]])
>>> sim = np.array([[2, 1, 2, 1], [1, 1, 2, 2]])
>>> tmc = ls.simulate.three_map_comparison(t0, t1, sim, factors=[1, 2])
>>> tmc.figure_of_merit()["overall"].round(3).tolist()
[0.333, 0.333]

ThreeMapComparison

ThreeMapComparison(tables: dict, factors: list, categories: list)

Three maps compared at several resolutions (lulcc's ThreeMapComparison).

Attributes:

Name	Type	Description
factors	list of int	The resolutions (block sides, in cells).
categories	list	The classes.
tables	dict	{factor: array}: lulcc's table, (n + 1)^2 rows (reference t0 class i within reference t1 class j, row i + j (n + 1), 0-based; index n the totals) by n + 1 columns (simulated class, then the total), as shares of the area.

table

table(factor=1) -> DataFrame

The three-dimensional table as a DataFrame: rows (reference t0, reference t1), columns simulated.

agreement_budget

agreement_budget() -> DataFrame

Sources of agreement and disagreement (lulcc's AgreementBudget), as shares of the area.

misses (change simulated as persistence), hits, wrong_hits (change simulated as change to another class), false_alarms (persistence simulated as change) and correct_rejections (persistence simulated correctly; not for a transition): a row per level (overall, a class, a transition) and resolution.

figure_of_merit

figure_of_merit() -> dict

The figure of merit (Pontius et al. 2011, equations 9-11) at each resolution.

Returns:

Type	Description
dict	overall: a Series by factor; category: a DataFrame (factor x class: change from the class); transition: {factor: DataFrame} (rows: from, columns: to).

compare_maps

compare_maps(reference, simulated, *, nodata=None) -> dict

MOLUSCE's comparison of two maps (DependenceCoef): kappa, agreement and dependence.

Returns:

Type	Description
dict	kappa, kappa_location, kappa_histogram and correctness (% of cells of the same class) from the cross table over the classes of both maps (MOLUSCE's validation); cramer (Cramer's V) and jiu (joint information uncertainty) from the cross table of the classes each map has; correlation (Pearson, of the class values).

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> a = np.array([[1, 1, 2, 2], [1, 3, 3, 2]])
>>> b = np.array([[1, 1, 2, 3], [1, 3, 3, 2]])
>>> round(ls.simulate.compare_maps(a, b)["kappa"], 4)
0.814

error_budget

error_budget(
    reference, simulated, *, iterations: int = 3, scale: int = 2, nodata=None
) -> DataFrame

MOLUSCE's error budget (Pontius 2002): agreement for levels of information, at coarser resolutions.

The share of agreement between the maps when the simulation is given no / medium / perfect information of quantity and of location (no_no, no_medium, medium_medium, medium_perfect, perfect_perfect), at the map's resolution and iterations - 1 times coarser, each time merging scale x scale cells (the edges left over dropped, as MOLUSCE does).

Returns:

Type	Description
DataFrame	A row per resolution (cells: the block side in cells).

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> a = np.repeat(np.array([[1, 1, 2, 2], [1, 3, 3, 2]]), 2, axis=0)
>>> b = np.repeat(np.array([[1, 1, 2, 3], [1, 3, 3, 2]]), 2, axis=0)
>>> ls.simulate.error_budget(a, b, iterations=2)["medium_medium"].round(3).tolist()
[0.875, 0.875]

futures_validation

futures_validation(reference, simulated, original=None, *, nodata=None) -> dict

FUTURES' validation (r.futures.validation): disagreement, kappa, kappa simulation, change metrics.

Parameters:

Name	Type	Description	Default
reference		The observed and the simulated maps (integer classes; FUTURES uses 0 undeveloped, 1 developed).	required
simulated		The observed and the simulated maps (integer classes; FUTURES uses 0 undeveloped, 1 developed).	required
original		The map at the start: adds kappa simulation (van Vliet et al. 2011) and, for 0/1 maps, the change metrics of Pontius et al. (2008).	None

Returns:

Type	Description
dict	classes; per class quantity and allocation disagreement (Pontius & Millones 2011) and their totals total_quantity, total_allocation; kappa; with the original: kappasimulation and (0/1 maps) misses, hits, false_alarms, null_successes, initially_developed (shares of the cells), figure_of_merit, producer and user accuracy. The cells without a value in any map are left out, as r.stats -n does.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> orig = np.array([[0, 0, 1], [0, 0, 1]])
>>> ref = np.array([[1, 0, 1], [1, 0, 1]])
>>> sim = np.array([[1, 1, 1], [0, 0, 1]])
>>> round(ls.simulate.futures_validation(ref, sim, orig)["figure_of_merit"], 4)
0.3333

futures_validation_grid

futures_validation_grid(
    reference,
    simulated,
    original=None,
    *,
    size: int | None = None,
    regions=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataFrame

FUTURES' validation per tile or per region (r.futures.gridvalidation).

Parameters:

Name	Type	Description	Default
reference		As futures_validation().	required
simulated		As futures_validation().	required
original		As futures_validation().	required
size	int | None	Tiles of size x size cells from the top-left corner.	None
regions		Or a map of region ids (one row per region).	None

Returns:

Type	Description
DataFrame	A row per tile (row, col and the centre x, y in map units when the maps are georeferenced) or region (region), the metrics of futures_validation() (None where they cannot be computed), and cells.

fuzzy_similarity

fuzzy_similarity(
    first,
    second,
    *,
    window: int = 5,
    decay: str = "exponential",
    nodata=None,
    n_threads: int | None = None,
) -> FuzzySimilarity

Fuzzy similarity of two categorical maps (Dinamica EGO's CalcReciprocalSimilarityMap).

At each cell with a class in first: the membership of that class in second around the cell (2^(-d/2) for the nearest cell of the class within the window x window square, d in cells; or 1 with decay="constant"; 0 if none), the smaller of it and the reverse membership where second has a class too; and the same for the cells of second.

Parameters:

Name	Type	Description	Default
first		Categorical maps (cells without a value take no part).	required
second		Categorical maps (cells without a value take no part).	required
window	int	Side of the square window, odd (1: only the cell itself).	5
decay	('exponential', 'constant')		"exponential"

Returns:

Type	Description
FuzzySimilarity

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> a = np.full((5, 5), np.nan); a[2, 1] = 1
>>> b = np.full((5, 5), np.nan); b[2, 3] = 1
>>> s = ls.simulate.fuzzy_similarity(a, b, window=5)
>>> s.first_mean, s.second_mean
(0.5, 0.5)

FuzzySimilarity dataclass

FuzzySimilarity(
    first: object,
    second: object,
    first_mean: float,
    second_mean: float,
    window: int,
    decay: str,
)

The fuzzy similarity of two maps (fuzzy_similarity()).

Attributes:

Name	Type	Description
first, second		The similarity at each cell with a class in the first (second) map (NaN elsewhere).
first_mean, second_mean		Their means (Dinamica's "first" and "second" similarity).
window, decay		The settings.

minimum property

minimum: float

The smaller mean: Dinamica's two-way similarity.

differences

differences(t0, t1, *, nodata=None)

The map of differences: each cell that changed between the maps, in its new class; NaN elsewhere.

similarity_of_differences

similarity_of_differences(
    t0,
    observed,
    simulated,
    *,
    windows=(1, 3, 5, 7, 9, 11),
    decay="exponential",
    nodata=None,
    n_threads: int | None = None,
) -> DataFrame

Dinamica EGO's validation: fuzzy similarity of the observed and the simulated change, per window.

The maps of differences() of the observed and the simulated end are compared with fuzzy_similarity().

Returns:

Type	Description
DataFrame	A row per window: simulated (the mean similarity over the simulated changes), observed (over the observed changes), minimum (the smaller: Dinamica's CalcSimilarityOfDifferences result) and the cells of each.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> t0 = np.ones((6, 6))
>>> obs = t0.copy(); obs[2, 2] = 2
>>> sim = t0.copy(); sim[2, 3] = 2
>>> ls.simulate.similarity_of_differences(t0, obs, sim, windows=[1, 3])["minimum"].round(4).tolist()
[0.0, 0.7071]

validate

validate(
    t0,
    observed,
    simulated,
    *,
    factors=(1, 2, 4, 8, 16),
    windows=(1, 3, 5, 7, 9, 11),
    decay: str = "exponential",
    budget_iterations: int = 4,
    metrics=None,
    level: str = "class",
    nodata=None,
    n_threads: int | None = None,
) -> Validation

Compare a simulated map with the observed one in every way the land change tools do.

Parameters:

Name	Type	Description	Default
t0		The map at the start of the simulated period.	required
observed		The observed and the simulated end.	required
simulated		The observed and the simulated end.	required
factors		Resolutions of the three-map comparison (block sides, cells).	(1, 2, 4, 8, 16)
windows		Windows of the fuzzy similarity of the differences.	(1, 3, 5, 7, 9, 11)
decay		Windows of the fuzzy similarity of the differences.	(1, 3, 5, 7, 9, 11)
budget_iterations	int	Resolutions of MOLUSCE's error budget (each twice the last).	4
metrics		Landscape metrics (ls.landscape_metrics names, at level) to put side by side.	None

Returns:

Type	Description
Validation

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> rng = np.random.default_rng(0)
>>> t0 = rng.integers(1, 4, (40, 40)).astype(float)
>>> obs = t0.copy(); obs[5:15, 5:15] = 3
>>> sim = t0.copy(); sim[7:17, 6:16] = 3
>>> v = ls.simulate.validate(t0, obs, sim, factors=[1, 2, 4], windows=[1, 3, 5])
>>> bool(v.summary()["similarity_5"] > v.summary()["similarity_1"])
True

Validation dataclass

Validation(
    three_map: object,
    similarity: DataFrame,
    kappa: dict,
    error_budget: DataFrame,
    futures: dict,
    metrics: DataFrame | None = None,
    settings: dict = dict(),
)

Every comparison of a simulated map with the observed one (validate()).

Attributes:

Name	Type	Description
three_map	ThreeMapComparison	Pontius et al. (2011) at several resolutions (agreement_budget(), figure_of_merit()).
similarity	DataFrame	Fuzzy similarity of the differences per window (Dinamica EGO).
kappa	dict	MOLUSCE's kappa, kappa location and histogram, correctness...
error_budget	DataFrame	MOLUSCE's error budget at coarser resolutions.
futures	dict	Quantity and allocation disagreement, kappa simulation (FUTURES).
metrics	DataFrame or None	Landscape metrics of the observed and the simulated map side by side.

summary

summary() -> Series

The headline numbers: figure of merit, kappas, disagreement, similarity at each window.

Calibration

calibrate

calibrate(
    model,
    start,
    observed,
    drivers=None,
    *,
    params: dict,
    targets: dict | None = None,
    method: str = "nelder-mead",
    replicates: int = 1,
    seed: int | None = None,
    max_evals: int = 50,
    grid_points: int = 5,
    x0: dict | None = None,
    n_workers: int | None = None,
    nodata=None,
    on=None,
    n_threads: int | None = None,
) -> Calibration

Search the settings of a model that best reproduce an observed map.

Parameters:

Name	Type	Description	Default
model		A Model, its settings named by params as paths into its pieces ("allocation.step_size", "allocation.seeds[2]", "rules.elasticity[1]"); or a function build(**params) returning a model.	required
start		The map the simulation starts from and the observed map at its end.	required
observed		The map the simulation starts from and the observed map at its end.	required
drivers		The drivers of the model.	None
params	dict	{name: (low, high)} (integers if both bounds are) or {name: [values]} (searched by method="grid").	required
targets	dict | None	{name: weight}: "figure_of_merit" (1 - FoM), "similarity" or "similarity_<window>" (1 - fuzzy similarity of the differences, window 5 by default), "kappa", "quantity", "allocation" (disagreement); "lsm": ([metrics], weight) or {"metrics": [...], "level": "landscape", "weight": w} (mean relative difference of the landscape metrics); or {name: f} / {name: (f, weight)} with f(observed, simulated) a loss. Default: figure of merit and similarity, equally weighted.	None
method	('nelder-mead', 'differential-evolution', 'grid', 'random')	The search (scipy's optimizers; "grid": every combination, each range cut into grid_points values; "random": max_evals Latin hypercube draws).	"nelder-mead"
replicates	int	Simulations per evaluation (their mean loss), seeds seed + r in every evaluation.	1
max_evals	int	Evaluations at most (not for "grid").	50
x0	dict | None	Starting settings for Nelder-Mead (default: the middle of the ranges).	None
on		ls.helpers(): evaluate the points of a "grid" or "random" search on this computer and on the helpers of its group (see parallel), with the same losses. Each point's model travels as its Terraflux graph (save()), and the targets must be the built-in ones.	None

Returns:

Type	Description
Calibration

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> S = ls.simulate
>>> rng = np.random.default_rng(0)
>>> d1, d2 = rng.random((30, 30)), rng.random((30, 30))
>>> land = np.ones((30, 30))
>>> demand = S.AreaDemand(np.array([[900, 0], [800, 100]]), classes=[1, 2])
>>> def build(w):
...     pot = S.RasterPotential({1: np.full((30, 30), 0.5), 2: w * d1 + (1 - w) * d2})
...     return S.Model(pot, demand, S.Ordered(order=[2, 1], stochastic=False))
>>> observed = build(0.3).run(land, seed=1).maps[-1]
>>> cal = S.calibrate(build, land, observed, params={"w": [0.0, 0.1, 0.2, 0.3, 0.4, 0.5]},
...                   targets={"figure_of_merit": 1}, method="grid")
>>> cal.best
{'w': 0.3}

Calibration dataclass

Calibration(
    best: dict,
    loss: float,
    history: DataFrame,
    model: object,
    simulation: object,
    settings: dict = dict(),
    _maps: tuple = None,
)

The result of calibrate().

Attributes:

Name	Type	Description
best	dict	The settings of lowest loss.
loss	float	Its loss.
history	DataFrame	Every evaluation: the settings, the loss and each target's part.
model	Model	The model with the best settings.
simulation	Simulation	Its first replicate.

validate

validate(**settings)

Every comparison of the best simulation with the observed map (validate()).

Presets and one-call projection

run

run(
    t0,
    t1,
    drivers=None,
    *,
    until=None,
    years=None,
    preset="dinamica",
    seed=None,
    validate=True,
    dynamic=None,
    rules=None,
    nodata=None,
    **options,
) -> Projection

Fit a land change model on two maps, check it on the observed period and project it.

1. The preset is fitted on t0 -> t1 with one step per year;
2. validate: the model simulates the observed period from t0 and its last map is compared with t1 (agreement());
3. until: the model projects from t1, a step per year.

Parameters:

Name	Type	Description	Default
t0		The earlier and later maps.	required
t1		The earlier and later maps.	required
drivers		The explanatory variables.	None
until		Last year of the projection (default: none, only the calibration).	None
years		(year of t0, year of t1); default: the maps' year/time/ band coordinate, else (0, 1).	None
preset		"dinamica" (Markov + weights of evidence + Expander/Patcher with the observed patch sizes), "plus" (LEAS + CARS) or "molusce" (its WoE and simulator, a step per observed period), or a function f(t0, t1, drivers, steps=..., **options) returning a Model.	'dinamica'
seed		Seed of every random draw (default: drawn, kept in settings).	None
validate		Simulate the observed period first and compare with t1.	True
dynamic		Dynamic variables and rules given to the model.	None
rules		Dynamic variables and rules given to the model.	None
options		The preset's settings (e.g. woe=30, patches=0.6, isometry=1.3 for "dinamica"; estimator=, seeds= for "plus"; bins= for "molusce").	{}

Returns:

Type	Description
Projection

Examples:

>>> import landschaft as ls
>>> lu = ls.load_raster(ls.simulate.example_data("sibuyan_landuse"))
>>> drivers = ls.load_raster(ls.simulate.example_data("sibuyan_drivers"))[[0, 10]]
>>> proj = ls.simulate.run(lu[0], lu[1], drivers, until=2015, seed=1, woe=250)   # maps of 1997 and 2011
>>> proj.maps["year"].values.tolist()
[2011, 2012, 2013, 2014, 2015]
>>> sorted(proj.report["stage"])
['fit', 'projection', 'validation']

Projection dataclass

Projection(
    maps: DataArray,
    calibration: object,
    validation: dict | None,
    model: object,
    report: DataFrame,
    settings: dict,
    simulation: object = None,
    _call: dict = None,
)

The result of run().

Attributes:

Name	Type	Description
maps	DataArray	(year, y, x): the later observed map and the projected years.
calibration	Simulation or None	The simulation of the observed period, from the earlier map.
validation	dict or None	agreement() of the calibration's last map with the later map.
model	Model	The fitted model.
report	DataFrame	What each stage did.
settings	dict	The choices (preset, years, seed, stage settings), enough to repeat the projection.

redo

redo(**changes) -> Projection

Run again from the same maps with some settings changed (proj.redo(seed=2, until=2060)).

validate

validate(**settings)

Every comparison of the calibration's last map with the later map (validate()).

to_dataset

to_dataset() -> Dataset

The projected maps (and the calibration's, on its own year axis) as one Dataset.

agreement

agreement(observed, simulated, start, *, nodata=None) -> dict

How a simulated map matches the observed one, cell by cell and in quantity.

Returns overall (share of cells of the same class), figure_of_merit (hits / (hits + misses + false alarms + wrong hits), Pontius et al. 2008), the hits, misses, false_alarms and wrong_hits (change simulated where it happened, to another class), and the quantity difference of each class (simulated - observed cells).

presets

The classic land change models, assembled from ls.simulate's pieces.

Each function returns a Model. The first four take the pieces the model needs (they name the allocation and its usual settings); dinamica, plus and molusce fit everything from two maps and the drivers, as those programs do::

import landschaft as ls

S = ls.simulate
model = S.presets.dinamica(lc2010, lc2020, drivers, steps=10)   # Markov + WoE + Expander/Patcher
sim = model.run(lc2010, drivers, seed=1)

clue_s

clue_s(potential, demand, *, rules=None, dynamic=None, **clue)

CLUE-S (Verburg et al. 2002) as lulcc runs it: CLUE allocation of an AreaDemand.

clue are CLUE's settings (jitter, scale, max_iter, max_diff, mean_diff); rules the Rules (elasticities, conversion matrix, neighbourhood...).

ordered

ordered(potential, demand, order, *, rules=None, dynamic=None, stochastic: bool = True)

The ordered procedure of Fuchs et al. (2013), lulcc's OrderedModel.

clumondo

clumondo(potential, demand, supply, *, rules=None, dynamic=None, **settings)

CLUMondo: land systems meeting a ServiceDemand with their Supply.

settings are CLUMondo's (order, deviation, speed, neighbourhood...); give the starting map or a LandSystemStack to run.

futures

futures(potential, demand, patch_sizes, *, rules=None, **settings)

FUTURES' patch growing (r.futures.simulation): PatchGrowing with a DevelopmentDemand.

observed_patches

observed_patches(t0, t1, transitions=None, *, nodata=None) -> DataFrame

The patches of each observed transition: expansions of existing patches and new ones.

The cells that went from i to j are grouped in patches (8 neighbours); a patch touching a cell of j in t0 is an expansion (Dinamica's Expander), the others are new patches (the Patcher).

Returns:

Type	Description
DataFrame	Per transition: cells changed, expansion (share of the cells in expansions), the number, mean and variance (cells) of the expansion and new patches.

dinamica

dinamica(
    t0,
    t1,
    drivers,
    *,
    steps: int = 1,
    woe=None,
    patches="observed",
    isometry: float = 1.5,
    prune_factor: float = 10.0,
    rules=None,
    dynamic=None,
    nodata=None,
)

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

Fitted on two maps: the rates of each transition per step (the steps-th root of the period's matrix), weights of evidence of every observed transition (the ranges found automatically, or woe: the WeightsOfEvidence settings, a number = the increment), and the share of the Expander, the mean and variance of the patch sizes of each transition measured on the observed change (observed_patches(); or patches: the ExpanderPatcher settings, a number = the Expander's share); isometry is not measurable and given (Dinamica's usual 1.5).

plus

plus(
    t0,
    t1,
    drivers,
    *,
    steps: int = 1,
    estimator=None,
    sample=0.05,
    demand=None,
    rules=None,
    dynamic=None,
    seed=None,
    nodata=None,
    **cars,
)

PLUS (Liang et al. 2021): LEAS (a random forest per class on the expansion areas) and CARS.

estimator is any scikit-learn classifier (default: a random forest of 50 trees, the paper's), fitted for each class that expanded between the maps on a sample of the cells (the paper's 5 %, drawn from the expansion cells and the others separately); the demand is the Markov chain of the two maps projected (PLUS's default) unless demand is given; cars are CARS's settings.

molusce

molusce(
    t0,
    t1,
    drivers,
    *,
    bins=None,
    iterations: int = 1,
    unit_cell: float = 1,
    dynamic=None,
    nodata=None,
)

MOLUSCE: its weights of evidence (MolusceWoE) and simulator (MolusceCA).

Each iteration repeats the change observed between the maps (their cross table); start the run from the later map, as MOLUSCE does.

Running models with Terraflux

Terraflux runs land change models as graphs of operators, with the same C++ kernels: Model.save writes a fitted model as its graph (with the data it runs on), Model.load reads it back, and terraflux runs it with the terraflux command, giving the same maps as Model.run with the same seed.

terraflux

terraflux(model, *, out=None, seed=None, exe=None, quiet: bool = True) -> DataArray

Run a Terraflux model with the terraflux command; the simulated maps (step, y, x).

Parameters:

Name	Type	Description	Default
model		The model file (Model.save writes one).	required
out		The folder of the outputs (default: output next to the model).	None
seed		Replace the model's seed.	None
exe		The terraflux executable (default: the TERRAFLUX variable, else the one on PATH).	None

Returns:

Type	Description
DataArray	The maps output: one map per step (1, 2...), the start not included.

Scenarios and replicates

scenarios

scenarios(
    model,
    start=None,
    drivers=None,
    *,
    scenarios=None,
    replicates: int = 10,
    seed=None,
    steps=None,
    years=None,
    keep=None,
    n_workers=None,
    n_threads=None,
    nodata=None,
    on=None,
) -> Ensemble

Run a model under several scenarios, several times each, in parallel.

Parameters:

Name	Type	Description	Default
model		A Model, or a Projection (its fitted model, from its later map with its drivers, over its years).	required
start		The starting map (default: the projection's later map).	None
drivers		The drivers (default: the projection's).	None
scenarios		{name: change}: a Model, a dict of the pieces to replace ({"demand": ..., "rules": ...}) or a function of the model; default one scenario, "baseline", the model as given.	None
replicates	int	Runs per scenario; replicate r uses the seed seed + r in every scenario (common random numbers: scenarios differ by their settings, not by their draws).	10
steps		Steps to simulate, set on the allocation (default: the model's own).	None
years		Year of the start (one step a year), or the years of all the maps.	None
keep		The years to keep (default all): with large maps, e.g. [2030, 2050].	None
n_workers		Runs at the same time (default: as many as the threads allow).	None
n_threads		Threads in all (default: all cores but one), shared by the workers.	None
on		ls.helpers(): run the simulations on this computer and on the helpers of its group (see parallel), with the same maps. Each scenario's model travels as its Terraflux graph (save()), so it must be made of the pieces that save() writes.	None

Returns:

Type	Description
Ensemble

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> S = ls.simulate
>>> land = np.ones((30, 30)); land[10:20, 10:20] = 2
>>> pot = S.RasterPotential({1: np.full((30, 30), 0.2), 2: np.random.default_rng(0).random((30, 30))})
>>> demand = S.AreaDemand(np.array([[800, 100], [750, 150], [700, 200]]), classes=[1, 2])
>>> fast = S.AreaDemand(np.array([[800, 100], [700, 200], [600, 300]]), classes=[1, 2])
>>> ens = S.scenarios(S.Model(pot, demand, S.CARS(step_size=10)), land, replicates=4, seed=1,
...                   scenarios={"trend": {}, "fast": {"demand": fast}}, years=2020)
>>> ens.probability(2, year=2022).sel(scenario="fast").values.sum().round()
np.float64(300.0)

Ensemble dataclass

Ensemble(
    scenarios: list,
    seeds: list,
    years: list,
    classes: list,
    codes: ndarray,
    steps: DataFrame,
    template: object = None,
)

The simulations of scenarios(): every scenario, replicate and year.

Attributes:

Name	Type	Description
scenarios	list of str
seeds	list of int	The seed of each replicate (the same in every scenario).
years	list	The year (or step) of each kept map, the start first.
classes	list	The classes of the maps.
steps	DataFrame	The step tables of every run, with scenario and replicate.

maps

maps(scenario=None, replicate: int = 0) -> DataArray

The maps of one run, (year, y, x).

probability

probability(cls, *, year=None, scenario=None) -> DataArray

Share of the replicates where each cell has the class: (scenario, year, y, x).

year and scenario select one (the dimension is dropped).

change_probability

change_probability(*, scenario=None) -> DataArray

Share of the replicates where each cell's class differs from the start: (scenario, year, y, x).

cover

cover() -> DataFrame

Cells of each class: a row per scenario, replicate and year.

metrics

metrics(
    metrics, *, years=None, level: str = "class", n_threads: int | None = None
) -> DataFrame

Landscape metrics (or any function of a map) of every run and year, in long format.

Parameters:

Name	Type	Description	Default
metrics		Names of ls.landscape_metrics metrics (["pland", "np", "ai"], at level), or a dict {name: function} of functions of a map (a DataArray) giving a number, a dict or a Series (e.g. an MSPA share or a HabitatGraph index).	required
years		The years to measure (default: all kept).	None
level	str	The level of the named metrics ("class", "landscape").	'class'

Returns:

Type	Description
DataFrame	scenario, replicate, year, metric, class (NaN at landscape level and for functions), value.

to_dataset

to_dataset() -> Dataset

Every map as one Dataset: landuse (scenario, replicate, year, y, x).

Example data

example_data

example_data(name: str = 'sibuyan_landuse') -> Path

Path of an example data set for land change modelling.

Parameters:

Name	Type	Description	Default
name	str	Plum Island Ecosystems, Massachusetts (Pontius & Parmentier 2014; 434 x 497 cells of ~100 m, EPSG:26986): "pie_landuse": land use in 1985, 1991 and 1999 (bands named by year; 1 forest, 2 built, 3 other; 255 = no data); "pie_drivers": elevation, slope and distance to built land in 1985. Sibuyan Island, Philippines (Verburg et al. 2002, the CLUE-S demonstration data, modified for teaching: not for other uses; 108 x 128 cells of 250 m, no CRS): "sibuyan_landuse": land use in 1997 and 2011 (1 forest, 2 coconut, 3 grassland, 4 rice, 5 other); "sibuyan_drivers": 13 drivers (distance to sea, population density, diorite, ultramafic rock, sediments, no erosion, moderate erosion, elevation, slope, aspect, distance to roads, to urban areas and to streams); "sibuyan_restrictions": current and proposed national park; "sibuyan_demand": CSV with the demand scenarios demand1 to demand4 (cells of each land use per year). Sao Lourenco River Basin, Mato Grosso, Brazil ("sao_lourenco"; SOS Pantanal 2015; 6372 x 6546 cells of 30 m, UTM 21S): land use and cover in 2002, 2008, 2010, 2012 and 2014 (2 pasture, 3 forest, 4 park savannah, 5 gramineous savannah, 7 anthropized vegetation, 8 wooded savannah, 9 water, 10 urban, 11 crop farming, 12 reforestation, 13 mining; 255 = no data). Downloaded (5 MB) and converted the first time it is asked for (the conversion needs ~2 GB of memory once).	'sibuyan_landuse'

Returns:

Type	Description
Path

Examples:

>>> import landschaft as ls
>>> lu = ls.load_raster(ls.simulate.example_data("sibuyan_landuse"))
>>> lu.shape
(2, 108, 128)