Taylor Diagram with observations and multi-case

lib/adf_config.py

@@ -237,20 +237,28 @@ def expand_references(self, config_dict):
         """
 
         #copy YAML config dictionary:
-        config_dict_copy = copy.copy(config_dict)
+        config_dict_copy = copy.deepcopy(config_dict)
 
-        #Loop through dictionary:
-        for key, value in config_dict_copy.items():
+        #Recursively expand references
+        self.__expand_refs(config_dict_copy)
 
-            #Skip non-strings (as they won't contain a keyword):
-            if not isinstance(value, str):
-                continue
+        #Update the original dict
+        config_dict.update(config_dict_copy)
 
-            #expand any keywords to their full values:
-            new_value = self.__expand_yaml_var_ref(value)
+    def __expand_refs(self, container):
 
-            #Set config variable to new, expanded value:
-            config_dict[key] = new_value
+        """
+        Recursive helper: expand keyword references in a nested dict or list, in place.
+        """
+
+        items = container.items() if isinstance(container, dict) else enumerate(container)
+
+        for key, value in items:
+            if isinstance(value, str):
+                #expand any keywords to their full values:
+                container[key] = self.__expand_yaml_var_ref(value)
+            elif isinstance(value, (dict, list)):
+                self.__expand_refs(value)
 
     #########
 

lib/adf_dataset.py

@@ -35,6 +35,17 @@
 #       Set AdfData.ref_nickname to that.
 #       Could be altered from "Obs" to be the data source label.
 
+# NOTE: Standard ADF workflow creates time series files with NCO.
+#       Climo files are then generated with create_climo_files.py
+#       Since neither of these apply units conversions (add_offset, scale_factor),
+#       the methods here default to applying them when loading 
+#       time series and climo files, using the kwarg apply_scaling.
+#       Regridded files are made with regrid_and_vert_interp[_2].py,
+#       which uses this module for loading climo files, so will apply
+#       scaling.
+#       Therefore the default on loading regridded files is to NOT
+#       apply scaling.
+
 class AdfData:
     """A class instantiated with an AdfDiag object.
        Methods provide means to load data.
@@ -197,9 +208,12 @@ def load_timeseries_da(self, case, variablename):
             return None
         return self.load_da(fils, variablename, add_offset=add_offset, scale_factor=scale_factor)
 
-    def load_reference_timeseries_da(self, field):
+    def load_reference_timeseries_da(self, field, apply_scaling=True):
         """Return a DataArray time series to be used as reference
           (aka baseline) for variable field.
+
+          apply_scaling: bool
+            If True, apply add_offset and scale_factor to data (if present).
         """
         fils = self.get_ref_timeseries_file(field)
         if not fils:
@@ -215,6 +229,10 @@ def load_reference_timeseries_da(self, field):
         else:
             add_offset, scale_factor = self.get_value_converters(self.ref_case_label, field)
 
+        if not apply_scaling:
+            add_offset = 0
+            scale_factor = 1
+
         return self.load_da(fils, field, add_offset=add_offset, scale_factor=scale_factor)
 
 
@@ -225,10 +243,38 @@ def load_reference_timeseries_da(self, field):
     #------------------
 
     # Test case(s)
-    def load_climo_da(self, case, variablename):
-        """Return DataArray from climo file"""
+    def load_climo_ds(self, case, variablename):
+        """Return Dataset from climo file; applies scale factor and offset to `variablename`."""
         add_offset, scale_factor = self.get_value_converters(case, variablename)
         fils = self.get_climo_file(case, variablename)
