#!/usr/bin/env python ''' run_1yr_benchmark.py: Driver script for creating benchmark plots. Run this script to generate benchmark comparisons between: (1) GCC (aka GEOS-Chem "Classic") vs. GCC (2) GCHP vs GCC (3) GCHP vs GCHP (4) GCHP vs GCC diff-of-diffs You can customize this by editing the following settings in the "Configurables" section below: (1) Edit the path variables so that they point to folders w/ model data (2) Edit the version strings for each benchmark simulation (3) Edit the switches that turn on/off creating of plots and tables (4) If necessary, edit labels for the dev and ref versions Calling sequence: ./run_1yr_benchmark.py Remarks: By default, matplotlib will try to open an X window for plotting. If you are running this script in an environment where you do not have an active X display (such as in a computational queue), then you will need to use these commands to disable the X-window functionality. import os os.environ['QT_QPA_PLATFORM']='offscreen' For more information, please see this issue posted at the ipython site: https://github.com/ipython/ipython/issues/10627 This issue might be fixed in matplotlib 3.0. ''' # ===================================================================== # Imports and global settings (you should not need to edit these) # ===================================================================== from calendar import monthrange import os from os.path import join import xarray as xr from gcpy import benchmark as bmk from gcpy.core import get_filepaths, skip_these_vars import pandas as pd import numpy as np import warnings # Tell matplotlib not to look for an X-window os.environ['QT_QPA_PLATFORM']='offscreen' # Suppress harmless run-time warnings (mostly about underflow in division) warnings.filterwarnings('ignore', category=RuntimeWarning) warnings.filterwarnings('ignore', category=UserWarning) # ===================================================================== # Configurables # ===================================================================== # Benchmark information (*MUST EDIT*) # Note: When doing GCHP vs GCC comparisions gchp_dev will be compared to # gcc_dev (not gcc_ref!). This ensures consistency in version names when # doing GCHP vs GCC diff-of-diffs (mps, 6/27/19) maindir = '/n/scratchlfs02/jacob_lab/msulprizio/GC/benchmarks/1yr/12.7.0/TransportTracers/GCC' gcc_ref_version = 'Run0-StdWetDep' gcc_dev_version = 'Run1-LuoWetDep' gchp_ref_version = 'gchp_ref_version_string' gchp_dev_version = 'gchp_dev_version_string' # Comparisons to run (edit as needed) # NOTE: For now, just develop GCC vs GCC and add others later (bmy, 10/11/19) gcc_vs_gcc = True gchp_vs_gcc = False gchp_vs_gchp = False gchp_vs_gcc_diff_of_diffs = False # Output to generate (edit as needed) # Plots/tables will be created in this order: plot_conc = False plot_emis = False emis_table = False plot_jvalues = False plot_aod = False mass_table = True budget_table = False OH_metrics = False # Start and end of the benchmark year (as numpy.datetime64 dates) bmk_start = np.datetime64('2016-01-01') bmk_end = np.datetime64('2017-01-01') # Data directories (edit as needed) gcc_vs_gcc_refdir = join(maindir, gcc_ref_version, 'OutputDir') gcc_vs_gcc_devdir = join(maindir, gcc_dev_version, 'OutputDir') gchp_vs_gcc_refdir = join(maindir, gcc_dev_version, 'OutputDir') gchp_vs_gcc_devdir = join(maindir, gchp_dev_version, 'OutputDir') gchp_vs_gchp_refdir = join(maindir, gchp_ref_version, 'OutputDir') gchp_vs_gchp_devdir = join(maindir, gchp_dev_version, 'OutputDir') # Plots directories (edit as needed) gcc_vs_gcc_plotsdir = join(maindir, gcc_dev_version, 'Plots') gchp_vs_gchp_plotsdir = join(maindir, gchp_dev_version, 'Plots/GCHP_version_comparison') gchp_vs_gcc_plotsdir = join(maindir, gchp_dev_version, 'Plots/GCHP_GCC_comparison') diff_of_diffs_plotsdir = join(maindir, gchp_dev_version, 'Plots/GCHP_GCC_diff_of_diffs') # Plot title strings (edit as needed) # For gchp_vs_gcc_refstr use gcc_dev_version, not ref (mps, 6/27/19) gcc_vs_gcc_refstr = '{}'.format(gcc_ref_version) gcc_vs_gcc_devstr = '{}'.format(gcc_dev_version) gchp_vs_gcc_refstr = '{}'.format(gcc_dev_version) gchp_vs_gcc_devstr = '{}'.format(gchp_dev_version) gchp_vs_gchp_refstr = '{}'.format(gchp_ref_version) gchp_vs_gchp_devstr = '{}'.format(gchp_dev_version) diff_of_diffs_refstr = '{} - {}'.format(gcc_dev_version, gcc_ref_version) diff_of_diffs_devstr = '{} - {}'.format(gchp_dev_version, gchp_ref_version) # ===================================================================== # The rest of these settings should not need to be changed # ===================================================================== ############################################################################### # Under development, comment these out for now ## Paths to restart files #gcc_vs_gcc_refrst = core.get_gcc_filepath(gcc_vs_gcc_refdir, 'Restart', # year=bmk_year, month=bmk_months) #gcc_vs_gcc_devrst = core.get_gcc_filepath(gcc_vs_gcc_devdir, 'Restart', # year=bmk_year, month=bmk_months) # ## Paths to budget data #gcc_vs_gcc_devbgt = core.get_gcc_filepath(gcc_vs_gcc_devdir, 'Budget', # year=bmk_year, month=bmk_months) # ## Paths to concentration after chemistry files #gcc_vs_gcc_refcac = core.get_gcc_filepath(gcc_vs_gcc_refdir, 'ConcAfterChem', # year=bmk_year, month=bmk_months) # #gcc_vs_gcc_devcac = core.get_gcc_filepath(gcc_vs_gcc_devdir, 'ConcAfterChem', # year=bmk_year, month=bmk_months) ############################################################################### # Monthly array of dates bmk_delta_1m = np.timedelta64(1, 'M') bmk_months = np.arange(bmk_start, bmk_end, step=bmk_delta_1m, dtype='datetime64[M]') # Get the benchmark year from the datetime64 object bmk_year = bmk_months[0].astype('datetime64[Y]').astype(int) + 1970 # Seasonal array of dates bmk_delta_3m = np.timedelta64(3, 'M') bmk_seasons = np.arange(bmk_start, bmk_end, step=bmk_delta_3m, dtype='datetime64[M]') bmk_nseasons = len(bmk_seasons) # Names for each season (e.g. Jan2016, Apr2016, Jul2016, Oct2016) bmk_seasons_names = ['Jan', 'Apr', 'Jul', 'Oct'] bmk_seasons_names = [v + str(bmk_year) for v in bmk_seasons_names] # Seconds in each month of the benchmark year sec_per_month = np.zeros(12) for m in range(1, 13): month_info = monthrange(bmk_year, m) sec_per_month[m-1] = month_info[1] * 86400.0 # Files that will contain lists of quantities that have significant # differences -- we need these for the benchmark approval forms. gcc_vs_gcc_sigdiff = {} for mon_yr_str in bmk_seasons_names: vstr = '{}_vs_{}.{}'.format( gcc_vs_gcc_refstr, gcc_vs_gcc_devstr, mon_yr_str) sigdiff_files = [ join(gcc_vs_gcc_plotsdir, '{}.sig_diffs_sfc.txt'.format(vstr)), join(gcc_vs_gcc_plotsdir, '{}.sig_diffs_500hpa.txt'.format(vstr)), join(gcc_vs_gcc_plotsdir, '{}.sig_diffs_zonalmean.txt'.format(vstr)), join(gcc_vs_gcc_plotsdir, '{}.sig_diffs_emissions.txt'.format(vstr))] gcc_vs_gcc_sigdiff[mon_yr_str] = sigdiff_files ############################################################################## # For now we are only developing the GCC vs. GCC 1-yr benchmark code # so leave this commented out (bmy, 10/15/19) #vstr = '{}_vs_{}'{}.format(gchp_vs_gcc_refstr, gchp_vs_gcc_devstr) #gchp_vs_gcc_sigdiff = [ # join(gchp_vs_gcc_plotsdir, '{}_sig_diffs_sfc.txt'.format(vstr)), # join(gchp_vs_gcc_plotsdir, '{}_sig_diffs_500hpa.txt'.format(vstr)), # join(gchp_vs_gcc_plotsdir, '{}_sig_diffs_zonalmean.txt'.format(vstr)), # join(gchp_vs_gcc_plotsdir, '{}_sig_diffs_emissions.txt'.format(vstr))] # #vstr = '{}_vs_{}'.format(gchp_vs_gchp_refstr, gchp_vs_gchp_devstr) #gchp_vs_gchp_sigdiff = [ # join(gchp_vs_gchp_plotsdir, '{}_sig_diffs_sfc.txt'.format(vstr)), # join(gchp_vs_gchp_plotsdir, '{}_sig_diffs_500hpa.txt'.format(vstr)), # join(gchp_vs_gchp_plotsdir, '{}_sig_diffs_zonalmean.txt'.format(vstr)), # join(gchp_vs_gchp_plotsdir, '{}_sig_diffs_emissions.txt').format(vstr)] ############################################################################## # ====================================================================== # Create GCC vs GCC benchmark plots and tables # ====================================================================== if gcc_vs_gcc: if plot_conc: # -------------------------------------------------------------- # GCC vs GCC Concentration plots # (includes lumped species and separates by category) # -------------------------------------------------------------- print('\n%%% Creating GCC vs. GCC concentration plots %%%') # File lists for emissions data (seasonal) gcc_vs_gcc_refspc = get_filepaths(gcc_vs_gcc_refdir, 'SpeciesConc', bmk_seasons, is_gcc=True) gcc_vs_gcc_devspc = get_filepaths(gcc_vs_gcc_devdir, 'SpeciesConc', bmk_seasons, is_gcc=True) # Create seasonal concentration plots for s in range(bmk_nseasons): mon_yr_str = bmk_seasons_names[s] sigdiff_files= gcc_vs_gcc_sigdiff[mon_yr_str] bmk.make_benchmark_conc_plots(gcc_vs_gcc_refspc[s], gcc_vs_gcc_refstr, gcc_vs_gcc_devspc[s], gcc_vs_gcc_devstr, dst=gcc_vs_gcc_plotsdir, subdst=mon_yr_str, overwrite=True, sigdiff_files=sigdiff_files) if plot_emis: # -------------------------------------------------------------- # GCC vs GCC emissions plots # -------------------------------------------------------------- print('\n%%% Creating GCC vs. GCC emissions plots %%%') # File lists for emissions data (seasonal) gcc_vs_gcc_refhco = get_filepaths(gcc_vs_gcc_refdir, 'Emissions', bmk_seasons, is_gcc=True) gcc_vs_gcc_devhco = get_filepaths(gcc_vs_gcc_devdir, 'Emissions', bmk_seasons, is_gcc=True) # Create seasonal emissions plots for s in range(bmk_nseasons): mon_yr_str = bmk_seasons_names[s] sigdiff_files= gcc_vs_gcc_sigdiff[mon_yr_str] bmk.make_benchmark_emis_plots(gcc_vs_gcc_refhco[s], gcc_vs_gcc_refstr, gcc_vs_gcc_devhco[s], gcc_vs_gcc_devstr, dst=gcc_vs_gcc_plotsdir, subdst=mon_yr_str, plot_by_benchmark_cat=False, plot_by_hco_cat=False, overwrite=True, sigdiff_files=sigdiff_files) if emis_table: # -------------------------------------------------------------- # GCC vs GCC tables of emission and inventory totals # -------------------------------------------------------------- print('\n%%% Creating GCC vs. GCC emissions and inventory tables %%%') # File lists for J-values data (monthly) gcc_vs_gcc_refhco = get_filepaths(gcc_vs_gcc_refdir, 'Emissions', bmk_months, is_gcc=True) gcc_vs_gcc_devhco = get_filepaths(gcc_vs_gcc_devdir, 'Emissions', bmk_months, is_gcc=True) # Create emission tables bmk.make_benchmark_emis_tables(gcc_vs_gcc_refhco, gcc_vs_gcc_refstr, gcc_vs_gcc_devhco, gcc_vs_gcc_devstr, dst=gcc_vs_gcc_plotsdir, interval=sec_per_month, overwrite=True) if plot_jvalues: # -------------------------------------------------------------- # GCC vs GCC J-value plots # -------------------------------------------------------------- print('\n%%% Creating GCC vs. GCC J-value plots %%%') # Paths to J-value data (seasonal) gcc_vs_gcc_refjv = get_filepaths(gcc_vs_gcc_refdir, 'JValues', bmk_seasons, is_gcc=True) gcc_vs_gcc_devjv = get_filepaths(gcc_vs_gcc_devdir, 'JValues', bmk_seasons, is_gcc=True) # Create seasonal J-values plots for s in range(bmk_nseasons): mon_yr_str = bmk_seasons_names[s] sigdiff_files= gcc_vs_gcc_sigdiff[mon_yr_str] bmk.make_benchmark_jvalue_plots(gcc_vs_gcc_refjv[s], gcc_vs_gcc_refstr, gcc_vs_gcc_devjv[s], gcc_vs_gcc_devstr, dst=gcc_vs_gcc_plotsdir, subdst=mon_yr_str, overwrite=True, sigdiff_files=sigdiff_files) if plot_aod: # -------------------------------------------------------------- # GCC vs. GCC column AOD plots # -------------------------------------------------------------- print('\n%%% Creating GCC vs. GCC column AOD plots %%%') ## Paths to aerosol optical depth data gcc_vs_gcc_refaod = get_filepaths(gcc_vs_gcc_refdir, 'Aerosols', bmk_seasons, is_gcc=True) gcc_vs_gcc_devaod = get_filepaths(gcc_vs_gcc_devdir, 'Aerosols', bmk_seasons, is_gcc=True) # Create seasonal column AOD plots for s in range(bmk_nseasons): mon_yr_str = bmk_seasons_names[s] sigdiff_files= gcc_vs_gcc_sigdiff[mon_yr_str] bmk.make_benchmark_aod_plots(gcc_vs_gcc_refaod[s], gcc_vs_gcc_refstr, gcc_vs_gcc_devaod[s], gcc_vs_gcc_devstr, dst=gcc_vs_gcc_plotsdir, subdst=mon_yr_str, overwrite=True, sigdiff_files=sigdiff_files) if mass_table: # -------------------------------------------------------------- # GCC vs. GCC global mass tables # -------------------------------------------------------------- print('\n%%% Creating GCC vs. GCC global mass tables %%%') gcc_rstfile = 'restarts/GEOSChem.Restart.20170101_0000z.nc4' gcc_vs_gcc_refrst = join(maindir, gcc_ref_version, gcc_rstfile) gcc_vs_gcc_devrst = join(maindir, gcc_dev_version, gcc_rstfile) bmk.make_benchmark_mass_tables(gcc_vs_gcc_refrst, gcc_vs_gcc_refstr, gcc_vs_gcc_devrst, gcc_vs_gcc_devstr, dst=gcc_vs_gcc_plotsdir, overwrite=True) if budget_table: # -------------------------------------------------------------- # GCC vs GCC budgets tables # -------------------------------------------------------------- print('\n%%% Creating GCC vs. GCC budget tables %%%') bmk.make_benchmark_budget_tables(gcc_vs_gcc_devbgt, gcc_vs_gcc_devstr, dst=gcc_vs_gcc_plotsdir, overwrite=True) if OH_metrics: # -------------------------------------------------------------- # GCC vs GCC Global mean OH, MCF Lifetime, CH4 Lifetime # -------------------------------------------------------------- print('\n%%% Creating GCC vs. GCC OH metrics %%%') gcc_vs_gcc_reflist = [gcc_vs_gcc_refcac, gcc_vs_gcc_refmet] gcc_vs_gcc_devlist = [gcc_vs_gcc_devcac, gcc_vs_gcc_devmet] bmk.make_benchmark_oh_metrics(gcc_vs_gcc_reflist, gcc_vs_gcc_refstr, gcc_vs_gcc_devlist, gcc_vs_gcc_devstr, dst=gcc_vs_gcc_plotsdir, overwrite=True) ############################################################################### # NOTE: For now, we are developing the GCC vs. GCC 1-yr benchmark plots. # Leave this commented out for now. ## ===================================================================== ## Create GCHP vs GCC benchmark plots and tables ## ===================================================================== # #if gchp_vs_gcc: # if plot_conc: # # Concentration plots # # (includes lumped species and separates by category) # print('\n%%% Creating GCHP vs. GCC concentration plots %%%') # bmk.make_benchmark_conc_plots(gchp_vs_gcc_refspc, # gchp_vs_gcc_refstr, # gchp_vs_gcc_devspc, # gchp_vs_gcc_devstr, # dst=gchp_vs_gcc_plotsdir, # overwrite=True, # sigdiff_files=gchp_vs_gcc_sigdiff) # # if plot_emis: # # Emissions plots # print('\n%%% Creating GCHP vs. GCC emissions plots %%%') # bmk.make_benchmark_emis_plots(gchp_vs_gcc_refhco, # gchp_vs_gcc_refstr, # gchp_vs_gcc_devhco, # gchp_vs_gcc_devstr, # dst=gchp_vs_gcc_plotsdir, # plot_by_benchmark_cat=True, # plot_by_hco_cat=True, # overwrite=True, # flip_dev=True, # sigdiff_files=gchp_vs_gcc_sigdiff) # # if emis_table: # # Tables of emissions and inventory totals # print('\n%%% Creating GCHP vs. GCC emissions and inventory tables %%%') # gchp_vs_gcc_reflist = [gchp_vs_gcc_refhco] # gchp_vs_gcc_devlist = [gchp_vs_gcc_devhco, gchp_vs_gcc_devmet] # bmk.make_benchmark_emis_tables(gchp_vs_gcc_reflist, # gchp_vs_gcc_refstr, # gchp_vs_gcc_devlist, # gchp_vs_gcc_devstr, # dst=gchp_vs_gcc_plotsdir, # overwrite=True) # # if plot_jvalues: # # Local noon J-values plots # print('\n%%% Creating GCHP vs. GCC J-value plots %%%') # bmk.make_benchmark_jvalue_plots(gchp_vs_gcc_refjv, # gchp_vs_gcc_refstr, # gchp_vs_gcc_devjv, # gchp_vs_gcc_devstr, # dst=gchp_vs_gcc_plotsdir, # overwrite=True, # sigdiff_files=gchp_vs_gcc_sigdiff) # # if plot_aod: # # Column AOD plots # print('\n%%% Creating GCHP vs. GCC column AOD plots %%%') # bmk.make_benchmark_aod_plots(gchp_vs_gcc_refaod, # gchp_vs_gcc_refstr, # gchp_vs_gcc_devaod, # gchp_vs_gcc_devstr, # dst=gchp_vs_gcc_plotsdir, # overwrite=True, # sigdiff_files=gchp_vs_gcc_sigdiff) # ## Under development, leave commented out for now (bmy, 9/5/19) ## if mass_table: ## # Global mass tables ## print('\n%%% Creating GCHP vs. GCC global mass tables %%%') ## gchp_vs_gcc_reflist = [gchp_vs_gcc_refrst] ## gchp_vs_gcc_devlist = [gchp_vs_gcc_devrst, gchp_vs_gcc_devmetinst] ## bmk.make_benchmark_mass_tables(gchp_vs_gcc_reflist, ## gchp_vs_gcc_refstr, ## gchp_vs_gcc_devlist, ## gchp_vs_gcc_devstr, ## dst=gchp_vs_gcc_plotsdir, ## overwrite=True) # # if OH_metrics: # # Global mean OH, MCF Lifetime, CH4 Lifetime # print('\n%%% Creating GCHP vs. GCC OH metrics %%%') # gchp_vs_gcc_reflist = [gchp_vs_gcc_refcac, gchp_vs_gcc_refmet] # gchp_vs_gcc_devlist = [gchp_vs_gcc_devcac, gchp_vs_gcc_devmet] # bmk.make_benchmark_oh_metrics(gchp_vs_gcc_reflist, # gchp_vs_gcc_refstr, # gchp_vs_gcc_devlist, # gchp_vs_gcc_devstr, # dst=gchp_vs_gcc_plotsdir, # overwrite=True) # ## ===================================================================== ## Create GCHP vs GCHP benchmark plots and tables ## ===================================================================== # #if gchp_vs_gchp: # # if plot_conc: # # Concentration plots # # (includes lumped species and separates by category) # print('\n%%% Creating GCHP vs. GCHP concentration plots %%%') # bmk.make_benchmark_conc_plots(gchp_vs_gchp_refspc, # gchp_vs_gchp_refstr, # gchp_vs_gchp_devspc, # gchp_vs_gchp_devstr, # dst=gchp_vs_gchp_plotsdir, # overwrite=True, # sigdiff_files=gchp_vs_gchp_sigdiff) # # if plot_emis: # # Emissions plots # print('\n%%% Creating GCHP vs. GCHP emissions plots %%%') # bmk.make_benchmark_emis_plots(gchp_vs_gchp_refhco, # gchp_vs_gchp_refstr, # gchp_vs_gchp_devhco, # gchp_vs_gchp_devstr, # dst=gchp_vs_gchp_plotsdir, # plot_by_benchmark_cat=True, # plot_by_hco_cat=True, # overwrite=True, # flip_ref=True, # flip_dev=True, # sigdiff_files=gchp_vs_gchp_sigdiff) # # if emis_table: # # Tables of emissions and inventory totals # print('\n%%% Creating GCHP vs. GCHP emissions and inventory tables %%%') # gchp_vs_gchp_reflist = [gchp_vs_gchp_refhco, gchp_vs_gchp_refmet] # gchp_vs_gchp_devlist = [gchp_vs_gchp_devhco, gchp_vs_gchp_devmet] # bmk.make_benchmark_emis_tables(gchp_vs_gchp_reflist, # gchp_vs_gchp_refstr, # gchp_vs_gchp_devlist, # gchp_vs_gchp_devstr, # dst=gchp_vs_gchp_plotsdir, # overwrite=True) # # if plot_jvalues: # # Local noon J-values plots # print('\n%%% Creating GCHP vs. GCHP J-value plots %%%') # bmk.make_benchmark_jvalue_plots(gchp_vs_gchp_refjv, # gchp_vs_gchp_refstr, # gchp_vs_gchp_devjv, # gchp_vs_gchp_devstr, # dst=gchp_vs_gchp_plotsdir, # overwrite=True, # sigdiff_files=gchp_vs_gchp_sigdiff) # # if plot_aod: # # Column AOD plots # print('\n%%% Creating GCHP vs. GCHP column AOD plots %%%') # bmk.make_benchmark_aod_plots(gchp_vs_gchp_refaod, # gchp_vs_gchp_refstr, # gchp_vs_gchp_devaod, # gchp_vs_gchp_devstr, # dst=gchp_vs_gchp_plotsdir, # overwrite=True, # sigdiff_files=gchp_vs_gchp_sigdiff) # ## Under development, leave commented out for now (bmy, 9/5/19) ## if mass_table: ## # Global mass tables ## print('\n%%% Creating GCHP vs. GCHP global mass tables %%%') ## gchp_vs_gchp_reflist = [gchp_vs_gchp_refrst, gchp_vs_gchp_refmetinst] ## gchp_vs_gchp_devlist = [gchp_vs_gchp_devrst, gchp_vs_gchp_devmetinst] ## bmk.make_benchmark_mass_tables(gchp_vs_gchp_reflist, ## gchp_vs_gchp_refstr, ## gchp_vs_gchp_devlist, ## gchp_vs_gchp_devstr, ## dst=gchp_vs_gchp_plotsdir, ## overwrite=True) # # if OH_metrics: # # Global mean OH, MCF Lifetime, CH4 Lifetime # print('\n%%% Creating GCHP vs. GCHP OH metrics %%%') # gchp_vs_gcc_reflist = [gchp_vs_gchp_refcac, gchp_vs_gchp_refmet] # gchp_vs_gcc_devlist = [gchp_vs_gchp_devcac, gchp_vs_gchp_devmet] # bmk.make_benchmark_oh_metrics(gchp_vs_gchp_reflist, # gchp_vs_gchp_refstr, # gchp_vs_gchp_devlist, # gchp_vs_gchp_devstr, # dst=gchp_vs_gchp_plotsdir, # overwrite=True) # ## ===================================================================== ## Create GCHP vs GCC difference of differences benchmark plots ## ===================================================================== # #if gchp_vs_gcc_diff_of_diffs: # # # NOTE: This can be expanded to differences beyond species # # concentrations by following how this is done for conc plots. # print('%\n%% Creating GCHP vs. GCC diff-of-diffs concentration plots %%%') # # # Get a list of variables that GCPy should not read # skip_vars = core.skip_these_vars() # # # Target output files # diff_of_diffs_refspc = './gcc_diffs_spc.nc4' # diff_of_diffs_devspc = './gchp_diffs_spc.nc4' # # # Create a ref file that contains GCC differences # gcc_ref = xr.open_dataset(gcc_vs_gcc_refspc, drop_variables=skip_vars) # gcc_dev = xr.open_dataset(gcc_vs_gcc_devspc, drop_variables=skip_vars) # with xr.set_options(keep_attrs=True): # gcc_diffs = gcc_dev - gcc_ref # for v in gcc_dev.data_vars.keys(): # # Ensure the gcc_diffs Dataset includes attributes # gcc_diffs[v].attrs = gcc_dev[v].attrs # gcc_diffs.to_netcdf(diff_of_diffs_refspc) # # # Create a dev file that contains GCHP differences. Include special # # handling if cubed sphere grid dimension names are different since they # # changed in MAPL v1.0.0. # gchp_ref = xr.open_dataset(gchp_vs_gchp_refspc, drop_variables=skip_vars) # gchp_dev = xr.open_dataset(gchp_vs_gchp_devspc, drop_variables=skip_vars) # refdims = gchp_ref.dims # devdims = gchp_dev.dims # if 'lat' in refdims and 'Xdim' in devdims: # gchp_ref_newdimnames = gchp_dev.copy() # for v in gchp_dev.data_vars.keys(): # if 'Xdim' in gchp_dev[v].dims: # gchp_ref_newdimnames[v].values = gchp_ref[v].values.reshape( # gchp_dev[v].values.shape) # # NOTE: the reverse conversion is gchp_dev[v].stack(lat=('nf','Ydim')).transpose('time','lev','lat','Xdim').values # gchp_ref = gchp_ref_newdimnames.copy() # with xr.set_options(keep_attrs=True): # gchp_diffs = gchp_dev.copy() # for v in gchp_dev.data_vars.keys(): # if 'Xdim' in gchp_dev[v].dims or 'lat' in gchp_dev[v].dims: # gchp_diffs[v] = gchp_dev[v] - gchp_ref[v] # # NOTE: The gchp_diffs Dataset is created without variable # # attributes; we have to reattach them # gchp_diffs[v].attrs = gchp_dev[v].attrs # gchp_diffs.to_netcdf(diff_of_diffs_devspc) # # # # Create diff-of-diff plots for species concentrations # # (includes lumped species and separates by category) # # # # NOTE: Since at the present time we are only printing out # # diff-of-diffs for concentration plots, we can take this # # call out of the "if plot_conc:" block. # bmk.make_benchmark_conc_plots(diff_of_diffs_refspc, # diff_of_diffs_refstr, # diff_of_diffs_devspc, # diff_of_diffs_devstr, # dst=diff_of_diffs_plotsdir, # overwrite=True, # use_cmap_RdBu=True) # ###############################################################################