<jupyter_start><jupyter_text>Compare two GEOS-Chem diagnostic datasetsIncludes GEOS-Chem Classic and GCHP, binary and netcdf formats. Datasets on different grids and at different grid resolutions can be compared directly following the example in this notebook. Plots may be displayed inline in this notebook or sent to PDF.  Tips for using Jupyter Notebook1. You can use unix commands such as ls and pwd2. Hover over a function and press Shift-Tab to see function arguments3. Access documentation with help(fnct) or fnct?, where fnct is function name4. See function source code with fnct??5. Use tab-complete to see what is possible, e.g. gcpy.cTAB6. Explore modules with dir(module), e.g. dir(gcpy) or dir(gcpy.core)  Import dependencies<jupyter_code>import os
import numpy as np
import xarray as xr
import xbpch
import gcpy

# Enable auto-reloading modules
%load_ext autoreload
%autoreload 2

%matplotlib inline
import warnings; warnings.filterwarnings("ignore")<jupyter_output><empty_output><jupyter_text>Choose comparison type<jupyter_code># Choose comparison type: 
#   'bpch' = bpch vs bpch
#   'nc' = nc vs nc
#.  'bpch_vs_nc' = bpch vs nc
# NOTE: The example in the existing settings of this notebook is for bpch vs nc

#comparison = 'bpch'
#comparison = 'nc'
comparison = 'bpch_vs_nc'<jupyter_output><empty_output><jupyter_text>Load data  Binary data (if using)<jupyter_code># Set bpch ref data is doing bpch vs nc or bpch vs bpch
if 'bpch' in comparison:
    refdir = '/path/to/ref/binary/file/'
    refbpfile = os.path.join(refdir,'trac_avg.geosfp_4x5_standard.201607010000')
    reftracerinfo = os.path.join(refdir,'tracerinfo.dat') 
    refdiaginfo = os.path.join(refdir,'diaginfo.dat')
    refdata = xbpch.open_bpchdataset(refbpfile, tracerinfo_file=reftracerinfo, diaginfo_file=refdiaginfo)
    
# Set bpch dev data if doing bpch vs bpch
if comparison == 'bpch':
    devdir = '/path/to/dev/binary/file'
    devbpfile = os.path.join(devdir,'trac_avg.geosfp_4x5_standard.201607010000')
    devtracerinfo = os.path.join(devdir,'tracerinfo.dat') 
    devdiaginfo = os.path.join(devdir,'diaginfo.dat')
    devdata = xbpch.open_bpchdataset(devbpfile, tracerinfo_file=reftracerinfo, diaginfo_file=refdiaginfo)<jupyter_output><empty_output><jupyter_text>Netcdf data (if using)<jupyter_code># Set netcdf ref data if doing nc vs nc
if comparison == 'nc':
    refdata = xr.open_dataset('/path/to/ref/netcdf/file')
    
# Set netcdf dev data if doing bpch vs nc
if 'nc' in comparison:
    devdata = xr.open_dataset('/path/to/dev/netcdf/file')<jupyter_output><empty_output><jupyter_text>Inspect datasets<jupyter_code>refdata
devdata<jupyter_output><empty_output><jupyter_text>If using binary data, process it further<jupyter_code># Reorder the dimensions in the array to always be time, lev, lat, lon (so same order as netcdf)
# NOTE: check dimensions in refdata in previous step. Make sure they are all here and there are no extra.
if 'bpch' in comparison:
    refdata = refdata.transpose('lev','lev_edge','lat', 'lon') 
if comparison == 'bpch':
    devdata = devdata.transpose('lev','lev_edge','lat', 'lon')     
# convert binary names to analogous netcdf diagnostic names
# NOTE: not necessary if comparing bpch vs bpch and not doing an official GEOS-Chem benchmark
if 'bpch' in comparison:
    refdata = gcpy.convert_bpch_names_to_netcdf_names(refdata)    
