Merge PR #397 (Fix errors in annual mean diff-of-diffs plots)

This merge brings PR #397 (Fix errors in annual mean GCHP vs
GC-Classic "diff-of-diffs" plots, by @yantosca) into the
GCPy 1.7.0 development stream.

In PR #397 We have done the following:

1. Updated functions compare_single_level and compare_zonal_mean so
   as to not take the time mean of data if we are generating
   diff-of-diff plots.

2. Simplified function get_diff_of_diffs (in gcpy/util.py) so that a
   single algorithm is used for both GC-Classic and GCHP.

3. Now take the time mean of absolute and fractional differences in
   get_diff_of_diffs

4. Added cosmetic changes (comments)

Signed-off-by: Bob Yantosca <yantosca@seas.harvard.edu>

Author: yantosca

CHANGELOG.md

@@ -52,12 +52,14 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 - Fixed the restart regridding for stretched GCHP when target lat/lon is exactly 0.0 in `gcpy/regrid_restart_file.py`
 - Fixed computation of the global AOD benchmark table caused by hardwired species names
 - Fixed error in `create_benchmark_emissions_table` where all species were assumed to be in Ref and Dev even if they were not
+- Fixed error that caused the GCHP vs GCC diff-of-diffs AnnualMean plots to be computed as a difference of means instead of a mean of differences
 
 ### Removed
 - Removed `PdfMerger()` from `compare_single_level` and `compare_zonal_mean`, it has been removed in pypdf >= 5.0.0
 - Removed `.load()` statements from xarray Datasets to improve performance
 - Removed `paths:spcdb_dir` YAML tag in benchmark configuration files
 - Removed `st_Ox` from `benchmark_categories.yml`; this species is no longer used in TransportTracers simulations
+- Removed special data handling for files generated with MAPL versions prior to 1.0.0 in function `get_diff_of_diffs` (located in `gcpy/util.py`)
 
 ## [1.6.2] - 2025-06-12
 ### Added

gcpy/benchmark/modules/benchmark_funcs.py