+        if not fils:
+            warnings.warn("\t    WARNING: Did not find climo file for case: "
+                          f"{case}, variable: {variablename}")
+            return None
+        ds = self.load_dataset(fils)
+        if ds is None:
+            return None
+        # xarray arithmetic drops attrs, so carry them across by hand -- otherwise
+        # the regridded files lose 'units' and the plotting scripts KeyError on it.
+        attrs = ds[variablename].attrs.copy()
+        ds[variablename] = ds[variablename] * scale_factor + add_offset
+        ds[variablename].attrs = attrs
+        if scale_factor != 1 or add_offset != 0:
+            new_unit = self.adf.variable_defaults.get(variablename, {}).get("new_unit")
+            if new_unit:
+                ds[variablename].attrs['units'] = new_unit
+                # int, not bool: netCDF4 cannot store a Python bool as an attribute
+                ds[variablename].attrs['transformed'] = 1
+        return ds
+
+    def load_climo_da(self, case, variablename, apply_scaling=True):
+        """Return DataArray from climo file"""
+        if not apply_scaling:
+            add_offset = 0
+            scale_factor = 1
+        else:
+            add_offset, scale_factor = self.get_value_converters(case, variablename)
+        fils = self.get_climo_file(case, variablename)
         return self.load_da(fils, variablename, add_offset=add_offset, scale_factor=scale_factor)
 
 
@@ -264,11 +310,46 @@ def get_climo_file(self, case, variablename):
 
 
     # Reference case (baseline/obs)
-    def load_reference_climo_da(self, case, variablename):
-        """Return DataArray from reference (aka baseline) climo file"""
+    def load_reference_climo_ds(self, case, variablename, apply_scaling=True):
+        """Return Dataset from reference climo file; applies scale factor and offset to `variablename`.
+        """
         add_offset, scale_factor = self.get_value_converters(case, variablename)
         fils = self.get_reference_climo_file(variablename)
-        return self.load_da(fils, variablename, add_offset=add_offset, scale_factor=scale_factor)
+        if not fils:
+            warnings.warn("\t    WARNING: Did not find reference climo file for "
+                          f"variable: {variablename}")
+            return None
+        ds = self.load_dataset(fils)
+        if ds is None:
+            return None
+        vname = self.ref_var_nam[variablename]  # name of variable in the reference data
+        # Check if already transformed (via attribute or units)
+        new_unit = self.adf.variable_defaults.get(variablename, {}).get('new_unit')
+        unit_match = new_unit is not None and ds[vname].attrs.get('units') == new_unit
+        if ds[vname].attrs.get('transformed', False) or unit_match:
+            apply_scaling = False
+        if not apply_scaling:
+            add_offset = 0
+            scale_factor = 1
+
+        attrs = ds[vname].attrs.copy()
+        ds[vname] = ds[vname] * scale_factor + add_offset
+        ds[vname].attrs = attrs
+        if scale_factor != 1 or add_offset != 0:
+            # int, not bool: netCDF4 cannot store a Python bool as an attribute
+            ds[vname].attrs['transformed'] = 1
+        return ds
+
+    def load_reference_climo_da(self, case, variablename, apply_scaling=True):
+        """Return DataArray from reference (aka baseline) climo file"""
+        fils = self.get_reference_climo_file(variablename)
+        vname = self.ref_var_nam[variablename]
+        if not apply_scaling:
+            add_offset = 0
+            scale_factor = 1
+        else:
+            add_offset, scale_factor = self.get_value_converters(case, variablename)
+        return self.load_da(fils, vname, add_offset=add_offset, scale_factor=scale_factor)
 
     def get_reference_climo_file(self, var):
         """Return a list of files to be used as reference (aka baseline) for variable var."""
@@ -294,7 +375,7 @@ def get_reference_climo_file(self, var):
     # Test case(s)
     def get_regrid_file(self, case, field):
         """Return list of test regridded files"""
-        model_rg_loc = Path(self.adf.get_basic_info("cam_regrid_loc", required=True))
+        model_rg_loc = Path(self.model_rgrid_loc)
         # rlbl = "reference label" = name of the reference data that defines the target grid
         rlbl = self.ref_labels[field]
         return sorted(model_rg_loc.glob(f"{rlbl}_{case}_{field}_regridded.nc"))
@@ -310,9 +391,13 @@ def load_regrid_dataset(self, case, field):
         return self.load_dataset(fils)
 
 
-    def load_regrid_da(self, case, field):
+    def load_regrid_da(self, case, field, apply_scaling=False):
         """Return a data array to be used as reference (aka baseline) for variable field."""
