refining taylor diagram and associated changes

Author: brianpm

lib/adf_dataset.py

@@ -240,10 +240,20 @@ def load_reference_climo_ds(self, case, variablename, apply_scaling=True):
         fils = self.get_reference_climo_file(variablename)
         ds = self.load_dataset(fils)
         vname = self.ref_var_nam[variablename]  # name of variable in the reference data
+        # Check if already transformed (via attribute or units)
+        # Check if already transformed (via attribute or units)
+        unit_match = ds[vname].attrs.get('units') == self.adf.variable_defaults[variablename].get('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:
+            ds[vname].attrs['transformed'] = True
         return ds
 
     def load_reference_climo_da(self, case, variablename, apply_scaling=True):
@@ -391,20 +401,24 @@ def load_da(self, fils, variablename, **kwargs):
         if ds is None:
             warnings.warn(f"\t    WARNING: Load failed for {variablename}")
             return None
-        da = (ds[variablename]).squeeze()
+        da = ds[variablename].squeeze()
         apply_scaling = kwargs.get('apply_scaling', True)
         if not apply_scaling:
             add_offset = 0
             scale_factor = 1
         else:
             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
+                    da.attrs['transformed'] = True
         return da
 
     # Get variable conversion defaults, if applicable
@@ -433,8 +447,4 @@ def get_value_converters(self, case, variablename):
                     add_offset = vres.get("add_offset", 0)
         return add_offset, scale_factor
 
-    #------------------
-
-
-
-    
+    #------------------

lib/adf_utils.py

@@ -126,7 +126,7 @@ def mask_land_or_ocean(arr, msk, use_nan=False):
         missing_value = -999.
     #End if
 
-    arr = xr.where(msk>=0.9,arr,missing_value)
+    arr = xr.where(msk>=0.9,arr,missing_value, keep_attrs=True)
     arr.attrs["missing_value"] = missing_value
     return(arr)
 
@@ -535,7 +535,7 @@ def lev_to_plev(data, ps, hyam, hybm, P0=100000., new_levels=None,
     #Temporary print statement to notify users to ignore warning messages.
     #This should be replaced by a debug-log stdout filter at some point:
     print("Please ignore the interpolation warnings that follow!")
-
+    
     #Apply GeoCAT hybrid->pressure interpolation:
     if new_levels is not None:
         data_interp = gcomp.interpolation.interp_hybrid_to_pressure(data, ps,
@@ -559,10 +559,22 @@ def lev_to_plev(data, ps, hyam, hybm, P0=100000., new_levels=None,
     #Rename vertical dimension back to "lev" in order to work with
     #the ADF plotting functions:
     data_interp_rename = data_interp.rename({"plev": "lev"})
+    attrs = data_interp_rename.attrs.copy()
+    lev_orig = data_interp_rename["lev"]
+    lev_orig_attrs = lev_orig.attrs.copy()
 
     #Convert vertical dimension to mb/hPa, if requested:
     if convert_to_mb:
-        data_interp_rename["lev"] = data_interp_rename["lev"] / 100.0
+        lev_new = lev_orig / 100.0
+        lev_new.attrs = lev_orig_attrs
+        lev_new.name = "lev"
+        lev_new.attrs["units"] = "hPa"
+        lev_new.attrs["history"] = f"converted to hPa by dividing by 100 in adf_utils.lev_to_plev"
+        data_interp_rename["lev"] = lev_new
+        data_interp_rename.attrs = attrs
+    else:
+        data_interp_rename.attrs['units'] = "Pa"
+        data_interp_rename.attrs['history'] = f"Interpolated using GeoCAT, assume units of Pa in adf_utils.lev_to_plev"
 
     return data_interp_rename
 
@@ -678,6 +690,7 @@ def plev_to_plev(data, new_levels=None, convert_to_mb=False):
 
     # convert vertical dimension to mb/hPa, if requested:
     if convert_to_mb:
+        ## DEAL WITH METADATA BETTER HERE
         output["lev"] = output["lev"] / 100.0
     #End if