.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/emissions/plot_emg.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_emissions_plot_emg.py: Calculate Emissions by Fitting a Modified Guassian ================================================== This example uses Dallas Fort Worth Texas as an illustrative application of this technique. For illustration, only 3 days are used. In practice, it should use as whole year or the ozone season to get a more robust estimate. Steps: 1. Create a 3km custom L3 product on the HRRR grid, 2. Retrieve HRRR winds at center of domain, 3. Rotate daily rasters so that wind blows from left to right, 4. Fit Eq 1 of Goldberg et al. (10.5194/acp-22-10875-2022; ~exponnorm), 5. Derive tau parameter and estimate emissions. The 4 day result 41 moles/s, which is in good agreement with the 2017 NEI[1,2] 60870 short-tons/yr (38 moles/s) and larger than the 2020 NEI[3,4] 47110 short- tons/yr (29 moles/s). These days were selected because of the reasonable agreement. In actual application, you should check the rotation of each day and apply to many days. The alignment of winds with plumes is imperfect due to meandering winds. [1] https://enviro.epa.gov/envirofacts/embed/nei?pType=TIER&pReport=county&pState=&pState=48&pPollutant=&pPollutant=NOX&pSector=&pYear=2017&pCounty=48439&pTier=&pWho=NEI [2] https://enviro.epa.gov/envirofacts/embed/nei?pType=TIER&pReport=county&pState=&pState=48&pPollutant=&pPollutant=NOX&pSector=&pYear=2017&pCounty=48113&pTier=&pWho=NEI [3] https://enviro.epa.gov/envirofacts/embed/nei?pType=TIER&pReport=county&pState=&pState=48&pPollutant=&pPollutant=NOX&pSector=&pYear=2020&pCounty=48439&pTier=&pWho=NEI [4] https://enviro.epa.gov/envirofacts/embed/nei?pType=TIER&pReport=county&pState=&pState=48&pPollutant=&pPollutant=NOX&pSector=&pYear=2020&pCounty=48113&pTier=&pWho=NEI .. GENERATED FROM PYTHON SOURCE LINES 29-35 .. code-block:: default import matplotlib.pyplot as plt import numpy as np import pyrsig import pandas as pd .. GENERATED FROM PYTHON SOURCE LINES 36-42 Define Location and Date Range ------------------------------ - Location is a single point - Date range can be a few days at a time - vbbox is the box to retrieve VCDs .. GENERATED FROM PYTHON SOURCE LINES 42-49 .. code-block:: default workdir = 'dfw' loc = (-96.95, 32.85) dates = pd.date_range('2023-05-01T00Z', '2023-05-04T00Z', freq='1d') vbbox = np.array(loc + loc) + np.array([-1, -1, 1, 1]) * 1.5 .. GENERATED FROM PYTHON SOURCE LINES 50-53 Get TropOMI NO2 Over vbbox -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 53-70 .. code-block:: default vrsig = pyrsig.RsigApi( bbox=vbbox, workdir=workdir, grid_kw='HRRR3K', gridfit=True ) vkey = 'tropomi.offl.no2.nitrogendioxide_tropospheric_column' ds = vrsig.to_ioapi(bdate=dates, key=vkey) # Make the output square by windowing in on the center n = min(ds.sizes['COL'], ds.sizes['ROW']) n = n - 1 if n % 2 == 0 else n dn = n / 2 mc = ds.COL.mean().astype('i') + 0.5 # Get nearest cell center mr = ds.ROW.mean().astype('i') + 0.5 # Get nearest cell center ds = ds.sel(ROW=slice(mr - dn, mr + dn), COL=slice(mc - dn, mc + dn)) # removing missing values ds = ds.where(lambda x: x > -9e30) .. GENERATED FROM PYTHON SOURCE LINES 71-76 Get Winds At Domain Center -------------------------- - Find the domain center - Create a 0.01 degree box around it - retrieve the HRRR 80m wind components (requires xdr) .. GENERATED FROM PYTHON SOURCE LINES 76-90 .. code-block:: default # Find the center lon/lat to query for wind speed llds = ds.isel(TSTEP=0, LAY=0).sel(ROW=ds.ROW.mean(), COL=ds.COL.mean(), method='nearest') wloc = float(llds.LONGITUDE), float(llds.LATITUDE) wbbox = np.array(wloc + wloc) + np.array([-1, -1, 1, 1]) * .005 wrsig = pyrsig.RsigApi(bbox=wbbox, workdir=workdir) wkey = 'hrrr.wind_80m' df = wrsig.to_dataframe(bdate=dates, key=wkey, unit_keys=False, parse_dates=True, backend='xdr') df = df.loc[df['time'].dt.tz_convert(None).isin(ds.TSTEP.data)] # only if a valid record is avail # Add wind components to the VCD dataset ds['wind_80m_u'] = ('TSTEP',), df.wind_80m_u.values, dict(units='m/s') ds['wind_80m_v'] = ('TSTEP',), df.wind_80m_v.values, dict(units='m/s') .. GENERATED FROM PYTHON SOURCE LINES 91-95 Fit Exponential Modified Gaussian --------------------------------- - Convert from molecules/cm2 to mole/m2 - Fit using u/v wind components and the width of cells in m .. GENERATED FROM PYTHON SOURCE LINES 95-104 .. code-block:: default molperm2 = ds.NO2 / 6.022e23 * 1e4 dxm = ds.XCELL emgout = pyrsig.emiss.fitemg(molperm2, ds.wind_80m_u, ds.wind_80m_v, dx=dxm) emgout['wind_80m_u'] = ('time',), ds.wind_80m_u.data emgout['wind_80m_v'] = ('time',), ds.wind_80m_v.data emgout.to_netcdf(f'{workdir}/emg.nc') .. GENERATED FROM PYTHON SOURCE LINES 105-109 Plot EMG Result and Derived Props --------------------------------- - Create a 3-panel plot (unrotated VCD mean, rotated VCD mean and VLD) - Fit using u/v wind components and the width of cells in m .. GENERATED FROM PYTHON SOURCE LINES 109-141 .. code-block:: default gskw = dict(left=0.05, right=0.95, bottom=0.3) fig, axx = plt.subplots(1, 2, figsize=(12, 6), gridspec_kw=gskw) axw = .9 / 2.1 ds.NO2.mean('TSTEP').plot(ax=axx[0]) emgout.ROTATED_VCD.mean('time').plot(ax=axx[1]) bbox = axx[1].get_position() bbox.y0 = 0.075 bbox.y1 = 0.225 axx[1].xaxis.set_major_formatter(plt.matplotlib.ticker.NullFormatter()) inax = fig.add_axes(bbox) emgout.vldhat.plot(ax=inax) emgout.vld.plot(ax=inax) inax.set_xticks(axx[1].get_xticks()) inax.set_xlim(axx[1].get_xlim()) tag = ', '.join([f'{k}={v:.1f}' for k, v in emgout.vldhat.attrs.items()]) tag = f'E={emgout.emis_no2:.0f} [{emgout.emis_no2.units}] = $\\alpha \\overline{{ws}} / x0$\n{tag}' fig.text(0.05, 0.075, tag, ha='left', va='bottom') figpath = f'{workdir}/emg.png' fig.savefig(figpath) # Requires geopandas to add county and primary roads try: import geopandas as gpd cntypath = 'https://www2.census.gov/geo/tiger/TIGER2023/COUNTY/tl_2023_us_county.zip' cnty = gpd.read_file(cntypath, bbox=tuple(vbbox)).to_crs(ds.crs_proj4) roadpath = 'https://www2.census.gov/geo/tiger/TIGER2023/PRIMARYROADS/tl_2023_us_primaryroads.zip' road = gpd.read_file(roadpath, bbox=tuple(vbbox)).to_crs(ds.crs_proj4) cnty.plot(facecolor='none', edgecolor='gray', linewidth=.6, ax=axx[0]) road.plot(facecolor='none', edgecolor='gray', linewidth=.3, ax=axx[0]) fig.savefig(figpath) except Exception as e: print('failed to add map', str(e)) .. image-sg:: /auto_examples/emissions/images/sphx_glr_plot_emg_001.png :alt: LAY = 0.9975 :srcset: /auto_examples/emissions/images/sphx_glr_plot_emg_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /opt/hostedtoolcache/Python/3.13.7/x64/lib/python3.13/site-packages/pyogrio/raw.py:198: RuntimeWarning: HTTP response code on https://www2.census.gov/geo/tiger/TIGER2023/COUNTY/tl_2023_us_county.zip: 0 return ogr_read( failed to add map '/vsizip/vsicurl/https://www2.census.gov/geo/tiger/TIGER2023/COUNTY/tl_2023_us_county.zip' does not exist in the file system, and is not recognized as a supported dataset name. .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 31.979 seconds) .. _sphx_glr_download_auto_examples_emissions_plot_emg.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_emg.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_emg.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_