.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/satellite/run_satellite_01_simpleomi.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_satellite_run_satellite_01_simpleomi.py: Simple OMI ========== This example compares CAMx to the Ozone Monitoring Instrument (OMI) for NO2 tropospheric vertical column densitities (tropNO2). OMI's tropNO2 is derived from a measurement along the path light from the sun and reflected off the surface of the earth. The original measurement is called a slant column density (SCD). The SCD is split into tropospheric and stratospheric parts, based on a model profile (a priori). The a priori and slant path are combined with a sensitivity vector (scattering weight) to convert the tropospheric SCD to tropNO2 [molecules/cm2]. CAMx can output mixing ratios at all levels, which can be converted to molecules/cm2 and summed to make a tropNO2 product. Note, however, that the CAMx tropNO2 inherently assumes a different model profile. A more complete comparison would adjust OMI to use CAMx's profile, or adjust CAMx to use OMI's profile. Each has strengths and weaknesses. *Reminder*: You must have already activated your python environment. .. GENERATED FROM PYTHON SOURCE LINES 22-24 Imports and Folders ''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 24-31 .. code-block:: Python import pyrsig import os os.makedirs('outputs', exist_ok=True) .. GENERATED FROM PYTHON SOURCE LINES 32-34 Configuration ''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 34-49 .. code-block:: Python cdate = '20160610' # Using first available CAMx date odate = '2017-06-10' # oddly Jun 10-11 of 2016 have no OMI over the US # Use dates to find avrg file and 3d met file cpath = f'../../camx/outputs/CAMx.v7.32.36.12.{cdate}.3D.avrg.grd01.nc' # If you used camx.tar.gz for inputs, update the met folder to metss mpath = f'../../camx/met/camx.3d.36km.{cdate}.nc' # Define output paths. outstem = f'outputs/OMINO2_CAMx.v7.32.36.12.{cdate}.3D.avrg.grd01' outpath = outstem + '.nc' tilepath = outstem + '_tile.png' scatterpath = outstem + '_scatter.png' .. GENERATED FROM PYTHON SOURCE LINES 50-52 Data Acquisition '''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 52-63 .. code-block:: Python # Open a CAMx 3D CONC file for NO2 cf = pyrsig.open_ioapi(cpath) # Open a 3d met file (pressure, z, temperature, humiditiy) and match avrg file time mkeys = ['pressure', 'z', 'temperature', 'humidity'] mf = pyrsig.open_ioapi(mpath)[mkeys].interp(TSTEP=cf.TSTEP) # Process OMI for the CAMx grid api = pyrsig.RsigApi(bbox=(-130, 20, -60, 60), workdir='outputs', grid_kw=cf.attrs) omids = api.to_ioapi('omi.l2.omno2.ColumnAmountNO2Trop', bdate=odate) .. GENERATED FROM PYTHON SOURCE LINES 64-70 Derive Conversion Factors ''''''''''''''''''''''''' - OMI tropospheric NO2 column (tropNO2) in molecules/cm2 - CAMx average mixing ratios (X) in ppmV of dry-air by layer (k) - Need dry-air intensity (I) molecules/cm2 by layer - CAMx tropNO2 = $\sum_k{X(k) * I(k)}$ .. GENERATED FROM PYTHON SOURCE LINES 70-88 .. code-block:: Python # Calculate layer depths from interface heights dz = mf['z'].copy() dz[:, 1:] = mf['z'].diff('LAY').data # Combine pressure, temperature, humidity, and dz to get dry-air intensity molecules_per_cm2 = ( mf['pressure'] * 100 # mb -> Pa / 8.314 / mf['temperature'] # Pa -> moles/m3 * (1 - mf['humidity'] / 1e6) # moles/m3 -> moles_dry/m3 * dz * 6.022e23 / 1e4 # moles_dry/m3 -> molecules_dry/m2 ) vcdno2 = (cf['NO2'] * 1e-6 * molecules_per_cm2.data).sum('LAY', keepdims=True) vcdno2.attrs.update( long_name='CAMx_tropNO2', units='molecules/cm2', var_desc='Simple sum of molecules' ) .. GENERATED FROM PYTHON SOURCE LINES 89-91 Organizing for output ''''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 91-103 .. code-block:: Python # Allowing OMI coordinates (time) to override CAMx. omids['CAMx_tropNO2'] = vcdno2.dims, vcdno2.data, vcdno2.attrs # Mask missing values from OMI. outf = omids.where(omids['COLUMNAMOUNTNO2'] > -9e36) # Make OMI naming more clear. outf = outf.rename(COLUMNAMOUNTNO2='OMI_tropNO2') outf['OMI_tropNO2'].attrs.update(long_name='OMI_tropNO2') # Save to disk pyrsig.cmaq.save_ioapi(outf.drop_vars('TFLAG'), outpath) .. GENERATED FROM PYTHON SOURCE LINES 104-106 Visualize Results ''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 106-117 .. code-block:: Python import matplotlib.pyplot as plt import matplotlib.colors as mc import pyrsig import pycno import numpy as np from scipy.stats import linregress # Average potential overlapping overpasses omids = pyrsig.open_ioapi(outpath).mean(('TSTEP', 'LAY')) .. GENERATED FROM PYTHON SOURCE LINES 118-120 Make a tile-plot '''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 120-127 .. code-block:: Python Z = omids[['OMI_tropNO2', 'CAMx_tropNO2']].to_dataarray() Z.attrs.update(long_name='tropNO2', units='molecules/cm2') fca = Z.plot(col='variable') pycno.cno(omids.crs_proj4).drawstates(ax=fca.axs) fca.fig.savefig(tilepath) .. image-sg:: /auto_examples/satellite/images/sphx_glr_run_satellite_01_simpleomi_001.png :alt: variable = OMI_tropNO2, variable = CAMx_tropNO2 :srcset: /auto_examples/satellite/images/sphx_glr_run_satellite_01_simpleomi_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 128-130 Make a scatter-plot ''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 130-148 .. code-block:: Python # Get values from OMI and CAMx separately omivals = omids['OMI_tropNO2'].to_numpy().ravel() camxvals = omids['CAMx_tropNO2'].to_numpy().ravel() bothvalid = ~(np.isnan(omivals) | np.isnan(camxvals)) # Perform least-square regression lr = linregress(omivals[bothvalid], camxvals[bothvalid]) lrlabel = f'y = {lr.slope:.3}x + {lr.intercept:.3} (r={lr.rvalue:.2f})' # Create a hexbin plot fig, ax = plt.subplots(figsize=(6, 6)) pc = ax.hexbin(x=omivals[bothvalid], y=camxvals[bothvalid], gridsize=50, norm=mc.LogNorm(20)) ax.set(xlabel='OMI_tropNO2', ylabel='CAMx_tropNO2', title='tropNO2 molecules/cm2') ax.axline((0, 0), slope=1, label='1:1', color='grey', linestyle='--') ax.axline((0, lr.intercept), slope=lr.slope, label=lrlabel, color='k') ax.legend() fig.savefig(scatterpath) .. image-sg:: /auto_examples/satellite/images/sphx_glr_run_satellite_01_simpleomi_002.png :alt: tropNO2 molecules/cm2 :srcset: /auto_examples/satellite/images/sphx_glr_run_satellite_01_simpleomi_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 12.816 seconds) .. _sphx_glr_download_auto_examples_satellite_run_satellite_01_simpleomi.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: run_satellite_01_simpleomi.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: run_satellite_01_simpleomi.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: run_satellite_01_simpleomi.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_