From ccbe4a224ad54fef40d0bc4a6eb42769f535e648 Mon Sep 17 00:00:00 2001
From: Ysobel <ysobelsims@gmail.com>
Date: Sun, 12 Jul 2026 14:58:04 +1000
Subject: [PATCH 1/5] Add lyzenga variants

---
 src/earthrs/processing/core.py  |   5 -
 src/earthrs/processing/depth.py | 215 ++++++++++++++++++++++++++++----
 tests/test_processing.py        |  83 ++++++++++--
 3 files changed, 263 insertions(+), 40 deletions(-)

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..04ae156 100644
--- a/src/earthrs/processing/depth.py
+++ b/src/earthrs/processing/depth.py
@@ -1,19 +1,22 @@
-"""Depth-correction processing routines."""
+"""Depth-correction processing routines.
+"""
 
 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,
@@ -33,35 +36,42 @@ 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,
+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, ...], variant: str
 ) -> Scene:
-    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
+    """Build the result `Scene` shared by all Lyzenga variants."""
+
     metadata = dict(scene.metadata)
     metadata["depth_method"] = f"lyzenga_{variant}"
     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}"),
@@ -71,6 +81,151 @@ def _lyzenga_depth(
     )
 
 
+def _lyzenga_1978(
+    scene: Scene,
+    *,
+    deep_water_radiance: Mapping[str, float] | float | None = None,
+    epsilon: float = 1e-6,
+) -> 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.
+    """
+
+    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, variant="1978")
+
+
+def _lyzenga_1981(
+    scene: Scene,
+    *,
+    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,
+) -> 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.
+    """
+
+    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, variant="1981")
+
+
+def _lyzenga_2006(
+    scene: Scene,
+    *,
+    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,
+) -> 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.
+    """
+
+    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, variant="2006")
+
+
+_LYZENGA_VARIANTS = {
+    "1978": _lyzenga_1978,
+    "1981": _lyzenga_1981,
+    "2006": _lyzenga_2006,
+}
+
+
+def _lyzenga_depth(
+    scene: Scene,
+    *,
+    depth: Any | None = None,
+    variant: str = "1981",
+    **kwargs: Any,
+) -> Scene:
+    """Dispatch to the requested Lyzenga variant (``"1978"``, ``"1981"``, ``"2006"``).
+
+    Lyzenga methods derive depth information from radiance; unlike Maritorena's
+    attenuation model, none of them take an already-known depth as input, so
+    `depth` is accepted (for API symmetry with other depth-correction methods) but
+    ignored.
+    """
+
+    _ = depth
+    if not isinstance(scene.data, dict):
+        raise TypeError("Lyzenga correction expects mapping-based scene data.")
+    variant_key = str(variant)
+    implementation = _LYZENGA_VARIANTS.get(variant_key)
+    if implementation is None:
+        raise ValueError(
+            f"Unknown Lyzenga variant '{variant}'. Expected one of "
+            f"{sorted(_LYZENGA_VARIANTS)}."
+        )
+    return implementation(scene, **kwargs)
+
+
 def _stumpf_depth(
     scene: Scene,
     *,
@@ -81,6 +236,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 +249,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 +277,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:
diff --git a/tests/test_processing.py b/tests/test_processing.py
index 49b2485..d2b52d4 100644
--- a/tests/test_processing.py
+++ b/tests/test_processing.py
@@ -51,31 +51,94 @@ def test_hedley_glint_removal_requires_mapping_data() -> None:
         remove_glint(scene, method="hedley")
 
 
-def test_lyzenga_depth_correction_default_variant() -> None:
+def test_lyzenga_1978_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="lyzenga", variant="1978")
 
-    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_1978"
+    assert result.history == ("depth_correct:lyzenga:1978",)
+
+
+def test_lyzenga_1978_subtracts_deep_water_radiance() -> None:
+    scene = Scene(data={"blue": [0.5]}, band_names=["blue"])
+
+    result = depth_correct(
+        scene, method="lyzenga", variant="1978", deep_water_radiance={"blue": 0.2}
+    )
+
+    assert result.data["blue"] == pytest.approx([math.log(0.3)])
+
+
+def test_lyzenga_1981_computes_depth_invariant_index_by_default() -> None:
+    scene = Scene(
+        data={"blue": [math.e], "green": [math.e**2]},
+        band_names=["blue", "green"],
+    )
+
+    result = depth_correct(scene, method="lyzenga", 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"] == "lyzenga_1981"
     assert result.history == ("depth_correct:lyzenga:1981",)
 
 
-def test_lyzenga_depth_correction_records_requested_variant() -> None:
-    scene = Scene(data={"blue": [0.1]}, band_names=["blue"])
+def test_lyzenga_1981_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="lyzenga")
 
-    result = depth_correct(scene, method="lyzenga", variant="2006")
 
+def test_lyzenga_1981_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="lyzenga", k_i=0.0, k_j=0.0)
+
+
+def test_lyzenga_2006_computes_depth_regression() -> None:
+    scene = Scene(
+        data={"blue": [math.e], "green": [math.e**2]},
+        band_names=["blue", "green"],
+    )
+
+    result = depth_correct(
+        scene,
+        method="lyzenga",
+        variant="2006",
+        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"] == "lyzenga_2006"
+    assert result.history == ("depth_correct:lyzenga:2006",)
 
 
-def test_lyzenga_depth_correction_divides_by_depth_when_given() -> None:
-    scene = Scene(data={"blue": [0.1]}, band_names=["blue"])
+def test_lyzenga_2006_requires_coefficients() -> None:
+    scene = Scene(data={"blue": [0.1], "green": [0.2]}, band_names=["blue", "green"])
+
+    with pytest.raises(ValueError):
+        depth_correct(scene, method="lyzenga", variant="2006")
 
-    result = depth_correct(scene, method="lyzenga", depth=[2.0])
 
-    assert result.data["blue"] == pytest.approx([-math.log(0.1) / 2.0])
+def test_lyzenga_2006_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="lyzenga", variant="2006", bands=(), coefficients={})
+
+
+def test_lyzenga_rejects_unknown_variant() -> None:
+    scene = Scene(data={"blue": [0.1]}, band_names=["blue"])
+
+    with pytest.raises(ValueError):
+        depth_correct(scene, method="lyzenga", variant="1999")
 
 
 def test_stumpf_depth_correction_adds_derived_band() -> None:

From cbe3f9d412ffc1a162bd98e2981be106a5d491f4 Mon Sep 17 00:00:00 2001
From: Ysobel <ysobelsims@gmail.com>
Date: Mon, 13 Jul 2026 08:07:13 +1000
Subject: [PATCH 2/5] Update docs

---
 docs/processing.md | 100 +++++++++++++++++++++++++++++++++++++++++++--
 docs/status.md     |   2 -
 2 files changed, 96 insertions(+), 6 deletions(-)

diff --git a/docs/processing.md b/docs/processing.md
index 7b32a5d..2123201 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**: `"lyzenga"` (three distinct variants, see below), `"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,97 @@ 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
+
+`depth_correct(scene, method="lyzenga", variant=...)` selects between three literature-accurate
+formulas — they are genuinely different algorithms, not the same transform under three names, so
+each has its own required keyword arguments.
+
+### `"1978"` — 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="lyzenga", variant="1978", 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)`.
+
+### `"1981"` — 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="lyzenga",
+    variant="1981",
+    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 `"1978"`.
+
+This is the default variant (`depth_correct(scene, method="lyzenga")` is equivalent to
+`variant="1981"`), so `k_i`/`k_j` must be supplied even when `variant` is omitted.
+
+### `"2006"` — 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="lyzenga",
+    variant="2006",
+    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 `"1978"`.
+
+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
 

From a4d8088a08e6a4416ba9f2403ac953b76345af28 Mon Sep 17 00:00:00 2001
From: Ysobel <ysobelsims@gmail.com>
Date: Mon, 13 Jul 2026 14:13:22 +1000
Subject: [PATCH 3/5] dont use separate variant arg

---
 docs/processing.md              | 31 ++++++-------
 src/earthrs/processing/depth.py | 80 +++++++++++++++------------------
 tests/test_processing.py        | 62 +++++++++++++------------
 3 files changed, 83 insertions(+), 90 deletions(-)

diff --git a/docs/processing.md b/docs/processing.md
index 2123201..934783d 100644
--- a/docs/processing.md
+++ b/docs/processing.md
@@ -19,7 +19,8 @@ register_glint_method("my_method", my_glint_method)
 ## Implemented methods
 
 - **Glint removal**: `"hedley"` (Hedley et al., 2005).
-- **Depth correction**: `"lyzenga"` (three distinct variants, see below), `"stumpf"` (ratio
+- **Depth correction**: `"lyzenga1978"`, `"lyzenga1981"`, `"lyzenga2006"` (three distinct
+  formulas, see below; plain `"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).
@@ -32,11 +33,12 @@ yet have its own registry.
 
 ## Lyzenga depth-correction variants
 
-`depth_correct(scene, method="lyzenga", variant=...)` selects between three literature-accurate
-formulas — they are genuinely different algorithms, not the same transform under three names, so
-each has its own required keyword arguments.
+The Lyzenga formula is registered once per publication year — `"lyzenga1978"`, `"lyzenga1981"`,
+`"lyzenga2006"` — since each is a genuinely different algorithm, not the same transform under
+three names, so each has its own required keyword arguments. Plain `method="lyzenga"` is an
+alias for `"lyzenga2006"`, the most recent and most widely used variant.
 
-### `"1978"` — per-band log-linearising transform
+### `"lyzenga1978"` — per-band log-linearising transform
 
 The single-band transform from Lyzenga (1978), eqs. (1) & (7). Replaces every spectral band with
 
@@ -49,7 +51,7 @@ through unchanged. Bands are not combined with each other, so this variant does
 bottom-reflectance variation — it's the foundational linearisation the other two variants build on.
 
 ```python
-depth_correct(scene, method="lyzenga", variant="1978", deep_water_radiance={"blue": 0.02})
+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
@@ -57,7 +59,7 @@ depth_correct(scene, method="lyzenga", variant="1978", deep_water_radiance={"blu
 - `epsilon` (default `1e-6`) — floor applied to `L_i - L_si` before the log, to avoid `log(0)`
   or `log(negative)`.
 
-### `"1981"` — two-band depth-invariant bottom index
+### `"lyzenga1981"` — two-band depth-invariant bottom index
 
 The classic two-band index from Lyzenga (1981), eq. (2). Adds a `lyzenga_dii` band:
 
@@ -73,8 +75,7 @@ water depth.
 ```python
 depth_correct(
     scene,
-    method="lyzenga",
-    variant="1981",
+    method="lyzenga1981",
     band_i="blue",
     band_j="green",
     k_i=0.056,
@@ -85,12 +86,9 @@ depth_correct(
 - `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 `"1978"`.
+- `deep_water_radiance`, `epsilon` — as in `"lyzenga1978"`.
 
-This is the default variant (`depth_correct(scene, method="lyzenga")` is equivalent to
-`variant="1981"`), so `k_i`/`k_j` must be supplied even when `variant` is omitted.
-
-### `"2006"` — empirical multi-band depth regression
+### `"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):
@@ -107,8 +105,7 @@ depth measurements (for example LIDAR) for the site being processed.
 ```python
 depth_correct(
     scene,
-    method="lyzenga",
-    variant="2006",
+    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,
@@ -119,7 +116,7 @@ depth_correct(
 - `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 `"1978"`.
+- `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/src/earthrs/processing/depth.py b/src/earthrs/processing/depth.py
index 04ae156..02d3740 100644
--- a/src/earthrs/processing/depth.py
+++ b/src/earthrs/processing/depth.py
@@ -27,7 +27,11 @@ 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. The Lyzenga
+    formula is registered under a distinct name per publication year --
+    ``"lyzenga1978"``, ``"lyzenga1981"``, ``"lyzenga2006"`` -- since each is a
+    genuinely different algorithm; plain ``"lyzenga"`` is an alias for
+    ``"lyzenga2006"``, the most recent and most widely used variant.
     """
 
     processor = _DEPTH_REGISTRY.get(method.lower())
@@ -36,6 +40,11 @@ def depth_correct(
     return processor(scene, depth=depth, **kwargs)
 
 
+def _require_mapping_data(scene: Scene, label: str) -> None:
+    if not isinstance(scene.data, dict):
+        raise TypeError(f"{label} correction expects mapping-based scene data.")
+
+
 def _band_param(
     bands: Iterable[str], value: Any, default: float
 ) -> dict[str, float]:
@@ -60,12 +69,12 @@ def _log_transform(values: Any, deep_water_radiance: float, epsilon: float) -> A
 
 
 def _finalise_lyzenga(
-    scene: Scene, *, data: dict[str, Any], band_names: tuple[str, ...], variant: str
+    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=data,
         crs=scene.crs,
@@ -74,7 +83,7 @@ def _finalise_lyzenga(
         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,
@@ -84,8 +93,10 @@ def _finalise_lyzenga(
 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.
 
@@ -94,6 +105,8 @@ def _lyzenga_1978(
     ``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():
@@ -101,18 +114,22 @@ def _lyzenga_1978(
             updated[band] = values
             continue
         updated[band] = _log_transform(values, dwr[band], epsilon)
-    return _finalise_lyzenga(scene, data=updated, band_names=scene.band_names, variant="1978")
+    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.
 
@@ -123,6 +140,8 @@ def _lyzenga_1981(
     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 "
@@ -141,17 +160,19 @@ def _lyzenga_1981(
     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, variant="1981")
+    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.
 
@@ -159,9 +180,12 @@ def _lyzenga_2006(
     ``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.
+    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:
@@ -188,42 +212,7 @@ def _lyzenga_2006(
     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, variant="2006")
-
-
-_LYZENGA_VARIANTS = {
-    "1978": _lyzenga_1978,
-    "1981": _lyzenga_1981,
-    "2006": _lyzenga_2006,
-}
-
-
-def _lyzenga_depth(
-    scene: Scene,
-    *,
-    depth: Any | None = None,
-    variant: str = "1981",
-    **kwargs: Any,
-) -> Scene:
-    """Dispatch to the requested Lyzenga variant (``"1978"``, ``"1981"``, ``"2006"``).
-
-    Lyzenga methods derive depth information from radiance; unlike Maritorena's
-    attenuation model, none of them take an already-known depth as input, so
-    `depth` is accepted (for API symmetry with other depth-correction methods) but
-    ignored.
-    """
-
-    _ = depth
-    if not isinstance(scene.data, dict):
-        raise TypeError("Lyzenga correction expects mapping-based scene data.")
-    variant_key = str(variant)
-    implementation = _LYZENGA_VARIANTS.get(variant_key)
-    if implementation is None:
-        raise ValueError(
-            f"Unknown Lyzenga variant '{variant}'. Expected one of "
-            f"{sorted(_LYZENGA_VARIANTS)}."
-        )
-    return implementation(scene, **kwargs)
+    return _finalise_lyzenga(scene, data=updated, band_names=band_names, method="lyzenga2006")
 
 
 def _stumpf_depth(
@@ -308,6 +297,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 d2b52d4..d84550c 100644
--- a/tests/test_processing.py
+++ b/tests/test_processing.py
@@ -51,56 +51,54 @@ def test_hedley_glint_removal_requires_mapping_data() -> None:
         remove_glint(scene, method="hedley")
 
 
-def test_lyzenga_1978_log_linearises_each_band() -> 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", variant="1978")
+    result = depth_correct(scene, method="lyzenga1978")
 
     assert result.data["blue"] == pytest.approx([math.log(0.1)])
     assert result.data["qa60"] == [1]
-    assert result.metadata["depth_method"] == "lyzenga_1978"
-    assert result.history == ("depth_correct:lyzenga:1978",)
+    assert result.metadata["depth_method"] == "lyzenga1978"
+    assert result.history == ("depth_correct:lyzenga1978",)
 
 
-def test_lyzenga_1978_subtracts_deep_water_radiance() -> None:
+def test_lyzenga1978_subtracts_deep_water_radiance() -> None:
     scene = Scene(data={"blue": [0.5]}, band_names=["blue"])
 
-    result = depth_correct(
-        scene, method="lyzenga", variant="1978", deep_water_radiance={"blue": 0.2}
-    )
+    result = depth_correct(scene, method="lyzenga1978", deep_water_radiance={"blue": 0.2})
 
     assert result.data["blue"] == pytest.approx([math.log(0.3)])
 
 
-def test_lyzenga_1981_computes_depth_invariant_index_by_default() -> None:
+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="lyzenga", k_i=3.0, k_j=4.0)
+    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"] == "lyzenga_1981"
-    assert result.history == ("depth_correct:lyzenga:1981",)
+    assert result.metadata["depth_method"] == "lyzenga1981"
+    assert result.history == ("depth_correct:lyzenga1981",)
 
 
-def test_lyzenga_1981_requires_attenuation_coefficients() -> None:
+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="lyzenga")
+        depth_correct(scene, method="lyzenga1981")
 
 
-def test_lyzenga_1981_rejects_both_attenuation_coefficients_zero() -> None:
+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="lyzenga", k_i=0.0, k_j=0.0)
+        depth_correct(scene, method="lyzenga1981", k_i=0.0, k_j=0.0)
 
 
-def test_lyzenga_2006_computes_depth_regression() -> None:
+def test_lyzenga2006_computes_depth_regression() -> None:
     scene = Scene(
         data={"blue": [math.e], "green": [math.e**2]},
         band_names=["blue", "green"],
@@ -108,37 +106,43 @@ def test_lyzenga_2006_computes_depth_regression() -> None:
 
     result = depth_correct(
         scene,
-        method="lyzenga",
-        variant="2006",
+        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"] == "lyzenga_2006"
-    assert result.history == ("depth_correct:lyzenga:2006",)
+    assert result.metadata["depth_method"] == "lyzenga2006"
+    assert result.history == ("depth_correct:lyzenga2006",)
 
 
-def test_lyzenga_2006_requires_coefficients() -> None:
+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="lyzenga", variant="2006")
+        depth_correct(scene, method="lyzenga2006")
 
 
-def test_lyzenga_2006_requires_at_least_one_band() -> None:
+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="lyzenga", variant="2006", bands=(), coefficients={})
+        depth_correct(scene, method="lyzenga2006", bands=(), coefficients={})
 
 
-def test_lyzenga_rejects_unknown_variant() -> None:
-    scene = Scene(data={"blue": [0.1]}, band_names=["blue"])
+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"],
+    )
 
-    with pytest.raises(ValueError):
-        depth_correct(scene, method="lyzenga", variant="1999")
+    result = depth_correct(
+        scene, method="lyzenga", coefficients={"blue": 1.0, "green": 0.5}, intercept=10.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:

From 8faf8faad74a3e67236aae98854dbe2a77f29daf Mon Sep 17 00:00:00 2001
From: Ysobel Sims <35280100+ysims@users.noreply.github.com>
Date: Mon, 13 Jul 2026 18:54:45 +1000
Subject: [PATCH 4/5] Apply suggestions from code review

Co-authored-by: Ysobel Sims <35280100+ysims@users.noreply.github.com>
---
 docs/processing.md              | 8 +++-----
 src/earthrs/processing/depth.py | 9 +++------
 2 files changed, 6 insertions(+), 11 deletions(-)

diff --git a/docs/processing.md b/docs/processing.md
index 934783d..9fa285a 100644
--- a/docs/processing.md
+++ b/docs/processing.md
@@ -19,8 +19,7 @@ register_glint_method("my_method", my_glint_method)
 ## Implemented methods
 
 - **Glint removal**: `"hedley"` (Hedley et al., 2005).
-- **Depth correction**: `"lyzenga1978"`, `"lyzenga1981"`, `"lyzenga2006"` (three distinct
-  formulas, see below; plain `"lyzenga"` is an alias for `"lyzenga2006"`), `"stumpf"` (ratio
+- **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).
@@ -33,9 +32,8 @@ yet have its own registry.
 
 ## Lyzenga depth-correction variants
 
-The Lyzenga formula is registered once per publication year — `"lyzenga1978"`, `"lyzenga1981"`,
-`"lyzenga2006"` — since each is a genuinely different algorithm, not the same transform under
-three names, so each has its own required keyword arguments. Plain `method="lyzenga"` is an
+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
diff --git a/src/earthrs/processing/depth.py b/src/earthrs/processing/depth.py
index 02d3740..6d7c9b7 100644
--- a/src/earthrs/processing/depth.py
+++ b/src/earthrs/processing/depth.py
@@ -1,5 +1,4 @@
-"""Depth-correction processing routines.
-"""
+"""Depth-correction processing routines."""
 
 from __future__ import annotations
 
@@ -27,10 +26,8 @@ def depth_correct(
 ) -> Scene:
     """Apply depth correction to raster data.
 
-    Method names are resolved from the depth-correction registry. The Lyzenga
-    formula is registered under a distinct name per publication year --
-    ``"lyzenga1978"``, ``"lyzenga1981"``, ``"lyzenga2006"`` -- since each is a
-    genuinely different algorithm; plain ``"lyzenga"`` is an alias for
+    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.
     """
 

From 0d7ea8159f764d9b3113386a27c06cba7d3a2a55 Mon Sep 17 00:00:00 2001
From: Ysobel <ysobelsims@gmail.com>
Date: Mon, 13 Jul 2026 18:56:38 +1000
Subject: [PATCH 5/5] lint

---
 src/earthrs/processing/depth.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/earthrs/processing/depth.py b/src/earthrs/processing/depth.py
index 6d7c9b7..3f0e4d2 100644
--- a/src/earthrs/processing/depth.py
+++ b/src/earthrs/processing/depth.py
@@ -26,9 +26,10 @@ def depth_correct(
 ) -> Scene:
     """Apply depth correction to raster data.
 
-    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.
+    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())