Habitat patterns¶

Structural classes of habitat pixels (MSPA, SPA), objects by size, parcellation and boundary contortion. Tutorial: habitat patterns.

MSPA¶

mspa ¶

Morphological Spatial Pattern Analysis (MSPA).

Soille, P. & Vogt, P. (2009). Morphological segmentation of binary patterns. Pattern Recognition Letters 30(4), 456-459. doi:10.1016/j.patrec.2008.10.015

mspa ¶

mspa(
    data,
    *,
    connectivity: int = 8,
    edge_width: int = 1,
    transition: bool = True,
    internal: bool = True,
    n_threads: int | None = None,
    validate: bool = True,
    output_dir=None,
    n_workers: int | None = None,
    max_memory: int | str | None = None,
    overwrite: bool = False,
)

Morphological Spatial Pattern Analysis, identical to GuidosToolbox's MSPA.

Parameters:

Name	Type	Description	Default
`data`		Map coded `0` = missing, `1` = background, `2` = foreground (`uint8`). Accepts a 2-D numpy array, a dask array or an `xarray.DataArray`. Arrays with more dimensions (e.g. `(time, y, x)`) are processed slice by slice over the last two (or `y`/`x`) dimensions; with dask the slices run in parallel. A `RasterCollection` (`load_raster` of a list, directory or glob) is processed one raster per thread, see `output_dir`.	required
`connectivity`	`int`	Foreground connectivity, 8 (default) or 4 (GTB "FGconn").	`8`
`edge_width`	`int`	Effective edge width in pixels, 1-255 (GTB "EdgeWidth").	`1`
`transition`	`bool`	Whether loops/bridges crossing edges or perforations are shown as transitions (GTB "Transition"). As in GuidosToolbox this only changes the colour table written by `write_raster()`, not the values.	`True`
`internal`	`bool`	Distinguish internal features with a +100 offset and label core-openings (100) and border-openings (220) (GTB "Intext").	`True`
`n_threads`	`int \| None`	Threads used inside one map; `1` runs sequentially. `None` uses all cores but one, except for dask-backed input where it defaults to 1 because dask already runs the slices in parallel.	`None`
`validate`	`bool`	Check the input like GuidosToolbox before processing.	`True`
`output_dir`		Collections only. With `output_dir` every result is written there as soon as it is ready (GTB naming) and a list of `BatchResult` is returned, so memory stays bounded; see `mspa_batch()`. Without it, the results are returned as a list of DataArrays. `n_workers` rasters run concurrently (default: all cores but one), within the `max_memory` budget (e.g. `"16GB"`).	`None`
`n_workers`		Collections only. With `output_dir` every result is written there as soon as it is ready (GTB naming) and a list of `BatchResult` is returned, so memory stays bounded; see `mspa_batch()`. Without it, the results are returned as a list of DataArrays. `n_workers` rasters run concurrently (default: all cores but one), within the `max_memory` budget (e.g. `"16GB"`).	`None`
`max_memory`		Collections only. With `output_dir` every result is written there as soon as it is ready (GTB naming) and a list of `BatchResult` is returned, so memory stays bounded; see `mspa_batch()`. Without it, the results are returned as a list of DataArrays. `n_workers` rasters run concurrently (default: all cores but one), within the `max_memory` budget (e.g. `"16GB"`).	`None`
`overwrite`		Collections only. With `output_dir` every result is written there as soon as it is ready (GTB naming) and a list of `BatchResult` is returned, so memory stays bounded; see `mspa_batch()`. Without it, the results are returned as a list of DataArrays. `n_workers` rasters run concurrently (default: all cores but one), within the `max_memory` budget (e.g. `"16GB"`).	`None`

Returns:

Type	Description
Same container type as ``data`` with the ``uint8`` MSPA classes; see
`MSPA_CLASSES`. DataArrays keep coordinates and georeferencing and
record the parameters in ``attrs``. For a collection, a list (see
``output_dir``).

mspa_stats ¶

mspa_stats(result) -> DataFrame

Area of every MSPA class.

Parameters:

Name	Type	Description	Default
`result`		An MSPA map returned by `mspa()` (DataArray or numpy array).	required

Returns:

