Terrain hydrology¶

Where water goes on a terrain: depressions filled so that every cell drains, the direction each cell drains to, how much area drains through it, the stream network, the height above the nearest drainage (HAND) and the wetness, stream power and sediment transport indices. They live in the ls.hydro module.

The directions and accumulations are those of WhiteboxTools 2.4 cell for cell (D8Pointer, DInfPointer, D8FlowAccumulation, DInfFlowAccumulation, FD8FlowAccumulation, ElevationAboveStream); the filling is its FillDepressions with water let out at the map's edge and at missing cells, which WhiteboxTools does not do, and the breaching its BreachDepressionsLeastCost with the channels that reach the edge cut. Flat areas drain by the method of Barnes, Lehman & Mulla (2014). The "mfd" method and distance_to_channel are the routing of the InVEST models (pygeoprocessing 2.4's flow_dir_mfd, flow_accumulation_mfd, distance_to_channel_mfd), cell for cell; its fill_pits gives the same DEM as fill.

import landschaft as ls

dem = ls.load_raster("srtm.tif")
filled = ls.hydro.fill(dem)
acc = ls.hydro.flow_accumulation(filled)                 # D8, in cells
rivers = ls.hydro.streams(acc, 1000)
height = ls.hydro.hand(filled, rivers)                   # height above the drainage
wetness = ls.hydro.twi(filled)                           # FD8 contributing area

Conditioning¶

fill ¶

