Now exclude "anchor" from open_dataset and open_mfdataset calls

According to xarray issues:
   pydata/xarray#3286
   pydata/xarray#1378

The open_mfdataset function has problems in creating a merged
dataset from multiple files in which variables have repeated
dimension names.  The easiest thing to do in this case is to
prevent such variables from being read in.

We now have added the drop_variables keyword to avoid reading
in the "anchor" variable in all calls to open_dataset and
open_mfdataset in both benchmark.py and core.py.  This variable is
only present in GCHP-created netCDF files using MAPL v1.0.0, which
is in GCHP 12.5.0 and later.

This commit should resolve GCPy issue #26:
   #26

Signed-off-by: Bob Yantosca <[email protected]>

Author: yantosca

gcpy/benchmark.py

@@ -32,6 +32,9 @@
 emission_spc = 'emission_species.json' 
 emission_inv = 'emission_inventories.json' 
 
+# List of variables that should not be read by xarray
+drop_varlist = ['anchor']
+
 
 def compare_single_level(refdata, refstr, devdata, devstr, varlist=None,
                          ilev=0, itime=0,  weightsdir=None, pdfname='',
@@ -2972,15 +2975,15 @@ def make_benchmark_conc_plots(ref, refstr, dev, devstr, dst='./1mo_benchmark',
 
     # Ref dataset
     try:
-        refds = xr.open_dataset(ref)
+        refds = xr.open_dataset(ref, drop_variables=dropvars)
     except FileNotFoundError:
         print('Could not find Ref file: {}'.format(ref))
         raise
     refds = core.add_lumped_species_to_dataset(refds, verbose=verbose)
 
     # Dev dataset
     try:
-        devds = xr.open_dataset(dev)
+        devds = xr.open_dataset(dev, drop_variables=dropvars)
     except FileNotFoundError:
         print('Could not find Dev file: {}!'.format(dev))
         raise
@@ -3213,14 +3216,14 @@ def make_benchmark_emis_plots(ref, refstr, dev, devstr,
 
     # Ref dataset
     try:
-        refds = xr.open_dataset(ref)
+        refds = xr.open_dataset(ref, drop_variables=dropvars)
     except FileNotFoundError:
         print('Could not find Ref file: {}'.format(ref))
         raise
 
     # Dev dataset
     try:
-        devds = xr.open_dataset(dev)
+        devds = xr.open_dataset(dev, drop_variables=dropvars)
     except FileNotFoundError:
         print('Could not find Dev file: {}'.format(dev))
         raise
@@ -3441,23 +3444,23 @@ def make_benchmark_emis_tables(reflist, refstr, devlist, devstr,
 
     # Ref
     if len(reflist) == 1:
-        refds = xr.open_dataset(reflist[0])
+        refds = xr.open_dataset(reflist[0], drop_variables=dropvars)
         assert gcc_area_name in list(refds.keys()),'Ref file {} does not contain area variable {}'.format(reflist[0], gcc_area_name)
 
     elif len(reflist) == 2:
-        refds = xr.open_dataset(reflist[0])
-        metrefds = xr.open_dataset(reflist[1])
+        refds = xr.open_dataset(reflist[0], drop_variables=dropvars)
+        metrefds = xr.open_dataset(reflist[1], drop_variables=dropvars)
         assert gchp_area_name in list(metrefds.keys()),'Ref met file {} does not contain area variable {}'.format(reflist[1], gchp_area_name)
         refds[gcc_area_name] = metrefds[gchp_area_name]
 
     # Dev
     if len(devlist) == 1:
-        devds = xr.open_dataset(devlist[0])
+        devds = xr.open_dataset(devlist[0], drop_variables=dropvars)
         assert gcc_area_name in list(refds.keys()),'Dev file {} does not contain area variable {}'.format(devlist[0], gcc_area_name)
 
     elif len(devlist) == 2:
-        devds = xr.open_dataset(devlist[0])
-        metdevds = xr.open_dataset(devlist[1])
+        devds = xr.open_dataset(devlist[0], drop_variables=dropvars)
+        metdevds = xr.open_dataset(devlist[1], drop_variables=dropvars)
         assert gchp_area_name in list(metdevds.keys()),'Dev met file {} does not contain area variable {}'.format(devlist[1], gchp_area_name)
         devds[gcc_area_name] = metdevds[gchp_area_name]
 
@@ -3604,14 +3607,14 @@ def make_benchmark_jvalue_plots(ref, refstr, dev, devstr,
 
     # Ref dataset
     try:
-        refds = xr.open_dataset(ref)
+        refds = xr.open_dataset(ref, drop_variables=dropvars)
     except FileNotFoundError:
         print('Could not find Ref file: {}'.format(ref))
         raise
 
     # Dev dataset
     try:
-        devds = xr.open_dataset(dev)
+        devds = xr.open_dataset(dev, drop_variables=dropvars)
     except FileNotFoundError:
         print('Could not find Dev file: {}'.format(dev))
         raise
@@ -3832,14 +3835,14 @@ def make_benchmark_aod_plots(ref, refstr, dev, devstr,
 
     # Read the Ref dataset
     try:
-        refds = xr.open_dataset(ref)
+        refds = xr.open_dataset(ref, drop_variables=dropvars)
     except FileNotFoundError:
         print('Could not find Ref file: {}'.format(ref))
         raise
 
     # Read the Dev dataset
     try:
-        devds = xr.open_dataset(dev)
+        devds = xr.open_dataset(dev, drop_variables=dropvars)
     except FileNotFoundError:
         print('Could not find Dev file: {}'.format(dev))
         raise
@@ -4054,14 +4057,14 @@ def make_benchmark_mass_tables(reflist, refstr, devlist, devstr,
 
     # Ref
     try:
-        refds = xr.open_mfdataset(reflist)
+        refds = xr.open_mfdataset(reflist, drop_variables=dropvars)
     except FileNotFoundError:
         print('Error opening Ref files: {}'.format(reflist))
         raise
 
     # Dev dataset
     try:
-        devds = xr.open_mfdataset(devlist)
+        devds = xr.open_mfdataset(devlist, drop_variables=dropvars)
     except FileNotFoundError:
         print('Error opening Dev files: {}!'.format(devlist))
         raise
@@ -4214,7 +4217,7 @@ def make_benchmark_budget_tables(devlist, devstr, dst='./1mo_benchmark',
 
     # Dev
     try:
-        devds = xr.open_mfdataset(devlist)
+        devds = xr.open_mfdataset(devlist, drop_variables=dropvars)
     except FileNotFoundError:
         print('Could not find one of the Dev files: {}'.format(devlist))
         raise

gcpy/core.py

@@ -17,6 +17,9 @@
 lumped_spc = 'lumped_species.json'
 bpch_to_nc_names = 'bpch_to_nc_names.json'
 
+# List of variables that should not be read by xarray
+dropvars = ['anchor']
+
 
 def open_dataset(filename, **kwargs):
     '''
@@ -74,7 +77,7 @@ def open_dataset(filename, **kwargs):
                          'pass a BPCH or netCDF file with extension '
                          '"bpch" or "nc"!'.format(file_extension))
 
-    return _opener(filename, **kwargs)
+    return _opener(filename, **kwargs, drop_variables=dropvars)
 
 
 def open_mfdataset(filenames, concat_dim='time', compat='no_conflicts',
@@ -156,19 +159,23 @@ def open_mfdataset(filenames, concat_dim='time', compat='no_conflicts',
                          '"bpch" or "nc" or "nc4"'.format(file_extension))
 
     return _opener(filenames, concat_dim=concat_dim, compat=compat,
-                   preprocess=preprocess, lock=lock, **kwargs)
+                   preprocess=preprocess, lock=lock,
+                   drop_variables=dropvars, **kwargs)
 
 
 def get_gcc_filepath(outputdir, collection, day, time):
     if collection == 'Emissions':
-        filepath = os.path.join(outputdir, 'HEMCO_diagnostics.{}{}.nc'.format(day,time))
+        filepath = os.path.join(outputdir,
+                                'HEMCO_diagnostics.{}{}.nc'.format(day,time))
     else:
-        filepath = os.path.join(outputdir, 'GEOSChem.{}.{}_{}z.nc4'.format(collection,day,time))
+        filepath = os.path.join(outputdir,
+                                'GEOSChem.{}.{}_{}z.nc4'.format(collection,day,time))
     return filepath
 
 
 def get_gchp_filepath(outputdir, collection, day, time):
-    filepath = os.path.join(outputdir, 'GCHP.{}.{}_{}z.nc4'.format(collection,day,time))
+    filepath = os.path.join(outputdir,
+                            'GCHP.{}.{}_{}z.nc4'.format(collection,day,time))
     return filepath