Skip to content

Instantly share code, notes, and snippets.

@barronh

barronh/emgsatemis.py

Created March 12, 2025 13:20
Show Gist options
  • Save barronh/9bffd91e041d36ac6ea257ce6ff70266 to your computer and use it in GitHub Desktop.
Save barronh/9bffd91e041d36ac6ea257ce6ff70266 to your computer and use it in GitHub Desktop.
Download ZIP
Satellite Derived Emissions from Exponentially Modified Gaussian
Raw
emgsatemis.py
__version__ = '0.2.0'
__doc__ = """EMG Satellite-derived Emissions
===============================
EMG Satellite-derived Emissions (emgsatemis) is designed to estimate emissions
from satellite data following many previous studies based on Bierle et al.
Functions:
- to_polar(x, y, degrees=False)
- to_cartesian(r, alpha, theta=0, degrees=False)
- rotate_array(a, theta, degrees=False)
- EMG(x, alpha, x0, sigma, mu, beta, return_parts=False)
- fitlinedens(y, dx, verbose=0, addyhat=True)
- rotate_and_fitlinedens(orgs, us, vs, degrees=False, dx=1, ...)
- get_gfs(date)
- get_gfs_uv(date, lon, lat, retfile=False)
- squareit(f)
- tropomi_fit(lockey, lon, lat, dates, ...)
- plot_fit(tropf, axx=None, lax=None, nacross=21)
Example Application
-------------------
```python
# %pip install -qq xarray pandas pyrsig netcdf4
# %pip install -qq geopandas # for optional shapefile support
import pandas as pd
import emgsatemis as emg
import os
# allows dates to be distinct
dates = pd.date_range('2019-04-01', '2019-09-30')
lockey = 'Houston'
locattrs = {'lat': 29.720621, 'lon': -95.1}
nacross = 21
alltropf = emg.get_tropomi(lockey, dates=dates, **locattrs, met=None)
fmissing = alltropf['DAILY_NO2'].isnull().mean(('ROW', 'COL', 'LAY'))
tropf = alltropf.sel(TSTEP=fmissing < 0.2)
tropf = emg.fit_tropomi(tropf, nacross=nacross)
fig = emg.plot_fit(tropf)
# wget -r https://www2.census.gov/geo/tiger/TIGER2023/CSA/tl_2023_us_csa.zip
shppath = 'www2.census.gov/geo/tiger/TIGER2023/CSA/tl_2023_us_csa.zip'
if os.path.exists(shppath):
import geopandas as gpd
bbox = tuple(
locattrs[k] + dl
for k, dl in [('lon', -5), ('lat', -5), ('lon', 5), ('lat', 5)]
)
csaf = gpd.read_file(shppath, bbox=bbox).to_crs(tropf.crs_proj4)
csaf.plot(facecolor='none', edgecolor='k', ax=fig.axes[0])
os.makedirs('figs', exist_ok=True)
fig.savefig(f'figs/{lockey}_EMG.png')
```
```bash
python emgsatemis.py Houston -95.1 29.720621 2019-04-01 2019-09-30
```
"""
def to_polar(x, y, degrees=False):
"""
Convert cartesian (x, y)) to polar (r, alpha)
https://en.wikipedia.org/wiki/List_of_common_coordinate_transformations
#From_Cartesian_coordinates
Arguments
---------
x, y : float or np.ndarray
x and y coordinates
degrees : bool
If True, then alpha returned is provided in degrees.
If False, then alpha returned is provided in radians.
Returns
-------
r, alpha : float or np.ndarray
radius and angle in polar coordinates
Example
-------
import numpy as np
from emgsatemis import to_polar
us = np.array([1, 0, -1, 0, 1, -1])
vs = np.array([0, 1, 0, -1, 1, -1])
ws, wd = to_polar(us, vs, degrees=True)
print(ws)
#[ 1. 1. 1. 1. 1.41421356 1.41421356]
print(wd)
# [ 0. 90. 180. -90. 45. -135.]
"""
import numpy as np
x = np.asarray(x)
y = np.asarray(y)
alpha = np.arctan2(y.astype('f'), x.astype('f'))
r = (x**2 + y**2)**.5
if degrees:
alpha = np.degrees(alpha)
return r, alpha
def to_cartesian(r, alpha, theta=0, degrees=False):
"""
Calculate the (x, y) with an optional rotation defined by theta.
https://en.wikipedia.org/wiki/Rotation_of_axes
Arguments
---------
r : float or np.ndarray
radius of polar coordinate
alpha : float or np.ndarray
angle (in radians) of polar coordinate
theta : float
angle to rotate for the cartesian plane
degrees : bool
If True, then alpha and theta are provided in degrees.
If False, then alpha and theta are provided in radians.
Returns
-------
x, y : np.ndarray
cartesian coordinates with optional rotation
Example
-------
import numpy as np
from emgsatemis import to_cartesian
x, y = to_cartesian(np.sqrt(2), np.radians(45))
x, y
# (1.0000000000000002, 1.0)
x, y = to_cartesian(1, 45, -45, degrees=True)
x, y
# (2.220446049250313e-16, 1.0)
"""
import numpy as np
alpha = np.asarray(alpha)
theta = np.asarray(theta)
if degrees:
alpha = np.radians(alpha)
theta = np.radians(theta)
x = r * np.cos(alpha)
y = r * np.sin(alpha)
xp = x * np.cos(theta) + y * np.sin(theta)
yp = -x * np.sin(theta) + y * np.cos(theta)
return xp, yp
def rotate_array(a, theta, degrees=False):
"""
Arguments
---------
a : np.ndarray
symmetric array (n, n)
theta : float
Direction (in degrees) to rotate symmetric array
degrees : bool
If True, then theta is in degrees
Returns
-------
b : np.ndarray
Rotated symmetric grid (n, n)
Example
-------
import numpy as np
from emgsatemis import rotate_array
tmpl = np.eye(5)
tmpl[:2] = 0
rots = np.zeros(tmpl.shape)
orgs = np.zeros(tmpl.shape)
for v in [1, -1]:
for u in [1, -1]:
a = tmpl[::v, ::u]
d = np.degrees(np.arctan2(v, u)) # v, u or y, x
b = rotate_array(a, -d, degrees=d)
orgs += a
rots += b
print(f'# u={u:+.0f};v={v:+.0f};d={d}')
print('# Sum Originals:')
print('#' + np.array2string(np.flipud(orgs)).replace('\n', '\n#'))
print('# Sum Rotated:')
print('#' + np.array2string(np.flipud(rots)).replace('\n', '\n#'))
# u=+1;v=+1;d=45.0
# u=-1;v=+1;d=135.0
# u=+1;v=-1;d=-45.0
# u=-1;v=-1;d=-135.0
# Sum Originals:
#[[1. 0. 0. 0. 1.]
# [0. 1. 0. 1. 0.]
# [0. 0. 4. 0. 0.]
# [0. 1. 0. 1. 0.]
# [1. 0. 0. 0. 1.]]
# Sum Rotated:
#[[0. 0. 0. 0. 0.]
# [0. 0. 0. 0. 0.]
# [0. 0. 4. 4. 4.]
# [0. 0. 0. 0. 0.]
# [0. 0. 0. 0. 0.]]
"""
import numpy as np
ny, nx = a.shape
assert (ny % 2 == 1)
assert (nx % 2 == 1)
n = a.shape[0] // 2
if degrees:
theta = np.radians(theta)
j, i = np.indices(a.shape)
y = j - n
x = i - n
r, alpha = to_polar(x, y)
# Calculate the (x', y') in rotation
xp, yp = to_cartesian(r, alpha, theta)
# xp[n, n] = 0
# yp[n, n] = 0
# Sample at the closest grid cell
ip = np.minimum(n, np.maximum(
-n, np.ma.masked_invalid(np.round(xp))
)).astype('i').filled(0) + n
jp = np.minimum(n, np.maximum(
-n, np.ma.masked_invalid(np.round(yp))
)).astype('i').filled(0) + n
return a[jp, ip]
def EMG(x, alpha, x0, sigma, mu, beta, return_parts=False):
"""
Eq 1 of Goldberg et al.[1] as described
e = 1/x0 exp(mu/x0 + sigma^2 / (2 x0^2) - x / x0)
G = F((x - mu)/sigma - sigma/x_o)
f = e * G
y = alpha f + beta
Note: F is the Gaussian cumulative distribution function"
[1] Goldberg et al. doi: 10.5194/acp-22-10875-2022
Arguments
---------
x : array
x0 is the e-folding distance downwind, representing the length scale
of the NO2 decay
alpha : float
Mass distributed across distance
sigma : float
sigma is the standard deviation of the Gaussian function, representing
the Gaussian smoothing length scale
mu : float
mu is the location of the apparent source relative to the assumed
pollution source center
beta : float
Background line density
return_parts : bool
If False (default), return final output.
If True, return parts for further analysis
Returns
-------
out : array or tuple
If False (default), return y (alpha * e * G + beta)
If True, return alpha, e, G, beta and alpha * e * G + beta
Example
-------
import numpy as np
from emgsatemis import EMG
# Approximation of data in Goldberg[1] Figure 8a inset
# g/y y/h * h molNOx/gNOx NO2/NOx [=] mol
approx_alpha = 62e9 / 8760 * 1.7 / 46 / 1.32
approx_x0 = 20e3 # iteratively identified
exparams = {
'alpha': approx_alpha, 'x0': approx_x0,
'sigma': 28e3, 'mu': -8000.,
'beta': 2.2,
}
x = np.arange(-75e3, 105e3, 5000)
y = EMG(x, **exparams)
print(np.array2string(np.round(x / 1000, 1).astype('i')))
print(np.array2string(np.round(y, 1)))
# [-75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10
# 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100]
# [2.3 2.3 2.3 2.4 2.5 2.6 2.7 2.9 3.1 3.3 3.6 3.8 4.1 4.3 4.4 4.5 4.6 4.6
# 4.6 4.4 4.3 4.1 3.9 3.7 3.5 3.3 3.1 2.9 2.8 2.7 2.6 2.5 2.4 2.4 2.4 2.3]
"""
import numpy as np
from scipy.stats import norm
e = 1 / x0 * np.exp((mu / x0) + (sigma**2 / (2 * x0**2)) - x / x0)
# e = np.where(x < mu, 0, e)
G = norm.cdf((x - mu) / sigma - sigma / x0)
pdf = e * G
# precisely equal to scipy equivalent
# from scipy.stats import exponnorm
# pdf = exponnorm(K=x0 / sigma, loc=mu, scale=sigma).pdf(x)
out = alpha * pdf + beta
if return_parts:
return alpha, e, G, beta, out
else:
return out
def fitlinedens(y, dx, verbose=0, addyhat=True):
"""
Arguments
---------
y : np.array
shape (n,) is the line density assuming focal point at n // 2
dx : float
Width in length units (consistent with units of b, eg, b [=] mol/m2
then dx [=] m)
addyhat : bool
If True (default), add the fit line (yhat) to the output
Returns
-------
out : dict
x, y, params (from fitlinedens), and (optionally) yhat
Example
-------
from emgsatemis import fitlinedens
from scipy.stats import exponnorm
# Approximation of data in Goldberg[1] Figure 8a inset
# g/y y/h * h molNOx/gNOx NO2/NOx [=] mol
approx_alpha = 62e9 / 8760 * 1.7 / 46 / 1.32
approx_x0 = 20e3 # iteratively identified
alpha = 198155
x0 = 20e3
sigma = 28e3
mu = -8000.
beta = 2.2
x = np.arange(-75e3, 156e3, 5e3)
f = exponnorm(K=approx_x0 / sigma, loc=mu, scale=sigma).pdf(x)
y = approx_alpha * f + beta
yerr = np.random.normal(size=y.size) * 0.1
yfit = emgsatemis.fitlinedens(y, dx=1e3)
print(np.array2string(np.round(x / 1000, 1).astype('i')))
print(np.array2string(np.round(y, 1)))
print(np.array2string(np.round(yfit['yhat'], 1)))
# [-75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10
# 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100
# 105 110 115 120 125 130 135 140 145 150 155]
# [2.3 2.3 2.3 2.4 2.5 2.6 2.7 2.9 3.1 3.3 3.6 3.8 4.1 4.3 4.4 4.5 4.6 4.6
# 4.6 4.4 4.3 4.1 3.9 3.7 3.5 3.3 3.1 2.9 2.8 2.7 2.6 2.5 2.4 2.4 2.4 2.3
# 2.3 2.3 2.3 2.2 2.2 2.2 2.2 2.2 2.2 2.2 2.2]
# [2.3 2.3 2.3 2.4 2.5 2.6 2.7 2.9 3.1 3.3 3.6 3.8 4.1 4.3 4.4 4.5 4.6 4.6
# 4.6 4.4 4.3 4.1 3.9 3.7 3.5 3.3 3.1 2.9 2.8 2.7 2.6 2.5 2.4 2.4 2.4 2.3
# 2.3 2.3 2.3 2.2 2.2 2.2 2.2 2.2 2.2 2.2 2.2]
"""
from scipy.optimize import curve_fit
import numpy as np
ny = y.size
assert (ny % 2 == 1)
n = ny // 2
x = (np.arange(ny)[:] - n) * dx
# fitting requires a reasonable starting guess
dist = np.abs(x[-1] - x[0])
sigma0 = dist / 3
x00 = sigma0 / 2
mu0 = -2 * np.float32(dx)
pdf_max = EMG(x, alpha=1, x0=x00, sigma=sigma0, mu=mu0, beta=0).max()
# should this be divided by (pdf_max * dx)?
alpha0 = (y.max() - y.min()) / pdf_max
beta0 = y.min()
guess = np.array([alpha0, x00, sigma0, mu0, beta0], dtype='f')
if verbose > 0:
print('p0', dict(zip('alpha x0 sigma mu beta'.split(), guess)))
lb = np.array([1e-1, 1e-1, 1e-1, x.min(), 0])
ub = np.array([1e10, x.max() * 4, x.max() * 4, x.max(), y.max()])
pfit, pcov = curve_fit(
EMG, x, y, p0=guess, bounds=(lb, ub), check_finite=False
)
pfit = dict(zip('alpha x0 sigma mu beta'.split(), pfit))
out = {'x': x, 'y': y, 'params': pfit}
if addyhat:
out['yhat'] = EMG(x, **pfit)
return out
def rotate_and_fitlinedens(
orgs, us, vs, degrees=False, dx=1, acrossplumen=None, verbose=0
):
"""
Arguments
---------
orgs : array-like
3d (t, row, col) where row=col
us : array-like
1d (t,) U component of wind
vs : array-like
1d (t,) V component of wind
degrees : bool
Perform calculations of wind direction in degrees
Returns
-------
fitresult : dict
yhat = best fit function
ws = wind speeds
wd = wind directions
wsmean = mean of wind speeds
rots = rotated rasters
"""
import numpy as np
rots = []
wss, wds = to_polar(us, vs, degrees=degrees)
for ws, wd, a in zip(wss, wds, orgs):
b = rotate_array(a, -wd, degrees=degrees)
rots.append(b)
rots = np.ma.stack(rots, axis=0)
rot = rots.mean(0)
ny, nx = rot.shape
assert (ny == nx)
n = ny // 2
if acrossplumen is None:
acrossplumen = n - 1
sliceacross = slice(
max(0, n - acrossplumen), min(ny, n + acrossplumen + 1)
)
y = rot[sliceacross, :].sum(0) * dx
out = fitlinedens(y, dx=dx, verbose=verbose)
out['yhat'] = EMG(out['x'], **out['params'])
out['ws'] = wss
out['wd'] = wds
out['wsmean'] = wss.mean()
out['rots'] = rots
out['tau_s'] = out['params']['x0'] / out['wsmean']
out['emis_ps'] = 1.32 * out['params']['alpha'] / out['tau_s']
return out
def get_gfs(date):
"""
Get file with u- and v-component of wind from 0.5 degree GFS simulation at
100m from the 12Z initialization at 6 hour forecast (18Z, ~13LST)
Arguments
---------
date : pandas date-like
Target date for file
Returns
-------
u, v : float, float
u- and v-component of wind from GFS 0.5 degree
Example
-------
gfsf = emgsatemis.get_gfs('2019-01-01')
print(*list(gfsf))
# Temperature_surface, u-component_of_wind_height_above_ground, v-comp...
# Plot US map with Surface Temperature and wind barbs
qm = gfsf['Temperature_surface'].plot()
thin = 5
x = gfsf['lon'][::thin]
y = gfsf['lat'][::thin]
u = gfsf['u-component_of_wind_height_above_ground'][0, 0, ::thin, ::thin]
v = gfsf['v-component_of_wind_height_above_ground'][0, 0, ::thin, ::thin]
qm.axes.quiver(x, y, u, v, scale=200)
#qm.axes.set(xlim=(255, 265), ylim=(25, 30))
"""
import os
import numpy as np
import xarray as xr
import pandas as pd
date = pd.to_datetime(date)
timeslice = '[0:1:0]'
latslice = '[40:1:140]'
lonslice = '[420:1:620]'
urlroot = 'https://www.ncei.noaa.gov/thredds/dodsC/'
groot = 'model-gfs-g4-anl-files-old'
filepath = f'{groot}/{date:%Y%m/%Y%m%d/gfsanl_4_%Y%m%d_1200_006.grb2}'
localpath = f'{filepath}.nc'
ukey = 'u-component_of_wind_height_above_ground'
vkey = 'v-component_of_wind_height_above_ground'
if not os.path.exists(localpath):
print(f'Retrieving {localpath}', end='...')
os.makedirs(os.path.dirname(localpath), exist_ok=True)
url = f'{urlroot}/{filepath}?{ukey}'
tmpf = xr.open_dataset(url)
hdim = tmpf.variables[ukey].dims[1]
url = (
f'{urlroot}/{filepath}?{hdim}'
)
tmpf = xr.open_dataset(url)
# find 100m agl
hidx = np.abs(tmpf[hdim][:].data - 100).argmin()
aglslice = f'[{hidx}:1:{hidx}]'
uvslice = timeslice + aglslice + latslice + lonslice
url = (
f'{urlroot}/{filepath}'
f'?lat{latslice},lon{lonslice},reftime,time{timeslice},'
f'{hdim}{aglslice},{ukey}{uvslice},{vkey}{uvslice}'
f',Temperature_surface{timeslice}{latslice}{lonslice}'
)
f = xr.open_dataset(url).sel(lat=slice(None, None, -1))
f.coords['lon'] = ('lon',), f.coords['lon'].data, f.lon.attrs
# make sure the result is 100m agl
assert f[hdim][0] == 100
f.to_netcdf(localpath)
print()
return xr.open_dataset(localpath)
def get_gfs_uv(date, lon, lat, retfile=False):
"""
Get u- and v-component of wind from 0.5 degree GFS simulation at 100m
from the 12Z initialization at 6 hour forecast (18Z, ~13LST)
Arguments
---------
date : pandas date-like
Target date for file
lon : float
degrees east of 0E
lat : float
degrees north of 0N
Returns
-------
u, v : float, float
u- and v-component of wind from GFS 0.5 degree
Example
-------
u, v = emgsatemis.get_gfs_uv('2019-01-01', -95.1, 29.720621)
print(u, v)
# -3.3206224 -5.3561645
"""
lon = lon % 360
gfsf = get_gfs(date)
llf = gfsf.sel(lon=lon, lat=lat, method='nearest')
u = float(llf['u-component_of_wind_height_above_ground'][0, 0].values)
v = float(llf['v-component_of_wind_height_above_ground'][0, 0].values)
return u, v
def squareit(f):
"""
squareit takes a file and returns a square version.
rotate_and_fitlinedens requires a square domain with an odd number of
ROW and COL.
Arguments
---------
f : xr.Dataset
dataset with ROW and COL dimensions. The shortest of the two will
be used to create an odd number square.
Returns
-------
sqf : xr.Dataset
f with ROW and COL zoomed in to the smallest common odd number.
"""
nr = f.sizes['ROW']
nc = f.sizes['COL']
n = min(nr, nc)
if n % 2 == 0:
n = n - 1
cb = (nc - n) // 2
rb = (nr - n) // 2
outf = f.isel(ROW=slice(rb, rb + n), COL=slice(cb, cb + n))
return outf
def get_tropomi(
lockey, lon, lat, dates, bboxwidth=3, gdnam='4US1', met='hrrr', verbose=0
):
"""
Uses pyrsig to download TropOMI around lon/lat and fit exponentially
modified gaussian to the rotated plume.
Arguments
---------
lockey : str
Name of the location (used for caching data)
lon : float
Longitude center of the location extract
lat : float
Latitude center of the location extract
dates : list
Iterable of pandas.Datatime
bboxwidth : float
Width of bounding box in degrees for initial retrieval of TropOMI
gdnam : str
Name of grid to use 12US1, 4US1, HRRR3K, 1US1 (initally covers large area, but
fit to actual domain)
met : str
hrrr or gfs0.5. hrrr uses wind components (u/v) from HRRR 3km
downloaded using pyrsig. gfs0.5 uses 100m from GFS 0.5 degree model
that is downloaded from NCEI over the US and then sampled.
verbose : int
Level of verbosity for debuging
Returns
-------
f : xarray.Dataset
Data variables include:
* DAILY_NO2: Tropospheric Vertical Column Density from TropOMI
* us, vs : u- and v-component of wind from HRRR
"""
import pyrsig
import numpy as np
import xarray as xr
import pandas as pd
tkey = 'tropomi.offl.no2.nitrogendioxide_tropospheric_column'
tropfs = []
dg = np.array([-1, -1, 1, 1]) * bboxwidth / 2.
gbbox = np.array([lon, lat, lon, lat]) + dg
if verbose > 0:
print('gridbbox', gbbox)
gapi = pyrsig.RsigApi(
bbox=gbbox, workdir=lockey, gridfit=True, grid_kw=gdnam
)
gapi.grid_kw['REGRID_AGGREGATE'] = 'daily'
if verbose > 0:
print('windbbox', gbbox)
for date in dates:
print(f'\rTropOMI: {date:%F}', end='.', flush=True)
try:
tmpf = gapi.to_ioapi(tkey, bdate=date, verbose=-100)
except Exception as e:
if e.errno == -51:
print('Empty file', end='.', flush=True)
else:
print(str(e), end='.', flush=True)
continue
tropfs.append(tmpf)
print()
if verbose > 0:
print('stacking files')
tropf = xr.concat(tropfs, dim='TSTEP')
if verbose > 0:
print('converting to mol/m2')
# use mol/m2 for smaller numbers
tropf['DAILY_NO2'] = tropf['DAILY_NO2'].where(
tropf['DAILY_NO2'] > -9e36
) / 6.022e23 * 1e4
tropf['DAILY_NO2'].attrs.update(units='mol/m2')
metf = get_met(lockey, lon, lat, pd.to_datetime(tropf.TSTEP))
tropf['us'] = metf['us']
tropf['vs'] = metf['vs']
tropf.attrs['met'] = met
tropf.attrs['lockey'] = lockey
now = pd.to_datetime('now', utc=True).strftime('%Y-%m-%dT%HZ')
desc = f"TropOMI NO2 from RSIG and winds from {met} extracted using"
desc += " emgsat.get_tropomi."
hist = tropf.attrs.get('history', '')
hist += f"{now}: emgsatemis v{__version__} get_tropomi;"
ref = "Goldberg et al. doi: 10.5194/acp-22-10875-2022 (and references"
ref += " therein)"
tropf.attrs['title'] = 'TropOMI with Winds'
tropf.attrs['history'] = hist
tropf.attrs['institution'] = 'unknown'
tropf.attrs['source'] = 'retrieved thru pyrsig via emgsatemis.get_tropomi;'
tropf.attrs['comment'] = desc
tropf.attrs['references'] = ref
return tropf
def get_met(workdir, lon, lat, dates, met='hrrr', verbose=0):
"""
Retrieve 80m or 100m winds from met data source (hrrr or gfs0.5).
Arguments
---------
workdir : str
Used as the folder to save out retreived winds
lon : float
Longitude (degrees_east) of winds to retrieve
lat : float
Latitude (degrees_north) of winds to retrieve
dates : iterable
List or series of dates to process
Results
-------
windf : xarray.Dataset
Dataset with u-component winds (us) and v-component winds (vs) from
either gfs0.5 (wind_100m) or hrrr (wind_80m)
"""
import numpy as np
import xarray as xr
import pyrsig
import pandas as pd
outf = xr.Dataset(
coords=dict(TSTEP=dates),
attrs=dict(source=met),
)
uvpfx = {'hrrr': 'wind_80m', 'gfs0.5': 'wind_100m'}[met]
_z = np.zeros(dates.shape, dtype='f')
usattrs = dict(units='m/s', long_name=f'{uvpfx}_u')
outf['us'] = ('TSTEP',), _z * np.nan, usattrs
vsattrs = dict(units='m/s', long_name=f'{uvpfx}_v')
outf['vs'] = ('TSTEP',), _z * np.nan, vsattrs
pbbox = [lon - .025, lat - .025, lon + .025, lat + .025]
papi = pyrsig.RsigApi(bbox=pbbox, workdir=workdir)
hkey = 'hrrr.wind_80m'
s5p_utc = pd.to_timedelta(13 - lon / 15, unit='h')
s5p_utc_b = s5p_utc - pd.to_timedelta('10800s')
s5p_utc_e = s5p_utc + pd.to_timedelta('1799s')
print(met, end=':')
for i, date in enumerate(dates):
print(f'\r{met}: {date:%F}', end='.', flush=True)
if met == 'hrrr':
hbdate = (date + s5p_utc_b)
hedate = (date + s5p_utc_e)
s5pquery = f'time > "{hbdate}Z" and time < "{hedate}Z"'
if verbose > 0:
print('wind-interval', s5pquery)
hrrrdf = papi.to_dataframe(
hkey, bdate=date, backend='xdr', unit_keys=False,
parse_dates=True, verbose=-100
)
uvdf = hrrrdf.query(s5pquery).mean(numeric_only=True)
outf['us'][i] = uvdf['wind_80m_u']
outf['vs'][i] = uvdf['wind_80m_v']
elif met == 'gfs0.5':
u, v = get_gfs_uv(date, lon, lat)
outf['us'][i] = u
outf['vs'][i] = v
elif met is not None:
raise KeyError(f'met ({met}) must be hrrr or gfs0.5')
if verbose > 0:
print('adding winds')
print()
return outf
def fit_tropomi(tropf, nacross, verbose=0):
"""
Arguments
---------
tropf : xr.Dataset
Data variables must include:
* DAILY_NO2: Tropospheric Vertical Column Density from TropOMI
* us, vs : u- and v-component of wind
nacross : int
When fitting, use the sum across a subset of rotated space.
verbose : int
Level of verbosity for debuging
Returns
-------
f : xarray.Dataset
* ROT_DAILY_NO2: same as DAILY_NO2, but rotated by daily wind dir
* ws, wd : wind speed (ws) and direction (wd) from u/v
* fit_x : distance from center pixel [m]
* fit_y : line density (mol/m) integrated from ROT_DAILY_NO2
* fit_yhat : best fit of fit_y (attrs have fit params)
* tau : lifetime of NO2 (x0 / ws) [s]
* emis : emissions estimate from alpha / tau [mol/s]
"""
import numpy as np
import pandas as pd
tropf = squareit(tropf)
if verbose > 0:
print('Rotating and fitting')
no2 = np.ma.masked_invalid(tropf['DAILY_NO2'][:, 0].data)
dx = tropf.XCELL
results = rotate_and_fitlinedens(
no2, tropf['us'], tropf['vs'],
degrees=True, dx=dx, acrossplumen=nacross
)
if verbose > 0:
print('Preparing output')
tropf['ROT_DAILY_NO2'] = (
tropf['DAILY_NO2'].dims, results['rots'][:, None], dict(units='mol/m2')
)
tropf['wd'] = ('TSTEP',), results['wd'], dict(units='degrees')
tropf['ws'] = ('TSTEP',), results['ws'], dict(units='m/s')
tropf['fit_x'] = ('COL',), results['x'], dict(units='m')
tropf['fit_y'] = ('COL',), results['y'], dict(units='mol/m')
tropf['fit_yhat'] = ('COL',), results['yhat'], dict(units='mol/m')
tropf['fit_yhat'].attrs.update(results['params'])
tropf['emis'] = (), results['emis_ps'], dict(units='molNO2/s')
tropf['tau'] = (), results['tau_s'], dict(units='s')
met = tropf.attrs.get('met', 'unknown')
now = pd.to_datetime('now', utc=True).strftime('%Y-%m-%dT%HZ')
desc = tropf.attrs.get('comment', '')
desc += "Rotated plume and best fit line added using emgsatemis."
desc += "fit_tropomi."
hist = tropf.attrs.get('history', '')
hist += f"{now}: emgsatemis v{__version__} fit_tropomi;"
ref = "Goldberg et al. doi: 10.5194/acp-22-10875-2022 (and references"
ref += " therein)"
tropf.attrs['title'] = 'TropOMI with Winds, Rotated Plume, and fit line.'
tropf.attrs['history'] += hist
source = tropf.attrs.get('source' , '')
tropf.attrs['source'] = 'fit via emgsatemis.fit_tropomi;'
tropf.attrs['comment'] = desc
tropf.attrs['references'] = ref
return tropf
def plot_fit(tropf, axx=None, lax=None, nacross=21):
"""
Create a plot with detailed information from fit_tropom:
1. Map of unrotated NO2.
2. 2-d plot of rotated NO2
3. Line-plot showing the 1-d integrated mass pdf and EMG fit line
Arguments
---------
tropf : xarray.Dataset
Must have DAILY_NO2 and ROT_DAILY_NO2 (from fit_tropomi)
axx : array
Array of axes shape (2,) where the first will be used for a map of
unrotated data, the second will be used for the rotated plume plot.
If None, fig, axx = plt.subplots(1, 2, figsize=(16, 6))
lax : matplotlib.axes.Axes
Axes to use for the line plot. If None, lax will be inset in axx[1]
Returns
-------
fig : matplotlib.figure.Figure
Figure object from axx[1]
"""
import matplotlib.pyplot as plt
if axx is None:
fig, axx = plt.subplots(1, 2, figsize=(16, 6))
qm = tropf['DAILY_NO2'].mean(('TSTEP', 'LAY')).plot(ax=axx[0])
lockey = tropf.attrs.get('lockey', 'Unknown')
x = tropf.COL.values
c = x.mean()
r = tropf.ROW.values.mean()
n = x.size
x0inset = 20
dxinset = 50
y0inset = tropf.ROW.min() + 2
ax = axx[1]
tropf['ROT_DAILY_NO2'].mean(('TSTEP', 'LAY')).plot(ax=ax, norm=qm.norm)
ax.axhline(r + nacross / 2, color='k', linestyle=':')
ax.axhline(r - nacross / 2, color='k', linestyle=':')
ax.scatter(c, r, marker='+', color='r')
axx[0].scatter(c, r, marker='+', color='r')
E = tropf['emis'] * 3600 * 8760 * 46e-9 # mol/s s/h h/y Gg/mol [=] Gg/yr
tauh = tropf['tau'] / 3600
wsmean = tropf["ws"].mean()
label = f'E={E:.0f}Gg/y;$\\tau$={tauh:.2}h;ws={wsmean:.2}'
fitlabel = r"""$\alpha$={alpha:.0f};x0={x0:.0f};$\sigma$={sigma:.0f}
$\beta$={beta:.2f};$\mu$={mu:.1f}""".format(**tropf.fit_yhat.attrs)
ax.text(
0.05, 0.75, label + '\n' + fitlabel, backgroundcolor='w',
transform=ax.transAxes
)
lrect = [
(x0inset - x.min() + 0.5) / n,
(y0inset - tropf.ROW.min()) / n, dxinset / n, .2
]
if lax is None:
lax = ax.inset_axes(lrect, frameon=True)
lax.plot(x, tropf['fit_y'], label='y', color='k')
lax.plot(x, tropf['fit_yhat'], label='fit', color='r')
lax.legend(loc='upper right', ncol=2)
xticks = ax.get_xticks()
lax.set(xticks=xticks, xticklabels=['' for xt in xticks])
lax.set(
ylim=(0, 10), facecolor='w', ylabel='mol/m',
xlim=(x0inset, x0inset + dxinset)
)
axx[0].set_title(f'{lockey} Original Mean')
axx[1].set_title(f'{lockey} {tropf.met} Rotated Mean')
xf = plt.np.array([-8, dxinset + 8]) + x0inset
yf = plt.np.array([y0inset, y0inset])
xf = [lrect[0] - 0.1, sum(lrect[::2]) + 0.1]
yf = [lrect[1] - 0.1, sum(lrect[1::2]) + 0.02]
ax.fill_between(
x=xf, y1=yf[:1]*2, y2=yf[-1:]*2, color='w', zorder=2,
transform=ax.transAxes
)
return ax.figure
if __name__ == '__main__':
import argparse
import pandas as pd
prsr = argparse.ArgumentParser(
formatter_class=argparse.RawDescriptionHelpFormatter
)
aa = prsr.add_argument
aa('lockey', help='Short name describing location')
aa('lon', help='degrees east of 0E', type=float)
aa('lat', help='degrees north of equator', type=float)
aa('date', help='Date to process')
aa('edate', default=None, help='If provided, process range (date, edate)')
aa('--nacross', default=21, help='Number of grid boxes')
aa('--bboxwidth', default=3, help='Width (deg) used to query tropomi')
aa('--gdnam', default='4US1', help='Grid to fit tropomi')
aa('--met', default='hrrr', help='Met model (hrrr or gfs0.5)')
aa('--verbose', default=0, help='Level of verbosity')
aa('--outpath', default=None, help='Output path')
prsr.epilog = """
Example
-------
python emgsatemis.py Houston -95.1 29.720621 2019-04-01 2019-09-30
"""
args = prsr.parse_args()
if args.edate is None:
args.edate = args.date
args.date = pd.to_datetime(args.date)
args.edate = pd.to_datetime(args.edate)
opts = {
k: v for k, v in vars(args).items()
if k not in ('date', 'edate', 'outpath', 'nacross', 'met')
}
opts['dates'] = pd.date_range(args.date, args.edate, freq='1d')
tropf = get_tropomi(met=args.met, **opts)
tropf = fit_tropomi(tropf, nacross=args.nacross, verbose=args.verbose)
if args.outpath is None:
args.outpath = f'{args.lockey}_{args.date:%F}_{args.edate:%F}.nc'
tropf.to_netcdf(args.outpath)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment