Vertical skill documentation and improved docstrings

docs/user-guide/data-structures.qmd

@@ -20,31 +20,37 @@ The primary data of ModelSkill are the data that needs to be compared: observati
 
 * `point`: 0D time series data
 * `track`: 0D time series data at moving locations (trajectories)
+* `vertical`: 1D depth profile time series at a fixed location
 * `grid`: gridded 2D data
 * `dfsu`: flexible mesh 2D data
 
-Point and track data are both `TimeSeries` objects, while grid and dfsu data are both `SpatialField` objects. `TimeSeries` objects are ready to be compared whereas data from `SpatialField` object needs to be *extracted* first (the extracted object will be of the `TimeSeries` type).
+Point, track and vertical data are all `TimeSeries` objects, while grid and dfsu data are both `SpatialField` objects. `TimeSeries` objects are ready to be compared whereas data from `SpatialField` object needs to be *extracted* first (the extracted object will be of the `TimeSeries` type).
 
 `TimeSeries` objects contains its data in an [](`xarray.Dataset`) with the actual data in the first DataArray and optional auxilliary data in the following DataArrays. The DataArrays have a `kind` attribute with either `observation` or `model`.
 
+Vertical data is stored in *long* format: a single `time` dimension in which each timestamp is repeated once per depth, with the depth stored in a non-index `z` coordinate.
+
 
 ## Comparer objects
 
 Comparer objects are results of a matching procedure (between observations and model results) or constructed directly from already matched data. A comparison of a *single* observation and one or more model results are stored in a [](`~modelskill.Comparer`) object. A comparison of *multiple* observations and one or more model results are stored in a [](`~modelskill.ComparerCollection`) object which is a collection of [](`~modelskill.Comparer`) objects.
 
-The matched data in a [](`~modelskill.Comparer`) is stored in an [](`xarray.Dataset`) which can be accessed via the `data` attribute. The Dataset has an attribute `gtype` which is a string describing the type of data (e.g. `point`, `track`). The first DataArray in the Dataset is the observation data, the next DataArrays are model result data and optionally additional DataArrays are auxilliarye data. Each of the DataArrays have a `kind` attribute with either `observation`, `model` or `aux`.
+The matched data in a [](`~modelskill.Comparer`) is stored in an [](`xarray.Dataset`) which can be accessed via the `data` attribute. The Dataset has an attribute `gtype` which is a string describing the type of data (e.g. `point`, `track`, `vertical`). The first DataArray in the Dataset is the observation data, the next DataArrays are model result data and optionally additional DataArrays are auxilliarye data. Each of the DataArrays have a `kind` attribute with either `observation`, `model` or `aux`.
 
 Both [](`~modelskill.Comparer`) and [](`~modelskill.ComparerCollection`) have a `plot` accessor for plotting the data (e.g. `cmp.plot.timeseries()` or `cmp.plot.scatter()`).
 
+A [](`~modelskill.Comparer`) with `gtype="vertical"` additionally has a `vertical` accessor (e.g. `cmp.vertical.skill()` or `cmp.vertical.plot.hovmoller()`) for depth-binned skill and profile plots.
+
 
 
 ## Skill objects
 
-Calling a skill method on a comparer object will return a skill object with skill scores (statistics) from comparing observation and model result data using different metrics (e.g. root mean square error). Two skill objects are currently implemented: [](`~modelskill.SkillTable`) and [](`~modelskill.skill_grid.SkillGrid`). The first is relevant for all ModelSkill users while the latter is relevant for users of the track data (e.g. MetOcean studies using satellite altimetry data).
+Calling a skill method on a comparer object will return a skill object with skill scores (statistics) from comparing observation and model result data using different metrics (e.g. root mean square error). Three skill objects are currently implemented: [](`~modelskill.SkillTable`), [](`~modelskill.skill_grid.SkillGrid`) and [](`~modelskill.skill_profile.SkillProfile`). The first is relevant for all ModelSkill users, while the second is relevant for users of track data (e.g. MetOcean studies using satellite altimetry data) and the third for users of vertical profile data.
 
 If `c` is a comparer object, then the following skill methods are available:
 
 * `c.skill()` -> `SkillTable`
 * `c.mean_skill()` -> `SkillTable`
 * `c.gridded_skill()` -> `SkillGrid`
+* `c.vertical.skill()` -> `SkillProfile` (only for vertical data)
 

docs/user-guide/matching.qmd

@@ -91,7 +91,7 @@ cc = ms.match([o1, o2], mr_dfsu, spatial_method='contained')
 
 ::: callout-note
 * Extraction of *track* data does not currently support the "contained" method.
-* Extraction of point data from 3D dfsu files is not yet fully supported. It is recommended to extract the data "offline" prior to using ModelSkill.
+* Extraction of *point* data from 3D dfsu files is not yet fully supported. It is recommended to extract the data "offline" prior to using ModelSkill. Extraction of whole *profiles* from 3D dfsu files is supported, see [matching of vertical profiles](#matching-of-vertical-profiles).
 :::
 
 
@@ -104,6 +104,43 @@ cc = ms.match([o1, o2], mr_nc, spatial_method='nearest')
 ```
 
 
+## Matching of vertical profiles
+
+A [](`~modelskill.VerticalObservation`) (e.g. a CTD cast) is matched with either a
+[](`~modelskill.VerticalModelResult`) - a column extracted beforehand - or a *3D*
+[](`~modelskill.DfsuModelResult`), in which case the column at the observation
+position is extracted as part of the matching.
+
+::: {.callout-note}
+Vertical profile support is still in development and the API may change.
+:::
+
+```{python}
+vo = ms.VerticalObservation("../data/vertical/VerticalProfile_obs1.dfs0",
+                            item="Salinity", z_item="z", x=657500, y=6553600)
+mr_3d = ms.model_result("../data/vertical/sigma_z_coast.dfsu", item="Salinity")
+cmp = ms.match(vo, mr_3d)
+cmp
+```
+
+Matching is done in two steps:
+
+1. The model profile is taken from the model time *closest* to each observation time. The tolerance defaults to half the median model time step, so observation times without a nearby model profile are dropped.
+2. The model values are interpolated *linearly in depth* onto the observation depths. Observation depths outside the model depth range are discarded - no extrapolation is performed.
+
+Note that the temporal step is a nearest-neighbour selection rather than the
+linear interpolation in time used for point and track data. This preserves each
+modelled profile as a whole instead of blending two profiles.
+
+The result is a [](`~modelskill.Comparer`) with `gtype="vertical"` in *long* format:
+one row per (time, depth) pair.
+
+::: callout-note
+* Extraction from a 3D dfsu file only supports `spatial_method="contained"` (the default); the other methods are not implemented for vertical profiles.
+* The `max_model_gap` argument is not yet supported for vertical matching.
+:::
+
+
 ## Event-based matching and handling of gaps
 
 If the model result data contains gaps either because only events are stored or because of missing data, the `max_model_gap` argument of the [](`~modelskill.match`) function can be used to specify the maximum allowed gap (in seconds) in the model result data. This will avoid interpolating model data over long gaps in the model result data!

docs/user-guide/overview.qmd

@@ -18,7 +18,7 @@ they will need to be defined and then matched in space and time with the `match(
 
 ### Define observations
 
-The observations can be defined as either a [`PointObservation`](`modelskill.PointObservation`) or a [`TrackObservation`](`modelskill.TrackObservation`) (a moving point).
+The observations can be defined as a [`PointObservation`](`modelskill.PointObservation`), a [`TrackObservation`](`modelskill.TrackObservation`) (a moving point) or a [`VerticalObservation`](`modelskill.VerticalObservation`) (a depth profile at a fixed position).
 
 ```{python}
 import modelskill as ms
@@ -34,8 +34,9 @@ A model result will either be a simple point/track like the observations, or spa
 
 * [`PointModelResult`](`modelskill.PointModelResult`) - a point result from a dfs0/nc file or a DataFrame
 * [`TrackModelResult`](`modelskill.TrackModelResult`) - a track result from a dfs0/nc file or a DataFrame
+* [`VerticalModelResult`](`modelskill.VerticalModelResult`) - a depth profile from a dfs0 file or a DataFrame
 * [`GridModelResult`](`modelskill.GridModelResult`) - a spatial field from a dfs2/nc file or a Xarray Dataset
-* [`DfsuModelResult`](`modelskill.DfsuModelResult`) - a spatial field from a dfsu file
+* [`DfsuModelResult`](`modelskill.DfsuModelResult`) - a spatial field from a dfsu file (2D or 3D)
 
 ```{python}
 mr1 = ms.PointModelResult("../data/model.dfs0", item="WL")
@@ -51,6 +52,24 @@ cc
 ```
 
 
+## Vertical profiles
+
+Model results can also be validated against vertical profiles, e.g. CTD casts of
+salinity or temperature. The observation is a
+[`VerticalObservation`](`modelskill.VerticalObservation`) and the model result is either a
+[`VerticalModelResult`](`modelskill.VerticalModelResult`) (a pre-extracted column) or a
+3D [`DfsuModelResult`](`modelskill.DfsuModelResult`), from which the column at the
+observation position is *extracted* automatically (see [matching](matching.qmd)).
+
+The result is an ordinary [`Comparer`](`modelskill.Comparer`) with `gtype="vertical"`,
+so `cmp.skill()` and all the usual plots work as before. In addition, a `vertical`
+accessor gives access to depth-aware skill and plots (see [plotting](plotting.qmd)).
+The comparer can be sliced by depth with `cmp.sel(z=-2)` or `cmp.sel(z=slice(-10, 0))`
+(see [selecting data](selecting-data.qmd)), and the profile can be collapsed to a
+single value per timestep with `cmp.vertical.mean()`, `.min()` or `.max()` (which
+return point comparers).
+
+
 ## Analysis
 
 Once the observations and model results are matched, the `Comparer` object can be used for analysis and plotting. 

docs/user-guide/plotting.qmd

@@ -72,6 +72,55 @@ cmp.plot.scatter();
 ```
 
 
+## Vertical profiles
+
+A [](`~modelskill.Comparer`) holding vertical profile data (`gtype="vertical"`)
+has an extra `vertical` accessor with depth-aware plots, available through
+[`cmp.vertical.plot`](`modelskill.comparison.VerticalPlotter`).
+
+The `profile` plot shows model and observed values against depth. It is meant
+for a single point in time, so select a time first:
+
+```{python}
+vo = ms.VerticalObservation("../data/vertical/VerticalProfile_obs1.dfs0",
+                            item="Salinity", z_item="z", x=657500, y=6553600)
+mr3d = ms.DfsuModelResult("../data/vertical/sigma_z_coast.dfsu", item="Salinity")
+vcmp = ms.match(vo, mr3d)
+
+vcmp.sel(time="2022-06-15").vertical.plot.profile();
+```
+
+The model is drawn as a line (the full modelled column) and the observations as
+markers. Pass `show_matched_model=True` to also mark the model values
+interpolated to the observation depths. `vertical.plot()` is a shorthand for
+`vertical.plot.profile()`.
+
+The `hovmoller` plot shows the whole comparison as a function of depth and time,
+with the model as a filled contour and the observations overlaid as markers:
+
+```{python}
+vcmp.vertical.plot.hovmoller(figsize=(12, 4));
+```
+
+With more than one model result, the model to contour must be named explicitly,
+e.g. `vcmp.vertical.plot.hovmoller(model="m1")`.
+
+Skill computed in depth bins with
+[`vertical.skill()`](`modelskill.comparison.VerticalAccessor.skill`) returns a
+[](`~modelskill.skill_profile.SkillProfile`), whose metrics can be plotted as a
+horizontal bar plot per depth bin:
+
+```{python}
+sk = vcmp.vertical.skill(bins=4)
+sk.rmse.plot();
+```
+
+The bins can also be given explicitly, e.g. `bins=[-20, -15, -10, -5, 0]`.
+
+Note that all the ordinary Comparer plots shown above work on vertical data too;
+they simply ignore the depth dimension.
+
+
 ## Taylor diagrams
 
 A Taylor diagram shows how well a model result matches an observation in terms of correlation, standard deviation and root mean squared error. The `taylor` plot can be accessed through the Comparer [`plot`](`modelskill.comparison.ComparerPlotter`) accessor or the ComparerCollection [`plot`](`modelskill.comparison.ComparerCollectionPlotter`) accessor:

docs/user-guide/selecting-data.qmd

@@ -53,13 +53,14 @@ cmp = ms.match(o, [m1, m2])
 
 The [](`~modelskill.Comparer.sel`) method allows you to select data based on specific criteria such as time, model name, or spatial area. It returns a new `Comparer` object with the selected data. This method is highly versatile and supports multiple selection parameters, which can be combined.
 
-**Syntax:** `Comparer.sel(model=None, time=None, area=None)`
+**Syntax:** `Comparer.sel(model=None, time=None, area=None, z=None)`
 
 | Parameter    | Type                        | Description                                             | Default |
 |--------------|-----------------------------|---------------------------------------------------------|---------|
 | `model`      | str, int, or list           | Model name or index. Selects specific models.           | None    |
 | `time`       | str, datetime, or slice     | Specific time or range for selection.                  | None    |
 | `area`       | list of float or Polygon    | Bounding box [x0, y0, x1, y1] or a polygon area filter. | None    |
+| `z`          | float or slice              | Depth or depth range (vertical data only).              | None    |
 
 **Example 1: Selecting data by time**
 ```{python}
@@ -84,6 +85,43 @@ cmp_area = cmp.sel(area=[4.0, 52.5, 5.0, 53.0])
 
 This filters the data within the bounding box defined by `[x0, y0, x1, y1]`.
 
+**Example 4: Selecting by depth**
+
+For vertical profile data (`gtype="vertical"`), `sel()` additionally accepts a
+`z` argument, given either as a single depth or as a slice:
+
+```{python}
+#| code-fold: true
+#| code-summary: "Vertical observation and model data"
+vo = ms.VerticalObservation("../data/vertical/VerticalProfile_obs1.dfs0",
+                            item="Salinity", z_item="z", x=657500, y=6553600)
+mr3d = ms.DfsuModelResult("../data/vertical/sigma_z_coast.dfsu", item="Salinity")
+vcmp = ms.match(vo, mr3d)
+```
+
+```{python}
+vcmp.sel(z=slice(-10, 0))
+```
+
+A slice keeps all observation depths within the range. A single value instead
+selects the nearest depth:
+
+```{python}
+vcmp.sel(z=-2)
+```
+
+Selections can be combined and given a descriptive name, which is convenient
+when comparing skill between depth ranges:
+
+```{python}
+cmp_surface = vcmp.sel(z=slice(-10, 0)).rename({vcmp.name: "surface"})
+cmp_deep = vcmp.sel(z=slice(None, -10)).rename({vcmp.name: "deep"})
+ms.ComparerCollection([cmp_surface, cmp_deep]).skill()
+```
+
+`z` also works on a [](`~modelskill.ComparerCollection`), and depth can equally
+be filtered with `query()`, e.g. `vcmp.query("z > -10")`.
+
 ### `where()` method
 
 The [](`~modelskill.Comparer.where`) method is used to filter data conditionally. It works similarly to `xarray`'s `where` method and returns a new `Comparer` object with values satisfying a given condition. Other values will be masked (set to `NaN`).
@@ -94,14 +132,14 @@ The [](`~modelskill.Comparer.where`) method is used to filter data conditionally
 |-----------|------------------------------|-----------------------------------------------|
 | `cond`    | bool, np.ndarray, or xr.DataArray | Condition to filter values (True or False). |
 
-**Example 4: Filtering data conditionally**
+**Example 5: Filtering data conditionally**
 ```{python}
 cmp.where(cmp.data.Observation > 3)
 ```
 
 This filters out any rows where the observation values are not greater than 3.
 
-**Example 5: Multiple conditions**
+**Example 6: Multiple conditions**
 ```{python}
 cmp.where((cmp.data.m1 < 2.9) & (cmp.data.Observation > 3))
 ```
@@ -120,7 +158,7 @@ The [](`~modelskill.Comparer.query`) method uses a [](`pandas.DataFrame.query`)-
 |-----------|--------|---------------------------------------------|
 | `query`   | str    | Query string for filtering data.           |
 
-**Example 6: Querying data**
+**Example 7: Querying data**
 ```{python}
 cmp.query("Observation > 3.0 and m1 < 2.9")
 ```
@@ -143,22 +181,22 @@ sk = cmp.skill(metrics=["rmse", "mae", "si"])
 sk
 ```
 
-**Example 7: Select model**
+**Example 8: Select model**
 ```{python}
 sk.sel(model='m1')
 ```
 
 Here, `sk` contains skill scores for all models, and `sk_m1` filters the results to include only model "m1". Observations can be selected in the same way.
 
-**Example 8: Querying skill scores**
+**Example 9: Querying skill scores**
 ```{python}
 sk_high_rmse = sk.query("rmse > 0.3")
 sk_high_rmse
 ```
 
 This filters the `SkillTable` to include only rows where the root mean square error (RMSE) exceeds 0.3.
 
-**Example 9: Accessing and visualizing specific metrics**
+**Example 10: Accessing and visualizing specific metrics**
 ```{python}
 sk_rmse = sk.rmse
 sk_rmse

src/modelskill/comparison/_vertical_comparison.py

@@ -23,7 +23,7 @@ def __init__(self, comparer):
         self.comparer = comparer
 
     def __call__(self, *args, **kwargs) -> matplotlib.axes.Axes:
-        """Plot scatter plot of modelled vs observed data"""
+        """Plot vertical profiles of model and observations (alias for `profile`)"""
         return self.profile(*args, **kwargs)
 
     def profile(
@@ -43,6 +43,9 @@ def profile(
             Matplotlib Axes to plot on (if None, a new figure and axes will be created).
         figsize : tuple, optional
             Size of the figure (only used if ax is None).
+        show_matched_model : bool, optional
+            Also show the model values interpolated to the observation depths
+            as markers, by default False (only the raw model profile is drawn).
 
         Returns
         -------
@@ -284,15 +287,27 @@ def _agg(self, agg_func: str = "mean"):
         return Comparer(ds, raw_mod_data=raw_mod_data)
 
     def mean(self):
-        """Aggregate the comparison vertically using a specified aggregation function."""
+        """Aggregate the profile to a single mean value per timestep.
+
+        Returns a point `Comparer`. The raw model profile is restricted to the
+        depth range covered by the observation before aggregating.
+        """
         return self._agg(agg_func="mean")
 
     def min(self):
-        """Aggregate the comparison vertically using a specified aggregation function."""
+        """Aggregate the profile to a single min value per timestep.
+
+        Returns a point `Comparer`. The raw model profile is restricted to the
+        depth range covered by the observation before aggregating.
+        """
         return self._agg(agg_func="min")
 
     def max(self):
-        """Aggregate the comparison vertically using a specified aggregation function."""
+        """Aggregate the profile to a single max value per timestep.
+
+        Returns a point `Comparer`. The raw model profile is restricted to the
+        depth range covered by the observation before aggregating.
+        """
         return self._agg(agg_func="max")
 
     def skill(