Fix Maritorena and Stumpf depth-correction formulas

src/earthrs/processing/depth.py

@@ -217,22 +217,40 @@ def _stumpf_depth(
     green_band: str = "green",
     m0: float = 0.0,
     m1: float = 1.0,
+    n: float = 1000.0,
+    epsilon: float = 1e-6,
     **_: Any,
 ) -> Scene:
-    """Stumpf (2003) band-ratio depth transform.
-
-    Adds a ``stumpf_depth`` band from ``m0 - m1 * ln(blue) / ln(green)``.
+    """Stumpf, Holderied & Sinclair (2003) band-ratio depth transform.
+
+    Adds a ``stumpf_depth`` band from eq. (1),
+    ``depth = m1 * ln(n * Rw_i) / ln(n * Rw_j) - m0``, where `Rw_i`/`Rw_j` are
+    `blue_band`/`green_band` reflectance (blue is the numerator band, ``i``; green
+    is the denominator band, ``j``, following the paper's ``lambda_i``/``lambda_j``
+    notation), and `m0`/`m1` are empirically tuned per-image
+    coefficients. `n` is a fixed scaling constant, chosen so that ``n * Rw`` stays
+    comfortably above 1 across the sensor's expected reflectance range: this keeps
+    ``ln(n * Rw)`` positive (avoiding a sign flip in the ratio) and well away from
+    zero (avoiding the numerical instability that taking logs of small fractional
+    reflectance values, close to zero, is prone to). The original paper uses
+    ``n = 1000``. `epsilon` floors ``n * Rw`` at ``1 + epsilon`` so the log stays
+    strictly positive even for zero or negative input reflectance.
+
+    Citation: Stumpf, R. P., Holderied, K., & Sinclair, M. (2003). "Determination
+    of water depth with high-resolution satellite imagery over variable bottom
+    types." Limnology and Oceanography, 48(1), 547-556, eq. (1).
     """
 
     _ = depth
     if not isinstance(scene.data, dict):
         raise TypeError("Stumpf correction expects mapping-based scene data.")
     blue = scene.data[blue_band]
     green = scene.data[green_band]
+    floor = 1.0 + epsilon
     stumpf = _map_binary(
         blue,
         green,
-        lambda b, g: m0 - m1 * (math.log(max(b, 1e-6)) / math.log(max(g, 1e-6))),
+        lambda b, g: m1 * (math.log(max(n * b, floor)) / math.log(max(n * g, floor))) - m0,
     )
     updated = dict(scene.data)
     updated["stumpf_depth"] = stumpf
@@ -257,22 +275,51 @@ def _maritorena_depth(
     scene: Scene,
     *,
     depth: Any | None = None,
-    attenuation: float = 0.1,
+    deep_water_reflectance: Mapping[str, float] | float | None = 0.0,
+    attenuation: Mapping[str, float] | float | None = 0.1,
     **_: Any,
 ) -> Scene:
-    """Single-parameter exponential attenuation correction given a known `depth`.
-
-    Scales each band by ``exp(attenuation * depth)``.
+    """Maritorena, Morel & Gentili (1994) two-flow shallow-water reflectance model.
+
+    The forward model for water-leaving reflectance ``Rw`` over a finite-depth
+    bottom is ``Rw(depth) = Rw_inf + (Rb - Rw_inf) * exp(-2 * K_d * depth)``, where
+    ``Rw_inf`` is the reflectance of optically-deep water (same band, no bottom
+    contribution), ``Rb`` is the bottom albedo/reflectance, ``K_d`` is the diffuse
+    attenuation coefficient, and the factor of 2 accounts for the two-way (down-
+    and up-welling) light path. Unlike the Lyzenga/Stumpf transforms in this
+    module, `depth` here is a *known* input rather than a quantity being
+    estimated, so this correction inverts the forward model to recover bottom
+    reflectance from measured water-leaving reflectance:
+    ``Rb = Rw_inf + (Rw_measured - Rw_inf) / exp(-2 * K_d * depth)``.
+
+    `deep_water_reflectance` supplies ``Rw_inf`` (per-band mapping or scalar,
+    default ``0.0``). `attenuation` supplies ``K_d`` (per-band mapping or scalar,
+    default ``0.1``), in units matched to `depth`'s units.
+
+    This model is also the basis of later semi-analytical shallow-water
+    inversions (e.g. Lee et al.).
+
+    Citation: Maritorena, S., Morel, A., & Gentili, B. (1994). "Diffuse
+    reflectance of oceanic shallow waters: influence of water depth and bottom
+    albedo." Limnology and Oceanography, 39(7), 1689-1703.
     """
 
     if not isinstance(scene.data, dict):
         raise TypeError("Maritorena correction expects mapping-based scene data.")
     if depth is None:
         raise ValueError("Maritorena correction requires `depth`.")
+    rw_inf = _band_param(scene.data.keys(), deep_water_reflectance, 0.0)
+    k_d = _band_param(scene.data.keys(), attenuation, 0.1)
     updated = {}
     for band, values in scene.data.items():
+        band_rw_inf = rw_inf[band]
+        band_k_d = k_d[band]
         updated[band] = _map_binary(
-            values, depth, lambda value, d: value * math.exp(attenuation * d)
+            values,
+            depth,
+            lambda rw, d, rw_inf=band_rw_inf, kd=band_k_d: (
+                rw_inf + (rw - rw_inf) / math.exp(-2.0 * kd * d)
+            ),
         )
     metadata = dict(scene.metadata)
     metadata["depth_method"] = "maritorena"

tests/test_processing.py