Type	Description
`DataFrame`	Indexed by class (Core, Islet, Perforation, Edge, Loop, Bridge, Branch, Background, Core-Opening, Border-Opening, Missing), with `pixels`, the percentage of the mapped area `pct_data` (all non-missing pixels) and, for habitat classes, the percentage of the habitat `pct_foreground`. Classes inside core-openings count with their class.

Examples:

>>> import landschaft as ls
>>> morph = ls.mspa(ls.load_raster(ls.example_data("input")), edge_width=2)
>>> stats = ls.mspa_stats(morph)
>>> list(stats.columns)
['pixels', 'pct_data', 'pct_foreground']
>>> round(float(stats.loc["Core", "pct_foreground"]), 1)
61.8

mspa_stats ¶

mspa_stats(result) -> DataFrame

Area of every MSPA class.

Parameters:

Name	Type	Description	Default
`result`		An MSPA map returned by `mspa()` (DataArray or numpy array).	required

Returns:

Type	Description
`DataFrame`	Indexed by class (Core, Islet, Perforation, Edge, Loop, Bridge, Branch, Background, Core-Opening, Border-Opening, Missing), with `pixels`, the percentage of the mapped area `pct_data` (all non-missing pixels) and, for habitat classes, the percentage of the habitat `pct_foreground`. Classes inside core-openings count with their class.

Examples:

>>> import landschaft as ls
>>> morph = ls.mspa(ls.load_raster(ls.example_data("input")), edge_width=2)
>>> stats = ls.mspa_stats(morph)
>>> list(stats.columns)
['pixels', 'pct_data', 'pct_foreground']
>>> round(float(stats.loc["Core", "pct_foreground"]), 1)
61.8

MSPA class codes¶

ls.MSPA_CLASSES maps each code of an MSPA map to its class:

Code	Class	Code	Class
17	Core	65	Loop
9	Islet	67	Loop in edge
5	Perforation	69	Loop in perforation
3	Edge	33	Bridge
1	Branch	35	Bridge in edge
0	Background	37	Bridge in perforation
100	Core-opening	220	Border-opening
129	Missing	+100	the class inside a core-opening (`internal=True`)

SPA¶

spa ¶

Simplified Pattern Analysis (GuidosToolbox Pattern -> Morphological -> SPA2/3/5/6).

Simplified alternatives to MSPA for the 8-connected foreground, with a fixed edge width of one pixel:

SPA2: contiguous foreground and small linear features (SLF: foreground removed by a 3x3 opening);
SPA3: core, core-opening and margin;
SPA5: core, edge, perforation, margin (and core-opening);
SPA6: SPA5 plus islets (foreground objects without core).

The maps and statistics are identical to GTB's.

SPAResult `dataclass` ¶

SPAResult(
    map: DataArray | ndarray,
    variant: int,
    stats: DataFrame,
    counts: dict[str, int],
    areas: dict[str, int],
    porosity: float,
)

Bases: Saveable

Result of spa().

Attributes:

Name	Type	Description
`map`	`DataArray \| ndarray`	`uint8` map in GTB's codes (see `SPA_CLASSES`).
`stats`	`DataFrame`	Foreground classes with area [pixels] and their share of the data area and of the foreground (GTB's SPA report).
`counts`	`dict[str, int]`	Number of 8-connected objects per class (core-openings: 4-connected) and of foreground objects, as GTB reports them.
`porosity`	`float`	`100 - contiguous / (contiguous + core-openings) * 100`.

Examples:

>>> import landschaft as ls
>>> result = ls.spa(ls.load_raster(ls.example_data("input")), variant=6)
>>> result.map.shape
(1000, 1000)

summary ¶

summary() -> str

Report in the layout of GTB's spaN.txt (areas and proportions).

spa ¶

spa(data, variant: int = 5, *, n_threads: int = 0) -> SPAResult

Simplified pattern analysis: the habitat in 2, 3, 5 or 6 classes.

A coarser, faster summary than mspa():

2: contiguous habitat (core, edge, perforation) and small linear features (margin);
3: core, edge (with perforation) and margin;
5: core, edge, perforation and margin;
6: core, edge, perforation, margin and islet.

Parameters:

Name	Type	Description	Default
`data`		2-D `uint8` habitat map (array, DataArray or raster path) coded 0 missing, 1 background, 2 habitat.	required
`variant`	`int`	2, 3, 5 (default) or 6.	`5`
`n_threads`	`int`	Threads (`None` / 0: all cores but one).	`0`

Returns:

Type	Description
`SPAResult`	The class map, class areas, object counts and porosity.

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))
>>> result = ls.spa(forest, variant=5)
>>> list(result.stats["class"])
['Core', 'Edge', 'Perforation', 'Margin']
>>> round(result.porosity, 2)
1.42
>>> print(result.summary())

