Surface metrics¶

Measures of the texture of a continuous raster (NDVI, elevation, canopy height, any gradient) that need no classification into patches: how much it varies, how its heights are distributed, how its trend can be removed. They are the areal surface texture parameters of ISO 25178, brought to landscape ecology as gradient surface metrics, and live in the ls.surface module.

The metrics follow the R package geodiv, bit for bit wherever geodiv is right; where it errs (Sq and the moments built on it, the seven-point slopes, the valleys, the summit curvature, the areas of the bearing area curve, the centre of its Fourier images, the fractal dimension's prisms, the choice of trend surface), Landschaft computes what the method intends, and each function's notes say what differs.

import landschaft as ls

ndvi = ls.surface.example_data("orforest")
z = ls.surface.remove_plane(ndvi)       # residuals from the best trend surface
ls.surface.metrics(z)                   # every metric in one table

All metrics at once¶

metrics ¶

metrics(
    x,
    metrics: str | Iterable[str] | None = None,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> Series

Surface metrics of a raster, as a Series indexed by metric.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (a single band; NaN or `nodata` = missing), usually with its trend removed (`remove_plane()`).	required
`metrics`	`str or list of str`	Names from `list_metrics()` (default: all). Metrics that share a computation are computed once.	`None`
`extent`	`tuple of float`	`(xmin, xmax, ymin, ymax)` of an array (default: unit cells).	`None`
`crs`	`optional`	CRS of an array.	`None`
`band`	`int`	Band of a multi-band file or DataArray (1-based).	`None`
`nodata`	`float`	Value that marks missing cells.	`None`

Returns:

Type	Description
`Series`

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> ls.surface.metrics(z, ["sa", "sq", "sku"]).round(4)
sa     0.0443
sq     0.0648
sku    6.4979
dtype: float64

list_metrics ¶

list_metrics() -> DataFrame

The surface metrics metrics() computes.

Returns:

Type	Description
`DataFrame`	Columns `metric`, `group` and `description`.

Examples:

>>> import landschaft as ls
>>> ls.surface.list_metrics().head(3)
  metric      group                 description
0     sa  amplitude           average roughness
1     sq  amplitude  root mean square roughness
2    ssk  amplitude         skewness (adjusted)

Amplitude¶

sa ¶

sa(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Average roughness (Sa): the mean absolute deviation from the mean height.

Sa = sum(|z - mean(z)|) / N over the N cells with a value.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (a single band; NaN or `nodata` = missing). Remove a trend first where there is one (`remove_plane()`).	required
`extent`	`tuple of float`	`(xmin, xmax, ymin, ymax)` of an array (default: unit cells).	`None`
`crs`	`optional`	CRS of an array.	`None`
`band`	`int`	Band of a multi-band file or DataArray (1-based).	`None`
`nodata`	`float`	Value that marks missing cells.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.sa(z), 6)
0.044294

sq ¶

sq(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Root mean square roughness (Sq): the root mean square of the heights about their mean, sqrt(sum((z - mean(z))^2) / N) (ISO 25178; geodiv uses the sample standard deviation, over N - 1).

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `sa()`).	required
`extent`		As in `sa()`.	`None`
`crs`		As in `sa()`.	`None`
`band`		As in `sa()`.	`None`
`nodata`		As in `sa()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.sq(z), 6)
0.064846

ssk ¶

ssk(
    x,
    adj: bool = True,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> float

Skewness of the height distribution (Ssk).

(sum((z - mean(z))^3) / N) / Sq^3 with Sq from sq(); with adj, multiplied by sqrt(N (N - 1)) / (N - 2) (the adjusted Fisher-Pearson coefficient). Negative: valleys dominate; positive: peaks dominate. (geodiv divides by the sample standard deviation cubed.)

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `sa()`).	required
`adj`	`bool`	Adjusted coefficient.	`True`
`extent`		As in `sa()`.	`None`
`crs`		As in `sa()`.	`None`
`band`		As in `sa()`.	`None`
`nodata`		As in `sa()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.ssk(z), 4), round(ls.surface.ssk(z, adj=False), 4)
(-1.9341, -1.9341)

sku ¶

sku(
    x,
    excess: bool = True,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> float

Kurtosis of the height distribution (Sku).

(sum((z - mean(z))^4) / N) / Sq^4 with Sq from sq(); with excess, minus 3 (0 for a normal distribution). (geodiv divides by the sample standard deviation to the fourth.)

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `sa()`).	required
`excess`	`bool`	The excess kurtosis.	`True`
`extent`		As in `sa()`.	`None`
`crs`		As in `sa()`.	`None`
`band`		As in `sa()`.	`None`
`nodata`		As in `sa()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.sku(z), 4), round(ls.surface.sku(z, excess=False), 4)
(6.4979, 9.4979)

sv ¶

sv(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Maximum valley depth (Sv): the absolute value of the lowest height.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `sa()`), usually with its mean or trend removed.	required
`extent`		As in `sa()`.	`None`
`crs`		As in `sa()`.	`None`
`band`		As in `sa()`.	`None`
`nodata`		As in `sa()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.sv(z), 6)
0.585464

sph ¶

sph(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Maximum peak height (Sp): the absolute value of the highest height.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `sa()`), usually with its mean or trend removed.	required
`extent`		As in `sa()`.	`None`
`crs`		As in `sa()`.	`None`
`band`		As in `sa()`.	`None`
`nodata`		As in `sa()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.sph(z), 6)
0.158592

smean ¶

smean(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Mean peak height: the mean of the positive heights.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `sa()`), with its mean or trend removed.	required
`extent`		As in `sa()`.	`None`
`crs`		As in `sa()`.	`None`
`band`		As in `sa()`.	`None`
`nodata`		As in `sa()`.	`None`

Returns:

Type	Description
`float`	NaN without positive heights.

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.smean(z), 6)
0.036899

Slopes and developed area¶

sdq ¶

sdq(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Root mean square slope (Sdq).

The root of the mean of (z[r, c] - z[r, c+1])^2 + (z[r, c] - z[r+1, c])^2, the differences with the right and lower neighbours in height units per cell, over the cells that have both.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (a single band; NaN or `nodata` = missing).	required
`extent`	`tuple of float`	`(xmin, xmax, ymin, ymax)` of an array (default: unit cells).	`None`
`crs`	`optional`	CRS of an array.	`None`
`band`	`int`	Band of a multi-band file or DataArray (1-based).	`None`
`nodata`	`float`	Value that marks missing cells.	`None`

Returns:

Type	Description
`float`

Notes

geodiv divides by the number of cells, those without both neighbours too.

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.sdq(z), 6)
0.04389

sdq6 ¶

sdq6(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Area root mean square slope (Sdq6), from seven-point derivatives.

The slopes in x and y from ISO 25178's sixth-order central difference, (z[+3] - 9 z[+2] + 45 z[+1] - 45 z[-1] + 9 z[-2] - z[-3]) / 60 per cell, at the cells three or more cells away from the edges; the root of the mean of their squared sums.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `sdq()`).	required
`extent`		As in `sdq()`.	`None`
`crs`		As in `sdq()`.	`None`
`band`		As in `sdq()`.	`None`
`nodata`		As in `sdq()`.	`None`

Returns:

Type	Description
`float`	NaN on rasters with fewer than 7 rows or columns.

Notes

geodiv's stencil has - 9 z[+2] twice in place of + 9 z[-2] (it does not cancel the height, so a high surface gives large values) and it takes the root of the sum, not of the mean.

Examples:

>>> import landschaft as ls
>>> dem = ls.surface.example_data("orelevation")
>>> round(ls.surface.sdq6(dem), 4)
10.6134

sdr ¶

sdr(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Developed interfacial area ratio (Sdr), in percent.

How much larger the surface is than its projection: (surface_area() - flat_area()) / flat_area() x 100, with the heights rescaled to 0-1 as geodiv.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `sdq()`).	required
`extent`		As in `sdq()`.	`None`
`crs`		As in `sdq()`.	`None`
`band`		As in `sdq()`.	`None`
`nodata`		As in `sdq()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.sdr(z), 4)
0.1733

surface_area ¶

surface_area(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Area of the surface, triangulated, in geodiv's units.

Every block of 2 x 2 cell centres is split into two triangles along each diagonal and the two areas are averaged; the sum over the blocks whose four values are all present. As geodiv, on the heights rescaled to 0-1 and cells 1 wide by their height-to-width ratio high, which makes the area independent of the heights' unit.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `sdq()`).	required
`extent`		As in `sdq()`.	`None`
`crs`		As in `sdq()`.	`None`
`band`		As in `sdq()`.	`None`
`nodata`		As in `sdq()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.surface_area(z), 2)
137137.26

flat_area ¶

flat_area(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Area of the surface's projection on a plane (geodiv's flatsa).

surface_area() of the same cells with every height 0.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `sdq()`).	required
`extent`		As in `sdq()`.	`None`
`crs`		As in `sdq()`.	`None`
`band`		As in `sdq()`.	`None`
`nodata`		As in `sdq()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.flat_area(z), 2)
136900.0

Peaks and valleys¶

find_peaks ¶

find_peaks(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> DataFrame

The peaks of a surface: cells higher than their 8 neighbours and than 0.

Only cells off the edge whose 8 neighbours all have a value can be peaks; heights are compared with 0, so remove the mean or trend first (remove_plane()).

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (a single band; NaN or `nodata` = missing).	required
`extent`	`tuple of float`	`(xmin, xmax, ymin, ymax)` of an array (default: unit cells).	`None`
`crs`	`optional`	CRS of an array.	`None`
`band`	`int`	Band of a multi-band file or DataArray (1-based).	`None`
`nodata`	`float`	Value that marks missing cells.	`None`

Returns:

Type	Description
`DataFrame`	One row per peak, in row order: the coordinates of the cell centre (`x`, `y`), its `value`, its `row` and `col` (from 0).

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> peaks = ls.surface.find_peaks(z)
>>> len(peaks), round(float(peaks.value.max()), 4)
(6488, 0.1586)

find_valleys ¶

find_valleys(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> DataFrame

The valleys of a surface: cells lower than their 8 neighbours and than 0.

As find_peaks(), downwards.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `find_peaks()`).	required
`extent`		As in `find_peaks()`.	`None`
`crs`		As in `find_peaks()`.	`None`
`band`		As in `find_peaks()`.	`None`
`nodata`		As in `find_peaks()`.	`None`

Returns:

Type	Description
`DataFrame`	As `find_peaks()`.

Notes

geodiv's findvalleys reads each cell's neighbours from the values arranged by columns, not by rows: on its normforest it finds 33041 "valleys" where there are 5167.

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> len(ls.surface.find_valleys(z))
5167

sds ¶

sds(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Density of summits (Sds): peaks per cell.

The number of find_peaks() over (rows - 1) (columns - 1).

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `find_peaks()`).	required
`extent`		As in `find_peaks()`.	`None`
`crs`		As in `find_peaks()`.	`None`
`band`		As in `find_peaks()`.	`None`
`nodata`		As in `find_peaks()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.sds(z), 6)
0.047392

ssc ¶

ssc(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Mean summit curvature (Ssc).

ISO 25178's -1 / (2n) times the sum of the second derivatives in x and y at the n find_peaks(), (z[x+1] - 2 z + z[x-1]) / dx^2 and (z[y+1] - 2 z + z[y-1]) / dy^2, with geodiv's units: heights rescaled to 0-1, dx = 1 and dy the ratio of the cell's height to its width. Positive: peaks are rounded upwards.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `find_peaks()`).	required
`extent`		As in `find_peaks()`.	`None`
`crs`		As in `find_peaks()`.	`None`
`band`		As in `find_peaks()`.	`None`
`nodata`		As in `find_peaks()`.	`None`

Returns:

Type	Description
`float`	NaN without peaks.

Notes

geodiv has the opposite sign and first differences divided by dx, not second differences divided by dx^2.

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.ssc(z), 6)
0.043845

s10z ¶

s10z(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Ten-point height (S10z).

(sum of the 5 highest peaks + sum of the absolute depths of the 5 deepest valleys) / 5, from find_peaks() and find_valleys() (fewer when there are fewer).

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `find_peaks()`).	required
`extent`		As in `find_peaks()`.	`None`
`crs`		As in `find_peaks()`.	`None`
`band`		As in `find_peaks()`.	`None`
`nodata`		As in `find_peaks()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> round(ls.surface.s10z(z), 6)
0.665448

Bearing area curve¶

bearing_area ¶

bearing_area(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> DataFrame

The bearing area (Abbott-Firestone) curve: the share of the surface at or above each height.

For every distinct height, from the highest, the fraction of cells whose height is at least that high (geodiv's curve is that of 1 - z, which puts rounding errors in the heights).

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (a single band; NaN or `nodata` = missing).	required
`extent`	`tuple of float`	`(xmin, xmax, ymin, ymax)` of an array (default: unit cells).	`None`
`crs`	`optional`	CRS of an array.	`None`
`band`	`int`	Band of a multi-band file or DataArray (1-based).	`None`
`nodata`	`float`	Value that marks missing cells.	`None`

Returns:

Type	Description
`DataFrame`	Columns `bearing` (0-1, increasing) and `height` (decreasing), one row per distinct height: the points of geodiv's `plot_ba_curve`.

Examples:

>>> import landschaft as ls
>>> curve = ls.surface.bearing_area(ls.surface.example_data("orelevation"))
>>> curve.iloc[[0, -1]]
      bearing  height
0    0.000015  1390.0
953  1.000000   433.0

height_ba ¶

height_ba(
    x,
    xval,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
)

Height of the bearing area curve at a bearing fraction.

The height above which a fraction xval of the surface lies: R's type-7 quantile of z at 1 - xval.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `bearing_area()`).	required
`xval`	`float or sequence of float`	Bearing fraction(s), 0 (the highest height) to 1 (the lowest).	required
`extent`		As in `bearing_area()`.	`None`
`crs`		As in `bearing_area()`.	`None`
`band`		As in `bearing_area()`.	`None`
`nodata`		As in `bearing_area()`.	`None`

Returns:

Type	Description
`float or ndarray`	One height per fraction.

Examples:

>>> import landschaft as ls
>>> dem = ls.surface.example_data("orelevation")
>>> ls.surface.height_ba(dem, [0, 0.5, 1]).tolist()
[1390.0, 884.0, 433.0]

sdc ¶

sdc(
    x,
    low: float,
    high: float,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> float

Areal material ratio height difference (Sdc).

The height difference between two bearing fractions: height_ba() at low minus at high.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `bearing_area()`).	required
`low`	`float`	Bearing fractions, 0 <= low < high <= 1.	required
`high`	`float`	Bearing fractions, 0 <= low < high <= 1.	required
`extent`		As in `bearing_area()`.	`None`
`crs`		As in `bearing_area()`.	`None`
`band`		As in `bearing_area()`.	`None`
`nodata`		As in `bearing_area()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> ls.surface.sdc(ls.surface.example_data("orelevation"), 0.1, 0.8)
460.0

sbi ¶

sbi(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Surface bearing index (Sbi): Sq divided by the height at 5 % bearing.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `bearing_area()`), with its mean or trend removed.	required
`extent`		As in `bearing_area()`.	`None`
`crs`		As in `bearing_area()`.	`None`
`band`		As in `bearing_area()`.	`None`
`nodata`		As in `bearing_area()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> round(ls.surface.sbi(ls.surface.example_data("normforest")), 6)
0.801444

svi ¶

svi(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Valley fluid retention index (Svi): the void volume of the valley zone.

The area between the bearing area curve and its top from 80 % to 100 % bearing (Simpson's rule on 500 intervals), in height units (geodiv's definition, corrected: its area depends on whether the curve dips below 0, which gives negative values on positive surfaces, and its Simpson's rule drops the last point).

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `bearing_area()`).	required
`extent`		As in `bearing_area()`.	`None`
`crs`		As in `bearing_area()`.	`None`
`band`		As in `bearing_area()`.	`None`
`nodata`		As in `bearing_area()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> round(ls.surface.svi(ls.surface.example_data("normforest")), 6)
0.012704

sci ¶

sci(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Core fluid retention index (Sci): the void volume of the core zone.

The area between the bearing area curve and its top from 5 % to 100 % bearing (1000 intervals) minus svi().

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `bearing_area()`).	required
`extent`		As in `bearing_area()`.	`None`
`crs`		As in `bearing_area()`.	`None`
`band`		As in `bearing_area()`.	`None`
`nodata`		As in `bearing_area()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> round(ls.surface.sci(ls.surface.example_data("normforest")), 6)
0.069002

sk ¶

sk(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Core roughness depth (Sk): the height range of the core line.

The core line is the least-squares line through the bearing area curve over the flattest 40 % of it (geodiv's find_flat: the curve's slopes by central differences at 10000 points, one-sided at the ends, their mean absolute value over each 40 % window); Sk is the difference of its heights at 0 % and 100 %.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `bearing_area()`).	required
`extent`		As in `bearing_area()`.	`None`
`crs`		As in `bearing_area()`.	`None`
`band`		As in `bearing_area()`.	`None`
`nodata`		As in `bearing_area()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> round(ls.surface.sk(ls.surface.example_data("normforest")), 6)
0.102369

spk ¶

spk(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Reduced peak height (Spk): the peaks above the core.

The height difference between the top of the bearing area curve and the curve at the bearing where the core line (sk()) starts.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `bearing_area()`).	required
`extent`		As in `bearing_area()`.	`None`
`crs`		As in `bearing_area()`.	`None`
`band`		As in `bearing_area()`.	`None`
`nodata`		As in `bearing_area()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> round(ls.surface.spk(ls.surface.example_data("normforest")), 6)
0.096705

svk ¶

svk(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Reduced valley depth (Svk): the valleys below the core.

The height difference between the bottom of the bearing area curve and the curve at the bearing where the core line (sk()) ends.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `bearing_area()`).	required
`extent`		As in `bearing_area()`.	`None`
`crs`		As in `bearing_area()`.	`None`
`band`		As in `bearing_area()`.	`None`
`nodata`		As in `bearing_area()`.	`None`

Returns:

Type	Description
`float`

Examples:

>>> import landschaft as ls
>>> round(ls.surface.svk(ls.surface.example_data("normforest")), 6)
0.544987

Autocorrelation, orientation and wavelength¶

aacf ¶

aacf(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> DataArray

Areal autocorrelation function of a surface.

As geodiv: the surface times a two-dimensional Hann window, its power spectrum by FFT, the inverse FFT of that, shifted to put zero lag in the centre and divided by its maximum. With missing cells, the largest square without them around the centre is used (and the image covers it).

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (a single band; NaN or `nodata` = missing).	required
`extent`	`tuple of float`	`(xmin, xmax, ymin, ymax)` of an array (default: unit cells).	`None`
`crs`	`optional`	CRS of an array.	`None`
`band`	`int`	Band of a multi-band file or DataArray (1-based).	`None`
`nodata`	`float`	Value that marks missing cells.	`None`

Returns:

Type	Description
`DataArray`	The autocorrelation at every lag, 1 at zero lag.

Examples:

>>> import landschaft as ls
>>> a = ls.surface.aacf(ls.surface.example_data("normforest"))
>>> a.shape, float(a.max())
((371, 371), 1.0)

scl ¶

scl(
    x,
    threshold=THRESHOLDS,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> ndarray

Correlation lengths: how far the autocorrelation takes to decay.

On the upper half of aacf() (with the zero-lag row), the distances (map units; metres on longitude/latitude rasters) from zero lag to the cells whose autocorrelation is at or below each threshold: the shortest (fastest decay, Sal for 0.2) and the longest (slowest decay).

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `aacf()`).	required
`threshold`	`float or sequence of float`	Autocorrelation value(s).	`(0.2, 1/e)`
`extent`		As in `aacf()`.	`None`
`crs`		As in `aacf()`.	`None`
`band`		As in `aacf()`.	`None`
`nodata`		As in `aacf()`.	`None`

Returns:

Type	Description
`ndarray`	The fastest decay length for each threshold, then the slowest ones.

Notes

Distances are exact in map units on projected rasters and haversine distances (terra's, radius 6378137 m) on longitude/latitude ones. geodiv measures from one row above zero lag (its upper half leaves that row out), with GDAL's proximity in single precision and the cell width in both directions, and fails with missing cells.

Examples:

>>> import landschaft as ls
>>> ls.surface.scl(ls.surface.example_data("orelevation")).round(1).tolist()
[4662.9, 3317.3, 6873.0, 6873.0]

stxr ¶

stxr(
    x,
    threshold=THRESHOLDS,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> ndarray

Texture aspect ratio (Str): fastest over slowest decay length.

From scl(), for each threshold: near 0 for a strongly oriented texture, near 1 for an isotropic one.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `aacf()`).	required
`threshold`	`float or sequence of float`	Autocorrelation value(s).	`(0.2, 1/e)`
`extent`		As in `aacf()`.	`None`
`crs`		As in `aacf()`.	`None`
`band`		As in `aacf()`.	`None`
`nodata`		As in `aacf()`.	`None`

Returns:

Type	Description
`ndarray`	One ratio per threshold.

Examples:

>>> import landschaft as ls
>>> ls.surface.stxr(ls.surface.example_data("orelevation")).round(4).tolist()
[0.6784, 0.4827]

std ¶

std(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> ndarray

Texture direction (Std) and texture direction index (Stdi).

From the amplitude of the surface's Fourier spectrum (upper half, within the largest circle around zero frequency that fits), summed over orientations 0 to 179.5 degrees in steps of 0.5 (+-0.25): Std is the angle with the largest sum, Stdi the mean sum over the largest (near 1: no dominant direction). Angles are those of the frequency vector: 0 for a texture that varies along x (stripes running north-south), 90 for one that varies along y.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `aacf()`).	required
`extent`		As in `aacf()`.	`None`
`crs`		As in `aacf()`.	`None`
`band`		As in `aacf()`.	`None`
`nodata`		As in `aacf()`.	`None`

Returns:

Type	Description
`ndarray`	`[std, stdi]` (degrees, ratio).

Notes

geodiv takes the angles from one row above and, on even widths, one column left of zero frequency, and includes zero frequency (the sum of the heights) in the 0 degree bin.

Examples:

>>> import landschaft as ls
>>> ls.surface.std(ls.surface.example_data("normforest")).round(4).tolist()
[90.0, 0.1937]

srw ¶

srw(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> ndarray

Dominant radial wavelength (Srw), its index (Srwi) and the half-power wavelength (Shw).

The amplitude of the surface's Fourier spectrum (upper half) summed over circles around zero frequency (up to 100, evenly spaced): Srw is the radius with the largest sum, Srwi the mean sum over the largest, Shw the radius that splits the summed amplitude in half (NaN on ties), measured from zero frequency (geodiv: from one row above it). Radii are distances in the spectrum in map units: a texture repeating k times across a raster of n cells has its peak at k cell sizes, a wavelength of n / k cells.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `aacf()`).	required
`extent`		As in `aacf()`.	`None`
`crs`		As in `aacf()`.	`None`
`band`		As in `aacf()`.	`None`
`nodata`		As in `aacf()`.	`None`

Returns:

Type	Description
`ndarray`	`[srw, srwi, shw]`.

Examples:

>>> import landschaft as ls
>>> ls.surface.srw(ls.surface.example_data("normforest")).round(4).tolist()
[409.6675, 0.5731, 1638.6698]

sfd ¶

sfd(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> float

Fractal dimension of the surface (Sfd), by triangular prisms.

Clarke's (1986) method as geodiv implements it: heights rescaled to 0-1, the area of triangular prisms over the central square of a power-of-two side at cell sizes 1, 2, 4...; 2 minus the slope of log(area) against log(cell area). geodiv's code has a typo in one prism side, rescales against 99999 / -99999 when the heights exceed them and reads past the raster (1-based indices): corrected here.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `aacf()`), without missing cells.	required
`extent`		As in `aacf()`.	`None`
`crs`		As in `aacf()`.	`None`
`band`		As in `aacf()`.	`None`
`nodata`		As in `aacf()`.	`None`

Returns:

Type	Description
`float`	NaN with missing cells.

Examples:

>>> import landschaft as ls
>>> round(ls.surface.sfd(ls.surface.example_data("orelevation")), 6)
2.000007

Texture images¶

window_metrics ¶

window_metrics(
    x,
    metrics: str | Iterable[str],
    window: int = 5,
    shape: str = "square",
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> Dataset

Surface metrics in a moving window: one texture image per metric.

Every cell gets the value of each metric over the window x window cells around it, computed exactly as metrics() computes it for a raster, on a raster that covers the window. Cells beyond the edges of the raster are missing in the windows that reach them.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (a single band; NaN or `nodata` = missing).	required
`metrics`	`str or list of str`	Names from `list_metrics()`.	required
`window`	`int`	Side of the window in cells, odd.	`5`
`shape`	`('square', 'circle')`	`"circle"`: only the cells within `window // 2` cells of the centre (the others are missing).	`"square"`
`extent`	`tuple of float`	`(xmin, xmax, ymin, ymax)` of an array (default: unit cells).	`None`
`crs`	`optional`	CRS of an array.	`None`
`band`	`int`	Band of a multi-band file or DataArray (1-based).	`None`
`nodata`	`float`	Value that marks missing cells.	`None`

Returns:

Type	Description
`Dataset`	One variable per metric on the grid of `x`.

Notes

Windows run in parallel (n_threads). sv, sph and sds on square windows take shortcuts (running extremes, a summed-area table of the peaks) that give the same values; the other metrics cost what they cost on a raster the size of the window, times the number of cells. The bearing-area metrics (sbi...svk) and the spatial ones (sal... sfd) are expensive per window (sk, spk and svk average 6000 runs of 4001 slopes each): use them on small rasters.

geodiv's own texture_image does not work with current terra (it reads the raster as one column), its focal_metrics transposes every window and its circular windows mask nothing (their radius is the window's side): here each window gets the metric computed on it as on a raster.

Examples:

>>> import landschaft as ls
>>> z = ls.surface.example_data("normforest")
>>> tex = ls.surface.window_metrics(z, ["sa", "sku", "sds"], window=7)
>>> list(tex.data_vars), tex["sa"].shape
(['sa', 'sku', 'sds'], (371, 371))

Terrain¶

Indices of a DEM's shape as the R package spatialEco computes them (over terra's slope and aspect), corrected where it errs.

terrain ¶

terrain(
    x,
    v: str = "slope",
    unit: str = "degrees",
    neighbors: int = 8,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Slope or aspect of a DEM, as terra's terrain.

Horn's method (neighbors=8: the 3 x 3 neighbourhood weighted 1-2-1) or the 4 direct neighbours; on longitude/latitude rasters with the cell sizes in metres (geodesics on WGS84) and elevation in metres. The edge cells are missing.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The DEM (a single band; NaN or `nodata` = missing).	required
`v`	`('slope', 'aspect')`	Slope (0 flat) or aspect (clockwise from north, where the slope faces).	`"slope"`
`unit`	`('degrees', 'radians')`		`"degrees"`
`neighbors`	`(8, 4)`		`8`
`extent`	`tuple of float`	`(xmin, xmax, ymin, ymax)` of an array (default: unit cells).	`None`
`crs`	`optional`	CRS of an array.	`None`
`band`	`int`	Band of a multi-band file or DataArray (1-based).	`None`
`nodata`	`float`	Value that marks missing cells.	`None`

Returns:

Type	Description
`DataArray`

Examples:

>>> import landschaft as ls
>>> dem = ls.surface.example_data("orelevation")
>>> slope = ls.surface.terrain(dem, "slope")
>>> round(float(slope.mean()), 2)
18.45

tpi ¶

tpi(
    x,
    scale=3,
    win: str = "rectangle",
    normalize: bool = False,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Topographic position index: the elevation minus the mean of the window.

Positive on ridges and hilltops, negative in valleys, near 0 on flats and even slopes. The window includes the cell itself; cells whose window reaches a missing cell or the edge are missing.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The DEM (see `terrain()`).	required
`scale`	`int or float`	`"rectangle"`: the window's side in cells (odd); `"circle"`: its radius in map units (at least two cells).	`3`
`win`	`('rectangle', 'circle')`		`"rectangle"`
`normalize`	`bool`	Divided by the window's standard deviation.	`False`
`extent`		As in `terrain()`.	`None`
`crs`		As in `terrain()`.	`None`
`band`		As in `terrain()`.	`None`
`nodata`		As in `terrain()`.	`None`

Returns:

Type	Description
`DataArray`

Notes

The window means (and standard deviations) come from compensated prefix sums along the rows, one step per window row: within about 1e-13 of the values' range of summing each window.

spatialEco's zero.correct option computes the same mean (and, with normalize, a standard deviation of 0 by a misplaced parenthesis): it is not offered.

Examples:

>>> import landschaft as ls
>>> dem = ls.surface.example_data("orelevation")
>>> t = ls.surface.tpi(dem, 5)
>>> round(float(t.min()), 1), round(float(t.max()), 1)
(-31.2, 34.4)

tri ¶

tri(
    x,
    s=3,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Terrain ruggedness index (Riley et al. 1999).

sqrt(sum((z - z_i)^2)) over the other cells of the window that have a value (so the edges get one too).

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The DEM (see `terrain()`).	required
`s`	`int or (int, int)`	The window in cells (odd).	`3`
`extent`		As in `terrain()`.	`None`
`crs`		As in `terrain()`.	`None`
`band`		As in `terrain()`.	`None`
`nodata`		As in `terrain()`.	`None`

Returns:

Type	Description
`DataArray`

Examples:

>>> import landschaft as ls
>>> t = ls.surface.tri(ls.surface.example_data("orelevation"))
>>> round(float(t.mean()), 2)
23.15

vrm ¶

vrm(
    x,
    s=3,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Vector ruggedness measure (Sappington et al. 2007).

1 minus the length of the mean of the unit normal vectors of the window's cells: 0 for a plane however steep, towards 1 for rugged terrain. Normals from terrain()'s slope and aspect.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The DEM (see `terrain()`).	required
`s`	`int or (int, int)`	The window in cells (odd).	`3`
`extent`		As in `terrain()`.	`None`
`crs`		As in `terrain()`.	`None`
`band`		As in `terrain()`.	`None`
`nodata`		As in `terrain()`.	`None`

Returns:

Type	Description
`DataArray`

Examples:

>>> import landschaft as ls
>>> v = ls.surface.vrm(ls.surface.example_data("orelevation"))
>>> round(float(v.mean()), 4)
0.006

srr ¶

srr(
    x,
    s=5,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Surface relief ratio (Pike & Wilson 1971): (mean - min) / (max - min) of the window.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The DEM (see `terrain()`).	required
`s`	`int or (int, int)`	The window in cells (odd).	`5`
`extent`		As in `terrain()`.	`None`
`crs`		As in `terrain()`.	`None`
`band`		As in `terrain()`.	`None`
`nodata`		As in `terrain()`.	`None`

Returns:

Type	Description
`DataArray`

Examples:

>>> import landschaft as ls
>>> r = ls.surface.srr(ls.surface.example_data("orelevation"))
>>> round(float(r.mean()), 3)
0.501

dissection ¶

dissection(
    x,
    s=5,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Martonne's modified dissection: (z - min) / (max - min) of the window.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The DEM (see `terrain()`).	required
`s`	`int or (int, int)`	The window in cells (odd).	`5`
`extent`		As in `terrain()`.	`None`
`crs`		As in `terrain()`.	`None`
`band`		As in `terrain()`.	`None`
`nodata`		As in `terrain()`.	`None`

Returns:

Type	Description
`DataArray`

Examples:

>>> import landschaft as ls
>>> d = ls.surface.dissection(ls.surface.example_data("orelevation"))
>>> round(float(d.mean()), 3)
0.502

curvature ¶

curvature(
    x,
    type: str = "planform",
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Surface curvature of a DEM.

"planform", "profile" and "total" (their sum) from Zevenbergen & Thorne's (1987) partial derivatives of the 3 x 3 neighbourhood, p, q (first) and r, s, t (second): planform -(q^2 r - 2pqs + p^2 t) / ((p^2 + q^2) sqrt(1 + p^2 + q^2)), profile -(p^2 r + 2pqs + q^2 t) / ((p^2 + q^2) (1 + p^2 + q^2)^(3/2)); "mcnab": the mean of the differences between the cell and its 8 neighbours, divided by its largest value; "bolstad": 10000 (z - mean of the 3 x 3) / 1000 / (1.5 cell width).

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The DEM, on a projected raster (see `terrain()`).	required
`type`	`('planform', 'profile', 'total', 'mcnab', 'bolstad')`		`"planform"`
`extent`		As in `terrain()`.	`None`
`crs`		As in `terrain()`.	`None`
`band`		As in `terrain()`.	`None`
`nodata`		As in `terrain()`.	`None`

Returns:

Type	Description
`DataArray`	Missing on the edges (and flat cells for the Zevenbergen & Thorne curvatures).

Notes

spatialEco halves the second derivatives (divides by 2 L^2), pads the edges with 0, rounds to 6 decimals, uses the cell width for both directions, divides McNab's sum by 9 and subtracts 1 (operator precedence), and fails for "bolstad": corrected here.

Examples:

>>> import landschaft as ls
>>> c = ls.surface.curvature(ls.surface.example_data("normforest"), "total")
>>> c.shape
(371, 371)

spherical_sd ¶

spherical_sd(
    x,
    d=3,
    variance: bool = False,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Spherical standard deviation of the slope and aspect (terrain heterogeneity).

The unit normal vectors of the cells (from terrain()) averaged over the window: sqrt(1 - R^2), R the length of the mean vector (variance: 1 - R^2).

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The DEM (see `terrain()`).	required
`d`	`int or (int, int)`	The window in cells (odd).	`3`
`variance`	`bool`		`False`
`extent`		As in `terrain()`.	`None`
`crs`		As in `terrain()`.	`None`
`band`		As in `terrain()`.	`None`
`nodata`		As in `terrain()`.	`None`

Returns:

Type	Description
`DataArray`

Examples:

>>> import landschaft as ls
>>> v = ls.surface.spherical_sd(ls.surface.example_data("orelevation"))
>>> round(float(v.mean()), 4)
0.0978

hsp ¶

hsp(
    x,
    min_scale=3,
    max_scale=27,
    inc=4,
    win: str = "rectangle",
    normalize: bool = False,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Hierarchical slope position (Murphy et al. 2010).

The sum over the scales min_scale, min_scale + inc, ... max_scale of 100 times the z-score of tpi() at that scale ((tpi - mean) / sd over the raster).

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The DEM (see `terrain()`).	required
`min_scale`	`default 3, 27, 4`	The scales: odd window sides in cells, or radii in map units for circles.	`3`
`max_scale`	`default 3, 27, 4`	The scales: odd window sides in cells, or radii in map units for circles.	`3`
`inc`	`default 3, 27, 4`	The scales: odd window sides in cells, or radii in map units for circles.	`3`
`win`	`('rectangle', 'circle')`		`"rectangle"`
`normalize`	`bool`	Rescaled to 0-1.	`False`
`extent`		As in `terrain()`.	`None`
`crs`		As in `terrain()`.	`None`
`band`		As in `terrain()`.	`None`
`nodata`		As in `terrain()`.	`None`

Returns:

Type	Description
`DataArray`

Notes

spatialEco computes tpi - mean / sd (operator precedence), not the z-score: corrected here.

Examples:

>>> import landschaft as ls
>>> h = ls.surface.hsp(ls.surface.example_data("orelevation"), 3, 11, 4, normalize=True)
>>> round(float(h.min()), 1), round(float(h.max()), 1)
(0.0, 1.0)

trasp ¶

trasp(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> DataArray

Transformed aspect (Roberts & Cooper 1989): (1 - cos(aspect - 30 deg)) / 2.

0 on north-northeast slopes (cool, wet), 1 on south-southwest ones (hot, dry). spatialEco uses 3.142 for pi.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The DEM (see `terrain()`).	required
`extent`		As in `terrain()`.	`None`
`crs`		As in `terrain()`.	`None`
`band`		As in `terrain()`.	`None`
`nodata`		As in `terrain()`.	`None`

Returns:

Type	Description
`DataArray`

Examples:

>>> import landschaft as ls
>>> t = ls.surface.trasp(ls.surface.example_data("orelevation"))
>>> round(float(t.min()), 3), round(float(t.max()), 3)
(0.0, 1.0)

hli ¶

hli(
    aspect,
    slope,
    latitude=None,
    direct: bool = False,
    scaled: bool = False,
    units: str = "degrees",
    hemisphere: str = "northern",
    force_hemisphere: bool = True,
    equation: int = 1,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Heat load index (McCune & Keon 2002).

Their equations 1-3 of potential direct incident radiation from slope, folded aspect and latitude: the aspect folded about the northeast-southwest line (hottest on southwest slopes; southern hemisphere: about northwest- southeast) or, with direct, about north-south. The equations give the logarithm of the radiation (MJ cm^-2 yr^-1); scaled gives the radiation.

Parameters:

Name	Type	Description	Default
`aspect`	`(str, path, DataArray or array_like)`	Aspect and slope rasters (see `terrain()`).	required
`slope`	`(str, path, DataArray or array_like)`	Aspect and slope rasters (see `terrain()`).	required
`latitude`	`float`	Latitude (degrees). Default: the centre of the raster (which needs a CRS).	`None`
`direct`	`bool`	Fold about north-south (direct radiation) instead of NE-SW.	`False`
`scaled`	`bool`	exp() of the equation.	`False`
`units`	`('degrees', 'radians')`	Units of aspect and slope.	`"degrees"`
`hemisphere`	`('northern', 'southern')`		`"northern"`
`force_hemisphere`	`bool`	Use `hemisphere`; otherwise the sign of the latitude decides.	`True`
`equation`	`(1, 2, 3)`	McCune & Keon's equation (3: latitudes 30-60 degrees).	`1`
`extent`		As in `terrain()` (for arrays).	`None`
`crs`		As in `terrain()` (for arrays).	`None`
`band`		As in `terrain()` (for arrays).	`None`
`nodata`		As in `terrain()` (for arrays).	`None`

Returns:

Type	Description
`DataArray`

Notes

spatialEco's southern hemisphere folds the aspect in degrees without the absolute values (then takes its cosine as radians), uses the signed latitude with the mirrored aspect, and with units="radians" converts the aspect again: corrected here (the southern hemisphere uses the absolute latitude).

Examples:

>>> import landschaft as ls
>>> dem = ls.surface.example_data("orelevation")
>>> h = ls.surface.hli(ls.surface.terrain(dem, "aspect"), ls.surface.terrain(dem, "slope"))
>>> round(float(h.mean()), 3)
-0.288

Landforms (geomorphons)¶

The 10 landforms and 498 ternary patterns of Jasiewicz & Stepinski (2013), from the lines of sight of every cell: GRASS GIS's r.geomorphon, bit for bit (four of its errors corrected).

geomorphons ¶

geomorphons(
    dem,
    *,
    search: float = 3,
    skip: float = 0,
    flat: float = 1.0,
    dist: float = 0.0,
    units: str = "cells",
    comparison: str = "anglev1",
    extended: bool = False,
    measures: bool = False,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    _grass_errors: bool = False,
    n_threads: int | None = None,
) -> Geomorphons

Landforms of a DEM by geomorphons (Jasiewicz & Stepinski 2013).

From each cell, along the 8 directions, the terrain is looked at from skip + 1 cells away up to the search radius: the highest angle above the horizon (zenith) and the lowest (nadir). A direction is higher (+) when its zenith angle wins, lower (-) when the nadir angle does, level when neither exceeds the flatness threshold flat. The numbers of + and - name the landform (flat, peak, ridge, shoulder, spur, slope, hollow, footslope, valley, pit), and the 8 signs, up to rotation and reflection, one of 498 patterns. Unlike slope and curvature, which see a 3 x 3 window, the forms follow the terrain at the scale of the search radius: a better covariate of the site's shape for land change, species distributions and soils.

The same forms, patterns and measures as GRASS GIS's r.geomorphon (8.4), bit for bit, but for four corrections (see Notes).

Parameters:

Name	Type	Description	Default
`dem`	`(str, path, DataArray or array_like)`	The DEM (projected; a single band; NaN or `nodata` = missing). It is read as float32, as GRASS reads it.	required
`search`	`float`	The outer search radius, in cells (or map units with `units="map"`): `r.geomorphon`'s `search`.	`3`
`skip`	`float`	The inner radius not looked at (`skip`), less than `search`.	`0`
`flat`	`float`	The flatness threshold, in degrees.	`1`
`dist`	`float`	The flatness distance (`dist`): beyond it the threshold falls as the angle of the height `tan(flat) x dist` (0: none; it must lie between `skip` and `search`).	`0`
`units`	`('cells', 'map')`	The unit of `search`, `skip` and `dist` (`"map"`: map units, `r.geomorphon -m`).	`"cells"`
`comparison`	`('anglev1', 'anglev2', 'anglev2_distance')`	How the zenith and nadir angles decide a direction (`r.geomorphon`'s `comparison`).	`"anglev1"`
`extended`	`bool`	The extended correction (`-e`): with a search radius above 10 cells, a ridge, peak, shoulder or footslope is looked at again within half the radius.	`False`
`measures`	`bool`	Also the measures of the geomorphon (intensity, exposition, range, variance, elongation, azimuth, extend, width).	`False`
`extent`		As in `terrain()`.	`None`
`crs`		As in `terrain()`.	`None`
`band`		As in `terrain()`.	`None`
`nodata`		As in `terrain()`.	`None`

Returns:

Type	Description
`Geomorphons`

Notes

Four things r.geomorphon 8.4 does not do as intended, here in the correct form: its lines of sight stop one row short of the DEM's last row for the cells of its bottom search rows (a buffer of rows that holds the last one out of reach); its binary positive and negative patterns add the direction numbers instead of setting their bits; its extended correction keeps the inner radius in the shorter, second search, though it sets it to 0; and the shape of a geomorphon whose first direction (NE) has no point is all zeros. Elsewhere the values are GRASS's, compared on synthetic and real DEMs.

Longitude / latitude DEMs are refused: GRASS measures their distances on the ellipsoid; reproject the DEM first.

References

Jasiewicz, J. & Stepinski, T. F. (2013). Geomorphons, a pattern recognition approach to classification and mapping of landforms. Geomorphology 182, 147-156.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> yy, xx = np.mgrid[0:60, 0:60] * 10.0
>>> hill = 300 * np.exp(-((yy - 300) ** 2 + (xx - 300) ** 2) / 2e4)   # a round hill
>>> g = ls.surface.geomorphons(hill, search=10, extent=(0, 600, 0, 600))
>>> ls.surface.GEOMORPHON_FORMS[int(g.forms[30, 30])]                 # its summit
'peak'
>>> g                                                                  # the hill's flanks: slopes and spurs
Geomorphons(60 x 60; slope 61%, spur 27%, footslope 8%)
>>> g.share().loc[[2, 8], ["name", "cells"]]
           name  cells
form
2          peak     21
8     footslope    274

Geomorphons `dataclass` ¶

Geomorphons(
    forms: DataArray,
    ternary: DataArray,
    positive: DataArray,
    negative: DataArray,
    intensity: DataArray | None = None,
    exposition: DataArray | None = None,
    range: DataArray | None = None,
    variance: DataArray | None = None,
    elongation: DataArray | None = None,
    azimuth: DataArray | None = None,
    extend: DataArray | None = None,
    width: DataArray | None = None,
)

Bases: Saveable

The geomorphons of a DEM (geomorphons()).

Attributes:

Name	Type	Description
`forms`	`DataArray`	The landform of each cell, 1 to 10 (`GEOMORPHON_FORMS`: flat, peak, ridge, shoulder, spur, slope, hollow, footslope, valley, pit); NaN where the DEM has no value and on its outer `skip + 1` cells.
`ternary`	`DataArray`	The ternary pattern: the 8 signs (higher, level, lower) of the directions as a base-3 number, the smallest of its rotations and reflections (498 distinct codes, 0 to 6560).
`positive, negative`	`DataArray`	The directions whose zenith angle exceeds the flatness threshold (whose nadir angle is below minus the threshold), as 8 bits, the smallest of their rotations (0 to 255).
`intensity, exposition, range, variance, elongation, azimuth, extend, width`	`DataArray or None`	With `measures=True`: the mean height of the cell above the points that decide its directions (`intensity`), minus the height of the farthest of them in elevation (`exposition`), their range and variance of elevation; the shape of the octagon they draw on the plane: its elongation, the azimuth of its long axis (degrees), its area relative to the octagon of the search radius (`extend`) and its width (map units).

share ¶

share() -> DataFrame

Cells and share of each landform: form (code), name, cells, share.

ls.surface.GEOMORPHON_FORMS maps the codes to the names: 1 flat, 2 peak, 3 ridge, 4 shoulder, 5 spur, 6 slope, 7 hollow, 8 footslope, 9 valley, 10 pit.

Hilltops (APP de topo de morro)¶

hilltop ¶

hilltop(
    dem,
    *,
    min_height: float = 100.0,
    min_slope: float = 25.0,
    slope: str = "mean",
    fraction: float = 2 / 3,
    plain=None,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> Hilltop

Hilltops as permanent protection areas (APP de topo de morro).

Brazil's Forest Code (Lei 12.651/2012, art. 4, IX) protects the tops of hills, mounts, mountains and ranges at least 100 m high and of mean slope above 25 degrees: the areas above the contour at 2/3 of the height, measured from the base, which is the plane of the adjacent plain or water body or, in undulating relief, the elevation of the nearest saddle point.

Every summit of the DEM is found with its base by sweeping the cells from the highest down and joining each to the regions of its 8 neighbours (a merge tree): where regions meet, the cell is a saddle, and the regions of lower summits end there (their base is the saddle: the key col of topographic prominence, so a range is one elevation above its lowest saddle and its peaks are hills of their own above theirs). A region that reaches a cell of plain gets that cell's elevation as its base; the highest summit of a connected area that meets neither gets its lowest elevation. A hill qualifies when its height is at least min_height and its slope above min_slope; its APP is the area joined to its summit through cells at least base + fraction x height high.

Parameters:

Name	Type	Description	Default
`dem`	`(str, path, DataArray or array_like)`	The DEM (metres), projected or in longitude / latitude.	required
`min_height`	`float`	The least height of a hill (metres; CONAMA 303/2002: 50).	`100`
`min_slope`	`float`	The slope a hill must exceed (degrees; CONAMA 303/2002: about 16.7, 30 %, along the line of steepest slope).	`25`
`slope`	`('mean', 'line')`	`"mean"`: the mean slope (Horn, degrees) of the hill's cells above its base (the law's "inclinação média"); `"line"`: the slope of the line from the summit to the base cell.	`"mean"`
`fraction`	`float`	Where the APP starts, as a fraction of the height above the base.	`2/3`
`plain`	`DataArray or array_like`	Plains and water bodies (true / non-zero): a hill that reaches them takes their elevation as its base.	`None`
`extent`		As in `terrain()`.	`None`
`crs`		As in `terrain()`.	`None`
`band`		As in `terrain()`.	`None`
`nodata`		As in `terrain()`.	`None`

Returns:

Type	Description
`Hilltop`

Notes

The law leaves room for interpretation and these are the choices made: the nearest saddle of an elevation is the col where its region meets higher terrain (the hill's prominence), so the highest peak of a range is measured from the range's base and lower peaks from their saddles; the mean slope is that of the cells of the hill above its base (its region in the merge tree); the APP is the part of that region above the threshold that is joined to the summit (a separate peak above the threshold, joined below it, is judged as a hill of its own). Plains are not inferred from the slope: give plain (flat floodplains, lakes, rivers, the sea). Cells on the DEM's edges have no slope and are left out of the mean. Results depend on the DEM's resolution and accuracy, and automatic delineation is a starting point for field and legal work.

Dinamica EGO's CalcHilltop follows CONAMA 303/2002 on a DEM cut into elevation bands: bases are the band of a local minimum reached downhill, and its slope criterion is the angle between the summit and the nearest cell of the base band that borders other features or the map's edge (5 degrees for an isolated cone of 21.8 degrees in a plain, measured with Dinamica EGO 8.13), so its results differ.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> yy, xx = np.mgrid[0:100, 0:100] * 10.0
>>> cone = 500 + np.clip(300 - 1.0 * np.hypot(yy - 500, xx - 500), 0, None)  # 300 m high, 45 degrees
>>> h = ls.surface.hilltop(cone, extent=(0, 1000, 0, 1000))
>>> top = h.hills.iloc[0]
>>> float(top.base), float(top.height), round(float(top.mean_slope)), bool(top.qualifies)
(500.0, 300.0, 44, True)
>>> int(np.nansum(h.app)) == int(np.sum(cone >= 700))   # the APP: above 500 + 2/3 x 300
True

Hilltop `dataclass` ¶

Hilltop(app: DataArray, hill: DataArray, hills: DataFrame)

Bases: Saveable

The hilltops of a DEM (hilltop()).

Attributes:

Name	Type	Description
`app`	`DataArray`	1 in the upper part of the hills that qualify (the APP), 0 elsewhere on the DEM, NaN where it has no value.
`hill`	`DataArray`	The id (`hills` row) of the hill whose upper part holds each cell (of the highest threshold, where upper parts nest), NaN elsewhere.
`hills`	`DataFrame`	One row per summit, highest first: `x`, `y` (the summit), `top`, `base`, `base_kind` (`"saddle"`: the key col to higher terrain; `"plain"`: the plain or water body it meets; `"lowest"`: the lowest cell of its connected area), `base_x`, `base_y` (the saddle, the plain cell met, or the lowest cell nearest the summit), `height`, `mean_slope` (degrees, of its cells above the base), `line_slope` (degrees, of the line from the summit to its foot: the nearest cell at or below its base), `area` (above the base, map units squared), `area_ha` (hectares; at the centre latitude for longitude / latitude), `qualifies`, `threshold` (the elevation from which the APP starts) and `app_cells`.

Spatial interaction¶

gravity ¶

gravity(
    sources,
    exponent: float = 2.0,
    *,
    mask=None,
    unit: str = "cells",
    method: str = "auto",
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Gravitational interaction potential of source cells.

Every cell gets the sum over the source cells of their mass over the distance to them raised to exponent: P_i = sum_k m_k / d_ik^beta, the gravity model of spatial interaction (a "pressure" of towns, markets or roads to use as a driver of land change). The sources are the cells of sources with a value (other than 0); the value is the mass. A source's own cell is at a distance of one cell. This is Dinamica EGO's CalcInteractionPotentialMap.

Parameters:

Name	Type	Description	Default
`sources`	`(str, path, DataArray or array_like)`	The masses (NaN or `nodata`: not a source).	required
`exponent`	`float`	The distance decay exponent beta (0: no decay).	`2`
`mask`	`(str, path, DataArray or array_like)`	A raster on the same grid; the potential is missing where it is missing. The sources under its missing cells still count.	`None`
`unit`	`('cells', 'map')`	Distances in cells (Dinamica) or in map units (cell width along rows, height along columns; a source's own cell at the square root of the cell area).	`"cells"`
`method`	`('auto', 'direct', 'fft')`	`"direct"`: the sum over every source-cell pair, in a fixed order (the same bits with any thread count); `"fft"`: the same sum as a convolution by FFT (rounding errors of about 1e-15 of the largest term, much faster with many sources); `"auto"`: direct up to 2000 sources (or 1e8 source-cell pairs).	`"auto"`
`extent`		As in `terrain()`.	`None`
`crs`		As in `terrain()`.	`None`
`band`		As in `terrain()`.	`None`
`nodata`		As in `terrain()`.	`None`

Returns:

Type	Description
`DataArray`	The potential (float64), 0 where no source reaches (no sources).

Notes

Dinamica computes each term and their sum in float32 (relative errors around 1e-7, growing with the number of sources) and scales the potential to the range of an integer cell type when asked for one; here the sum is in float64 and never scaled. Its distances are always in cells.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> towns = np.full((3, 4), np.nan)
>>> towns[0, 0], towns[2, 3] = 100, 50
>>> p = ls.surface.gravity(towns)
>>> round(float(p[0, 0]), 4), round(float(p[1, 1]), 4)
(103.8462, 60.0)

Trend surfaces¶

fit_plane ¶

fit_plane(
    x,
    order: int,
    *,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

The least-squares polynomial surface of a given order.

A polynomial in the cell-centre coordinates, x^i y^j with i + j <= order, fitted to the cells with a value (coordinates scaled to [-1, 1] and a Householder QR, as spatial::surf.ls) and evaluated in every cell.

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (a single band; NaN or `nodata` = missing).	required
`order`	`int`	Order of the polynomial, 0 (the mean) to 6.	required
`extent`	`tuple of float`	`(xmin, xmax, ymin, ymax)` of an array (default: unit cells).	`None`
`crs`	`optional`	CRS of an array.	`None`
`band`	`int`	Band of a multi-band file or DataArray (1-based).	`None`
`nodata`	`float`	Value that marks missing cells.	`None`

Returns:

Type	Description
`DataArray`	The fitted surface on the grid of `x` (also where `x` is missing).

Examples:

>>> import landschaft as ls
>>> ndvi = ls.surface.example_data("orforest")
>>> trend = ls.surface.fit_plane(ndvi, 1)
>>> round(float(trend.min()), 4), round(float(trend.max()), 4)
(0.6652, 0.6979)

best_fit_plane ¶

best_fit_plane(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> DataArray

The polynomial surface of order 0 to 3 that fits a surface best.

Of the four least-squares surfaces (fit_plane()), the one with the least Bayesian information criterion, N log(RSS / N) + p log(N) (N cells with a value, p coefficients), the lowest order on ties. The order is in attrs["order"].

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `fit_plane()`).	required
`extent`		As in `fit_plane()`.	`None`
`crs`		As in `fit_plane()`.	`None`
`band`		As in `fit_plane()`.	`None`
`nodata`		As in `fit_plane()`.	`None`

Returns:

Type	Description
`DataArray`

Notes

geodiv's bestfitplane takes the smallest absolute mean residual, but every least-squares fit with a constant term has residuals that average to zero: its choice rests on rounding errors (about 1e-17). To remove a trend of a chosen order, subtract fit_plane().

Examples:

>>> import landschaft as ls
>>> ndvi = ls.surface.example_data("orforest")
>>> ls.surface.best_fit_plane(ndvi).attrs["order"]
2

remove_plane ¶

remove_plane(
    x, *, extent=None, crs=None, band=None, nodata=None, n_threads: int | None = None
) -> DataArray

A surface minus its best-fitting polynomial surface (geodiv's remove_plane).

The residuals from best_fit_plane(): positive above the trend, negative below; missing cells stay missing. The order removed is in attrs["order"].

Parameters:

Name	Type	Description	Default
`x`	`(str, path, DataArray or array_like)`	The surface (see `fit_plane()`).	required
`extent`		As in `fit_plane()`.	`None`
`crs`		As in `fit_plane()`.	`None`
`band`		As in `fit_plane()`.	`None`
`nodata`		As in `fit_plane()`.	`None`

Returns:

Type	Description
`DataArray`

Examples:

>>> import landschaft as ls
>>> ndvi = ls.surface.example_data("orforest")
>>> flat = ls.surface.remove_plane(ndvi)
>>> flat.attrs["order"], round(float(flat.mean()), 12)
(2, 0.0)

Example data¶

example_data ¶

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

Path of an example surface.

Parameters:

Name	Type	Description	Default
`name`	`('normforest', 'orforest', 'orelevation')`	Three 371 x 371 rasters of about 30 m cells (longitude/latitude) over a part of south-western Oregon, from the geodiv package: `"orforest"`, the mean NDVI of June-August 2000 (Landsat path 45, row 30); `"normforest"`, the same with its best-fitting trend surface removed (the input the surface metrics expect); `"orelevation"`, SRTM elevation (m).	`"normforest"`

Returns:

Type	Description
`Path`

Examples:

>>> import landschaft as ls
>>> ndvi = ls.load_raster(ls.surface.example_data("orforest"))
>>> ndvi.shape, round(float(ndvi.mean()), 3)
((371, 371), 0.682)

Surface metrics¶

All metrics at once¶

metrics ¶

list_metrics ¶

Amplitude¶

sa ¶

sq ¶

ssk ¶

sku ¶

sv ¶

sph ¶

smean ¶

Slopes and developed area¶

sdq ¶

sdq6 ¶

sdr ¶

surface_area ¶

flat_area ¶

Peaks and valleys¶

find_peaks ¶

find_valleys ¶

sds ¶

ssc ¶

s10z ¶

Bearing area curve¶

bearing_area ¶

height_ba ¶

sdc ¶

sbi ¶

svi ¶

sci ¶

sk ¶

spk ¶

svk ¶

Autocorrelation, orientation and wavelength¶

aacf ¶

scl ¶

stxr ¶

std ¶

srw ¶

sfd ¶

Texture images¶

window_metrics ¶

Terrain¶

terrain ¶

tpi ¶

tri ¶

vrm ¶

srr ¶

dissection ¶

curvature ¶

spherical_sd ¶

hsp ¶

trasp ¶

hli ¶

Landforms (geomorphons)¶

geomorphons ¶

Geomorphons dataclass ¶

share ¶

Hilltops (APP de topo de morro)¶

hilltop ¶

Hilltop dataclass ¶

Spatial interaction¶

gravity ¶

Trend surfaces¶

fit_plane ¶

best_fit_plane ¶

remove_plane ¶

Example data¶

example_data ¶

Geomorphons `dataclass` ¶

Hilltop `dataclass` ¶