@@ -187,26 +187,114 @@ def test_stumpf_depth_correction_adds_derived_band() -> None:
 
     result = depth_correct(scene, method="stumpf")
 
-    expected = 0.0 - 1.0 * (math.log(0.1) / math.log(0.2))
+    n = 1000.0
+    expected = 1.0 * (math.log(n * 0.1) / math.log(n * 0.2)) - 0.0
     assert result.data["stumpf_depth"] == pytest.approx([expected])
     assert "stumpf_depth" in result.band_names
     assert result.metadata["depth_method"] == "stumpf"
 
 
+def test_stumpf_depth_matches_hand_computed_values() -> None:
+    # Choose blue/green so that n * Rw is exactly e and e**2 respectively, giving
+    # ln(n * blue) == 1 and ln(n * green) == 2, i.e. a ratio of exactly 0.5.
+    n = 1000.0
+    scene = Scene(
+        data={"blue": [math.e / n], "green": [math.e**2 / n]},
+        band_names=["blue", "green"],
+    )
+
+    default_result = depth_correct(scene, method="stumpf")
+    assert default_result.data["stumpf_depth"] == pytest.approx([0.5])
+
+    tuned_result = depth_correct(scene, method="stumpf", m0=1.0, m1=2.0)
+    assert tuned_result.data["stumpf_depth"] == pytest.approx([0.0])
+
+
+def test_stumpf_depth_ratio_increases_with_relative_blue_reflectance() -> None:
+    # Deeper water attenuates the longer green wavelength more than blue, so the
+    # blue/green ratio -- and hence the Stumpf depth estimate -- should increase
+    # as green reflectance drops relative to blue (eq. 1's sign convention).
+    shallow = Scene(data={"blue": [0.05], "green": [0.05]}, band_names=["blue", "green"])
+    deep = Scene(data={"blue": [0.05], "green": [0.01]}, band_names=["blue", "green"])
+
+    shallow_depth = depth_correct(shallow, method="stumpf").data["stumpf_depth"][0]
+    deep_depth = depth_correct(deep, method="stumpf").data["stumpf_depth"][0]
+
+    assert deep_depth > shallow_depth
+
+
+def test_stumpf_depth_accepts_custom_scaling_constant() -> None:
+    scene = Scene(data={"blue": [0.1], "green": [0.2]}, band_names=["blue", "green"])
+
+    result = depth_correct(scene, method="stumpf", n=500.0)
+
+    expected = 1.0 * (math.log(500.0 * 0.1) / math.log(500.0 * 0.2)) - 0.0
+    assert result.data["stumpf_depth"] == pytest.approx([expected])
+
+
 def test_maritorena_depth_correction_requires_depth() -> None:
     scene = Scene(data={"blue": [1.0]}, band_names=["blue"])
 
     with pytest.raises(ValueError):
         depth_correct(scene, method="maritorena")
 
 
-def test_maritorena_depth_correction_applies_exponential_attenuation() -> None:
-    scene = Scene(data={"blue": [1.0]}, band_names=["blue"])
+def test_maritorena_depth_correction_recovers_bottom_reflectance_round_trip() -> None:
+    # Run the forward Maritorena, Morel & Gentili (1994) model by hand to derive a
+    # synthetic measured Rw at a known depth, then check the correction inverts it
+    # back to the original bottom reflectance Rb.
+    rb = 0.18
+    rw_inf = 0.02
+    k_d = 0.15
+    depth = 3.5
+    measured_rw = rw_inf + (rb - rw_inf) * math.exp(-2.0 * k_d * depth)
 
-    result = depth_correct(scene, method="maritorena", depth=[2.0], attenuation=0.1)
+    scene = Scene(data={"blue": [measured_rw]}, band_names=["blue"])
+    result = depth_correct(
+        scene,
+        method="maritorena",
+        depth=[depth],
+        deep_water_reflectance=rw_inf,
+        attenuation=k_d,
+    )
 
-    assert result.data["blue"] == pytest.approx([math.exp(0.2)])
+    assert result.data["blue"] == pytest.approx([rb])
     assert result.metadata["depth_method"] == "maritorena"
+    assert result.history == ("depth_correct:maritorena",)
+
+
+def test_maritorena_depth_correction_zero_depth_recovers_reflectance_exactly() -> None:
+    scene = Scene(data={"blue": [0.42]}, band_names=["blue"])
+
+    result = depth_correct(
+        scene, method="maritorena", depth=[0.0], deep_water_reflectance=0.05, attenuation=0.2
+    )
+
+    assert result.data["blue"] == pytest.approx([0.42])
+
+
+def test_maritorena_depth_correction_supports_per_band_parameters() -> None:
+    rb_blue, rb_green = 0.2, 0.1
+    rw_inf = {"blue": 0.03, "green": 0.05}
+    k_d = {"blue": 0.1, "green": 0.2}
+    depth = 4.0
+    measured_blue = rw_inf["blue"] + (rb_blue - rw_inf["blue"]) * math.exp(
+        -2.0 * k_d["blue"] * depth
+    )
+    measured_green = rw_inf["green"] + (rb_green - rw_inf["green"]) * math.exp(
+        -2.0 * k_d["green"] * depth
+    )
+
+    scene = Scene(
+        data={"blue": [measured_blue], "green": [measured_green]},
+        band_names=["blue", "green"],
+    )
+    result = depth_correct(
+        scene, method="maritorena", depth=[depth], deep_water_reflectance=rw_inf, attenuation=k_d
+    )
+
+    assert result.data["blue"] == pytest.approx([rb_blue])
+    assert result.data["green"] == pytest.approx([rb_green])
 
 
 def test_sentinel2_qa60_cloud_mask_checks_bits_10_and_11() -> None: