Distance & restoration¶

Distances to the habitat edge, resistance maps, accumulated costs and least-cost paths, the connectivity gain of restoring a corridor, the zones of influence of habitat objects, and two inputs of restoration plans: the original vegetation of converted land and ecological intactness. Tutorial: distance and restoration.

Distance to the habitat¶

euclidean_distance ¶

euclidean_distance(
    data: DataArray | ndarray,
    *,
    pixel_size: float | tuple[float, float] | None = None,
    threads: int = 0,
    n_threads: int | None = None,
) -> EuclideanDistanceResult

Euclidean distance to the habitat edge, inside and outside the habitat (GTB's batch_eucldist).

Every foreground pixel gets its straight-line distance to the nearest background pixel (positive), every background pixel its distance to the nearest foreground (negative). The border pixels are replicated so that the map edge does not act as an edge.

Parameters:

Name	Type	Description	Default
`data`	`DataArray \| ndarray`	Map coded 0 missing, 1 background, 2 foreground (DataArray, array or path).	required
`pixel_size`	`float \| tuple[float, float] \| None`	Pixel size (e.g. in metres) for the statistics in map units; read from a georeferenced DataArray when omitted.	`None`
`threads`	`int`	Threads (0 or `None`: all cores but one).	`0`
`n_threads`	`int`	Threads (0 or `None`: all cores but one).	`0`

Returns:

Type	Description
`EuclideanDistanceResult`

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))   # 0 missing, 1 background, 2 forest
>>> distance = ls.euclidean_distance(forest, pixel_size=25)
>>> round(distance.mean_fg_distance, 2), round(distance.max_bg_distance, 2)
(6.01, 33.54)
>>> distance.stats[["class", "pixel_count", "mean_distance_units"]]
                 class  pixel_count  mean_distance_units
0  Foreground (Forest)       428490               150.37
1  Background (Matrix)       571240               125.76

EuclideanDistanceResult `dataclass` ¶

EuclideanDistanceResult(
    signed_distance: DataArray | ndarray,
    byte_map: DataArray | ndarray,
    mean_fg_distance: float,
    max_fg_distance: float,
    mean_bg_distance: float,
    max_bg_distance: float,
    fg_pixel_count: int,
    bg_pixel_count: int,
    stats: DataFrame,
    pixel_size: float | None = None,
)

Bases: Saveable

Result of euclidean_distance().

Attributes:

Name	Type	Description
`signed_distance`	`DataArray \| ndarray`	Distance of every pixel in pixels: positive in the foreground, negative in the background, NaN on missing pixels.
`byte_map`	`DataArray \| ndarray`	GTB's `uint8` display map of the distances (see `dist_colormap()`).
`mean_fg_distance, max_fg_distance`		Mean and maximum distance of the foreground pixels, in pixels.
`mean_bg_distance, max_bg_distance`		The same for the background (absolute values).
`fg_pixel_count, bg_pixel_count`		Number of foreground and background pixels.
`stats`	`DataFrame`	The same statistics as a table, in pixels and in map units.
`pixel_size`	`float \| None`	Pixel size used for the map units (`None`: unknown).

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))   # 0 missing, 1 background, 2 forest
>>> print(ls.euclidean_distance(forest, pixel_size=25).summary())
================================================================
GuidosToolbox (GTB) - Euclidean Distance Transform
...

summary ¶

summary() -> str

Text summary of the distance statistics.

Returns:

Type	Description
`str`

distance_to_resistance ¶

distance_to_resistance(data, *, n_threads: int | None = None)

Resistance growing with the distance to the foreground (GTB Setup Tools -> Distance -> Resistance).

data is coded 0 missing, 1 background, 2 foreground. Background pixels get round(distance + 2) capped at 100 (so 3 next to an object), foreground 2 and missing 0. The distance is IDL's MORPH_DISTANCE as GTB uses it (quasi-Euclidean; see influence_zones() for the remaining differences to IDL).

Parameters:

Name	Type	Description	Default
`data`		Map coded 0 missing, 1 background, 2 foreground.	required
`n_threads`	`int \| None`	Threads (0 or `None`: all cores but one).	`None`

Returns:

Type	Description
`DataArray or ndarray`	`uint8` resistance map (see `resistance_colormap()`).

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))   # 0 missing, 1 background, 2 forest
>>> resistance = ls.distance_to_resistance(forest)
>>> int(resistance.min()), int(resistance.max())
(0, 38)

Restoration status¶

restoration_status ¶

restoration_status(data: DataArray | ndarray | str | Path) -> RestorationStatusResult

How connected the habitat is and how much reconnecting it could add (GTB's Restoration Status, batch_rss).

The habitat objects (8-connected foreground components) of areas A_i give:

ECA (equivalent connected area) = sqrt(sum of A_i²), the size of a single object with the same connectivity;
DOC (degree of coherence) = ECA / total habitat area, in %;
restoration potential = 100 - DOC;
RAC (relative area of habitat) = habitat area / valid area, in %;
APS (average patch size) and CNOA (critical new object area).

Parameters:

Name	Type	Description	Default
`data`	`DataArray \| ndarray \| str \| Path`	Map coded 0 missing, 1 background (or resistance 1-100), 2 foreground (DataArray, array or path).	required

Returns:

Type	Description
`RestorationStatusResult`

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))   # 0 missing, 1 background, 2 forest
>>> status = ls.restoration_status(forest)
>>> status.nr_comp, round(status.doc, 1), round(status.rest_pot, 1)
(2850, 51.6, 48.4)

RestorationStatusResult `dataclass` ¶

RestorationStatusResult(
    data_area: int,
    nr_comp: int,
    eca_max: int,
    larg_obj: int,
    aps: float,
    cnoa: int,
    eca: float,
    doc: float,
    rest_pot: float,
    rac: float,
)

Bases: Saveable

Result of restoration_status().

Attributes:

Name	Type	Description
`data_area`	`int`	Valid pixels (all but the missing ones).
`nr_comp`	`int`	Number of habitat objects (8-connected).
`eca_max`	`int`	Total habitat area, pixels.
`larg_obj`	`int`	Area of the largest object, pixels.
`aps`	`float`	Average patch size, eca_max / nr_comp.
`cnoa`	`int`	Critical new object area (GTB's formula), floor(2 b c² / (b² - c²) + 1) with b = eca_max and c = eca.
`eca`	`float`	Equivalent connected area, sqrt(sum of A_i²), pixels.
`doc`	`float`	Degree of coherence, ECA / eca_max in %.
`rest_pot`	`float`	Restoration potential, 100 - doc.
`rac`	`float`	Relative area of habitat, eca_max / data_area in %.

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))   # 0 missing, 1 background, 2 forest
>>> status = ls.restoration_status(forest)
>>> status.to_dataframe()[["NR_OBJ", "COH[%]", "REST_POT[%]"]]
   NR_OBJ  COH[%]  REST_POT[%]
0    2850   51.64        48.36

to_dict ¶

to_dict() -> dict[str, Any]

The metrics as a dictionary.

Returns:

Type	Description
`dict`

to_dataframe ¶

to_dataframe() -> DataFrame

The metrics as one row of GTB's batch_rss spreadsheet.

Returns:

Type	Description
`DataFrame`

summary ¶

summary() -> str

Text summary of the restoration status.

Returns:

Type	Description
`str`

Costs, paths and corridors¶

cost_distance ¶

cost_distance(
    resistance: DataArray | ndarray | str | Path,
    source_mask: DataArray | ndarray | str | Path,
) -> DataArray | ndarray

Accumulated cost of travelling from the source cells across a resistance surface (8 neighbours).

A step costs the mean resistance of the two cells times the step length (1, or sqrt(2) on the diagonals).

Parameters:

Name	Type	Description	Default
`resistance`	`DataArray \| ndarray \| str \| Path`	Resistance surface (1-100; 0 = impassable / NoData): DataArray, array or path.	required
`source_mask`	`DataArray \| ndarray \| str \| Path`	Cells > 0 (or `True`) are the sources, at cost 0.	required

Returns:

Type	Description
`DataArray or ndarray`	Accumulated cost; `inf` where the sources cannot be reached.

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))   # 0 missing, 1 background, 2 forest
>>> import numpy as np
>>> resistance = ls.distance_to_resistance(forest)
>>> objects = ls.accounting(forest)                     # the two largest forest objects
>>> ids = objects.ids.values
>>> marker = np.zeros(forest.shape, dtype="uint8")
>>> marker[ids == objects.largest[0][0]] = 1            # object A
>>> marker[ids == objects.largest[1][0]] = 2            # object B
>>> marker[forest.values == 0] = 3                      # missing: not crossed
>>> cost = ls.cost_distance(resistance, marker == 1)
>>> round(float(cost.values[marker == 2].min()), 1)      # cheapest cost to reach B
370.4

least_cost_path ¶

least_cost_path(
    resistance: DataArray | ndarray | str | Path,
    marker_a: DataArray | ndarray | str | Path,
    marker_b: DataArray | ndarray | str | Path,
    dist_map: DataArray | ndarray | None = None,
) -> LeastCostPathResult

The least-cost path between two areas A and B.

The costs from A and from B (cost_distance()) are added; the path is the cells of minimal total cost.

Parameters:

Name	Type	Description	Default
`resistance`	`DataArray \| ndarray \| str \| Path`	Resistance surface (1-100; 0 = impassable).	required
`marker_a`	`DataArray \| ndarray \| str \| Path`	Masks of areas A and B.	required
`marker_b`	`DataArray \| ndarray \| str \| Path`	Masks of areas A and B.	required
`dist_map`	`DataArray \| ndarray \| None`	Optional distance map (e.g. `euclidean_distance()`) to report the mean distance of the restored cells.	`None`

Returns:

Type	Description
`LeastCostPathResult`

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))   # 0 missing, 1 background, 2 forest
>>> import numpy as np
>>> resistance = ls.distance_to_resistance(forest)
>>> objects = ls.accounting(forest)                     # the two largest forest objects
>>> ids = objects.ids.values
>>> marker = np.zeros(forest.shape, dtype="uint8")
>>> marker[ids == objects.largest[0][0]] = 1            # object A
>>> marker[ids == objects.largest[1][0]] = 2            # object B
>>> marker[forest.values == 0] = 3                      # missing: not crossed
>>> path = ls.least_cost_path(resistance, marker == 1, marker == 2)
>>> path.path_length, path.rest_pixels, path.expense
(115, 113, 323.0)

LeastCostPathResult `dataclass` ¶

LeastCostPathResult(
    combined_cost: DataArray | ndarray,
    lcp_mask: DataArray | ndarray,
    path_length: int,
    rest_pixels: int,
    expense: float,
    avg_distance: float,
)

Bases: Saveable

Result of least_cost_path().

Attributes:

Name	Type	Description
`combined_cost`	`DataArray or ndarray`	Cost from A plus cost from B.
`lcp_mask`	`DataArray or ndarray`	1 on the path, 0 elsewhere.
`path_length`	`int`	Length of the path, pixels.
`rest_pixels`	`int`	Pixels of the path in the background: what would have to be restored.
`expense`	`float`	Total resistance of those pixels.
`avg_distance`	`float`	Mean distance of those pixels (with `dist_map`; else 0).

restoration_corridor ¶

restoration_corridor(
    data: DataArray | ndarray | str | Path,
    resistance: DataArray | ndarray | str | Path,
    marker_a: DataArray | ndarray,
    marker_b: DataArray | ndarray,
    dist_map: DataArray | ndarray | None = None,
) -> RestorationCorridorResult

Connectivity gain of restoring the least-cost corridor between two areas (GTB's Restoration Planner).

The corridor of least_cost_path() is turned into habitat and the landscape's ECA and degree of coherence (restoration_status()) are compared before and after.

Parameters:

Name	Type	Description	Default
`data`	`DataArray \| ndarray \| str \| Path`	Habitat map coded 0 missing, 1 background, 2 foreground.	required
`resistance`	`DataArray \| ndarray \| str \| Path`	Resistance surface (1-100; 0 = impassable).	required
`marker_a`	`DataArray \| ndarray`	Masks of areas A and B.	required
`marker_b`	`DataArray \| ndarray`	Masks of areas A and B.	required
`dist_map`	`DataArray \| ndarray \| None`	Optional distance map, as in `least_cost_path()`.	`None`

Returns:

Type	Description
`RestorationCorridorResult`

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))   # 0 missing, 1 background, 2 forest
>>> import numpy as np
>>> resistance = ls.distance_to_resistance(forest)
>>> objects = ls.accounting(forest)                     # the two largest forest objects
>>> ids = objects.ids.values
>>> marker = np.zeros(forest.shape, dtype="uint8")
>>> marker[ids == objects.largest[0][0]] = 1            # object A
>>> marker[ids == objects.largest[1][0]] = 2            # object B
>>> marker[forest.values == 0] = 3                      # missing: not crossed
>>> corridor = ls.restoration_corridor(forest, resistance, marker == 1, marker == 2)
>>> corridor.rest_pixels, round(corridor.delta_eca, 2), round(corridor.delta_doc, 2)
(113, 31784.28, 7.41)

RestorationCorridorResult `dataclass` ¶

RestorationCorridorResult(
    lcp_result: LeastCostPathResult,
    size_a: int,
    size_b: int,
    rest_pixels: int,
    expense: float,
    efficiency: float,
    eca_orig: float,
    eca_new: float,
    delta_eca: float,
    doc_orig: float,
    doc_new: float,
    delta_doc: float,
)

Bases: Saveable

Result of restoration_corridor().

Attributes:

Name	Type	Description
`lcp_result`	`LeastCostPathResult`	Underlying least-cost path result.
`size_a`	`int`	Pixel area of source component A.
`size_b`	`int`	Pixel area of target component B.
`rest_pixels`	`int`	Number of restored pixels needed to establish the corridor.
`expense`	`float`	Accumulated resistance cost of restored pixels.
`efficiency`	`float`	Restoration efficiency, delta_eca / expense.
`eca_orig`	`float`	Landscape ECA before restoration.
`eca_new`	`float`	Landscape ECA after establishing the corridor.
`delta_eca`	`float`	Net increase in ECA, eca_new - eca_orig.
`doc_orig`	`float`	Landscape Degree of Coherence (%) before restoration.
`doc_new`	`float`	Landscape Degree of Coherence (%) after restoration.
`delta_doc`	`float`	Net increase in Degree of Coherence (% points).

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))   # 0 missing, 1 background, 2 forest
>>> import numpy as np
>>> resistance = ls.distance_to_resistance(forest)
>>> objects = ls.accounting(forest)                     # the two largest forest objects
>>> ids = objects.ids.values
>>> marker = np.zeros(forest.shape, dtype="uint8")
>>> marker[ids == objects.largest[0][0]] = 1            # object A
>>> marker[ids == objects.largest[1][0]] = 2            # object B
>>> marker[forest.values == 0] = 3                      # missing: not crossed
>>> corridor = ls.restoration_corridor(forest, resistance, marker == 1, marker == 2)
>>> corridor.to_dataframe()[["REST_PIX", "COH_ORIG", "COH_NEW"]]
   REST_PIX  COH_ORIG  COH_NEW
