Explicit context information for unit conversion (#1227)

### What kind of change does this PR introduce?

* Disables hydro context
* Enable hydro context for specific operations
* Add `infer_context` function to encapsulate the logic determining
context when not explicit (e.g. generic indices, datachecks, etc.)
* Adjust tests

### Does this PR introduce a breaking change?
Things that worked out of the box by assuming the hydro context applied
to every xclim operation will break. A bit more care needs to be taken
when writing indices, especially generic ones, to handle unit conversion
correctly.

HISTORY.rst

@@ -12,10 +12,15 @@ New features and enhancements
 * ``xclim.core.units.convert_units_to`` can now perform automatic conversions based on the standard name of the input when needed. (:issue:`1205`, :pull:`1206`).
     - Conversion from amount (thickness) to flux (rate), using ``amount2rate`` and ``rate2amount``.
     - Conversion from amount to thickness for liquid water quantities, using the new ``amount2lwethickness`` and ``lwethickness2amount``. This is similar to the implicit transformations enabled by the "hydro" unit context.
+* ``xclim.core.units.convert_units_to`` can now perform automatic conversions based on the standard name of the input when needed. (:issue:`1205`, :pull:`1206`).
+    - Conversion from amount (thickness) to flux (rate), using ``amount2rate`` and ``rate2amount``.
+    - Conversion from amount to thickness for liquid water quantities, using the new ``amount2lwethickness`` and ``lwethickness2amount``. This is similar to the implicit transformations enabled by the "hydro" unit context.
+    - Passing ``context='infer'`` will activate the "hydro" context if the source or the target are DataArrays with a standard name that is compatible, as decided by the new ``xclim.core.units.infer_context`` function.
 
 Breaking changes
 ^^^^^^^^^^^^^^^^
 * Rewrite of ``xclim.core.calendar.time_bnds``. It should now be more resilient and versatile, but all ``cftime_*`` and ``cfindex_*`` functions were removed. (:issue:`74`, :pull:`1207`).
+* `hydro` context is not always enabled, as it led to unwanted unit conversions. Unit conversion operations now need to explicitly declare the `hydro` context to support conversions from `kg / m2 /s` to `mm/day`. (:issue:`1208`, :pull:`1227`).
 * Many previously deprecated indices and indicators have been removed from `xclim` (:pull:`1228`), with replacement indices/indicators suggested as follows:
     - ``xclim.indicators.atmos.fire_weather_indexes`` → ``xclim.indicators.atmos.cffwis_indices``
     - ``xclim.indices.freshet_start`` → ``xclim.indices.first_day_temperature_above``
@@ -38,6 +43,7 @@ Bug fixes
 Internal changes
 ^^^^^^^^^^^^^^^^
 * Minor adjustments to GitHub Actions workflows (newest Ubuntu images, updated actions version, better CI triggering). (:pull:`1221`).
+* Pint units `context` added to various operations, tests and `Indicator` attributes. (:issue:`1208`, :pull:`1227`).
 * Updated article from Alavoine & Grenier (2022) within documentation. Many article reference URLs have been updated to use HTTPS where possible. (:issue:`1246`, :pull:`1247`).
 * Added relevant variable dataflag checks for potential evaporation, convective precipitation, and air pressure at sea level. (:pull:`1241`).
 * Documentation restructured to include `ReadMe` page (as `About`) with some minor changes to documentation titles. (:pull:`1233`).

docs/notebooks/example.yml

@@ -16,6 +16,7 @@ indicators:
     base: rx1day
     cf_attrs:
       long_name: Highest 1-day precipitation amount
+    context: hydro
   RX5day_canopy:
     base: max_n_day_precipitation_amount
     cf_attrs:
@@ -25,12 +26,14 @@ indicators:
     parameters:
       freq: QS-DEC
       window: 5
+    context: hydro
   R75pdays:
     base: days_over_precip_thresh
     parameters:
       pr_per:
         description: Daily 75th percentile of wet day precipitation flux.
       thresh: 1 mm/day
+    context: hydro
   fd:
     compute: count_occurrences
     input:
@@ -60,10 +63,12 @@ indicators:
     parameters:
       perc: 95
     references: climdex
+    context: hydro
   R99p:
     base: .R95p
     cf_attrs:
       - var_name: R99p
       - var_name: R99p_days
     parameters:
       perc: 99
+    context: hydro

docs/notebooks/sdba.ipynb

@@ -452,6 +452,12 @@
     "dref = open_dataset(\n",
     "    \"sdba/ahccd_1950-2013.nc\", chunks={\"location\": 1}, drop_variables=[\"lat\", \"lon\"]\n",
     ").sel(time=slice(\"1981\", \"2010\"))\n",
+    "\n",
+    "# Fix the standard name of the `pr` variable.\n",
+    "# This allows the convert_units_to below to infer the correct CF transformation (precip rate to flux)\n",
+    "# see the \"Unit handling\" notebook\n",
+    "dref.pr.attrs[\"standard_name\"] = \"lwe_precipitation_rate\"\n",
+    "\n",
     "dref = dref.assign(\n",
     "    tasmax=convert_units_to(dref.tasmax, \"K\"),\n",
     "    pr=convert_units_to(dref.pr, \"kg m-2 s-1\"),\n",

setup.cfg

@@ -3,13 +3,13 @@ current_version = 0.39.12-beta
 commit = True
 tag = False
 parse = (?P<major>\d+)\.(?P<minor>\d+).(?P<patch>\d+)(\-(?P<release>[a-z]+))?
-serialize = 
+serialize =
 	{major}.{minor}.{patch}-{release}
 	{major}.{minor}.{patch}
 
 [bumpversion:part:release]
 optional_value = gamma
-values = 
+values =
 	beta
 	gamma
 
@@ -26,14 +26,14 @@ relative_files = True
 omit = */tests/*.py
 
 [flake8]
-exclude = 
+exclude =
 	.git,
 	docs,
 	build,
 	.eggs,
 max-line-length = 88
 max-complexity = 12
-ignore = 
+ignore =
 	C901
 	E203
 	E231
@@ -43,12 +43,12 @@ ignore =
 	F403
 	W503
 	W504
-per-file-ignores = 
+per-file-ignores =
 	tests/*:E402
-rst-directives = 
+rst-directives =
 	bibliography
 	autolink-skip
-rst-roles = 
+rst-roles =
 	doc,
 	mod,
 	py:attr,
@@ -67,7 +67,7 @@ rst-roles =
 	cite:p
 	cite:t
 	cite:ts
-extend-ignore = 
+extend-ignore =
 	RST399,
 	RST201,
 	RST203,
@@ -87,19 +87,19 @@ addopts =
 	--maxprocesses=6
 	--dist=loadscope
 norecursedirs = docs/notebooks/*
-filterwarnings = 
+filterwarnings =
 	ignore::UserWarning
 testpaths = xclim/testing/tests xclim/testing/tests/test_sdba
 usefixtures = xdoctest_namespace
 doctest_optionflags = NORMALIZE_WHITESPACE IGNORE_EXCEPTION_DETAIL NUMBER ELLIPSIS
-markers = 
+markers =
 	slow: marks tests as slow (deselect with '-m "not slow"')
 	requires_docs: mark tests that can only be run with documentation present
 
 [pycodestyle]
 count = False
 exclude = xclim/testing/tests
-ignore = 
+ignore =
 	E226,
 	E402,
 	E501,

tox.ini

@@ -36,9 +36,9 @@ commands =
 [testenv:docs]
 description = Build the documentation with makefile under {basepython}
 setenv =
+    !notebooks: SKIP_NOTEBOOKS = 1
     PYTHONPATH = {toxinidir}
     READTHEDOCS = 1
-    !notebooks: SKIP_NOTEBOOKS = 1
 commands =
     make docs
 allowlist_externals =
@@ -88,7 +88,10 @@ setenv =
     COV_CORE_SOURCE =
     PYTEST_ADDOPTS = "--color=yes"
     PYTHONPATH = {toxinidir}
-passenv = CI GITHUB_* LD_LIBRARY_PATH
+passenv =
+    CI
+    GITHUB_*
+    LD_LIBRARY_PATH
 extras = dev
 deps =
     coverage: coveralls

xclim/core/dataflags.py

@@ -19,7 +19,7 @@
 from ..indices.run_length import suspicious_run
 from .calendar import climatological_mean_doy, within_bnds_doy
 from .formatting import update_xclim_history
-from .units import convert_units_to, declare_units, str2pint
+from .units import convert_units_to, declare_units, infer_context, str2pint, units
 from .utils import (
     VARIABLES,
     InputKind,
@@ -322,7 +322,7 @@ def very_large_precipitation_events(
     >>> rate = "300 mm d-1"
     >>> flagged = very_large_precipitation_events(ds.pr, thresh=rate)
     """
-    thresh_converted = convert_units_to(thresh, da)
+    thresh_converted = convert_units_to(thresh, da, context="hydro")
     very_large_events = _sanitize_attrs(da > thresh_converted)
     description = f"Precipitation events in excess of {thresh} for {da.name}."
     very_large_events.attrs["description"] = description
@@ -365,7 +365,9 @@ def values_op_thresh_repeating_for_n_or_more_days(
     ...     ds.pr, n=days, thresh=units, op=comparison
     ... )
     """
-    thresh = convert_units_to(thresh, da)
+    thresh = convert_units_to(
+        thresh, da, context=infer_context(standard_name=da.attrs.get("standard_name"))
+    )
 
     repetitions = _sanitize_attrs(suspicious_run(da, window=n, op=op, thresh=thresh))
     description = (

xclim/core/units.py

@@ -30,6 +30,7 @@
     "check_units",
     "convert_units_to",
     "declare_units",
+    "infer_context",
     "infer_sampling_units",
     "lwethickness2amount",
     "pint_multiply",
@@ -104,7 +105,6 @@
     lambda ureg, x: x * (1000 * ureg.kg / ureg.m**3),
 )
 units.add_context(hydro)
-units.enable_contexts(hydro)
 
 
 CF_CONVERSIONS = safe_load(open_text("xclim.data", "variables.yml"))["conversions"]
@@ -293,6 +293,9 @@ def convert_units_to(
 ) -> Convertible:
     """Convert a mathematical expression into a value with the same units as a DataArray.
 
+    If the dimensionalities of source and target units differ, automatic CF conversions
+    will be applied when possible. See :py:func:`xclim.core.units.cf_conversion`.
+
     Parameters
     ----------
     source : str or xr.DataArray or units.Unit or units.Quantity
@@ -301,6 +304,9 @@ def convert_units_to(
         Target array of values to which units must conform.
     context : str, optional
         The unit definition context. Default: None.
+        If "infer", it will be inferred with :py:func:`xclim.core.units.infer_context` using
+        the standard name from the `source` or, if none is found, from the `target`.
+        This means that the 'hydro' context could be activated if any one of the standard names allows it.
 
     Returns
     -------
@@ -309,7 +315,17 @@ def convert_units_to(
         The outputted type is always similar to `source` initial type.
         Attributes are preserved unless an automatic CF conversion is performed,
         in which case only the new `standard_name` appears in the result.
+
+    See Also
+    --------
+    cf_conversion
+    amount2rate
+    rate2amount
+    amount2lwethickness
+    lwethickness2amount
     """
+    context = context or "none"
+
     # Target units
     if isinstance(target, units.Unit):
         target_unit = target
@@ -320,13 +336,25 @@ def convert_units_to(
             "target must be either a pint Unit or a xarray DataArray."
         )
 
+    if context == "infer":
+        ctxs = []
+        if isinstance(source, xr.DataArray):
+            ctxs.append(infer_context(source.attrs.get("standard_name")))
+        if isinstance(target, xr.DataArray):
+            ctxs.append(infer_context(target.attrs.get("standard_name")))
+        # If any one of the target or source is compatible with the "hydro" context, use it.
+        if "hydro" in ctxs:
+            context = "hydro"
+        else:
+            context = "none"
+
     if isinstance(source, str):
         q = str2pint(source)
         # Return magnitude of converted quantity. This is going to fail if units are not compatible.
-        return q.to(target_unit).m
+        return q.to(target_unit, context).m
 
     if isinstance(source, units.Quantity):
-        return source.to(target_unit).m
+        return source.to(target_unit, context).m
 
     if isinstance(source, xr.DataArray):
         source_unit = units2pint(source)
@@ -856,6 +884,8 @@ def check_units(val: str | int | float | None, dim: str | None) -> None:
     if dim is None or val is None:
         return
 
+    context = infer_context(dimension=dim)
+
     if str(val).startswith("UNSET "):
         warnings.warn(
             "This index calculation will soon require user-specified thresholds.",
@@ -886,11 +916,12 @@ def check_units(val: str | int | float | None, dim: str | None) -> None:
         return
 
     # Check if there is a transformation available
-    start = pint.util.to_units_container(val_dim)
-    end = pint.util.to_units_container(expected)
-    graph = units._active_ctx.graph  # noqa
-    if pint.util.find_shortest_path(graph, start, end):
-        return
+    with units.context(context):
+        start = pint.util.to_units_container(val_dim)
+        end = pint.util.to_units_container(expected)
+        graph = units._active_ctx.graph  # noqa
+        if pint.util.find_shortest_path(graph, start, end):
+            return
 
     raise ValidationError(
         f"Data units {val_units} are not compatible with requested {dim}."
@@ -986,3 +1017,42 @@ def ensure_delta(unit: str = None):
     if "degree_Rankine" in u._units:
         delta_unit = pint2cfunits(u / units2pint("°R") * units2pint("delta_degF"))
     return delta_unit
+
+
+def infer_context(standard_name=None, dimension=None):
+    """Return units context based on either the variable's standard name or the pint dimension.
+
+    Valid standard names for the hydro context are those including the terms "rainfall", "lwe" (liquid water equivalent) and
+    "precipitation". The latter is technically incorrect, as any phase of precipitation could be referenced.
+    Standard names for evapotranspiration, evaporation and canopy water amounts are also associated with the hydro context.
+
+    Parameters
+    ----------
+    standard_name: str
+      CF-Convention standard name.
+    dimension: str
+      Pint dimension, e.g. '[time]'.
+
+    Returns
+    -------
+    str
+      "hydro" if variable is a liquid water flux, otherwise "none"
+    """
+    csn = (
+        (
+            standard_name
+            in [
+                "water_potential_evapotranspiration_flux",
+                "canopy_water_amount",
+                "water_evaporation_amount",
+            ]
+            or "rainfall" in standard_name
+            or "lwe" in standard_name
+            or "precipitation" in standard_name
+        )
+        if standard_name is not None
+        else False
+    )
+    cdim = (dimension == "[precipitation]") if dimension is not None else False
+
+    return "hydro" if csn or cdim else "none"

xclim/data/anuclim.yml

@@ -43,6 +43,7 @@ indicators:
       standard_name: lwe_thickness_of_precipitation_amount
       cell_methods: "time: sum"
       units: mm
+    context: hydro
   P13_PrecipWettestPeriod:
     allowed_periods: [A]
     src_freq: ['D', '7D', 'M']
@@ -53,6 +54,7 @@ indicators:
       units: mm
     parameters:
       op: wettest
+    context: hydro
   P14_PrecipDriestPeriod:
     allowed_periods: [A]
     src_freq: ['D', '7D', 'M']
@@ -63,6 +65,7 @@ indicators:
       units: mm
     parameters:
       op: driest
+    context: hydro
   P15_PrecipSeasonality:
     allowed_periods: [A]
     src_freq: ['D', '7D', 'M']

xclim/data/cf.yml

@@ -36,6 +36,7 @@ indicators:
       reducer: max
       threshold: 1 mm day-1
     references: ETCCDI
+    context: hydro
   cddcoldTT:
     cf_attrs:
     - cell_methods: 'time: sum over days'
@@ -346,6 +347,7 @@ indicators:
       reducer: max
       threshold: 1 mm day-1
     references: ETCCDI
+    context: hydro
   ddgtTT:
     cf_attrs:
     - cell_methods: 'time: sum over days'