@@ -1060,21 +1060,33 @@ def make_benchmark_conc_plots(
         devmetds = reader(devmet, drop_variables=SKIP_THESE_VARS)
 
     # Determine if doing diff-of-diffs
-    diff_of_diffs = False
-    if second_ref is not None and second_dev is not None:
-        diff_of_diffs = True
+    diff_of_diffs = second_ref is not None and second_dev is not None
 
-
-    # Open second datasets if passed as arguments (used for diff of diffs)
-    # Regrid to same horz grid resolution if two refs or two devs do not match.
     if diff_of_diffs:
+
+        # --------------------------------------------------------------
+        # %%%%% We are plotting diff-of-diffs %%%%%
+        #
+        # Open the second Ref and Dev datasets, if they have been
+        # passed as keyword arguments, as these are needed for the
+        # diff-of-diffs plots.  Regrid to the same # horizontal grid
+        # resolution if the two Refs or two Devs are not on the same
+        # grid.
+        #
+        # Also, do not take the time mean (e.g. Annual Mean) of the
+        # datasets, as this will compute the the difference of means.
+        # Instead, we need to compute the mean of differences.  This
+        # will be done in the plotting functions compare_single_level
+        # and compare_zonal_mean.
+        # --------------------------------------------------------------
         second_refds = reader(second_ref, drop_variables=SKIP_THESE_VARS)
         second_devds = reader(second_dev, drop_variables=SKIP_THESE_VARS)
-
-        print('\nPrinting second_refds (dev of ref for diff-of-diffs)\n')
-        print(second_refds)
-        print('\nPrinting second_devds (dev of dev for diff-of-diffs)\n')
-        print(second_devds)
+        
+        if verbose:
+            print('\nPrinting second_refds (dev of ref for diff-of-diffs)\n')
+            print(second_refds)
+            print('\nPrinting second_devds (dev of dev for diff-of-diffs)\n')
+            print(second_devds)
 
         # Only regrid if ref grid resolutions differ.
         # Assume only may differ if both cubed-sphere.
@@ -1088,6 +1100,7 @@ def make_benchmark_conc_plots(
         if "Xdim" in devds.dims and "Xdim" in second_devds.dims:
             if devds['Xdim'].size != second_devds['Xdim'].size:
                 regrid_dev = True
+
         if regrid_ref or regrid_dev:
             # Assume regridding C24 to C48 to compute the difference at C48
             regridfile=os.path.join(weightsdir,'regrid_weights_c24_to_c48.nc')
@@ -1118,17 +1131,20 @@ def make_benchmark_conc_plots(
                 print('\nRegrid complete\n')
 
     else:
+
+        # --------------------------------------------------------------
+        # %%%%% We are not plotting diff-of-diffs %%%%%
+        #
+        # We do not need the second Ref and Dev data.
+        # We can also compute the time mean (e.g. Annual Mean) here.
+        # --------------------------------------------------------------
         second_refds = None
         second_devds = None
-
-    # Compute the annual mean of datasets (if necessary)
-    if time_mean:
-        refds = dataset_mean(refds)
-        devds = dataset_mean(devds)
-        refmetds = dataset_mean(refmetds)
-        devmetds = dataset_mean(devmetds)
-        second_ref = dataset_mean(second_refds)
-        second_dev = dataset_mean(second_devds)
+        if time_mean:
+            refds = dataset_mean(refds)
+            devds = dataset_mean(devds)
+            refmetds = dataset_mean(refmetds)
+            devmetds = dataset_mean(devmetds)
 
     # Create regridding files if necessary while not in parallel loop
     [_ for _ in create_regridders(refds, devds, weightsdir=weightsdir)]

gcpy/plot/compare_single_level.py

@@ -205,17 +205,18 @@ def compare_single_level(
     # Prepare diff-of-diffs datasets if needed
     if diff_of_diffs:
 
-#        # If needed, use fake time dim in case dates are different
-#        # in datasets.  This needs more work for case of single versus
-#        # multiple times.
-#        aligned_time = [np.datetime64('2000-01-01')] * refdata.dims['time']
-#        refdata = refdata.assign_coords({'time': aligned_time})
-#        devdata = devdata.assign_coords({'time': aligned_time})
-#        second_ref = second_ref.assign_coords({'time': aligned_time})
-#        second_dev = second_dev.assign_coords({'time': aligned_time})
+        ## If needed, use fake time dim in case dates are different
+        ## in datasets.  This needs more work for case of single versus
+        ## multiple times.
+        #aligned_time = [np.datetime64('2000-01-01')] * refdata.dims['time']
+        #refdata = refdata.assign_coords({'time': aligned_time})
+        #devdata = devdata.assign_coords({'time': aligned_time})
+        #second_ref = second_ref.assign_coords({'time': aligned_time})
+        #second_dev = second_dev.assign_coords({'time': aligned_time})
 
         refdata, fracrefdata = get_diff_of_diffs(refdata, second_ref)
         devdata, fracdevdata = get_diff_of_diffs(devdata, second_dev)
+
         frac_refstr = 'GCC_dev / GCC_ref'
         frac_devstr = 'GCHP_dev / GCHP_ref'
     # If no varlist is passed, plot all (surface only for 3D)

gcpy/plot/compare_zonal_mean.py

@@ -212,13 +212,13 @@ def compare_zonal_mean(
     # Prepare diff-of-diffs datasets if needed
     if diff_of_diffs:
 
-#        # If needed, use fake time dim in case dates are different in datasets.
-#        # This needs more work for case of single versus multiple times.
-#        aligned_time = np.datetime64('2000-01-01')
-#        refdata = refdata.assign_coords({'time' : [aligned_time]})
-#        devdata = devdata.assign_coords({'time' : [aligned_time]})
-#        second_ref = second_ref.assign_coords({'time' : [aligned_time]})
-#        second_dev = second_dev.assign_coords({'time' : [aligned_time]})
+        ## If needed, use fake time dim in case dates are different in datasets.
+        ## This needs more work for case of single versus multiple times.
+        #aligned_time = np.datetime64('2000-01-01')
+        #refdata = refdata.assign_coords({'time' : [aligned_time]})
+        #devdata = devdata.assign_coords({'time' : [aligned_time]})
+        #second_ref = second_ref.assign_coords({'time' : [aligned_time]})
+        #second_dev = second_dev.assign_coords({'time' : [aligned_time]})
 
         refdata, fracrefdata = get_diff_of_diffs(refdata, second_ref)
         devdata, fracdevdata = get_diff_of_diffs(devdata, second_dev)

gcpy/util.py

@@ -634,70 +634,61 @@ def get_diff_of_diffs(
         dev
 ):
     """
-    Generate datasets containing differences between two datasets
+    Generate datasets containing differences between two datasets.
 
     Args:
-        ref: xarray Dataset
-            The "Reference" (aka "Ref") dataset.
-        dev: xarray Dataset
-            The "Development" (aka "Dev") dataset
+    ref       : xr.Dataset : The Ref (aka "Reference") dataset
+    dev       : xr.Dataset : The Dev" (aka "Development") dataset
 
-    Returns:
-         absdiffs: xarray Dataset
-            Dataset containing dev-ref values
-         fracdiffs: xarray Dataset
-            Dataset containing dev/ref values
-    """
+    Returns
+    absdiffs  : xr.Dataset : Absolute differences (Dev - Ref)
+    fracdiffs : xr.Dataset : Fractional differences (Dev / Ref)
+    """
+    # ------------------------------------------------------------------
+    # Throw an error unless Ref and Dev are on the same grid
+    # ------------------------------------------------------------------
+    if not ref.sizes == dev.sizes:
+        msg = "Diff-of-diffs plot supports only identical grid types "
+        msg += "(lat/lon or cubed-sphere) within each Ref & Dev pair!"
+        raise ValueError(msg)
 
-    # get diff of diffs datasets for 2 datasets
-    # limit each pair to be the same type of output (GEOS-Chem Classic or GCHP)
-    # and same resolution / extent
+    # ------------------------------------------------------------------
+    # Include only the common fields in Ref and Dev
+    # ------------------------------------------------------------------
     vardict = compare_varnames(ref, dev, quiet=True)
     varlist = vardict["commonvars"]
-    # Select only common fields between the Ref and Dev datasets
     ref = ref[varlist]
     dev = dev[varlist]
-    if 'nf' not in ref.dims and 'nf' not in dev.dims:
+
+    # ------------------------------------------------------------------
+    # Compute absolute (Dev-Ref) and fractional (Dev/Ref) differences
+    # ------------------------------------------------------------------
+    with xr.set_options(keep_attrs=True):
+
         # if the coords do not align then set time dimensions equal
         try:
             xr.align(dev, ref, join='exact')
         except BaseException:
             ref.coords["time"] = dev.coords["time"]
-        with xr.set_options(keep_attrs=True):
-            absdiffs = dev - ref
-            fracdiffs = dev / ref
-            for var in dev.data_vars.keys():
-                # Ensure the diffs Dataset includes attributes
-                absdiffs[var].attrs = dev[var].attrs
-                fracdiffs[var].attrs = dev[var].attrs
-    elif 'nf' in ref.dims and 'nf' in dev.dims:
-
-        # Include special handling if cubed sphere grid dimension names are different
-        # since they changed in MAPL v1.0.0.
-        if "lat" in ref.dims and "Xdim" in dev.dims:
-            ref_newdimnames = dev.copy()
-            for var in dev.data_vars.keys():
-                if "Xdim" in dev[var].dims:
-                    ref_newdimnames[var].values = ref[var].values.reshape(
-                        dev[var].values.shape)
-                # NOTE: the reverse conversion is gchp_dev[v].stack(lat=("nf","Ydim")).transpose(
-                # "time","lev","lat","Xdim").values
 
-        with xr.set_options(keep_attrs=True):
-            absdiffs = dev.copy()
-            fracdiffs = dev.copy()
-            for var in dev.data_vars.keys():
-                if "Xdim" in dev[var].dims or "lat" in dev[var].dims:
-                    absdiffs[var].values = dev[var].values - ref[var].values
-                    fracdiffs[var].values = dev[var].values / ref[var].values
-                    # NOTE: The diffs Datasets are created without variable
-                    # attributes; we have to reattach them
-                    absdiffs[var].attrs = dev[var].attrs
-                    fracdiffs[var].attrs = dev[var].attrs
-    else:
-        print('Diff-of-diffs plot supports only identical grid types (lat/lon or cubed-sphere)' + \
-              ' within each dataset pair')
-        raise ValueError
+        # Compute absolute and fractional differences
+        # Also make sure to preserve variable attributes
+        absdiffs = dev.copy()
+        fracdiffs = dev.copy()
+        for var in varlist:
+            absdiffs[var].values = dev[var].values - ref[var].values
+            fracdiffs[var].values = dev[var].values / ref[var].values
+            absdiffs[var].attrs = dev[var].attrs
+            fracdiffs[var].attrs = dev[var].attrs
+
+        # If there are more than one time point, then take the
+        # time mean (e.g. Annual Mean) of differences
+        if "time" in absdiffs.coords:
+            if len(absdiffs.coords["time"]) > 1:
+                absdiffs = dataset_mean(absdiffs)
+        if "time" in fracdiffs.coords:
+            if len(fracdiffs.coords["time"]) > 1:
+                fracdiffs = dataset_mean(fracdiffs)
 
     return absdiffs, fracdiffs