API Reference

Spectral Indices

indices.py

Spectral index computation from multispectral xarray DataArrays.

Supported indices

MNDWI Modified Normalised Difference Water Index (Xu 2006)
= (Green - SWIR) / (Green + SWIR)
NDWI Normalised Difference Water Index (McFeeters 1996)
= (Green - NIR) / (Green + NIR)
NDVI Normalised Difference Vegetation Index
= (NIR - Red) / (NIR + Red)
NDTI Normalised Difference Turbidity Index
= (Red - Green) / (Red + Green)
AWEIsh Automated Water Extraction Index Shadow (Feyisa et al. 2014)
= Blue + 2.5*Green - 1.5*(NIR + SWIR) - 0.25*SWIR2
AWEInsh Automated Water Extraction Index No Shadow (Feyisa et al. 2014)
= Blue + 2.5*Green - 1.5*(NIR + SWIR)

Functions

compute_mndwi() : Single MNDWI computation
compute_ndwi() : Single NDWI computation
compute_ndvi() : Single NDVI computation
compute_ndti() : Single NDTI computation
compute_aweish() : Single AWEIsh computation
compute_aweinsh() : Single AWEInsh computation
compute_indices() : MNDWI, NDVI, NDTI (+ optionally AWEIsh, AWEInsh for WCT)
compute_water_indices() : All water indices (MNDWI, NDWI, AWEIsh, AWEInsh)

MNDWI is recommended for most wetland applications. AWEIsh and AWEInsh may perform better in areas with strong topographic shadow (mountains, deep valleys) or confusion with built-up surfaces.

wetlandmapper.indices.compute_aweinsh(ds: Dataset | DataArray, green_band: str = 'green', nir_band: str = 'nir', swir_band: str = 'swir') → DataArray[source]

Compute the Automated Water Extraction Index without shadow suppression.

AWEInsh = 4*(Green - SWIR1) - (0.25*NIR + 2.75*SWIR1)

A simplified water index that omits SWIR2 and is computationally lighter than AWEIsh. Less effective at suppressing shadow, but robust in areas affected by atmospheric haze or aerosol loading.

Parameters:

ds (xr.Dataset or xr.DataArray) – Multi-spectral dataset. Input bands should be surface reflectance in the range [0, 1].
green_band (str) – Default "green".
nir_band (str) – Default "nir".
swir_band (str) – SWIR1. Default "swir".

Returns:

AWEInsh values. A threshold of 0.0 separates water from non-water.

Return type:

xr.DataArray

Notes

The 0.0001 DN scaling factor from the original paper is omitted here because this function expects surface reflectance in [0, 1].

References

Feyisa, G. L., Meilby, H., Fensholt, R., & Proud, S. R. (2014). Remote Sensing of Environment, 140, 23-35. https://doi.org/10.1016/j.rse.2013.08.029

wetlandmapper.indices.compute_aweish(ds: Dataset | DataArray, blue_band: str = 'blue', green_band: str = 'green', nir_band: str = 'nir', swir_band: str = 'swir', swir2_band: str = 'swir2') → DataArray[source]

Compute the Automated Water Extraction Index with shadow suppression.

AWEIsh = Blue + 2.5*Green - 1.5*(NIR + SWIR1) - 0.25*SWIR2

Designed to suppress topographic shadow and built-up surface confusion that affects simpler water indices such as MNDWI. Positive values indicate surface water. Preferred over MNDWI in mountainous terrain or urban areas.

Parameters:

ds (xr.Dataset or xr.DataArray) – Multi-spectral dataset. Input bands should be surface reflectance in the range [0, 1] (not raw DN values).
blue_band (str) – Default "blue". Landsat 8/9: "SR_B2"; Sentinel-2: "B2".
green_band (str) – Default "green".
nir_band (str) – Default "nir".
swir_band (str) – SWIR1. Default "swir".
swir2_band (str) – SWIR2. Default "swir2". Landsat 8/9: "SR_B7"; Sentinel-2: "B12".

Returns:

AWEIsh values. A threshold of 0.0 separates water (positive) from non-water (negative).

Return type:

xr.DataArray

Notes

The original Feyisa et al. (2014) formula was defined for raw Landsat DN values and includes a 0.0001 scaling factor. This implementation operates on surface reflectance already scaled to [0, 1], so the constant is omitted.

References

Feyisa, G. L., Meilby, H., Fensholt, R., & Proud, S. R. (2014). Automated Water Extraction Index: A new technique for surface water mapping using Landsat imagery. Remote Sensing of Environment, 140, 23-35. https://doi.org/10.1016/j.rse.2013.08.029

wetlandmapper.indices.compute_indices(ds: Dataset | DataArray, green_band: str = 'green', red_band: str = 'red', nir_band: str = 'nir', swir_band: str = 'swir', swir2_band: str = 'swir2', blue_band: str = 'blue', include_awei: bool = False) → Dataset[source]

Compute spectral indices and return as a Dataset.

Always computes MNDWI, NDVI, and NDTI (required for WCT classification). Optionally adds AWEIsh and AWEInsh.

Parameters:

ds (xr.Dataset or xr.DataArray) – Multi-spectral image (single-date or seasonal composite).
green_band (str) – Band names within ds. Defaults match the wetlandmapper harmonised naming convention.
red_band (str) – Band names within ds. Defaults match the wetlandmapper harmonised naming convention.
nir_band (str) – Band names within ds. Defaults match the wetlandmapper harmonised naming convention.
swir_band (str) – Band names within ds. Defaults match the wetlandmapper harmonised naming convention.
swir2_band (str) – SWIR2 band name. Required only when include_awei=True. Default "swir2".
blue_band (str) – Blue band name. Required only for AWEIsh when include_awei=True. Default "blue".
include_awei (bool) – If True, also compute AWEIsh and AWEInsh and include them in the returned Dataset. Default False for backward compatibility.

Returns:

Dataset with variables MNDWI, NDVI, NDTI (always), and optionally AWEIsh, AWEInsh.

Return type:

xr.Dataset

Examples

Standard WCT workflow:

>>> indices = compute_indices(ds, green_band="B3", red_band="B4",
...                           nir_band="B5", swir_band="B6")

Including AWEI for shadow-affected scenes:

>>> indices = compute_indices(ds, green_band="B3", red_band="B4",
...                           nir_band="B5", swir_band="B6",
...                           swir2_band="B7", blue_band="B2",
...                           include_awei=True)

wetlandmapper.indices.compute_mndwi(ds: Dataset | DataArray, green_band: str = 'green', swir_band: str = 'swir') → DataArray[source]

Compute the Modified Normalised Difference Water Index (MNDWI).

MNDWI = (Green - SWIR) / (Green + SWIR)

Values range from -1 to +1. Positive values indicate surface water.

Parameters:

ds (xr.Dataset or xr.DataArray) – Multi-spectral dataset containing at least a green and a SWIR band.
green_band (str) – Name of the green band (default: "green"). Common alternatives: "B3" (Landsat 8/9), "B03" (Sentinel-2).
swir_band (str) – Name of the SWIR band (default: "swir"). Common alternatives: "B6" (Landsat 8/9 SWIR1), "B11" (Sentinel-2).

Returns:

MNDWI values in the range [-1, 1], preserving input dimensions.

Return type:

xr.DataArray

References

Xu, H. (2006). Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. International Journal of Remote Sensing, 27(14), 3025–3033. https://doi.org/10.1080/01431160600589179

wetlandmapper.indices.compute_ndti(ds: Dataset | DataArray, red_band: str = 'red', green_band: str = 'green') → DataArray[source]

Compute the Normalised Difference Turbidity Index (NDTI).

NDTI = (Red - Green) / (Red + Green)

Higher values indicate more turbid (suspended-sediment-rich) water.

Parameters:

ds (xr.Dataset or xr.DataArray)
red_band (str) – Name of the red band (default: "red").
green_band (str) – Name of the green band (default: "green").

Returns:

NDTI values in the range [-1, 1].

Return type:

xr.DataArray

wetlandmapper.indices.compute_ndvi(ds: Dataset | DataArray, nir_band: str = 'nir', red_band: str = 'red') → DataArray[source]

Compute the Normalised Difference Vegetation Index (NDVI).

NDVI = (NIR - Red) / (NIR + Red)

Parameters:

ds (xr.Dataset or xr.DataArray)
nir_band (str) – Name of the NIR band (default: "nir"). Common alternatives: "B5" (Landsat 8/9), "B08" (Sentinel-2).
red_band (str) – Name of the red band (default: "red"). Common alternatives: "B4" (Landsat 8/9), "B04" (Sentinel-2).

Returns:

NDVI values in the range [-1, 1].

Return type:

xr.DataArray

wetlandmapper.indices.compute_ndwi(ds: Dataset | DataArray, green_band: str = 'green', nir_band: str = 'nir') → DataArray[source]

Compute the Normalised Difference Water Index (NDWI).

NDWI = (Green - NIR) / (Green + NIR)

Values range from -1 to +1. Positive values indicate surface water. Less sensitive to built-up land and vegetation than MNDWI.

Parameters:

ds (xr.Dataset or xr.DataArray) – Multi-spectral dataset containing at least a green and a NIR band.
green_band (str) – Name of the green band (default: "green"). Common alternatives: "B3" (Landsat 8/9), "B03" (Sentinel-2).
nir_band (str) – Name of the NIR band (default: "nir"). Common alternatives: "B5" (Landsat 8/9), "B08" (Sentinel-2).

Returns:

NDWI values in the range [-1, 1], preserving input dimensions.

Return type:

xr.DataArray

References

McFeeters, S. K. (1996). The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features. International Journal of Remote Sensing, 17(7), 1425–1432. https://doi.org/10.1080/01431169608948714

wetlandmapper.indices.compute_water_indices(ds: Dataset | DataArray, green_band: str = 'green', red_band: str = 'red', nir_band: str = 'nir', swir_band: str = 'swir', swir2_band: str = 'swir2', blue_band: str = 'blue') → Dataset[source]

Compute all available water indices and return as a Dataset.

Computes MNDWI, NDWI, AWEIsh, and AWEInsh for comprehensive water detection. Useful for comparing different water indices or selecting the best performing one for a specific study area.

Parameters:

ds (xr.Dataset or xr.DataArray) – Multi-spectral image (single-date or time series).
green_band (str) – Band names within ds.
red_band (str) – Band names within ds.
nir_band (str) – Band names within ds.
swir_band (str) – Band names within ds.
swir2_band (str) – SWIR2 band name. Default "swir2".
blue_band (str) – Blue band name. Default "blue".

Returns:

Dataset with variables MNDWI, NDWI, AWEIsh, AWEInsh.

Return type:

xr.Dataset

Examples

>>> water_indices = compute_water_indices(
...     ds,
...     blue_band="B2", green_band="B3", red_band="B4",
...     nir_band="B5", swir_band="B6", swir2_band="B7",  # Landsat 8/9
... )
>>> water_indices
<xarray.Dataset>
Dimensions:  (time: ..., y: ..., x: ...)
Data variables:
    MNDWI    ...
    NDWI     ...
    AWEIsh   ...
    AWEInsh  ...

Wetland Dynamics

Wetland dynamics classification and temporal aggregation utilities.

wetlandmapper.dynamics.aggregate_time(da: DataArray | Dataset, freq: str = 'annual', method: str = 'median') → DataArray | Dataset[source]

Temporally aggregate a multi-temporal xarray object before classification.

Reduces a time series to one composite per chosen period by computing a pixel-wise statistic within each period. Useful for:

Dynamics: produce annual composites from all available scenes rather than using every raw overpass.
WCT: produce monthly or seasonal composites, then classify each with classify_wct() / classify_wct_ema().

Parameters:

da (xr.DataArray or xr.Dataset) – Input data with a time dimension. Accepts both GEE-fetched and locally constructed objects.
freq ({"annual", "monthly", "seasonal", "all"}) –
Aggregation period:

"annual"
One composite per calendar year (resampled to year-end).

"monthly"
One composite per calendar month.

"seasonal"
One composite per meteorological season per year: DJF (Dec–Jan–Feb), MAM (Mar–Apr–May), JJA (Jun–Jul–Aug), SON (Sep–Oct–Nov). Uses pandas quarterly resampling anchored to December.

"all"
No aggregation — returns da unchanged.
method ({"median", "mean", "max", "min"}) – Pixel-wise statistic computed within each period. Default "median".

Returns:

Same type as da with a reduced time dimension (one step per period). For "seasonal", the time coordinate is labelled with the first day of each quarter (e.g., 2003-12-01 for DJF 2004).

Return type:

xr.DataArray or xr.Dataset

Raises:

ValueError – If freq or method is not one of the valid options.

Examples

Produce annual MNDWI composites from a dense time series:

>>> from wetlandmapper.dynamics import aggregate_time
>>> mndwi_annual = aggregate_time(mndwi_ts, freq="annual")
>>> dynamics = classify_dynamics(mndwi_annual, nYear=3)

Produce seasonal composites for WCT classification:

>>> from wetlandmapper import compute_indices, classify_wct_ema
>>> from wetlandmapper.dynamics import aggregate_time
>>> indices_ts = fetch(aoi, "2010-01-01", "2023-12-31",
...                    index=["MNDWI","NDVI","NDTI"])
>>> seasonal = aggregate_time(indices_ts, freq="seasonal")
>>> # Classify each season independently
>>> for t in seasonal.time:
...     wct = classify_wct_ema(seasonal.sel(time=t))

Notes

Pixels that are NaN (masked by cloud or no-data) in all scenes within a period remain NaN in the composite. The statistic is computed ignoring NaNs (skipna=True is the xarray default).

wetlandmapper.dynamics.classify_dynamics(water_index: DataArray, nYear: int = 3, thresholdWet: float = 25.0, thresholdPersis: float = 75.0, water_threshold: float = 0.0, nan_policy: str = 'total', min_valid_obs: int | None = None, mndwi_threshold: float | None = None) → DataArray[source]

Classify wetland pixels into six mutually exclusive temporal dynamics classes.

Each pixel is assigned to exactly one class based on three temporal summary statistics derived from a multi-year water index stack:

Overall wet frequency W% (% of years above water_threshold)
Historic wet count W_historic (first nYear years)
Recent wet count W_recent (last nYear years)

Class priority (strictly exclusive — a pixel receives exactly one class):

Each priority is applied only to unclassified pixels (code = 0), preventing any pixel from receiving more than one class code even when multiple conditions are simultaneously true (e.g. a pixel that is both Persistent and Intensifying will be classified as Persistent only).

Parameters:

water_index (xr.DataArray) – Multi-temporal water index time series with a time dimension. Accepts MNDWI, AWEIsh, AWEInsh, NDWI, or any index where positive values indicate surface water.
nYear (int) – Length of the historic and recent windows in years. Default 3.
thresholdWet (float) – Minimum wet frequency (%) for a pixel to be any wetland class. Default 25.
thresholdPersis (float) – Wet frequency (%) above which a pixel is Persistent. Must be greater than thresholdWet. Default 75.
water_threshold (float) – Index value above which a pixel is counted as wet in a given year. Default 0.0 (positive MNDWI = water-dominated; Xu 2006).
nan_policy ({"total", "valid"}) –
Denominator used when computing wet frequency.

"total" (default)
Denominator = total number of time steps. NaN pixels count as dry. Reproduces the original Singh & Sinha (2022) method. Appropriate when NaN values are rare or randomly distributed.

"valid"
Denominator = per-pixel count of non-NaN observations. Wet frequency is the fraction of cloud-free years that were wet. The historic/recent windows are also normalised by their own valid counts so that delta is expressed as a fraction in [−1, +1]. Use when cloud masking produces substantial or spatially clustered NaN values.
min_valid_obs (int, optional) – (Only active when nan_policy=”valid”.) Minimum number of non-NaN observations required for classification. Pixels with fewer valid observations are set to NaN in the output. Default None (no minimum enforced).
mndwi_threshold (float, optional) – Deprecated. Use water_threshold instead.

Returns:

Integer DataArray (dtype int8) of shape (y, x) with class codes from DYNAMICS_CLASSES. Guaranteed to contain only valid class codes — no additive artefacts.

Return type:

xr.DataArray

Raises:

ValueError – If water_index lacks a time dimension, if nYear * 2 > n_time, if thresholds are out of range, or if nan_policy is not one of the accepted values.

References

Singh, M. & Sinha, R. (2022). Remote Sensing Letters, 13(1), 1–13. https://doi.org/10.1080/2150704X.2021.1980919

wetlandmapper.dynamics.compute_wet_frequency(water_index: DataArray, water_threshold: float = 0.0, nan_policy: str = 'total', mndwi_threshold: float | None = None) → DataArray[source]

Return the pixel-wise wet frequency (%) across the full time series.

Parameters:

water_index (xr.DataArray) – Multi-temporal water index with a time dimension.
water_threshold (float) – Index value above which a pixel is counted as wet. Default 0.0.
nan_policy ({"total", "valid"}) – "total" (default): denominator = total time steps; NaN = dry. "valid": denominator = per-pixel non-NaN count.
mndwi_threshold (float, optional) – Deprecated. Use water_threshold instead.

Returns:

Wet frequency in percent (0–100), shape (y, x).

Return type:

xr.DataArray

Wetland Cover Types (WCT)

wct.py

Two classification approaches are provided, both based on:: Singh, M., Allaka, S., Gupta, P. K., Patel, J. G., & Sinha, R. (2022). Deriving wetland-cover types (WCTs) from integration of multispectral indices based on Earth observation data. Environmental Monitoring and Assessment, 194(12), 878. https://doi.org/10.1007/s10661-022-10541-7

Functions

classify_wct_ema(indices, n_parts=4)

Original EMA paper method. Each index is independently binned into n_parts equal intervals over its positive [0, 1] range (default 4, giving quartile step 0.25), producing a level 0-4 per index. Every pixel then carries a three-way combination code w_level × v_level × t_level (e.g. “w2v3t1”) and the WCT class is read from a pre-built 5×5×5 lookup table. No priority ordering: each unique combination has exactly one class.

Critical implication: a pixel with negative MNDWI (w0) but high NDVI (v3 or v4) is assigned WCT 4 (Emergent/Floating Veg) because dense aquatic vegetation completely suppresses the water signal in MNDWI.

classify_wct(indices, thresholds=None)

Improved continuous-threshold method. Same five WCT classes, but boundaries are free-floating values that can be tuned per sensor, season, or region without changing the bin structure.

WCT class codes (shared by both methods)

1 Open Clear Water High MNDWI, Low NDVI, Low NDTI 2 Turbid / Sediment-laden High MNDWI, Low NDVI, High NDTI 3 Submerged Aquatic Veg. High MNDWI, Mod. NDVI, Low NDTI 4 Emergent / Floating Veg. Any MNDWI, High NDVI, Low NDTI ← includes w0 5 Moist / Waterlogged Soil Low MNDWI, Low NDVI, Mod. NDTI 0 Non-wetland / Dry

wetlandmapper.wct.build_ema_lookup_table(n_parts: int = 4) → ndarray[source]

Build the (n_parts+1)³ lookup table mapping index levels → WCT class.

Every entry table[w, v, t] is the WCT class code (0–5) for the combination of MNDWI level w, NDVI level v, and NDTI level t.

Level coding (same for all three indices):: 0 — negative value (dry / non-wetland) 1 — [0, 1/n_parts) Low 2 — [1/n_parts, 2/n_parts) Moderate-Low … n — [(n-1)/n_parts, 1.0] High (n = n_parts)

WCT class rules (applied in this order; later entries overwrite earlier):

WCT MNDWI level NDVI level NDTI level Cover type — ———– ———- ———- ———-

5 w == 1 v <= 1 t <= 2 Moist / Waterlogged Soil 3 w >= 2 1 <= v <= 2 t <= 1 Submerged Aquatic Veg. 2 w >= 2 v <= 1 t >= 2 Turbid / Sediment Water 4 w >= 0 v >= 3 t <= 1 Emergent / Floating Veg. (*) 1 w >= 3 v <= 1 t <= 1 Open Clear Water

(*) WCT 4 matches even w == 0 (negative MNDWI): dense aquatic vegetation completely suppresses the MNDWI water signal, so the pixel reads as dry on MNDWI but shows high NDVI. The combination code w0v3t1 (for example) is unambiguously aquatic-vegetated, not dry land.

Parameters:: n_parts (int) – Number of equal bins over [0, 1]. Default 4.
Returns:: Lookup table indexed as [mndwi_level, ndvi_level, ndti_level].
Return type:: np.ndarray, shape (n_parts+1, n_parts+1, n_parts+1), dtype int8

wetlandmapper.wct.classify_wct(indices: Dataset, thresholds: dict[str, float] | None = None) → DataArray[source]

Classify Wetland Cover Types using continuous, adjustable spectral thresholds.

An improvement of classify_wct_ema(): boundaries are free-floating values (not locked to multiples of 1/n_parts), enabling sub-quartile calibration for different sensors, seasons, or geographic regions. The five WCT classes and their spectral logic are identical to the EMA method; only the boundary representation differs.

WCT 4 (Emergent/Floating Veg.) is assigned whenever NDVI exceeds ndvi_veg_high regardless of the MNDWI sign, matching the EMA method’s behaviour for vegetation-masked pixels.

Parameters:

indices (xr.Dataset) – Dataset with variables "MNDWI", "NDVI", "NDTI".
thresholds (dict, optional) –
Override any subset of the defaults (from Singh et al. 2022):

"mndwi_water" (0.0) MNDWI above which pixel is open water "mndwi_moist" (-0.2) MNDWI lower bound for moist-soil zone "ndvi_veg_high" (0.2) NDVI above which emergent veg. present "ndvi_veg_low" (0.05) NDVI lower bound for submerged veg. "ndti_turbid" (0.0) NDTI above which water is turbid

Returns:

Integer raster of WCT codes (dtype int8). Name: "wetland_cover_type".

Return type:

xr.DataArray

References

Singh et al. (2022). Environmental Monitoring and Assessment, 194(12), 878.

wetlandmapper.wct.classify_wct_ema(indices: Dataset, n_parts: int = 4) → Dataset[source]

Classify Wetland Cover Types using the original Singh et al. (2022) EMA method.

Each pixel is independently binned into a level (0–n_parts) for MNDWI, NDVI, and NDTI. The three levels together form a combination code:

w_level v_level t_level → WCT class ——- ——- ——- ———

0 4 0 → 4 (Emergent veg — water signal masked) 3 0 0 → 1 (Open clear water) 2 0 3 → 2 (Turbid water) 2 2 0 → 3 (Submerged aquatic veg.) 1 0 1 → 5 (Moist soil)

The mapping is a pre-built lookup table (see build_ema_lookup_table()). There is no priority ordering — each combination has exactly one class.

A vegetation-masked pixel (negative MNDWI but high NDVI) is classified as WCT 4 (Emergent / Floating Vegetation) even though its MNDWI alone would suggest dry land, because the dense canopy fully attenuates the water signal.

Parameters:

indices (xr.Dataset) – Dataset with variables "MNDWI", "NDVI", "NDTI". Typically produced by wetlandmapper.compute_indices().
n_parts (int) – Number of equal bins over the positive [0, 1] range. Default 4 (step 0.25; levels 1–4 at boundaries 0.25 / 0.50 / 0.75 / 1.00).

Returns:

Dataset with two variables:

"wetland_cover_type": Integer raster of WCT class codes (dtype int8). CRS preserved from input if rioxarray is available.
"combination_code": Integer raster encoding the three index levels that produced each WCT class. Hundreds digit = MNDWI level, tens = NDVI, units = NDTI (e.g. 401 → MNDWI High / NDVI Negative / NDTI Low → WCT 1). Useful for inspecting which index fractions drove the classification at each pixel.

Return type:

xr.Dataset

Notes

The combination code is distinct from the class: two pixels with the same WCT class may have different combination codes, revealing how far each index was from the bin boundaries.

Examples

>>> from wetlandmapper import compute_indices, classify_wct_ema
>>> indices = compute_indices(ds, green_band="B3", red_band="B4",
...                           nir_band="B5",  swir_band="B6")
>>> wct = classify_wct_ema(indices)
>>> wct.attrs['classification_method']
'EMA-combination-lookup'

References

Singh et al. (2022). Environmental Monitoring and Assessment, 194(12), 878. https://doi.org/10.1007/s10661-022-10541-7

Terrain Analysis

terrain.py