0       113     51.64    59.06

to_dataframe ¶

to_dataframe() -> DataFrame

The metrics as one row of GTB's Restoration Planner spreadsheet.

Returns:

Type	Description
`DataFrame`

isochrone_map ¶

isochrone_map(resistance, marker) -> IsochroneResult

Accumulated cost from object A (GTB Restoration Planner -> Isochrone Map A).

Parameters:

Name	Type	Description	Default
`resistance`		`uint8` resistance map: 2 = foreground (cost 1), 3-100 = background resistance (GTB's Restoration Planner coding; 0 costs nothing, as in GTB).	required
`marker`		`uint8` map of the same size: 0 background, 1 object A (a single 8-connected object), 3 missing (not crossed); 2 is ignored.	required

Returns:

Type	Description
`IsochroneResult`

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))   # 0 missing, 1 background, 2 forest
>>> import numpy as np
>>> resistance = ls.distance_to_resistance(forest)
>>> objects = ls.accounting(forest)                     # the two largest forest objects
>>> ids = objects.ids.values
>>> marker = np.zeros(forest.shape, dtype="uint8")
>>> marker[ids == objects.largest[0][0]] = 1            # object A
>>> marker[ids == objects.largest[1][0]] = 2            # object B
>>> marker[forest.values == 0] = 3                      # missing: not crossed
>>> only_a = np.where(marker == 2, 0, marker).astype("uint8")
>>> rings = ls.isochrone_map(resistance, only_a)
>>> rings.min_cost, rings.max_cost
(2, 1166)

isochrone_map_ab ¶

isochrone_map_ab(
    resistance, marker, *, n_threads: int | None = None
) -> IsochroneResult

Cheapest route between objects A and B (GTB Restoration Planner -> Isochrone Map AB).

marker holds 1 = object A, 2 = object B (each a single 8-connected object, not touching each other), 3 = missing, 0 = background; see isochrone_map() for resistance.

Parameters:

Name	Type	Description	Default
`resistance`		`uint8` resistance map, as in `isochrone_map()`.	required
`marker`		`uint8` map: 1 object A, 2 object B, 3 missing, 0 background.	required
`n_threads`	`int \| None`	Threads (0 or `None`: all cores but one).	`None`

Returns:

Type	Description
`IsochroneResult`	With the least-cost path in `lcp`.

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))   # 0 missing, 1 background, 2 forest
>>> import numpy as np
>>> resistance = ls.distance_to_resistance(forest)
>>> objects = ls.accounting(forest)                     # the two largest forest objects
>>> ids = objects.ids.values
>>> marker = np.zeros(forest.shape, dtype="uint8")
>>> marker[ids == objects.largest[0][0]] = 1            # object A
>>> marker[ids == objects.largest[1][0]] = 2            # object B
>>> marker[forest.values == 0] = 3                      # missing: not crossed
>>> route = ls.isochrone_map_ab(resistance, marker)
>>> route.min_cost, int(np.asarray(route.lcp).sum())
(210, 215)

IsochroneResult `dataclass` ¶

IsochroneResult(
    map: DataArray | ndarray,
    cost: DataArray | ndarray,
    lcp: DataArray | ndarray | None,
    min_cost: int,
    max_cost: int,
)

Bases: Saveable

Result of isochrone_map() / isochrone_map_ab().

Attributes:

Name	Type	Description
`map`	`DataArray \| ndarray`	`uint8` isochrone map in GTB's codes: 0-100 cost scaled between `min_cost` and `max_cost`, 101 markers, 102 missing, 103 least-cost path (AB), 104 unreachable.
`cost`	`DataArray \| ndarray`	Accumulated cost as GTB rounds it (`int64`; unreachable/missing: -2^31). For map AB, the cost of the cheapest A-B route through each pixel.
`lcp`	`DataArray \| ndarray \| None`	Least-cost path mask (map AB), else `None`: the pixels of minimal A-B cost, thinned to a 1-pixel wide path as GTB's `ggeo` does.

Zones of influence¶

influence_zones ¶

influence_zones(
    data,
    min_area: int = 5000,
    *,
    core_buffer: int = 0,
    buffer: int = 0,
    n_threads: int = 0,
) -> InfluenceZonesResult

Zones of influence of the foreground objects (GTB Distance -> Influence Zones).

Parameters:

Name	Type	Description	Default
`data`		2-D `uint8` map (array, DataArray or raster path) coded 0 = missing, 1 = background, 2 = foreground.	required
`min_area`	`int`	Minimum object area [pixels] (GTB default 5000). Smaller objects are not objects; they are shown as code 176.	`5000`
`core_buffer`	`int`	Reduce the objects to their core: drop the foreground within `core_buffer` pixels of the object edge before measuring areas (GTB "CoreZone"; 0 = off).	`0`
`buffer`	`int`	Width [pixels] of a background buffer drawn around the objects (0 = none).	`0`
`n_threads`	`int`	Threads (0 or `None`: all cores but one).	`0`

Returns:

Type	Description
`InfluenceZonesResult`

Notes

Experimental: the method is GTB's, but its distance transform (IDL's MORPH_DISTANCE) is not yet reproduced exactly, so zone boundaries can differ in detail.

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))   # 0 missing, 1 background, 2 forest
>>> zones = ls.influence_zones(forest, min_area=5000)
>>> zones.n_objects
10
>>> zones.object_area[:3]
array([214811,   7927,   7398])

InfluenceZonesResult `dataclass` ¶

InfluenceZonesResult(
    classes: DataArray | ndarray,
    ids: DataArray | ndarray,
    zones: DataArray | ndarray,
    object_area: ndarray,
    min_area: int,
    core_buffer: int,
    buffer: int,
)

Bases: Saveable

Result of influence_zones().

Attributes:

Name	Type	Description
`classes`	`DataArray \| ndarray`	`uint8` map in GTB's viewport codes (see `INFLUENCE_CLASSES`).
`ids`	`DataArray \| ndarray`	`int32` object id on the object pixels (holes included); 0 zones, -1 holes, -2 foreground that is not part of an object, -3 watershed lines, -4 missing, -5 buffer boundary, -6 buffer.
`zones`	`DataArray \| ndarray`	`int32` zone (watershed basin) of every pixel, 0 on the lines.
`object_area`	`ndarray`	Area of object `k` (without its holes) at index `k - 1`.

n_objects `property` ¶

n_objects: int

Number of objects (of at least min_area pixels).

proximity ¶

proximity(
    data,
    min_area: int = 5000,
    *,
    core_buffer: int = 0,
    prox_max: int = 0,
    n_threads: int = 0,
) -> ProximityResult

Proximity between neighbouring objects (GTB Distance -> Proximity).

Parameters are those of influence_zones(), with prox_max: the watershed lines with a proximity up to prox_max (or the minimum proximity, if larger) are marked 165 (GTB default 0).

Parameters:

Name	Type	Description	Default
`data`		Map coded 0 missing, 1 background, 2 foreground (array, DataArray or path).	required
`min_area`	`int`	Minimum object area in pixels (GTB default 5000).	`5000`
`core_buffer`	`int`	Reduce the objects to their core first (see `influence_zones()`).	`0`
`prox_max`	`int`	Watershed lines with a proximity up to this value are marked 165.	`0`
`n_threads`	`int`	Threads (0 or `None`: all cores but one).	`0`

Returns:

Type	Description
`ProximityResult`

Notes

Experimental, as influence_zones().

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))   # 0 missing, 1 background, 2 forest
>>> near = ls.proximity(forest, min_area=5000)
>>> near.n_objects, near.prox_min
(10, 2)
>>> near.lines.columns.tolist()
['x', 'y', 'nbr', 'clength', 'cag', 'ids', 'areas']

ProximityResult `dataclass` ¶

ProximityResult(
    classes: DataArray | ndarray,
    ids: DataArray | ndarray,
    zones: DataArray | ndarray,
    distance: DataArray | ndarray,
    cag: DataArray | ndarray,
    object_area: ndarray,
    lines: DataFrame,
    prox_min: int | None,
    prox_min_count: int,
    min_area: int,
    core_buffer: int,
    prox_max: int,
    source: str = "",
)

Bases: Saveable

Result of proximity().

Attributes:

