Merge pull request #383 from DHI/spatial-interp-method

Spatial interp method

modelskill/connection.py

@@ -183,17 +183,15 @@ def _parse_observation(self, obs):
         else:
             raise ValueError(f"Unknown observation type {type(obs)}")
 
-    def extract(self, max_model_gap: Optional[float] = None) -> Optional[Comparer]:
+    def extract(self, **kwargs) -> Optional[Comparer]:
         """Extract model results at times and positions of observation.
 
         Returns
         -------
         Comparer
             A comparer object for further analysis and plotting.
         """
-        comparer = _single_obs_compare(
-            obs=self.obs, mod=self.modelresults, max_model_gap=max_model_gap
-        )
+        comparer = _single_obs_compare(obs=self.obs, mod=self.modelresults, **kwargs)
         if comparer.n_points == 0:
             warnings.warn(f"No overlapping data was found for {comparer.name}!")
             return None

modelskill/matching.py

@@ -6,6 +6,7 @@
 from typing import (
     Dict,
     Iterable,
+    Collection,
     List,
     Literal,
     Optional,
@@ -162,6 +163,7 @@ def match(
     mod_item: Optional[IdxOrNameTypes] = None,
     gtype: Optional[GeometryTypes] = None,
     max_model_gap: Optional[float] = None,
+    spatial_method: Optional[str] = None,
 ) -> Comparer:
     ...
 
@@ -175,6 +177,7 @@ def match(
     mod_item: Optional[IdxOrNameTypes] = None,
     gtype: Optional[GeometryTypes] = None,
     max_model_gap: Optional[float] = None,
+    spatial_method: Optional[str] = None,
 ) -> ComparerCollection:
     ...
 
@@ -187,6 +190,7 @@ def match(
     mod_item=None,
     gtype=None,
     max_model_gap=None,
+    spatial_method: Optional[str] = None,
 ):
     """Match observation and model result data in space and time
 
@@ -213,6 +217,13 @@ def match(
     max_model_gap : (float, optional)
         Maximum time gap (s) in the model result (e.g. for event-based
         model results), by default None
+    spatial_method : str, optional
+        For Dfsu- and GridModelResult, spatial interpolation/selection method.
+
+        - For DfsuModelResult, one of: 'contained' (=isel), 'nearest',
+        'inverse_distance' (with 5 nearest points), by default "inverse_distance".
+        - For GridModelResult, passed to xarray.interp() as method argument,
+        by default 'linear'.
 
     Returns
     -------
@@ -234,11 +245,12 @@ def match(
             mod_item=mod_item,
             gtype=gtype,
             max_model_gap=max_model_gap,
+            spatial_method=spatial_method,
         )
 
-    assert isinstance(obs, Iterable)
+    assert isinstance(obs, Collection)
 
-    if len(obs) > 1 and isinstance(mod, Iterable) and len(mod) > 1:
+    if len(obs) > 1 and isinstance(mod, Collection) and len(mod) > 1:
         if not all(isinstance(m, (DfsuModelResult, GridModelResult)) for m in mod):
             raise ValueError(
                 """
@@ -260,6 +272,7 @@ def match(
             mod_item=mod_item,
             gtype=gtype,
             max_model_gap=max_model_gap,
+            spatial_method=spatial_method,
         )
         for o in obs
     ]
@@ -303,13 +316,14 @@ def _single_obs_compare(
     mod_item: Optional[int | str] = None,
     gtype: Optional[GeometryTypes] = None,
     max_model_gap: Optional[float] = None,
+    spatial_method: Optional[str] = None,
 ) -> Comparer:
     """Compare a single observation with multiple models"""
     obs = _parse_single_obs(obs, obs_item, gtype=gtype)
 
     mods = _parse_models(mod, mod_item, gtype=gtype)
 
-    raw_mod_data = {m.name: m.extract(obs) for m in mods}
+    raw_mod_data = {m.name: m.extract(obs, spatial_method) for m in mods}
     matched_data = match_space_time(obs, raw_mod_data, max_model_gap)
     matched_data.attrs["weight"] = obs.weight
 

modelskill/model/_base.py

@@ -73,12 +73,18 @@ def _validate_overlap_in_time(time: pd.DatetimeIndex, observation: Observation)
 
 class SpatialField(Protocol):
     def extract(
-        self, observation: PointObservation | TrackObservation
+        self,
+        observation: PointObservation | TrackObservation,
+        spatial_method: Optional[str] = None,
     ) -> PointModelResult | TrackModelResult:
         ...
 
-    def extract_point(self, observation: PointObservation) -> PointModelResult:
+    def _extract_point(
+        self, observation: PointObservation, spatial_method: Optional[str] = None
+    ) -> PointModelResult:
         ...
 
-    def extract_track(self, observation: TrackObservation) -> TrackModelResult:
+    def _extract_track(
+        self, observation: TrackObservation, spatial_method: Optional[str] = None
+    ) -> TrackModelResult:
         ...

modelskill/model/dfsu.py

@@ -98,32 +98,65 @@ def time(self) -> pd.DatetimeIndex:
     def _in_domain(self, x: float, y: float) -> bool:
         return self.data.geometry.contains([x, y])  # type: ignore
 
-    def extract(self, observation: Observation) -> PointModelResult | TrackModelResult:
+    def extract(
+        self, observation: Observation, spatial_method: Optional[str] = None
+    ) -> PointModelResult | TrackModelResult:
         """Extract ModelResult at observation positions
 
+        Note: this method is typically not called directly, but by the match() method.
+
         Parameters
         ----------
         observation : <PointObservation> or <TrackObservation>
             positions (and times) at which modelresult should be extracted
+        spatial_method : Optional[str], optional
+            spatial selection/interpolation method, 'contained' (=isel),
+            'nearest', 'inverse_distance' (with 5 nearest points),
+            by default None = 'inverse_distance'
 
         Returns
         -------
         PointModelResult or TrackModelResult
-            extracted modelresult
+            extracted modelresult with the same geometry as the observation
         """
+        method = self._parse_spatial_method(spatial_method)
+
         _validate_overlap_in_time(self.time, observation)
         if isinstance(observation, PointObservation):
-            return self.extract_point(observation)
+            return self._extract_point(observation, spatial_method=method)
         elif isinstance(observation, TrackObservation):
-            return self.extract_track(observation)
+            return self._extract_track(observation, spatial_method=method)
         else:
             raise NotImplementedError(
                 f"Extraction from {type(self.data)} to {type(observation)} is not implemented."
             )
 
-    def extract_point(self, observation: PointObservation) -> PointModelResult:
-        """Spatially extract a PointModelResult from a DfsuModelResult (when data is a Dfsu object),
-        given a PointObservation. No time interpolation is done!"""
+    @staticmethod
+    def _parse_spatial_method(method):
+        if method == "isel":
+            method = "contained"
+        if method == "IDW":
+            method = "inverse_distance"
+        if method is None:
+            return None  # decide default later
+        elif method not in ["nearest", "contained", "inverse_distance"]:
+            raise ValueError(
+                f"spatial_method for Dfsu must be 'nearest', 'contained', or 'inverse_distance'. Not {method}."
+            )
+        return method
+
+    def _extract_point(
+        self, observation: PointObservation, spatial_method: Optional[str] = None
+    ) -> PointModelResult:
+        """Spatially extract a PointModelResult from a DfsuModelResult
+        given a PointObservation. No time interpolation is done!
+
+        Note: 'inverse_distance' method uses 5 nearest points and is the default.
+        """
+
+        method = spatial_method or "inverse_distance"
+        assert method in ["nearest", "contained", "inverse_distance"]
+        n_nearest = 5 if method == "inverse_distance" else 1
 
         assert isinstance(
             self.data, (mikeio.dfsu.Dfsu2DH, mikeio.DataArray, mikeio.Dataset)
@@ -135,19 +168,36 @@ def extract_point(self, observation: PointObservation) -> PointModelResult:
                 f"PointObservation '{observation.name}' ({x}, {y}) outside model domain!"
             )
 
-        # TODO: interp2d
-        xy = np.atleast_2d([x, y])
-        elemids = self.data.geometry.find_index(coords=xy)
-        if isinstance(self.data, mikeio.dfsu.Dfsu2DH):
-            ds_model = self.data.read(elements=elemids, items=self.sel_items.all)
-            aux_items = self.sel_items.aux
-        elif isinstance(self.data, mikeio.Dataset):
-            ds_model = self.data[self.sel_items.all].isel(element=elemids)
-            aux_items = self.sel_items.aux
-        elif isinstance(self.data, mikeio.DataArray):
-            da = self.data.isel(element=elemids)
-            ds_model = mikeio.Dataset({da.name: da})
-            aux_items = None
+        if method == "contained":
+            xy = np.atleast_2d([x, y])
+            elemids = self.data.geometry.find_index(coords=xy)
+            if isinstance(self.data, mikeio.dfsu.Dfsu2DH):
+                ds_model = self.data.read(elements=elemids, items=self.sel_items.all)
+            elif isinstance(self.data, mikeio.Dataset):
+                ds_model = self.data[self.sel_items.all].isel(element=elemids)
+            elif isinstance(self.data, mikeio.DataArray):
+                da = self.data.isel(element=elemids)
+                ds_model = mikeio.Dataset({da.name: da})
+        else:
+            if isinstance(self.data, mikeio.dfsu.Dfsu2DH):
+                elemids = self.data.geometry.find_nearest_elements(
+                    x, y, n_nearest=n_nearest
+                )
+                ds = self.data.read(elements=elemids, items=self.sel_items.all)
+                ds_model = (
+                    ds.interp(x=x, y=y, n_nearest=n_nearest) if n_nearest > 1 else ds
+                )
+            elif isinstance(self.data, mikeio.Dataset):
+                ds_model = self.data[self.sel_items.all].interp(
+                    x=x, y=y, n_nearest=n_nearest
+                )
+            elif isinstance(self.data, mikeio.DataArray):
+                da = self.data.interp(x=x, y=y, n_nearest=n_nearest)
+                ds_model = mikeio.Dataset({da.name: da})
+
+        aux_items = (
+            None if isinstance(self.data, mikeio.DataArray) else self.sel_items.aux
+        )
 
         assert isinstance(ds_model, mikeio.Dataset)
 
