Check quantity units in match()

docs/examples/Metocean_track_comparison.qmd

@@ -36,7 +36,7 @@ Csv files have no metadata on which quantity it contains, we add this manually,
 
 ```{python}
 o1 = ms.TrackObservation(df, item="surface_elevation", name='alti',
-                         quantity=ms.Quantity(name="Surface Elevation", unit="meter")) 
+                         quantity=ms.Quantity(name="Surface Elevation", unit="m")) 
 o1
 ```
 
@@ -67,7 +67,7 @@ mr
 df = pd.read_csv('../data/altimetry_NorthSea_20171027.csv',
                   index_col=0, parse_dates=True)
 o1 = ms.TrackObservation(df, item=2, name='alti',
-                         quantity=ms.Quantity(name="Surface Elevation", unit="meter"))
+                         quantity=ms.Quantity(name="Surface Elevation", unit="m"))
 o1
 ```
 
@@ -91,7 +91,7 @@ mr = ms.model_result('../data/NorthSeaHD_and_windspeed.dfsu',
 df = pd.read_csv('../data/altimetry_NorthSea_20171027.csv',
                  index_col=0, parse_dates=True)
 o1 = ms.TrackObservation(df, item=2, name='alti',
-                         quantity=ms.Quantity(name="Surface Elevation", unit="meter"))
+                         quantity=ms.Quantity(name="Surface Elevation", unit="m"))
 cmp = ms.match(o1, mr)
 cmp
 ```

src/modelskill/matching.py

@@ -204,6 +204,7 @@ def match(
     spatial_method: str | None = None,
     spatial_tolerance: float = 1e-3,
     obs_no_overlap: Literal["ignore", "error", "warn"] = "error",
+    check_quantity: Literal["error", "ignore"] = "ignore",
 ) -> Comparer: ...
 
 
@@ -216,6 +217,7 @@ def match(
     spatial_method: str | None = None,
     spatial_tolerance: float = 1e-3,
     obs_no_overlap: Literal["ignore", "error", "warn"] = "error",
+    check_quantity: Literal["error", "ignore"] = "ignore",
 ) -> ComparerCollection: ...
 
 
@@ -227,6 +229,7 @@ def match(
     spatial_method: str | None = None,
     spatial_tolerance: float = 1e-3,
     obs_no_overlap: Literal["ignore", "error", "warn"] = "error",
+    check_quantity: Literal["error", "ignore"] = "ignore",
 ):
     """Match observation and model result data in space and time
 
@@ -260,6 +263,10 @@ def match(
         and TrackObservation, by default 1e-3.
     obs_no_overlap: str, optional
         How to handle observations with no overlap with model results. One of: 'ignore', 'error', 'warn', by default 'error'.
+    check_quantity: str, optional
+        How to handle a model result whose unit differs from the observation's,
+        see [](`~modelskill.Quantity.is_compatible`). One of: 'error', 'ignore',
+        by default 'ignore'.
 
     Returns
     -------
@@ -273,6 +280,11 @@ def match(
     from_matched - Create a Comparer from observation and model results that are already matched
     """
 
+    if check_quantity not in ("error", "ignore"):
+        raise ValueError(
+            f"check_quantity must be 'error' or 'ignore', not {check_quantity!r}"
+        )
+
     if isinstance(obs, get_args(ObsInputType)):
         return _match_single_obs(
             obs,
@@ -281,6 +293,7 @@ def match(
             spatial_method=spatial_method,
             spatial_tolerance=spatial_tolerance,
             obs_no_overlap=obs_no_overlap,
+            check_quantity=check_quantity,
         )
 
     if isinstance(obs, Collection):
@@ -322,6 +335,7 @@ def match(
             spatial_method=spatial_method,
             spatial_tolerance=spatial_tolerance,
             obs_no_overlap=obs_no_overlap,
+            check_quantity=check_quantity,
         )
         for o in obs
     ]
@@ -339,6 +353,7 @@ def _match_single_obs(
     spatial_method: str | None,
     spatial_tolerance: float,
     obs_no_overlap: Literal["ignore", "error", "warn"],
+    check_quantity: Literal["error", "ignore"],
 ) -> Comparer | None:
     if isinstance(mod, get_args(MRInputType)):
         models: list = [mod]
@@ -364,6 +379,14 @@ def _match_single_obs(
         else:
             matching_obs = m
 
+        if check_quantity == "error" and not obs.quantity.is_compatible(
+            matching_obs.quantity
+        ):
+            raise ValueError(
+                f"Model result '{m.name}' has quantity {matching_obs.quantity}, "
+                f"which is not compatible with {obs.quantity} of observation "
+                f"'{obs.name}'."
+            )
         raw_mod_data[m.name] = matching_obs
 
     matched_data = _match_space_time(

src/modelskill/quantity.py

@@ -58,29 +58,33 @@ def __repr__(self):
             # hide is_directional if False to avoid clutter
             return f"Quantity(name='{self.name}', unit='{self.unit}')"
 
-    def is_compatible(self, other) -> bool:
+    def is_compatible(self, other: "Quantity") -> bool:
         """Check if the quantity is compatible with another quantity
 
+        Two quantities are compatible when their units agree. Names are not
+        compared, since the same physical quantity is often named differently
+        in observation and model, e.g. "Water Level" and "Surface Elevation".
+        Units are compared as written, so "meter" and "m" differ. An undefined
+        quantity, or one without a unit, is compatible with any other.
+
         Examples
         --------
         ```{python}
-        wl = Quantity(name="Water Level", unit="meter")
-        ws = Quantity(name="Wind Speed", unit="meter per second")
+        wl = Quantity(name="Water Level", unit="m")
+        ws = Quantity(name="Wind Speed", unit="m/s")
         wl.is_compatible(ws)
         ```
         ```{python}
-        uq = Quantity(name="Undefined", unit="Undefined")
-        wl.is_compatible(uq)
+        wl.is_compatible(Quantity(name="Surface Elevation", unit="m"))
+        ```
+        ```{python}
+        wl.is_compatible(Quantity.undefined())
         ```
         """
-
-        if self == other:
+        if _is_undefined_unit(self.unit) or _is_undefined_unit(other.unit):
             return True
 
-        if (self.name == "Undefined") or (other.name == "Undefined"):
-            return True
-
-        return False
+        return self.unit == other.unit
 
     @staticmethod
     def undefined() -> "Quantity":
@@ -185,7 +189,14 @@ def from_mikeio_eum_name(type_name: str) -> "Quantity":
                 raise ValueError(
                     f"{type_name=} is not recognized as a known type. Please create a Quantity(name='{type_name}' unit='<FILL IN UNIT>')"
                 )
-        unit = etype.units[0].name
+        unit = etype.units[0].short_name
         is_directional = unit == "degree"
         warnings.warn(f"{unit=} was automatically set for {type_name=}")
         return Quantity(name=type_name, unit=unit, is_directional=is_directional)
+
+
+def _is_undefined_unit(unit: str) -> bool:
+    # "" from Quantity.undefined() and from res1d/EPANET results, which carry no
+    # unit; "undefined" from mikeio items of EUM type Undefined; "Undefined" from
+    # earlier modelskill versions
+    return unit in ("", "undefined", "Undefined")

tests/test_match.py

@@ -7,6 +7,7 @@
 import modelskill as ms
 from modelskill.comparison._comparison import ItemSelection
 from modelskill.model.dfsu import DfsuModelResult
+
 try:
     from modelskill.network import _make_basic_network
 except ImportError:
@@ -1064,7 +1065,8 @@ def test_network_match_multi_obs_multi_model_comprehensive(
 def test_network_match_error_non_node_observation(network_mr, point_obs_error):
     """Test that non-NodeObservation raises appropriate error"""
     with pytest.raises(
-        TypeError, match="NetworkModelResult supports NodeObservation and ReachObservation"
+        TypeError,
+        match="NetworkModelResult supports NodeObservation and ReachObservation",
     ):
         ms.match(point_obs_error, network_mr)
 
@@ -1076,3 +1078,39 @@ def test_match_nodeobs_with_other_result(node_obs1, mr1):
         match="Extraction from .* to <class 'modelskill.obs.NodeObservation'> is not implemented.",
     ):
         ms.match(node_obs1, mr1)
+
+
+@pytest.fixture
+def wl_obs_discharge_model():
+    time = pd.date_range("2020-01-01", periods=5, freq="h")
+    obs = ms.PointObservation(
+        pd.Series(np.arange(5.0), index=time, name="obs"),
+        x=0.0,
+        y=0.0,
+        quantity=ms.Quantity(name="Water Level", unit="m"),
+    )
+    mod = ms.PointModelResult(
+        pd.Series(np.arange(5.0), index=time, name="mod"),
+        x=0.0,
+        y=0.0,
+        quantity=ms.Quantity(name="Discharge", unit="m^3/s"),
+    )
+    return obs, mod
+
+
+def test_match_raises_on_incompatible_quantity(wl_obs_discharge_model):
+    obs, mod = wl_obs_discharge_model
+    with pytest.raises(ValueError, match="not compatible"):
+        ms.match(obs, mod, check_quantity="error")
+
+
+def test_match_check_quantity_ignore(wl_obs_discharge_model):
+    obs, mod = wl_obs_discharge_model
+    cmp = ms.match(obs, mod, check_quantity="ignore")
+    assert cmp.n_points == 5
+
+
+def test_match_check_quantity_invalid_value(wl_obs_discharge_model):
+    obs, mod = wl_obs_discharge_model
+    with pytest.raises(ValueError, match="check_quantity"):
+        ms.match(obs, mod, check_quantity="warn")

tests/test_pointcompare.py

@@ -22,7 +22,7 @@ def drogden():
         item=0,
         x=355568.0,
         y=6156863.0,
-        quantity=ms.Quantity("Water Level", unit="meter"),
+        quantity=ms.Quantity("Water Level", unit="m"),
     )
 
 
@@ -149,7 +149,7 @@ def test_weighted_score(modelresult_oresund_WL):
         x=355568.0,
         y=6156863.0,
         quantity=ms.Quantity(
-            "Water Level", unit="meter"
+            "Water Level", unit="m"
         ),  # not sure if this is relevant in this test
     )
 