SPAResult `dataclass` ¶

SPAResult(
    map: DataArray | ndarray,
    variant: int,
    stats: DataFrame,
    counts: dict[str, int],
    areas: dict[str, int],
    porosity: float,
)

Bases: Saveable

Result of spa().

Attributes:

Name	Type	Description
`map`	`DataArray \| ndarray`	`uint8` map in GTB's codes (see `SPA_CLASSES`).
`stats`	`DataFrame`	Foreground classes with area [pixels] and their share of the data area and of the foreground (GTB's SPA report).
`counts`	`dict[str, int]`	Number of 8-connected objects per class (core-openings: 4-connected) and of foreground objects, as GTB reports them.
`porosity`	`float`	`100 - contiguous / (contiguous + core-openings) * 100`.

Examples:

>>> import landschaft as ls
>>> result = ls.spa(ls.load_raster(ls.example_data("input")), variant=6)
>>> result.map.shape
(1000, 1000)

summary ¶

summary() -> str

Report in the layout of GTB's spaN.txt (areas and proportions).

Objects by size¶

accounting ¶

Accounting: foreground objects grouped into size classes.

Reimplements GuidosToolbox's Image Analysis → Objects → Accounting (and its batch version): every foreground object is labelled, measured and assigned to one of up to six size classes. The maps and numbers are identical to the GTB outputs <name>_acc.tif, <name>_acc_ids.tif, <name>_acc_pixels.tif and <name>_acc.txt/.csv.

Reference

Vogt, P. & Riitters, K. (2017). GuidosToolbox: universal digital image object analysis. European Journal of Remote Sensing, 50(1), 352-361. doi:10.1080/22797254.2017.1330650

AccountingResult `dataclass` ¶

AccountingResult(
    classes: DataArray | ndarray,
    ids: DataArray | ndarray,
    area: DataArray | ndarray,
    stats: DataFrame,
    object_area: ndarray,
    largest: tuple[tuple[int, ...], ...],
    thresholds: tuple[int, ...],
    connectivity: int,
    pixel_resolution: float,
    highlighted: bool,
)

Bases: Saveable

Result of accounting().

Attributes:

Name	Type	Description
`classes`	`DataArray \| ndarray`	`uint8` map in GTB's colour codes (see `ACCOUNTING_CLASSES`).
`ids`	`DataArray \| ndarray`	`int32` object ids `1..n_objects` (numbered like GTB: raster order starting from the bottom row); 0 background, -1 missing, -3/-4 special background.
`area`	`DataArray \| ndarray`	`int32` map of each object's area in pixels, same codes as `ids`.
`stats`	`DataFrame`	One row per size class: range, colour, number of objects, area and their percentages of all objects / of the total foreground area.
`object_area`	`ndarray`	Area of object `k` at index `k - 1`.
`largest`	`tuple[tuple[int, ...], ...]`	Ids of the largest, second and third largest objects (objects of equal area share a rank).
`highlighted`	`bool`	Whether the three largest objects are painted pink in `classes`.

Examples:

>>> import landschaft as ls
>>> result = ls.accounting(ls.load_raster(ls.example_data("input")))
>>> result.largest[0]
(1,)
>>> int(result.object_area[0])
214811

median_area `property` ¶

median_area: int

Median object area as GTB reports it (the upper middle value for an even count).

std_area `property` ¶

std_area: float

Sample standard deviation of the object areas (NaN for a single object).

hectare_factor `property` ¶

hectare_factor: float

Hectares per pixel for pixel_resolution metres.

summary ¶

summary() -> str

Report in the layout of GTB's acc.txt.

accounting ¶

accounting(
    data: DataArray | ndarray | str | Path,
    thresholds: Sequence[int] | str = DEFAULT_THRESHOLDS,
    *,
    connectivity: int = 8,
    highlight_largest: bool = False,
    pixel_resolution: float = 25.0,
    n_threads: int = 0,
) -> AccountingResult

Group the foreground objects into size classes (GTB Accounting).

Parameters:

Name	Type	Description	Default
`data`	`DataArray \| ndarray \| str \| Path`	2-D `uint8` map (array, DataArray or raster path) coded 0 = missing, 1 = background, 2 = foreground, 3/4 = special background.	required
`thresholds`	`Sequence[int] \| str`	Upper bounds [pixels] of the size classes, 1 to 5 non-decreasing values >= 1 (repeats are dropped, as in GTB). Class 1 holds objects of up to `thresholds[0]` pixels, class k those in `(thresholds[k-2], thresholds[k-1]]` and the last class the larger ones. The default is GTB's batch setting; `"auto"` uses the classes GTB's interactive Accounting proposes (0.1, 1, 10, 50 and 90% of the range of the 8-connected object areas).	`DEFAULT_THRESHOLDS`
`connectivity`	`int`	Foreground connectivity, 8 (default) or 4.	`8`
`highlight_largest`	`bool`	Paint the three largest objects pink (code 80), GTB's "big3pink". Only applied when there are at least three objects, as in GTB.	`False`
`pixel_resolution`	`float`	Pixel size in metres, used for the hectare/acre factors of the report.	`25.0`
`n_threads`	`int`	OpenMP threads (0: all cores but one).	`0`

Returns:

Type	Description
`AccountingResult`

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))
>>> result = ls.accounting(forest, thresholds=(50, 250, 1000, 5000, 20000))
>>> result.n_objects
2850
>>> result.stats[["size_class", "range", "objects"]].head(2)
   size_class     range  objects