@@ -165,9 +215,22 @@ def extract_point(self, observation: PointObservation) -> PointModelResult:
             aux_items=aux_items,
         )
 
-    def extract_track(self, observation: TrackObservation) -> TrackModelResult:
+    def _extract_track(
+        self, observation: TrackObservation, spatial_method: Optional[str] = None
+    ) -> TrackModelResult:
         """Extract a TrackModelResult from a DfsuModelResult (when data is a Dfsu object),
-        given a TrackObservation."""
+        given a TrackObservation.
+
+        Wraps MIKEIO's extract_track method (which has the default method='nearest').
+
+        MIKE IO's extract_track, inverse_distance method, uses 5 nearest points.
+        """
+        method = spatial_method or "inverse_distance"
+        if method == "contained":
+            raise NotImplementedError(
+                "spatial method 'contained' (=isel) not implemented for track extraction in MIKE IO"
+            )
+        assert method in ["nearest", "inverse_distance"]
 
         assert isinstance(
             self.data, (mikeio.dfsu.Dfsu2DH, mikeio.DataArray, mikeio.Dataset)
@@ -176,16 +239,18 @@ def extract_track(self, observation: TrackObservation) -> TrackModelResult:
         track = observation.data.to_dataframe()
 
         if isinstance(self.data, mikeio.DataArray):
-            ds_model = self.data.extract_track(track=track)
+            ds_model = self.data.extract_track(track=track, method=method)
             ds_model.rename({self.data.name: self.name}, inplace=True)
             aux_items = None
         else:
             if isinstance(self.data, mikeio.dfsu.Dfsu2DH):
                 ds_model = self.data.extract_track(
-                    track=track, items=self.sel_items.all
+                    track=track, items=self.sel_items.all, method=method
                 )
             elif isinstance(self.data, mikeio.Dataset):
-                ds_model = self.data[self.sel_items.all].extract_track(track=track)
+                ds_model = self.data[self.sel_items.all].extract_track(
+                    track=track, method=method
+                )
             ds_model.rename({self.sel_items.values: self.name}, inplace=True)
             aux_items = self.sel_items.aux
 

modelskill/model/dummy.py

@@ -1,5 +1,5 @@
 from __future__ import annotations
-from typing import Literal
+from typing import Literal, Optional
 
 import pandas as pd
 
@@ -54,26 +54,27 @@ def __init__(
         self.strategy = strategy
 
     def extract(
-        self, observation: PointObservation | TrackObservation
+        self,
+        observation: PointObservation | TrackObservation,
+        spatial_method: Optional[str] = None,
     ) -> PointModelResult | TrackModelResult:
+        if spatial_method is not None:
+            raise NotImplementedError(
+                "spatial interpolation not possible when matching point model results with point observations"
+            )
+
+        da = observation.data[observation.name].copy()
+        if self.strategy == "mean":
+            da[:] = da.mean()
+        else:
+            da[:] = self.data
+
         if isinstance(observation, PointObservation):
-            da = observation.data[observation.name].copy()
-            if self.strategy == "mean":
-                da[:] = da.mean()
-            else:
-                da[:] = self.data
-            pmr = PointModelResult(
+            return PointModelResult(
                 data=da, x=observation.x, y=observation.y, name=self.name
             )
-            return pmr
-
-        if isinstance(observation, TrackObservation):
-            da = observation.data[observation.name].copy()
-            if self.strategy == "mean":
-                da[:] = da.mean()
-            else:
-                da[:] = self.data
 
+        elif isinstance(observation, TrackObservation):
             data = pd.DataFrame(
                 {
                     "x": observation.x,
@@ -82,5 +83,8 @@ def extract(
                 },
                 index=da.time,
             )
-            tmr = TrackModelResult(data=data, name=self.name)
-            return tmr
+            return TrackModelResult(data=data, name=self.name)
+        else:
+            raise ValueError(
+                f"observation must be a PointObservation or TrackObservation not {type(observation)}"
+            )