Name	Type	Description
`classes, ids, zones, object_area`		As in `InfluenceZonesResult` (objects are all coded 12; lines at or below `prox_max` 165, at the minimum proximity 193).
`distance`	`DataArray \| ndarray`	`int32` proximity map: `round(2 * distance to the nearest object)`, on lines between more than two zones scaled by their number.
`cag`	`DataArray \| ndarray`	`int64` connected area gain on the watershed lines: `round(sqrt(2 * sum_{i<j} a_i * a_j))` over the adjacent objects.
`lines`	`DataFrame`	One row per watershed pixel between zones (GTB's proximity CSV): 1-based column/row, number of adjacent zones, connector length, CAG and the adjacent zone ids with their object areas.
`prox_min, prox_min_count`		Minimum proximity on the watershed lines and the number of pixels with it.

n_objects `property` ¶

n_objects: int

Number of objects (of at least min_area pixels).

csv_text ¶

csv_text(source: str | None = None) -> str

GTB's proximity CSV (<name>_proximity_<min_area>_<core_buffer>_<prox_max>.csv) as text.

Returns:

Type	Description
`str`

to_csv ¶

to_csv(path: str | Path, source: str | None = None) -> Path

Write GTB's proximity CSV.

Parameters:

Name	Type	Description	Default
`path`	`str \| Path`	Output file.	required

Original vegetation¶

original_vegetation ¶

original_vegetation(
    land_cover,
    natural: Iterable[str],
    anthropic: Iterable[str],
    *,
    past=None,
    ecoregions=None,
    failsafe: Iterable[str] | None = None,
    threshold: float = 0.95,
    restore_to: str = "original",
    n_threads: int | None = None,
) -> OriginalVegetation

The original vegetation of converted land, and what restoration brings back.

Every pixel holds shares of land-cover classes: natural ones (forest, grassland, wetland...), anthropic ones (crops, pasture...) and others that are neither (water, rock). Restoring a pixel's anthropic share brings back natural vegetation; this function says which, as plangea prepares it for its restoration model (Strassburg et al. 2020):

the anthropic share is split among the natural classes in the proportions they had in the pixel in the past;
where the past has no natural cover, or more than threshold of the pixel is anthropic today, in the proportions of the pixel's ecoregions (equal proportions without them);
where those are all zero too, in the mean proportions of the failsafe classes over the pixels with anthropic land (over all the land with restore_to="current").

The original cover of a class is then its current share plus its part of the anthropic share.

Parameters:

Name	Type	Description	Default
`land_cover`		The current share (0 to 1) of every class in each pixel: a `(band, y, x)` raster whose band names are the classes, a mapping class -> raster, a list of named rasters or files. NaN counts as 0.	required
`natural`	`Iterable[str]`	The natural and the anthropic classes. The other layers are ignored.	required
`anthropic`	`Iterable[str]`	The natural and the anthropic classes. The other layers are ignored.	required
`past`		The past share of every natural class (same forms; NaN counts as 0). By default the current shares.	`None`
`ecoregions`		For every natural class, its proportion in each pixel's ecoregion (same forms). By default equal proportions.	`None`
`failsafe`	`Iterable[str] \| None`	The classes the failsafe fills (default: every natural class).	`None`
`threshold`	`float`	Above this anthropic share, the ecoregion's proportions are used even where there is a past cover.	`0.95`
`restore_to`	`str`	`"original"`: restored land becomes the original cover's natural classes in their proportions (restoration alone); `"current"`: the current natural classes of the pixel where it has some, else the original ones (plangea's restoration with conversion).	`'original'`
`n_threads`	`int \| None`	Threads (default: all cores but one).	`None`

Returns:

Type	Description
`OriginalVegetation`	The original cover, the restored proportions, the source of each pixel's proportions and the anthropic share.

Notes

The arithmetic is plangea's, bit for bit (harmonize_oa), with two corrections: without ecoregion maps, plangea gives every natural class the pixel's total natural share divided by the number of classes instead of equal proportions (its own comment says "flat among natural types"), which leaves most of a converted pixel to no class; and pixels the failsafe fills get 0 for the classes it does not fill, where plangea can stop with NA.

Examples:

>>> import numpy as np
>>> import xarray as xr
>>> import landschaft as ls
>>> cover = xr.DataArray(np.array([[[0.6, 0.0]], [[0.2, 0.1]], [[0.2, 0.9]]]),
...                      dims=("band", "y", "x"), coords={"band": ["forest", "savanna", "crops"]})
>>> past = xr.DataArray(np.array([[[0.9, 0.0]], [[0.1, 0.0]]]),
...                     dims=("band", "y", "x"), coords={"band": ["forest", "savanna"]})
>>> veg = ls.original_vegetation(cover, ["forest", "savanna"], ["crops"], past=past)
>>> veg.original.sel(band="forest").values.round(3)
array([[0.78, 0.45]])
>>> veg.source.values
array([[0., 1.]])

OriginalVegetation `dataclass` ¶

OriginalVegetation(
    original: DataArray,
    restore: DataArray,
    source: DataArray,
    anthropic: DataArray,
    failsafe: dict = dict(),
)

Bases: Saveable

Result of original_vegetation().

Attributes:

Name	Type	Description
`original`	`DataArray`	The share of every natural class in each pixel's original cover (`(band, y, x)`, one band per natural class): its current share plus the part of the anthropic share that was that class. NaN outside the land (no natural or anthropic cover).
`restore`	`DataArray`	The shares of the natural classes that restored land in each pixel becomes (they sum to 1).
`source`	`DataArray`	Where each pixel's proportions come from: 0 its past cover, 1 its ecoregion, 2 the failsafe (see `VEGETATION_SOURCES`).
`anthropic`	`DataArray`	The anthropic share of each pixel: the land restoration can work on.
`failsafe`	`dict`	The failsafe proportions of the natural classes.

