From 8b64c4d41a3bbce291cd29c15c850929d34b9e39 Mon Sep 17 00:00:00 2001 From: Ysobel Date: Tue, 28 Jul 2026 20:52:15 +1000 Subject: [PATCH] Fix Maritorena and Stumpf depth-correction formulas _maritorena_depth was previously a single-parameter exponential radiance scaler with no relation to any published Maritorena model. Replace it with the inverse of the Maritorena, Morel & Gentili (1994) two-flow shallow-water reflectance model: given known depth, deep-water reflectance (Rw_inf), and diffuse attenuation (K_d, repurposing the existing `attenuation` parameter), recover bottom reflectance Rb from measured water-leaving reflectance Rw via Rb = Rw_inf + (Rw - Rw_inf) / exp(-2 * K_d * depth), the algebraic inverse of the forward model Rw = Rw_inf + (Rb - Rw_inf) * exp(-2 * K_d * depth). Both Rw_inf and K_d now accept per-band mappings via the existing `_band_param` helper, matching the Lyzenga functions' convention. _stumpf_depth previously computed `m0 - m1 * ln(blue) / ln(green)` directly on raw reflectance, which both flips the sign convention of Stumpf, Holderied & Sinclair (2003) eq. (1) (`m1 * ratio - m0`, not `m0 - m1 * ratio`) and omits the paper's fixed scaling constant `n` (commonly 1000), which keeps `ln(n * Rw)` positive and away from the numerical instability of taking logs of small fractional reflectance values near zero. This is a genuine formula bug, not just a re-parameterisation: fix both the sign and add `n` (default 1000.0, matching the original paper) as a keyword. Behavioural change: existing callers of method="maritorena" or method="stumpf" with default parameters will see different numeric output, which is expected -- the old formulas did not implement the cited papers. Co-Authored-By: Claude Sonnet 5 --- src/earthrs/processing/depth.py | 65 +++++++++++++++++++--- tests/test_processing.py | 98 +++++++++++++++++++++++++++++++-- 2 files changed, 149 insertions(+), 14 deletions(-) diff --git a/src/earthrs/processing/depth.py b/src/earthrs/processing/depth.py index 8785057..d828ef0 100644 --- a/src/earthrs/processing/depth.py +++ b/src/earthrs/processing/depth.py @@ -217,11 +217,28 @@ 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 @@ -229,10 +246,11 @@ def _stumpf_depth( 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" diff --git a/tests/test_processing.py b/tests/test_processing.py index 11fb4c3..d8d5530 100644 --- a/tests/test_processing.py +++ b/tests/test_processing.py @@ -157,12 +157,51 @@ 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"]) @@ -170,13 +209,62 @@ def test_maritorena_depth_correction_requires_depth() -> None: 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: