pyrsig.emiss package

Module contents

pyrsig.emiss.divergence(c, dy, dx, order=4, withdiag=True)[source]

calculate gradient of c in w-e and s-n directions and optionally in the ne-sw and se-nw diagonals

Parameters:

  • c (array) – 2-dimensional array of column densities

  • dy (float or array) – pixel edge length in the y direction

  • dx (float or array) – pixel edge length in the x direction

Returns:

dcdx, dcdy[, dcdr, dcds] – tuple of 2 or 4 divergence

Return type:

tuple

Notes

Adapted from https://github.com/Kang-Sun-CfA/Oversampling_matlab/blob

/v0.2/popy.py#L1711C1-L1732C1

pyrsig.emiss.emg(x, alpha, x0, sigma, mu, beta, return_parts=False)[source]

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

Parameters:

  • 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 – If False (default), return y (alpha * e * G + beta) If True, return alpha, e, G, beta and alpha * e * G + beta

Return type:

array or tuple

Example

import numpy as np from pyrsig.emiss._emg 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]

pyrsig.emiss.fitemg(vcds, us=None, vs=None, degrees=False, dx=1, nxplume=None, verbose=0)[source]
Parameters:

  • vcds (xr.DataArray) – Vertical column densities 3d (t, y, x) where (x, y) must be square (x == y). Units seem to work best when mol/m2. The area unit must match the length unit of dx.

  • us (array-like) – 1d (t,) U component of wind in the x direction at center x/y

  • vs (array-like) – 1d (t,) V component of wind in the y direction at center x/y

  • dx (float) – edge length in units commensurate with vcds (typically m)

  • nxplume (int) – Number of cells to sum in the cross-plume direction after rotation when calculating the line density (mol/m) along the plume. Defaults to the half width minus 1.

  • degrees (bool) – Perform calculations of wind direction in degrees (default=True)

Returns:

fitresult – vldhat: best fit function ws: wind speeds wd: wind directions wsmean: mean of wind speeds rots: rotated rasters

Return type:

dict

pyrsig.emiss.fitemg1d(vld, dx, verbose=0, addvldhat=True)[source]
Parameters:

  • vld (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)

  • addvldhat (bool) – If True (default), add the fit line (vldhat) to the output

Returns:

out – x, vld, params (from fitemg1d), and (optionally) vldhat

Return type:

dict

Example

from pyrsig.emiss import fitemg1d 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 = fitemg1d(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[‘vldhat’], 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]