Feat/calculate nrd

configs/config.yaml

@@ -36,13 +36,15 @@ preprocessing:
     nodata: 0
 
 pad_profile:
-  scanning_angle: True
-  limit_N_points: 400
-  limit_flight_agl: 800
-  z0: 0.0
-  dz: 1.0
-  nlayers: 60
-  ground_margin: 0.1
+  cos_theta:
+    scanning_angle: True
+    limit_N_points: 400
+    limit_flight_agl: 800
+  compute_ni_n:
+    z0: 0.0
+    dz: 1.0
+    nlayers: 60
+    ground_margin: 0.1
 
 dtm:
   download:

lidar_for_fuel/main_pad_profile.py

@@ -131,13 +131,13 @@ def main_on_one_tile(filename):
             input_filename=os.path.join(input_dir, filename),
             srid=config.io.spatial_reference,
             deviation_days=config.commons.filter_date.deviation_days,
-            scanning_angle=config.pad_profile.scanning_angle,
-            limit_N_points=config.pad_profile.limit_N_points,
-            limit_flight_agl=config.pad_profile.limit_flight_agl,
-            z0=config.pad_profile.z0,
-            dz=config.pad_profile.dz,
-            nlayers=config.pad_profile.nlayers,
-            ground_margin=config.pad_profile.ground_margin,
+            scanning_angle=config.pad_profile.cos_theta.scanning_angle,
+            limit_N_points=config.pad_profile.cos_theta.limit_N_points,
+            limit_flight_agl=config.pad_profile.cos_theta.limit_flight_agl,
+            z0=config.pad_profile.compute_ni_n.z0,
+            dz=config.pad_profile.compute_ni_n.dz,
+            nlayers=config.pad_profile.compute_ni_n.nlayers,
+            ground_margin=config.pad_profile.compute_ni_n.ground_margin,
         )
 
     if initial_las_filename:

lidar_for_fuel/pad_profile/calculate_pad_profile.py

@@ -10,6 +10,7 @@
 from lidar_for_fuel.commons.filter_points_by_date import filter_by_date
 from lidar_for_fuel.pad_profile.compute_cos_theta import compute_cos_theta
 from lidar_for_fuel.pad_profile.compute_ni_n import compute_ni_n
+from lidar_for_fuel.pad_profile.compute_nrd import compute_nrd
 
 logger = logging.getLogger(__name__)
 
