diff --git a/docs/processing.md b/docs/processing.md index 7b32a5d..9fa285a 100644 --- a/docs/processing.md +++ b/docs/processing.md @@ -19,10 +19,8 @@ register_glint_method("my_method", my_glint_method) ## Implemented methods - **Glint removal**: `"hedley"` (Hedley et al., 2005). -- **Depth correction**: `"lyzenga"` (Lyzenga-style log transform; the `variant` keyword accepts - `"1978"`/`"1981"`/`"2006"` but currently applies the same transform for all three — see - [Project status](status.md)), `"stumpf"` (ratio transform), `"maritorena"` (single-parameter - exponential attenuation). +- **Depth correction**: `"lyzenga1978"`, `"lyzenga1981"`, `"lyzenga2006"` (the default `"lyzenga"` is an alias for `"lyzenga2006"`), `"stumpf"` (ratio + transform, Stumpf et al., 2003), `"maritorena"` (single-parameter exponential attenuation). - **Cloud masking**: `"sentinel2_qa60"`, `"sentinel2_scl"`, `"landsat_qa_pixel"`, `"probability"` (threshold on a cloud-probability layer), `"user_mask"` (pass through a caller-supplied mask). `cloud_mask(scene, method="auto")` resolves one of these automatically from `scene.sensor` and @@ -31,3 +29,92 @@ register_glint_method("my_method", my_glint_method) `atmospheric_correction` is **intentionally deferred** — it raises `NotImplementedError` until a backend is implemented in a follow-up change. Unlike glint/depth/cloud processing, it does not yet have its own registry. + +## Lyzenga depth-correction variants + +There are three Lyzenga algorithms for depth correction - `"lyzenga1978"`, `"lyzenga1981"`, +`"lyzenga2006"`. `method="lyzenga"` is an +alias for `"lyzenga2006"`, the most recent and most widely used variant. + +### `"lyzenga1978"` — per-band log-linearising transform + +The single-band transform from Lyzenga (1978), eqs. (1) & (7). Replaces every spectral band with + +```text +X_i = ln(L_i - L_si) +``` + +which is approximately linear in water depth. Non-spectral bands (`qa60`, `scl`, `qa_pixel`) pass +through unchanged. Bands are not combined with each other, so this variant does not correct for +bottom-reflectance variation — it's the foundational linearisation the other two variants build on. + +```python +depth_correct(scene, method="lyzenga1978", deep_water_radiance={"blue": 0.02}) +``` + +- `deep_water_radiance` (mapping or scalar, default `0.0` per band) — the deep-water radiance + `L_si` to subtract before taking the log. +- `epsilon` (default `1e-6`) — floor applied to `L_i - L_si` before the log, to avoid `log(0)` + or `log(negative)`. + +### `"lyzenga1981"` — two-band depth-invariant bottom index + +The classic two-band index from Lyzenga (1981), eq. (2). Adds a `lyzenga_dii` band: + +```text +Y_i = [K_j * X_i - K_i * X_j] / sqrt(K_i^2 + K_j^2) +``` + +where `X = ln(L - L_s)` and `K_i`/`K_j` are the water-attenuation coefficients for the two chosen +bands. The result is (ideally) independent of water depth, so it's a bottom-type index rather than +a depth estimate — it removes depth dependence so bottom features can be compared across varying +water depth. + +```python +depth_correct( + scene, + method="lyzenga1981", + band_i="blue", + band_j="green", + k_i=0.056, + k_j=0.074, +) +``` + +- `band_i`, `band_j` (default `"blue"`, `"green"`) — the two bands to combine. +- `k_i`, `k_j` (**required**) — water-attenuation coefficients for `band_i`/`band_j`. Raises + `ValueError` if either is missing, or if both are zero (undefined). +- `deep_water_radiance`, `epsilon` — as in `"lyzenga1978"`. + +### `"lyzenga2006"` — empirical multi-band depth regression + +The regression-based estimator from Lyzenga et al. (2006), eq. (9). Adds a `lyzenga_depth` band +containing an actual depth estimate (in the same units as the calibration data, typically metres): + +```text +depth = h_0 - sum(h_j * X_j) +``` + +where `X_j = ln(L_j - L_sj)` and `h_j` are empirically fitted coefficients (`h_0` is the +intercept). Unlike the 1978/1981 variants, this produces a direct depth estimate rather than a +depth-invariant index — the coefficients are normally obtained by regressing against ground-truth +depth measurements (for example LIDAR) for the site being processed. + +```python +depth_correct( + scene, + method="lyzenga2006", # or method="lyzenga" — it's an alias for this variant + bands=("blue", "green"), + coefficients={"blue": -17.42, "green": 26.7}, + intercept=17.84, +) +``` + +- `bands` (default `("blue", "green")`) — the bands to combine; must be non-empty. +- `coefficients` (**required**) — mapping of band name to its fitted regression weight `h_j`. + Raises `ValueError` if missing, or if any band in `bands` has no entry. +- `intercept` (default `0.0`) — the `h_0` term. +- `deep_water_radiance`, `epsilon` — as in `"lyzenga1978"`. + +None of the three variants take a known `depth` as input — the whole point is to derive depth (or +a depth-invariant index) from radiance, unlike `"maritorena"` which requires one. diff --git a/docs/status.md b/docs/status.md index 866a70b..20a234c 100644 --- a/docs/status.md +++ b/docs/status.md @@ -11,8 +11,6 @@ scaffolded but not yet built — see [Getting started](getting-started.md) and data the `Scene` docstring describes as the typical case. - Spectral indices treat only Python `None` as a missing value; there is no nodata/fill-value (for example `-9999`) handling yet. -- `depth_correct(..., method="lyzenga", variant=...)` records the requested variant in scene - metadata/history but does not yet apply distinct 1978/1981/2006 formulas. ## Scaffolded, not yet implemented diff --git a/src/earthrs/processing/core.py b/src/earthrs/processing/core.py index a754ce5..e9d1668 100644 --- a/src/earthrs/processing/core.py +++ b/src/earthrs/processing/core.py @@ -2,7 +2,6 @@ from __future__ import annotations -import math from typing import Any from earthrs.scene import Scene @@ -117,10 +116,6 @@ def _map_binary(left: Any, right: Any, func) -> Any: return func(float(left), float(right)) -def _safe_log(value: float) -> float: - return math.log(max(value, 1e-6)) - - def _flatten_numeric(values: Any) -> list[float]: if isinstance(values, (list, tuple)): out: list[float] = [] diff --git a/src/earthrs/processing/depth.py b/src/earthrs/processing/depth.py index d948d39..3f0e4d2 100644 --- a/src/earthrs/processing/depth.py +++ b/src/earthrs/processing/depth.py @@ -3,17 +3,19 @@ from __future__ import annotations import math +from collections.abc import Iterable, Mapping from typing import Any from earthrs.processing.core import ( _DEPTH_REGISTRY, _map_binary, _map_unary, - _safe_log, register_depth_method, ) from earthrs.scene import Scene +_NON_SPECTRAL_BANDS = {"qa60", "scl", "qa_pixel"} + def depth_correct( scene: Scene, @@ -24,7 +26,10 @@ def depth_correct( ) -> Scene: """Apply depth correction to raster data. - Method names are resolved from the depth-correction registry. + Method names are resolved from the depth-correction registry. There are three + Lyzenga formulations, each denoted using the year of publication -- + ``"lyzenga1978"``, ``"lyzenga1981"``, ``"lyzenga2006"``. ``"lyzenga"`` is an + alias for ``"lyzenga2006"``, the most recent and most widely used variant. """ processor = _DEPTH_REGISTRY.get(method.lower()) @@ -33,44 +38,181 @@ def depth_correct( return processor(scene, depth=depth, **kwargs) -def _lyzenga_depth( - scene: Scene, - *, - depth: Any | None = None, - variant: str = "1981", - epsilon: float = 1e-6, - **_: Any, -) -> Scene: +def _require_mapping_data(scene: Scene, label: str) -> None: if not isinstance(scene.data, dict): - raise TypeError("Lyzenga correction expects mapping-based scene data.") - updated = {} - for band, values in scene.data.items(): - if band in {"qa60", "scl", "qa_pixel"}: - updated[band] = values - continue - corrected = _map_unary(values, lambda value: -math.log(max(value, epsilon))) - if depth is not None: - corrected = _map_binary( - corrected, depth, lambda value, depth_value: value / max(depth_value, epsilon) - ) - updated[band] = corrected + raise TypeError(f"{label} correction expects mapping-based scene data.") + + +def _band_param( + bands: Iterable[str], value: Any, default: float +) -> dict[str, float]: + """Normalise a per-band parameter to a ``{band: value}`` mapping. + + Accepts ``None`` (every band gets `default`), a single scalar (every band gets + that value), or a mapping (missing bands fall back to `default`). + """ + + band_list = list(bands) + if value is None: + return dict.fromkeys(band_list, default) + if isinstance(value, Mapping): + return {band: float(value.get(band, default)) for band in band_list} + return dict.fromkeys(band_list, float(value)) + + +def _log_transform(values: Any, deep_water_radiance: float, epsilon: float) -> Any: + """Apply the Lyzenga log-linearising transform ``X = ln(L - L_s)`` elementwise.""" + + return _map_unary(values, lambda v: math.log(max(v - deep_water_radiance, epsilon))) + + +def _finalise_lyzenga( + scene: Scene, *, data: dict[str, Any], band_names: tuple[str, ...], method: str +) -> Scene: + """Build the result `Scene` shared by all Lyzenga variants.""" + metadata = dict(scene.metadata) - metadata["depth_method"] = f"lyzenga_{variant}" + metadata["depth_method"] = method return Scene( - data=updated, + data=data, crs=scene.crs, transform=scene.transform, metadata=metadata, - band_names=scene.band_names, + band_names=band_names, acquisition_time=scene.acquisition_time, sensor=scene.sensor, - history=(*scene.history, f"depth_correct:lyzenga:{variant}"), + history=(*scene.history, f"depth_correct:{method}"), masks=dict(scene.masks), cloud_mask=scene.cloud_mask, cloud_probability=scene.cloud_probability, ) +def _lyzenga_1978( + scene: Scene, + *, + depth: Any | None = None, + deep_water_radiance: Mapping[str, float] | float | None = None, + epsilon: float = 1e-6, + **_: Any, +) -> Scene: + """Lyzenga (1978) per-band log-linearising transform. + + Replaces each spectral band with ``X_i = ln(L_i - L_si)`` (eqs. 1 & 7), which + is approximately linear in water depth. Non-spectral bands (``qa60``, ``scl``, + ``qa_pixel``) pass through unchanged. Bands are not combined with each other. + """ + + _ = depth + _require_mapping_data(scene, "Lyzenga 1978") + dwr = _band_param(scene.data.keys(), deep_water_radiance, 0.0) + updated = {} + for band, values in scene.data.items(): + if band in _NON_SPECTRAL_BANDS: + updated[band] = values + continue + updated[band] = _log_transform(values, dwr[band], epsilon) + return _finalise_lyzenga( + scene, data=updated, band_names=scene.band_names, method="lyzenga1978" + ) + + +def _lyzenga_1981( + scene: Scene, + *, + depth: Any | None = None, + band_i: str = "blue", + band_j: str = "green", + k_i: float | None = None, + k_j: float | None = None, + deep_water_radiance: Mapping[str, float] | float | None = None, + epsilon: float = 1e-6, + **_: Any, +) -> Scene: + """Lyzenga (1981) two-band depth-invariant bottom index. + + Adds a ``lyzenga_dii`` band computed from eq. (2), + ``Y_i = [K_j * X_i - K_i * X_j] / sqrt(K_i^2 + K_j^2)``, where + ``X = ln(L - L_s)`` and `k_i`/`k_j` are the water-attenuation coefficients for + `band_i`/`band_j`. The result is (ideally) independent of water depth, so it can + be used as a bottom-type index rather than an estimate of depth itself. + """ + + _ = depth + _require_mapping_data(scene, "Lyzenga 1981") + if k_i is None or k_j is None: + raise ValueError( + "Lyzenga 1981 depth-invariant index requires water-attenuation " + "coefficients `k_i` and `k_j`." + ) + missing = [band for band in (band_i, band_j) if band not in scene.data] + if missing: + raise ValueError(f"Lyzenga 1981 requires bands {missing} in scene data.") + denominator = math.sqrt(k_i**2 + k_j**2) + if denominator == 0: + raise ValueError("Lyzenga 1981 requires `k_i` and `k_j` not both be zero.") + dwr = _band_param((band_i, band_j), deep_water_radiance, 0.0) + x_i = _log_transform(scene.data[band_i], dwr[band_i], epsilon) + x_j = _log_transform(scene.data[band_j], dwr[band_j], epsilon) + dii = _map_binary(x_i, x_j, lambda xi, xj: (k_j * xi - k_i * xj) / denominator) + updated = dict(scene.data) + updated["lyzenga_dii"] = dii + band_names = (*scene.band_names, "lyzenga_dii") + return _finalise_lyzenga(scene, data=updated, band_names=band_names, method="lyzenga1981") + + +def _lyzenga_2006( + scene: Scene, + *, + depth: Any | None = None, + bands: tuple[str, ...] = ("blue", "green"), + coefficients: Mapping[str, float] | None = None, + intercept: float = 0.0, + deep_water_radiance: Mapping[str, float] | float | None = None, + epsilon: float = 1e-6, + **_: Any, +) -> Scene: + """Lyzenga et al. (2006) empirical multi-band depth regression. + + Adds a ``lyzenga_depth`` band computed from eq. (9), + ``depth = h_0 - sum(h_j * X_j)``, where ``X_j = ln(L_j - L_sj)`` and `h_j` are + empirically fitted `coefficients` (`intercept` supplies ``h_0``). Unlike the + 1978/1981 variants, this produces a direct depth estimate rather than a + depth-invariant index. This is the most recent and most widely used Lyzenga + variant, so it is also registered as the plain ``"lyzenga"`` method. + """ + + _ = depth + _require_mapping_data(scene, "Lyzenga 2006") + if not bands: + raise ValueError("Lyzenga 2006 depth regression requires at least one band.") + if coefficients is None: + raise ValueError( + "Lyzenga 2006 depth regression requires `coefficients` (empirically " + "fitted h_j weights per band, see eq. (9))." + ) + missing_bands = [band for band in bands if band not in scene.data] + if missing_bands: + raise ValueError(f"Lyzenga 2006 requires bands {missing_bands} in scene data.") + missing_coeff = [band for band in bands if band not in coefficients] + if missing_coeff: + raise ValueError(f"Missing regression coefficients for bands {missing_coeff}.") + dwr = _band_param(bands, deep_water_radiance, 0.0) + depth_values: Any = None + for band in bands: + x_band = _log_transform(scene.data[band], dwr[band], epsilon) + term = _map_unary(x_band, lambda xv, h=coefficients[band]: -h * xv) + if depth_values is None: + depth_values = term + else: + depth_values = _map_binary(depth_values, term, lambda a, b: a + b) + depth_values = _map_unary(depth_values, lambda v: v + intercept) + updated = dict(scene.data) + updated["lyzenga_depth"] = depth_values + band_names = (*scene.band_names, "lyzenga_depth") + return _finalise_lyzenga(scene, data=updated, band_names=band_names, method="lyzenga2006") + + def _stumpf_depth( scene: Scene, *, @@ -81,6 +223,11 @@ def _stumpf_depth( m1: float = 1.0, **_: Any, ) -> Scene: + """Stumpf (2003) band-ratio depth transform. + + Adds a ``stumpf_depth`` band from ``m0 - m1 * ln(blue) / ln(green)``. + """ + _ = depth if not isinstance(scene.data, dict): raise TypeError("Stumpf correction expects mapping-based scene data.") @@ -89,7 +236,7 @@ def _stumpf_depth( stumpf = _map_binary( blue, green, - lambda b, g: m0 - m1 * (_safe_log(max(b, 1e-6)) / _safe_log(max(g, 1e-6))), + lambda b, g: m0 - m1 * (math.log(max(b, 1e-6)) / math.log(max(g, 1e-6))), ) updated = dict(scene.data) updated["stumpf_depth"] = stumpf @@ -117,6 +264,11 @@ def _maritorena_depth( attenuation: float = 0.1, **_: Any, ) -> Scene: + """Single-parameter exponential attenuation correction given a known `depth`. + + Scales each band by ``exp(attenuation * depth)``. + """ + if not isinstance(scene.data, dict): raise TypeError("Maritorena correction expects mapping-based scene data.") if depth is None: @@ -143,6 +295,9 @@ def _maritorena_depth( ) -register_depth_method("lyzenga", _lyzenga_depth) +register_depth_method("lyzenga1978", _lyzenga_1978) +register_depth_method("lyzenga1981", _lyzenga_1981) +register_depth_method("lyzenga2006", _lyzenga_2006) +register_depth_method("lyzenga", _lyzenga_2006) # alias: most recent/popular variant register_depth_method("stumpf", _stumpf_depth) register_depth_method("maritorena", _maritorena_depth) diff --git a/tests/test_processing.py b/tests/test_processing.py index 49b2485..d84550c 100644 --- a/tests/test_processing.py +++ b/tests/test_processing.py @@ -51,31 +51,98 @@ def test_hedley_glint_removal_requires_mapping_data() -> None: remove_glint(scene, method="hedley") -def test_lyzenga_depth_correction_default_variant() -> None: +def test_lyzenga1978_log_linearises_each_band() -> None: scene = Scene(data={"blue": [0.1], "qa60": [1]}, band_names=["blue", "qa60"]) - result = depth_correct(scene, method="lyzenga") + result = depth_correct(scene, method="lyzenga1978") - assert result.data["blue"] == pytest.approx([-math.log(0.1)]) + assert result.data["blue"] == pytest.approx([math.log(0.1)]) assert result.data["qa60"] == [1] - assert result.metadata["depth_method"] == "lyzenga_1981" - assert result.history == ("depth_correct:lyzenga:1981",) + assert result.metadata["depth_method"] == "lyzenga1978" + assert result.history == ("depth_correct:lyzenga1978",) -def test_lyzenga_depth_correction_records_requested_variant() -> None: - scene = Scene(data={"blue": [0.1]}, band_names=["blue"]) +def test_lyzenga1978_subtracts_deep_water_radiance() -> None: + scene = Scene(data={"blue": [0.5]}, band_names=["blue"]) - result = depth_correct(scene, method="lyzenga", variant="2006") + result = depth_correct(scene, method="lyzenga1978", deep_water_radiance={"blue": 0.2}) - assert result.metadata["depth_method"] == "lyzenga_2006" + assert result.data["blue"] == pytest.approx([math.log(0.3)]) -def test_lyzenga_depth_correction_divides_by_depth_when_given() -> None: - scene = Scene(data={"blue": [0.1]}, band_names=["blue"]) +def test_lyzenga1981_computes_depth_invariant_index() -> None: + scene = Scene( + data={"blue": [math.e], "green": [math.e**2]}, + band_names=["blue", "green"], + ) + + result = depth_correct(scene, method="lyzenga1981", k_i=3.0, k_j=4.0) + + assert result.data["lyzenga_dii"] == pytest.approx([-0.4]) + assert "lyzenga_dii" in result.band_names + assert result.metadata["depth_method"] == "lyzenga1981" + assert result.history == ("depth_correct:lyzenga1981",) + + +def test_lyzenga1981_requires_attenuation_coefficients() -> None: + scene = Scene(data={"blue": [0.1], "green": [0.2]}, band_names=["blue", "green"]) + + with pytest.raises(ValueError): + depth_correct(scene, method="lyzenga1981") + + +def test_lyzenga1981_rejects_both_attenuation_coefficients_zero() -> None: + scene = Scene(data={"blue": [0.1], "green": [0.2]}, band_names=["blue", "green"]) + + with pytest.raises(ValueError): + depth_correct(scene, method="lyzenga1981", k_i=0.0, k_j=0.0) + + +def test_lyzenga2006_computes_depth_regression() -> None: + scene = Scene( + data={"blue": [math.e], "green": [math.e**2]}, + band_names=["blue", "green"], + ) - result = depth_correct(scene, method="lyzenga", depth=[2.0]) + result = depth_correct( + scene, + method="lyzenga2006", + coefficients={"blue": 1.0, "green": 0.5}, + intercept=10.0, + ) + + assert result.data["lyzenga_depth"] == pytest.approx([8.0]) + assert "lyzenga_depth" in result.band_names + assert result.metadata["depth_method"] == "lyzenga2006" + assert result.history == ("depth_correct:lyzenga2006",) + + +def test_lyzenga2006_requires_coefficients() -> None: + scene = Scene(data={"blue": [0.1], "green": [0.2]}, band_names=["blue", "green"]) + + with pytest.raises(ValueError): + depth_correct(scene, method="lyzenga2006") + + +def test_lyzenga2006_requires_at_least_one_band() -> None: + scene = Scene(data={"blue": [0.1], "green": [0.2]}, band_names=["blue", "green"]) + + with pytest.raises(ValueError): + depth_correct(scene, method="lyzenga2006", bands=(), coefficients={}) + + +def test_bare_lyzenga_method_is_an_alias_for_lyzenga2006() -> None: + scene = Scene( + data={"blue": [math.e], "green": [math.e**2]}, + band_names=["blue", "green"], + ) + + result = depth_correct( + scene, method="lyzenga", coefficients={"blue": 1.0, "green": 0.5}, intercept=10.0 + ) - assert result.data["blue"] == pytest.approx([-math.log(0.1) / 2.0]) + assert result.data["lyzenga_depth"] == pytest.approx([8.0]) + assert result.metadata["depth_method"] == "lyzenga2006" def test_stumpf_depth_correction_adds_derived_band() -> None: