Vertical skill

src/modelskill/__init__.py

@@ -36,6 +36,7 @@
 from .model import (
     PointModelResult,
     TrackModelResult,
+    VerticalModelResult,
     GridModelResult,
     DfsuModelResult,
     DummyModelResult,
@@ -44,6 +45,7 @@
     observation,
     PointObservation,
     TrackObservation,
+    VerticalObservation,
 )
 from .matching import from_matched, match
 from .configuration import from_config
@@ -91,12 +93,14 @@ def load(filename: Union[str, Path]) -> Comparer | ComparerCollection:
     "model_result",
     "PointModelResult",
     "TrackModelResult",
+    "VerticalModelResult",
     "GridModelResult",
     "DfsuModelResult",
     "DummyModelResult",
     "observation",
     "PointObservation",
     "TrackObservation",
+    "VerticalObservation",
     "TimeSeries",
     "match",
     "from_matched",

src/modelskill/comparison/_comparer_plotter.py

@@ -746,7 +746,14 @@ def taylor(
         df = df.rename(columns={"_std_obs": "obs_std", "_std_mod": "std"})
 
         pts = [
-            TaylorPoint(name=r.model, obs_std=r.obs_std, std=r.std, cc=r.cc, marker=marker, marker_size=marker_size)
+            TaylorPoint(
+                name=r.model,
+                obs_std=r.obs_std,
+                std=r.std,
+                cc=r.cc,
+                marker=marker,
+                marker_size=marker_size,
+            )
             for r in df.itertuples()
         ]
 

src/modelskill/comparison/_comparison.py

@@ -25,7 +25,7 @@
 from .. import Quantity
 from ..types import GeometryType
 from ..obs import PointObservation, TrackObservation
-from ..model import PointModelResult, TrackModelResult
+from ..model import PointModelResult, TrackModelResult, VerticalModelResult
 from ..timeseries._timeseries import _validate_data_var_name
 from ._comparer_plotter import ComparerPlotter
 from ..metrics import _parse_metric
@@ -444,7 +444,10 @@ class Comparer:
     def __init__(
         self,
         matched_data: xr.Dataset,
-        raw_mod_data: dict[str, PointModelResult | TrackModelResult] | None = None,
+        raw_mod_data: dict[
+            str, PointModelResult | TrackModelResult | VerticalModelResult
+        ]
+        | None = None,
     ) -> None:
         self.data = _parse_dataset(matched_data)
         self.raw_mod_data = (
@@ -464,7 +467,9 @@ def __init__(
     @staticmethod
     def from_matched_data(
         data: xr.Dataset | pd.DataFrame,
-        raw_mod_data: Optional[Dict[str, PointModelResult | TrackModelResult]] = None,
+        raw_mod_data: Optional[
+            Dict[str, PointModelResult | TrackModelResult | VerticalModelResult]
+        ] = None,
         obs_item: str | int | None = None,
         mod_items: Optional[Iterable[str | int]] = None,
         aux_items: Optional[Iterable[str | int]] = None,
@@ -734,7 +739,9 @@ def _to_observation(self) -> PointObservation | TrackObservation:
         else:
             raise NotImplementedError(f"Unknown gtype: {self.gtype}")
 
-    def _to_model(self) -> list[PointModelResult | TrackModelResult]:
+    def _to_model(
+        self,
+    ) -> list[PointModelResult | TrackModelResult | VerticalModelResult]:
         mods = list(self.raw_mod_data.values())
         return mods
 
@@ -1257,8 +1264,13 @@ def load(filename: Union[str, Path]) -> "Comparer":
         if data.gtype == "track":
             return Comparer(matched_data=data)
 
+        if data.gtype == "vertical":
+            return Comparer(matched_data=data)  # FIXME: consider during Phase3
+
         if data.gtype == "point":
-            raw_mod_data: Dict[str, PointModelResult | TrackModelResult] = {}
+            raw_mod_data: Dict[
+                str, PointModelResult | TrackModelResult | VerticalModelResult
+            ] = {}
 
             for var in data.data_vars:
                 var_name = str(var)

src/modelskill/matching.py

@@ -28,7 +28,8 @@
 from .model.dummy import DummyModelResult
 from .model.grid import GridModelResult
 from .model.track import TrackModelResult
-from .obs import Observation, PointObservation, TrackObservation
+from .model.vertical import VerticalModelResult
+from .obs import Observation, PointObservation, TrackObservation, VerticalObservation
 from .timeseries import TimeSeries
 from .types import Period
 
@@ -40,6 +41,7 @@
     GridModelResult,
     DfsuModelResult,
     TrackModelResult,
+    VerticalModelResult,
     DummyModelResult,
 ]
 MRInputType = Union[
@@ -56,7 +58,7 @@
     TimeSeries,
     MRTypes,
 ]
-ObsTypes = Union[PointObservation, TrackObservation]
+ObsTypes = Union[PointObservation, TrackObservation, VerticalObservation]
 ObsInputType = Union[
     str,
     Path,
@@ -351,7 +353,9 @@ def _get_global_start_end(idxs: Iterable[pd.DatetimeIndex]) -> Period:
 
 def _match_space_time(
     observation: Observation,
-    raw_mod_data: Mapping[str, PointModelResult | TrackModelResult],
+    raw_mod_data: Mapping[
+        str, PointModelResult | TrackModelResult | VerticalModelResult
+    ],
     max_model_gap: float | None,
     spatial_tolerance: float,
     obs_no_overlap: Literal["ignore", "error", "warn"],
@@ -375,6 +379,9 @@ def _match_space_time(
                 )
             case PointModelResult() as pmr, PointObservation():
                 aligned = pmr.align(observation, max_gap=max_model_gap)
+            case VerticalModelResult(), VerticalObservation():
+                raise NotImplementedError("Vertical matching not implemented yet!")
+                # aligned = vmr.align(observation, max_gap=max_model_gap)
             case _:
                 raise TypeError(
                     f"Matching not implemented for model type {type(mr)} and observation type {type(observation)}"

src/modelskill/model/__init__.py

@@ -18,13 +18,15 @@
 from .factory import model_result
 from .point import PointModelResult
 from .track import TrackModelResult
+from .vertical import VerticalModelResult
 from .dfsu import DfsuModelResult
 from .grid import GridModelResult
 from .dummy import DummyModelResult
 
 __all__ = [
     "PointModelResult",
     "TrackModelResult",
+    "VerticalModelResult",
     "DfsuModelResult",
     "GridModelResult",
     "model_result",

src/modelskill/model/_base.py

@@ -10,9 +10,10 @@
 if TYPE_CHECKING:
     from .point import PointModelResult
     from .track import TrackModelResult
+    from .vertical import VerticalModelResult
 
 from ..utils import _get_name
-from ..obs import Observation, PointObservation, TrackObservation
+from ..obs import Observation, PointObservation, TrackObservation, VerticalObservation
 
 
 @dataclass
@@ -75,9 +76,9 @@ def _validate_overlap_in_time(time: pd.DatetimeIndex, observation: Observation)
 class SpatialField(Protocol):
     def extract(
         self,
-        observation: PointObservation | TrackObservation,
+        observation: PointObservation | TrackObservation | VerticalObservation,
         spatial_method: Optional[str] = None,
-    ) -> PointModelResult | TrackModelResult: ...
+    ) -> PointModelResult | TrackModelResult | VerticalModelResult: ...
 
     def _extract_point(
         self, observation: PointObservation, spatial_method: Optional[str] = None
@@ -86,3 +87,7 @@ def _extract_point(
     def _extract_track(
         self, observation: TrackObservation, spatial_method: Optional[str] = None
     ) -> TrackModelResult: ...
+
+    def _extract_vertical(
+        self, observation: VerticalObservation, spatial_method: Optional[str] = None
+    ) -> VerticalModelResult: ...

src/modelskill/model/dfsu.py

@@ -13,7 +13,8 @@
 from ..utils import _get_idx
 from .point import PointModelResult
 from .track import TrackModelResult
-from ..obs import Observation, PointObservation, TrackObservation
+from .vertical import VerticalModelResult
+from ..obs import Observation, PointObservation, TrackObservation, VerticalObservation
 
 
 class DfsuModelResult(SpatialField):
@@ -73,7 +74,9 @@ def __init__(
             item = data.name
             self.sel_items = SelectedItems(values=data.name, aux=[])
             data = mikeio.Dataset({data.name: data})
-        elif isinstance(data, (mikeio.dfsu.Dfsu2DH, mikeio.dfsu.Dfsu3D)):
+        elif isinstance(
+            data, (mikeio.dfsu.Dfsu2DH, mikeio.dfsu.Dfsu3D, mikeio.Dataset)
+        ):
             item_names = [i.name for i in data.items]
             idx = _get_idx(x=item, valid_names=item_names)
             item_info = data.items[idx]
@@ -113,14 +116,14 @@ def _in_domain(self, x: float, y: float) -> bool:
 
     def extract(
         self, observation: Observation, spatial_method: Optional[str] = None
-    ) -> PointModelResult | TrackModelResult:
+    ) -> PointModelResult | TrackModelResult | VerticalModelResult:
         """Extract ModelResult at observation positions
 
         Note: this method is typically not called directly, but by the match() method.
 
         Parameters
         ----------
-        observation : <PointObservation> or <TrackObservation>
+        observation : <PointObservation>, <TrackObservation>, or <VerticalObservation>
             positions (and times) at which modelresult should be extracted
         spatial_method : Optional[str], optional
             spatial selection/interpolation method, 'contained' (=isel),
@@ -129,7 +132,7 @@ def extract(
 
         Returns
         -------
-        PointModelResult or TrackModelResult
+        PointModelResult, TrackModelResult, or VerticalModelResult
             extracted modelresult with the same geometry as the observation
         """
         method = self._parse_spatial_method(spatial_method)
@@ -139,6 +142,8 @@ def extract(
             return self._extract_point(observation, spatial_method=method)
         elif isinstance(observation, TrackObservation):
             return self._extract_track(observation, spatial_method=method)
+        elif isinstance(observation, VerticalObservation):
+            return self._extract_vertical(observation, spatial_method=method)
         else:
             raise NotImplementedError(
                 f"Extraction from {type(self.data)} to {type(observation)} is not implemented."
@@ -162,6 +167,76 @@ def _parse_spatial_method(method: str | None) -> str | None:
         else:
             return METHOD_MAP[method]
 
+    def _extract_vertical(
+        self,
+        observation: VerticalObservation,
+        spatial_method: Optional[str] = None,
+    ) -> VerticalModelResult:
+        x, y = observation.x, observation.y
+
+        if not self._in_domain(x, y):
+            raise ValueError(
+                f"PointObservation '{observation.name}' (x={x}, y={y}) outside model domain!"
+            )
+
+        elemids = self.data.geometry.find_index(x=x, y=y)
+        if len(elemids) == 0:
+            raise ValueError(f"No elements found at observation location ({x}, {y})")
+
+        method = spatial_method or "contained"
+
+        # Check if 3D for dfsu and dataset
+        if not hasattr(self.data, "n_layers") and not hasattr(
+            self.data.geometry, "n_layers"
+        ):
+            raise ValueError("Cannot extract vertical profile from 2D dfsu file.")
+
+        if method == "contained":
+            if isinstance(self.data, mikeio.Dataset):
+                ds_column = self.data.sel(x=x, y=y)
+            elif isinstance(self.data, mikeio.dfsu.Dfsu3D):
+                # FIXME: open with specific element instead...make sure mikeio fix is in place before
+                # ds_column = self.data.read(elements=elemids) # when bug is fixed in mikeio
+                ds_column = self.data.read(items=self.sel_items.all).sel(x=x, y=y)
+            else:
+                raise ValueError(
+                    "Unsupported data type for vertical profile extraction."
+                )
+        else:
+            raise NotImplementedError(
+                "Only spatial_method='contained' is currently implemented for vertical profile extraction from DfsuModelResult. "
+            )
+
+        # get layer depth info
+        layer_boundaries = ds_column.geometry.calc_ze(ds_column._zn)
+
+        item_name = self.sel_items.values
+
+        # Flatten ds_column (time, layers) to 1D format for ProfileModelResult
+        # Need to repeat time for each layer and flatten z and values
+        n_layers = ds_column.geometry.n_layers
+
+        # Create flattened arrays # Repeat each timestamp n_layers times
+        time_flat = np.repeat(ds_column.time, n_layers)
+        z_flat = layer_boundaries.flatten()  # Flatten z-coordinates
+        item_values_1d = ds_column[item_name].to_numpy().flatten()  # Flatten item?
+        # aux_items_flat = {aux_item: ds_column[aux_item].to_numpy().flatten() for aux_item in self.sel_items.aux}
+
+        # Create DataFrame in long format
+        df = pd.DataFrame({"time": time_flat, "z": z_flat, self.name: item_values_1d})
+        # Set time as index for ProfileModelResult parsing
+        df = df.set_index("time")
+        return VerticalModelResult(
+            data=df,
+            item=self.name,
+            z_item="z",
+            x=ds_column.geometry.element_coordinates[0, 0],
+            y=ds_column.geometry.element_coordinates[0, 1],
+            name=self.name,
+            quantity=self.quantity,
+            aux_items=self.sel_items.aux,
+        )
+
     def _extract_point(
         self, observation: PointObservation, spatial_method: Optional[str] = None
     ) -> PointModelResult:

src/modelskill/model/dummy.py

@@ -6,7 +6,7 @@
 
 from modelskill.model.point import PointModelResult
 from modelskill.model.track import TrackModelResult
-from modelskill.obs import PointObservation, TrackObservation
+from modelskill.obs import PointObservation, TrackObservation, VerticalObservation
 
 
 @dataclass
@@ -49,14 +49,18 @@ def __post_init__(self):
 
     def extract(
         self,
-        observation: PointObservation | TrackObservation,
+        observation: PointObservation | TrackObservation | VerticalObservation,
         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"
             )
 
+        if isinstance(observation, VerticalObservation):
+            raise NotImplementedError(
+                "DummyModelResult does not support VerticalObservation yet"
+            )
         da = observation.data[observation.name].copy()
         if self.strategy == "mean":
             da[:] = da.mean()

src/modelskill/model/grid.py

@@ -12,7 +12,7 @@
 from ..quantity import Quantity
 from .point import PointModelResult
 from .track import TrackModelResult
-from ..obs import PointObservation, TrackObservation
+from ..obs import PointObservation, TrackObservation, VerticalObservation
 
 
 class GridModelResult(SpatialField):
@@ -124,7 +124,7 @@ def _in_domain(self, x: float, y: float) -> bool:
 
     def extract(
         self,
-        observation: PointObservation | TrackObservation,
+        observation: PointObservation | TrackObservation | VerticalObservation,
         spatial_method: Optional[str] = None,
     ) -> PointModelResult | TrackModelResult:
         """Extract ModelResult at observation positions
@@ -144,6 +144,10 @@ def extract(
         PointModelResult or TrackModelResult
             extracted modelresult
         """
+        if isinstance(observation, VerticalObservation):
+            raise NotImplementedError(
+                "Extraction of VerticalObservation from GridModelResult is not implemented yet."
+            )
         _validate_overlap_in_time(self.time, observation)
         if isinstance(observation, PointObservation):
             return self._extract_point(observation, spatial_method)