classes `property` ¶

classes: list[str]

The natural classes, in the order of the bands.

Ecological intactness¶

intactness ¶

Ecological intactness from the human footprint (Beyer et al. 2020).

A cell's habitat quality is q = exp(-gamma HFI) for its human footprint index HFI, and its intactness Q_i = sum_j (q_i q_j)^z w_ij / sum_j w_ij: the quality of the cell and of its surroundings within a radius, weighted by w = exp(beta d) at a distance d (km). The computation and its parameters follow plangea's multi_calc_it (International Institute for Sustainability; GPL >= 3), which uses it as a benefit of restoration: restoring a share x of a cell lowers its footprint by a reduction per unit, and intactness_gradient() gives what that adds to the total intactness.

References

Beyer, H. L., Venter, O., Grantham, H. S. & Watson, J. E. M. (2020). Substantial losses in ecoregion intactness highlight urgency of globally coordinated action. Conservation Letters 13, e12692. doi:10.1111/conl.12692

intactness ¶

intactness(
    hfi,
    radius: int = 3,
    *,
    gamma: float = INTACTNESS_GAMMA,
    beta: float = -0.2,
    z: float = 0.5,
    cell_size=None,
    reduction=None,
    restored=None,
    n_threads: int | None = None,
) -> DataArray

Ecological intactness of every cell from the human footprint (Beyer et al. 2020).

The habitat quality of a cell is q = exp(-gamma HFI), from 1 where people leave no footprint down towards 0; its intactness

.. math:: Q_i = \frac{\sum_j (q_i q_j)^z w_{ij}}{\sum_j w_{ij}}

over the cells j within radius cells of it (itself included), with weights w = exp(beta d) that fall with the distance d in km: a cell is intact when it and its surroundings are of high quality.

Parameters:

