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
@@ -175,6 +175,7 @@ class ItemSelection:
     aux: Sequence[str]
     x: Optional[str] = None
     y: Optional[str] = None
+    z: Optional[str] = None
 
     def __post_init__(self) -> None:
         # check that obs, model and aux are unique, and that they are not overlapping
@@ -189,6 +190,8 @@ def all(self) -> Sequence[str]:
             res.append(self.x)
         if self.y is not None:
             res.append(self.y)
+        if self.z is not None:
+            res.append(self.z)
         return res
 
     @staticmethod
@@ -199,6 +202,7 @@ def parse(
         aux_items: Optional[Iterable[str | int]] = None,
         x_item: str | int | None = None,
         y_item: str | int | None = None,
+        z_item: str | int | None = None,
     ) -> ItemSelection:
         """Parse item selection.
 
@@ -214,26 +218,30 @@ def parse(
             raise ValueError("data must contain at least two items")
         obs_name = _get_name(obs_item, items) if obs_item else items[0]
 
-        # Check existance of items and convert to names
-
+        # Check existence of items and convert to names
         if aux_items is not None:
             if isinstance(aux_items, (str, int)):
                 aux_items = [aux_items]
             aux_names = [_get_name(a, items) for a in aux_items]
         else:
             aux_names = []
+
+        x_name = _get_name(x_item, items) if x_item is not None else None
+        y_name = _get_name(y_item, items) if y_item is not None else None
+        z_name = _get_name(z_item, items) if z_item is not None else None
+
         if mod_items is not None:
             if isinstance(mod_items, (str, int)):
                 mod_items = [mod_items]
             mod_names = [_get_name(m, items) for m in mod_items]
         else:
+            # Add remaining items as model items
             mod_names = [
-                item for item in items if item not in aux_names and item != obs_name
+                item
+                for item in items
+                if item not in [x_name, y_name, z_name] + aux_names + [obs_name]
             ]
 
-        x_name = _get_name(x_item, items) if x_item is not None else None
-        y_name = _get_name(y_item, items) if y_item is not None else None
-
         if len(mod_names) == 0:
             raise ValueError("no model items were found! Must be at least one")
         if obs_name in mod_names:
@@ -242,7 +250,7 @@ def parse(
             raise ValueError("observation item must not be an auxiliary item")
 
         return ItemSelection(
-            obs=obs_name, model=mod_names, aux=aux_names, x=x_name, y=y_name
+            obs=obs_name, model=mod_names, aux=aux_names, x=x_name, y=y_name, z=z_name
         )
 
 
@@ -302,18 +310,23 @@ def _matched_data_to_xarray(
     z: Optional[float] = None,
     x_item: str | int | None = None,
     y_item: str | int | None = None,
+    z_item: str | int | None = None,
     quantity: Optional[Quantity] = None,
 ) -> xr.Dataset:
     """Convert matched data to accepted xarray.Dataset format"""
     assert isinstance(df, pd.DataFrame)
-    cols = list(df.columns)
-    items = ItemSelection.parse(cols, obs_item, mod_items, aux_items, x_item, y_item)
 
+    cols = list(df.columns)
+    items = ItemSelection.parse(
+        cols, obs_item, mod_items, aux_items, x_item, y_item, z_item
+    )
     # check that x and x_item is not both specified (same for y and y_item)
     if (x is not None) and (x_item is not None):
         raise ValueError("x and x_item cannot both be specified")
     if (y is not None) and (y_item is not None):
         raise ValueError("y and y_item cannot both be specified")
+    if (z is not None) and (z_item is not None):
+        raise ValueError("z and z_item cannot both be specified")
 
     if x is not None:
         try:
@@ -329,7 +342,6 @@ def _matched_data_to_xarray(
             raise TypeError(
                 f"y must be scalar, not {type(y)}; if y is a coordinate variable, use y_item"
             )
-
     # check that items.obs and items.model are numeric
     if not np.issubdtype(df[items.obs].dtype, np.number):
         raise ValueError(
@@ -370,16 +382,21 @@ def _matched_data_to_xarray(
     else:
         ds.coords["y"] = np.nan
 
-    # No z-item so far (relevant for ProfileObservation)
-    if z is not None:
+    if z_item is not None:
+        ds = ds.rename({items.z: "z"}).set_coords("z")
+    elif z is not None:
         ds.coords["z"] = z
+    elif "z" in ds.data_vars:
+        ds = ds.set_coords("z")
+    else:
+        pass  # z is optional, if not provided, we don't add it as a coordinate
 
-    if ds.coords["x"].size == 1:
+    if "z" in ds.coords:
+        ds.attrs["gtype"] = str(GeometryType.VERTICAL)
+    elif ds.coords["x"].size == 1:
         ds.attrs["gtype"] = str(GeometryType.POINT)
     else:
         ds.attrs["gtype"] = str(GeometryType.TRACK)
-    # TODO
-    # ds.attrs["gtype"] = str(GeometryType.PROFILE)
 
     if quantity is None:
         q = Quantity.undefined()
@@ -444,7 +461,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 +484,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,
@@ -475,6 +497,7 @@ def from_matched_data(
         z: Optional[float] = None,
         x_item: str | int | None = None,
         y_item: str | int | None = None,
+        z_item: str | int | None = None,
         quantity: Optional[Quantity] = None,
     ) -> "Comparer":
         """Initialize from compared data"""
@@ -491,9 +514,11 @@ def from_matched_data(
                 z=z,
                 x_item=x_item,
                 y_item=y_item,
+                z_item=z_item,
                 quantity=quantity,
             )
             data.attrs["weight"] = weight
+            # FIXME: Consider changes needed for vertical
         return Comparer(matched_data=data, raw_mod_data=raw_mod_data)
 
     def __repr__(self):
@@ -734,7 +759,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
 
@@ -1193,6 +1220,7 @@ def to_dataframe(self) -> pd.DataFrame:
         """Convert matched data to pandas DataFrame
 
         Include x, y coordinates only if gtype=track
+        Include z coordinate only if gtype=vertical
 
         Returns
         -------
@@ -1208,6 +1236,16 @@ def to_dataframe(self) -> pd.DataFrame:
             # make sure that x, y cols are first
             cols = ["x", "y"] + [c for c in df.columns if c not in ["x", "y"]]
             return df[cols]
+        elif self.gtype == str(GeometryType.VERTICAL):
+            df = self.data.drop_vars(["x", "y"]).to_dataframe()
+            # make sure that Observations is first and z is last
+            cols = (
+                ["Observation"]
+                + [c for c in df.columns if c not in ["Observation", "z"]]
+                + ["z"]
+            )
+            return df[cols]
+
         else:
             raise NotImplementedError(f"Unknown gtype: {self.gtype}")
 
@@ -1221,6 +1259,8 @@ def save(self, filename: Union[str, Path]) -> None:
         """
         ds = self.data
 
+        # FIXME: Consider changes needed for vertical
+
         # add self.raw_mod_data to ds with prefix 'raw_' to avoid name conflicts
         # an alternative strategy would be to use NetCDF groups
         # https://docs.xarray.dev/en/stable/user-guide/io.html#groups
@@ -1257,8 +1297,14 @@ def load(filename: Union[str, Path]) -> "Comparer":
         if data.gtype == "track":
             return Comparer(matched_data=data)
 
+        if data.gtype == "vertical":
+            # FIXME: consider during Phase3
+            return Comparer(matched_data=data)
+
         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,
@@ -85,6 +87,7 @@ def from_matched(
     z: Optional[float] = None,
     x_item: str | int | None = None,
     y_item: str | int | None = None,
+    z_item: str | int | None = None,
 ) -> Comparer:
     """Create a Comparer from data that is already matched (aligned).
 
@@ -113,6 +116,8 @@ def from_matched(
         Name of x item, only relevant for track data
     y_item: [str, int], optional
         Name of y item, only relevant for track data
+    z_item: [str, int], optional
+        Name of z item, only relevant for vertical data for which it must be provided
 
     Returns
     -------
@@ -165,6 +170,7 @@ def from_matched(
         z=z,
         x_item=x_item,
         y_item=y_item,
+        z_item=z_item,
         quantity=quantity,
     )
 
@@ -351,7 +357,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 +383,8 @@ def _match_space_time(
                 )
             case PointModelResult() as pmr, PointObservation():
                 aligned = pmr.align(observation, max_gap=max_model_gap)
+            case VerticalModelResult() as vmr, VerticalObservation():
+                aligned = vmr.align(observation)
             case _:
                 raise TypeError(
                     f"Matching not implemented for model type {type(mr)} and observation type {type(observation)}"
@@ -386,7 +396,6 @@ def _match_space_time(
             raise ValueError(
                 f"Aux variables are not allowed to have identical names. Choose either aux from obs or model. Overlapping: {overlapping}"
             )
-
         for dv in aligned:
             data[dv] = aligned[dv]
 

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",