0           1    1 - 50     2668
1           2  51 - 250      129
>>> _ = ls.write_raster(result.classes, "size_classes.tif")

Classes proposed from the range of object sizes:

>>> auto = ls.accounting(forest, thresholds="auto")

AccountingResult `dataclass` ¶

AccountingResult(
    classes: DataArray | ndarray,
    ids: DataArray | ndarray,
    area: DataArray | ndarray,
    stats: DataFrame,
    object_area: ndarray,
    largest: tuple[tuple[int, ...], ...],
    thresholds: tuple[int, ...],
    connectivity: int,
    pixel_resolution: float,
    highlighted: bool,
)

Bases: Saveable

Result of accounting().

Attributes:

Name	Type	Description
`classes`	`DataArray \| ndarray`	`uint8` map in GTB's colour codes (see `ACCOUNTING_CLASSES`).
`ids`	`DataArray \| ndarray`	`int32` object ids `1..n_objects` (numbered like GTB: raster order starting from the bottom row); 0 background, -1 missing, -3/-4 special background.
`area`	`DataArray \| ndarray`	`int32` map of each object's area in pixels, same codes as `ids`.
`stats`	`DataFrame`	One row per size class: range, colour, number of objects, area and their percentages of all objects / of the total foreground area.
`object_area`	`ndarray`	Area of object `k` at index `k - 1`.
`largest`	`tuple[tuple[int, ...], ...]`	Ids of the largest, second and third largest objects (objects of equal area share a rank).
`highlighted`	`bool`	Whether the three largest objects are painted pink in `classes`.

Examples:

>>> import landschaft as ls
>>> result = ls.accounting(ls.load_raster(ls.example_data("input")))
>>> result.largest[0]
(1,)
>>> int(result.object_area[0])
214811

median_area `property` ¶

median_area: int

Median object area as GTB reports it (the upper middle value for an even count).

std_area `property` ¶

std_area: float

Sample standard deviation of the object areas (NaN for a single object).

hectare_factor `property` ¶

hectare_factor: float

Hectares per pixel for pixel_resolution metres.

summary ¶

summary() -> str

Report in the layout of GTB's acc.txt.

Parcels and shape¶

parcellation ¶

parcellation(
    data: DataArray | ndarray | str | Path,
    connectivity: int = 8,
    return_labels: bool = False,
    *,
    n_threads: int | None = None,
) -> ParcellationResult

How finely the landscape is divided into parcels.

