Spatial prioritization

Two ways of setting conservation priorities, side by side: ranking every cell of a landscape (below), and selecting the best set of planning units by optimization, with targets, budgets, spatial penalties and constraints, portfolios, multi-action plans and the evaluation of the solutions (Optimization onwards). Both read the same rasters and share their pre-processing, and each can start or check the other.

Ranking the cells of a landscape from the least to the most important for conservation, given where a set of features (species, habitats, ecosystem services, carbon...) occur: the core of systematic conservation planning. The ranking balances all features at once, weighted, with costs and hierarchical masks (for example, protected areas kept to the end), and its performance curves tell how much of every feature the best 10 %, 20 %... of the landscape hold. It lives in the ls.prioritize module.

The method is Zonation 5 (Moilanen et al. 2022): the cells are sorted by their marginal loss, the loss to the features of removing a cell at the point where it is removed, with the core-area (CAZ, driven by the features that lose the largest fraction) or additive-benefit (ABF) rules, and sorted again until the order is consistent. Landschaft gives the rank maps and curves of the z5 program (v2.4) bit for bit, faster, with the curves of the ranking it returns (Zonation's own report those of the order before its last sort). The tutorial walks through a project, and the course chapter Conservation priorities explains the method.

import landschaft as ls

res = ls.prioritize.zonation(["sp1.tif", "sp2.tif", "sp3.tif"], weights=[2, 1, 1],
                             cost="land_price.tif", hierarchic_mask="protected.tif")
res.rank                      # 0 = removed first ... ~1 = highest priority
res.feature_curves            # fraction of every feature kept, by fraction removed
res.write("zonation_output")  # rankmap.tif and Zonation's CSV files

Ranking

zonation

zonation(
    features,
    *,
    weights=None,
    rule: str = "CAZ2",
    rule_param: float | None = None,
    abf_z=None,
    names=None,
    weight_groups=None,
    group_weights: dict | None = None,
    groups=None,
    thresholds=None,
    zero_mode: str = "all",
    keep_zeros=None,
    extract=None,
    cost=None,
    cost_in_ranking: bool = True,
    hierarchic_mask=None,
    analysis_area=None,
    area: int = 0,
    solution=None,
    bits: int = 8,
    max_iterations: int = 500,
    ranking_error: float = 0.01,
    smoothing_factor: float = 0.0,
    stop_on_plateau: bool = False,
    curves_interval: float = 0.0001,
    seed: int | None = None,
    condition=None,
    renormalize=False,
    connectivity=None,
    self_multiply=False,
    interaction=None,
    interaction2=None,
    matrix_connectivity=None,
    matrices: dict | None = None,
    retention=None,
    preprocessing: bool = False,
    transformed=False,
    impact_groups=None,
    impact_bins=None,
    impact=None,
    cumulative_impact=None,
    nodata=None,
    storage=None,
    n_threads: int | None = None,
) -> ZonationResult

Rank the cells of a landscape for conservation, as Zonation 5 does.

Every feature (species distribution, habitat, ecosystem service...) is normalized to sum 1 and quantized to bits bits; the cells are then ordered from the least to the most important by their marginal loss: the loss to the features, weighted, of removing a cell at the point where it is removed, with the other cells removed before it. Zonation iterates stable sorts of the cells by that loss until the order is consistent (the mean ranking error below ranking_error %). The result is the same, bit for bit, as the z5 program (v2.4) on one thread.

Parameters:

Name	Type	Description	Default
features		Raster files (every band is a feature), 2-D arrays or DataArrays, or a list of them; or a 3-D array/DataArray (features, rows, columns), NumPy or dask. All on the same grid.	required
weights		One weight per feature (Zonation's -w); negative weights mark features to avoid. None: all 1.	None
rule	str	Marginal loss rule: "CAZ2" (default), "CAZ1", "CAZMAX" (core-area Zonation: the loss is driven by the features that lose the largest fraction), "CAZP"/"CAZPW" (p-norms), "ABF" (additive benefit function, benefit R^z of what remains), "RAND" (random order) or "LOAD" (the order of solution, to evaluate a given ranking). A parameter may follow the name ("ABF 0.25").	'CAZ2'
rule_param	float | None	ABF exponent z (0.25) or CAZP/CAZPW exponent p (1).	None
abf_z		ABF exponent per feature (None entries use rule_param).	None
names		Feature names for the tables.	None
weight_groups		Weight group of each feature (-1: none) and the total weight of each group ({group: weight}): the weights of a group are rescaled to sum (in absolute value) to its weight (Zonation's -W).	None
group_weights		Weight group of each feature (-1: none) and the total weight of each group ({group: weight}): the weights of a group are rescaled to sum (in absolute value) to its weight (Zonation's -W).	None
groups		Output group of each feature (-1: none), summarized in group_curves.	None
thresholds		Per feature, values below become 0.	None
zero_mode	str	"all" (zeros are data), "none" (zeros are NoData) or "dynamic" (per feature, from keep_zeros).	'all'
extract		Per feature, a value to extract from a categorical raster (cells with it keep it, the others become NoData; None: no extraction).	None
cost		Cost of every cell. With cost_in_ranking the losses are divided by the relative cost (-X); otherwise the cost is only followed in the curves (-x). Cells of cost 0 are removed last.	None
hierarchic_mask		Integer levels 0-255 per cell: lower levels are removed first.	None
analysis_area		Raster of analysis areas: cells outside (0, or other than area when area > 0) are left out.	None
area		Raster of analysis areas: cells outside (0, or other than area when area > 0) are left out.	None
solution		The ranking to evaluate with rule="LOAD".	None
bits	int	Quantization of the feature values: 8 (z5) or 16 (z5_16bit).	8
max_iterations	int	The iteration: at most max_iterations sorts, until the mean ranking error is at most ranking_error % and no longer decreasing; sort by a running mix of the losses (smoothing_factor of the previous ones); stop when the error stops decreasing after 5 iterations.	500
ranking_error	int	The iteration: at most max_iterations sorts, until the mean ranking error is at most ranking_error % and no longer decreasing; sort by a running mix of the losses (smoothing_factor of the previous ones); stop when the error stops decreasing after 5 iterations.	500
smoothing_factor	int	The iteration: at most max_iterations sorts, until the mean ranking error is at most ranking_error % and no longer decreasing; sort by a running mix of the losses (smoothing_factor of the previous ones); stop when the error stops decreasing after 5 iterations.	500
stop_on_plateau	int	The iteration: at most max_iterations sorts, until the mean ranking error is at most ranking_error % and no longer decreasing; sort by a running mix of the losses (smoothing_factor of the previous ones); stop when the error stops decreasing after 5 iterations.	500
curves_interval	float	Fraction of the cells between two rows of the curves (0.0001: 10 000 rows).	0.0001
seed	int | None	Seed of RAND.	None
condition		Per feature, a raster of the condition of its habitat, 0-1 (None: none; Zonation's -c): the quantized values are multiplied by it (truncated), and the feature's initial level becomes the conditioned fraction.	None
renormalize		Per feature (or for all), rescale a conditioned feature to level 1.	False
connectivity		Per feature, a Kernel (None: none; Zonation's -d): the feature is replaced by its convolution with the kernel (its connectivity), kept where it had data, and quantized again.	None
self_multiply		Per feature (or for all): multiply the smoothed feature by the feature.	False
interaction		Per feature, an Interaction (-i, -I): the feature weighted by (or reduced by) the connectivity of a resource. The second set is applied after the matrix transforms.	None
interaction2		Per feature, an Interaction (-i, -I): the feature weighted by (or reduced by) the connectivity of a resource. The second set is applied after the matrix transforms.	None
matrix_connectivity		Per feature, (kernel, key) (-m): the features sharing a key are smoothed and mixed by the matrix matrices[key] (row: source feature, column: target feature, in the features' order), then multiplied by the target feature.	None
matrices		Per feature, (kernel, key) (-m): the features sharing a key are smoothed and mixed by the matrix matrices[key] (row: source feature, column: target feature, in the features' order), then multiplied by the target feature.	None
retention		Per feature (Zonation's -r): None; a number, the fraction of the feature retained everywhere; a raster of fractions retained per cell; or ("absolute", raster), the amounts retained per cell. What is retained is taken out of the feature (it stays whatever the ranking), so its curve ends at the retained fraction.	None
preprocessing	bool	Also compute Zonation's pre-processing maps (-b): per cell, the sum of the features' original values (richness) and of their normalized values (range-size rarity, rsr), for the positive, negative and zero weights and each output group, and weighted with weights.	False
transformed		Per feature (or True for all), return the feature as ranked (-f).	False
impact_groups		Per feature, its impact group (-1: none), for the impact scenarios.	None
impact_bins		Per feature, the value classes its impacts are reported in: None (none), 0 (ten equal classes up to its maximum) or a list of class limits in the feature's units.	None
impact		An ImpactScenario (-s: land-use change around impact sources) and/or a CumulativeScenario (-e: direct habitat loss and cumulative pressures). The features of the impact groups lose what the scenarios take before the ranking; the impacts are in impacts.	None
cumulative_impact		An ImpactScenario (-s: land-use change around impact sources) and/or a CumulativeScenario (-e: direct habitat loss and cumulative pressures). The features of the impact groups lose what the scenarios take before the ranking; the impacts are in impacts.	None
nodata		NoData value of inputs without one (NaN is always NoData).	None
storage		A directory for temporary files: the large arrays (the quantized features, the rows of the ranking and its work arrays) are then memory-mapped files that the operating system pages to disk when they do not fit in memory, for maps larger than the RAM (out-of-core). The files are deleted when no longer needed. None: in memory.	None
n_threads	int | None	Threads (default: all cores but one).	None

Returns:

Type	Description
ZonationResult

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> rng = np.random.default_rng(0)
>>> species = rng.random((3, 20, 30)) ** 3          # three features
>>> res = ls.prioritize.zonation(species, weights=[2, 1, 1])
>>> res.rank.shape
(20, 30)
>>> float(res.rank.max()) < 1.0
True
>>> round(float(res.feature_curves.iloc[0, 1]), 4)  # nothing removed yet
1.0

ZonationResult dataclass

ZonationResult(
    rank: object,
    feature_curves: DataFrame,
    summary_curves: DataFrame,
    density_curves: DataFrame,
    group_curves: DataFrame | None,
    features: DataFrame,
    iterations: int,
    ranking_error: float,
    auc: float,
    history: DataFrame,
    order: ndarray,
    info: dict = dict(),
    preprocessing: dict = dict(),
    transformed: dict = dict(),
    hierarchic_mask_ranks: dict = dict(),
    impacts: dict = dict(),
    _ctx: dict = dict(),
)

The result of zonation().

Attributes:

Name	Type	Description
rank	DataArray	Priority rank of every cell with data, from 0 (first removed, lowest priority) to (n - 1) / n (last removed, highest priority); NaN elsewhere.
feature_curves	DataFrame	Fraction remaining of every feature (columns) as the lowest-ranked cells are removed (rank = fraction of the cells removed); Zonation's feature_curves.csv.
summary_curves, density_curves	DataFrame	Mean, weighted mean, maximum and minimum of the curves (by the sign of the weights), remaining area and cost; and the loss rates of the means.
group_curves	DataFrame or None	Mean, weighted mean and minimum per output group.
features	DataFrame	Per feature: input number, name, weight, the sum of its values, cells with values above zero, cells without data, quantization error (%).
iterations	int	Sorting iterations run (0 for RAND and LOAD).
ranking_error	float	Mean ranking error of the last iteration (%).
auc	float	Area under the mean curve of the positively weighted features.
history	DataFrame	Ranking error (%), number of cells out of order and mean rank change (%) per iteration.
order	ndarray	Indices (row-major) of the cells in removal order.
preprocessing	dict	With preprocessing=True: Zonation's pre-processing maps (richness.pos, rsr.pos, ... .weighted; see zonation()).
transformed	dict	With transformed: the feature layers as ranked (after quantization, condition and retention), normalized to sum 1, by feature name.
hierarchic_mask_ranks	dict	With a hierarchic mask: the lowest rank of each of its levels.
impacts	dict	With impact scenarios: zones (per area, group, source, zone, feature and bin: the features' distribution in the zone, its loss and cells), cumulative (per area, group, phase and feature) and rasters (the impact of each group).

area_statistics

area_statistics(areas)

Statistics of areas of the landscape: their size, mean rank, the share of every feature they hold and how rich they are relative to their size (Zonation's area_statistics_external).

Parameters:

Name	Type	Description	Default
areas		Raster of area numbers (0 or NoData: none) on the features' grid.	required

Returns:

Type	Description
AreaStatistics

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> x = np.random.default_rng(0).random((3, 20, 30))
>>> res = ls.prioritize.zonation(x)
>>> zones = np.repeat([[1, 2, 3]], 20 * 10, axis=0).reshape(20, 30)
>>> st = res.area_statistics(zones)
>>> st.areas[["area", "size"]].values.tolist()
[[1, 200], [2, 200], [3, 200]]

area_statistics_rank_range

area_statistics_rank_range(min_rank: float, max_rank: float = 1.0)

area_statistics() of the cells ranked in [min_rank, max_rank] (Zonation's area_statistics_rank_range); the mask is in the result.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> res = ls.prioritize.zonation(np.random.default_rng(0).random((3, 20, 30)))
>>> top = res.area_statistics_rank_range(0.9)      # about the best 10 %
>>> int(top.areas["size"].iloc[0])                  # 540/600 in float32 is just below 0.9
59

area_statistics_separate_areas

area_statistics_separate_areas(
    min_rank: float, max_rank: float = 1.0, min_area_size: int = 1
)

area_statistics() of the separate (8-connected) areas of cells ranked in [min_rank, max_rank] of at least min_area_size cells (Zonation's area_statistics_separate_areas).

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> res = ls.prioritize.zonation(np.random.default_rng(0).random((3, 20, 30)))
>>> st = res.area_statistics_separate_areas(0.8, 1.0, min_area_size=3)
>>> bool((st.areas["size"] >= 3).all())
True

solution_comparison

solution_comparison(
    other, min_rank: float, max_rank: float, min_value: float, max_value: float
)

Compare the cells ranked in [min_rank, max_rank] with the cells of another solution with values in [min_value, max_value] (Zonation's solution_comparison).

Returns:

Type	Description
dict	both, only_this, only_other, neither (cell counts) and map (1 both, 2 only this, 3 only the other, 4 neither, 0 no data).

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> res = ls.prioritize.zonation(np.random.default_rng(0).random((3, 20, 30)))
>>> cmp = res.solution_comparison(res.rank, 0.5, 1.0, 0.5, 1.0)   # itself
>>> cmp["both"], cmp["only_this"], cmp["only_other"]
(300, 0, 0)

histogram

histogram(min_rank: float, max_rank: float = 1.0, bins: int = 10) -> dict

How many features keep which fraction in the cells ranked in [min_rank, max_rank]: per weight class and output group, the count of features per bin of coverage (Zonation's histogram).

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> res = ls.prioritize.zonation(np.random.default_rng(0).random((3, 20, 30)))
>>> h = res.histogram(0.9, bins=10)["pos"]
>>> int(h["count"].sum())
3

flexibility_map

flexibility_map(max_error: float = 2.0, *, up: bool = True, down: bool = True) -> dict

How far each cell could move in the ranking before its marginal loss, at the other position, differs from that of the cell there by max_error % (Zonation's flexibility_map): up, down and their sum both, as fractions of the cells.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> res = ls.prioritize.zonation(np.random.default_rng(0).random((3, 20, 30)))
>>> f = res.flexibility_map(5.0)
>>> sorted(f)
['both', 'down', 'up']

postprocess

postprocess(calls, directory) -> list

Run Zonation's post-processing calls and write its files into directory.

Parameters:

Name	Type	Description	Default
calls		A post-processing file (its paths are relative to it), its text, or a list of (name, [parameters]): area_statistics_external(file, id), area_statistics_rank_range(min_rank, max_rank, id), area_statistics_separate_areas(min_rank, max_rank, min_area_size, id), solution_comparison(min_rank, max_rank, file, min_value, max_value, id), histogram(min_rank, max_rank, bins, id), flexibility_map(max_error, both, up, down, id).	required
directory		Output directory.	required

Returns:

Type	Description
list of pathlib.Path	The files written.

Examples:

>>> import tempfile
>>> import numpy as np
>>> import landschaft as ls
>>> res = ls.prioritize.zonation(np.random.default_rng(0).random((3, 20, 30)))
>>> files = res.postprocess('histogram(0.9, 1.0, 10, "top")', tempfile.mkdtemp())
>>> [f.name for f in files]
['histogram.pos.top.txt']

save

save(directory) -> Path

Save the result as Zonation 5 does (the same as write()).

Examples:

>>> import tempfile
>>> import numpy as np
>>> import landschaft as ls
>>> res = ls.prioritize.zonation(np.random.default_rng(0).random((2, 10, 12)))
>>> (res.save(tempfile.mkdtemp()) / "rankmap.tif").exists()
True

write

write(directory) -> Path

Write the rank map and the tables as Zonation 5 does.

Parameters:

Name	Type	Description	Default
directory		Output directory (created if needed): rankmap.tif, feature_curves.csv, summary_curves.csv, feature_density_curves.csv, group_curves.csv (with output groups) and features_info.csv, in Zonation's formats.	required

Returns:

Type	Description
Path	The directory.

Examples:

>>> import tempfile
>>> import numpy as np
>>> import landschaft as ls
>>> res = ls.prioritize.zonation(np.random.default_rng(0).random((2, 10, 12)))
>>> out = res.write(tempfile.mkdtemp())
>>> sorted(p.name for p in out.iterdir())
['feature_curves.csv', 'feature_density_curves.csv', 'features_info.csv', 'rankmap.tif', 'summary_curves.csv']

Large maps

The features are read in blocks of rows (files, NumPy, xarray or dask arrays; dask stacks are computed block by block, never whole) and quantized to one or two bytes per value, kept as runs of NoData, zeros and values: a feature costs about its occurrences, however many cells it is zero in, so thousands of species ranges fit where Zonation stores every feature in every cell. During the ranking every cell's occurrences are a compact row (a bitmap or layer gaps, whichever is smaller) and the work arrays take about 130 bytes per cell with data.

For maps larger than the memory, storage puts the large arrays in memory-mapped temporary files in a directory (on a fast disk), which the operating system pages as needed; the reads and writes of the ranking are in sequence or nearly so. The result is the same bit for bit; the files are deleted when no longer needed.

import dask.array as da

stack = da.stack([da.from_zarr(f"ranges.zarr/{s}") for s in species])  # (features, y, x)
res = ls.prioritize.zonation(stack, storage="D:/scratch", n_threads=16)
res.write("global_ranking")

The ranking is limited to 2^32 - 1 cells with data and 65 535 features.

Zonation projects

Existing Zonation 5 projects run as they are: the settings file, its feature list and link files (weight groups, conditions, retention), with the switches of the z5 command line.

run_project

run_project(settings, flags='', *, mode: str | None = None, output=None, **overrides)

Run a Zonation 5 project (see read_project()).

Parameters:

Name	Type	Description	Default
settings		As read_project().	required
flags		As read_project().	required
mode		As read_project().	required
output		Directory to write Zonation's output files into (optional).	None
**overrides		Other zonation() arguments (bits=16 for z5_16bit, n_threads, storage...).	{}

Returns:

Type	Description
ZonationResult

Examples:

>>> import tempfile
>>> from pathlib import Path
>>> import numpy as np, rasterio
>>> import landschaft as ls
>>> d = Path(tempfile.mkdtemp())
>>> for i in range(2):
...     with rasterio.open(d / f"sp{i}.tif", "w", driver="GTiff", height=5, width=6, count=1,
...                        dtype="float32", crs="EPSG:32633",
...                        transform=rasterio.transform.from_origin(0, 500, 100, 100)) as f:
...         f.write(np.random.default_rng(i).random((5, 6), dtype=np.float32), 1)
>>> _ = (d / "features.txt").write_text("filename\nsp0.tif\nsp1.tif\n")
>>> _ = (d / "settings.z5").write_text("feature list file = features.txt\n")
>>> res = ls.prioritize.run_project(d / "settings.z5", mode="CAZ2")
>>> res.rank.shape
(5, 6)

read_project

read_project(settings, flags='', *, mode: str | None = None) -> dict

The zonation() arguments of a Zonation 5 project.

Parameters:

Name	Type	Description	Default
settings		The project's settings file (.z5).	required
flags		The z5 command-line switches that enable its parts, as a string or list: -w weights, -W weight groups, -g output groups, -t thresholds, -c conditions, -r retention, -x/-X cost (tracked / in the ranking), -h hierarchic mask, -a analysis area (--area=N), -b pre-processing maps, -f transformed layers, -d/-i/-I/-m connectivity, -s/-e impact scenarios, -p post-processing (run by run_project() with an output).	''
mode	str | None	The marginal loss rule (z5's --mode), e.g. "CAZ2" or "ABF 0.25".	None

Returns:

Type	Description
dict	Keyword arguments of zonation() (features included; with -p, postprocessing_file too, which run_project() runs).

Examples:

>>> import tempfile
>>> from pathlib import Path
>>> import numpy as np, rasterio
>>> import landschaft as ls
>>> d = Path(tempfile.mkdtemp())
>>> for i in range(2):
...     with rasterio.open(d / f"sp{i}.tif", "w", driver="GTiff", height=5, width=6, count=1,
...                        dtype="float32", crs="EPSG:32633",
...                        transform=rasterio.transform.from_origin(0, 500, 100, 100)) as f:
...         f.write(np.random.default_rng(i).random((5, 6), dtype=np.float32), 1)
>>> _ = (d / "features.txt").write_text("weight filename\n2.0 sp0.tif\n1.0 sp1.tif\n")
>>> _ = (d / "settings.z5").write_text("feature list file = features.txt\n")
>>> kw = ls.prioritize.read_project(d / "settings.z5", "-w", mode="ABF")
>>> kw["weights"], kw["rule"]
([2.0, 1.0], 'ABF')

Connectivity

Dispersal kernels for the connectivity transforms of the features (the connectivity, interaction and matrix_connectivity arguments of zonation): a feature smoothed by the dispersal of its species, weighted by the connectivity of a resource, or mixed with others by a matrix.

Kernel dataclass

Kernel(shape: str, distance: float, cut: float = 1.0)

A dispersal kernel of Zonation 5's connectivity transforms.

Parameters:

Name	Type	Description	Default
shape	str	"uniform" (1 within distance), "triangular" (1 - d / distance within distance), "negexp" (exp(-2 d / distance): distance is the mean dispersal distance) or "gaussian" (normal density of standard deviation distance).	required
distance	float	In cells.	required
cut	float	negexp and gaussian: the kernel is 0 beyond distance * cut.	1.0

Examples:

>>> import landschaft as ls
>>> k = ls.prioritize.Kernel("negexp", 10, cut=3)
>>> k.image().shape       # the kernel image Zonation convolves with
(128, 128)

size

size(max_distance: int = 2 ** 31 - 1) -> int

The side of the kernel image (getKernel's getFilterSize).

image

image(max_distance: int = 2 ** 31 - 1) -> ndarray

The kernel image (float64, not normalized).

Interaction dataclass

Interaction(
    kernel: Kernel,
    resource: object,
    positive: bool = True,
    beta: float = 1.0,
    normalize: bool = True,
)

An interaction connectivity transform (Zonation's -i/-I): a feature weighted by the connectivity of a resource (positive) or reduced by it.

Parameters:

Name	Type	Description	Default
kernel	Kernel	Dispersal kernel from the resource.	required
resource	object	Raster of the resource (path, DataArray or array).	required
positive	bool	True: the feature is multiplied by the smoothed resource; False: beta times the feature's share of the product is removed.	True
beta	float	Strength of the negative interaction.	1.0
normalize	bool	Negative interactions: keep the feature's level (False: its level becomes what remains).	True

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> rng = np.random.default_rng(0)
>>> species = rng.random((2, 30, 40))
>>> food = rng.random((30, 40))          # a resource of species 0
>>> ia = ls.prioritize.Interaction(ls.prioritize.Kernel("negexp", 3, 3), food)
>>> res = ls.prioritize.zonation(species, interaction=[ia, None])
>>> res.iterations > 0
True

convolve

convolve(image: ndarray, kernel: Kernel, workers: int = 1) -> ndarray

Convolve a raster with a kernel normalized to sum 1, as Zonation 5 does.

FFT convolution in blocks of the kernel's size, the blocks' results added as float32 in Zonation's order: the smoothing of its connectivity transforms.

Parameters:

Name	Type	Description	Default
image	ndarray	2-D array (rows, columns); cells outside the data should be 0.	required
kernel	Kernel	A Kernel.	required
workers	int	FFT threads.	1

Returns:

Type	Description
ndarray	float32, the image's shape.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> img = np.zeros((50, 60))
>>> img[25, 30] = 1.0
>>> out = ls.prioritize.convolve(img, ls.prioritize.Kernel("uniform", 2))
>>> round(float(out.sum()), 5), int((out > 1e-6).sum())   # mass kept, spread over the disc
(1.0, 13)

Post-processing

Analyses of a ranking (methods of ZonationResult): statistics of areas (external zones, the top fraction, separate high-priority areas), comparison with another solution, histograms of the features' coverage and the flexibility map; ZonationResult.postprocess runs Zonation's post-processing files.

AreaStatistics dataclass

AreaStatistics(areas: DataFrame, coverage: dict, mask: ndarray | None = None)

Bases: Saveable

Statistics of areas of a ranking (area_statistics()).

Attributes:

Name	Type	Description
areas	DataFrame	Per area (and weight class or output group, column features): size in cells, mean rank, the sum of the features' normalized values (ds) and weighted (wds), relative to what the same area would hold on average (rfd, wfd), and the cost.
coverage	dict	Per weight class / group, a DataFrame: the fraction of every feature in every area.
mask	ndarray or None	The areas (0: outside), when they come from the ranking.

Impact assessment

Scenarios of land-use change and pressures, evaluated for the features and applied to them before the ranking (the impact, cumulative_impact, impact_groups and impact_bins arguments of zonation): the loss around impact sources, per footprint and buffer zone, and the direct habitat loss and cumulative pressures of every impact group.

ImpactScenario dataclass

ImpactScenario(
    sources: object,
    impacts: list,
    name: str = "scenario",
    subareas: object = None,
    area_unit: str = "",
    impact_unit: str = "",
    area_factor: float = 1.0,
    rasters: bool = True,
)

An impact assessment scenario (Zonation's -s): land-use change at sources.

Parameters:

Name	Type	Description	Default
sources	object	Raster of the impact sources (integer values; 0: none).	required
impacts	list	The ImpactSource of every group and source.	required
name	str	Project id (the output directory's name).	'scenario'
subareas	object	Raster of subareas (integers, 0: none) reported apart.	None
area_unit	str	Names of the units of area and of the features' values, and the area of a cell in the area unit.	''
impact_unit	str	Names of the units of area and of the features' values, and the area of a cell in the area unit.	''
area_factor	str	Names of the units of area and of the features' values, and the area of a cell in the area unit.	''
rasters	bool	Write the impact rasters of the groups.	True

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> P = ls.prioritize
>>> species = np.random.default_rng(0).random((2, 30, 40))
>>> road = np.zeros((30, 40), np.uint16)
>>> road[15, :] = 1                                   # source 1: a planned road
>>> sc = P.ImpactScenario(road, [P.ImpactSource(group=1, source=1, footprint=100,
...                                             buffers=((3, 50),))])   # 50 % loss up to 3 cells
>>> res = P.zonation(species, impact_groups=[1, -1], impact=sc)
>>> zones = res.impacts["zones"]
>>> sorted(zones["zone"].unique().tolist())          # footprint and buffer ring
[0, 1]

ImpactSource dataclass

ImpactSource(group: int, source: int, footprint: float, buffers: tuple = ())

An impact source of an impact group (a line of Zonation's impact buffers file).

Parameters:

Name	Type	Description	Default
group	int	Impact group (the features with this impact_groups entry).	required
source	int	Value of the source's cells in the impact source raster.	required
footprint	float	Loss (%) in the source's own cells.	required
buffers	tuple	(radius, loss %) of successive rings around the source (radius in cells, increasing).	()

CumulativeScenario dataclass

CumulativeScenario(
    name: str = "cumulative",
    subareas: object = None,
    area_unit: str = "",
    impact_unit: str = "",
    area_factor: float = 1.0,
    rasters: bool = True,
    direct: str | None = None,
    direct_rasters: list = list(),
    direct_matrix: object = None,
    cumulative: str | None = None,
    pressures: list = list(),
    cumulative_matrix: object = None,
    aggregation: str = "add",
    functions: dict = dict(),
)

A cumulative impact assessment scenario (Zonation's -e).

Parameters:

Name	Type	Description	Default
name	str	As ImpactScenario.	'cumulative'
subareas	str	As ImpactScenario.	'cumulative'
area_unit	str	As ImpactScenario.	'cumulative'
impact_unit	str	As ImpactScenario.	'cumulative'
area_factor	str	As ImpactScenario.	'cumulative'
rasters	str	As ImpactScenario.	'cumulative'
direct	str | None	Direct habitat loss: None, "sum" (the weighted sum of the loss rasters, at most 1) or "product" (1 - the product of what the rasters leave).	None
direct_rasters	list	Loss rasters (fractions) and the weights (impact groups x rasters).	list()
direct_matrix	list	Loss rasters (fractions) and the weights (impact groups x rasters).	list()
cumulative	str | None	Cumulative pressures: None, "variant1" (the weighted sum of the pressures through the group's response function) or "variant2" (each pressure through its response function, then aggregated).	None
pressures	list	The Pressure rasters.	list()
cumulative_matrix	object	variant1: weights (impact groups x pressures).	None
aggregation	str	variant2: "add", "multiply" or "max".	'add'
functions	dict	Response functions, {(group, pressure): [(x, loss), ...]}: pressure None for the group's default (1-based pressure numbers otherwise).	dict()

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> P = ls.prioritize
>>> species = np.random.default_rng(0).random((2, 30, 40))
>>> pressure = np.linspace(0, 1, 30 * 40, dtype=np.float32).reshape(30, 40)
>>> sc = P.CumulativeScenario(cumulative="variant1", pressures=[P.Pressure(pressure)],
...                           cumulative_matrix=np.array([[1.0]]),
...                           functions={(1, None): [(0, 0), (1, 0.5)]})   # up to 50 % loss
>>> res = P.zonation(species, impact_groups=[1, 1], cumulative_impact=sc)
>>> sorted(res.impacts["cumulative"]["phase"].unique().tolist())
['cumulative']

Pressure dataclass

Pressure(
    raster: object, normalize: bool = False, offset: float = 0.0, divisor: float = 1.0
)

A pressure raster of a cumulative impact scenario.

Optimization

A problem of planning units (raster cells, polygons, H3 cells or table rows), features, an objective, targets, penalties and constraints, compiled to a mixed-integer linear model and solved to a guaranteed gap. The grammar and the models are prioritizr's (Hanson et al.): compile() gives prioritizr 8.0.6's model bit for bit and the solvers reach the same optima as prioritizr on Gurobi. HiGHS is built in, with Landschaft's own engines in front of it for the problems it finds hard (covering problems of millions of units, budgeted objectives, boundary penalties by minimum cuts, multi-action plans); Gurobi is used when installed and asked for. Coming from R or prioripy? The migration guide lists the equivalent of every function.

import landschaft as ls

P = ls.prioritize
cost = ls.load_raster(P.example_data("cost"))
features = ls.load_raster(P.example_data("features"))
protected = ls.load_raster(P.example_data("protected"))
p = (P.problem(cost, features)
       .add_min_set_objective()
       .add_relative_targets(0.2)
       .add_locked_in_constraints(protected)
       .add_boundary_penalties(0.0005)
       .add_binary_decisions()
       .add_highs_solver(gap=0.05))
sol = p.solve()
sol.status, round(sol.gap, 3)

problem

problem(
    x, features=None, *, cost_column=None, rij=None, contribution=None, zones=None
) -> Problem

A conservation planning problem (prioritizr's problem()).

Parameters:

Name	Type	Description	Default
x		The planning units and their cost: a raster (the cells with a cost are the planning units), polygons (a GeoDataFrame) or a table, both with a cost_column (the rows with a cost are the planning units).	required
features		With a raster: the features' rasters on the same grid (a (band, y, x) DataArray, a Dataset or a list of rasters). With polygons: the features' rasters, summed over every polygon weighted by the exact share of every cell it covers (as exactextractr, which prioritizr uses). With a table: the names of its feature columns, or, with rij, a table of the features (id, name).	None
cost_column		The cost column of a table.	None
rij		Marxan-style amounts (pu, species, amount) with a table of planning units that has an id column (with zones, a zone column too).	None
contribution		Management zones: every planning unit goes to at most one zone, each with its own cost (x: a mapping zone -> cost raster, a raster with a band per zone, or cost_column a mapping zone -> column). The features are given once, and contribution says how much of them counts in each zone: a mapping zone -> a number, one number per feature, or a raster (a column with tables) of factors per unit (zones not listed count fully). Instead, features may be a mapping zone -> that zone's own features (prioritizr's zones()).	None
zones		With rij and a zone column: the zones table (id, name), in the order of the cost columns.	None

Examples:

>>> import landschaft as ls
>>> cost = ls.load_raster(ls.prioritize.example_data("cost"))
>>> features = ls.load_raster(ls.prioritize.example_data("features"))
>>> p = (ls.prioritize.problem(cost, features)
...        .add_min_set_objective()
...        .add_relative_targets(0.3)
...        .add_binary_decisions())
>>> sol = p.solve()
>>> sol.status
'optimal'
>>> bool(sol.coverage["met"].all())
True

Problem dataclass

Problem(
    units: Units,
    objective: tuple | None = None,
    targets: tuple | None = None,
    locked_in: object = None,
    locked_out: object = None,
    bounds: tuple = (),
    decisions: tuple = ("binary",),
    solver: tuple = ("highs", {}),
    weights: object = None,
    subregional: tuple | None = None,
    spatial: tuple = (),
    portfolio: tuple | None = None,
    zone_targets: tuple = (),
    mandatory: bool = False,
)

A conservation planning problem (prioritizr's problem()): see problem().

Every add_*_objective that takes a budget takes one number (the most the selected units may cost); with zones, also one budget per zone (a list in the zones' order, or a mapping zone -> budget).

zones property

zones: tuple[str, ...]

The names of the management zones (empty without zones).

add_min_set_objective

add_min_set_objective() -> Problem

Minimum set: meet every target at the least cost.

add_min_shortfall_objective

add_min_shortfall_objective(budget: float) -> Problem

Minimum shortfall: within the budget, the least sum of the targets' shortfalls, each as a share of its target (weighted by add_feature_weights()).

add_min_largest_shortfall_objective

add_min_largest_shortfall_objective(budget: float) -> Problem

Minimum largest shortfall: within the budget, the least of the largest shortfall of any target, as a share of the target.

add_max_utility_objective

add_max_utility_objective(budget: float) -> Problem

Maximum utility: within the budget, the most of the features' amounts (weighted); no targets.

add_max_features_objective

add_max_features_objective(budget: float) -> Problem

Maximum features: within the budget, the most targets met (weighted).

add_max_wtd_sum_objective

add_max_wtd_sum_objective(budget: float) -> Problem

add_max_utility_objective() under prioritizr 9's name (the model is 8.0.6's, with the small cost term that breaks ties, which prioritizr 9 dropped).

add_max_n_targets_met_objective

add_max_n_targets_met_objective(budget: float) -> Problem

add_max_features_objective() under prioritizr 9's name (the model is 8.0.6's, with the small cost term that breaks ties, which prioritizr 9 dropped).

add_max_cover_objective

add_max_cover_objective(budget: float) -> Problem

Maximum cover: within the budget, the most features held at all (any amount above zero, weighted); no targets.

add_max_phylo_div_objective

add_max_phylo_div_objective(budget: float, tree) -> Problem

Maximum phylogenetic diversity: within the budget, the most total length of the branches of a phylogeny with a target met below them (Faith 1992). tree: an ls.bio tree (read_tree()) or a Newick string, whose tips are the features.

add_max_phylo_end_objective

add_max_phylo_end_objective(budget: float, tree) -> Problem

Maximum phylogenetic endemism: as add_max_phylo_div_objective(), every branch's length divided by the total amount of the features below it (Rosauer et al. 2009).

add_min_extinction_risk_objective

add_min_extinction_risk_objective(budget: float, z=0.25, segments: int = 16) -> Problem

Minimum extinction risk (our extension): within the budget, the least sum over the features (weighted) of the risk the species-area relationship gives, 1 - (held / total)^z (Pereira & Daily 2006; z = 0.25 by default, one value or one per feature), as a piecewise-linear function through segments + 1 points of the curve (exact at them, slightly above in between). No targets.

add_min_penalties_objective

add_min_penalties_objective(budget: float) -> Problem

Minimum penalties: meet the targets within the budget with the least penalties (the penalties come in a later phase: until then any solution that meets them is optimal).

add_relative_targets

add_relative_targets(targets, zone=None) -> Problem

Targets as shares (0-1) of every feature's total amount (one value or one per feature).

With zones, the targets are met by the sum over the zones of the features' amounts in them, each counted by its zone's contribution: a share of the feature's total (the base features of the contributions; with every zone's own features, the sum of its totals in the zones, as prioritizr). zone: targets of one zone (shares of the feature's total in that zone) or of a list of zones, added to the others (targets of the same zones replace the earlier ones). A matrix (features x zones) gives every zone its own targets (prioritizr's form).

add_absolute_targets

add_absolute_targets(targets, zone=None) -> Problem

Targets as amounts of every feature (one value or one per feature); zone as in add_relative_targets().

add_manual_targets

add_manual_targets(targets: DataFrame) -> Problem

Targets from a table, one row per target: feature (name), type ("relative" or "absolute"), target and, optionally, sense (">=", the default, "<=" or "="). A feature may have several targets, or none. With zones, an optional zone column: a zone or a list of zones whose amounts the target adds up (default: every zone; a relative target as in add_relative_targets()).

add_subregional_targets

add_subregional_targets(regions, targets: DataFrame) -> Problem

Targets within sub-regions (an extension: prioritizr has none), added to the other targets: one row per target, region, feature (name), target, optionally type ("relative", the default, or "absolute") and sense (">=", "<=", "="). A relative target is a share of the feature's total in the region, over every unit with or without a cost, as the other relative targets are. With zones, an optional zone column (a zone or a list of zones, default every zone), as in add_relative_targets().

Parameters:

Name	Type	Description	Default
regions		The region of every unit: with raster units, a raster on their grid (NaN: no region); with tables and polygons, a column name or one value per row (missing: no region).	required
targets	DataFrame	The targets (a DataFrame or a dict of columns).	required

add_loglinear_targets

add_loglinear_targets(
    lower_bound_amount: float,
    lower_bound_target: float,
    upper_bound_amount: float,
    upper_bound_target: float,
    cap_amount: float | None = None,
    cap_target: float | None = None,
    abundances=None,
) -> Problem

Targets that fall log-linearly with a feature's total amount (Rodrigues et al. 2004): lower_bound_target (a share) up to lower_bound_amount, upper_bound_target from upper_bound_amount, linear in log10(amount) between them; with a cap, features of at least cap_amount get the absolute target cap_target. abundances: the amounts (default: every feature's total).

add_feature_weights

add_feature_weights(weights) -> Problem

Weights of the features (of the targets, when there are targets) in the objectives that add them up: min shortfall, max utility, max features, max cover (ignored, with a warning, by min set and min largest shortfall). With zones and no targets (max utility, max cover): one weight per feature (the same in every zone) or a features x zones matrix.

add_boundary_penalties

add_boundary_penalties(
    penalty: float, edge_factor: float = 0.5, data=None, zones=None
) -> Problem

Penalise the boundary of the selection (more compact solutions): penalty per unit of exposed boundary length. edge_factor: the share of the boundary at the edge of the study area (not shared with another planning unit) that counts. data: a boundary matrix of every unit (shared lengths, every unit's perimeter on the diagonal); by default from the planning units' raster grid or polygons.

With zones: edge_factor one or one per zone, and zones how much the boundary between two zones counts: "same" (the default, prioritizr's identity matrix: the boundary of every zone's units), "any", or a symmetric zones x zones matrix (values from -1 to 1; the diagonal for units of the same zone, the other cells for neighbours in different zones; a DataFrame by zone names works too).

add_connectivity_penalties

add_connectivity_penalties(penalty: float, data, zones=None) -> Problem

Reward selecting connected units: penalty x the connectivity of every pair selected together (a symmetric matrix of every unit; the diagonal rewards a unit alone). With zones, zones as in add_boundary_penalties() (how much a connection between two zones counts).

add_asym_connectivity_penalties

add_asym_connectivity_penalties(penalty: float, data, zones=None) -> Problem

Penalise asymmetric connections that the selection breaks: penalty x the connectivity from a selected unit to one not selected (a matrix of every unit, rows the sources). With zones, zones as in add_boundary_penalties() (it need not be symmetric).

add_linear_penalties

add_linear_penalties(penalty: float, data) -> Problem

Penalise selecting units by a value of each (penalty x value): a raster on the units' grid, a column name, or one value per unit. With zones, penalty one or one per zone (a list or a mapping) and data the same values in every zone or a mapping zone -> values (zones not given: 0) or a raster with a layer per zone.

add_linear_constraints

add_linear_constraints(threshold: float, sense: str, data) -> Problem

The sum over the selected units of a value of each (as add_linear_penalties()) must be >=, <= or = the threshold. With zones, the sum over every zone (data as in add_linear_penalties()).

add_contiguity_constraints

add_contiguity_constraints(data=None, zones=None) -> Problem

The selected units must form one contiguous patch. data: an adjacency matrix of every unit; by default units that share an edge (raster cells, 4 directions) or that touch (polygons).

The model is exact: a spanning arborescence of the selected units, with a flow that rules out cycles. prioritizr's formulation (a tree whose edges always run from a unit to one of a lower index) is a restriction: it misses patches in which a unit has no selected neighbour of a higher index, so it can return a costlier patch than the best, or call a feasible problem infeasible. On HiGHS the problem is solved in three steps: a first patch from that restriction, a bound from directed cuts, then the exact MIP from the first patch. Contiguity stays the hardest spatial rule: give it a time limit on more than a few hundred units.

With zones, zones is a binary symmetric zones x zones matrix: the zones linked through it form a cluster whose units (in any of its zones) must form one patch; a zone with 0 on the diagonal is free. By default ("same") the units of every zone form one patch.

add_neighbors_constraints

add_neighbors_constraints(k: int, clamp: bool = True, data=None, zones=None) -> Problem

Every selected unit must have at least k selected neighbours (clamp: at most as many as it has). data: an adjacency matrix of every unit, by default as add_contiguity_constraints(). With zones, k one or one per zone, and zones a binary zones x zones matrix of the zones whose units count as neighbours of a unit in a zone (default "same": only its own zone's).

add_neighbor_constraints

add_neighbor_constraints(k: int, clamp: bool = True, data=None, zones=None) -> Problem

add_neighbors_constraints() under the name of prioritizr's newer versions.

add_feature_contiguity_constraints

add_feature_contiguity_constraints(data=None) -> Problem

Every feature's units (where it is counted towards its targets) must form one contiguous patch (the expanded formulation: a column per feature and unit; the exact model of add_contiguity_constraints() for every feature). data: an adjacency matrix of every unit for all features, a list of one per feature, or None (as add_contiguity_constraints()).

add_corridor_constraints

add_corridor_constraints(
    terminals, resistance=None, penalty: float = 1.0, data=None
) -> Problem

Corridors (an extension, after prioripy): the terminal units of every network are selected and joined by selected units (a flow from the first terminal to each of the others).

Parameters:

Name	Type	Description	Default
terminals		The terminals of one network (a raster on the units' grid, a column name or one value per unit; > 0 marks a terminal), or a list of networks.	required
resistance		Optional: values per unit (as terminals; an edge gets the mean of its two units) or an ls.cost Transition (raster units; its connections are the edges, each with the reciprocal of its conductance). The flow's penalty x resistance is added to the objective, so corridors run through low resistance.	None
data		An adjacency matrix of every unit (default: as add_contiguity_constraints()).	None

add_locked_in_constraints

add_locked_in_constraints(locked) -> Problem

Planning units that must be selected: a raster (1 = locked), a column name, a boolean array over the planning units or their (0-based) indices. With zones, a mapping zone -> any of these: the units locked into each zone (a unit in one zone at most, and only where it has a cost).

add_locked_out_constraints

add_locked_out_constraints(locked) -> Problem

Planning units that must not be selected (as add_locked_in_constraints(); with zones, a mapping zone -> the units kept out of that zone).

add_mandatory_allocation_constraints

add_mandatory_allocation_constraints() -> Problem

With zones: every planning unit goes to a zone (by default a unit may go to none).

add_manual_bounded_constraints

add_manual_bounded_constraints(
    data: DataFrame | None = None, *, lower=None, upper=None, zone=None
) -> Problem

The selected share of given planning units must lie between a lower and an upper bound (prioritizr's add_manual_bounded_constraints; prioripy's add_continuous_bounds).

Parameters:

Name	Type	Description	Default
data	DataFrame | None	prioritizr's table: a pu column (the 1-based index of the unit among all units: the cell of the raster, row-major, or the row of the table, with or without cost; units without cost are skipped), lower and upper. A unit listed twice keeps its last bounds.	None
lower		The same as values instead of a table: a number for every planning unit, a raster on the planning units' grid, a column of their table, or one value per planning unit (or per cell or row). A unit with NaN in both is left alone; a NaN (or a side not given) takes the decisions' own bound (0, or 1 or the semicontinuous limit).	None
upper		The same as values instead of a table: a number for every planning unit, a raster on the planning units' grid, a column of their table, or one value per planning unit (or per cell or row). A unit with NaN in both is left alone; a NaN (or a side not given) takes the decisions' own bound (0, or 1 or the semicontinuous limit).	None
zone		With zones: the zone of the values (the table takes a zone column instead). A unit without a cost in a zone stays out of it (a lower bound above 0 there is an error).	None

Notes

The bounds overwrite the units' bounds in the order the manual constraints were added, after add_locked_in_constraints() and add_locked_out_constraints(). They must lie within the decisions' bounds, and with binary decisions be 0 or 1 (prioritizr accepts a share there and the solver rounds it): both raise an error when the problem is compiled.

add_manual_locked_constraints

add_manual_locked_constraints(
    data: DataFrame | None = None, *, status=None, zone=None
) -> Problem

Fix the selected share of given planning units (prioritizr's add_manual_locked_constraints): a status of 1 locks a unit in, 0 out, a share in between (proportion decisions) fixes it.

Parameters:

Name	Type	Description	Default
data	DataFrame | None	prioritizr's table: pu (as in add_manual_bounded_constraints()) and status.	None
status		The same as values: a number, a raster on the planning units' grid, a column of their table or one value per planning unit; NaN leaves a unit alone.	None
zone		With zones: as in add_manual_bounded_constraints().	None

add_binary_decisions

add_binary_decisions() -> Problem

Every planning unit is selected or not.

add_proportion_decisions

add_proportion_decisions() -> Problem

Any share (0-1) of every planning unit can be selected (a linear program).

add_semicontinuous_decisions

add_semicontinuous_decisions(upper_limit: float) -> Problem

Any share of every planning unit up to upper_limit can be selected.

add_cuts_portfolio

add_cuts_portfolio(number_solutions: int = 10) -> Problem

The best number_solutions distinct solutions, one solve after another, each forbidding the earlier ones (binary decisions).

add_top_portfolio

add_top_portfolio(number_solutions: int = 10) -> Problem

As add_cuts_portfolio(), every solve to optimality (gap 0).

add_gap_portfolio

add_gap_portfolio(number_solutions: int = 10, pool_gap: float = 0.1) -> Problem

Up to number_solutions distinct solutions whose objective is within pool_gap (relative) of the optimum's.

add_shuffle_portfolio

add_shuffle_portfolio(
    number_solutions: int = 10, remove_duplicates: bool = True, seed: int | None = None
) -> Problem

number_solutions solves of the model with its planning units (and other columns) in random orders: different solutions of equal quality when there are ties (remove_duplicates: keep distinct ones).

add_highs_solver

add_highs_solver(
    gap: float = 0.1,
    time_limit: float | None = None,
    n_threads: int | None = None,
    presolve: str = "choose",
    verbose: bool = False,
) -> Problem

Solve on HiGHS (built into landschaft), to a relative gap (0.1 = 10%).

add_gurobi_solver

add_gurobi_solver(
    gap: float = 0.1,
    time_limit: float | None = None,
    n_threads: int | None = None,
    verbose: bool = False,
) -> Problem

Solve on Gurobi (needs the gurobipy package and a Gurobi licence).

add_anneal_solver

add_anneal_solver(
    num_reps: int = 10,
    num_iters: int | None = None,
    start_temp: float = -1.0,
    cool_fac: float = 0.0,
    num_temp: int = 10000,
    spf=(1.0, 2.0),
    itimp_type: int = 0,
    time_limit: float | None = None,
    seed: int = 1,
    n_threads: int | None = None,
) -> Problem

Solve on Marxan's simulated annealing (marxan()'s engine), for minimum-set problems that are hard for a MILP solver: boundary, connectivity and linear penalties, targets, locked units and other ">=" rows on the planning units, binary decisions.

The problem is Marxan's exactly (prioritizr's boundary penalty is Marxan's BLM, its edge factor Marxan's fixed boundaries); every run of annealing and iterative improvement is then repaired to meet every target and trimmed of the units the targets do not need, and the solution is the run with the lowest objective. Its bound is the linear relaxation without the boundary terms (a lower bound on the optimum) and status is "optimal" only when the gap closes, "heuristic" otherwise ("time_limit" if the annealing was cut short).

Parameters:

Name	Type	Description	Default
num_reps	int	Runs (in parallel), from different random reserves.	10
num_iters	int | None	Annealing iterations per run (default: 1000 per planning unit, at least a million); with a time_limit, the runs take as many as fit (the schedule spread over their share of it).	None
start_temp	float	Marxan's schedule: start_temp < 0 for Connolly's (the default) or the adaptive (-2) one, else the initial temperature and its cooling factor; num_temp temperature decreases.	-1.0
cool_fac	float	Marxan's schedule: start_temp < 0 for Connolly's (the default) or the adaptive (-2) one, else the initial temperature and its cooling factor; num_temp temperature decreases.	-1.0
num_temp	float	Marxan's schedule: start_temp < 0 for Connolly's (the default) or the adaptive (-2) one, else the initial temperature and its cooling factor; num_temp temperature decreases.	-1.0
spf		Every feature's penalty factor during the annealing (Marxan's SPF): a shortfall costs spf x the cost of meeting the target alone. Low values let the annealing roam and the repair meets the targets anyway; several values: the runs take them in turn (run r the r % len(spf)-th), the best run wins.	(1.0, 2.0)
itimp_type	int	Marxan's iterative improvement after the annealing (0 normal, 3 twice).	0
time_limit	float | None	Seconds for the runs: every run's annealing gets its share (the runs in waves of n_threads) and cools over it, then improvement and repair.	None
seed	int	The same seed gives the same solution on any number of threads.	1
n_threads	int | None	Runs in parallel (default: all cores but one).	None

compile

compile() -> OptimizationProblem

The model prioritizr's compile() gives for this problem.

solve

solve(start=None)

Solve the problem on its solver (HiGHS unless another was added): a Solution.

Parameters:

Name	Type	Description	Default
start		A solution to start from: one value per planning unit, or a Zonation ranking to build one from: "zonation" (Zonation's CAZ2 on the problem's features, with the cost raster if the units are cells), a ZonationResult or a rank raster. The units go in from the highest rank down until every target is met (locked-in units always, locked-out never); the solver then drops those not needed and goes on from there.	None

Solution dataclass

Solution(
    selection: object,
    values: ndarray,
    coverage: DataFrame,
    objective: float = nan,
    bound: float = nan,
    gap: float = nan,
    status: str = "",
    solver: str = "",
    seconds: float = 0.0,
    zone_names: tuple = (),
    zones: object = None,
    _cost: ndarray = None,
)

Bases: Saveable

A solution of a Problem.

Attributes:

Name	Type	Description
selection	object	The selected share of every planning unit (1 / 0 with binary decisions): a raster on the planning units' grid (NaN outside the planning units) or a copy of their table with a solution column. With zones: a raster with a layer per zone (dimension zone; NaN where a unit has no cost in the zone), or a column solution_<zone> per zone and a zone column with the name of the unit's zone.
values	ndarray	The same, one value per planning unit (in their order); with zones, planning units x zones.
zones	object	With zones: the zone of every planning unit as a raster of zone numbers (1, 2, ... in the order of zone_names; 0: none; NaN outside the planning units; the names in attrs["zones"]), or a column of zone names. With shares (proportion decisions), the zone with the largest share.
coverage	DataFrame	Every target's feature, sense, the feature's total amount, the target, the amount held and the shortfall, absolute and relative, and whether it is met (prioritizr's eval_target_coverage_summary()); without targets (max utility, max cover), every feature's amount held.
objective, bound, gap		The objective of the solution, the solver's bound on the optimum and their relative gap.
status	str	"optimal" (within the gap asked for), "time_limit" (the gap asked for not reached within the time limit: with a limit, models too wide for HiGHS's MIP to keep it stop at the engine's own solution, before the limit), "infeasible", ...
solver, seconds		The solver used and the time it took.

cost property

cost: float

Total cost of the selection (over every zone).

zone_costs property

zone_costs: dict

With zones: the cost of the selection in every zone.

Portfolio dataclass

Portfolio(solutions: list, selection: object, frequency: object, summary: DataFrame)

Bases: Saveable

Several solutions of a problem (add_cuts_portfolio() and the others).

Attributes:

Name	Type	Description
solutions	list	The solutions, best first (cuts, top, gap) or in the order found (shuffle).
selection	object	Every solution's selection: a raster with one band per solution, or the planning units' table with columns solution_1, solution_2...
frequency	object	The share of the solutions that select every planning unit (a raster or a table column, as selection): a simple irreplaceability.
summary	DataFrame	Objective, cost, status and gap of every solution.

OptimizationProblem dataclass

OptimizationProblem(
    A: csr_matrix,
    obj: ndarray,
    lb: ndarray,
    ub: ndarray,
    rhs: ndarray,
    sense: tuple[str, ...],
    vtype: tuple[str, ...],
    modelsense: str,
    row_ids: tuple[str, ...],
    col_ids: tuple[str, ...],
    connected: tuple = (),
    corridors: tuple = (),
)

A compiled model: modelsense of obj @ x subject to A x (sense) rhs, lb <= x <= ub and x of type vtype ("B" binary, "C" continuous, "I" integer); the layout of prioritizr's compile().

ncell property

ncell: int

Non-zeros of the constraint matrix.

row_bounds

row_bounds() -> tuple[ndarray, ndarray]

Lower and upper bounds of every row (-inf / inf when open).

h3_units

h3_units(cost: DataArray, resolution: int, n_threads: int | None = None)

H3 hexagons as planning units: the hexagons of an H3 resolution that cover a cost raster, with their cost.

Pass them to problem() as polygons (cost_column="cost"): the features' amounts are then summed over every hexagon like any polygon's.

Parameters:

Name	Type	Description	Default
cost	DataArray	The cost raster (with a CRS). Every hexagon that covers a cell with a cost is a planning unit; its cost is the sum of the costs of the cells it covers, each weighted by the exact share of the cell inside it.	required
resolution	int	The H3 resolution (0-15; 7 is about 5 km², 8 about 0.7 km²).	required
n_threads	int | None	Threads for the sums (default: all cores but one).	None

Returns:

Type	Description
GeoDataFrame	One hexagon per row (in the raster's CRS, its six vertices reprojected): h3 (the cell's index) and cost.

Examples:

>>> import landschaft as ls
>>> hexes = ls.prioritize.h3_units(cost, 8)
>>> p = ls.prioritize.problem(hexes, features, cost_column="cost")

example_data

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

Path of an example data set for prioritization: Serra Verde, an imaginary coastal region of 60 x 90 cells of 1 km (SIRGAS 2000 / UTM 23S). Mountains rise in the north-west; the coastal plain has three cities, a highway and two rivers; the ocean in the south-east has no data, which leaves 4825 planning units on land. The data are synthetic, made to behave like real planning data.

Parameters:

Name	Type	Description	Default
name	('cost', 'features', 'protected', 'urban', 'condition', 'regions', 'carbon', 'terminals', 'threats', 'elevation', 'lines', 'land_cover', 'past_cover', 'footprint', 'ecoregions', 'carbon_potential')	"cost": land cost (thousand USD per km²), cheap in the mountains, dear near cities and roads; "features": habitat suitability (0-1) of eight features, a stack with the bands jaguar (widespread), lowland_tapir (river valleys), harpy_eagle (remote mountain forest), maned_wolf (western grasslands), golden_lion_tamarin (two small lowland fragments: rare), neotropical_otter (the rivers), araucaria_forest (cold highlands) and mangrove (a thin coastal strip: rare); "protected": three existing protected areas (1); "urban": the cities (1); "condition": habitat condition (0 degraded, 1 pristine); "regions": three administrative regions (1-3); "carbon": carbon stock (t/ha); "terminals": the two mountain parks to link by a corridor (1); "threats": deforestation and fire intensity (0-1); "elevation": relative elevation (0-1); "lines": the rivers and the roads (1), for maps. For restoration plans (restoration_problem()): "land_cover": the share (0-1) of six classes in every cell, the bands forest, savanna, mangrove (natural), pasture, cropland (anthropic) and urban; "past_cover": the share of the three natural classes in the past, before most of the conversion; "footprint": the human footprint index (0-50); "ecoregions": four ecoregions (1 coast, 2 lowland forest, 3 montane forest, 4 savanna); "carbon_potential": the carbon restored vegetation would hold after 30 years (t/ha).	"cost"

Returns:

Type	Description
Path

Examples:

>>> import landschaft as ls
>>> features = ls.load_raster(ls.prioritize.example_data("features"))
>>> features.shape
(8, 60, 90)
>>> [str(b) for b in features.band.values][:3]
['jaguar', 'lowland_tapir', 'harpy_eagle']

Management zones

With zones (prioritizr's zones), every planning unit goes to at most one of several zones, each with its own cost and its own contribution to the features: give problem the costs of the zones (a mapping zone -> raster, or cost columns) and the contributions, or every zone's own features. Targets are summed over the zones unless a zone is given. The solution's zones is the zone of every unit; these two helpers convert between such maps and stacks of zone layers.

category_layer

category_layer(stack) -> DataArray

The zone of every cell of a stack of zone layers (prioritizr's category_layer()): the number (1, 2, ...) of the layer with the largest value above 0, 0 where none is, NaN where every layer is NaN. stack: a raster with a layer per zone (a solution's selection) or a mapping zone -> raster.

Examples:

>>> import landschaft as ls
>>> cost = ls.load_raster(ls.prioritize.example_data("cost"))
>>> features = ls.load_raster(ls.prioritize.example_data("features"))
>>> p = (ls.prioritize.problem({"protect": cost, "use": cost * 0.4}, features, contribution={"use": 0.5})
...        .add_min_set_objective().add_relative_targets(0.3))
>>> sol = p.solve()
>>> zones = ls.prioritize.category_layer(sol.selection)
>>> zones.attrs["zones"]
{1: 'protect', 2: 'use'}
>>> bool((zones == sol.zones).sum() == (sol.zones >= 0).sum())
True

binary_stack

binary_stack(categories: DataArray, zones=None) -> dict

One 0 / 1 raster per zone from a raster of zone numbers (prioritizr's binary_stack()): a mapping zone -> raster, ready for add_locked_in_constraints() (for example an existing zoning plan). zones: the zones' names in the order of their numbers (default: those of a solution's zones map, else "1", "2", ...); cells with 0 or NaN are in no zone.

Examples:

>>> import landschaft as ls
>>> cost = ls.load_raster(ls.prioritize.example_data("cost"))
>>> features = ls.load_raster(ls.prioritize.example_data("features"))
>>> protected = ls.load_raster(ls.prioritize.example_data("protected"))
>>> p = (ls.prioritize.problem({"protect": cost, "use": cost * 0.4}, features, contribution={"use": 0.5})
...        .add_min_set_objective().add_relative_targets(0.3))
>>> existing = protected.fillna(0)              # the protected areas: zone 1 (protect)
>>> locked = ls.prioritize.binary_stack(existing, zones=p.zones)
>>> sorted(locked)
['protect', 'use']
>>> sol = p.add_locked_in_constraints(locked).solve()
>>> bool((sol.zones.values[existing.values == 1] == 1).all())
True

Spatial relations

The matrices that penalties and constraints read: the length of the boundary units share, which units are adjacent, the connectivity between them from a raster or from distances.

boundary_matrix

boundary_matrix(x) -> csc_matrix

The length of the boundary every two units share, and every unit's perimeter on the diagonal (prioritizr's boundary_matrix()): raster cells with data share xres with the cells above and below and yres with those left and right, perimeter 2 (xres + yres); polygons share the length of their common paths, perimeters their own, rounded to 8 decimals.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> bm = P.boundary_matrix(cost)
>>> bm.shape
(5400, 5400)
>>> float(bm[1000, 1000])
4000.0

adjacency_matrix

adjacency_matrix(x, directions: int = 4) -> csc_matrix

1 for every two adjacent units (prioritizr's adjacency_matrix()): raster cells with data that share an edge (directions 4) or also a corner (8); polygons that intersect, touching at a corner included.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> adj = P.adjacency_matrix(cost)
>>> int(adj.sum(axis=1).max())
4

connectivity_matrix

connectivity_matrix(x, y) -> csc_matrix

Connectivity between adjacent units (prioritizr's connectivity_matrix()): the length of their shared boundary times the mean of y in the two (y: a raster on the units' grid; with polygons a raster summed over each, or a column name).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> cm = P.connectivity_matrix(cost, features[0])
>>> cm.shape
(5400, 5400)

proximity_matrix

proximity_matrix(x, distance: float) -> csc_matrix

1 for every two units within distance of each other (prioritizr's proximity_matrix()): raster cell centres, or polygons by their nearest points.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> pm = P.proximity_matrix(cost, 1500)
>>> pm.shape
(5400, 5400)
>>> int(pm.getrow(2000).nnz)
8

decay_connectivity

decay_connectivity(
    distances, median: float, min_probability: float = 0.001
) -> csc_matrix

Connectivity from distances between units as a dispersal probability (Conefor's negative exponential): exp(-ln(2) d / median), 0.5 at the median dispersal distance, dropped below min_probability (sparse), 0 on the diagonal and for unreachable pairs (inf, NaN).

distances: units x units (every unit), e.g. least-cost distances from cost_distance() between unit_points(), commute distances, effective resistances (circuit) or Euclidean ones.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> from scipy.spatial.distance import cdist
>>> pts = P.unit_points(cost)
>>> dc = P.decay_connectivity(cdist(pts, pts), median=3000)
>>> dc.shape
(5400, 5400)
>>> round(float(dc[0, 1]), 3)
0.794

unit_points

unit_points(x) -> ndarray

A point per unit, (n, 2) coordinates: raster cell centres (every cell), or a point inside every polygon (shapely's point on surface), to measure distances between units with cost_distance() or commute_distance().

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> pts = P.unit_points(cost)
>>> pts.shape
(5400, 2)

Evaluating solutions

The figures of a solution (prioritizr's eval_*, equal bit for bit) and the importance of every selected unit; the landscape the solution makes (patches, MSPA) and its species-area extinction risk. They take a Solution, its selection or plain values per planning unit.

eval_n_summary

eval_n_summary(p, solution) -> DataFrame

The number of selected planning units (the sum of the solution's values).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> P.eval_n_summary(p, sol)
   summary      n
0  overall  534.0

eval_cost_summary

eval_cost_summary(p, solution) -> DataFrame

The solution's total cost.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> round(P.eval_cost_summary(p, sol)["cost"].item(), 1) == round(sol.cost, 1)
True

eval_boundary_summary

eval_boundary_summary(
    p, solution, edge_factor: float = 0.5, data=None, zones=None
) -> DataFrame

The solution's boundary length: the edges between selected and other units, plus edge_factor x the selected units' edges on the outside of the planning units.

data: the boundary matrix (every unit's shared lengths and its perimeter on the diagonal; default boundary_matrix() of the planning units). With zones: edge_factor one or one per zone and zones the zones matrix (as add_boundary_penalties()); the rows are the boundary overall (with the zones matrix) and of every zone's units alone.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> P.eval_boundary_summary(p, sol)
   summary  boundary
0  overall  382500.0

eval_connectivity_summary

eval_connectivity_summary(p, solution, data, zones=None) -> DataFrame

The connectivity among the selected units: the sum of data[i, j] x the solution's values of units i and j (a symmetric matrix over the units). With zones, zones the zones matrix (default the identity): the rows are the connectivity overall and within every zone.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> cm = P.connectivity_matrix(cost, features[0])
>>> round(P.eval_connectivity_summary(p, sol, cm)["connectivity"].item(), 2)
1276348.0

eval_asym_connectivity_summary

eval_asym_connectivity_summary(p, solution, data, zones=None) -> DataFrame

As eval_connectivity_summary() for an asymmetric matrix (its diagonal counted once).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> import scipy.sparse as sp
>>> cm = P.connectivity_matrix(cost, features[0])
>>> asym = sp.triu(cm).tocsr()
>>> round(P.eval_asym_connectivity_summary(p, sol, asym)["asym_connectivity"].item(), 2)
638174.0

eval_feature_representation_summary

eval_feature_representation_summary(p, solution) -> DataFrame

Every feature's total amount and the amount the solution holds, absolute and relative.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> rep = P.eval_feature_representation_summary(p, sol)
>>> rep[["feature", "relative_held"]].round(3)
               feature  relative_held
0               jaguar          0.200
1        lowland_tapir          0.201
2          harpy_eagle          0.299
3           maned_wolf          0.205
4  golden_lion_tamarin          0.204
5    neotropical_otter          0.201
6     araucaria_forest          0.428
7             mangrove          0.210

eval_target_coverage_summary

eval_target_coverage_summary(
    p, solution, include_sense: bool = False, include_zone: bool | None = None
) -> DataFrame

Every target's feature, whether it is met, the feature's total, and the target, the amount held and the shortfall, absolute and relative (shares of the total; 0 when the total is 0).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> cov = P.eval_target_coverage_summary(p, sol)
>>> bool(cov["met"].all())
True

eval_replacement_importance

eval_replacement_importance(
    p, solution, rescale: bool = True, n_threads: int | None = None, on=None
)

The replacement cost of every selected unit (Cabeza & Moilanen 2006): how much the objective worsens when the unit is left out and the problem solved again with the other selected units kept (infinite when nothing can replace it), rescaled to 0.01-1. One solve per selected unit, in parallel (n_threads, default all cores but one). on=ls.helpers() solves them on the helpers of this computer's group too (parallel).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> rc = P.eval_replacement_importance(p, sol)
>>> float(rc.max())
1.0

eval_ferrier_importance

eval_ferrier_importance(p, solution)

Ferrier et al. (2000)'s irreplaceability of the selected units for every feature and in total (prioritizr's version: the problem needs one ">=" target per feature, amounts >= 0).

Returns a raster with a band per feature and total (or the units' table with those columns): 0 for the units not selected.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> fi = P.eval_ferrier_importance(p, sol)
>>> [str(b) for b in fi.band.values][-2:]
['mangrove', 'total']
>>> float(fi.sel(band="total").max()) > 0
True

eval_rare_richness_importance

eval_rare_richness_importance(p, solution, rescale: bool = True)

The rarity-weighted richness of the selected units (Albuquerque & Beier 2015): the sum over the features of the unit's amount / the feature's largest amount / its total, rescaled to 0.01-1.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> rwr = P.eval_rare_richness_importance(p, sol)
>>> float(rwr.max())
1.0

eval_rank_importance

eval_rank_importance(
    p,
    solution,
    n: int | None = None,
    budgets=None,
    rescale: bool = True,
    by_zone: bool = True,
)

Rank importance (Jung et al. 2021, prioritizr's eval_rank_importance): the problem is solved with increasing budgets (the min set objective becomes min shortfall), only the solution's units allowed and each step keeping the units of the step before; a unit's rank is how early it enters (n = the first budget, ..., 1 = the last), rescaled to 0.01-1. n budgets evenly spaced between the locked-in units' cost and the solution's, or the budgets themselves.

Returns the per-unit ranks (a raster or the units' table) and a table of every budget's objective, status and gap. With zones, a layer per zone, and by_zone: a budget per zone (every zone's cost in equal steps), else one budget over every zone; budgets may then be steps x zones.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> ranks, steps = P.eval_rank_importance(p, sol, n=4)
>>> steps[["budget", "status"]].round(1)
    budget   status
0   5898.5  optimal
1  11797.0  optimal
2  17695.4  optimal
3  23593.9  optimal

eval_extinction_risk

eval_extinction_risk(p, solution, z=0.25) -> DataFrame

Every feature's extinction risk by the species-area relationship, 1 - (held / total)^z (z = 0.25 by default, one value or one per feature), and their sum (the expected number of features lost).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> risk = P.eval_extinction_risk(p, sol)
>>> round(float(risk["extinction_risk"].iloc[-1]), 3)
2.423

eval_landscape_metrics

eval_landscape_metrics(
    solution, level: str | None = "class", metric=None, zone=None, **kw
) -> DataFrame

Landscape metrics of a raster solution (landscape_metrics() of the map of the selected cells: class 1 selected, 0 not, outside the planning units no data): e.g. the number, sizes and shapes of the selected patches. Class-level and patch-level rows of the selected class only.

With zones, the map of the zones (a class per zone, 0 for the units in none; class and patch rows of the zones, with a zone column), or with zone that zone's units against the others.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> P.eval_landscape_metrics(sol, metric=["np", "ca"])[["metric", "value"]]
  metric    value
0     ca  53400.0
1     np     13.0

eval_mspa

eval_mspa(solution, zone=None, **kw)

MSPA (mspa()) of a raster solution: the selected cells as foreground, the other planning units as background, the cells outside them missing. Returns (the MSPA map, its class areas). With zones, give the zone whose units are the foreground.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> mspa, areas = P.eval_mspa(sol)
>>> int(areas.loc["Core", "pixels"]) >= 0
True

Areas of a priority surface

The separate areas of the top fraction of any priority surface (a Zonation ranking, an importance, a selection frequency) and what each holds, and the overlap of the top fractions of two surfaces (Zonation's landscape identification and comparison; on a Zonation ranking, the same areas as ZonationResult.area_statistics_separate_areas).

eval_landscape_identification

eval_landscape_identification(
    surface,
    top_fraction: float = 0.1,
    *,
    min_size: int = 1,
    merge_distance: int = 0,
    connectivity: int = 8,
    problem=None,
    features=None,
    cost=None,
) -> LandscapeAreas

Landscape identification: the separate areas of the top fraction of a priority surface, and what each holds (Zonation's area statistics of separate areas, over any surface).

Parameters:

Name	Type	Description	Default
surface		The priority surface: a raster (a ranking, an importance, a selection frequency...), a ZonationResult (its rank), a Portfolio (its selection frequency), a Solution (its selection) or a 2-D array; a raster of several bands needs one chosen, except an eval_ferrier_importance() result (its total).	required
top_fraction	float	The share of the surface's cells taken from the top, by Zonation's rule (see Notes).	0.1
min_size	int	Areas of fewer cells are dropped.	1
merge_distance	int	Areas this many cells apart (dilations of the top cells) are merged into one; the area keeps only its own top cells (prioripy's merge_distance).	0
connectivity	int	8 (Zonation's) or 4.	8
problem		A problem on raster planning units of the surface's grid: the features' amounts and the cost in every area, the shares of the features' totals (those relative targets use).	None
features		Without a problem: feature rasters (bands) on the surface's grid and a cost raster.	None
cost		Without a problem: feature rasters (bands) on the surface's grid and a cost raster.	None

Returns:

Type	Description
LandscapeAreas

Notes

The cells are ranked by position as Zonation does (position_rank(); ties: the earlier cell in row order ranks lower) and the top fraction f is the cells ranked at or above 1 - f: on 600 cells, f = 0.1 takes 59 cells (rank 540/600 is just below 0.9 in single precision), as area_statistics_rank_range() does. Areas are numbered by their first cell in row order; the numbers of dropped areas are left out (as Zonation's).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> ferrier = P.eval_ferrier_importance(p, p.add_highs_solver(gap=0.01).solve())
>>> lsi = P.eval_landscape_identification(ferrier, 0.1, min_size=5, problem=p)
>>> len(lsi.table) > 0, bool((lsi.table["cells"] >= 5).all())
(True, True)

eval_landscape_comparison

eval_landscape_comparison(
    a, b, fractions=(0.05, 0.1, 0.2, 0.3, 0.5)
) -> LandscapeComparison

Landscape comparison: how much the top fractions of two priority surfaces overlap (Zonation's solution comparison at several fractions, prioripy's eval_landscape_comparison).

Parameters:

Name	Type	Description	Default
a		Two surfaces on the same grid, of any kind eval_landscape_identification() takes. Cells outside either surface are left out of both.	required
b		Two surfaces on the same grid, of any kind eval_landscape_identification() takes. Cells outside either surface are left out of both.	required
fractions		The top fractions compared (Zonation's rule: see eval_landscape_identification()).	(0.05, 0.1, 0.2, 0.3, 0.5)

Returns:

Type	Description
LandscapeComparison

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> features = ls.load_raster(P.example_data("features"))
>>> cost = ls.load_raster(P.example_data("cost"))
>>> caz = P.zonation(features, rule="CAZ1")
>>> abf = P.zonation(features, rule="ABF")
>>> cmp = P.eval_landscape_comparison(caz, abf, fractions=[0.1, 0.3])
>>> cmp.table[["fraction", "both"]].shape
(2, 2)

position_rank

position_rank(values) -> ndarray

Zonation's ranks of a surface by position: the i-th lowest of the N finite cells gets float32(i * (1 / N)); ties go to the earlier cell in row order first; NaN stays NaN.

On a Zonation rank map this returns the rank map itself, bit for bit.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> ls.prioritize.position_rank(np.array([[3.0, 1.0], [np.nan, 1.0]])).tolist()
[[0.6666666865348816, 0.0], [nan, 0.3333333432674408]]

LandscapeAreas dataclass

LandscapeAreas(areas: object, table: DataFrame, top: object, top_fraction: float = 0.1)

Bases: Saveable

The areas of the top fraction of a priority surface (eval_landscape_identification()).

Attributes:

Name	Type	Description
areas	object	The area of every cell (0: not in an area, NaN outside the surface): a raster like the surface (an array for an array).
table	DataFrame	Per area: area (its number), cells, mean_value (the surface's mean), cost (with a cost) and, per feature, the amount in the area (<feature>) and its share of the feature's total (<feature>_share).
top	object	The cells of the top fraction (1, 0, NaN outside), before areas under min_size were dropped.

LandscapeComparison dataclass

LandscapeComparison(table: DataFrame, maps: object)

Bases: Saveable

Two priority surfaces compared at several top fractions (eval_landscape_comparison()).

Attributes:

Name	Type	Description
table	DataFrame	Per fraction: the cells in the top of a, of b, of both, of only one and of neither, the Jaccard index (both / either) and the share of a's top in b's.
maps	object	Per fraction (band), Zonation's solution comparison codes: 1 in both tops, 2 only in a's, 3 only in b's, 4 in neither, NaN outside.

Irreplaceability (ProxIRR)

ProxIRR (Baisero et al. 2022): for every target, the share of the feature's spare amount a unit holds (alpha), and per unit the chance it is needed for at least one feature (beta) and the number of features it is essential for plus beta (gamma). The bits are those of the CRAN package proxirr.

eval_proxirr

eval_proxirr(p, triage: bool = False, n_threads: int | None = None) -> ProxIRR

ProxIRR (Baisero et al. 2022) of every planning unit for the problem's targets, as the CRAN package proxirr computes it.

For every >= target of the problem, local is the amount in the unit, global the total its relative targets are shares of (the whole raster or table, the region for sub-regional targets) and target its amount; alpha, beta and gamma follow (see the module's notes).

Parameters:

Name	Type	Description	Default
p		A Problem with targets (relative, absolute, manual, log-linear or sub-regional); its objective, constraints and penalties play no part.	required
triage	bool	proxirr's triage: a feature whose target exceeds what is left counts for nothing (alpha 0) instead of making every unit holding it essential.	False
n_threads	int | None	Threads (default: all cores but one).	None

Returns:

Type	Description
ProxIRR

Notes

A feature absent from the landscape (global amount 0) gives alpha 0 (proxirr stops with an error). Targets of sense <= or = are left out.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_relative_targets(0.3)
>>> irr = P.eval_proxirr(p)
>>> float(irr.beta.max()) <= 1.0, float(irr.gamma.max()) >= float(irr.beta.max())
(True, True)

eval_proxirr_marxan

eval_proxirr_marxan(
    pu: DataFrame,
    spec: DataFrame,
    puvspr: DataFrame,
    triage: bool = False,
    n_threads: int | None = None,
) -> ProxIRR

ProxIRR of a Marxan problem's planning units (proxirr's marxan_alphas, marxan_betas and marxan_gammas), from the tables marxan() takes.

Parameters:

Name	Type	Description	Default
pu	DataFrame	Marxan's pu.dat (id), spec.dat (id and prop (a share of the feature's total, used when the column is there) or target) and puvspr.dat (species, pu, amount).	required
spec	DataFrame	Marxan's pu.dat (id), spec.dat (id and prop (a share of the feature's total, used when the column is there) or target) and puvspr.dat (species, pu, amount).	required
puvspr	DataFrame	Marxan's pu.dat (id), spec.dat (id and prop (a share of the feature's total, used when the column is there) or target) and puvspr.dat (species, pu, amount).	required
triage	bool	As eval_proxirr().	False
n_threads	bool	As eval_proxirr().	False

Returns:

Type	Description
ProxIRR	alpha: puvspr with global, target and alpha; beta / gamma: pu with the column; targets: spec with global and target.

Notes

proxirr 0.5 looks a feature's global amount and target up by its position in spec rather than by its id, so with ids other than 1, 2, ... in order every feature gets another's; and a feature of puvspr missing from spec stops marxan_betas with an error. Here features are looked up by id, and those without a target (NaN alphas) are left out of beta and gamma.

Examples:

>>> import pandas as pd
>>> import landschaft as ls
>>> pu = pd.DataFrame({"id": [1, 2, 3], "cost": [1.0, 1.0, 1.0]})
>>> spec = pd.DataFrame({"id": [10, 20], "target": [4.0, 1.0]})
>>> puvspr = pd.DataFrame({"species": [10, 10, 20, 20], "pu": [1, 2, 2, 3], "amount": [3.0, 2.0, 1.0, 1.0]})
>>> irr = ls.prioritize.eval_proxirr_marxan(pu, spec, puvspr)
>>> irr.alpha["alpha"].tolist()
[1.0, 1.0, 1.0, 1.0]
>>> irr.gamma["gamma"].tolist()
[1.0, 2.0, 1.0]

ProxIRR dataclass

ProxIRR(beta: object, gamma: object, alpha: object, targets: DataFrame)

Bases: Saveable

ProxIRR of every planning unit (eval_proxirr(), eval_proxirr_marxan()).

Attributes:

Name	Type	Description
beta, gamma		Per planning unit, like a solution's selection: a raster (NaN outside the units) or the units' table with a beta / gamma column (from Marxan tables: the pu table with the column, as proxirr's marxan_betas and marxan_gammas).
alpha	object	Per (target, unit): a sparse matrix, targets x planning units (from Marxan tables: the puvspr table with global, target and alpha columns, as marxan_alphas).
targets	DataFrame	Every target's feature, global amount and amount (from Marxan tables: spec with global and target).

Calibration and trade-offs

Candidate solutions over a trade-off (a penalty's weight, cost thresholds, budgets, any parameter), the choice of one by TOPSIS or the knee of the curve, and solving again as a problem changes.

calibrate_penalty

calibrate_penalty(
    p,
    penalty: str = "boundary",
    values=None,
    n: int = 9,
    data=None,
    edge_factor: float = 0.5,
    n_threads: int | None = None,
    on=None,
) -> TradeOffs

The blended approach: the problem solved with a range of penalty weights (values; default n values evenly spaced in log scale over 1/100 to 100 times the Cohon et al. penalty), every candidate's cost and criterion (boundary length, connectivity, ...), and the candidates TOPSIS and the knee of the curve choose. p: the problem without the penalty. Solves in parallel. on=ls.helpers() solves them on the helpers of this computer's group too (parallel).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> c, f = P.coarsen(cost, features, 6)
>>> p = P.problem(c, f).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> cal = P.calibrate_penalty(p.add_highs_solver(gap=0.05), "boundary", values=[0.001, 0.01, 0.1])
>>> cal.table[["value", "status"]]
   value   status
0  0.001  optimal
1  0.010  optimal
2  0.100  optimal
>>> best = cal.best("topsis")

calibrate_budget

calibrate_budget(
    p,
    penalty: str = "boundary",
    thresholds=None,
    inflation=None,
    data=None,
    edge_factor: float = 0.5,
    n_threads: int | None = None,
    on=None,
) -> TradeOffs

The hierarchical approach: the least cost solution first (p, a min set problem without the penalty, solved to optimality), then the least penalty within cost thresholds above its cost (thresholds, or the cost x (1 + inflation); default inflation 1e-5 to 4 in 9 steps, the thresholds rounded up as in prioritizr's tutorial). Solves in parallel.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> c, f = P.coarsen(cost, features, 6)
>>> p = P.problem(c, f).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> cal = P.calibrate_budget(p.add_highs_solver(gap=0.05), "boundary", inflation=[0.1, 0.5, 1.0])
>>> len(cal.table)
3
>>> bool((cal.table["cost"] <= cal.table["value"] + 1e-6).all())
True

cohon_penalty

cohon_penalty(
    p, penalty: str = "boundary", data=None, edge_factor: float = 0.5, gap: float = 0.0
)

The penalty of Cohon et al. (1979): the trade-off rate between the solution with the least cost (p as given, a min set problem without the penalty) and the one with the least penalty (min penalties within a budget above every unit's cost): |cost difference / criterion difference|, rounded to 5 decimals. Returns (penalty, a table of the two ideal solutions' cost and criterion).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> c, f = P.coarsen(cost, features, 6)
>>> p = P.problem(c, f).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> value, ideals = P.cohon_penalty(p.add_highs_solver(gap=0.05))
>>> value > 0
True
>>> list(ideals.columns)
['solution', 'cost', 'boundary']

pareto_front

pareto_front(
    p, budgets=None, n: int = 10, n_threads: int | None = None, on=None
) -> TradeOffs

The trade-off between cost and representation: the least total shortfall of the targets (min shortfall objective: the sum of the shortfalls as shares of their targets) within n budgets evenly spaced up to the least-cost solution's cost (above the cost of the locked-in units, which every budget includes), or the budgets given. p: a problem with targets (its objective is replaced). The table adds met (targets met) and efficient (no other candidate is cheaper and not worse, or better and not dearer).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> front = P.pareto_front(p, n=4)
>>> front.table[["met", "efficient"]]
   met  efficient
0    1       True
1    3       True
2    2       True
3    8       True

calibrate

calibrate(
    make,
    values,
    criterion,
    name: str = "criterion",
    n_threads: int | None = None,
    on=None,
) -> TradeOffs

Any trade-off parameter: make(value) gives the problem for every value (a feature's weight, a target, a penalty of any kind...), criterion(problem, solution) the quantity traded off against the cost (smaller is better). The candidates' table as calibrate_penalty()'s (TOPSIS, knee), solved in parallel. on=ls.helpers() solves them on the helpers of this computer's group too (parallel).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> def make(w):
...     return P.problem(cost, features).add_max_utility_objective(5000).add_feature_weights([w, 1, 1, 1, 1, 1, 1, 1])
>>> def lost(q, s):
...     rep = P.eval_feature_representation_summary(q, s)
...     return 1 - rep["relative_held"].iloc[1]
>>> cal = P.calibrate(make, [1, 10], lost, name="tapir_lost")
>>> list(cal.table.columns[:3])
['value', 'cost', 'tapir_lost']

TradeOffs dataclass

TradeOffs(table: DataFrame, solutions: list, criterion: str = '')

Bases: Saveable

Candidate solutions over a trade-off parameter.

Attributes:

Name	Type	Description
table	DataFrame	One row per candidate: value (the penalty, the cost threshold or the budget), cost, the criterion traded off against it (boundary, connectivity, shortfall...), the objective, status and gap; topsis and knee mark the candidates those methods choose.
solutions	list	The candidates' Solution\ s, in the table's order.

best

best(method: str = 'topsis')

The candidate the method chooses ("topsis" with equal weights, or "knee").

topsis

topsis(decision, weights=None, impacts=None) -> DataFrame

TOPSIS (as R's topsis package): the alternatives (rows of decision) scored by their closeness to the ideal (impacts: "+" larger is better, "-" smaller is better, per column; equal weights by default). Returns score (0-1, larger is better) and rank (1 = best).

Examples:

>>> import landschaft as ls
>>> ls.prioritize.topsis([[1.0, 9.0], [4.0, 4.0], [9.0, 1.0]], impacts=["-", "-"])["rank"].tolist()
[2.5, 1.0, 2.5]

knee_point

knee_point(x, y) -> int

The knee (elbow) of a trade-off curve: the point farthest from the straight line between its ends, both axes scaled to 0-1 (the points in the order of x). Returns its position.

Examples:

>>> import landschaft as ls
>>> ls.prioritize.knee_point([0, 1, 2, 3, 4, 5], [10, 4, 2, 1.5, 1.2, 1.0])
2

mean_sd

mean_sd(decision, impacts=None) -> DataFrame

plangea's choice of weights (post_choose_weights): every criterion (a column of decision) rescaled to 0-1 over the alternatives (a constant one to 1), the "-" ones (smaller is better) inverted, and every alternative scored by the mean of its criteria divided by their standard deviation (good on all of them at once). Returns mean, sd, score (larger is better) and rank (1 = best). plangea took the mean after adding the deviation as a column (the mean of the criteria and the deviation); here it is the criteria's.

Examples:

>>> import landschaft as ls
>>> ls.prioritize.mean_sd([[1.0, 9.0], [4.0, 6.0], [9.0, 1.0]], impacts=["+", "+"])["rank"].tolist()
[2.5, 1.0, 2.5]

iterative_solve

iterative_solve(p, update, steps: int, lock_previous: bool = False)

Solve, change the problem, solve again: update(step, solution, problem) gives the next problem (for instance from the landscape a ls.simulate step leaves, with the solution's cells protected by protected_cells() in the simulation's rules), or None to stop. With lock_previous, every solution's units stay selected in the next problem (protected areas that accumulate). Returns the solutions and a table (step, cost, units, status).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> def update(step, sol, q):
...     return q.add_relative_targets(0.2 + 0.05 * (step + 1))
>>> sols, table = P.iterative_solve(p, update, 3, lock_previous=True)
>>> table["step"].tolist()
[0, 1, 2]
>>> bool(table["cost"].is_monotonic_increasing)
True

protected_cells

protected_cells(solution)

The selected cells of a raster solution as a boolean raster (for ls.simulate's restrictions).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> mask = P.protected_cells(sol)
>>> mask.dtype.name
'bool'

Ranking and optimization together

Zonation's feature transforms as the features of a problem, the two paradigms compared at the same cost, features from species distribution models, a report of a solution, and problems solved on coarser blocks.

zonation_features

zonation_features(
    features,
    *,
    names=None,
    condition=None,
    renormalize=False,
    connectivity=None,
    self_multiply=False,
    interaction=None,
    interaction2=None,
    matrix_connectivity=None,
    matrices: dict | None = None,
    retention=None,
    impact_groups=None,
    impact=None,
    cumulative_impact=None,
    bits: int = 16,
    nodata=None,
    n_threads: int | None = None,
)

The features as Zonation 5 transforms them before ranking (condition, connectivity kernels, interactions, matrix connectivity, retention, impacts: see zonation() for each), back in the features' own units, for an optimization problem (problem(cost, zonation_features(...))).

The values carry Zonation's quantization (bits 16: a relative error of about 1e-5 of the feature's largest value). Returns a stack (band = the features' names), NaN outside the data.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> features = ls.load_raster(P.example_data("features"))
>>> condition = ls.load_raster(P.example_data("condition"))
>>> z = P.zonation_features(features, condition=[condition] * 8)
>>> z.shape
(8, 60, 90)

compare_zonation

compare_zonation(
    p,
    solution=None,
    zonation=None,
    rule: str = "CAZ2",
    n_threads: int | None = None,
    **zonation_kw,
) -> ZonationComparison

Compare an optimization solution of p (solved if not given) with a Zonation ranking of the same features and costs (run with rule and zonation_kw if not given: a ZonationResult or a rank raster) cut at the solution's cost. Raster planning units.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> c = P.compare_zonation(p, sol)
>>> c.summary["method"].tolist()
['milp', 'zonation']
>>> bool(c.table["milp_met"].all())
True

ZonationComparison dataclass

ZonationComparison(
    table: DataFrame,
    summary: DataFrame,
    selection: DataArray,
    solution: object,
    zonation: object,
)

Bases: Saveable

A Zonation ranking and an optimization solution with the same budget.

Attributes:

Name	Type	Description
table	DataFrame	Every feature's share held by the solution (milp) and by Zonation's best cells up to the same cost (zonation), with its target and whether each meets it (with targets).
summary	DataFrame	Cost, units and boundary of each.
selection	DataArray	The two selections (band milp and zonation: 1 selected, 0 not, NaN outside the units).

sdm_features

sdm_features(results, env, *, threshold=None, ensemble: str = 'weighted') -> DataArray

Features from species distribution models (ls.sdm): every result's prediction on env (its ensemble when it has several models, else its model), as a stack for problem().

Parameters:

Name	Type	Description	Default
results		{name: SDMResult} or a list of them (named species_1, ...).	required
env		The predictors (as for predict()).	required
threshold		None: the suitability (0-1) is the amount; a number, or "max_tss" (every model's maximum-TSS threshold, their mean for an ensemble): presence (1) or absence (0).	None

Examples:

>>> import landschaft as ls
>>> from sklearn.linear_model import LogisticRegression
>>> occ = ls.load_table(ls.sdm.example_data("bradypus"))
>>> env = ls.load_raster(ls.sdm.example_data("bioclim"))
>>> data = ls.sdm.prepare(occ, env, clean=False, background=500, seed=1)
>>> res = ls.sdm.fit(data, {"glm": LogisticRegression(max_iter=1000)})
>>> f = ls.prioritize.sdm_features({"bradypus": res}, env, threshold="max_tss")
>>> [str(b) for b in f.band.values]
['bradypus']

report

report(
    p,
    solution=None,
    *,
    importance: bool = True,
    zonation: bool = False,
    extinction_z=None,
    title: str = "Prioritization",
    n_threads: int | None = None,
) -> PrioritizationReport

Solve p (unless a solution is given) and evaluate the solution: the number of units, cost, boundary (raster and polygon units), the features' representation and the targets' coverage; with importance the Ferrier irreplaceability of the selected units (problems with one target per feature) or else their rarity-weighted richness; with extinction_z the species-area extinction risk; with zonation the comparison with a Zonation ranking at the same cost (raster units). With zones, every summary has a row overall and one per zone (no irreplaceability, rarity or Zonation comparison: they need a single zone).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> sol = p.solve()
>>> r = P.report(p, sol)
>>> sorted(r.tables)
['representation', 'summary', 'targets']
>>> page = ls.plot(r, show=False).save("report.html")

PrioritizationReport dataclass

PrioritizationReport(
    solution: object,
    tables: dict,
    maps: dict,
    comparison: ZonationComparison | None = None,
    title: str = "Prioritization",
)

Bases: Saveable

A solution and its evaluation in one object: ls.plot(report) draws it, ls.plot(report).save( "report.html") writes it as a page that works offline, ls.save(report, folder) its tables and maps.

coarsen

coarsen(cost: DataArray, features, factor: int)

A problem's rasters on blocks of factor x factor cells (for continental problems: solve the coarse problem, then refine() the selection): every block's cost and amounts are the sums of its cells' with a cost (the fine planning units; a block without any has none), so the refined selection holds what the coarse one does. Returns (cost, features).

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> c, f = P.coarsen(cost, features, 4)
>>> c.shape, f.shape
((15, 23), (8, 15, 23))

refine

refine(selection, like: DataArray) -> DataArray

A selection on coarse blocks (a Solution of a coarsen()-ed problem, or its selection raster) back on the fine grid of like (the fine cost raster): every cell takes its block's value, NaN where like has none.

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> c, f = P.coarsen(cost, features, 4)
>>> coarse = P.problem(c, f).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions().solve()
>>> fine = P.refine(coarse, cost)
>>> fine.shape
(60, 90)

Multi-action planning

Which units to monitor and which actions to take against which threats in them, so that the features recover (or are kept) to their targets at the least cost, or recover the most within a budget: prioriactions's models (Salgado-Rojas et al. 2020), bit for bit, from tables or rasters.

actions_data

actions_data(
    pu: DataFrame,
    features: DataFrame,
    dist_features: DataFrame,
    threats: DataFrame,
    dist_threats: DataFrame,
    sensitivity: DataFrame | None = None,
    boundary: DataFrame | None = None,
) -> ActionsData

The inputs of a multi-action problem, checked and completed as prioriactions's inputData().

Parameters:

Name	Type	Description	Default
pu	DataFrame	Planning units: id, monitoring_cost and optionally status (0 free, 2 locked in, 3 locked out).	required
features	DataFrame	id, target_recovery, optionally target_conservation (default 0) and name.	required
dist_features	DataFrame	The features' amounts in the units: pu, feature, amount.	required
threats	DataFrame	id, optionally name and blm_actions (the weight of the boundary of the actions against it; default 0).	required
dist_threats	DataFrame	The threats in the units and the actions against them: pu, threat, amount (the threat's intensity), action_cost and optionally status (0, 2 the action locked in, 3 locked out).	required
sensitivity	DataFrame | None	Which features each threat harms: feature, threat and optionally the response's delta1 to delta4 (default: every feature to every threat, deltas 0, the threat's highest amount, 0 and 1).	None
boundary	DataFrame | None	id1, id2, boundary (shared lengths).	None

Notes

Rows with amount 0 are dropped; features without units, threats without units and threats that harm no feature are dropped with a warning; numbers are rounded to 3 decimals (R's round). With no status in dist_threats the actions are free (prioriactions resets the units' status there instead), and the actions in locked-out units are locked out.

Examples:

>>> import pandas as pd
>>> import landschaft as ls
>>> pu = pd.DataFrame({"id": [1, 2, 3, 4], "monitoring_cost": [2.0, 1.0, 3.0, 1.0]})
>>> features = pd.DataFrame({"id": [1, 2], "target_recovery": [1.0, 0.5]})
>>> dist_features = pd.DataFrame({"pu": [1, 2, 3, 3, 4], "feature": [1, 1, 1, 2, 2], "amount": 1.0})
>>> threats = pd.DataFrame({"id": [1]})
>>> dist_threats = pd.DataFrame({"pu": [1, 2, 3, 4], "threat": 1, "amount": 1.0, "action_cost": [1.0, 4.0, 2.0, 1.0]})
>>> d = ls.prioritize.actions_data(pu, features, dist_features, threats, dist_threats)
>>> sol = d.problem("minimizeCosts").solve()
>>> sol.objective
5.0
>>> sol.action_table()
   pu  threat  solution
0   1       1       1.0
1   2       1       0.0
2   3       1       0.0
3   4       1       1.0

actions_raster_data

actions_raster_data(
    monitoring_cost: DataArray,
    features,
    threats,
    action_cost=1.0,
    *,
    target_recovery=0.0,
    target_conservation=0.0,
    relative: bool = True,
    sensitivity: DataFrame | None = None,
    blm_actions=0.0,
    locked_in=None,
    locked_out=None,
    boundary: bool = False,
    n_threads: int | None = None,
) -> ActionsData

The inputs of a multi-action problem from rasters (see actions_data() for the tables).

Parameters:

Name	Type	Description	Default
monitoring_cost	DataArray	The cost of monitoring every cell: the planning units are the cells with a cost (not NaN), in the order of the cells (ids 1, 2, ...).	required
features		The features' amounts: a stack (band, y, x), a Dataset or a list of rasters on the same grid (the layers' names are the features' names; NaN = 0).	required
threats		The threats' intensity, the same way: a threat occurs in a cell with a value > 0, and its action there removes it.	required
action_cost		The cost of acting against every threat: a number, a raster, or one of them per threat (a sequence in the threats' order, a mapping by name or a stack); it needs a value wherever the threat occurs.	1.0
target_recovery		The features' targets: a number, a sequence in the features' order or a mapping by name. With relative (the default) shares of the amount the feature has where a threat harms it (recovery) or where none does (conservation); otherwise absolute amounts of benefit.	0.0
target_conservation		The features' targets: a number, a sequence in the features' order or a mapping by name. With relative (the default) shares of the amount the feature has where a threat harms it (recovery) or where none does (conservation); otherwise absolute amounts of benefit.	0.0
sensitivity	DataFrame | None	Which features each threat harms (feature, threat by name or 1-based position, and optionally delta1 to delta4); default every feature to every threat.	None
blm_actions		The weight of the boundary of every threat's actions (a number, sequence or mapping).	0.0
locked_in		Rasters of the cells that must (or must not) be monitored (non-zero, not NaN).	None
locked_out		Rasters of the cells that must (or must not) be monitored (non-zero, not NaN).	None
boundary	bool	Give the problem the cells' shared edges as boundary data (for blm and blm_actions).	False

Returns:

Type	Description
ActionsData	With grid set: the solution's maps().

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> threats = ls.load_raster(P.example_data("threats"))
>>> d = P.actions_raster_data(cost, features, threats, action_cost=[20, 5], target_recovery=0.2)
>>> d.n_units, d.n_actions
(4825, 2838)
>>> sol = d.problem("minimizeCosts").add_highs_solver(gap=0.01).solve()
>>> sol.status
'optimal'
>>> [str(b) for b in sol.maps().band.values]
['monitoring', 'deforestation', 'fire']

ActionsData dataclass

ActionsData(
    pu: DataFrame,
    features: DataFrame,
    dist_features: DataFrame,
    threats: DataFrame,
    dist_threats: DataFrame,
    sensitivity: DataFrame,
    boundary: DataFrame | None = None,
    grid: object = None,
)

Checked inputs of a multi-action problem (prioriactions's inputData()): see actions_data().

The tables keep their columns plus internal_id (pu, features, threats: positions from 0) and internal_pu, internal_feature, internal_threat, internal_id1, internal_id2.

problem

problem(
    model: str = "minimizeCosts",
    budget: float = 0.0,
    blm: float = 0.0,
    curve: int = 1,
    segments: int = 3,
) -> ActionsProblem

The optimization problem (prioriactions's problem()).

Parameters:

Name	Type	Description	Default
model	str	"minimizeCosts" ("min_costs"): meet the features' recovery and conservation targets at the least cost of monitoring and actions (plus blm x the boundary); "maximizeBenefits" ("max_benefits"): the most benefit within budget.	'minimizeCosts'
budget	float	The most the monitoring and the actions may cost (maximizeBenefits).	0.0
blm	float	The weight of the selected units' boundary (needs boundary data); every threat's blm_actions weighs the boundary of its actions.	0.0
curve	int	The features' response to the actions: 1 linear, 2 quadratic, 3 cubic (benefit = b^curve, needs threat amounts of 0 or 1).	1
segments	int	Pieces of the curve's linear approximation (1 to 5): pieces of length round(1 / segments, 1).	3

potential_benefit

potential_benefit() -> DataFrame

Every feature's distribution (units, threatened units) and the most benefit the actions can give it: conservation (amount in the units without its threats), recovery and total (prioriactions's getPotentialBenefit(), the conservation benefit as the amount).

ActionsProblem dataclass

ActionsProblem(
    data: ActionsData,
    model: str = "minimizeCosts",
    budget: float = 0.0,
    blm: float = 0.0,
    curve: int = 1,
    segments: int = 3,
    solver: tuple = ("highs", {"gap": 0.0}),
    dependencies: tuple = (),
    portfolio: tuple | None = None,
)

A multi-action problem (prioriactions's problem()): see problem().

compile

compile() -> OptimizationProblem

prioriactions's model (getModelInfo()'s matrices), bit for bit; with curve 2 or 3 the power constraints are not in it (Gurobi's general constraints in prioriactions).

add_directional_dependencies

add_directional_dependencies(
    flow: DataFrame, penalty: float = 1.0, hard: bool = False
) -> ActionsProblem

Dependencies between planning units, such as upstream and downstream along a river (an extension of prioriactions, after prioripy): monitoring a target_pu without its source_pu adds strength x penalty to the cost (or takes it from the benefit), or, with hard, is not allowed.

Parameters:

Name	Type	Description	Default
flow	DataFrame	source_pu, target_pu (the units' ids) and optionally strength (default 1).	required
penalty	float	Multiplies every dependency's strength.	1.0
hard	bool	Make the dependencies constraints: a target is monitored only with its source.	False

add_cuts_portfolio

add_cuts_portfolio(number_solutions: int = 10) -> ActionsProblem

The best number_solutions distinct solutions (differing in a unit monitored or an action taken), one solve after another, each forbidding the earlier ones.

add_top_portfolio

add_top_portfolio(number_solutions: int = 10) -> ActionsProblem

As add_cuts_portfolio(), every solve to optimality (gap 0).

add_gap_portfolio

add_gap_portfolio(number_solutions: int = 10, pool_gap: float = 0.1) -> ActionsProblem

Up to number_solutions distinct solutions within pool_gap (relative) of the optimum.

add_shuffle_portfolio

add_shuffle_portfolio(
    number_solutions: int = 10, remove_duplicates: bool = True, seed: int | None = None
) -> ActionsProblem

number_solutions solves of the model with its units and actions in random orders.

add_highs_solver

add_highs_solver(
    gap: float = 0.0,
    time_limit: float | None = None,
    n_threads: int | None = None,
    verbose: bool = False,
    engine: bool = True,
) -> ActionsProblem

Solve on HiGHS to a relative gap (prioriactions's gap_limit, default 0), behind landschaft's closure engine (engine=False: HiGHS's MIP alone; see solve()).

add_gurobi_solver

add_gurobi_solver(
    gap: float = 0.0,
    time_limit: float | None = None,
    n_threads: int | None = None,
    verbose: bool = False,
) -> ActionsProblem

Solve on Gurobi (needs gurobipy and a licence).

solve

solve(start=None)

Solve the problem.

On HiGHS (the default), the equivalent model (activation one row per action, the benefit columns substituted) first goes to landschaft's closure engine (a Lagrangian bound as strong as the linear relaxation, by a closed form or minimum cuts, and repaired selections), then, while the gap is wider than asked, to HiGHS's MIP from its solution with the units and actions the bound proves fixed. start: a solution to start from (an ActionsSolution or the units then the actions, 0/1). With a portfolio, an ActionsPortfolio (solved on HiGHS's MIP or Gurobi: the no-good cuts do not fit the engine).

ActionsSolution dataclass

ActionsSolution(
    data: ActionsData,
    units: ndarray,
    actions: ndarray,
    local_benefit: DataFrame,
    objective: float = nan,
    bound: float = nan,
    gap: float = nan,
    status: str = "",
    solver: str = "",
    seconds: float = 0.0,
    blm: float = 0.0,
    x: ndarray = None,
    dependencies: tuple = (),
)

Bases: Saveable

A solution of an ActionsProblem.

Attributes:

Name	Type	Description
units, actions		The monitored units (0/1, in the order of data.pu) and the actions taken (0/1, in the order of data.dist_threats): prioriactions's sol_monitoring and sol_actions.
local_benefit	DataFrame	Every (unit, feature)'s benefit, by conservation or recovery (getSolutionBenefit(type = "local")).
objective, bound, gap, status, solver, seconds		As Solution's.

maps

maps()

The monitored cells and every threat's actions as rasters (band monitoring, then the threats' names; NaN outside the planning units), for problems from actions_raster_data().

dependency_penalty

dependency_penalty() -> float

The directional dependencies' penalty: the weights of the targets monitored without their sources (0 without soft dependencies).

monitoring

monitoring() -> DataFrame

id and solution of every planning unit.

action_table

action_table() -> DataFrame

pu, threat and solution of every action.

cost

cost() -> DataFrame

The monitoring cost and every threat's actions' cost (getCost()).

connectivity_penalty

connectivity_penalty() -> DataFrame | None

The boundary of the monitored units and of every threat's actions (getConnectivityPenalty(); None without boundary data).

action_summary

action_summary(format: str = 'wide') -> DataFrame

Every unit's actions (by threat id), whether it conserves a feature without actions (conservation) and whether it is monitored only for connectivity (connectivity): prioriactions's getActions(). format="long": one row per (pu, action).

benefit

benefit() -> DataFrame

Every feature's benefit, by conservation and recovery (getSolutionBenefit(); the conservation benefit as the amount in the monitored units).

ActionsPortfolio dataclass

ActionsPortfolio(solutions: list, summary: DataFrame, frequency: dict)

Bases: Saveable

Several solutions of an ActionsProblem (add_cuts_portfolio() and the others).

Attributes:

Name	Type	Description
solutions	list	The ActionsSolution\ s, best first (cuts, top, gap) or in the order found (shuffle).
summary	DataFrame	Objective, cost (monitoring and actions), status and gap of every solution.
frequency	dict	The share of the solutions that monitor every unit and take every action: monitoring (id, frequency) and actions (pu, threat, frequency).

maps

maps()

The frequencies as rasters (band monitoring, then the threats), for problems from actions_raster_data().

Restoration planning

How much of every converted cell to restore, for the most reduction of the species' extinction risk (and of the ecosystems' risk of collapse), carbon and ecological intactness within area targets and a budget, or the cheapest plan for a benefit: the restoration problem of plangea (Strassburg et al. 2020), with its conservation (conservable=) and conversion (convertible=) forms. The species-area curves are concave, so the optimum is found and proven, by simplicial decomposition (linear programs over the cells, mixed by a small concave problem), from thousands of cells to millions; plangea's own procedure (benefit per cost per cell, the target in steps) is method="steps". Targets per region come from projections of land demand (add_projected_change), shares of the land (relative=True), other values (add_linear_target) or the ecosystems' risk of collapse (ecosystem_targets); tradeoff gives the plans of a growing target and the gradient of priorities, compare_weights ranks sets of weights. The inputs come from ls.original_vegetation, ls.sdm.area_of_habitat and ls.intactness. Tutorial: restoration priorities.

restoration_problem

restoration_problem(
    restorable=None,
    *,
    conservable=None,
    convertible=None,
    limit=1.0,
    threat=1.0,
    cell_area=None,
) -> RestorationProblem

A restoration problem: how much of every converted cell to restore (and of natural land to conserve or convert).

The planning units are the cells with land an action can take; the decision is the share x_i of each, from 0 to the share that can be times limit. Then add benefits (add_extinction_risk_benefit(), add_linear_benefit(), add_intactness_benefit()), targets (add_area_target()), a budget, and solve.

Parameters:

Name	Type	Description	Default
restorable		The share of every cell that can be restored: a raster, or a OriginalVegetation (its anthropic share).	None
conservable		The share of every cell's natural land that conservation can keep (plangea's C): without the plan, threat of it is lost (plangea: all of it). The benefits count what conserving saves.	None
convertible		The share of every cell's natural land that can be converted (plangea's D, a demand for land met where it costs the least benefit). With restorable too, the plan restores some cells and converts others (plangea's R + D). Not with conservable.	None
limit		plangea's bound limit: a number or a raster multiplying the shares, or a mapping action -> limit.	1.0
threat		The share of the conservable land lost without conservation, in [0, 1]: a number or a raster (default 1, plangea's).	1.0
cell_area		The cells' area in km²: a number, a raster, or None (from the grid, see cell_area_km2()).	None

Returns:

Type	Description
RestorationProblem

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> rng = np.random.default_rng(1)
>>> anthropic = rng.random((20, 30)) * 0.6
>>> carbon = rng.random((20, 30)) * 100
>>> p = (ls.prioritize.restoration_problem(anthropic, cell_area=1.0)
...        .add_linear_benefit(carbon, name="carbon")
...        .add_area_target(50.0))
>>> sol = p.solve()
>>> round(sol.area, 6), sol.status
(50.0, 'optimal')

RestorationProblem dataclass

RestorationProblem(
    template: DataArray | None,
    shape: tuple,
    index: ndarray,
    upper: ndarray,
    area: ndarray,
    benefits: tuple = (),
    targets: tuple = (),
    cost: ndarray | None = None,
    cost_layer: ndarray | None = None,
    budget: float | None = None,
    objective: tuple = ("max_benefit",),
    solver: tuple = ("highs", {}),
    action: ndarray | None = None,
    natural: ndarray | None = None,
    threat: ndarray | None = None,
    lower: ndarray | None = None,
)

A restoration problem (see restoration_problem()), built up by its add_* methods; every method returns a new problem.

n property

n: int

The number of decision columns: the planning units of every action.

actions property

actions: tuple

The actions of the plan, of "restore", "conserve" and "convert".

add_extinction_risk_benefit

add_extinction_risk_benefit(
    areas,
    z=PLANGEA_Z,
    weight: float = 1.0,
    restoration_discount=0.0,
    name: str | None = None,
) -> RestorationProblem

The reduction of the species' extinction risk (or of ecosystems' risk of collapse).

A species with habitat A of a potential P has the risk 1 - (A / P)^z (0 at A >= P); restoring a share x_i of a cell gives it g_i x_i more habitat, g the gain of area_of_habitat() (or of ecosystem_areas()). Conserving a share of a cell's natural land keeps the habitat that would be lost (its loss times the threat); converting it takes that habitat away. The benefit is the sum over the species of the risk without the plan minus the risk with it (negative when conversion raises the risk).

Parameters:

Name	Type	Description	Default
areas		A HabitatAreas with its gain (area_of_habitat given the original vegetation) for restoration and its loss (loss=True) for conservation and conversion, on the problem's grid; the species that are not valid are left out.	required
z		The exponent of the species-area relationship, in (0, 1): a number for every species (default 0.25, plangea's), one per species (in the order of areas.table) or a mapping species -> z (species left out of it get 0.25).	PLANGEA_Z
weight	float	The benefit's weight in the objective.	1.0
restoration_discount		The share of the habitat value that restored land does NOT deliver, in [0, 1] (0: restored land is habitat in full, as Strassburg et al. 2020 assume): the gain becomes g (1 - discount). A number, one per species (in the order of areas.table) or a mapping species -> discount. It is not a financial discount rate, and only restoration has it.	0.0
name	str | None	The benefit's name (default "extinction_risk", numbered when repeated).	None

add_linear_benefit

add_linear_benefit(
    values, weight: float = 1.0, name: str | None = None, per: str = "area"
) -> RestorationProblem

A benefit proportional to the area acted on, such as carbon.

values: a raster on the problem's grid, a number or one value per unit, or a mapping action -> values (0 for the actions left out); per km² (per="area") or per cell in full (per="cell") of the action. A loss is negative (the carbon stock of converted land, say); conservation's values count the land the plan saves from loss (they are multiplied by the threat).

add_intactness_benefit

add_intactness_benefit(
    hfi,
    reduction=None,
    radius: int = 3,
    weight: float = 1.0,
    *,
    increase=None,
    gamma: float | None = None,
    beta: float = -0.2,
    z: float = 0.5,
    cell_size=None,
    form: str = "linear",
    name: str | None = None,
) -> RestorationProblem

Ecological intactness (Beyer et al. 2020): restoring a share x of a cell lowers its human footprint by reduction x, converting it raises the footprint by increase x (and natural land lost for want of conservation by increase times the threat, which conserving spares); the benefit is the total intactness gained. hfi, reduction and increase on the problem's grid (or numbers); the other parameters as intactness().

form="linear": linearised at the plan of nothing (intactness_gradient(); with method="steps" it is recomputed at every step, as plangea does). form="exact": the total intactness itself, which is not concave (a sum of exponentials of the footprint): the exact method improves the plan by successive concave minorants (the tangent of every pair's exponential, which never overestimates; each step solved exactly, the benefit never falls) to a local optimum, and bounds the best plan by the secant of every pair over the footprints the plans can reach. The status is "optimal" when the two meet within the tolerance, else "local" (the bound is reported). The weight must be at least 0.

add_area_target

add_area_target(
    target,
    sense: str = "<=",
    regions=None,
    action: str | None = None,
    relative: bool = False,
) -> RestorationProblem

The area acted on (km²): in the whole landscape, or in every region.

Parameters:

Name	Type	Description	Default
target		A number (km²), or with regions a mapping region -> target, or a table with columns region, target and optionally sense.	required
sense	str	"<=" (at most: the most valuable cells for that area), ">=" or "=".	'<='
regions		A raster of region codes on the problem's grid (NaN: no region) or one code per unit.	None
action	str | None	The area of "restore", "conserve" or "convert", of "all" of them, or the "net" change of natural land (restored minus converted). Default: "net" when the plan restores and converts, "all" otherwise.	None
relative	bool	The targets are shares (0-1) of the area that can be acted on (in the landscape or in every region): plangea's flat targets.	False

add_linear_target

add_linear_target(
    values,
    target,
    sense: str = "<=",
    regions=None,
    per: str = "area",
    relative: bool = False,
) -> RestorationProblem

A target on a value other than the area (production, opportunity cost, carbon ...): the sum of values over the plan, in the landscape or in every region.

values: a raster, a number, one value per unit or a mapping action -> values; per km² (per="area") or per cell in full (per="cell") of the action. target, sense, regions and relative as add_area_target() (plangea's targets of other types).

add_projected_change

add_projected_change(
    change,
    regions,
    sense: str = "=",
    restrict: bool = True,
    compensate: bool | None = None,
) -> RestorationProblem

Projected changes of natural land per region (plangea's econometric targets): positive, land spared for restoration; negative, a demand for land to convert.

Parameters:

Name	Type	Description	Default
change		A mapping region -> change (km²) or a table with columns region, change and optionally sense.	required
regions		A raster of region codes on the problem's grid or one code per unit.	required
sense	str	Of every region's row ("=": the change as projected).	'='
restrict	bool	A plan that restores and converts acts in every region in the projected direction only: restoring where land is spared, converting where it is demanded, neither where the change is 0 (what plangea's harmonize_sr_r_or_c_bvals meant to do). Its rows are the net change.	True
compensate	bool | None	A plan that only restores cannot convert: with compensate (the default, plangea's harmonize_sr_sparable) the land demanded elsewhere is taken from the restoration, every region with spared land restoring its share of the net spared land; without it, the regions restore what they spare. A plan that only converts converts what the regions demand.	None

add_budget

add_budget(budget: float, cost) -> RestorationProblem

Spend at most budget, cost per km² acted on (a raster, a number, one value per unit or a mapping action -> cost; the opportunity cost of the land, say).

add_max_benefit_objective

add_max_benefit_objective() -> RestorationProblem

The most weighted benefit within the targets and the budget (the default).

add_min_cost_objective

add_min_cost_objective(benefit: float, cost=None) -> RestorationProblem

The cheapest plan whose weighted benefit reaches benefit; cost per km² acted on (or the budget's, without it).

add_highs_solver

add_highs_solver(
    tolerance: float = 1e-07,
    time_limit: float | None = None,
    n_threads: int | None = None,
    verbose: bool = False,
) -> RestorationProblem

Solve the linear programs on HiGHS (built in). tolerance: the relative gap between the bound and the plan's benefit (or cost) at which the exact method stops.

add_gurobi_solver

add_gurobi_solver(
    tolerance: float = 1e-07,
    time_limit: float | None = None,
    n_threads: int | None = None,
    verbose: bool = False,
) -> RestorationProblem

Solve the linear programs on Gurobi (needs gurobipy and a licence).

solve

solve(method: str = 'exact', steps: int = 1) -> RestorationSolution

Solve the problem.

method="exact": the optimum within the solver's tolerance (relative, between the plan's benefit and a bound on the best one). Simplicial decomposition: every vertex of the set of plans is the solution of a linear program over the units with the gradient of the benefits as its objective and only the targets and the budget as rows; the plan is the best mix of the vertices found (a small concave problem, by SLSQP); the bound is F(x) + grad F(x) . (v - x), v the next vertex (valid because the benefits are concave). The minimum cost: a bisection on the budget, each step the most benefit for it; the bound is the largest budget proven too small.

method="steps": plangea's procedure, for comparison: every cell scored by its weighted benefits (each layer scaled by 1 / 10^trunc(log10(max))) per cost, the area targets met in steps equal parts, the extinction-risk and intactness layers recomputed after every part. It is a heuristic: its benefit is at most the exact one, and depends on steps. Area targets only (no budget, no minimum cost).

evaluate

evaluate(
    plan=None, *, restored=None, conserved=None, converted=None
) -> RestorationSolution

The benefits, cost and area of a plan: the shares of every column (plan: a raster on the grid or one value per unit when the problem has one action, or one value per column), or the shares restored, conserved and converted (rasters or one value per unit; 0 for those left out).

uniform

uniform(area: float | None = None) -> RestorationSolution

plangea's uniform scenario: every unit acted on by the same share of what it can be, for a total area (default: the first landscape area target). For a plan of one action.

tradeoff

tradeoff(
    levels=(0.25, 0.5, 0.75, 1.0),
    *,
    nested: bool = True,
    targets: str = "landscape",
    method: str = "exact",
    steps: int = 1,
    threshold: float = 0.0,
) -> RestorationTradeOff

The plans for a growing target: the trade-off curve, and the priority of every cell.

Parameters:

Name	Type	Description	Default
levels		Shares of the targets (in (0, 1], increasing): plangea's gradient steps.	(0.25, 0.5, 0.75, 1.0)
nested	bool	Every level's plan contains the previous one (a priority order: the exact best plan given the earlier levels, with their shares as the least shares), or (False) every level is the exact optimum on its own (the frontier: plans need not contain each other).	True
targets	str	The target rows that grow: the landscape's ("landscape", plangea's default; the regions' targets hold at every level) or "all".	'landscape'
method	str	"exact" or "steps" (plangea's procedure, nested by construction: the levels follow one another, steps parts each, the layers recomputed from the whole plan so far).	'exact'
threshold	float	A cell counts in a level when its share is above it (the gradient).	0.0

Returns:

Type	Description
RestorationTradeOff	The table (one row per level: the target, area, cost, benefits), the plans and the gradient (plangea's post_gradient_output): for restoration and conservation, the number of levels whose plan acts on the cell (the most for the first chosen); for conversion, the levels plus one minus that number (the least for the first converted).

compare_weights

compare_weights(
    weights, *, method: str = "exact", steps: int = 1
) -> RestorationTradeOff

Solve with every set of the benefits' weights and rank the sets by plangea's metric (post_choose_weights): every benefit's value (and the cost, inverted) rescaled to 0-1 over the sets, the mean over them divided by their standard deviation (see mean_sd()); the best set does well on every benefit at once.

weights: a list of mappings benefit name -> weight (the others keep theirs). Returns the table (w.<benefit>, the values, cost, mean, sd, score, rank) and the plans; best() is the first in rank.

RestorationSolution dataclass

RestorationSolution(
    restored: DataArray | None,
    values: ndarray,
    area: float = 0.0,
    cost: float = nan,
    benefits: DataFrame = None,
    benefit: float = 0.0,
    objective: float = nan,
    bound: float = nan,
    status: str = "",
    method: str = "",
    rounds: int = 0,
    seconds: float = 0.0,
    conserved: DataArray | None = None,
    converted: DataArray | None = None,
    areas: dict = dict(),
)

Bases: Saveable

A restoration plan.

Attributes:

Name	Type	Description
restored, conserved, converted		The share of every cell restored, conserved and converted (NaN outside the planning units; None for an action the problem does not have).
values	ndarray	The shares, one value per column (the units of every action, block after block: restore, conserve, convert).
area, areas		The area acted on (km²) and the area of every action.
cost	float	Its cost (NaN without a cost).
benefits	DataFrame	One row per benefit: its value, weight and weighted value.
benefit	float	The weighted benefit.
objective, bound		The objective (the weighted benefit, or the cost with the minimum-cost objective) and the bound on the optimum the exact method proved (NaN for method="steps").
status, method, rounds, seconds		"optimal" (within the tolerance), "time_limit", "iteration_limit", "heuristic" (steps) ...; the method, the iterations (vertices or bisection steps; steps for method="steps"), the time.

RestorationTradeOff dataclass

RestorationTradeOff(table: DataFrame, solutions: list, gradients: dict = dict())

Bases: Saveable

Plans over a trade-off (tradeoff()) or over sets of weights (compare_weights()).

Attributes:

Name	Type	Description
table	DataFrame	One row per plan.
solutions	list	The RestorationSolution of every row.
gradients	dict	The trade-off's gradient per action (an empty mapping for weights).

gradient property

gradient: DataArray

The gradient of a plan of one action (see tradeoff()).

best

best() -> RestorationSolution

The plan ranked first (compare_weights), or the knee of the curve of benefit against area.

ecosystem_areas

ecosystem_areas(land_cover, natural, anthropic, ecoregions, *, cell_area=None)

The natural area of every ecosystem (ecoregion) now and before conversion, and what restoring a cell gives it: the input of add_extinction_risk_benefit() for the risk of ecosystem collapse (plangea's ec).

Parameters:

Name	Type	Description	Default
land_cover		The share of every class in each cell ((band, y, x), band names the classes).	required
natural		The natural and the anthropic classes.	required
anthropic		The natural and the anthropic classes.	required
ecoregions		A raster of ecoregion codes (NaN: none).	required
cell_area		As restoration_problem().	None

Returns:

Type	Description
HabitatAreas	One row per ecoregion (species its code): current the natural area, potential the natural plus the anthropic area, the gain (the cell's area per unit restored, in the cells with anthropic land) and the loss (the cell's area per unit of natural land lost, in the cells with natural land).

Notes

plangea sums the natural area of an ecoregion only over the cells of its problem (those with anthropic land), which leaves out the ecoregion's untouched land and overstates its loss; here every cell of the ecoregion counts.

Examples:

>>> import numpy as np
>>> import xarray as xr
>>> import landschaft as ls
>>> cover = xr.DataArray(np.array([[[0.8, 0.2, 0.0, 0.5]], [[0.2, 0.8, 1.0, 0.5]]]),
...                      dims=("band", "y", "x"), coords={"band": ["forest", "crops"]})
>>> eco = np.array([[1.0, 1.0, 2.0, 2.0]])
>>> ea = ls.prioritize.ecosystem_areas(cover, ["forest"], ["crops"], eco, cell_area=1.0)
>>> ea.table[["species", "current", "potential"]]
  species  current  potential
0       1      1.0        2.0
1       2      0.5        2.0

ecosystem_targets

ecosystem_targets(
    areas,
    target: float,
    *,
    minimum: float = 0.0,
    z: float = PLANGEA_Z,
    half_life: float = 0.0,
    sense: str = "<=",
    method: str = "exact",
    steps: int = 10,
) -> DataFrame

Share a restoration target among ecosystems for the least risk of collapse (plangea's ecosystemic sub-region targets, harmonize_sr_ecosys).

Every ecosystem e with natural area N_e of a potential T_e (ecosystem_areas()) gets a_e of the target, between its least share and its converted area T_e - N_e, for the most total of ((N_e + a_e) / T_e)^z, each ecosystem's marginal value damped by exp(half_life (N_e + a_e) / T_e) (plangea's "yukawa" factor; 0: none).

Parameters:

Name	Type	Description	Default
areas		A HabitatAreas of ecosystems (ecosystem_areas()).	required
target	float	The area to share (km²).	required
minimum	float	Every ecosystem's least share, as a part (0-1) of an equal split of the target (plangea's subregion_ecosystems_minimum_par), up to its converted area.	0.0
z	float	The exponent of the species-area relationship.	PLANGEA_Z
half_life	float	The damping exponent (below 1 - z, so the marginal value falls with the area).	0.0
sense	str	The sense of the targets returned (plangea's subregion_ecosystemic_targets_sense).	'<='
method	str	"exact": the optimum (the marginal values equal among the ecosystems not at a bound, by bisection). "steps": plangea's procedure (the target in steps parts, each a linear program on the marginal values, recomputed between them).	'exact'

Returns:

Type	Description
DataFrame	region (the ecosystem's code) and target, sense: a table for add_area_target() with the ecosystem raster as regions.

Examples:

>>> import numpy as np
>>> import xarray as xr
>>> import landschaft as ls
>>> cover = xr.DataArray(np.array([[[0.8, 0.2, 0.1, 0.5]], [[0.2, 0.8, 0.9, 0.5]]]),
...                      dims=("band", "y", "x"), coords={"band": ["forest", "crops"]})
>>> eco = np.array([[1.0, 1.0, 2.0, 2.0]])
>>> ea = ls.prioritize.ecosystem_areas(cover, ["forest"], ["crops"], eco, cell_area=1.0)
>>> ls.prioritize.ecosystem_targets(ea, 1.0).round(4)
   region  target sense
0     1.0     0.3    <=
1     2.0     0.7    <=

eval_risk_reduction

eval_risk_reduction(
    problem: RestorationProblem,
    solution,
    z_sd: float | None = None,
    *,
    benefit: str | None = None,
    method: str = "correlated",
    step: float = 1e-06,
) -> RiskReduction

The reduction of the extinction risk a plan brings (plangea's multi_exrisk_aggregate), with the uncertainty of the species-area exponent z.

Parameters:

Name	Type	Description	Default
problem	RestorationProblem	The restoration problem.	required
solution		A RestorationSolution, or the share restored per unit (or a raster).	required
z_sd	float | None	The standard deviation of z (plangea: 0.1): a number for every species, one per species (in the order of the benefit's species) or a mapping species -> deviation (0 for the others); None: no uncertainty.	None
benefit	str | None	The extinction-risk benefit to evaluate (default: the first).	None
method	str	How the species' sensitivities to z add up. "correlated" (the default): the species' z err together (one z for all of them, or per-species values from one source), so the total's deviation is |sum_s d reduction_s / dz_s sd_s| (the delta method on the total). "independent" sums them in quadrature and "linear" sums their absolute values, plangea's two forms (they treat the species' z as separate parameters).	'correlated'
step	float	The forward difference of z (plangea's 1e-6).	1e-06

Returns:

Type	Description
RiskReduction

Examples:

>>> import numpy as np
>>> import xarray as xr
>>> import landschaft as ls
>>> cover = xr.DataArray(np.array([[[0.6, 0.0, 0.3]], [[0.2, 0.1, 0.0]], [[0.2, 0.9, 0.7]]]),
...                      dims=("band", "y", "x"), coords={"band": ["forest", "savanna", "crops"]})
>>> ranges = xr.DataArray(np.ones((2, 1, 3)), dims=("band", "y", "x"), coords={"band": ["toucan", "rhea"]})
>>> veg = ls.original_vegetation(cover, ["forest", "savanna"], ["crops"])
>>> aoh = ls.sdm.area_of_habitat(ranges, cover, ["forest", "savanna"], original=veg, cell_area=1.0,
...                              suitability={"toucan": ["forest"], "rhea": ["savanna"]})
>>> p = (ls.prioritize.restoration_problem(veg, cell_area=1.0)
...        .add_extinction_risk_benefit(aoh).add_area_target(0.8))
>>> sol = p.solve()
>>> rr = ls.prioritize.eval_risk_reduction(p, sol, z_sd=0.1)
>>> round(rr.total, 4), round(rr.sd, 4)
(0.269, 0.0871)

RiskReduction dataclass

RiskReduction(table: DataFrame, total: float = 0.0, sd: float = nan, method: str = '')

Bases: Saveable

Result of eval_risk_reduction(): per species, its habitat now (current), without the plan (baseline: less than now when conservation's baseline loses land) and gained by it, its risk before and after; the total reduction (negative when conversion raises the risk) and its standard deviation from the uncertainty of z.

Simulated annealing

Marxan 4's runs, unit by unit as the program, and its annealing as a solver of any problem (Problem.add_anneal_solver) for scales where a gap matters less than an answer.

marxan

marxan(
    pu: DataFrame,
    spec: DataFrame,
    puvspr: DataFrame,
    bound: DataFrame | None = None,
    *,
    blm: float = 0.0,
    num_reps: int = 10,
    num_iters: int = 1000000,
    start_temp: float = -1.0,
    cool_fac: float = 0.0,
    num_temp: int = 10000,
    run_mode: int = 1,
    heur_type: int = 1,
    itimp_type: int = 0,
    prop: float = 0.0,
    misslevel: float = 1.0,
    cost_thresh: float = 0.0,
    thresh_pen1: float = 0.0,
    thresh_pen2: float = 0.0,
    connectivity_in: bool = False,
    seed: int = 1,
    n_threads: int | None = None,
) -> MarxanResult

Run Marxan: num_reps independent runs of annealing and iterative improvement (or the greedy heuristics) minimising cost + blm x boundary + penalties for the features' shortfalls.

Parameters:

Name	Type	Description	Default
pu	DataFrame	Planning units: id, cost and optionally status (0 available, 1 in the initial reserve, 2 locked in, 3 locked out). Marxan's pu.dat.	required
spec	DataFrame	Features: id, and target (an amount), or prop (a share of the total), spf (the feature penalty factor, default 1) and optionally targetocc (occurrences). Marxan's spec.dat.	required
puvspr	DataFrame	Amount of every feature in every planning unit: species, pu, amount (one row per pair). Marxan's puvspr.dat.	required
bound	DataFrame | None	Boundaries: id1, id2, boundary; a row with id1 == id2 is a unit's boundary with the outside (its fixed boundary). Marxan's bound.dat.	None
blm	float		0.0
num_reps	float		0.0
num_iters	float		0.0
start_temp	float		0.0
cool_fac	float		0.0
num_temp	float		0.0
run_mode	float		0.0
heur_type	float		0.0
itimp_type	float		0.0
prop	float		0.0
misslevel	float	Marxan's input.dat options (BLM, NUMREPS, NUMITNS, STARTTEMP: < 0 for Connolly's (-1) or the adaptive (-2) schedule, COOLFAC, NUMTEMP, RUNMODE 0, 1, 2 or 6, HEURTYPE 0 or 1, ITIMPTYPE 0 or 3, PROP, MISSLEVEL, COSTTHRESH, THRESHPEN1, THRESHPEN2).	1.0
cost_thresh	float	Marxan's input.dat options (BLM, NUMREPS, NUMITNS, STARTTEMP: < 0 for Connolly's (-1) or the adaptive (-2) schedule, COOLFAC, NUMTEMP, RUNMODE 0, 1, 2 or 6, HEURTYPE 0 or 1, ITIMPTYPE 0 or 3, PROP, MISSLEVEL, COSTTHRESH, THRESHPEN1, THRESHPEN2).	1.0
thresh_pen1	float	Marxan's input.dat options (BLM, NUMREPS, NUMITNS, STARTTEMP: < 0 for Connolly's (-1) or the adaptive (-2) schedule, COOLFAC, NUMTEMP, RUNMODE 0, 1, 2 or 6, HEURTYPE 0 or 1, ITIMPTYPE 0 or 3, PROP, MISSLEVEL, COSTTHRESH, THRESHPEN1, THRESHPEN2).	1.0
thresh_pen2	float	Marxan's input.dat options (BLM, NUMREPS, NUMITNS, STARTTEMP: < 0 for Connolly's (-1) or the adaptive (-2) schedule, COOLFAC, NUMTEMP, RUNMODE 0, 1, 2 or 6, HEURTYPE 0 or 1, ITIMPTYPE 0 or 3, PROP, MISSLEVEL, COSTTHRESH, THRESHPEN1, THRESHPEN2).	1.0
connectivity_in	bool	Marxan's CONNECTIVITYIN: count the boundaries inside the reserve instead of at its edge.	False
seed	int	RANDSEED: the same seed gives the same runs on any number of threads.	1
n_threads	int | None	Runs in parallel (default: all cores but one).	None

Returns:

Type	Description
MarxanResult

Examples:

>>> import landschaft as ls
>>> P = ls.prioritize
>>> cost = ls.load_raster(P.example_data("cost"))
>>> features = ls.load_raster(P.example_data("features"))
>>> p = P.problem(cost, features).add_min_set_objective().add_relative_targets(0.2).add_binary_decisions()
>>> pu = p.units
>>> import numpy as np, pandas as pd
>>> rij = pu.rij.tocoo()
>>> res = P.marxan(pd.DataFrame({"id": np.arange(1, pu.n + 1), "cost": pu.cost, "status": 0}),
...                pd.DataFrame({"id": np.arange(1, 9), "prop": 0.2, "spf": 10.0}),
...                pd.DataFrame({"species": rij.row + 1, "pu": rij.col + 1, "amount": rij.data}),
...                num_reps=4, num_iters=20000, seed=1)
>>> len(res.summary)
4
>>> list(res.solutions.columns[:3])
['planning_unit', 'solution_1', 'solution_2']

MarxanResult dataclass

MarxanResult(
    solutions: DataFrame,
    best: DataFrame,
    summary: DataFrame,
    penalties: DataFrame,
    frequency: DataFrame,
    best_run: int = 1,
    unmet: int = 0,
)

Bases: Saveable

The runs of marxan().

Attributes:

Name	Type	Description
solutions	DataFrame	Planning units (planning_unit: their ids) and every run's selection (solution_1, ...: 1 selected, 0 not; locked units 1 or 0).
best	DataFrame	The best run's selection (planning_unit, solution): the lowest score, the first run on ties.
summary	DataFrame	Every run's figures, as Marxan's _sum file: Run_Number, Score, Cost, Planning_Units, Connectivity (boundary length), Penalty, Shortfall, Missing_Values (features below misslevel of their target) and MPM (the lowest share of a target met).
penalties	DataFrame	Every feature's spid, its target (after the proportion targets) and its penalty (Marxan's greedy estimate of the cost of meeting its target alone).
frequency	DataFrame	planning_unit and frequency: the share of the runs that select each unit (Marxan's summed solution over the runs).
best_run	int	The best run's number (from 1).