@@ -69,15 +70,16 @@ def pad_metrics_core(
             (m). Default 0.1.
 
     Returns:
-        tuple[float, np.ndarray, np.ndarray, np.ndarray] | None: `(cos_theta, ni, n,
-        min_layer)`, or None if a quality guard fails.
+        tuple[float, np.ndarray, np.ndarray, np.ndarray, np.ndarray] | None: `(cos_theta, ni, n,
+        min_layer, nrd)`, or None if a quality guard fails.
             cos_theta: scan angle factor (1.0 if `scanning_angle=False`).
             ni: vegetation/ground hit count per stratum.
             n: cumulative entering-ray count per stratum (non-decreasing). A
                 point below the ground margin still counts toward every
                 stratum's `n` via the cumulative sum, since a ray ending there
                 had to pass through every stratum above it on the way down.
             min_layer: lower height bound of each stratum (m), pre-ground-margin-shift.
+            ndr: fractions of incoming rays intercepted for each "NRD" stratum..
     """
     # # Step 1:
     # Filter points by ±deviation_days around the most densely sampled calendar day.
@@ -130,4 +132,8 @@ def pad_metrics_core(
     # Then, get number of returns intercepted and "pulses" entering in each strata
     Ni, N, min_layer = compute_ni_n(h_abg, veg_ground_points, z0, dz, nlayers, ground_margin)
 
-    return cos_theta, Ni, N, min_layer
+    # # Step 6:
+    # Compute the fractions of incoming rays intercepted for each "NRD" stratum.
+    NRD = compute_nrd(Ni, N)
+
+    return cos_theta, Ni, N, min_layer, NRD

lidar_for_fuel/pad_profile/compute_nrd.py

@@ -0,0 +1,39 @@
+"""
+Calculate the fractions of incoming rays intercepted for each "NRD" stratum.
+
+"""
+import numpy as np
+
+
+def compute_nrd(
+    Ni: np.ndarray,
+    N: np.ndarray,
+) -> np.ndarray:
+    """Calculate the fractions of incoming rays intercepted for each "NRD" stratum.
+
+    Args:
+        Ni (np.ndarray): each of length `nlayers`, vegetation/ground hit count per stratum.
+        N (np.ndarray): each of length `nlayers`, cumulative entering-ray count per stratum (non-decreasing).
+
+    Returns:
+        np.ndarray: NRD values, corrected for boundary cases (NRD=0 or NRD=1)
+            per equations 23-24 of Pimont et al. 2018.
+    """
+    # General case:
+    # Compute the fraction of intercepted rays per stratum.
+    # /!\ np.errstate(invalid="ignore") suppresses the RuntimeWarning triggered by 0/0 (when N=0 and Ni=0).
+    with np.errstate(invalid="ignore"):
+        nrd = Ni / N
+
+    # When N=0 then Ni=0 (no incoming rays → no interception):
+    # 0/0 yields NaN, replaced by 0 (null intercepted fraction).
+    nrd[np.isnan(nrd)] = 0.0
+
+    # Boundary cases NRD=0 or NRD=1:
+    # Laplace Bayesian correction (Pimont et al. 2018, eq. 23-24).
+    # These extreme values make the PAD computation unstable.
+    # => (Ni+1)/(N+2) pulls the extremes slightly toward the center without affecting intermediate values.
+    i_nrdc = (nrd == 0) | (nrd == 1)
+    nrd[i_nrdc] = (Ni[i_nrdc] + 1) / (N[i_nrdc] + 2)
+
+    return nrd

test/pad_profile/test_calculate_pad_profile.py

@@ -47,7 +47,7 @@ def test_pad_metrics_core_returns_cos_theta_between_0_and_1():
 
 def test_pad_metrics_core_output_format():
     """Verify the output format: `None` when there are too few points, otherwise a
-    `(cos_theta, ni, n, min_layer)` tuple shaped for the default 60 strata."""
+    `(cos_theta, ni, n, min_layer, nrd)` tuple shaped for the default 60 strata."""
     n = 5
     gpstime = np.zeros(n, dtype=np.float64)
     points = _points(n, gpstime)
@@ -62,17 +62,17 @@ def test_pad_metrics_core_output_format():
     )
     assert result_none is None
 
-    # Enough points -> (cos_theta, ni, n, min_layer)
+    # Enough points -> (cos_theta, ni, n, min_layer, nrd)
     result = pad_metrics_core(
         **points,
         **_DEFAULT_PARAMS,
         scanning_angle=False,
         limit_N_points=1,
         deviation_days=np.inf,
     )
-    cos_theta, ni, n_per_stratum, min_layer = result
+    cos_theta, ni, n_per_stratum, min_layer, nrd = result
     assert isinstance(cos_theta, float)
-    assert len(ni) == len(n_per_stratum) == len(min_layer) == 60
+    assert len(ni) == len(n_per_stratum) == len(min_layer) == len(nrd) == 60
     assert np.issubdtype(ni.dtype, np.integer)
     assert np.issubdtype(n_per_stratum.dtype, np.integer)
     assert min_layer.dtype == np.float64
@@ -81,6 +81,8 @@ def test_pad_metrics_core_output_format():
     # full point total via the cumulative sum (see compute_ni_n).
     np.testing.assert_array_equal(ni, np.zeros(60, dtype=ni.dtype))
     np.testing.assert_array_equal(n_per_stratum, np.full(60, n, dtype=n_per_stratum.dtype))