-        add_offset, scale_factor = self.get_value_converters(case, field)
+        if not apply_scaling:
+            add_offset = 0
+            scale_factor = 1
+        else:
+            add_offset, scale_factor = self.get_value_converters(case, field)
         fils = self.get_regrid_file(case, field)
         if not fils:
             warnings.warn("\t    WARNING: Did not find regrid file(s) for case: "
@@ -331,7 +416,7 @@ def get_ref_regrid_file(self, case, field):
             else:
                 fils = []
         else:
-            model_rg_loc = Path(self.adf.get_basic_info("cam_regrid_loc", required=True))
+            model_rg_loc = Path(self.model_rgrid_loc)
             fils = sorted(model_rg_loc.glob(f"{case}_{field}_baseline.nc"))
         return fils
 
@@ -346,9 +431,13 @@ def load_reference_regrid_dataset(self, case, field):
         return self.load_dataset(fils)
 
 
-    def load_reference_regrid_da(self, case, field):
+    def load_reference_regrid_da(self, case, field, apply_scaling=False):
         """Return a data array to be used as reference (aka baseline) for variable field."""
-        add_offset, scale_factor = self.get_value_converters(case, field)
+        if not apply_scaling:
+            add_offset = 0
+            scale_factor = 1
+        else:
+            add_offset, scale_factor = self.get_value_converters(case, field)
         fils = self.get_ref_regrid_file(case, field)
         if not fils:
             warnings.warn("\t    WARNING: Did not find regridded file(s) for case: "
@@ -360,13 +449,9 @@ def load_reference_regrid_da(self, case, field):
             field = self.ref_var_nam[field]
         return self.load_da(fils, field, add_offset=add_offset, scale_factor=scale_factor)
 
-    #------------------
-
-
+    #---------------------------
     # DataSet and DataArray load
     #---------------------------
-
-    # Load DataSet
     def load_dataset(self, fils):
         """Return xarray DataSet from file(s)"""
         if len(fils) == 0:
@@ -384,7 +469,6 @@ def load_dataset(self, fils):
             warnings.warn("\t    WARNING: invalid data on load_dataset")
         return ds
 
-    # Load DataArray
     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)
@@ -394,12 +478,17 @@ def load_da(self, fils, variablename, **kwargs):
         da = ds[variablename].squeeze()
         scale_factor = kwargs.get('scale_factor', 1)
         add_offset = kwargs.get('add_offset', 0)
+        attrs = da.attrs.copy()
         da = da * scale_factor + add_offset
-        if variablename in self.adf.variable_defaults:
-            vres = self.adf.variable_defaults[variablename]
-            da.attrs['units'] = vres.get("new_unit", da.attrs.get('units', 'none'))
-        else:
-            da.attrs['units'] = 'none'
+        da.attrs = attrs
+
+        if scale_factor != 1 or add_offset != 0:
+            if variablename in self.adf.variable_defaults:
+                new_unit = self.adf.variable_defaults[variablename].get("new_unit")
+                if new_unit:
+                    da.attrs['units'] = new_unit
+                    # int, not bool: netCDF4 cannot store a Python bool as an attribute
+                    da.attrs['transformed'] = 1
         return da
 
     # Get variable conversion defaults, if applicable
@@ -429,3 +518,4 @@ def get_value_converters(self, case, variablename):
         return add_offset, scale_factor
 
     #------------------
+

lib/adf_info.py

@@ -125,10 +125,9 @@ def __init__(self, config_file, debug=False):
                     emsg += f" {self.__num_cases} entries, instead it has {len(conf_val)}"
                     self.end_diag_fail(emsg)
             else:
-                #If not a list, then convert it to one:
-                self.__cam_climo_info[conf_var] = [conf_val]
-            #End if
-        #End for
+                # If not a list, replicate the scalar value for each case
+                self.__cam_climo_info[conf_var] = [conf_val] * self.__num_cases
+        # End for
 
         #Initialize ADF variable list:
         self.__diag_var_list = self.read_config_var('diag_var_list', required=True)
@@ -371,12 +370,14 @@ def __init__(self, config_file, debug=False):
 
         #Make lists of None to be iterated over for case_names
         if syears is None:
-            syears = [None]*len(case_names)
-        #End if
-        if eyears is None:
-            eyears = [None]*len(case_names)
-        #End if
+            syears = [None] * self.__num_cases
+        elif not isinstance(syears, list):
+            syears = [syears] * self.__num_cases
 
+        if eyears is None:
+            eyears = [None] * self.__num_cases
+        elif not isinstance(eyears, list):
+            eyears = [eyears] * self.__num_cases
         #Extract cam history files location:
         cam_hist_locs = self.get_cam_info('cam_hist_loc')
 
@@ -447,13 +448,13 @@ def __init__(self, config_file, debug=False):
             hist_str_case = hist_str[case_idx]
             if any(cam_hist_locs):
                 #Grab first possible hist string, just looking for years of run
-                hist_str = hist_str_case[0]
+                hist_str_use = hist_str_case[0]
 
                 #Get climo years for verification or assignment if missing
                 starting_location = Path(cam_hist_locs[case_idx])
                 print(f"\tChecking history files in '{starting_location}'")
 
-                file_list = sorted(starting_location.glob('*'+hist_str+'.*.nc'))
+                file_list = sorted(starting_location.glob('*'+hist_str_use+'.*.nc'))
 
                 #Check if the history file location exists
                 if not starting_location.is_dir():
@@ -466,7 +467,7 @@ def __init__(self, config_file, debug=False):
                     self.end_diag_fail(emsg)
 
                 #Check if there are any history files
-                file_list = sorted(starting_location.glob('*'+hist_str+'.*.nc'))
+                file_list = sorted(starting_location.glob('*'+hist_str_use+'.*.nc'))
                 if len(file_list) == 0:
                     msg = "Checking history files:\n"
                     msg += f"\tThere are no history files in '{starting_location}'."
@@ -485,7 +486,7 @@ def __init__(self, config_file, debug=False):
                 #Since the last part always includes the time range, grab that with last index (2)
                 #NOTE: this is based off the current CAM file name structure in the form:
                 #  $CASE.cam.h#.YYYY<other date info>.nc
-                case_climo_yrs = [int(str(i).partition(f"{hist_str}.")[2][0:4]) for i in file_list]
+                case_climo_yrs = [int(str(i).partition(f"{hist_str_use}.")[2][0:4]) for i in file_list]
                 if not case_climo_yrs:
                     msg = f"\t ERROR: No climo years found in {cam_hist_locs[case_idx]}, "
                     raise AdfError(msg)
@@ -619,15 +620,34 @@ def __init__(self, config_file, debug=False):
 
     def hist_str_to_list(self, conf_var, conf_val):
         """
-        Make hist_str a nested list [ncases,nfiles] of the given value(s)
+        Normalizes hist_str input into a nested list [ncases][nfiles].
         """
-        if isinstance(conf_val, list):
-            hist_str = conf_val
-        else:  # one case, one hist str
-            hist_str = [
-                conf_val
-            ]
-        self.__cam_climo_info[conf_var] = [hist_str]
+        n = self.__num_cases
+        result = None
+
+        # 1. Handle Single String input: "h0" -> [["h0"], ["h0"], ...]
+        if isinstance(conf_val, str):
+            result = [[conf_val] for _ in range(n)]
+
+        elif isinstance(conf_val, list):
+            # 2. Check if it's already a nested list: [["h0"], ["h0"]]
+            # We check the first element to see if it's a list.
+            if len(conf_val) == n and all(isinstance(i, list) for i in conf_val):
+                result = conf_val
+
+            # 3. Check if it's a list of strings matching N cases: ["h0", "h1"]
+            elif len(conf_val) == n and all(isinstance(i, str) for i in conf_val):
+                result = [[i] for i in conf_val]
+
+            # 4. Otherwise, treat it as a single set of files for ALL cases: ["h0", "h1"]
+            else:
+                # We wrap the list and multiply it
+                result = [conf_val for _ in range(n)]
+
+        if result is None:
+            raise ValueError(f"Invalid format for {conf_var}: {conf_val}")
+
+        self.__cam_climo_info[conf_var] = result
     #########
 
     # Create property needed to return "user" name to user: