Ecosystem services¶

Models of InVEST (the Natural Capital Project's natcap.invest, Apache-2.0), reimplemented in C++ from land use maps and parameter tables: their maps, bit for bit where the arithmetic permits. A simulated land use map (ls.simulate) or a series of scenarios goes straight in.

Carbon storage and sequestration¶

carbon ¶

carbon(
    lulc,
    pools,
    *,
    future=None,
    price: float | None = None,
    discount_rate: float | None = None,
    rate_change: float = 0.0,
    year: int | None = None,
    future_year: int | None = None,
    nodata=None,
    n_threads: int | None = None,
) -> Carbon

Carbon storage of a landscape and its change between two maps (InVEST's carbon model).

Each land use class holds a density of carbon (tons per hectare) in four pools: aboveground and belowground biomass, soil and dead organic matter. The storage of a map is the sum of its pools cell by cell; with a future map, the change is the sequestration (future minus baseline), and with a price, its net present value.

The maps are those of InVEST's carbon model (natcap.invest 3.20), bit for bit: float32 densities, the pools summed in their order, the change and its value in float32; the landscape totals agree to the last digits.

Parameters:

Name	Type	Description	Default
`lulc`		The baseline land use map (integer codes; a projected raster in metres).	required
`pools`		The densities by class: a DataFrame with `lucode` (or indexed by the codes) and `c_above`, `c_below`, `c_soil`, `c_dead` (Mg C / ha), or a `dict` `{code: {"c_above": ..., ...}}`. Every class of the maps must be in it.	required
`future`		The alternate (future) land use map, on the same grid.	`None`
`price`	`float \| None`	Valuation of the change (with `future`): the price of a ton of carbon and the annual discount rate and change of the price, in percent.	`None`
`discount_rate`	`float \| None`	Valuation of the change (with `future`): the price of a ton of carbon and the annual discount rate and change of the price, in percent.	`None`
`rate_change`	`float \| None`	Valuation of the change (with `future`): the price of a ton of carbon and the annual discount rate and change of the price, in percent.	`None`
`year`	`int \| None`	The years of the two maps (needed with `price`).	`None`
`future_year`	`int \| None`	The years of the two maps (needed with `price`).	`None`
`nodata`		Value marking missing cells in the maps (besides their own nodata).	`None`

Returns:

Type	Description
`Carbon`

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> from affine import Affine
>>> import landschaft as ls
>>> lulc = ls.write_raster(np.array([[1, 1, 2], [2, 3, 3]]), "lulc.tif", crs="EPSG:31983",
...                        transform=Affine(30, 0, 500000, 0, -30, 8000000))   # 30 m cells
>>> pools = pd.DataFrame({"lucode": [1, 2, 3], "c_above": [140, 35, 5], "c_below": [40, 10, 2],
...                       "c_soil": [90, 80, 60], "c_dead": [12, 3, 0]})
>>> later = np.array([[1, 2, 2], [2, 3, 3]])
>>> c = ls.services.carbon(ls.load_raster(lulc), pools, future=later)
>>> c.storage.values[0, 0], c.change.values[0, 1]
(np.float32(282.0), np.float32(-154.0))
>>> c.totals
                      label  total units
0   Baseline Carbon Storage  85.86     t
1  Alternate Carbon Storage  72.00     t
2  Change in Carbon Storage -13.86     t

Carbon `dataclass` ¶

Carbon(
    storage: DataArray,
    pools: dict,
    totals: DataFrame,
    storage_future: DataArray | None = None,
    pools_future: dict = dict(),
    change: DataArray | None = None,
    npv: DataArray | None = None,
    valuation: float | None = None,
)

Bases: Saveable

Carbon stored in a landscape, and its change (carbon()).

Maps are densities in metric tons of carbon per hectare (float32, as InVEST writes them), NaN where the land use map has no value.

Attributes:

Name	Type	Description
`storage`	`DataArray`	Total carbon of the baseline map: the sum of the four pools.
`pools`	`dict`	Each pool's map (`c_above`, `c_below`, `c_soil`, `c_dead`).
`storage_future, pools_future`		The same for the alternate (future) map, if given.
`change`	`DataArray or None`	Sequestration: future minus baseline storage.
`npv`	`DataArray or None`	Net present value of the change (currency per hectare), with valuation.
`totals`	`DataFrame`	The landscape totals (`label`, `total`, `units`): tons of carbon (densities times the cell areas) and currency.
`valuation`	`float or None`	The factor turning a change in carbon into its net present value.

valuation_constant ¶

valuation_constant(
    price: float, discount_rate: float, rate_change: float, year: int, future_year: int
) -> float

The factor turning a change in carbon into its net present value (InVEST's).

The change of storage between the two years, spread evenly over them, valued at price per ton rising by rate_change percent a year and discounted at discount_rate percent a year: price / n * sum_{t=0}^{n-1} 1 / ((1 + d/100) (1 + r/100))^t, n the years.

Examples:

>>> import landschaft as ls
>>> round(ls.services.valuation_constant(43.0, 7.0, 0.0, 2000, 2030), 6)
19.031333

Annual water yield and hydropower¶

water_yield ¶

water_yield(
    lulc,
    biophysical,
    precipitation,
    eto,
    depth,
    pawc,
    *,
    seasonality: float,
    watersheds=None,
    subwatersheds=None,
    demand=None,
    valuation=None,
    nodata=None,
    n_threads: int | None = None,
) -> WaterYield

Annual water yield, consumption and hydropower value (InVEST's annual water yield model).

The water balance of each cell on the Budyko curve in Fu and Zhang's form: actual evapotranspiration over precipitation AET/P = 1 + PET/P - (1 + (PET/P)^w)^(1/w), at most PET/P, with PET = Kc ET0 and w = Z AWC / P + 1.25 (at most 5), AWC = min(root depth, soil depth) x PAWC; classes that are not vegetation (water, wetlands, towns) evaporate min(Kc ET0, P). The yield is P - AET. With watersheds, the means over each polygon and the yield volume; with a demand table, the consumption and the realized supply; with a valuation table, the hydropower energy and its net present value.

The maps are InVEST's (natcap.invest 3.20) bit for bit, in the same float32 (or float64) arithmetic; the watershed means agree to float32 precision (InVEST adds the cells in float32).

Parameters:

Name	Type	Description	Default
`lulc`		Land use map (integer codes; projected, in metres).	required
`biophysical`		A DataFrame by class (`lucode`): `lulc_veg` (1 vegetation, 0 not), `root_depth` (mm) and `kc` (crop coefficient).	required
`precipitation`		Annual precipitation and reference evapotranspiration (mm), on the grid of `lulc`. Cells with no rain are left out, as InVEST does.	required
`eto`		Annual precipitation and reference evapotranspiration (mm), on the grid of `lulc`. Cells with no rain are left out, as InVEST does.	required
`depth`		Depth to the root restricting layer (mm).	required
`pawc`		Plant available water content (fraction, 0 to 1).	required
`seasonality`	`float`	The Z parameter (about 1 to 30): the seasonal distribution of rain.	required
`watersheds`		Polygons (GeoDataFrame) with a `ws_id` (`subws_id`) column, in the maps' CRS.	`None`
`subwatersheds`		Polygons (GeoDataFrame) with a `ws_id` (`subws_id`) column, in the maps' CRS.	`None`
`demand`		A DataFrame by class (`lucode`, `demand`): consumptive use, m3 per year and cell.	`None`
`valuation`		A DataFrame by watershed (`ws_id`): `efficiency`, `fraction`, `height` (m), `kw_price`, `cost` (per year), `time_span` (years) and `discount` (percent per year). Needs the demand table.	`None`
`nodata`		Value marking missing cells in the land use map.	`None`

Returns:

Type	Description
`WaterYield`

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> from affine import Affine
>>> import landschaft as ls
>>> codes = np.array([[1, 1, 2], [1, 2, 3]])                 # forest, pasture, water
>>> lulc = ls.load_raster(ls.write_raster(codes, "wy_lulc.tif", crs="EPSG:31983",
...                                       transform=Affine(30, 0, 500000, 0, -30, 8000000)))
>>> bio = pd.DataFrame({"lucode": [1, 2, 3], "lulc_veg": [1, 1, 0], "root_depth": [5000, 1500, 1],
...                     "kc": [1.0, 0.65, 1.05]})
>>> rain, et0 = lulc.copy(data=np.full((2, 3), 1500.0)), lulc.copy(data=np.full((2, 3), 1200.0))
>>> soil, pawc = lulc.copy(data=np.full((2, 3), 2000.0)), lulc.copy(data=np.full((2, 3), 0.2))
>>> wy = ls.services.water_yield(lulc, bio, rain, et0, soil, pawc, seasonality=5.0)
>>> wy.wyield.values.round(1)                                # mm per year
array([[582.6, 582.6, 864.3],
       [582.6, 864.3, 240. ]])

WaterYield `dataclass` ¶

WaterYield(
    wyield: DataArray,
    aet: DataArray,
    pet: DataArray,
    fractp: DataArray,
    demand: DataArray | None = None,
    watersheds: object = None,
    subwatersheds: object = None,
)

Bases: Saveable

Annual water yield of a landscape (water_yield()).

Maps are in millimetres per year, NaN where an input has no value (or there is no rain, as InVEST leaves those cells out of the balance); float32, or float64 where InVEST computes in float64 (float64 inputs).

Attributes:

Name	Type	Description
`wyield`	`DataArray`	Water yield: precipitation minus actual evapotranspiration.
`aet`	`DataArray`	Actual evapotranspiration.
`pet`	`DataArray`	Potential evapotranspiration: reference evapotranspiration times Kc.
`fractp`	`DataArray`	Actual evapotranspiration as a fraction of precipitation (float32).
`demand`	`DataArray or None`	Consumptive water use per cell (m3 per year), with a demand table.
`watersheds, subwatersheds`	`GeoDataFrame or None`	The polygons with `precip_mn`, `PET_mn`, `AET_mn`, `wyield_mn` (means over the cells with values, mm), `wyield_vol` (the mean yield over the polygon's area, m3), with demand `consum_vol` (m3), `consum_mn` (m3 per cell), `rsupply_vl` (realized supply, yield minus consumption, m3) and `rsupply_mn` (m3 per cell), and with valuation (watersheds) `hp_energy` (kWh per year) and `hp_val` (net present value).

Sediment delivery ratio¶

sdr ¶

sdr(
    dem,
    erosivity,
    erodibility,
    lulc,
    biophysical,
    watersheds,
    *,
    threshold: float,
    k: float = 2.0,
    ic_0: float = 0.5,
    sdr_max: float = 0.8,
    l_max: float = 122.0,
    method: str = "mfd",
    drainage=None,
    nodata=None,
    n_threads: int | None = None,
) -> SDR

Sediment export, deposition and retention (InVEST's sediment delivery ratio model).

The soil loss of each cell by the revised universal soil loss equation, USLE = R K LS C P (the LS factor of Desmet & Govers 1996 from the flow accumulation), and the share of it that reaches a stream, the sediment delivery ratio SDR = SDR_max / (1 + exp((IC_0 - IC) / k)) of the index of connectivity IC = log10(D_up / D_dn) (Borselli et al. 2008): D_up = W S sqrt(A) upslope (W the C factor, S the slope, averaged over the contributing area A) and D_dn = sum d / (W S) down to the stream. The sediment that does not reach a stream is deposited downslope; the retention is what the vegetation and practices avoid, plus that deposition.

The rasters are InVEST's (natcap.invest 3.20) bit for bit, in its float32 and float64 types: the DEM filled, routed (pygeoprocessing's multiple flow directions or D8) and traced into streams as pygeoprocessing does, its percent slope, and the deposition; the watershed totals agree to float32 precision.

Parameters:

Name	Type	Description	Default
`dem`		Elevations (projected, in metres).	required
`erosivity`		Rainfall erosivity R (MJ mm / (ha h yr)) and soil erodibility K (t ha h / (ha MJ mm)), on the DEM's grid.	required
`erodibility`		Rainfall erosivity R (MJ mm / (ha h yr)) and soil erodibility K (t ha h / (ha MJ mm)), on the DEM's grid.	required
`lulc`		Land use map (integer codes).	required
`biophysical`		A DataFrame by class (`lucode`): `usle_c` and `usle_p`.	required
`watersheds`		Polygons (GeoDataFrame): cells outside them are left out (as InVEST masks the inputs to them), and the totals are given per polygon.	required
`threshold`	`float`	Flow accumulation (cells) from which a cell is a stream.	required
`k`	`float`	Borselli's calibration parameters of the SDR.	`2.0`
`ic_0`	`float`	Borselli's calibration parameters of the SDR.	`2.0`
`sdr_max`	`float`	Borselli's calibration parameters of the SDR.	`2.0`
`l_max`	`float`	The largest slope length factor L.	`122.0`
`method`	`str`	`"mfd"` (multiple flow directions, the default) or `"d8"`.	`'mfd'`
`drainage`		Cells that drain to a stream otherwise (1), such as roads or ditches.	`None`
`nodata`		Value marking missing cells in the land use map.	`None`

Returns:

Type	Description
`SDR`

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> import geopandas as gpd
>>> from affine import Affine
>>> from shapely.geometry import box
>>> import landschaft as ls
>>> y, x = np.mgrid[:40, :50]
>>> z = (100 + 0.5 * x + 8 * np.sin(y / 5.0)).astype("float32")      # a slope with ridges
>>> t = Affine(30, 0, 500000, 0, -30, 8000000)
>>> dem = ls.load_raster(ls.write_raster(z, "sdr_dem.tif", crs="EPSG:31983", transform=t))
>>> lulc = dem.copy(data=np.where(x < 25, 1, 2).astype("int16"))    # forest, crops
>>> r, k = dem.copy(data=np.full(z.shape, 4000, "float32")), dem.copy(data=np.full(z.shape, 0.03, "float32"))
>>> bio = pd.DataFrame({"lucode": [1, 2], "usle_c": [0.003, 0.25], "usle_p": [1.0, 0.8]})
>>> ws = gpd.GeoDataFrame({"ws_id": [1]}, geometry=[box(*dem.rio.bounds())], crs=dem.rio.crs)
>>> s = ls.services.sdr(dem, r, k, lulc, bio, ws, threshold=50)
>>> crops, forest = np.nanmean(s.usle.values[:, 30:]), np.nanmean(s.usle.values[:, :20])
>>> bool(crops > 10 * forest)                     # crops lose far more soil (t / ha / yr)
True
>>> list(s.watersheds.columns)
['ws_id', 'geometry', 'usle_tot', 'sed_export', 'sed_dep', 'avoid_exp', 'avoid_eros']

SDR `dataclass` ¶

SDR(
    usle: DataArray,
    sed_export: DataArray,
    sed_deposition: DataArray,
    avoided_export: DataArray,
    avoided_erosion: DataArray,
    rkls: DataArray,
    sdr: DataArray,
    ic: DataArray,
    stream: DataArray,
    drains_to_stream: DataArray,
    watersheds: object,
    intermediate: dict = dict(),
)

Bases: Saveable

Sediment export and retention of a landscape (sdr()).

Maps in tons per hectare and year (float32), NaN where a cell has no value; usle, rkls and avoided_erosion have none on the streams, as in InVEST.

Attributes:

Name	Type	Description
`usle`	`DataArray`	Soil loss, RKLS x C x P.
`sed_export`	`DataArray`	The soil loss that reaches the streams: USLE x SDR.
`sed_deposition`	`DataArray`	Sediment from upslope that the cell traps.
`avoided_export`	`DataArray`	Total retention: avoided erosion x SDR plus the deposition.
`avoided_erosion`	`DataArray`	Erosion the vegetation and practices avoid: RKLS - USLE.
`rkls`	`DataArray`	Soil loss of bare soil, R x K x LS.
`sdr`	`DataArray`	Sediment delivery ratio (1 on the streams).
`ic`	`DataArray`	Index of connectivity, log10(D_up / D_dn).
`stream`	`DataArray`	The streams (1), with the drainage given (0 elsewhere).
`drains_to_stream`	`DataArray`	1 where the flow reaches a stream, 0 where it does not (no D_dn).
`watersheds`	`GeoDataFrame`	The watersheds with `usle_tot`, `sed_export`, `sed_dep`, `avoid_exp` and `avoid_eros` (tons per year).
`intermediate`	`dict`	The other rasters of the chain (`slope`, `flow_accumulation`, `ls`, `w`, `cp`, `w_bar`, `s_bar`, `d_up`, `d_dn`, `e_prime`, `flux`...), in InVEST's types.

Habitat quality and rarity¶

habitat_quality ¶

habitat_quality(
    lulc,
    threats,
    sensitivity,
    threat_maps,
    *,
    half_saturation: float,
    future=None,
    future_threat_maps=None,
    base=None,
    base_threat_maps=None,
    access=None,
    nodata=None,
    n_threads: int | None = None,
) -> HabitatQuality

Habitat degradation, quality and rarity (InVEST's habitat quality model).

Each threat (crops, roads, towns...) degrades the habitat around it: its impact decays linearly or exponentially with the distance to the nearest cell of the threat, up to its maximum distance (exponentially: to 5 % there), and weighs by its relative weight, by the sensitivity of the cell's land use class to it and by the cell's access (legal protection). The degradation D of a cell is the sum over the threats; its quality is h (1 - D^z / (D^z + k^z)), h the habitat suitability of its class, k the half-saturation constant and z = 2.5. With a base map, the rarity of each class of the current (and future) map: 1 - n_x / (n_base + n_x), the cell counts in the two maps.

The maps are InVEST's (natcap.invest 3.20, which takes the distance to the nearest threat cell since 3.12): the same distance transform in float32 and the same float32 expressions, bit for bit (the exponential decay through numpy's float32 exp, as InVEST computes it).

Parameters:

Name	Type	Description	Default
`lulc`		The current land use map (integer codes; projected, in metres).	required
`threats`		A DataFrame with a row per threat: `threat` (its name), `max_dist` (metres), `weight` (relative) and `decay` (`"linear"` or `"exponential"`).	required
`sensitivity`		A DataFrame by land use class (`lucode`): `habitat` (0 to 1) and a column per threat with the sensitivity of the class to it (0 to 1).	required
`threat_maps`		The current map of each threat, `{threat: raster}` (or a 3-D DataArray with a band per threat), on the grid of `lulc`: cells other than 0 and nodata are the threat (values from 0 to 1).	required
`half_saturation`	`float`	The half-saturation constant k (often half the highest degradation).	required
`future`		A future land use map and its threat maps.	`None`
`future_threat_maps`		A future land use map and its threat maps.	`None`
`base`		A base (reference) land use map, for the rarity, and its threat maps (optional: without them the base map has no degradation nor quality).	`None`
`base_threat_maps`		A base (reference) land use map, for the rarity, and its threat maps (optional: without them the base map has no degradation nor quality).	`None`
`access`		Accessibility to threats, from 0 (inaccessible) to 1: a raster, or polygons (GeoDataFrame) with an `access` column (1 outside them).	`None`
`nodata`		Value marking missing cells in the land use maps.	`None`

Returns:

Type	Description
`HabitatQuality`

Examples:

>>> import numpy as np
>>> import pandas as pd
>>> from affine import Affine
>>> import landschaft as ls
>>> lulc = np.ones((40, 40), int)
>>> lulc[:, 20:] = 2                                       # forest (1), pasture (2)
>>> road = np.zeros((40, 40))
>>> road[:, 30] = 1
>>> lulc = ls.load_raster(ls.write_raster(lulc, "hq_lulc.tif", crs="EPSG:31983",
...                                       transform=Affine(30, 0, 500000, 0, -30, 8000000)))
>>> threats = pd.DataFrame({"threat": ["road"], "max_dist": [600.0], "weight": [1.0],
...                         "decay": ["linear"]})
>>> sens = pd.DataFrame({"lucode": [1, 2], "habitat": [1.0, 0.0], "road": [0.8, 0.0]})
>>> hq = ls.services.habitat_quality(lulc, threats, sens, {"road": lulc.copy(data=road)},
...                                  half_saturation=0.1)
>>> q = hq.quality["current"].values
>>> q[0, 0] > q[0, 19], q[0, 25]                           # less quality near the road
(np.True_, np.float32(0.0))

HabitatQuality `dataclass` ¶

HabitatQuality(
    degradation: dict,
    quality: dict,
    threats: dict,
    rarity: dict = dict(),
    rarity_table: DataFrame | None = None,
)

Bases: Saveable

Habitat degradation, quality and rarity (habitat_quality()).

Maps are float32 (as InVEST writes them), NaN where the land use map has no value. degradation, quality, threats and rarity are dicts by scenario: "current" and, when given, "future" and "base".

Attributes:

Name	Type	Description
`degradation`	`dict`	Total degradation of each scenario: the sum over the threats.
`quality`	`dict`	Habitat quality: `h (1 - D^z / (D^z + k^z))`, z = 2.5.
`threats`	`dict`	Each scenario's degradation by threat (a dict by threat name).
`rarity`	`dict`	Rarity of the current and future map against the base map: by cell (`rarity[scenario]`), and by class in `rarity_table`.
`rarity_table`	`DataFrame`	`scenario`, `lucode`, `rarity`: 1 - n_x / (n_base + n_x), the cell counts of the class in the scenario and the base map (0 for a class the base map lacks).

Ecosystem services¶

Carbon storage and sequestration¶

carbon ¶

Carbon dataclass ¶

valuation_constant ¶

Annual water yield and hydropower¶

water_yield ¶

WaterYield dataclass ¶

Sediment delivery ratio¶

sdr ¶

SDR dataclass ¶

Habitat quality and rarity¶

habitat_quality ¶

HabitatQuality dataclass ¶

Carbon `dataclass` ¶

WaterYield `dataclass` ¶

SDR `dataclass` ¶

HabitatQuality `dataclass` ¶