Dimensionless Units

data/dimensionless_mappings.yaml

@@ -0,0 +1,5 @@
+# In general:
+# model_variable_name:
+#   cmor_unit_string: pint_friendly_SI_units
+so:
+  "0.001": g/kg

examples/sample.yaml

@@ -14,6 +14,7 @@ pymorize:
   fixed_jobs: 12
   # minimum_jobs: 8
   # maximum_jobs: 30
+  dimensionless_mapping_table: ../data/dimensionless_mappings.yaml
 rules:
   - name: paul_example_rule
     description: "You can put some text here"

src/pymorize/cmorizer.py

@@ -2,6 +2,7 @@
 
 import dask
 import pandas as pd
+import xarray as xr
 import questionary
 import yaml
 from dask.distributed import Client
@@ -12,9 +13,14 @@
 from prefect_dask import DaskTaskRunner
 from rich.progress import track
 
-from .data_request import (DataRequest, DataRequestTable, DataRequestVariable,
-                           IgnoreTableFiles)
+from .data_request import (
+    DataRequest,
+    DataRequestTable,
+    DataRequestVariable,
+    IgnoreTableFiles,
+)
 from .filecache import fc
+from .units import handle_unit_conversion
 from .logging import logger
 from .pipeline import Pipeline
 from .rule import Rule
@@ -52,6 +58,7 @@ def __init__(
         self._post_init_create_data_request()
         self._post_init_populate_rules_with_tables()
         self._post_init_data_request_variables()
+        self._post_init_read_dimensionless_unit_mappings()
 
     def _post_init_configure_dask(self):
         """
@@ -148,6 +155,36 @@ def _post_init_data_request_variables(self):
         self._rules_expand_drvs()
         self._rules_depluralize_drvs()
 
+    def _post_init_read_dimensionless_unit_mappings(self):
+        """
+        Reads the dimensionless unit mappings from a configuration file and
+        updates the rules with these mappings.
+
+        This method reads the dimensionless unit mappings from a file specified
+        in the configuration. If the file is not specified or does not exist,
+        an empty dictionary is used. The mappings are then added to each rule
+        in the `rules` attribute.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        None
+        """
+        pymorize_cfg = self._pymorize_cfg
+        unit_map_file = pymorize_cfg.get("dimensionless_mapping_table", None)
+        if unit_map_file is None:
+            logger.warning("No dimensionless unit mappings file specified!")
+            dimensionless_unit_mappings = {}
+        else:
+            with open(unit_map_file, "r") as f:
+                dimensionless_unit_mappings = yaml.safe_load(f)
+        # Add to rules:
+        for rule in self.rules:
+            rule.dimensionless_unit_mappings = dimensionless_unit_mappings
+
     def find_matching_rule(
         self, data_request_variable: DataRequestVariable
     ) -> Rule or None:
@@ -225,6 +262,7 @@ def validate(self):
         # self._check_rules_for_table()
         # self._check_rules_for_output_dir()
         self._check_is_subperiod()
+        self._check_units()
 
     def _check_is_subperiod(self):
         logger.info("checking frequency in netcdf file and in table...")
@@ -257,6 +295,25 @@ def _check_is_subperiod(self):
         if errors:
             for err in errors:
                 logger.error(err)
+            raise errors[0]
+
+    def _check_units(self):
+        errors = []
+        for rule in self.rules:
+            for input_collection in rule.inputs:
+                try:
+                    filename = input_collection.files[0]
+                except IndexError:
+                    break
+                da = xr.open_dataarray(filename, use_cftime=True)
+                try:
+                    handle_unit_conversion(da, rule)
+                except ValueError as e:
+                    msg = f"Error while converting units for {filename}: {e}"
+                    logger.error(msg)
+                    errors.append(e)
+        if errors:
+            raise errors[0]
 
     @classmethod
     def from_dict(cls, data):
@@ -285,6 +342,7 @@ def from_dict(cls, data):
         instance._post_init_populate_rules_with_tables()
         instance._post_init_create_data_request()
         instance._post_init_data_request_variables()
+        instance._post_init_read_dimensionless_unit_mappings()
         return instance
 
     def add_rule(self, rule):

src/pymorize/units.py

@@ -5,15 +5,10 @@
 In case the units in model files differ from CMIP Tables, this module attempts to
 convert them automatically.
 
-In case of missing units in either model files or CMIP Tables, this module can
-not convert from a dimentionless base to something with dimension. Dealing with
-such thing have to done with `action` section in the Rules module on a per
-variable basis.
-
-Additionally, the cmip frequencies are mapped here. The CMIP6 frequency
-names and corresponding number of days are available as a dictionary in the
-``CMIP_FREQUENCIES`` variable. Assignment of these frequencies to the unit registry
-can be done with the ``assign_frequency_to_unit_registry`` function.
+Conversion to-or-from a dimensionless quantity is ambiguous. In this case,
+provide a mapping of what this dimensionless quantity represents and that
+is used for the conversion. `data/dimensionless_mappings.yaml` contains some
+examples on how the mapping is written.
 """
 
 import re
@@ -29,17 +24,12 @@
 import xarray as xr
 from chemicals import periodic_table
 
-from .frequency import CMIP_FREQUENCIES
 from .logging import logger
 from .rule import Rule
 
 ureg = pint_xarray.unit_registry
 
-
-def assign_frequency_to_unit_registry():
-    """Assign the CMIP6 frequencies to the unit registry."""
-    for freq_name, days in CMIP_FREQUENCIES.items():
-        ureg.define(f"{freq_name} = {days} * d")
+ureg.define("practical_salinity_unit = g/kg = psu = PSU")
 
 
 def handle_chemicals(
@@ -112,15 +102,28 @@ def handle_unit_conversion(da: xr.DataArray, rule: Rule) -> xr.DataArray:
     model_unit = rule.get("model_unit")
     from_unit = da.attrs.get("units")
     if model_unit is not None:
-        logger.debug(
-            f"using user defined unit ({model_unit}) instead of ({from_unit}) from DataArray "
+        logger.info(
+            f"using user defined unit ({model_unit}) instead of ({from_unit}) from the original file"
         )
         from_unit = model_unit
     handle_chemicals(from_unit)
     handle_chemicals(to_unit)
     new_da = da.pint.quantify(from_unit)
-    logger.debug(f"Converting units: {from_unit} -> {to_unit}")
-    new_da = new_da.pint.to(to_unit).pint.dequantify()
+    dimless = rule.get("dimensionless_unit_mappings", {})
+    cmor_variable = rule.get("cmor_variable")
+    if cmor_variable in dimless:
+        _to_unit = dimless[cmor_variable][to_unit]
+    else:
+        _to_unit = to_unit
+    if _to_unit == to_unit:
+        logger.info(
+            f"Converting units: ({da.name} -> {cmor_variable}) {from_unit} -> {to_unit}"
+        )
+    else:
+        logger.info(
+            f"Converting units: ({da.name} -> {cmor_variable}) {from_unit} -> {_to_unit} ({to_unit})"
+        )
+    new_da = new_da.pint.to(_to_unit).pint.dequantify()
     if new_da.attrs.get("units") != to_unit:
         logger.debug(
             "Pint auto-unit attribute setter different from requested unit string, setting manually."

tests/unit/test_dimensionless_units.py

@@ -0,0 +1,33 @@
+import numpy as np
+import xarray as xr
+
+from pymorize.data_request import DataRequestVariable
+from pymorize.rule import Rule
+from pymorize.units import handle_unit_conversion
+
+
+def test_units_with_psu():
+    da = xr.DataArray(10, name="sos", attrs={"units": "psu"})
+    rule = Rule(
+        cmor_variable="sos",
+        dimensionless_unit_mappings={"sos": {"0.001": "g/kg"}},
+    )
+    drv = DataRequestVariable(
+        variable_id="sos",
+        unit="0.001",
+        description="Salinity of the ocean",
+        time_method="MEAN",
+        table="Omon",
+        frequency="Monthly",
+        realms="Ocean",
+        standard_name="salinity_ocean",
+        cell_methods="area: mean where sea",
+        cell_measures="area: areacello",
+    )
+    rule.data_request_variable = drv
+    new_da = handle_unit_conversion(da, rule)
+    assert new_da.attrs.get("units") == "0.001"
+    # Check the magnitude of the data
+    # 1 g/kg is approximately 1 psu, so the values should be 10
+    # after conversion:
+    assert np.allclose(new_da.values, 10)