From c9639b17ad98eee1c0921f3843575ac69db807fe Mon Sep 17 00:00:00 2001 From: "Timothy W. Hilton" Date: Sat, 15 Aug 2026 16:20:59 +1200 Subject: [PATCH] Add bounded satellite-index forward fill Provide an opt-in per-pixel forward-fill method for EVI, LSWI, or caller-selected satellite indices after timestamp merging. Preserve leading and over-age gaps as missing through an optional observation-age limit, and add regression coverage for independent pixel histories and bounded carry-forward behavior. Runs much faster than LOWESS. --- pyVPRM/VPRM.py | 50 +++++++++++++++++ tests/test_satellite_index_forward_fill.py | 63 ++++++++++++++++++++++ 2 files changed, 113 insertions(+) create mode 100644 tests/test_satellite_index_forward_fill.py diff --git a/pyVPRM/VPRM.py b/pyVPRM/VPRM.py index 6fbb34e..4ca1725 100644 --- a/pyVPRM/VPRM.py +++ b/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( diff --git a/tests/test_satellite_index_forward_fill.py b/tests/test_satellite_index_forward_fill.py new file mode 100644 index 0000000..78e4df9 --- /dev/null +++ b/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