Untangled spatial and time resampling

Author: AntonioTorga

.gitignore

@@ -8,4 +8,5 @@ SUBPLAN-insitu-readers.md
 /.ruff_cache/
 **/DATOS-INSITU/*
 CLAUDE.md
-test_data/results/
+test_data/results/
+USAGE.md

ClimateGraph/data/data.py

@@ -4,7 +4,6 @@
 
 import cartopy.crs as ccrs
 import numpy as np
-import pandas as pd
 import xarray as xr
 
 from ClimateGraph.reader import Reader
@@ -118,7 +117,6 @@ def __init__(
         self._path = None
         self._bbox = None  # minlon, minlat, maxlon, maxlat
         self._dims = None
-        self.resampled = None
 
         # Provided by user
         self.name = name
@@ -145,7 +143,6 @@ def copy(self):
         new._geom = self._geom
         new._bbox = self._bbox  # minlon, minlat, maxlon, maxlat
         new._dims = self._dims
-        new.resampled = self.resampled
 
         return new
 
@@ -180,7 +177,6 @@ def obj(self, _obj: xr.Dataset):
         self._bbox = None
         self._geom = None
         self._dims = None
-        self.resampled = None
 
     @property
     def geom(self):
@@ -402,28 +398,27 @@ def resample_vars(
         other: "Data",
         vars: str | list[str],
         radius_of_influence: int = 10000,
-        time_tolerance: str | None = "30min",
         engine: str | ResampleEngine = "pyresample",
         engine_kwargs: dict | None = None,
-    ) -> xr.DataArray | xr.Dataset:
-        """Project ``other``'s vars onto ``self``'s geometry (and time axis).
+    ) -> xr.Dataset:
+        """Project ``other``'s vars onto ``self``'s spatial geometry.
 
-        A single responsibility: spatial resampling with the requested engine, plus
-        a nearest-neighbour time *alignment* of ``other`` onto ``self``'s time axis.
-        It does NOT convert units or filter/resample time — those belong to the
-        caller (``get_var`` / the plot's ``change_unit`` and ``time_resampling``).
+        Purely spatial and stateless: each variable of ``other`` is reprojected
+        onto ``self``'s geometry with the requested engine, keeping ``other``'s own
+        time axis and any extra dims (z, member, …) untouched, and inheriting
+        ``self``'s spatial coordinates (site / latitude / longitude / region / …).
+        It does NOT align time, convert units, or cache — those are the caller's
+        concern. Each call is independent, so two sources with different time
+        extents resampled onto the same target never clobber each other.
 
         Parameters
         ----------
         other : Data
-            Other data object to resample into the "self" geometry.
+            Data whose vars are reprojected onto ``self``'s geometry.
         vars : str | list[str]
-            Variable name or list of names to resample (kept in their source units).
+            Variable name or list of names to resample.
         radius_of_influence : int, optional
             Radius length in meters to use for resampling. by default 10000
-        time_tolerance : str | None, optional
-            Pandas-style timedelta used as the tolerance when snapping ``other``'s
-            time axis onto ``self``'s via nearest-neighbour reindex. Default ``"30min"``.
         engine : str | ResampleEngine, optional
             Resample backend name or instance. Default ``"pyresample"``.
         engine_kwargs : dict | None, optional
@@ -432,40 +427,37 @@ def resample_vars(
 
         Returns
         -------
-        xr.DataArray | xr.Dataset
-            Resampled data on ``self``'s geometry, aligned onto ``self``'s time axis.
+        xr.Dataset
+            One variable per input var, keyed ``"{var}__{other.name}"``, on
+            ``self``'s geometry and carrying ``other``'s time / extra dims.
         """
         if isinstance(vars, str):
             vars = [vars]
 
         resample_engine = get_engine(engine, **(engine_kwargs or {}))
-        src_geom = other.geom
-        dst_geom = self.geom
-
         info = resample_engine.prepare(
-            src_geom,
-            dst_geom,
+            other.geom,
+            self.geom,
             radius_of_influence=radius_of_influence,
         )
-        _resample = resample_engine.make_resampler(info, dst_geom.shape)
-
-        if self.resampled is None:
-            self.resampled = self.obj.drop_vars(list(self.obj.data_vars))
-
-        new_vars = []
+        _resample = resample_engine.make_resampler(info, self.geom.shape)
+
+        # Target geometry straight from self.obj (no data vars, no time axis): the
+        # geom-dim ORDER follows self.obj so it matches self.geom.shape, and the
+        # spatial coords (those depending only on the geom dims) are reattached to
+        # the result — this is how the source inherits site/region/lat/lon.
+        dst_geom_dims = [d for d in self.obj.dims if d in set(self.geom_dims)]
+        dst_sizes = {d: self.obj.sizes[d] for d in dst_geom_dims}
+        dst_coords = {
+            name: coord
+            for name, coord in self.obj.coords.items()
+            if set(coord.dims) <= set(self.geom_dims)
+        }
+
+        resampled_vars = {}
         for var in vars:
-            var_dst_dims = self.get_var(var).sizes
-            var_src = other.get_var(var)
-
-            if time_tolerance is not None and "time" in var_src.dims:
-                var_src = var_src.reindex(
-                    time=self.resampled["time"],
-                    method="nearest",
-                    tolerance=pd.Timedelta(time_tolerance),
-                )
-
+            var_src = other.get_var(var)  # keeps other's own time + extra dims
             src_geom_dims = [d for d in var_src.dims if d in set(other.geom_dims)]
-            dst_geom_dims = [d for d in var_dst_dims if d in set(self.geom_dims)]
             resampled = xr.apply_ufunc(
                 _resample,
                 var_src,
@@ -474,30 +466,22 @@ def resample_vars(
                 vectorize=True,
                 dask="parallelized",
                 output_dtypes=[var_src.dtype],
-                dask_gufunc_kwargs={
-                    "output_sizes": {
-                        name: value
-                        for name, value in var_dst_dims.items()
-                        if name in set(self.geom_dims)
-                    }
-                },
+                dask_gufunc_kwargs={"output_sizes": dst_sizes},
             )
+            resampled_vars[f"{var}__{other.name}"] = resampled.assign_coords(dst_coords)
 
-            new_name = f"{var}__{other.name}"
-            self.resampled[new_name] = resampled
-            new_vars.append(new_name)
-
+        result = xr.Dataset(resampled_vars)
         _record(
-            self.resampled,
-            f"spatially resampled {new_vars} from {other.name} "
+            result,
+            f"spatially resampled {list(resampled_vars)} from {other.name} "
             f"({type(other).__name__}) to {self.name} ({type(self).__name__})",
         )
 
         save_to = self.reader_kwargs.get("save_resampled_to")
         if save_to:
             target = Path(save_to)
             target.parent.mkdir(parents=True, exist_ok=True)
-            self.resampled.to_netcdf(target)
-            _record(self.resampled, f"saved resampled data to {target}")
+            result.to_netcdf(target)
+            _record(result, f"saved resampled data to {target}")
 
-        return self.resampled[new_vars]
+        return result

ClimateGraph/domain/domain.py

@@ -122,11 +122,10 @@ def apply(self, data: "Data") -> "Data":
     def _resample(self, data: "Data") -> "Data":
         """_resample Reproject ``data`` onto the resample target's geometry.
 
-        Resampling runs on the REAL target so its ``.resampled`` cache accumulates
-        and ``save_resampled_to`` fires as designed. The returned value is a cheap
-        target-topology wrapper wearing the SOURCE's identity (name + vars); building
-        it never disturbs the target's cache because the ``obj`` setter resets only
-        the wrapper's own caches.
+        ``resample_vars`` is purely spatial and stateless, so the reprojected
+        result keeps ``data``'s own time axis and inherits the target's spatial
+        coords (site / region / lat / lon). The returned value is a cheap
+        target-topology wrapper wearing the SOURCE's identity (name + vars).
 
         Parameters
         ----------

ClimateGraph/plot/plots.py

@@ -5,6 +5,7 @@
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 import numpy as np
+import xarray as xr
 from pydantic import BaseModel, ConfigDict, Field, field_validator, model_validator
 
 mpl.use("Agg")
@@ -324,6 +325,13 @@ def _plot_one(self, time_interval: str | None):
                     reduced[ds_name] = da
 
                 x_var, y_var = reduced[x_name], reduced[y_name]
+                # Scatter pairs the two datasets element-wise along `dimension`.
+                # resample_vars is now purely spatial (it no longer snaps time), so
+                # align the two on the shared axis here — the correct home for the
+                # time-alignment concern. inner-join keeps only matching coords.
+                if dimension in x_var.dims and dimension in y_var.dims:
+                    x_var, y_var = xr.align(x_var, y_var, join="inner")
+
                 figure = plt.figure(**self.figure_kwargs())
                 ax = figure.add_subplot(1, 1, 1)
                 min_val, max_val = (

test_data/configs/template_custom_timeseries.yaml

@@ -0,0 +1,84 @@
+# ============================================================================
+# TEMPLATE — custom timeseries (`type: custom` + `series` primitives)
+# ----------------------------------------------------------------------------
+# A "custom" plot overlays rendering primitives on shared axes. A `series` with
+# `x: time` is a line over time — space is reduced away automatically. Copy a
+# block below and adapt `dataset` / `var` / `time` to your data.
+#
+# Runs as-is on the bundled WRF (grid) + DMC (stations) samples.
+# ============================================================================
+analysis:
+  output_path: "./test_data/results/template-custom-ts/"
+  debug: false
+
+data:
+  DMC:                                   # in-situ stations (PointSurface)
+    path: "./test_data/data/dmc*.nc"
+    topology: PointSurface
+    reader: DMC
+    vars:
+      Temperatura: {name: temperatura, unit: degC}
+  WRF:                                   # model grid (RegularGrid)
+    path: "./test_data/data/wrf-20*.nc"
+    topology: RegularGrid
+    reader: wrf
+    vars:
+      Temperatura: {name: T2, unit: kelvin}
+
+plots:
+
+  # 1) SIMPLEST — one line, one dataset.
+  single:
+    type: custom
+    time: "1/1/2019 - 1/2/2019"
+    vars: [Temperatura]                  # one figure per var in the list
+    subplots:
+      - type: series
+        dataset: DMC
+        x: time
+        label: "Stations"
+
+  # 2) COMPARISON — two datasets, same var, on one figure.
+  #    Plot-level `vars` as a {var: unit} map converts every series to that unit
+  #    (WRF kelvin -> degC) so they share a y-axis. `grid: true` adds faded
+  #    dashed reference lines at each tick.
+  obs_vs_model:
+    type: custom
+    time: "1/1/2019 - 1/2/2019"
+    timestep: D                          # daily-mean before plotting (optional)
+    vars: {Temperatura: degC}
+    grid: true
+    subplots:
+      - type: series
+        dataset: DMC
+        x: time
+        label: "Observations"
+        color: black
+        linewidth: 2
+      - type: series
+        dataset: WRF
+        x: time
+        label: "WRF model"
+        color: tab:red
+
+  # 3) TWO PERIODS on one figure — per-subplot `time` override.
+  #    Omit plot-level `vars` and give each subplot its own `var` + `time`; the
+  #    plot renders once (no fan-out) with both windows overlaid.
+  two_periods:
+    type: custom
+    subplots:
+      - type: series
+        dataset: DMC
+        var: Temperatura
+        x: time
+        time: "1/1/2019 - 1/2/2019"
+        label: "January"
+        color: tab:blue
+      - type: series
+        dataset: DMC
+        var: Temperatura
+        x: time
+        time: "1/2/2019 - 28/2/2019"
+        label: "February"
+        color: tab:orange
+    title: "DMC temperature — January vs February"

tests/test_data_module/test_data_methods.py

@@ -94,12 +94,34 @@ def test_resample_grid_to_point(self, regular_grid_data, point_surface_data):
         # Result should live on the destination (point) geometry.
         assert "site" in result.dims
 
-    def test_resample_caches_resampled(self, regular_grid_data, point_surface_data):
-        assert point_surface_data.resampled is None
-        point_surface_data.resample_vars(
+    def test_resample_is_stateless(self, regular_grid_data, point_surface_data):
+        # No cache: resampling never mutates the target, and each call is a fresh
+        # independent Dataset (so different-time sources can't clobber each other).
+        assert not hasattr(point_surface_data, "resampled")
+        target_before = point_surface_data.obj
+        r1 = point_surface_data.resample_vars(
             regular_grid_data, "Temperatura", radius_of_influence=500_000
         )
-        assert point_surface_data.resampled is not None
+        r2 = point_surface_data.resample_vars(
+            regular_grid_data, "Temperatura", radius_of_influence=500_000
+        )
+        assert point_surface_data.obj is target_before  # target untouched
+        assert r1 is not r2
+
+    def test_resample_keeps_source_time_and_inherits_target_coords(
+        self, regular_grid_data, point_surface_data
+    ):
+        # The source keeps its OWN time axis (no snapping onto the target), and
+        # inherits the target's spatial coords (site + region), which is what makes
+        # region/site attribute filtering on a resampled grid possible.
+        result = point_surface_data.resample_vars(
+            regular_grid_data, "Temperatura", radius_of_influence=500_000
+        )
+        src_time = regular_grid_data.obj["time"]
+        assert result.sizes["time"] == src_time.sizes["time"]
+        assert result["time"].equals(src_time)
+        assert "region" in result.coords  # inherited from the point target
+        assert result.sizes["site"] == point_surface_data.obj.sizes["site"]
 
     def test_resample_4d_grid_preserves_z(self, point_surface_data):
         """Resampling a 4D (time, z, y, x) RegularGrid should broadcast across z
@@ -173,13 +195,11 @@ def test_cached_state_is_propagated(self, regular_grid_data):
         _ = regular_grid_data.bbox
         _ = regular_grid_data.dims
         regular_grid_data._set_geom()
-        regular_grid_data.resampled = "sentinel"
 
         copy = regular_grid_data.copy()
         assert copy._bbox == regular_grid_data._bbox
         assert copy._dims == regular_grid_data._dims
         assert copy._geom is regular_grid_data._geom
-        assert copy.resampled == "sentinel"
 
     def test_obj_setter_on_copy_does_not_mutate_original(
         self, regular_grid_data, make_regular_grid

tests/test_plot/test_plot_runs.py

@@ -336,3 +336,29 @@ def test_runs_with_grid_and_point(
         )
         plot.plot()
         assert _outputs(tmp_output_dir, "sc")
+
+    def test_aligns_datasets_with_mismatched_time(
+        self, point_surface_data, tmp_output_dir
+    ):
+        # resample_vars no longer snaps time, so Scatter aligns the two series
+        # itself. Pair a 6-step dataset with a 4-step one: without the inner-join
+        # align the reduced arrays would be different lengths and scatter would
+        # raise; with it, they pair on the overlapping timestamps.
+        short = point_surface_data.copy()
+        short.obj = point_surface_data.obj.isel(time=slice(0, 4))
+
+        cfg = ScatterConfig(
+            type="scatter",
+            data=["A", "B"],
+            time=TIME_INTERVAL,
+            vars={"Temperatura": "degC"},
+        )
+        plot = Scatter(
+            name="sc_align",
+            plot_config=cfg,
+            data_registry={"A": point_surface_data, "B": short},
+            domain_registry={},
+            output_path=tmp_output_dir,
+        )
+        plot.plot()  # must not raise on the length mismatch
+        assert _outputs(tmp_output_dir, "sc_align")