stilt.particles.transport_error#

stilt.particles.transport_error(particles, error_particles, flux, *, transforms=(), receptor=None, directory=None, levels=20, length_scale=356.0, percentile=1.0, noise_splits=16, background=None)[source]#

Return the transport-error variance of a modeled enhancement (Lin and Gerbig, 2005).

Parameters:

  • particles (DataFrame) – Unperturbed particle table of a receptor (sim.particles).

  • error_particles (DataFrame | Sequence[DataFrame]) – Particle table of a wind-error variant of the same receptor, or one table per realization for a variant with realizations: N (see Examples). A table is one simulation’s particles: one receptor and one realization.

  • flux (DataArray) – Surface flux field (see particles.stilt.enhancement).

  • transforms (Sequence[Any]) – The footprint’s particle transforms (sim.variant.footprint.transforms), applied to both tables so the error is weighted like the footprint.

  • receptor (stilt.receptor.Receptor | None) – The receptor (sim.receptor), for the transforms.

  • directory (str | Path | None) – Where a kernel table’s relative path starts (project.directory).

  • levels (int | Sequence[float]) – Release-height levels to compute statistics on: a number of equal-width bins between the lowest and highest release height, or bin edges in meters. When the particles have no more distinct release heights than levels (a multipoint receptor), each height is a level. A point receptor is one level.

  • length_scale (float | None) – Vertical e-folding length of the error correlation between levels, in meters (X-STILT’s value). None treats the levels as uncorrelated. It has no effect for a point receptor.

  • percentile (float) – In each level, drop particles above this quantile of the enhancement before taking the variance. 1.0 keeps every particle, as Lin and Gerbig do. X-STILT uses 0.99 to limit the few particles that cross a point source. Means always use every particle.

  • noise_splits (int) – Number of random half-splits of the unperturbed particles used to estimate noise. Fewer than 2 skips it and gives NaN.

  • background (DataArray | None) – Background field sampled at each particle’s endpoint (see background()). Wind errors move the endpoints as well as the surface contact, so with a field the statistics are of the modeled mole fraction per particle, as X-STILT computes them.

Return type:

TransportError

Notes

Per level l with n_l of N particles, dvar_l is the change in the variance of the per-particle enhancement under the perturbation, s_l = sign(dvar_l) sqrt(|dvar_l|), and the column variance is Σ_ij w_i w_j s_i s_j exp(-|h_i - h_j| / L) with w_l = n_l / N and h the level heights. For one level it is dvar itself. The transforms are applied to the particles first, so the level statistics are already column-weighted and the level weights are particle counts.

With several error realizations, each level’s mean_err and var_err are averaged over them before dvar is formed, so the perturbed side’s sampling noise falls as 1/sqrt(N). The unperturbed side is the same particles in every realization, so its noise does not fall: with N realizations the null spread of variance is sqrt((1 + 1/N) / 2) times the single-realization noise, which tends to 1/sqrt(2). noise includes that factor. More realizations therefore give at most a sqrt(2) tighter estimate and cannot resolve a case that one realization leaves unresolved.

The signal is the extra spread the perturbation adds. It is small when the wind error decorrelates quickly (HYSPLIT decorrelates it with the distance a particle travels as well as with time), and when turbulence already spreads the particles widely, as on a convective afternoon. Compare a single value with noise. Over many receptors, combine the signed variance values with a median instead of clipping each at zero.

Examples

For the variant hrrr-err with realizations: 4, select one receptor’s simulations and group their particles by realization. Grouping by variant gives one table, which counts as one realization and gives a wrong variance. Without the receptor, each realization’s table holds every receptor’s particles.

>>> from stilt.particles import transport_error
>>> sim = project.simulation(rid, "hrrr")
>>> sims = project.simulations
>>> ensemble = sims[
...     (sims.variant == "hrrr-err") & (sims.receptor == rid)
... ]
>>> by_realization = project.particles(ensemble).groupby(
...     "realization"
... )
>>> error_particles = [rows for _, rows in by_realization]
>>> result = transport_error(sim.particles, error_particles, flux)
>>> result.realizations
4

sim.transport_error(realizations, flux) builds this list from the realizations’ simulations. It also passes transforms, receptor, and directory, which a variant with footprint transforms needs and this call leaves out.