if comparison == 'bpch':
    devdata = gcpy.convert_bpch_names_to_netcdf_names(devdata)
# Change the units for concentrations. If ppb or ppbv then change to mol/mol. Will need to do this for
# other binary diagnostic units if they are not the same as netcdf. Ideally put the mapping in a dictionary.
conc_units = ['ppb','ppbv']
if 'bpch' in comparison:
    for var in refdata.data_vars:
        if refdata[var].units.strip() in conc_units:
            refdata[var].attrs['units'] = 'mol mol-1 dry'
            refdata[var].values = refdata[var].values * 1e-9
if comparison == 'bpch':
    for var in refdata.data_vars:
        if devdata[var].units.strip() in conc_units:
            devdata[var].attrs['units'] = 'mol mol-1 dry'
            devdata[var].values = devdata[var].values * 1e-9<jupyter_output><empty_output><jupyter_text>Inspect common variables in ref and devIf you are comparing binary and netcdf diagnostics then you can expect a large number of variables with different names.<jupyter_code>[cmnvars, cmnvarsOther, cmnvars2D, cmnvars3D] = gcpy.compare_varnames(refdata, devdata)<jupyter_output>161 common variables
686 variables in ref only (skip)
   Variable names: ['ProdSO2fromDMSandOH', 'ProdSO2fromDMSandNO3', 'ProdSO2fromDMS', 'ProdMSAfromDMS', 'ProdSO4fromGasPhase', 'ProdSO4fromH2O2inCloud', 'ProdSO4fromO3inCloud', 'PL_SUL_S_SO4o2', 'ProdSO4fromO3inSeaSalt', 'PL_SUL_S_LOH', 'PL_SUL_S_LNO3', 'ProdSO4fromOxidationOnDust', 'ProdNITfromHNO3uptakeOnDust', 'ProdSO4fromUptakeOfH2SO4g', 'LossHNO3onSeaSalt', 'ProdSO4fromHOBrInCloud', 'ProdSO4fromSRO3', 'ProdSO4fromSRHObr', 'ProdSO4fromO3s', 'PL_SUL_S_HPLUS', 'EmisDST1_Natural', 'EmisDST2_Natural', 'EmisDST3_Natural', 'EmisDST4_Natural', 'BC_ANTH_BLKC', 'BC_BIOF_BLKC', 'BC_BIOB_BLKC', 'EmisOC_Anthro', 'OC_BIOF_ORGC', 'OC_BIOB_ORGC', 'ProdBCPIfromBCPO', 'OC_BIOG_ORGC', 'ProdOCPIfromOCPO', 'EmisSALA_Natural', 'EmisSALC_Natural', 'EmisACET_Monoterp', 'EmisACET_MethylBut', 'EmisACET_DirectBio', 'EmisACET_Ocean', 'EmisDMS_Ocean', 'EmisSO2_Aircraft', 'SO2_AN_S_SO2', 'SO2_BIOB_SO2', 'SO2_BIOF_SO2', 'EmisSO2_EVOL', 'EmisSO2_NVOL', 'EmisSO[...]<jupyter_text>Set and ref and dev strings for including in plots<jupyter_code>refstr='GCC binary, 12.2.0' # EDIT AS NEEDED
devstr='GCHP netcdf, 12.2.0' # EDIT AS NEEDED
    
print('Ref: {}'.format(refstr))
print('Dev: {}'.format(devstr))<jupyter_output>Ref: GCC binary, 12.2.0
Dev: GCHP netcdf, 12.2.0<jupyter_text>Compare species concentrations  Define list of variables to plot<jupyter_code>#------------ Below are a few examples ---------------

# Example 1: all species
#desc='spc_conc'
#varlist = [k for k in cmnvars if 'SpeciesConc_' in k]

# Example 2: Subset of species
#desc='Dust'
#varlist = [k for k in cmnvars if 'DST' in k and 'SpeciesConc' in k]

# Example 3: one species
desc='Ozone'
varlist = ['SpeciesConc_O3']

print('{} variables found'.format(len(varlist)))
varlist<jupyter_output>1 variables found<jupyter_text>Inspect units to make sure they match (they might not if binary versus netcdf)<jupyter_code>for var in varlist:
    print(var)
    print('   ref units: {}'.format(refdata[var].units))
    print('   dev units: {}'.format(devdata[var].units))
# If units do not match, the following code may be helpful:
#
# Print all unique units in bpch datasets
# unitstrs = []
# for var in refdata.data_vars:
#     if refdata[var].units.strip() not in unitstrs:
#         unitstrs.append(refdata[var].units.strip())
# unitstrs<jupyter_output><empty_output><jupyter_text>Plot inline in notebookIf you use this function in ipython then import matplotlib.pyplot as plt and do plt.show() to create pop-up with figure.<jupyter_code># Try hovering over function name and pressing shift-tab to see optional arguments lst 
# NOTE: cmpres argument is only for when ref and dev grid resolutions differ
gcpy.compare_single_level( refdata, refstr, devdata, devstr, cmpres='1x1.25', varlist=['SpeciesConc_DST3'])
# Zonal mean - quick plot
gcpy.compare_zonal_mean( refdata, refstr, devdata, devstr, varlist=['SpeciesConc_DST3'])<jupyter_output>Reuse existing file: ./conservative_4x5_1x1.25.nc
Reuse existing file: ./conservative_c24_1x1.25_0.nc
Reuse existing file: ./conservative_c24_1x1.25_1.nc
Reuse existing file: ./conservative_c24_1x1.25_2.nc
Reuse existing file: ./conservative_c24_1x1.25_3.nc
Reuse existing file: ./conservative_c24_1x1.25_4.nc
Reuse existing file: ./conservative_c24_1x1.25_5.nc<jupyter_text>Send plots to PDFIf you want to a create a pdf with more than one species then redefine varlist prior to passing to plot functions.<jupyter_code># This directory must already exist
plotsdir = '/path/to/plots/directory'
# Surface - PDF
pdfname = os.path.join(plotsdir,'{}_surface_{}_diag.pdf'.format(desc,comparison))
gcpy.compare_single_level( refdata, refstr, devdata, devstr, varlist=varlist, pdfname=pdfname )
gcpy.add_bookmarks_to_pdf( pdfname, varlist, remove_prefix='SpeciesConc_')
# Zonal mean - PDF
pdfname = os.path.join(plotsdir,'{}_zonalmean_{}_diag.pdf'.format(desc, comparison))
gcpy.compare_zonal_mean( refdata, refstr, devdata, devstr, varlist=varlist, pdfname=pdfname )
gcpy.add_bookmarks_to_pdf( pdfname, varlist, remove_prefix='SpeciesConc_')<jupyter_output>Reuse existing file: ./conservative_4x5_1x1.25.nc
Reuse existing file: ./conservative_c24_1x1.25_0.nc
Reuse existing file: ./conservative_c24_1x1.25_1.nc
Reuse existing file: ./conservative_c24_1x1.25_2.nc
Reuse existing file: ./conservative_c24_1x1.25_3.nc
Reuse existing file: ./conservative_c24_1x1.25_4.nc
Reuse existing file: ./conservative_c24_1x1.25_5.nc

Creating /Users/lizzielundgren/gc/data/12.2.0_test/plots/Ozone_zonalmean_bpch_vs_nc_diag.pdf for 1 variables
0
{'fraction_non_alphanumeric': 0.10920747209019564, 'fraction_numerical': 0.024096385542168676, 'mean_word_length': 5.542299349240781, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 44, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '474118', 'license': 'NCSA', 'n_tokens_mistral': 3446, 'n_tokens_neox': 3266, 'n_words': 850, 'path': 'examples/diagnostics/compare_diagnostics.ipynb', 'repo_name': 'sdeastham/gcpy', 'size': '9047'}