Every connected region of one value (e.g. habitat and background, or the classes of a land-cover map) is a parcel. Per class and for the whole map: the number of parcels, the average parcel size (APS), the area-weighted average parcel size (AWAPS, the size of the parcel a random pixel falls in), the division index (probability that two random pixels are in different parcels) and the degree of parcellation (PARC).

Parameters:

Name	Type	Description	Default
`data`	`DataArray \| ndarray \| str \| Path`	Class raster (array, DataArray or path; integer classes, 0 = NoData).	required
`connectivity`	`int`	Pixels of one parcel touch by a side (4) or also by a corner (8).	`8`
`return_labels`	`bool`	Also return a map of the parcels.	`False`
`n_threads`	`int \| None`	Threads (`None` / 0: all cores but one).	`None`

Returns:

Type	Description
`ParcellationResult`	`stats` per class, `global_stats` for the map and, with `return_labels`, `parcel_labels`.

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))
>>> result = ls.parcellation(forest)
>>> list(result.stats.columns)
['Class', 'Count', 'Area[pixels]', 'APS', 'AWAPS', 'AWAPS/data', 'DIVISION', 'PARC[%]']
>>> result.global_stats["total_parcels"]
4020
>>> labels = ls.parcellation(forest, return_labels=True).parcel_labels

ParcellationResult `dataclass` ¶

ParcellationResult(
    stats: DataFrame,
    global_stats: dict[str, Any],
    parcel_labels: DataArray | ndarray | None = None,
)

Bases: Saveable

Results from landscape parcellation analysis.

Attributes:

Name	Type	Description
`stats`	`DataFrame`	Class-level parcellation metrics (Count, Area, APS, AWAPS, DIVISION, PARC).
`global_stats`	`dict[str, Any]`	Landscape-level metrics across all valid classes.
`parcel_labels`	`(DataArray or ndarray, optional)`	Labeled parcels raster (if return_labels=True was requested).

Examples:

>>> import landschaft as ls
>>> result = ls.parcellation(ls.load_raster(ls.example_data("input")))
>>> round(result.global_stats["division"], 3)
0.689

summary ¶

summary() -> str

Return a formatted string summary of parcellation metrics.

contortion ¶

contortion(
    data: DataArray | ndarray | str | Path,
    min_area: int = 1,
    return_display: bool = False,
) -> DataArray | ndarray

How convoluted the boundary of every habitat object is.

For each habitat object, a measure of how much its boundary bends relative to its size: compact objects score low, elongated or dendritic ones high.

Parameters:

Name	Type	Description	Default
`data`	`DataArray \| ndarray \| str \| Path`	Habitat map (array, DataArray or path) coded 0 missing, 1 background, 2 habitat.	required
`min_area`	`int`	Objects smaller than this (pixels), or without an interior (no 3 x 3 core), are left out and coded 2.	`1`
`return_display`	`bool`	Rescale the values to 3-255 for display (the byte scale of the colour table `contortion_colormap()`).	`False`

Returns:

Type	Description
`DataArray or ndarray`	0 missing, 1 background, 2 left-out objects, 3 and above the contortion of each object.

Examples:

>>> import landschaft as ls
>>> forest = ls.load_raster(ls.example_data("input"))
>>> values = ls.contortion(forest, min_area=100)
>>> int(values.max())
42

Habitat patterns¶

MSPA¶

mspa ¶

mspa ¶

mspa_stats ¶

mspa_stats ¶

MSPA class codes¶

SPA¶

spa ¶

SPAResult dataclass ¶

summary ¶

spa ¶

SPAResult dataclass ¶

summary ¶

Objects by size¶

accounting ¶

AccountingResult dataclass ¶

median_area property ¶

std_area property ¶

hectare_factor property ¶

summary ¶

accounting ¶

AccountingResult dataclass ¶

median_area property ¶

std_area property ¶

hectare_factor property ¶

summary ¶

Parcels and shape¶

parcellation ¶

ParcellationResult dataclass ¶

summary ¶

contortion ¶

SPAResult `dataclass` ¶

SPAResult `dataclass` ¶

AccountingResult `dataclass` ¶

median_area `property` ¶

std_area `property` ¶

hectare_factor `property` ¶

AccountingResult `dataclass` ¶

median_area `property` ¶

std_area `property` ¶

hectare_factor `property` ¶

ParcellationResult `dataclass` ¶