gchp_vs_gcc_comparison.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare GCHP and GCC NetCDF Diagnostics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview of this Notebook\n",
    "\n",
    "* Import dependencies\n",
    "* Define functions: utility, regridding, and plotting (will eventually be in gcpy)\n",
    "* Define data sources\n",
    "* Compare restart files\n",
    "* Store grid areas for normalization\n",
    "* Interactively inspect species concentrations\n",
    "* Create PDFs comparing other diagnostic collections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import dependencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import cubedsphere as cs\n",
    "import xesmf as xe\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "from matplotlib import ticker\n",
    "from matplotlib.backends.backend_pdf import PdfPages\n",
    "from matplotlib.colors import ListedColormap\n",
    "from cartopy import crs\n",
    "from cartopy.mpl.geoaxes import GeoAxes\n",
    "import random\n",
    "\n",
    "%matplotlib inline\n",
    "import warnings; warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "# Also define colormap. Colormap file source: https://bitbucket.org/gcst/gcpy\n",
    "rgb_WhGrYlRd = np.genfromtxt('/n/home08/elundgren/GC/python/WhGrYlRd.txt',delimiter=' ')\n",
    "WhGrYlRd = ListedColormap(rgb_WhGrYlRd/255.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define general utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_gcc_filepath(outputdir, collection, day, time):\n",
    "    if collection == 'Emissions':\n",
    "        filepath = os.path.join(outputdir, 'HEMCO_diagnostics.{}{}.nc'.format(day,time))\n",
    "    else:\n",
    "        filepath = os.path.join(outputdir, 'GEOSChem.{}.{}_{}z.nc4'.format(collection,day,time))\n",
    "    return filepath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_gchp_filepath(outputdir, collection, day, time):\n",
    "    filepath = os.path.join(outputdir, 'GCHP.{}.{}_{}z.nc4'.format(collection,day,time))\n",
    "    return filepath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_paths( gcc_path, gchp_path):\n",
    "    if not os.path.exists(gcc_path):\n",
    "        print('ERROR! GCC path not exist: {}'.format(gcc_path))\n",
    "    else:\n",
    "        print('GCC path exists: {}'.format(gcc_path))\n",
    "    if not os.path.exists(gchp_path):\n",
    "        print('ERROR! GCHP file does not exist: {}'.format(gchp_path))\n",
    "    else:\n",
    "        print('GCHP path exists: {}'.format(gchp_path))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def compare_varnames(dgcc, dgchp):\n",
    "    \n",
    "    # Find common variables in collection by generating lists and list overlap\n",
    "    gchpvars = [k for k in dgchp.data_vars.keys()]\n",
    "    gccvars= [k for k in dgcc.data_vars.keys()]\n",
    "    commonvars = sorted(list(set(gchpvars).intersection(set(gccvars))))\n",
    "    gcconly = [v for v in gccvars if v not in gchpvars]\n",
    "    gchponly = [v for v in gchpvars if v not in gccvars]\n",
    "    dimmismatch = [v for v in commonvars if dgcc[v].ndim != dgchp[v].ndim]\n",
    "    commonvars2D = [v for v in commonvars if dgcc[v].ndim == 3]\n",
    "    commonvars3D = [v for v in commonvars if dgcc[v].ndim == 4]\n",
    "    \n",
    "    # Print information on common and mismatching variables, as well as dimensions\n",
    "    print('{} common variables ({} are 2-dim and {} are 3-dim)'.format(len(commonvars), len(commonvars2D), len(commonvars3D)))\n",
    "    if len(gcconly) > 0:\n",
    "        print('{} variables in gcc only (skip)'.format(len(gcconly)))\n",
    "        print('   Variable names: {}'.format(gcconly))\n",
    "    if len(gchponly) > 0:\n",
    "        print('{} variables in gchp only (skip)'.format(len(gchponly)))\n",
    "        print('   Variable names: {}'.format(gchponly))\n",
    "    if len(dimmismatch) > 0:\n",
    "        print('{} common variables have different dimensions'.format(len(dimmismatch)))\n",
    "        print('   Variable names: {}'.format(dimmismatch))\n",
    "        \n",
    "    return [commonvars, commonvars2D, commonvars3D]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_stats(dgcc, dgchp, varname):\n",
    "    ds_gcc_var = dgcc[varname]\n",
    "    ds_gchp_var = dgchp[varname]\n",
    "    print('Variable array type: {}'.format(type(ds_gcc_var)))\n",
    "    print('Variable values array type: {}'.format(type(ds_gcc_var.values)))\n",
    "    print('Array shapes:')\n",
    "    print('    GCC:  {}'.format(ds_gcc_var.shape))\n",
    "    print('    GCHP: {}'.format(ds_gchp_var.shape))\n",
    "    print('Global Mean:')\n",
    "    print('    GCC:  {}'.format(np.round(ds_gcc_var.values.mean(),20)))\n",
    "    print('    GCHP: {}'.format(np.round(ds_gchp_var.values.mean(),20)))\n",
    "    print('Global Min:')\n",
    "    print('    GCC:  {}'.format(np.round(ds_gcc_var.values.min(),20)))\n",
    "    print('    GCHP: {}'.format(np.round(ds_gchp_var.values.min(),20)))\n",
    "    print('Global Max:')\n",
    "    print('    GCC:  {}'.format(np.round(ds_gcc_var.values.max(),20)))\n",
    "    print('    GCHP: {}'.format(np.round(ds_gchp_var.values.max(),20)))\n",
    "    print('Global Sum:')\n",
    "    print('    GCC:  {}'.format(np.round(ds_gcc_var.values.sum(),20)))\n",
    "    print('    GCHP: {}'.format(np.round(ds_gchp_var.values.sum(),20)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_n_random_vars(varlist,N):\n",
    "    sublist = [varlist[i] for i in sorted(random.sample(range(len(varlist)), N))]\n",
    "    return sublist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define functions for regridding\n",
    "\n",
    "These functions utilize the xESMF package and the cubedsphere package developed by Jaiwei Zhuang, a graduate student at Harvard University. More specifically I use the cubedsphere package CSGrid class and the xESMF package Regridder class. The latter serves as a wrapper for the ESMF/ESMPy packages. These tools enable simple conservative regridding between lat/lon and cubed sphere horizontal grids as specifically defined by the NASA Global Modeling and Assimilation Office (GMAO). See http://xesmf.readthedocs.io/en/latest/ and https://github.com/JiaweiZhuang/cubedsphere for more information. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_grid_LL(llres):\n",
    "    [dlat,dlon] = list(map(float, llres.split('x')))\n",
    "    lon_b = np.linspace(-180 - dlon/2, 180 - dlon/2, int(360/dlon) + 1, endpoint=True)\n",
    "    lat_b = np.linspace(-90 - dlat/2, 90 + dlat/2, int(180/dlat) + 2, endpoint=True).clip(-90,90)\n",
    "    lat = (lat_b[1:] + lat_b[:-1]) / 2\n",
    "    lon = (lon_b[1:] + lon_b[:-1]) / 2\n",
    "    llgrid = {'lat': lat, \n",
    "              'lon': lon, \n",
    "              'lat_b': lat_b, \n",
    "              'lon_b': lon_b}\n",
    "    return llgrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_grid_CS(csres):\n",
    "    csgrid = cs.csgrid_GMAO(csres)\n",
    "    csgrid_list = [None]*6\n",
    "    for i in range(6):\n",
    "        csgrid_list[i] = {'lat': csgrid['lat'][i], \n",
    "                          'lon': csgrid['lon'][i],\n",
    "                          'lat_b': csgrid['lat_b'][i], \n",
    "                          'lon_b': csgrid['lon_b'][i]}\n",
    "    return [csgrid, csgrid_list]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_regridder_C2L( csres_in, llres_out, weightsdir='.', reuse_weights=False,  ):\n",
    "    csgrid, csgrid_list = make_grid_CS(csres_in)\n",
    "    llgrid = make_grid_LL(llres_out)\n",
    "    regridder_list = []\n",
    "    for i in range(6):\n",
    "        weightsfile = os.path.join(weightsdir, 'conservative_c{}_{}_{}.nc'.format(str(csres_in), llres_out, str(i)))\n",
    "        regridder = xe.Regridder(csgrid_list[i], llgrid, method='conservative', filename=weightsfile, reuse_weights=reuse_weights)\n",
    "        regridder_list.append(regridder)\n",
    "    return regridder_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_regridder_L2L( llres_in, llres_out, weightsdir='.', reuse_weights=False ):\n",
    "    llgrid_in = make_grid_LL(llres_in)\n",
    "    llgrid_out = make_grid_LL(llres_out)\n",
    "    weightsfile = os.path.join(weightsdir,'conservative_{}_{}.nc'.format(llres_in, llres_out))\n",
    "    regridder = xe.Regridder(llgrid_in, llgrid_out, method='conservative', filename=weightsfile, reuse_weights=reuse_weights)\n",
    "    return regridder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define functions for plotting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function to plot and compare single vertical level of GCHP and GCC\n",
    "\n",
    "Notes: This function regrids both GCC and GCHP on-the-fly using the pre-generated weights. You can use it to plot interactively or to generate a multi-page pdf of plots. It takes a list of variable names to plot for a single collection only. You can plot for any level and any time slice in the file. You can also input GCHP and GCC output at any resolution and it will be regridded to a common lat-lon resolution for the comparison, showing plots for both native resolutions as well as the regridded. Some diagnostics, such as wet loss, need to be normalized by area for a comparison that makes sense. You can do this by passing the appropriate arguments. If comparing restart files you must flip the veritical axis of the GCHP restart file, also configurable by passing an argument. This is not necessary if comparing output diagnostic files. By default the colorbars for the GCHP and GCC plots will have the same range, but you can turn this feature off. Also by default the colorbar of the fractional difference between the model outputs will be limited to +/-2, but you can change this as well via the passed parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_single_level(dgcc, dgchp, varlist=None, weightsdir='.', ilev=0, itime=0, llres_raw='4x5', \n",
    "                         llres_cmp='1x1.25', savepdf=False, pdfname='gchp_vs_gcc_map.pdf', \n",
    "                         match_cbar=True, full_ratio_range=False, normalize_by_area=False, \n",
    "                         area1=None, area2=None, check_units=True, flip_vert=False):\n",
    "    \n",
    "    # If no varlist is passed, plot all (surface only for 3D)\n",
    "    if varlist == None:\n",
    "        [varlist, commonvars2D, commonvars3D] = compare_varnames(dgcc, dgchp)\n",
    "        print('Plotting all common variables (surface only if 3D)')\n",
    "    n_var = len(varlist)\n",
    "    \n",
    "    # Get lat-lon grids and regridder. Assume regridding weights have already been generated\n",
    "    llgrid_raw = make_grid_LL(llres_raw)\n",
    "    llgrid_cmp = make_grid_LL(llres_cmp)\n",
    "    ll_regridder = make_regridder_L2L(llres_raw, llres_cmp, weightsdir=weightsdir, reuse_weights=True)\n",
    "\n",
    "    # Get cubed sphere grid and regridder\n",
    "    csres = dgchp['lon'].size\n",
    "    [csgrid, csgrid_list] = make_grid_CS(csres)\n",
    "    cs_regridder_list = make_regridder_C2L(csres, llres_cmp, weightsdir=weightsdir, reuse_weights=True)\n",
    "\n",
    "    # Get lat/lon extents\n",
    "    [minlon_raw, maxlon_raw] = [min(llgrid_raw['lon_b']), max(llgrid_raw['lon_b'])]\n",
    "    [minlat_raw, maxlat_raw] = [min(llgrid_raw['lat_b']), max(llgrid_raw['lat_b'])]\n",
    "    [minlon_cmp, maxlon_cmp] = [min(llgrid_cmp['lon_b']), max(llgrid_cmp['lon_b'])]\n",
    "    [minlat_cmp, maxlat_cmp] = [min(llgrid_cmp['lat_b']), max(llgrid_cmp['lat_b'])]\n",
    "\n",
    "    # Create pdf (if saving)\n",
    "    if savepdf:\n",
    "        print('\\nCreating {} for {} variables'.format(pdfname,n_var))\n",
    "        pdf = PdfPages(pdfname)\n",
    "\n",
    "    # Loop over variables\n",
    "    for ivar in range(n_var):\n",
    "        if savepdf: print('{} '.format(ivar), end='')\n",
    "        varname = varlist[ivar]\n",
    "        \n",
    "        # Do some checks: dimensions and units\n",
    "        varndim = dgchp[varname].ndim\n",
    "        varndim2 = dgcc[varname].ndim\n",
    "        if check_units: assert varndim == varndim2, 'GCHP and GCC dimensions do not agree for {}!'.format(varname)\n",
    "        units_raw = dgchp[varname].units\n",
    "        units2 = dgcc[varname].units\n",
    "        if check_units: assert units_raw == units2, 'GCHP and GCC units do not match for {}!'.format(varname)\n",
    "            \n",
    "        # if normalizing by area, adjust units to be per m2, and adjust title string\n",
    "        units = units_raw\n",
    "        subtitle_extra = ''\n",
    "                 \n",
    "        # Slice the data\n",
    "        if varndim == 4: \n",
    "            if 'ilev' in dgcc[varname].dims:\n",
    "                ds1 = dgcc[varname].isel(time=itime,ilev=ilev)\n",
    "            else:\n",
    "                ds1 = dgcc[varname].isel(time=itime,lev=ilev)\n",
    "            if flip_vert: ds2 = dgchp[varname].isel(time=itime,lev=71-ilev)\n",
    "            else: ds2 = dgchp[varname].isel(time=itime,lev=ilev)\n",
    "        elif varndim == 3: \n",
    "            ds1 = dgcc[varname].isel(time=itime)\n",
    "            ds2 = dgchp[varname].isel(time=itime)\n",
    "            \n",
    "        # if normalizing by area, transform on the native grid and adjust units and subtitle string\n",
    "        exclude_list = ['WetLossConvFrac','Prod_','Loss_']\n",
    "        if normalize_by_area and not any(s in varname for s in exclude_list):\n",
    "            ds1.values = ds1.values / area1\n",
    "            ds2.values = ds2.values / area2\n",
    "            units = '{} m-2'.format(units_raw)\n",
    "            subtitle_extra = ', Normalized by Area'\n",
    "            \n",
    "        # Regrid the slices\n",
    "        csdata = ds2.data.reshape(6,csres,csres)\n",
    "        gchp_ll = np.zeros([llgrid_cmp['lat'].size, llgrid_cmp['lon'].size])\n",
    "        for i in range(6):\n",
    "            regridder = cs_regridder_list[i]\n",
    "            gchp_ll += regridder(csdata[i])\n",
    "        gcc_ll = ll_regridder(ds1)\n",
    "        \n",
    "        # Get min and max for colorbar limits\n",
    "        vmin_gchp = np.min([csdata.min(), gchp_ll.min()])\n",
    "        vmin_gcc = np.min([ds1.values.min(), gcc_ll.values.min()])\n",
    "        vmin_cmn = np.min([vmin_gchp, vmin_gcc])\n",
    "        vmax_gchp = np.max([csdata.max(), gchp_ll.max()])\n",
    "        vmax_gcc = np.max([ds1.values.max(), gcc_ll.values.max()])\n",
    "        vmax_cmn = np.max([vmax_gchp, vmax_gcc])\n",
    "        if match_cbar: [vmin, vmax] = [vmin_cmn, vmax_cmn]\n",
    "        \n",
    "        # Create 3x2 figure\n",
    "        figs, ((ax0, ax1), (ax2, ax3), (ax4, ax5)) = plt.subplots(3, 2, \n",
    "                                                                  figsize=[12,14], \n",
    "                                                                  subplot_kw={'projection': crs.PlateCarree()})\n",
    "        # Give the figure a title\n",
    "        offset = 0.96\n",
    "        fontsize=25\n",
    "        if varndim == 4:\n",
    "            if ilev == 0: levstr = 'Surface'\n",
    "            elif ilev == 22: levstr = '500 hPa'\n",
    "            else: levstr = 'Level ' +  str(ilev-1)\n",
    "            figs.suptitle('{}, {}'.format(varname,levstr), fontsize=fontsize, y=offset)\n",
    "        elif varndim == 3: \n",
    "            figs.suptitle('{}'.format(varname), fontsize=fontsize, y=offset)\n",
    "        else:\n",
    "            print('varndim is 2 for {}! Must be 3 or 4.'.format(varname))\n",
    "            \n",
    "        # Set bounds of auto-tick range\n",
    "        \n",
    "            \n",
    "        # Subplot (0,0): GCHP raw\n",
    "        ax0.coastlines()\n",
    "        if not match_cbar: [vmin, vmax] = [vmin_gchp, vmax_gchp]        \n",
    "        masked_csdata = np.ma.masked_where(np.abs(csgrid['lon'] - 180) < 2, csdata) # based on cubedsphere plotCS_quick_raw\n",
    "        for i in range(6):\n",
    "            plot0 = ax0.pcolormesh(csgrid['lon_b'][i,:,:], csgrid['lat_b'][i,:,:], masked_csdata[i,:,:], \n",
    "                                   cmap=WhGrYlRd,vmin=vmin, vmax=vmax)\n",
    "        ax0.set_title('GCHP Raw{}\\nC{}'.format(subtitle_extra,str(csres)))\n",
    "        cb = plt.colorbar(plot0, ax=ax0, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.1 or (vmax-vmin) > 100:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)            \n",
    "        \n",
    "        # Subplot (0,1): GCHP regridded\n",
    "        ax1.coastlines()\n",
    "        if not match_cbar: [vmin, vmax] = [vmin_gchp, vmax_gchp]\n",
    "        plot1 = ax1.imshow(gchp_ll, extent=(minlon_cmp, maxlon_cmp, minlat_cmp, maxlat_cmp), \n",
    "                           cmap=WhGrYlRd,vmin=vmin, vmax=vmax)\n",
    "        ax1.set_title('GCHP Regridded\\n{}'.format(llres_cmp))\n",
    "        cb = plt.colorbar(plot1, ax=ax1, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.1 or (vmax-vmin) > 100:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)\n",
    "        \n",
    "        # Subplot (1,0): GCC raw\n",
    "        ax2.coastlines()\n",
    "        if not match_cbar: [vmin, vmax] = [vmin_gcc, vmax_gcc]\n",
    "        plot2 = ax2.imshow(ds1, extent=(minlon_raw, maxlon_raw, minlat_raw, maxlat_raw), \n",
    "                           cmap=WhGrYlRd,vmin=vmin, vmax=vmax)\n",
    "        ax2.set_title('GCC Raw{}\\n{}'.format(subtitle_extra,llres_raw)) \n",
    "        cb = plt.colorbar(plot2, ax=ax2, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.1 or (vmax-vmin) > 100:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)\n",
    "        \n",
    "        # Subplot (1,1): GCC regridded\n",
    "        ax3.coastlines()\n",
    "        if not match_cbar: [vmin, vmax] = [vmin_gcc, vmax_gcc]\n",
    "        plot3 = ax3.imshow(gcc_ll, extent=(minlon_cmp, maxlon_cmp, minlat_cmp, maxlat_cmp), \n",
    "                           cmap=WhGrYlRd,vmin=vmin, vmax=vmax)\n",
    "        ax3.set_title('GCC Regridded\\n{}'.format(llres_cmp))\n",
    "        cb = plt.colorbar(plot3, ax=ax3, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.1 or (vmax-vmin) > 100:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)\n",
    "            \n",
    "        # Subplot (2,0): Difference\n",
    "        gc_absdiff = gchp_ll - gcc_ll\n",
    "        diffabsmax = max([np.abs(gc_absdiff.min()), np.abs(gc_absdiff.max())])\n",
    "        [vmin, vmax] = [-diffabsmax, diffabsmax]\n",
    "        ax4.coastlines()\n",
    "        gc_absdiff.plot.imshow\n",
    "        plot4 = ax4.imshow(gc_absdiff, cmap='RdBu_r', extent=(minlon_cmp, maxlon_cmp, minlat_cmp, maxlat_cmp), \n",
    "                           vmin=vmin, vmax=vmax)\n",
    "        ax4.set_title('Difference\\n(GCHP - GCC)')\n",
    "        cb = plt.colorbar(plot4, ax=ax4, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.1 or (vmax-vmin) > 100:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)\n",
    "        \n",
    "        # Subplot (2,1): Fractional Difference (restrict to +/-2)\n",
    "        gc_fracdiff = (gchp_ll - gcc_ll) / gcc_ll\n",
    "        if full_ratio_range: [vmin, vmax] = [None, None]\n",
    "        else: [vmin, vmax] = [-2, 2]\n",
    "        ax5.coastlines()\n",
    "        plot5 = ax5.imshow(gc_fracdiff, vmin=vmin, vmax=vmax, cmap='RdBu_r', \n",
    "                           extent=(minlon_cmp, maxlon_cmp, minlat_cmp, maxlat_cmp))\n",
    "        ax5.set_title('Fractional Difference\\n(GCHP-GCC)/GCC')\n",
    "        cb = plt.colorbar(plot5, ax=ax5, orientation='horizontal', pad=0.10)\n",
    "        cb.set_clim(vmin=vmin, vmax=vmax)\n",
    "        cb.set_label('unitless')      \n",
    "            \n",
    "        if savepdf:    \n",
    "            pdf.savefig(figs)\n",
    "            plt.close(figs)\n",
    "            \n",
    "    if savepdf: pdf.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function to plot and compare zonal means of GCHP and GCC\n",
    "\n",
    "Note: Zonal mean is defined as the mean value across a constant latitudinal band. Calculating zonal mean from cubed sphere and comparing with GCC is not straight-forward. This function therefore regrids both model outputs to a common lat-lon grid resolution prior to plotting. Many of the same features available for plotting a single level above are also available for this function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare_zonal_mean(dgcc, dgchp, varlist=None, weightsdir='.', itime=0, llres_raw='4x5', \n",
    "                       llres_cmp='1x1.25', savepdf=False, pdfname='gchp_vs_gcc_map.pdf', match_cbar=True, \n",
    "                       full_ratio_range=False, normalize_by_area=False, area1=None, area2=None,\n",
    "                      flip_vert=False):\n",
    "\n",
    "    # If no varlist is passed, plot all 3D variables in the dataset\n",
    "    if varlist == None:\n",
    "        [commonvars, commonvars2D, varlist] = compare_varnames(dgcc, dgchp)\n",
    "        print('Plotting all 3D variables')\n",
    "    n_var = len(varlist)\n",
    "    \n",
    "    # Get lat-lon grids and regridder. Assume regridding weights have already been generated\n",
    "    llgrid_raw = make_grid_LL(llres_raw)\n",
    "    llgrid_cmp = make_grid_LL(llres_cmp)\n",
    "    ll_regridder = make_regridder_L2L(llres_raw, llres_cmp, weightsdir=weightsdir, reuse_weights=True)\n",
    "\n",
    "    # Get cubed sphere grid and regridder\n",
    "    csres = dgchp['lon'].size\n",
    "    [csgrid, csgrid_list] = make_grid_CS(csres)\n",
    "    cs_regridder_list = make_regridder_C2L(csres, llres_cmp, weightsdir=weightsdir, reuse_weights=True)\n",
    "    \n",
    "    # Universal plot setup\n",
    "    xtick_positions = np.arange(-90,91,30)\n",
    "    xticklabels = ['{}$\\degree$'.format(x) for x in xtick_positions]\n",
    "    ytick_positions = np.arange(0,61,20)\n",
    "    yticklabels = [str(y) for y in ytick_positions]\n",
    "    \n",
    "    # Create pdf (if saving)\n",
    "    if savepdf:\n",
    "        print('\\nCreating {} for {} variables'.format(pdfname, n_var))\n",
    "        pdf = PdfPages(pdfname)\n",
    "\n",
    "    # Loop over variables\n",
    "    for ivar in range(n_var):\n",
    "        if savepdf: print('{} '.format(ivar), end='')\n",
    "        varname = varlist[ivar]\n",
    "        \n",
    "        # Do some checks: dimensions and units\n",
    "        varndim = dgchp[varname].ndim\n",
    "        varndim2 = dgcc[varname].ndim\n",
    "        if 'ilev' in dgcc[varname].dims: nlev = 73\n",
    "        else: nlev = 72\n",
    "        assert varndim == varndim2, 'GCHP and GCC dimensions do not agree for {}!'.format(varname)\n",
    "        units_raw = dgchp[varname].units\n",
    "        units2 = dgcc[varname].units\n",
    "        assert units_raw == units2, 'GCHP and GCC units do not match for {}!'.format(varname)\n",
    "        \n",
    "        # Set plot extent\n",
    "        extent=(-90,90,0,nlev)\n",
    "        \n",
    "        # if normalizing by area, adjust units to be per m2, and adjust title string\n",
    "        units = units_raw\n",
    "        subtitle_extra = ''\n",
    "         \n",
    "        # Slice the data\n",
    "        ds1 = dgcc[varname].isel(time=itime)\n",
    "        ds2 = dgchp[varname].isel(time=itime)\n",
    "\n",
    "        # if normalizing by area, transform on the native grid and adjust units and subtitle string\n",
    "        exclude_list = ['WetLossConvFrac','Prod_','Loss_']\n",
    "        if normalize_by_area and not any(s in varname for s in exclude_list):\n",
    "            ds1.values = ds1.values / area1.values[np.newaxis,:,:]\n",
    "            ds2.values = ds2.values / area2.values[np.newaxis,:,:]\n",
    "            units = '{} m-2'.format(units_raw)\n",
    "            subtitle_extra = ', Normalized by Area'\n",
    "            \n",
    "        # Regrid the slices\n",
    "        if flip_vert: ds2.data = ds2.data[::-1,:,:]\n",
    "        csdata = ds2.data.reshape(nlev,6,csres,csres).swapaxes(0,1)\n",
    "        gchp_ll = np.zeros([nlev, llgrid_cmp['lat'].size, llgrid_cmp['lon'].size])\n",
    "        for i in range(6):\n",
    "            regridder = cs_regridder_list[i]\n",
    "            gchp_ll += regridder(csdata[i])\n",
    "        gcc_ll = ll_regridder(ds1)\n",
    "        \n",
    "        # Calculate zonal mean of the regridded data\n",
    "        gchp_zm = gchp_ll.mean(axis=2)\n",
    "        gcc_zm = gcc_ll.mean(axis=2)\n",
    "            \n",
    "        # Get min and max for colorbar limits\n",
    "        [vmin_gchp, vmax_gchp] = [gchp_zm.min(), gchp_zm.max()]\n",
    "        [vmin_gcc, vmax_gcc] = [gcc_zm.min(), gcc_zm.max()]\n",
    "        vmin_cmn = np.min([vmin_gchp, vmin_gcc])\n",
    "        vmax_cmn = np.max([vmax_gchp, vmax_gcc])\n",
    "        if match_cbar: [vmin, vmax] = [vmin_cmn, vmax_cmn]\n",
    "        \n",
    "        # Create 2x2 figure\n",
    "        figs, ((ax0, ax1), (ax2, ax3)) = plt.subplots(2, 2, figsize=[12,12], \n",
    "                                                      subplot_kw={'projection': crs.PlateCarree()})\n",
    "        # Give the page a title\n",
    "        offset = 0.96\n",
    "        fontsize=25\n",
    "        figs.suptitle('{}, Zonal Mean'.format(varname), fontsize=fontsize, y=offset)\n",
    "\n",
    "        # Subplot 0: GCHP regridded\n",
    "        if not match_cbar: [vmin, vmax] = [vmin_gchp, vmax_gchp]\n",
    "        plot0 = ax0.imshow(gchp_zm, cmap=WhGrYlRd, extent=extent, vmin=vmin, vmax=vmax)\n",
    "        ax0.set_title('GCHP Regridded{}\\n{}'.format(subtitle_extra, llres_cmp))\n",
    "        ax0.set_aspect('auto')\n",
    "        ax0.set_xticks(xtick_positions)\n",
    "        ax0.set_xticklabels(xticklabels)\n",
    "        ax0.set_yticks(ytick_positions)\n",
    "        ax0.set_yticklabels(yticklabels)\n",
    "        cb = plt.colorbar(plot0, ax=ax0, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.001 or (vmax-vmin) > 1000:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)\n",
    "        \n",
    "        # Subplot 1: GCC regridded\n",
    "        if not match_cbar: [vmin, vmax] = [vmin_gcc, vmax_gcc]\n",
    "        plot1 = ax1.imshow(gcc_zm, cmap=WhGrYlRd, extent=extent, vmin=vmin, vmax=vmax)\n",
    "        ax1.set_title('GCC Regridded{}\\n{}'.format(subtitle_extra, llres_cmp))\n",
    "        ax1.set_aspect('auto')\n",
    "        ax1.set_xticks(xtick_positions)\n",
    "        ax1.set_xticklabels(xticklabels)\n",
    "        ax1.set_yticks(ytick_positions)\n",
    "        ax1.set_yticklabels(yticklabels)\n",
    "        cb = plt.colorbar(plot1, ax=ax1, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.001 or (vmax-vmin) > 1000:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)\n",
    "            \n",
    "        # Subplot 2: Difference\n",
    "        gc_absdiff = gchp_zm - gcc_zm\n",
    "        diffabsmax = max([np.abs(gc_absdiff.min()), np.abs(gc_absdiff.max())])\n",
    "        [vmin, vmax] = [-diffabsmax, diffabsmax]\n",
    "        plot2 = ax2.imshow(gc_absdiff, cmap='RdBu_r', extent=extent, vmin=vmin, vmax=vmax)\n",
    "        ax2.set_title('Difference\\n(GCHP - GCC)')\n",
    "        ax2.set_aspect('auto')\n",
    "        ax2.set_xticks(xtick_positions)\n",
    "        ax2.set_xticklabels(xticklabels)\n",
    "        ax2.set_yticks(ytick_positions)\n",
    "        ax2.set_yticklabels(yticklabels)\n",
    "        cb = plt.colorbar(plot2, ax=ax2, orientation='horizontal', pad=0.10)\n",
    "        if (vmax-vmin) < 0.001 or (vmax-vmin) > 1000:\n",
    "            cb.locator = ticker.MaxNLocator(nbins=4)\n",
    "            cb.update_ticks()\n",
    "        cb.set_label(units)\n",
    "        \n",
    "        # Subplot 3: Fractional Difference (restrict to +/-2)\n",
    "        gc_fracdiff = (gchp_zm - gcc_zm) / gcc_zm\n",
    "        if full_ratio_range: [vmin, vmax] = [None, None]\n",
    "        else: [vmin, vmax] = [-2, 2]\n",
    "        plot3 = ax3.imshow(gc_fracdiff, vmin=vmin, vmax=vmax, cmap='RdBu_r', extent=extent)\n",
    "        ax3.set_title('Fractional Difference\\n(GCHP-GCC)/GCC')\n",
    "        ax3.set_aspect('auto')\n",
    "        ax3.set_xticks(xtick_positions)\n",
    "        ax3.set_xticklabels(xticklabels)\n",
    "        ax3.set_yticks(ytick_positions)\n",
    "        ax3.set_yticklabels(yticklabels)\n",
    "        cb = plt.colorbar(plot3, ax=ax3, orientation='horizontal', pad=0.10)\n",
    "        cb.set_clim(vmin=vmin, vmax=vmax)\n",
    "        cb.set_label('unitless')      \n",
    "            \n",
    "        if savepdf:    \n",
    "            pdf.savefig(figs)\n",
    "            plt.close(figs)\n",
    "            \n",
    "    if savepdf: pdf.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Directories with Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "datadir = '/n/home08/elundgren/GC/testruns/12.1.0_gchp_dev'\n",
    "gccdatadir = os.path.join(datadir,'geosfp_4x5_standard_1day')\n",
    "gchpdatadir = os.path.join(datadir,'gchp_standard/c48_1day/data')\n",
    "check_paths(gccdatadir, gchpdatadir)\n",
    "\n",
    "# Also set directory to store plots and where to store regridding weights\n",
    "plotsdir = os.path.join(datadir,'plots', '1day')\n",
    "weightsdir = '/n/home08/elundgren/GC/python/regrid_weights'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print GCC netcdf filenames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GCC_files = [k for k in os.listdir(gccdatadir) if '.nc' in k]\n",
    "for k in GCC_files:\n",
    "    print(k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print GCHP netcdf filenames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GCHP_files = [k for k in os.listdir(gchpdatadir) if '.nc' in k]\n",
    "for k in GCHP_files:\n",
    "    print(k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Print GCHP netcdf filenames one up from data directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GCHP_files = [k for k in os.listdir(os.path.join(gchpdatadir,'..')) if '.nc' in k]\n",
    "for k in GCHP_files:\n",
    "    print(k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create Regridding Weights\n",
    "\n",
    "The xESMF Regridder object can create weights for regridding from one grid resolution to another, or it can be provided weights that are pre-generated to save time. Generating weights therefore only needs to be done once for each combination of input and output grid type and resolution. This can be done by calling the make_regridder functions and passing reuse_weights=False. Once the weights exist, all subsequent calls to make a regridder can be passed reuse_weights=True along with the location of the pre-existing weights.\n",
    "\n",
    "Weights are created on the fly in the plotting functions above. However, you need to specify the location where you want weights stored. Edit and run the cell below to set that. The cell beyond is for example purposes only and therefore is commented out. Uncomment and run if you want specific weights for use outside of the plotting functions, or to play around with seeing the regridder object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#gcc_res = '4x5'  # GEOS-Chem Classic run resolution\n",
    "#gchp_res = 48    # GCHP run resolution\n",
    "#compare_res = '1x1.25'  # Lat-lon resolution for comparison\n",
    "#regridder_C2L_list = make_regridder_C2L(gchp_res, compare_res, \n",
    "#                                        weightsdir=weightsdir, reuse_weights=False)\n",
    "#regridder_L2L = make_regridder_L2L(gcc_res, compare_res, \n",
    "#                                   weightsdir=weightsdir, reuse_weights=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare Restart Files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#gcc_restart = os.path.join(gccdatadir,'GEOSChem_restart.201607010000.nc')\n",
    "#gchp_res=48\n",
    "#gchp_restart = os.path.join(gchpdatadir,'..','..','initial_GEOSChem_rst.c{}_standard.nc'.format(gchp_res))\n",
    "#check_paths(gcc_restart, gchp_restart)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#ds_gcc = xr.open_dataset(gcc_restart)\n",
    "#ds_gchp = xr.open_dataset(gchp_restart)\n",
    "#[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#get_stats(ds_gcc, ds_gchp, 'SPC_O3')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Comparing restarts requires flipping GCHP in the vertical. Do this by passing flip_vert=True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compare_zonal_mean(ds_gcc, ds_gchp, varlist=['SPC_DHDN'], weightsdir=weightsdir, flip_vert=True )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#compare_single_level(ds_gcc, ds_gchp, varlist=['SPC_DHDN'], ilev=0, weightsdir=weightsdir, flip_vert=True )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Store Grid Areas for Normalizing by Area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'StateMet'\n",
    "gcc_day = '20160701' # YYYYMMDD start date of time-averaging period\n",
    "gcc_time='0000'      # HHmm: GCC filename shows start time of time-averaging period\n",
    "\n",
    "gchp_collection = 'StateMet_avg'\n",
    "gchp_day = '20160701' # YYYYMMDD\n",
    "gchp_time='1200'      # HHmm: GCHP filename shows mid-point time of time-averaging period\n",
    "\n",
    "gcc_met = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_met = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_met, gchp_met)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc_met = xr.open_dataset(gcc_met)\n",
    "ds_gchp_met = xr.open_dataset(gchp_met)\n",
    "varname = 'Met_AREAM2'\n",
    "gcc_area = ds_gcc_met[varname].isel(lev=0,time=0)\n",
    "gchp_area = ds_gchp_met[varname].isel(lev=0,time=0)\n",
    "print('GCC native grid area array shape: {}'.format(gcc_area.shape))\n",
    "print('GCHP native grid area array shape: {}'.format(gchp_area.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Species concentration quick plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'SpeciesConc'\n",
    "gcc_day = '20160701' # YYYYMMDD start date of time-averaging period\n",
    "gcc_time='0000'      # HHmm: GCC filename shows start time of time-averaging period\n",
    "\n",
    "gchp_collection = 'SpeciesConc_avg'\n",
    "gchp_day = '20160701' # YYYYMMDD\n",
    "gchp_time='1200'      # HHmm: GCHP filename shows mid-point time of time-averaging period\n",
    "\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varname='SpeciesConc_BrSALC'\n",
    "#get_stats(ds_gcc, ds_gchp, varname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compare_single_level(ds_gcc, ds_gchp, varlist=[varname], weightsdir=weightsdir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#compare_zonal_mean(ds_gcc, ds_gchp, varlist=[varname], weightsdir=weightsdir )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## All diagnostics collections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GCHP_files = [k for k in os.listdir(gchpdatadir) if '.nc' in k]\n",
    "for k in GCHP_files:\n",
    "    print(k)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Instantaneous"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_day = '20160702' # YYYYMMDD start date of time-averaging period\n",
    "gcc_time='0000' \n",
    "gchp_day = '20160702' # YYYYMMDD\n",
    "gchp_time='0000' "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'SpeciesConc_inst'\n",
    "gchp_collection = 'SpeciesConc_inst'\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=commonvars\n",
    "pdfname = os.path.join(plotsdir,'{}_Surface.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir,\n",
    "                     savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Time-averaged"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_day = '20160701' # YYYYMMDD start date of time-averaging period\n",
    "gcc_time='0000' \n",
    "gchp_day = '20160701' # YYYYMMDD\n",
    "gchp_time='1200' "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### SpeciesConc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'SpeciesConc'\n",
    "gchp_collection = 'SpeciesConc_avg'\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=commonvars\n",
    "pdfname = os.path.join(plotsdir,'{}_Surface.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir,\n",
    "                     savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Budget\n",
    "Units are kg/s so plotting requires normalization by area."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'Budget'\n",
    "gchp_collection = gcc_collection\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=commonvars\n",
    "pdfname = os.path.join(plotsdir,'{}.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname, normalize_by_area=True, area1=gcc_area,\n",
    "                    area2=gchp_area)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### CloudConvFlux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'CloudConvFlux'\n",
    "gchp_collection = gcc_collection\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samplevars = get_n_random_vars(commonvars,30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=samplevars\n",
    "pdfname = os.path.join(plotsdir,'{}_Surface_RandomSample.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean_RandomSample.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir,\n",
    "                     savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ConcAfterChem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'ConcAfterChem'\n",
    "gchp_collection = gcc_collection\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=commonvars\n",
    "pdfname = os.path.join(plotsdir,'{}_Surface.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir,\n",
    "                     savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### JValues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'JValues'\n",
    "gchp_collection = gcc_collection\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samplevars = get_n_random_vars(commonvars,30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=samplevars\n",
    "pdfname = os.path.join(plotsdir,'{}_Surface_RandomSample.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean_RandomSample.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir,\n",
    "                     savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Aerosols"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'Aerosols'\n",
    "gchp_collection = gcc_collection\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=commonvars\n",
    "pdfname = os.path.join(plotsdir,'{}_Surface.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir,\n",
    "                     savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### AerosolMass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'AerosolMass'\n",
    "gchp_collection = gcc_collection\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=commonvars\n",
    "pdfname = os.path.join(plotsdir,'{}_Surface.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir,\n",
    "                     savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Emissions\n",
    "GCHP emissions diagnostics are vertically inverted relative to GC classic and therefore must be flipped for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'Emissions'\n",
    "gchp_collection = gcc_collection\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=commonvars2D\n",
    "pdfname = os.path.join(plotsdir,'{}_2D.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname, flip_vert=True )\n",
    "varlist=commonvars3D\n",
    "pdfname = os.path.join(plotsdir,'{}_3D_Surface.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname, flip_vert=True )\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir,\n",
    "                     savepdf=True, pdfname=pdfname, flip_vert=True )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### DryDep"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'DryDep'\n",
    "gchp_collection = gcc_collection\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samplevars = get_n_random_vars(commonvars,30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=samplevars\n",
    "pdfname = os.path.join(plotsdir,'{}_RandomSample.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### StateMet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'StateMet'\n",
    "gchp_collection = 'StateMet_avg'\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=commonvars2D\n",
    "pdfname = os.path.join(plotsdir,'{}_2D.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )\n",
    "varlist=commonvars3D\n",
    "pdfname = os.path.join(plotsdir,'{}_3D_Surface.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### StateChm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'StateChm'\n",
    "gchp_collection = gcc_collection\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=commonvars\n",
    "pdfname = os.path.join(plotsdir,'{}_Surface.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir,\n",
    "                     savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ProdLoss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'ProdLoss'\n",
    "gchp_collection = gcc_collection\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=commonvars\n",
    "pdfname = os.path.join(plotsdir,'{}_Surface.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname, normalize_by_area=True, area1=gcc_area,\n",
    "                     area2=gchp_area)\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### WetLossLS\n",
    "Units are kg/s and therefore data must be normalized by area for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'WetLossLS'\n",
    "gchp_collection = gcc_collection\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samplevars = get_n_random_vars(commonvars,30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=samplevars\n",
    "pdfname = os.path.join(plotsdir,'{}_Surface_RandomSample.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname, normalize_by_area=True, area1=gcc_area,\n",
    "                    area2=gchp_area)\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean_RandomSample.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname, normalize_by_area=True, area1=gcc_area,\n",
    "                    area2=gchp_area)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### WetLossConv\n",
    "Units are kg/s and data must therefore be normalized by area for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'WetLossConv'\n",
    "gchp_collection = gcc_collection\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samplevars = get_n_random_vars(commonvars,30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=samplevars\n",
    "pdfname = os.path.join(plotsdir,'{}_Surface_RandomSample.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname, normalize_by_area=True, area1=gcc_area,\n",
    "                    area2=gchp_area)\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean_RandomSample.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname, normalize_by_area=True, area1=gcc_area,\n",
    "                    area2=gchp_area)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### LevelEdgeDiags"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gcc_collection = 'LevelEdgeDiags'\n",
    "gchp_collection = gcc_collection\n",
    "gcc_spc = get_gcc_filepath(gccdatadir, gcc_collection, gcc_day, gcc_time)\n",
    "gchp_spc = get_gchp_filepath(gchpdatadir, gchp_collection, gchp_day, gchp_time)\n",
    "check_paths(gcc_spc, gchp_spc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_gcc = xr.open_dataset(gcc_spc)\n",
    "ds_gchp = xr.open_dataset(gchp_spc)\n",
    "[commonvars, commonvars2D, commonvars3D] = compare_varnames(ds_gcc, ds_gchp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "varlist=commonvars\n",
    "pdfname = os.path.join(plotsdir,'{}_Surface.pdf'.format(gchp_collection))\n",
    "compare_single_level(ds_gcc, ds_gchp, varlist=varlist, ilev=0, weightsdir=weightsdir, \n",
    "                     savepdf=True, pdfname=pdfname )\n",
    "pdfname = os.path.join(plotsdir,'{}_ZonalMean.pdf'.format(gchp_collection))\n",
    "compare_zonal_mean(ds_gcc, ds_gchp, varlist=varlist, weightsdir=weightsdir,\n",
    "                     savepdf=True, pdfname=pdfname )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}