@@ -176,7 +176,7 @@ def test_weighted_score(modelresult_oresund_WL):
         x=355568.0,
         y=6156863.0,
         quantity=ms.Quantity(
-            "Water Level", unit="meter"
+            "Water Level", unit="m"
         ),  # not sure if this is relevant in this test
         weight=0.1,
     )

tests/test_quantity.py

@@ -12,14 +12,13 @@ def test_from_EUMType_string():
         # mikeio.EUMType.Significant_wave_height
         q = ms.Quantity.from_mikeio_eum_name("Significant_wave_height")
 
-    assert q.unit == "meter"
+    assert q.unit == "m"
 
     with pytest.warns(match="unit"):
         # mikeio.EUMType.Discharge
         q = ms.Quantity.from_mikeio_eum_name("Discharge")
 
-    # TODO should this be "meter^3 per second"? or m3/s?
-    assert q.unit == "meter_pow_3_per_sec"
+    assert q.unit == "m^3/s"
 
 
 def test_unknown_quantity_raises_error():
@@ -45,3 +44,33 @@ def test_from_cf_attrs_incomplete():
     q = ms.Quantity.from_cf_attrs({"long_name": "Wind speed"})
     assert q.name == ""
     assert q.unit == ""
+
+
+def test_is_compatible_ignores_name():
+    wl = ms.Quantity(name="Water Level", unit="m")
+    assert wl.is_compatible(ms.Quantity(name="Surface Elevation", unit="m"))
+
+
+def test_is_compatible_different_units():
+    wl = ms.Quantity(name="Water Level", unit="m")
+    assert not wl.is_compatible(ms.Quantity(name="Discharge", unit="m^3/s"))
+
+
+def test_is_compatible_compares_units_as_written():
+    wl = ms.Quantity(name="Water Level", unit="meter")
+    assert not wl.is_compatible(ms.Quantity(name="Water Level", unit="m"))
+
+
+@pytest.mark.parametrize(
+    "other",
+    [
+        ms.Quantity.undefined(),
+        ms.Quantity(name="Undefined", unit="Undefined"),
+        ms.Quantity(name="Undefined", unit="undefined"),
+        ms.Quantity(name="Pressure", unit=""),
+    ],
+)
+def test_is_compatible_without_unit(other):
+    q = ms.Quantity(name="Pressure", unit="MetresWater")
+    assert q.is_compatible(other)
+    assert other.is_compatible(q)