.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/comparisons/plot_cmaqtropomi.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_comparisons_plot_cmaqtropomi.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_comparisons_plot_cmaqtropomi.py:


Pairing CMAQ with TropOMI Columns
=================================

This example performs a simple comparison of CMAQ to TropOMI. TropOMI retrieves
total column HCHO at overpass time. The CMAQ PHOTDIAG1 outputs the column
integral at each hour.

- This analysis ignores satellite differential sensitivity, since CMAQ PHOTDIAG1
  outputs the whole column with perfect sensitivity.
  
- CMAQ column outputs include all molecules up to the model top defined by VGTOP
  When VGTOP is 5000Pa, this includes some, but not all stratosphere. Further,
  CMAQ does not have many stratosphere reactions/species. Comparisons could be
  improved (e.g., with cmaqsatproc)

.. GENERATED FROM PYTHON SOURCE LINES 19-22

Import libraries
----------------
And tell xarray to keep variable attributes

.. GENERATED FROM PYTHON SOURCE LINES 22-30

.. code-block:: default

    import matplotlib.pyplot as plt
    import pyrsig
    import pycno
    import pandas as pd
    import xarray as xr

    xr.set_options(keep_attrs=True)





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    <xarray.core.options.set_options object at 0x7f3185a143b0>



.. GENERATED FROM PYTHON SOURCE LINES 31-34

Define the focus
----------------
- area/time/domain of interest

.. GENERATED FROM PYTHON SOURCE LINES 34-42

.. code-block:: default


    bbox = (-130, 20, -55, 60)
    start = '2019-06-01T00Z'
    end = '2019-06-07T00Z'
    gdnam = '12US1'
    spc = 'NO2'  # or NO2
    # qtmpl = 'CCTM_PHOTDIAG1_%Y%m%d.nc'








.. GENERATED FROM PYTHON SOURCE LINES 43-46

Set species options
-------------------
- area/time/domain of interest

.. GENERATED FROM PYTHON SOURCE LINES 46-58

.. code-block:: default


    tdkey = {
        'HCHO': 'tropomi.offl.hcho.formaldehyde_tropospheric_vertical_column',
        'NO2': 'tropomi.offl.no2.nitrogendioxide_tropospheric_column'
    }[spc]
    cdkey = {
        'HCHO': 'cmaq.equates.conus.integrated.HCHO_COLUMN',
        'NO2': 'cmaq.equates.conus.integrated.NO2_COLUMN',
    }[spc]
    tvkey = spc
    qvkey = f'{spc}_COLUMN'








.. GENERATED FROM PYTHON SOURCE LINES 59-62

Collect Daily Datasets
----------------------


.. GENERATED FROM PYTHON SOURCE LINES 62-81

.. code-block:: default

    api = pyrsig.RsigApi(bbox=bbox, grid_kw=gdnam, workdir=gdnam)
    dates = pd.date_range(start, end)
    dt = pd.to_timedelta('1d') - pd.to_timedelta('1s')
    todss = []
    cqdss = []
    for bdate in dates:
      edate = bdate + dt
      # cqdss.append(pyrsig.open_ioapi(bdate.strftime(qtmpl))[[qvkey]])
      cqdss.append(
          api.to_ioapi(cdkey, bdate=bdate, edate=edate).where(
              lambda x: x > -9.999e+36
          )
      )
      todss.append(
          api.to_ioapi(tdkey, bdate=bdate, edate=edate).where(
              lambda x: x > -9.999e+36
          )
      )








.. GENERATED FROM PYTHON SOURCE LINES 82-86

Mask CMAQ like TropOMI
----------------------
- implicitly applies overpass time
- implicitly cloud filter (assuming similar)

.. GENERATED FROM PYTHON SOURCE LINES 86-99

.. code-block:: default

    tods = xr.concat(todss, dim='TSTEP')
    cqds = xr.concat(cqdss, dim='TSTEP')
    # sample CMAQ when TropOMI is valid
    tcqds = cqds.where(~tods[tvkey].isnull().isel(LAY=0))
    mcqds = tcqds.mean('TSTEP', keepdims=True)
    mtods = tods.mean('TSTEP', keepdims=True)
    # CMAQ stores NO2_COLUMN and HCHO column as petamolecules/cm2
    mtods[f'{tvkey}_CMAQ'] = mcqds[qvkey].dims, mcqds[qvkey].data * 1e15 
    mtods[f'{tvkey}_CMAQ'].attrs.update(
        units='molecules/cm**2',
        long_name=qvkey.ljust(16),
        var_desc=f'CMAQ total column {qvkey}'.ljust(80)
    )







.. GENERATED FROM PYTHON SOURCE LINES 100-103

Output TropOMI with CMAQ
------------------------


.. GENERATED FROM PYTHON SOURCE LINES 103-108

.. code-block:: default

    mtods.to_netcdf(f'TropOMI_and_CMAQ_AVG_{tvkey}_{start}_{end}_{gdnam}.nc')
    fig, axx = plt.subplots(1, 2, figsize=(12, 4))
    qm = mtods[tvkey].plot(ax=axx[0], vmin=0, vmax=3e15, extend='max')
    mtods[f'{tvkey}_CMAQ'].plot(ax=axx[1], norm=qm.norm, extend='max')
    pycno.cno(mtods.crs_proj4).drawstates(ax=axx)
    fig.savefig(f'TropOMI_and_CMAQ_AVG_{tvkey}_{start}_{end}_{gdnam}.png')


.. image-sg:: /auto_examples/comparisons/images/sphx_glr_plot_cmaqtropomi_001.png
   :alt: LAY = 0.9975, LAY = 0.9975
   :srcset: /auto_examples/comparisons/images/sphx_glr_plot_cmaqtropomi_001.png
   :class: sphx-glr-single-img






.. rst-class:: sphx-glr-timing

   **Total running time of the script:** ( 1 minutes  49.765 seconds)


.. _sphx_glr_download_auto_examples_comparisons_plot_cmaqtropomi.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example




    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_cmaqtropomi.py <plot_cmaqtropomi.py>`

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_cmaqtropomi.ipynb <plot_cmaqtropomi.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_