fill(
    dem,
    *,
    drain_flats: bool = False,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Fill the depressions of a DEM, so that water can flow from every cell to the edge.

Every depression (a cell or group of cells lower than all its surroundings) is raised to the elevation where it would spill over: the lowest DEM that keeps the original elevations everywhere else and drains to the map's edge or to a missing cell. The Priority-Flood algorithm of Barnes, Lehman & Mulla (2014), the same DEM as WhiteboxTools' FillDepressions.

Parameters:

Name	Type	Description	Default
`dem`		The elevations (DataArray, array or file); NaN or `nodata` mark missing cells.	required
`drain_flats`	`bool`	Raise each filled cell by the smallest representable step above the cell it drains to (Priority-Flood+epsilon), so the filled areas slope towards their outlet. The steps are tiny: keep the result in float64. Without it, the filled areas are flat and `flow_direction()` drains them by the method of Barnes et al. (2014), which gives more natural flow paths.	`False`
`extent`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`crs`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`band`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`nodata`		For arrays without georeferencing, or to pick a band and a missing value.	`None`

Returns:

Type	Description
`DataArray`	The filled DEM, float64, on the DEM's grid.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> dem = np.array([[5, 5, 5, 5], [5, 1, 2, 5], [5, 2, 3, 4], [5, 5, 5, 5]], float)
>>> ls.hydro.fill(dem).values
array([[5., 5., 5., 5.],
       [5., 4., 4., 5.],
       [5., 4., 4., 4.],
       [5., 5., 5., 5.]])

breach ¶

breach(
    dem,
    *,
    max_length: int = 100,
    max_cost: float | None = None,
    min_dist: bool = True,
    flat_increment: float | None = None,
    fill_rest: bool = True,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Remove the depressions of a DEM by cutting channels through their rims.

Breaching often changes the terrain less than filling (fill()): a road embankment across a valley or a false rim of a few cells is cut, not the whole valley flooded behind it (Lindsay 2016). Each pit (a cell with no lower neighbour), from the lowest, is first raised to just below its lowest neighbour, then joined to the nearest cell it can drain to (lower than the channel, missing, or the map's edge) by the channel of least cost, cut to descend by flat_increment per cell. WhiteboxTools' BreachDepressionsLeastCost (after Lindsay & Dhun 2015), cell for cell, except that a channel that reaches the map's edge is cut (WhiteboxTools leaves it uncut).

Parameters:

Name	Type	Description	Default
`dem`		The elevations; NaN or `nodata` mark missing cells.	required
`max_length`	`int`	The longest channel, in cells.	`100`
`max_cost`	`float \| None`	The costliest channel to cut (`None`: any); pits that would need more are left (and filled with `fill_rest`).	`None`
`min_dist`	`bool`	Count the channel's length in its cost (the lowering times the distance), so shorter channels are preferred; otherwise only the lowering counts.	`True`
`flat_increment`	`float \| None`	The drop per cell along a channel (default: from the elevations and the cell size, as WhiteboxTools chooses it). The steps are small: keep the result in float64.	`None`
`fill_rest`	`bool`	Fill the depressions that no channel could solve (`fill()`).	`True`
`extent`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`crs`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`band`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`nodata`		For arrays without georeferencing, or to pick a band and a missing value.	`None`

Returns:

Type	Description
`DataArray`	The breached DEM, float64.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> dem = np.array([[9, 9, 9, 9, 9], [9, 5, 5, 5, 9], [9, 5, 1, 6, 4], [9, 5, 5, 5, 9], [9, 9, 9, 9, 9]], float)
>>> b = ls.hydro.breach(dem, flat_increment=0.01)
>>> float(b.values[2, 2])                       # the pit, raised below its lowest neighbour
4.99
>>> int((b.values != dem).sum())                # and one rim cell cut on the way to the 4 (east)
2
>>> bool(ls.hydro.flow_direction(b, flats=False).values[2, 2] > 0)
True

Flow¶

flow_direction ¶

flow_direction(
    dem,
    *,
    method: str = "d8",
    flats: bool = True,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

The direction each cell drains to.

Parameters:

Name	Type	Description	Default
`dem`		The elevations, best with its depressions filled (`fill()`).	required
`method`	`('d8', 'dinf', 'mfd')`	`"d8"` (O'Callaghan & Mark 1984): the neighbour of steepest descent (drop over distance; of equal slopes the first from north-east clockwise), as an ESRI code (`D8_CODES`: 1 east, 2 south-east, 4 south ... 128 north-east), 0 where no neighbour is lower. `"dinf"` (Tarboton 1997): the azimuth of steepest descent over the eight triangular facets around the cell, in degrees clockwise from north, -1 where nothing is lower. `"mfd"`: the multiple flow directions of pygeoprocessing (the routing of the InVEST models), each lower neighbour's share in proportion to the drop (over sqrt 2 on diagonals), plateaus draining to their nearest outlets, packed as eight 4-bit weights (`(code >> 4 * k) & 15` for neighbour `k`: 0 east, 1 north-east, 2 north ... 7 south-east), 0 where nothing flows.	`"d8"`
`flats`	`bool`	Drain the flat areas towards lower terrain and away from higher terrain (Barnes, Lehman & Mulla 2014), so that every cell of a filled DEM has a direction; `False` leaves them without (as WhiteboxTools does).	`True`
`extent`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`crs`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`band`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`nodata`		For arrays without georeferencing, or to pick a band and a missing value.	`None`

Returns:

Type	Description
`DataArray`	The codes or azimuths (float64), NaN where the DEM is missing.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> dem = np.array([[3, 3, 3], [3, 2, 3], [3, 1, 3]], float)
>>> ls.hydro.flow_direction(dem).values[1, 1]   # to the south
np.float64(4.0)

flow_accumulation ¶

flow_accumulation(
    dem,
    *,
    method: str = "d8",
    units: str = "cells",
    weights=None,
    exponent: float | None = None,
    flats: bool = True,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

The area that drains through each cell (the cell itself included).

Parameters:

Name	Type	Description	Default
`dem`		The elevations, best with its depressions filled (`fill()`).	required
`method`	`('d8', 'dinf', 'fd8', 'quinn', 'mfd')`	How a cell passes its flow on: all to its steepest neighbour (`"d8"`), split between the two cells of the D-infinity azimuth (`"dinf"`, Tarboton 1997), or among all lower neighbours in proportion to slope^exponent (`"fd8"`, Freeman 1991, exponent 1.1) or to slope times the contour length (0.5 cardinal, 0.354 diagonal; `"quinn"`, Quinn et al. 1991, exponent 1), or pygeoprocessing's multiple flow directions (`"mfd"`: the routing of the InVEST models, see `flow_direction()`; `flats` and `exponent` do not apply). The multiple-flow methods avoid the parallel lines of D8 on hillslopes.	`"d8"`
`units`	`('cells', 'area', 'sca')`	Number of cells, area (cells times the cell area, in map units squared; square metres for longitude / latitude), or specific contributing area (the area per unit contour width: area divided by the mean cell side, as WhiteboxTools does).	`"cells"`
`weights`		A raster of what each cell contributes instead of 1 (runoff, sediment...).	`None`
`exponent`	`float \| None`	The exponent of the multiple-flow methods (default 1.1 for `"fd8"`, 1 for `"quinn"`).	`None`
`flats`	`bool`	Drain flat areas (see `flow_direction()`); without it their cells keep what they receive.	`True`
`extent`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`crs`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`band`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`nodata`		For arrays without georeferencing, or to pick a band and a missing value.	`None`

Returns:

Type	Description
`DataArray`	float64, NaN where the DEM is missing.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> dem = np.array([[4, 4, 4], [3, 3, 3], [2, 1, 2]], float)
>>> ls.hydro.flow_accumulation(dem).values
array([[1., 1., 1.],
       [2., 2., 2.],
       [1., 9., 1.]])

streams ¶

streams(
    accumulation, threshold: float, *, trace: float | None = None, dem=None
) -> DataArray

The stream cells: those whose flow accumulation is at least threshold.

Parameters:

Name	Type	Description	Default
`accumulation`		A flow accumulation (`flow_accumulation()`).	required
`threshold`	`float`	The smallest accumulation of a stream, in its units.	required
`trace`	`float \| None`	With `dem`: pygeoprocessing's stream extraction from the `"mfd"` accumulation, as InVEST's sediment model runs it (`trace=0.7`): the cells of at least `threshold` traced up from the outlets, and tributaries down to `trace` x `threshold` kept where they lead to such a stream (with divergent multiple flow, the cells above the threshold do not always join).	`None`
`dem`		The (filled) DEM whose `"mfd"` directions are traced.	`None`

Returns:

Type	Description
`DataArray`	1 on the streams, 0 elsewhere, NaN where the accumulation is missing.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> acc = ls.hydro.flow_accumulation(np.array([[4, 4, 4], [3, 3, 3], [2, 1, 2]], float))
>>> ls.hydro.streams(acc, 4).values
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 1., 0.]])

distance_to_channel ¶

distance_to_channel(
    dem,
    streams,
    *,
    weights=None,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

The mean distance water travels from each cell to a stream, along pygeoprocessing's multiple flow directions (as InVEST's sediment delivery model measures its downslope component).

Each cell's distance is the mean over the neighbours its flow reaches (weighted by its shares) of their distance plus the step to them; only the flow that reaches a stream counts.

Parameters:

Name	Description	Default
`dem`	The elevations, with depressions filled (`fill()`).	required
`streams`	A raster, 1 on the stream cells.	required
`weights`	The length of each step as a raster (the value of the cell stepped onto), instead of 1 (sides) and sqrt 2 (diagonals) cell lengths.	`None`
`extent`	For arrays without georeferencing, or to pick a band and a missing value.	`None`
`crs`	For arrays without georeferencing, or to pick a band and a missing value.	`None`
`band`	For arrays without georeferencing, or to pick a band and a missing value.	`None`
`nodata`	For arrays without georeferencing, or to pick a band and a missing value.	`None`

Returns:

Type	Description
`DataArray`	In cell lengths (or in the units of the weights); 0 on the streams, NaN where no flow reaches a stream.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> dem = np.add.outer(np.arange(4.0, 0, -1), np.zeros(3))      # falling to the south
>>> rivers = np.zeros((4, 3)); rivers[3] = 1
>>> ls.hydro.distance_to_channel(dem, rivers).values[:, 1].round(3)   # part flows diagonally
array([3.628, 2.433, 1.237, 0.   ])

ls.hydro.D8_CODES maps the D8 codes to their directions: 1 east, 2 south-east, 4 south, 8 south-west, 16 west, 32 north-west, 64 north, 128 north-east (0: no lower neighbour).

Stream links and watersheds¶

The division of a stream network into links (between junctions and outlets), the area that drains to each seed and outlets moved onto the streams: the delineation steps of the watershed model (see ls.hydro.swat.watershed).

stream_links ¶

stream_links(
    dem,
    streams,
    *,
    outlets=None,
    snap: float | None = None,
    flats: bool = True,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> tuple[DataArray, DataFrame]

The links of a stream network: the stretches between sources, junctions and outlets.

A link starts at a source (a stream cell nothing on the streams drains into), just below a junction (two or more streams meeting) or just below an outlet, and ends where the next one starts or where the stream leaves the map. Links are numbered 1, 2, ... in the raster order of their first cells.

Parameters:

Name	Type	Description	Default
`dem`		The elevations, best with depressions filled (`fill()`): the D8 directions (`flow_direction()`) are taken from it.	required
`streams`		1 on the stream cells (`streams()`), on the DEM's grid.	required
`outlets`		Points where a link must end (gauges, dams, the watershed outlet): (x, y) pairs, an (n, 2) array or a GeoDataFrame; each snaps to the nearest stream cell (`snap_outlets()`, within `snap`).	`None`
`snap`	`float \| None`	The largest move of an outlet, in map units (default: 300 m).	`None`
`flats`	`bool`	Drain flat areas (as `flow_direction()`).	`True`
`extent`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`crs`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`band`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`nodata`		For arrays without georeferencing, or to pick a band and a missing value.	`None`

Returns:

Name	Type	Description
`links`	`DataArray`	The link of each stream cell (1..n), 0 elsewhere on the DEM, NaN where it is missing.
`table`	`DataFrame`	By link (index `link`): `downstream` (the link its last cell drains into, 0 for none), `outlet` (whether it ends at one of `outlets`), the first (`start_row`, `start_col`) and last cell (`end_row`, `end_col`), `cells`, `length` (along the D8 path, from the first cell's centre to the centre of the cell below the last one; a last cell that drains nowhere counts one cell width; metres for projected and unprojected grids), the elevations `z_start` and `z_end` (of the first cell and of the cell below the last one, or of the last one when nothing is below), and `order` (Strahler).

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> dem = np.array([[5, 4, 5], [5, 3, 5], [5, 2, 5], [5, 1, 5]], float)
>>> st = ls.hydro.streams(ls.hydro.flow_accumulation(dem), 2)
>>> links, table = ls.hydro.stream_links(dem, st)
>>> links.values
array([[0., 0., 0.],
       [0., 1., 0.],
       [0., 1., 0.],
       [0., 1., 0.]])

watersheds ¶

watersheds(
    dem,
    seeds,
    *,
    flats: bool = True,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

The area draining to each seed: every cell labelled with the first seed below it.

Parameters:

Name	Type	Description	Default
`dem`		The elevations, best with depressions filled (`fill()`).	required
`seeds`		Positive integer labels on the cells water drains to (stream links from `stream_links()`, outlet cells, lakes), 0 or NaN elsewhere; seed cells keep their own label.	required
`flats`	`bool`	As `stream_links()`.	`True`
`extent`	`bool`	As `stream_links()`.	`True`
`crs`	`bool`	As `stream_links()`.	`True`
`band`	`bool`	As `stream_links()`.	`True`
`nodata`	`bool`	As `stream_links()`.	`True`

Returns:

Type	Description
`DataArray`	The label of the seed each cell drains to, 0 where its path reaches none (it leaves the map or ends in a pit first), NaN where the DEM is missing.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> dem = np.array([[3, 2, 3], [3, 1, 3]], float)
>>> ls.hydro.watersheds(dem, np.array([[0, 0, 0], [0, 7, 0]])).values
array([[7., 7., 7.],
       [7., 7., 7.]])

snap_outlets ¶

snap_outlets(
    points, streams, *, distance: float | None = None, accumulation=None
) -> ndarray

The stream cell each outlet point stands for.

Each point moves to the stream cell nearest to it (cell centres, in map units), within distance; of equally near cells, the one with the largest accumulation (then the first in raster order).

Parameters:

Name	Type	Description	Default
`points`		(x, y) pairs, an (n, 2) array or a GeoDataFrame of points (map coordinates of the stream raster).	required
`streams`		1 on the streams (`streams()`).	required
`distance`	`float \| None`	The largest move, in map units (default: 300 m, QSWAT+'s snap threshold, or 3 cells for unprojected grids).	`None`
`accumulation`		A flow accumulation to break ties (optional).	`None`

Returns:

Type	Description
`ndarray`	(n, 2) row and column of the cell of each point.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> dem = np.array([[5, 4, 5], [5, 3, 5], [5, 2, 5], [5, 1, 5]], float)
>>> st = ls.hydro.streams(ls.hydro.flow_accumulation(dem), 2)

a point half a cell west of the stream, in the raster's x, y (cell centres 0.5, 1.5...):

>>> ls.hydro.snap_outlets([(0.9, 1.5)], st, distance=1.0)
array([[2, 1]])

Height above the nearest drainage¶

hand ¶

hand(
    dem,
    streams=None,
    *,
    threshold: float | None = None,
    flats: bool = True,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Height Above the Nearest Drainage (HAND; Rennó et al. 2008, Nobre et al. 2011).

The height of each cell above the stream cell where its water reaches the drainage, following the D8 directions downhill: 0 on the streams, small on floodplains and valley floors, large on hilltops. A terrain descriptor of soil water and flood susceptibility developed at INPE.

Parameters:

Name	Type	Description	Default
`dem`		The elevations, with depressions filled (`fill()`), so that every cell reaches a stream.	required
`streams`		The drainage: a raster, non-zero on the stream cells (`streams()`).	`None`
`threshold`	`float \| None`	Or the drainage as the cells of D8 accumulation (in cells) of at least this.	`None`
`flats`	`bool`	Drain flat areas (see `flow_direction()`).	`True`
`extent`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`crs`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`band`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`nodata`		For arrays without georeferencing, or to pick a band and a missing value.	`None`

Returns:

Type	Description
`DataArray`	The height in the DEM's units, NaN where the DEM is missing or the cell reaches no stream.

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> dem = np.array([[9, 9, 9], [6, 5, 6], [4, 3, 4]], float)
>>> ls.hydro.hand(dem, threshold=3).values
array([[6., 6., 6.],
       [3., 2., 3.],
       [1., 0., 1.]])

Topographic indices¶

twi ¶

twi(
    dem,
    *,
    method: str = "fd8",
    min_slope: float = 0.0,
    sca=None,
    slope=None,
    flats: bool = True,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Topographic wetness index: ln(a / tan(slope)) (Beven & Kirkby 1979).

High where much area drains through a gentle slope (valley floors, hollows): a proxy of soil moisture and saturation used in species distribution and land change models.

Parameters:

Name	Type	Description	Default
`dem`		The elevations, with depressions filled (`fill()`).	required
`method`	`str`	The flow accumulation of the specific contributing area `a` (see `flow_accumulation()`): `"fd8"` (default), `"quinn"`, `"dinf"`, `"d8"`.	`'fd8'`
`min_slope`	`float`	A floor for tan(slope), so flat cells get a value (e.g. 0.001); with 0 they are NaN, as in WhiteboxTools.	`0.0`
`sca`		A specific contributing area raster to use instead.	`None`
`slope`		A slope raster in degrees to use instead of the gradient of the DEM.	`None`
`flats`	`bool`	Drain flat areas (see `flow_direction()`).	`True`
`extent`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`crs`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`band`		For arrays without georeferencing, or to pick a band and a missing value.	`None`
`nodata`		For arrays without georeferencing, or to pick a band and a missing value.	`None`

Returns:

Type	Description
`DataArray`

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> dem = np.add.outer(np.arange(5.0, 0, -1), np.zeros(4))   # a plane falling to the south
>>> w = ls.hydro.twi(dem)
>>> bool(w.values[4, 1] > w.values[0, 1])                     # wetter downslope
True

spi ¶

spi(
    dem,
    *,
    method: str = "fd8",
    sca=None,
    slope=None,
    flats: bool = True,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Stream power index: a · tan(slope) (Moore, Grayson & Ladson 1991), the erosive power of overland flow.

Parameters:

Name	Description	Default
`dem`	As in `twi()`.	required
`method`	As in `twi()`.	required
`sca`	As in `twi()`.	required
`slope`	As in `twi()`.	required
`flats`	As in `twi()`.	required
`extent`	As in `twi()`.	required
`crs`	As in `twi()`.	required
`band`	As in `twi()`.	required
`nodata`	As in `twi()`.	required

Returns:

Type	Description
`DataArray`

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> dem = np.add.outer(np.arange(5.0, 0, -1), np.zeros(4))
>>> bool(ls.hydro.spi(dem).values[3, 1] > ls.hydro.spi(dem).values[1, 1])
True

sti ¶

sti(
    dem,
    *,
    method: str = "fd8",
    sca_exponent: float = 0.4,
    slope_exponent: float = 1.3,
    sca=None,
    slope=None,
    flats: bool = True,
    extent=None,
    crs=None,
    band=None,
    nodata=None,
    n_threads: int | None = None,
) -> DataArray

Sediment transport index, the LS factor of the USLE from the contributing area (Moore & Burch 1986): (m + 1) (a / 22.13)^m (sin(slope) / 0.0896)^n.

Parameters:

Name	Type	Description	Default
`dem`		As in `twi()`.	required
`method`		As in `twi()`.	required
`sca`		As in `twi()`.	required
`slope`		As in `twi()`.	required
`flats`		As in `twi()`.	required
`extent`		As in `twi()`.	required
`crs`		As in `twi()`.	required
`band`		As in `twi()`.	required
`nodata`		As in `twi()`.	required
`sca_exponent`	`float`	m (0.4) and n (1.3).	`0.4`
`slope_exponent`	`float`	m (0.4) and n (1.3).	`0.4`

Returns:

Type	Description
`DataArray`

Examples:

>>> import numpy as np
>>> import landschaft as ls
>>> dem = np.add.outer(np.arange(5.0, 0, -1), np.zeros(4))
>>> float(ls.hydro.sti(dem).values[0, 0]) > 0
True

Sediment connectivity¶

The index of connectivity of Borselli et al. (2008) and Cavalli et al. (2013), as InVEST's sediment delivery ratio model computes it (ls.services.sdr).

connectivity_index ¶

connectivity_index(
    dem,
    weight,
    *,
    streams=None,
    threshold: float | None = None,
    method: str = "mfd",
    n_threads: int | None = None,
) -> DataArray

The index of connectivity of Borselli et al. (2008), as InVEST's SDR model computes it.

IC = log10(D_up / D_dn): D_up = W S sqrt(A), the weight W and the slope S (m/m) averaged over the contributing area A of each cell, and D_dn = sum d_i / (W_i S_i) along the flow down to a stream. Higher values: cells better connected to the streams. The DEM is filled and routed by pygeoprocessing's algorithms (multiple flow directions or D8), the slope is pygeoprocessing's percent slope; slopes are kept between 0.005 and 1 and the weights from 0.001 (Cavalli et al. 2013), as in InVEST.

Parameters:

Name	Type	Description	Default
`dem`		Elevations (projected, in metres).	required
`weight`		The impedance W, usually the USLE C factor of the land cover (0 to 1), on the DEM's grid.	required
`streams`		The stream cells (1), on the DEM's grid; or `threshold`.	`None`
`threshold`	`float \| None`	Flow accumulation (cells) from which a cell is a stream (traced as InVEST).	`None`
`method`	`str`	`"mfd"` (default) or `"d8"`.	`'mfd'`

Returns:

Type	Description
`DataArray`	The index (float32), NaN on the streams and where the flow reaches none.

Examples:

>>> import numpy as np
>>> from affine import Affine
>>> import landschaft as ls
>>> y, x = np.mgrid[:30, :40]
>>> z = (50 + 0.8 * x + 3 * np.cos(y / 4.0)).astype("float32")
>>> dem = ls.load_raster(ls.write_raster(z, "ic_dem.tif", crs="EPSG:31983",
...                                      transform=Affine(30, 0, 500000, 0, -30, 8000000)))
>>> c = dem.copy(data=np.where(x > 20, 0.2, 0.01))                   # crops uphill of forest
>>> ic = ls.hydro.connectivity_index(dem, c, threshold=30)
>>> bool(np.nanmean(ic.values[:, 25:]) > np.nanmean(ic.values[:, 5:15]))
True