Apply the variable defaults and the recorded time interval to AMWG tables (#423)

lib/adf_dataset.py

@@ -213,18 +213,16 @@ def load_timeseries_dataset(self, fils):
             ds = xr.open_dataset(sfil, decode_times=False)
         if ds is None:
             warnings.warn("\t    WARNING: invalid data on load_dataset")
-        # assign time to midpoint of interval (even if it is already)
-        if 'time_bnds' in ds:
-            t = ds['time_bnds'].mean(dim='nbnd')
-            t.attrs = ds['time'].attrs
-            ds = ds.assign_coords({'time':t})
-        elif 'time_bounds' in ds:
-            t = ds['time_bounds'].mean(dim='nbnd')
-            t.attrs = ds['time'].attrs
-            ds = ds.assign_coords({'time':t})
-        else:
+            return ds
+        # Assign time to the midpoint of the interval each step covers.  The
+        # shared helper reads the bounds the file names for itself, so a file
+        # calling them something other than 'time_bnds' is handled too, and it
+        # hands back the dataset it was given when the file records no bounds:
+        fixed = utils.use_time_bounds_midpoint(ds)
+        if fixed is ds:
             warnings.warn("\t    INFO: Timeseries file does not have time bounds info.")
-        return xr.decode_cf(ds)
+        # End if
+        return xr.decode_cf(fixed)
 
     def load_timeseries_da(self, case, variablename):
         """Return DataArray from time series file(s).
@@ -236,7 +234,13 @@ def load_timeseries_da(self, case, variablename):
             warnings.warn("\t    WARNING: Did not find case time series file(s), "
                           f"variable: {variablename}")
             return None
-        return self.load_da(fils, variablename, add_offset=add_offset, scale_factor=scale_factor)
+        return self.load_da(
+            fils,
+            variablename,
+            use_time_bounds=True,
+            add_offset=add_offset,
+            scale_factor=scale_factor,
+        )
 
     def load_reference_timeseries_da(self, field, apply_scaling=True):
         """Return a DataArray time series to be used as reference
@@ -263,7 +267,13 @@ def load_reference_timeseries_da(self, field, apply_scaling=True):
             add_offset = 0
             scale_factor = 1
 
-        return self.load_da(fils, field, add_offset=add_offset, scale_factor=scale_factor)
+        return self.load_da(
+            fils,
+            field,
+            use_time_bounds=True,
+            add_offset=add_offset,
+            scale_factor=scale_factor,
+        )
 
 
     #------------------
@@ -529,8 +539,15 @@ def _regrid_converters(self, fils, file_field, case, field, apply_scaling):
     #---------------------------
     # DataSet and DataArray load
     #---------------------------
-    def load_dataset(self, fils):
-        """Return xarray DataSet from file(s)"""
+    def load_dataset(self, fils, use_time_bounds=False):
+        """Return xarray DataSet from file(s).
+
+        `use_time_bounds` moves the time coordinate to the midpoint of the
+        interval each step covers, which is what a time series wants.  It is
+        off by default: climatology and regridded files carry a time
+        coordinate of month numbers, and turning that into dates would change
+        the files the ADF writes and reads back.
+        """
         if len(fils) == 0:
             warnings.warn("\t    WARNING: Input file list is empty.")
             return None
@@ -544,11 +561,20 @@ def load_dataset(self, fils):
             ds = xr.open_dataset(sfil)
         if ds is None:
             warnings.warn("\t    WARNING: invalid data on load_dataset")
+            return ds
+        if use_time_bounds:
+            # Time stamps that name one end of an averaging interval put steps
+            # in the wrong year, so use what the file records about it:
+            ds = utils.use_time_bounds_midpoint(ds)
+        # End if
         return ds
 
-    def load_da(self, fils, variablename, **kwargs):
-        """Return xarray DataArray from file(s) w/ optional scale factor, offset, new units."""
-        ds = self.load_dataset(fils)
+    def load_da(self, fils, variablename, use_time_bounds=False, **kwargs):
+        """Return xarray DataArray from file(s) w/ optional scale factor, offset, new units.
+
+        `use_time_bounds` is passed to `load_dataset`; see there.
+        """
+        ds = self.load_dataset(fils, use_time_bounds=use_time_bounds)
         if ds is None:
             warnings.warn(f"\t    WARNING: Load failed for {variablename}")
             return None

lib/adf_info.py

@@ -1008,22 +1008,9 @@ def ts_patterns(var):
         else:
             cam_ts_data = xr.open_mfdataset(ts_files, decode_times=True, combine='by_coords')
 
-        #Average time dimension over time bounds, if bounds exist:
-        if 'time_bnds' in cam_ts_data:
-            time_bounds_name = 'time_bnds'
-        elif 'time_bounds' in cam_ts_data:
-            time_bounds_name = 'time_bounds'
-        else:
-            time_bounds_name = None
-
-        if time_bounds_name:
-            time = cam_ts_data['time']
-            #NOTE: force `load` here b/c if dask & time is cftime, throws a NotImplementedError:
-            time = xr.DataArray(cam_ts_data[time_bounds_name].load().mean(dim='nbnd').values,
-                                dims=time.dims, attrs=time.attrs)
-            cam_ts_data['time'] = time
-            cam_ts_data.assign_coords(time=time)
-            cam_ts_data = xr.decode_cf(cam_ts_data)
+        # Use the interval each step covers rather than its stamp, so that the
+        # years found here are the years the data actually covers:
+        cam_ts_data = utils.use_time_bounds_midpoint(cam_ts_data)
 
         #Extract first and last years from dataset:
         syr = int(cam_ts_data.time[0].dt.year.values)
@@ -1036,4 +1023,4 @@ def ts_patterns(var):
 
 #++++++++++++++++++++
 #End Class definition
-#++++++++++++++++++++
+#++++++++++++++++++++

lib/adf_utils.py

@@ -6,6 +6,10 @@
 find_ts_files(), select_ts_files(), ts_files_overlap(), ts_file_span(),
 as_hist_str_list(), pick_hist_str()
     re-exported from adf_file_utils; time series file discovery
+plain_text_units()
+    render a unit string from the variable defaults as plain text
+use_time_bounds_midpoint()
+    set the time coordinate to the midpoint of the interval each step covers
 load_dataset()
     generalized load dataset method used for plotting/analysis functions
 mask_land_or_ocean(arr, msk, use_nan=False)
@@ -90,7 +94,120 @@ def my_formatwarning(msg, *args, **kwargs):
 #HELPER FUNCTIONS
 #################
 
-def load_dataset(fils):
+# Simple LaTeX seen in the units of the variable defaults, and what each means
+# as plain text.  The defaults were written for plot labels, where matplotlib
+# renders them; a table written as text or HTML shows them as they are.
+_LATEX_UNIT_REPLACEMENTS = (
+    ("\\mu", "u"),
+    ("^{", "^"),
+    ("}", ""),
+    ("$", ""),
+    ("\\,", " "),
+    ("\\", ""),
+)
+
+
+def plain_text_units(unit):
+    """
+    Render a unit string as plain text.
+
+    The units in the variable defaults are written for plot labels, where
+    matplotlib renders the LaTeX in them, so ``mm d$^{-1}$`` appears as mm
+    d^-1 with a proper superscript.  Written into a table -- a csv file, or a
+    web page -- the same string appears exactly as it is typed, which is not
+    something to put in front of a reader.
+
+    Parameters
+    ----------
+    unit : str
+        the unit as the variable defaults give it
+
+    Returns
+    -------
+    str
+        the same unit with its markup resolved: ``mm d$^{-1}$`` becomes
+        ``mm d^-1``.  A unit with no markup is returned unchanged, so this is
+        safe to apply to anything and does nothing to units that are already
+        plain.
+    """
+    if not isinstance(unit, str) or ("$" not in unit and "\\" not in unit):
+        return unit
+    # End if
+    for markup, plain in _LATEX_UNIT_REPLACEMENTS:
+        unit = unit.replace(markup, plain)
+    # End for
+    return unit.strip()
+
+
+def use_time_bounds_midpoint(ds, time_name="time"):
+    """
+    Set the time coordinate to the midpoint of the interval each step covers.
+
+    CAM stamps a monthly average with one end of the interval it covers, and
+    which end depends on the model version: an older CAM h0 file stamps January
+    with February 1st.  Anything that then asks which year a step belongs to --
+    selecting years, averaging by year, grouping by month -- puts that step in
+    the wrong place, which silently drops a year from an annual mean and shifts
+    a seasonal cycle by a month.  The interval itself is not in doubt: the file
+    records it, so the midpoint of the recorded interval is used instead.
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+        dataset to correct
+    time_name : str, optional
+        name of the time coordinate; defaults to ``"time"``
+
+    Returns
+    -------
+    xr.Dataset
+        a dataset whose time coordinate is the midpoint of its bounds, or the
+        dataset unchanged when the file does not say what its bounds are.
+
+    Notes
+    -----
+    The bounds variable is taken from the ``bounds`` attribute of the time
+    coordinate, which is where CF says it belongs and is authoritative when it
+    is there.  Only if that attribute is missing, or names something the file
+    does not contain, are the conventional names ``time_bnds`` and
+    ``time_bounds`` tried.  Files that record no bounds are returned untouched
+    -- there is then nothing better to go on than the stamp itself.
+
+    Doing this where files are opened, rather than in each script, is
+    deliberate: a script that does not know about the stamping convention
+    should not have to.
+    """
+    if time_name not in ds.variables:
+        return ds
+    # End if
+
+    # What the file itself says its bounds are, then the conventional names:
+    candidates = [ds[time_name].attrs.get("bounds"), "time_bnds", "time_bounds"]
+    bounds_name = next(
+        (name for name in candidates if name and name in ds.variables), None
+    )
+    if bounds_name is None:
+        return ds
+    # End if
+
+    bounds = ds[bounds_name]
+    # The bounds are (time, 2), but the second dimension is called 'nbnd' by
+    # CAM and other names elsewhere, so take whichever one is not time:
+    other_dims = [dim for dim in bounds.dims if dim != time_name]
+    if len(other_dims) != 1:
+        # Not a shape this can make sense of, so leave the file alone:
+        return ds
+    # End if
+
+    # load() first: averaging cftime under dask raises NotImplementedError.
+    midpoint = bounds.load().mean(dim=other_dims[0])
+    attrs = ds[time_name].attrs
+    ds = ds.assign_coords({time_name: midpoint})
+    ds[time_name].attrs = attrs
+    return ds
+
+
+def load_dataset(fils, use_time_bounds=False):
     """
     This method exists to get an xarray Dataset from input file information that can be passed into the plotting methods.
 
@@ -106,15 +223,27 @@ def load_dataset(fils):
     Notes
     -----
     When just one entry is provided, use `open_dataset`, otherwise `open_mfdatset`
+
+    Pass ``use_time_bounds=True`` for time series files, so that steps stamped
+    at one end of their averaging interval are counted in the right month and
+    year.  See `use_time_bounds_midpoint`.
     """
     if len(fils) == 0:
         warnings.warn(f"\t    WARNING: Input file list is empty.")
         return None
     elif len(fils) > 1:
-        return xr.open_mfdataset(fils, combine='by_coords')
+        ds = xr.open_mfdataset(fils, combine="by_coords")
     else:
-        return xr.open_dataset(fils[0])
-    #End if
+        ds = xr.open_dataset(fils[0])
+    # End if
+    if use_time_bounds:
+        # Time stamps that name one end of an averaging interval put steps in
+        # the wrong year, so use what the file records about it instead.  Off
+        # by default: this function also reads climatology files, whose time
+        # coordinate is month numbers and should stay that way.
+        ds = use_time_bounds_midpoint(ds)
+    # End if
+    return ds
 #End def