8/8 Add bounded forward filling for satellite indices

pyVPRM/VPRM.py

@@ -582,6 +582,56 @@ def sort_and_merge_by_timestamp(self, min_length_snow_period=21):
                 self.sat_imgs.sat_img[sat_ind] = da.where(np.isfinite(da))
         return
 
+    def forward_fill_satellite_indices(self, keys=("evi", "lswi"), max_age_days=None):
+        """Forward-fill missing satellite-index values along composite time.
+
+        Parameters
+        ----------
+        keys : sequence of str, default=("evi", "lswi")
+            Satellite-index variables to fill. Each must use the current
+            ``time_key`` dimension.
+        max_age_days : float, optional
+            Maximum permitted age in days of a carried-forward observation.
+            Values with no earlier valid observation, or whose earlier value
+            exceeds this age, remain missing. ``None`` permits an unlimited
+            forward fill.
+
+        Returns
+        -------
+        None
+            Selected index variables are updated in :attr:`sat_imgs`.
+
+        Raises
+        ------
+        ValueError
+            If ``max_age_days`` is negative or a requested index is absent.
+
+        Notes
+        -----
+        ``sort_and_merge_by_timestamp`` represents composite time as elapsed
+        days from the first image. Call this method after merging satellite
+        images and before calculating index-derived VPRM statistics.
+        """
+        if max_age_days is not None and max_age_days < 0:
+            raise ValueError("max_age_days must be non-negative or None.")
+
+        satellite_dataset = self.sat_imgs.sat_img
+        time_values = satellite_dataset[self.time_key]
+        for key in keys:
+            if key not in satellite_dataset:
+                raise ValueError("Satellite index is unavailable: {}.".format(key))
+            index_values = satellite_dataset[key]
+            filled_values = index_values.ffill(dim=self.time_key)
+            if max_age_days is not None:
+                last_valid_time = xr.where(
+                    index_values.notnull(), time_values, np.nan
+                ).ffill(dim=self.time_key)
+                filled_values = filled_values.where(
+                    (time_values - last_valid_time) <= max_age_days
+                )
+            satellite_dataset[key] = filled_values
+        return
+
     def clip_to_box(self, sat_to_crop):
         bounds = sat_to_crop.sat_img.rio.bounds()
         self.sat_imgs.sat_img = self.sat_imgs.sat_img.rio.clip_box(

tests/test_satellite_index_forward_fill.py

@@ -0,0 +1,63 @@
+"""Tests for bounded forward filling of satellite indices."""
+
+import numpy as np
+import xarray as xr
+
+from pyVPRM.VPRM import vprm_preprocessor
+from pyVPRM.sat_managers.base_manager import satellite_data_manager
+
+
+def test_forward_fill_satellite_indices_respects_maximum_age():
+    """Fill only gaps supported by a sufficiently recent observation.
+
+    Returns
+    -------
+    None
+        The test verifies per-pixel EVI and LSWI filling, preservation of
+        leading gaps, and rejection of observations older than the limit.
+    """
+    preprocessor = object.__new__(vprm_preprocessor)
+    preprocessor.time_key = "time"
+    satellite_indices = xr.Dataset(
+        {
+            "evi": (
+                ("time", "y", "x"),
+                np.array(
+                    [
+                        [[0.2, np.nan], [np.nan, np.nan]],
+                        [[np.nan, 0.4], [np.nan, np.nan]],
+                        [[np.nan, np.nan], [np.nan, np.nan]],
+                        [[0.8, np.nan], [np.nan, np.nan]],
+                    ]
+                ),
+            ),
+            "lswi": (
+                ("time", "y", "x"),
+                np.array(
+                    [
+                        [[0.1, np.nan], [np.nan, np.nan]],
+                        [[np.nan, 0.3], [np.nan, np.nan]],
+                        [[np.nan, np.nan], [np.nan, np.nan]],
+                        [[0.7, np.nan], [np.nan, np.nan]],
+                    ]
+                ),
+            ),
+        },
+        coords={
+            "time": [0.0, 8.0, 16.0, 32.0],
+            "y": [200.0, 100.0],
+            "x": [100.0, 200.0],
+        },
+    )
+    preprocessor.sat_imgs = satellite_data_manager(sat_img=satellite_indices)
+
+    preprocessor.forward_fill_satellite_indices(max_age_days=12)
+
+    filled_evi = preprocessor.sat_imgs.sat_img["evi"]
+    filled_lswi = preprocessor.sat_imgs.sat_img["lswi"]
+    assert filled_evi.sel(time=8.0, y=200.0, x=100.0).item() == 0.2
+    assert np.isnan(filled_evi.sel(time=16.0, y=200.0, x=100.0).item())
+    assert filled_evi.sel(time=16.0, y=200.0, x=200.0).item() == 0.4
+    assert np.isnan(filled_evi.sel(time=32.0, y=200.0, x=200.0).item())
+    assert np.isnan(filled_evi.sel(time=0.0, y=200.0, x=200.0).item())
+    assert filled_lswi.sel(time=8.0, y=200.0, x=100.0).item() == 0.1