+    # ni=0, n=5 for every stratum -> NRD=0 -> Laplace correction: (0+1)/(5+2) = 1/7
+    np.testing.assert_allclose(nrd, np.full(60, 1 / 7))
 
 
 def test_pad_metrics_core_scanning_angle_false_returns_one():
@@ -96,5 +98,5 @@ def test_pad_metrics_core_scanning_angle_false_returns_one():
         deviation_days=np.inf,
     )
 
-    cos_theta, _, _, _ = result
+    cos_theta, _, _, _, _ = result
     assert cos_theta == 1.0

test/pad_profile/test_compute_nrd.py

@@ -0,0 +1,72 @@
+import numpy as np
+import pytest
+
+from lidar_for_fuel.pad_profile.compute_nrd import compute_nrd
+
+# ── NRD values for various Ni/N combinations ──────────────────────────────────
+
+
+@pytest.mark.parametrize(
+    "Ni,N,expected",
+    [
+        (
+            # Intermediate values: NRD in (0, 1) -> no Laplace correction applied
+            np.array([2, 3], dtype=float),
+            np.array([5, 6], dtype=float),
+            np.array([2 / 5, 3 / 6]),
+        ),
+        (
+            # NRD=0 (Ni=0, N>0): Laplace correction -> (0+1) / (N+2)
+            np.array([0], dtype=float),
+            np.array([5], dtype=float),
+            np.array([1 / 7]),
+        ),
+        (
+            # NRD=1 (Ni=N, N>0): Laplace correction -> (N+1) / (N+2)
+            np.array([5], dtype=float),
+            np.array([5], dtype=float),
+            np.array([6 / 7]),
+        ),
+        (
+            # Ni=0, N=0: 0/0 -> NaN -> set to 0 -> Laplace correction -> 1/2
+            np.array([0], dtype=float),
+            np.array([0], dtype=float),
+            np.array([0.5]),
+        ),
+        (
+            # Mixed: covers NaN, NRD=0, intermediate, and NRD=1 in one array
+            np.array([0, 0, 2, 5], dtype=float),
+            np.array([0, 5, 5, 5], dtype=float),
+            np.array([0.5, 1 / 7, 0.4, 6 / 7]),
+        ),
+    ],
+    ids=[
+        "intermediate_no_correction",
+        "nrd_zero_laplace_correction",
+        "nrd_one_laplace_correction",
+        "nan_case_ni0_n0",
+        "mixed_all_cases",
+    ],
+)
+def test_compute_nrd(Ni, N, expected):
+    result = compute_nrd(Ni, N)
+    np.testing.assert_allclose(result, expected)
+
+
+# ── invariants ────────────────────────────────────────────────────────────────
+
+
+def test_output_strictly_between_zero_and_one():
+    rng = np.random.default_rng(42)
+    N = rng.integers(0, 20, size=100).astype(float)
+    Ni = np.array([rng.integers(0, int(n) + 1) for n in N], dtype=float)
+    result = compute_nrd(Ni, N)
+    assert np.all(result > 0)
+    assert np.all(result < 1)
+
+
+def test_no_nan_in_output():
+    Ni = np.array([0, 0, 3], dtype=float)
+    N = np.array([0, 5, 5], dtype=float)
+    result = compute_nrd(Ni, N)
+    assert not np.any(np.isnan(result))

test/pad_profile/test_main_pad_profile.py

@@ -34,7 +34,7 @@ def test_pad_profile_one_tile_real_las_returns_coherent_output_values():
         ground_margin=0.1,
     )
 
-    cos_theta, ni, n, min_layer = result
+    cos_theta, ni, n, min_layer, nrd = result
     assert isinstance(cos_theta, (float, int)), "Expected a numeric cos_theta value"
     assert 0.0 <= float(cos_theta) <= 1.0
-    assert len(ni) == len(n) == len(min_layer)
+    assert len(ni) == len(n) == len(min_layer) == len(nrd)