Name	Type	Description	Default
`hfi`		The human footprint index (a raster, NaN outside the area).	required
`radius`	`int`	The radius of the neighbourhood, in cells (plangea's `n_pixels`): the cells whose centres are within `radius` cell widths.	`3`
`gamma`	`float`	Beyer et al.'s parameters: gamma = -ln(0.2) / 4 (a footprint of 4 gives a quality of 0.2), beta = -0.2 per km, z = 0.5.	`INTACTNESS_GAMMA`
`beta`	`float`	Beyer et al.'s parameters: gamma = -ln(0.2) / 4 (a footprint of 4 gives a quality of 0.2), beta = -0.2 per km, z = 0.5.	`INTACTNESS_GAMMA`
`z`	`float`	Beyer et al.'s parameters: gamma = -ln(0.2) / 4 (a footprint of 4 gives a quality of 0.2), beta = -0.2 per km, z = 0.5.	`INTACTNESS_GAMMA`
`cell_size`		The width (and height) of a cell in km, or `(width, height)`; by default from a projected grid (its unit converted to km).	`None`
`reduction`		Restoration: the footprint a cell loses per unit of its area restored (a raster or a number) and the share of each cell restored. The footprint becomes HFI - reduction x restored, at least 0.	`None`
`restored`		Restoration: the footprint a cell loses per unit of its area restored (a raster or a number) and the share of each cell restored. The footprint becomes HFI - reduction x restored, at least 0.	`None`
`n_threads`	`int \| None`	Threads (default: all cores but one).	`None`

Returns:

Type	Description
`DataArray`	Q on the footprint's grid, 0 to 1 (NaN where the footprint is).

Notes

The values are plangea's multi_calc_it (with aggregate = TRUE), bit for bit, except at the edges: plangea leaves out the cells closer than the radius to an edge (and one more row and column at the end) and gives them the mean of the others; here their window is the part inside the raster, as it is around cells without a footprint.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> hfi = np.zeros((5, 5))
>>> hfi[:, 3:] = 20.0                      # a city on the right
>>> q = ls.intactness(hfi, radius=1, cell_size=1.0)
>>> q.values[2].round(3)
array([1.   , 1.   , 0.812, 0.004, 0.   ])

intactness_gradient ¶

intactness_gradient(
    hfi,
    reduction,
    radius: int = 3,
    *,
    gamma: float = INTACTNESS_GAMMA,
    beta: float = -0.2,
    z: float = 0.5,
    cell_size=None,
    restored=None,
    n_threads: int | None = None,
) -> DataArray

What restoring each cell adds to the total ecological intactness.

Restoring a share x of a cell lowers its human footprint by reduction x; this is the derivative of the sum of intactness() over all cells with respect to x: the cell's own intactness rises, and so does that of every cell whose neighbourhood holds it,

.. math:: \frac{\partial \sum_i Q_i}{\partial x_k} = z\,\gamma\,a_k \Big(Q_k + q_k^z \sum_{i} \frac{q_i^z w_{ik}}{W_i}\Big)

(a the reduction, W_i the sum of the weights of i's window), 0 where the footprint is already 0. It is the linear benefit of restoration for intactness in a restoration plan.

Parameters:

Name	Description	Default
`hfi`	As `intactness()`.	required
`radius`	As `intactness()`.	required
`gamma`	As `intactness()`.	required
`beta`	As `intactness()`.	required
`z`	As `intactness()`.	required
`cell_size`	As `intactness()`.	required
`n_threads`	As `intactness()`.	required
`reduction`	The footprint a cell loses per unit of its area restored (a raster or a number): in plangea, the anthropic classes' shares weighted by their footprint (`hfi_weights`).	required
`restored`	Shares already restored (the derivative at that point).	`None`

Returns:

Type	Description
`DataArray`	The derivative on the footprint's grid.

Notes

plangea's derivative changes only the cell's own factor q_i (by a finite difference), leaving out the cell inside its own window and its effect on its neighbours' intactness; this is the complete one (the derivative of the total intactness).

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> hfi = np.full((5, 5), 4.0)
>>> g = ls.intactness_gradient(hfi, reduction=2.0, radius=1, cell_size=1.0)
>>> round(float(g.values[2, 2]), 4)
0.1609

intactness_gradient ¶

intactness_gradient(
    hfi,
    reduction,
    radius: int = 3,
    *,
    gamma: float = INTACTNESS_GAMMA,
    beta: float = -0.2,
    z: float = 0.5,
    cell_size=None,
    restored=None,
    n_threads: int | None = None,
) -> DataArray

What restoring each cell adds to the total ecological intactness.

Restoring a share x of a cell lowers its human footprint by reduction x; this is the derivative of the sum of intactness() over all cells with respect to x: the cell's own intactness rises, and so does that of every cell whose neighbourhood holds it,

.. math:: \frac{\partial \sum_i Q_i}{\partial x_k} = z\,\gamma\,a_k \Big(Q_k + q_k^z \sum_{i} \frac{q_i^z w_{ik}}{W_i}\Big)

(a the reduction, W_i the sum of the weights of i's window), 0 where the footprint is already 0. It is the linear benefit of restoration for intactness in a restoration plan.

Parameters:

Name	Description	Default
`hfi`	As `intactness()`.	required
`radius`	As `intactness()`.	required
`gamma`	As `intactness()`.	required
`beta`	As `intactness()`.	required
`z`	As `intactness()`.	required
`cell_size`	As `intactness()`.	required
`n_threads`	As `intactness()`.	required
`reduction`	The footprint a cell loses per unit of its area restored (a raster or a number): in plangea, the anthropic classes' shares weighted by their footprint (`hfi_weights`).	required
`restored`	Shares already restored (the derivative at that point).	`None`

Returns:

Type	Description
`DataArray`	The derivative on the footprint's grid.

Notes

plangea's derivative changes only the cell's own factor q_i (by a finite difference), leaving out the cell inside its own window and its effect on its neighbours' intactness; this is the complete one (the derivative of the total intactness).

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> hfi = np.full((5, 5), 4.0)
>>> g = ls.intactness_gradient(hfi, reduction=2.0, radius=1, cell_size=1.0)
>>> round(float(g.values[2, 2]), 4)
0.1609

Distance & restoration¶

Distance to the habitat¶

euclidean_distance ¶

EuclideanDistanceResult dataclass ¶

summary ¶

distance_to_resistance ¶

Restoration status¶

restoration_status ¶

RestorationStatusResult dataclass ¶

to_dict ¶

to_dataframe ¶

summary ¶

Costs, paths and corridors¶

cost_distance ¶

least_cost_path ¶

LeastCostPathResult dataclass ¶

restoration_corridor ¶

RestorationCorridorResult dataclass ¶

to_dataframe ¶

isochrone_map ¶

isochrone_map_ab ¶

IsochroneResult dataclass ¶

Zones of influence¶

influence_zones ¶

InfluenceZonesResult dataclass ¶

n_objects property ¶

proximity ¶

ProximityResult dataclass ¶

n_objects property ¶

csv_text ¶

to_csv ¶

Original vegetation¶

original_vegetation ¶

OriginalVegetation dataclass ¶

classes property ¶

Ecological intactness¶

intactness ¶

intactness ¶

intactness_gradient ¶

intactness_gradient ¶

EuclideanDistanceResult `dataclass` ¶

RestorationStatusResult `dataclass` ¶

LeastCostPathResult `dataclass` ¶

RestorationCorridorResult `dataclass` ¶

IsochroneResult `dataclass` ¶

InfluenceZonesResult `dataclass` ¶

n_objects `property` ¶

ProximityResult `dataclass` ¶

n_objects `property` ¶

OriginalVegetation `dataclass` ¶

classes `property` ¶