DEM-based terrain analysis for masking high-altitude artefacts in wetland mapping.

Glaciers, permanent snowpacks, and steep mountain terrain produce strong positive wetness index signals that are not wetlands. This module provides tools to identify and mask these artefacts using a Digital Elevation Model (DEM) by characterising local terrain flatness — the key topographic criterion that separates true wetlands (flat, low-lying areas) from high-altitude false positives.

Three flatness metrics are provided, each serving a different purpose

compute_slope(dem): Most broadly applicable. A wetland pixel almost never sits on a slope >5°. Computed from numpy gradient on the DEM with approximate metre conversion for geographic coordinates. Independent of window size.
compute_tpi(dem, window): Adds discrimination between flat plateau (near-zero TPI, high elevation) and valley bottom (negative TPI). Slope alone can’t separate a flat glacier from a flat floodplain. TPI can, if you combine it with an elevation ceiling.
compute_local_range(dem, window): Directly replicates the GEE rolling-window approach. Most intuitive to threshold (e.g. “< 30 m variation in a 5×5 window”). Maximum minus minimum elevation within an NxN window. A low local range indicates flat terrain.
map_dem_depressions(raw_dem, filled_dem, …): Depression mapping from raw and pit-filled DEM using integer division and binary reclassification (depression=1, non-depression=0), following the protocol described by Sinha et al. (2017, Current Science).
mask_terrain_artifacts(wetness, dem, …): Combines any or all three metrics plus an elevation ceiling. Emits a warning if the mask retains <10% of pixels (likely thresholds are too strict).

Notes

All functions operate on xarray DataArrays with spatial dimensions (y, x) or (lat, lon). CRS is preserved where rioxarray is available. No GEE dependency — these functions are designed for use on locally downloaded data (from wetlandmapper.gee.fetch() or any other source).

For server-side DEM masking within GEE (without downloading the DEM), use the dem_mask parameter of wetlandmapper.gee.fetch().

wetlandmapper.terrain.compute_local_range(dem: DataArray, window: int = 5) → DataArray[source]

Compute local elevation range (max - min) within a rolling window.

A low local range indicates flat terrain. This directly replicates the rolling-window approach used in GEE scripts to retain only flat neighbourhoods for wetland mapping.

Parameters:

dem (xr.DataArray) – Digital Elevation Model with spatial dimensions (y, x).
window (int) – Square window size in pixels. Default 5 (5 x 5 window).

Returns:

Local elevation range in the same units as dem (metres). Name: "local_range".

Return type:

xr.DataArray

wetlandmapper.terrain.compute_slope(dem: DataArray, units: str = 'degrees') → DataArray[source]

Compute terrain slope from a DEM DataArray.

Uses central-difference gradient (numpy.gradient) applied to the spatial dimensions. The result represents the steepness of terrain at each pixel — low values indicate flat ground suitable for wetlands.

Parameters:

dem (xr.DataArray) – Digital Elevation Model with spatial dimensions (y, x) or (lat, lon). Units should be metres.
units ({"degrees", "radians", "percent"}) – Output slope units. Default "degrees".

Returns:

Slope values, same spatial shape as dem. Name: "slope".

Return type:

xr.DataArray

wetlandmapper.terrain.compute_tpi(dem: DataArray, window: int = 5) → DataArray[source]

Compute the Topographic Position Index (TPI).

TPI = elevation - focal_mean(elevation, window x window)

Positive TPI indicates a pixel is higher than its surroundings (hilltop, ridge); negative TPI indicates a valley or depression; near-zero TPI indicates flat terrain or a mid-slope position.

Parameters:

dem (xr.DataArray) – Digital Elevation Model with spatial dimensions (y, x) or (lat, lon).
window (int) – Size of the square focal window in pixels. Odd numbers produce a symmetric neighbourhood. Default 5 (5 x 5 = 25-pixel window).

Returns:

TPI values in the same elevation units as dem (usually metres). Name: "TPI".

Return type:

xr.DataArray

wetlandmapper.terrain.map_dem_depressions(raw_dem: DataArray, filled_dem: DataArray, *, require_integer: bool = True, apply_cleanup: bool = True, cleanup_window: int = 3, min_neighbours: int = 2) → DataArray[source]

Map topographic depressions from raw and pit-filled DEMs.

This implements a depression protocol used in floodplain wetland mapping:

Integer-divide raw_dem / filled_dem.
Pixels with value 1 are unchanged terrain (no pit).
Pixels with value 0 are depressions (raw < filled).
Reclassify to a binary mask: depression=1, non-depression=0.

Parameters:

raw_dem (xr.DataArray) – Original (unfilled) DEM.
filled_dem (xr.DataArray) – Pit-filled DEM created from raw_dem.
require_integer (bool) – If True (default), both DEMs must have integer dtype.
apply_cleanup (bool) – If True (default), remove isolated one-pixel/very small speckles using a neighbourhood-count filter.
cleanup_window (int) – Square rolling window size for cleanup. Default 3.
min_neighbours (int) – Minimum number of depression pixels (including self) within cleanup_window to retain a depression pixel. Default 2.

Returns:

Binary depression mask with values {0, 1}. Name: "depression_mask".

Return type:

xr.DataArray

Notes

This method performs best in low-relief floodplains and may be less reliable in rugged terrain. Residual speckle can still occur and may need additional post-processing for specific study areas.

References

Sinha, R., Saxena, S., & Singh, M. (2017). Protocols for Riverine Wetland Mapping and Classification Using Remote Sensing and GIS. Current Science, 112(7), 1544-1552. http://www.jstor.org/stable/24912702

wetlandmapper.terrain.mask_terrain_artifacts(wetness: DataArray | Dataset, dem: DataArray, max_slope: float | None = 5.0, max_tpi: float | None = None, max_local_range: float | None = None, local_range_window: int = 5, tpi_window: int = 5, max_elevation: float | None = None, invert: bool = False) → DataArray | Dataset[source]

Mask wetness data using terrain flatness and elevation filters.

Parameters:

wetness (xr.DataArray or xr.Dataset) – Wetness or index data to mask. Can include a time dimension.
dem (xr.DataArray) – Digital Elevation Model with spatial dimensions matching wetness.
max_slope (float or None) – Maximum slope in degrees. Pixels steeper than this are masked. Default 5.0. Set to None to disable slope filtering.
max_tpi (float or None) – Maximum absolute TPI (metres). Default None (disabled).
max_local_range (float or None) – Maximum local elevation range (metres). Default None (disabled).
local_range_window (int) – Window size for local range. Default 5.
tpi_window (int) – Window size for TPI. Default 5.
max_elevation (float or None) – Absolute elevation ceiling (metres). Pixels above this are masked. Default None.
invert (bool) – If True, invert the mask so excluded terrain is retained.

Returns:

wetness masked by terrain suitability. The type matches the input.

Return type:

xr.DataArray or xr.Dataset

Plotting

plotting.py

All functions: - Respect geographic coordinates when present (x/y or lon/lat),

so axes ticks show real coordinate values rather than pixel indices.

Detect y-axis direction (ascending vs descending) and set origin accordingly so images are never vertically flipped.
Place the class legend outside the plot area to avoid obscuring data.
Return (fig, ax) so callers can further customise or save.

Coordinate conventions

wetlandmapper.gee.fetch(): Returns arrays with dims (time, y, x) where y is descending (north → south, standard raster convention).
wetlandmapper.gee.fetch_xee(): After the built-in sortby fix, also returns y descending. If you have an older array with lat/lon dims, rename them first: da = da.rename({"lat": "y", "lon": "x"}).

wetlandmapper.plotting.plot_dynamics(dynamics, ax=None, title: str = 'Wetland Dynamics', figsize: tuple = (8, 7), add_colorbar: bool = True, legend_loc: str = 'outside right', savepath: str | None = None, dpi: int = 150)[source]

Plot a wetland dynamics classification raster.

Parameters:

dynamics (xr.DataArray) – Output of wetlandmapper.classify_dynamics(). Spatial coordinates (x/y or lon/lat) are used for axis ticks when present.
ax (matplotlib.axes.Axes, optional) – Axes to draw into. Created if not provided.
title (str) – Plot title.
figsize (tuple) – Figure size in inches.
add_colorbar (bool) – Add a class legend.
legend_loc (str) – Legend placement. Use "outside right" (default), "outside bottom", or any standard Matplotlib loc string (e.g. "lower right"). "outside right" / "outside bottom" never overlap the data.
savepath (str, optional) – If given, save the figure to this path (PNG, PDF, TIFF, etc.).
dpi (int) – Resolution for saved figure. Default 150.

Returns:

fig, ax

Return type:

matplotlib Figure and Axes

wetlandmapper.plotting.plot_index(da, index_name: str = 'Index', ax=None, figsize: tuple = (8, 7), vmin: float = -1.0, vmax: float = 1.0, cmap: str = 'RdYlGn', time_step: int | None = None, savepath: str | None = None, dpi: int = 150)[source]

Plot a single spectral index (MNDWI, NDVI, or NDTI).

Parameters:

da (xr.DataArray) – 2-D or 3-D (time, y, x) index DataArray.
index_name (str) – Used in the title and colorbar label.
time_step (int, optional) – Select this time index (0-based) when da has a time dim. Defaults to the temporal mean if not provided.
savepath – Same as plot_dynamics().
dpi – Same as plot_dynamics().

Return type:

fig, ax

wetlandmapper.plotting.plot_wct(wct, ax=None, title: str = 'Wetland Cover Types', figsize: tuple = (8, 7), add_colorbar: bool = True, legend_loc: str = 'outside right', savepath: str | None = None, dpi: int = 150)[source]

Plot a Wetland Cover Type classification raster.

Parameters:

wct (xr.DataArray) – Output of wetlandmapper.classify_wct() or wetlandmapper.classify_wct_ema().
ax – Same as plot_dynamics().
title – Same as plot_dynamics().
figsize – Same as plot_dynamics().
add_colorbar – Same as plot_dynamics().
legend_loc – Same as plot_dynamics().
savepath – Same as plot_dynamics().
dpi – Same as plot_dynamics().

Return type:

fig, ax

wetlandmapper.plotting.plot_wet_frequency(mndwi, ax=None, figsize: tuple = (8, 7), mndwi_threshold: float = 0.0, savepath: str | None = None, dpi: int = 150)[source]

Plot wet frequency (%) derived from an MNDWI time series.

Parameters:

mndwi (xr.DataArray) – Multi-temporal MNDWI with a time dimension.
mndwi_threshold (float) – Pixels with MNDWI above this value are counted as wet.
savepath – Same as plot_dynamics().
dpi – Same as plot_dynamics().

Return type:

fig, ax

Google Earth Engine

gee.py

Supports all Landsat missions (4, 5, 7, 8, 9), Sentinel-2, and MODIS (Terra/Aqua), including a "LandsatAll" option that automatically merges available missions for any requested date range.

Two retrieval functions

fetch(aoi, …): Downloads data immediately via GEE’s getDownloadURL API. Best for small-to-medium AOIs (up to ~100 km²).
fetch_xee(aoi, …): Opens the GEE collection as a lazy Dask-backed xarray via xee. No data is transferred until .compute() is called. Suitable for large AOIs and long time series.

Temporal aggregation

Both functions accept a temporal_aggregation parameter:

"all" — every individual scene (default) "annual" — one median composite per calendar year "monthly" — one median composite per calendar month "seasonal" — one composite per meteorological season (DJF/MAM/JJA/SON)

Sensors and date coverage

"Landsat" is an alias for "Landsat8" for backward compatibility.

Landsat 7 SLC-off note

The Scan Line Corrector (SLC) on Landsat 7 ETM+ failed on 2003-05-31. Images acquired after this date have wedge-shaped data gaps covering roughly 22 % of each scene. Use use_slc_off=False (default) to exclude these images and use only the good-quality 1999–2003 record. Set use_slc_off=True to include post-failure images (useful when other sensors have no coverage, e.g. 1999–2012 before Landsat 8).

MODIS note

MODIS provides 500m resolution surface reflectance composites (8-day intervals). Use "MODISAll" to automatically merge Terra and Aqua for continuous coverage. MODIS has coarser resolution but longer temporal record (2000–present) compared to Landsat. Suitable for regional-scale studies where 500m resolution is adequate. Band mapping differs from Landsat — MODIS stores Red as Band 1, NIR as Band 2, Blue as Band 3, Green as Band 4. Cloud masking uses StateQA bits 0-1 (cloud state) and bit 2 (shadow). Scale factor is 0.0001 with no offset (unlike Landsat C02 L2’s -0.2 offset). AWEIsh and AWEInsh are computed server-side for MODIS since all required bands are available.

DEM masking

Server-side DEM masking using Copernicus GLO-30 can be activated with the dem_mask parameter in fetch(). This applies terrain filters directly in GEE using ee.Terrain.slope() and reduceNeighborhood, avoiding the need to download the DEM. The snow question is addressed by min_temp_c in climate-adaptive mode (ERA5 precipitation includes snowfall; filtering to months ≥5°C removes cold months where precipitation is frozen) and by max_elevation_m (above the local glaciation line, always mask regardless of other criteria).

Requirements

earthengine-api : always required
rasterio : required by fetch()
xee + dask : required by fetch_xee()

Install: pip install wetlandmapper[gee] Authenticate once: earthengine authenticate

wetlandmapper.gee.authenticate() → None[source]: Run interactive GEE authentication (opens browser; only needed once).

wetlandmapper.gee.fetch(aoi: dict | str, start: str, end: str, sensor: str = 'Landsat8', index: str | Sequence[str] = 'MNDWI', custom_indices: dict[str, str] | None = None, scale: int = 30, max_cloud_cover: float = 20.0, temporal_aggregation: str = 'all', use_slc_off: bool = False, project: str | None = None, climate_adaptive: bool = False, min_precip_mm: float = 20.0, min_temp_c: float = 5.0, hydroperiod_months: int = 1, wetness_index: str = 'MNDWI', wetness_threshold: float = 0.0, dem_mask: bool = False, max_slope_deg: float | None = 5.0, max_tpi_m: float | None = None, tpi_window_px: int = 5, max_local_range_m: float | None = None, local_range_window_px: int = 5, max_elevation_m: float | None = None, months: list[int] | None = None, reduction_method: str = 'median', percentile: float = 50.0) → DataArray | Dataset[source]

Retrieve spectral indices from GEE as an xarray object (immediate download).

Fetches an image collection, applies cloud masking and surface-reflectance scaling, computes MNDWI / NDVI / NDTI server-side, optionally composites by time period, and downloads the result.

Parameters:

aoi (dict, str, or Path) –
Area of interest. Accepts:
- GeoJSON dict — plain geometry, Feature, or FeatureCollection.
- Shapefile path — str / Path to a .shp (or any format readable by geopandas). Multiple features are dissolved into one boundary. Requires pip install geopandas.
- GeoJSON file path — a .geojson / .json file.
start (str) – Start date ISO 8601, e.g. "2000-01-01".
end (str) – End date (inclusive) ISO 8601, e.g. "2023-12-31".
sensor (str) – Satellite sensor. One of: "Landsat4", "Landsat5", "Landsat7", "Landsat8" (default), "Landsat9", "LandsatAll", "Sentinel2", "MODIS_Terra", "MODIS_Aqua", "MODISAll". "Landsat" is an alias for "Landsat8".
index (str or list of str) – One or more of "MNDWI", "NDWI", "NDVI", "NDTI", "AWEIsh", "AWEInsh". Single str → DataArray; list → Dataset (both with time dim).
custom_indices (dict[str, str], optional) –
User-defined index formulas evaluated server-side via ee.Image.expression. Dictionary keys are output band names and values are expression strings. Available formula symbols are blue, green, red, nir, swir (alias swir1), swir2, and qa.

Example: {"NDSI": "(green - swir) / (green + swir)"}. To request a custom index, include its name in index.
scale (int) – Spatial resolution in metres. Default 30 (Landsat native).
max_cloud_cover (float) – Maximum cloud cover (%) per image. Default 20.
temporal_aggregation ({"all", "annual", "monthly", "seasonal"}) – Server-side temporal compositing. Default "all" (every scene). Using "annual" or "monthly" greatly reduces download volume.
use_slc_off (bool) – Include Landsat 7 SLC-off images (acquired after 2003-05-31)? Default False. Only relevant when sensor is "Landsat7" or "LandsatAll".
project (str, optional) – GEE cloud project ID.
climate_adaptive (bool) – If True, replace the standard temporal composite with a climate-adaptive annual composite guided by ERA5-Land precipitation and temperature. When enabled, temporal_aggregation is ignored (output is always one image per year). Default False.
months (list of int, optional) – Restrict compositing to these calendar months (1=Jan, …, 12=Dec). Example: [6, 7, 8] for June–August only. Default None (all months).
reduction_method ({"median", "mean", "percentile"}) – Collection reducer applied within each temporal aggregation window. Default "median".
percentile (float) –
Percentile to use when reduction_method="percentile". Must be between 0 and 100 inclusive. Default 50.0.

When climate_adaptive=True, the algorithm:
1. Builds monthly Landsat composites internally.
2. Joins with ERA5-Land monthly precipitation and 2m temperature.
3. Filters months where precip >= min_precip_mm AND temp >= min_temp_c (excludes dry season and snow months).
4. For each year, selects per-pixel values from the month with peak precipitation using qualityMosaic (captures maximum wetness rather than an arbitrary median).
5. Masks pixels wet for fewer than hydroperiod_months months per year on average (removes transient waterlogging).
min_precip_mm (float) – Minimum monthly precipitation (mm) to include a month in the composite window. Used only when climate_adaptive=True. Default 20 mm.
min_temp_c (float) – Minimum monthly mean 2m temperature (degrees C) to include a month. Excludes frozen-ground and snow months. Used only when climate_adaptive=True. Default 5 degrees C.
hydroperiod_months (int) – Minimum number of months per year a pixel must be wet (index above wetness_threshold) on average across the full record to be retained. Pixels below this are masked as transient waterlogging. Used only when climate_adaptive=True. Default 1. Increase to 2-3 for stricter wetland delineation.
wetness_index (str) – Which index band to use as the wetness indicator for hydroperiod counting and qualityMosaic selection. Must be one of the bands in index. Default "MNDWI".
wetness_threshold (float) – Index value above which a pixel is counted as wet for the hydroperiod calculation. Default 0.0.
dem_mask (bool) – If True, apply a server-side terrain flatness mask using the Copernicus GLO-30 DEM before compositing. Masks out glaciers, snowpacks, and steep mountain terrain that produce false wetness signals. Default False.
max_slope_deg (float or None) – Maximum terrain slope (degrees) to retain when dem_mask=True. Default 5.0. Set to None to disable slope filtering.
max_tpi_m (float or None) – Maximum absolute TPI (metres) when dem_mask=True. Default None (disabled).
tpi_window_px (int) – Focal window radius in pixels for TPI. Default 5.
max_local_range_m (float or None) – Maximum local elevation range (metres) in the rolling window when dem_mask=True. Default None (disabled).
local_range_window_px (int) – Window radius for local elevation range. Default 5.
max_elevation_m (float or None) – Absolute elevation ceiling (metres). Pixels above this elevation are always masked when dem_mask=True. Default None.

Returns:

xr.DataArray – DataArray with dims (time, y, x) when a single index is requested.
xr.Dataset – Dataset with one variable per index, dims (time, y, x).

Notes

Time steps where no valid cloud-free pixels exist are skipped with a UserWarning rather than raising an error. The returned object will have fewer time steps than requested periods in such cases.

Examples

Long-record annual MNDWI for dynamics using all available Landsat missions:

>>> mndwi = fetch(aoi, "1984-01-01", "2023-12-31",
...               sensor="LandsatAll",
...               temporal_aggregation="annual")
>>> dynamics = classify_dynamics(mndwi, nYear=3)

Post-monsoon WCT composite from Landsat 5 era:

>>> indices = fetch(aoi, "2005-10-01", "2005-12-31",
...                 sensor="Landsat5",
...                 index=["MNDWI", "NDVI", "NDTI"])

Annual MNDWI from MODIS for regional-scale analysis:

>>> mndwi_modis = fetch(aoi, "2000-01-01", "2023-12-31",
...                     sensor="MODISAll",
...                     temporal_aggregation="annual",
...                     scale=500)

wetlandmapper.gee.fetch_xee(aoi: dict | str, start: str, end: str, sensor: str = 'Landsat8', index: str | Sequence[str] = 'MNDWI', custom_indices: dict[str, str] | None = None, scale: int = 30, max_cloud_cover: float = 20.0, temporal_aggregation: str = 'all', use_slc_off: bool = False, project: str | None = None, climate_adaptive: bool = False, min_precip_mm: float = 20.0, min_temp_c: float = 5.0, hydroperiod_months: int = 1, wetness_index: str = 'MNDWI', wetness_threshold: float = 0.0, dem_mask: bool = False, max_slope_deg: float | None = 5.0, max_tpi_m: float | None = None, tpi_window_px: int = 5, max_local_range_m: float | None = None, local_range_window_px: int = 5, max_elevation_m: float | None = None, chunks: dict | None = None, months: list[int] | None = None, reduction_method: str = 'median', percentile: float = 50.0) → DataArray | Dataset[source]

Retrieve spectral indices from GEE as a lazy Dask-backed xarray via xee.

Parameters:

aoi (dict, str, or Path) – GeoJSON dict, shapefile path, or GeoJSON file path.
start (str) – ISO 8601 date strings.
end (str) – ISO 8601 date strings.
sensor (str) – Sensor key — see fetch() for the full list. Default "Landsat8".
index (str or list of str) – One or more of "MNDWI", "NDWI", "NDVI", "NDTI", "AWEIsh", "AWEInsh". Single str → DataArray; list → Dataset.
custom_indices (dict[str, str], optional) –
User-defined index formulas evaluated server-side via ee.Image.expression. Dictionary keys are output band names and values are expression strings. Available formula symbols are blue, green, red, nir, swir (alias swir1), swir2, and qa.

Example: {"NDSI": "(green - swir) / (green + swir)"}. To request a custom index, include its name in index.
scale (int) – Pixel resolution in metres. Default 30.
max_cloud_cover (float) – Maximum per-image cloud cover (%). Default 20.
temporal_aggregation ({"all", "annual", "monthly", "seasonal"}) –
Server-side compositing mode. Default "all".

Warning

"all" passes the raw scene collection to xee with no compositing. Older xee versions may return integer time indices (0, 1, 2 …) instead of real timestamps for this mode. This function detects the issue, emits a UserWarning, and patches the time coordinate automatically by querying system:time_start from GEE. Use "annual" or "monthly" to avoid the extra round-trip.
use_slc_off (bool) – Include Landsat 7 SLC-off scenes. Default False.
project (str, optional) – GEE cloud project ID.
climate_adaptive
min_precip_mm
min_temp_c
hydroperiod_months

:param : :param wetness_index: :param wetness_threshold: :param dem_mask: :param max_slope_deg: :param : :param max_tpi_m: :param tpi_window_px: :param max_local_range_m: :param : :param local_range_window_px: Same semantics as fetch(). :param max_elevation_m: Same semantics as fetch(). :param chunks: Dask chunk sizes, e.g. {"time": 1, "lon": 512, "lat": 512}. :type chunks: dict, optional :param months: Restrict compositing to these calendar months (1=Jan, …, 12=Dec).

Example: [6, 7, 8] for June–August only. Default None (all months).

Parameters:

reduction_method ({"median", "mean", "percentile"}) – Collection reducer applied within each temporal aggregation window. Default "median".
percentile (float) – Percentile to use when reduction_method="percentile". Must be between 0 and 100 inclusive. Default 50.0.

Returns:

Lazy object with dims (time, lat, lon).

Return type:

xr.DataArray or xr.Dataset

Notes

After calling .compute(), orient dimensions with:

da = da.rename({"lat": "y", "lon": "x"}).transpose("time", "y", "x")

Examples

>>> mndwi_lazy = fetch_xee(
...     aoi, "1984-01-01", "2023-12-31",
...     sensor="LandsatAll", temporal_aggregation="annual",
... )
>>> mndwi = (
...     mndwi_lazy
...     .rename({"lat": "y", "lon": "x"})
...     .transpose("time", "y", "x")
...     .compute()
... )

wetlandmapper.gee.init(project: str | None = None) → None[source]

Initialise the GEE Python client.

Parameters:: project (str, optional) – GEE cloud project ID. Required for accounts created after 2023.