Plume Background#

Background takes the background from a model field at the trajectory endpoints. For a satellite swath, the swath itself holds a background. Soundings that the city’s plume did not reach saw clean air at the same time, with the same instrument and the same biases. To find those soundings, run particles forward from the city over the hours before the overpass, see where they are when the satellite passes, and draw the plume around them. plume_polygon() draws the plume and plume_background() picks the background soundings. This follows X-STILT’s fit.kde.plume and calc.bg.upwind (method M3 in Wu et al., 2018).

Forward runs#

A forward run is a run with a positive n_hours. Everything else stays the same. The receptor is where and when the particles are released, and the particle table records where they go.

A plume is what the whole city puts out, so spread the release over a box around the city, just above the surface. Repeat it every half hour over the hours before the overpass, so that air of every age is included. X-STILT’s defaults are a 0.3° box, 10 m above ground, and releases from ten hours before the overpass up to the overpass, each followed twelve hours forward with 1000 particles. In PYSTILT each release is a MultiPointReceptor, and jitter_points() spreads its points over the box:

import pandas as pd
import stilt
from stilt.observations import jitter_points
from stilt.receptors import MultiPointReceptor

site = (-111.9, 40.75)                       # Salt Lake City
overpass = pd.Timestamp("2023-10-19 19:42")  # from the sounding times
half = 0.15
box = [(site[0] - half, site[1] - half), (site[0] + half, site[1] - half),
       (site[0] + half, site[1] + half), (site[0] - half, site[1] + half)]
points = jitter_points(box, 100)

receptors = [
    MultiPointReceptor(t, [p[0] for p in points], [p[1] for p in points],
                       [10.0] * len(points))
    for t in pd.date_range(overpass - pd.Timedelta(hours=10), overpass, freq="30min")
]

forward = stilt.Model(
    project="./forward_12h",     # keep forward runs in a project of their own
    receptors=receptors,
    mets={"hrrr": met},          # as in the quickstart
    n_hours=12,                  # positive: forward in time
    numpar=1000,
)
forward.run()

No footprint is needed. The plume only needs the particle tables, which hold each particle’s position at every output step along with its datetime.

The plume#

Pool the particle rows that fall within a few minutes of the overpass, across all the forward runs, and outline them:

from stilt.observations import plume_polygon

window = (overpass - pd.Timedelta(minutes=3), overpass + pd.Timedelta(minutes=3))
rows = []
for sim in forward.simulations:
    p = sim.trajectories.data
    rows.append(p[p["datetime"].between(*window)])
particles = pd.concat(rows)

plume = plume_polygon(particles["long"], particles["lati"])
plume.polygon      # shapely Polygon in longitude and latitude
plume.density      # the kernel density it was cut from, (lat, lon), max 1

The plume is where the 2-D kernel density of these positions is at least a tenth of its maximum. threshold sets that fraction. A higher value gives a tighter plume, and a lower value a wider one. bandwidth sets the smoothing, 0.1° in longitude and 0.15° in latitude by default. As in R’s kde2d, the kernel’s standard deviation is a quarter of the bandwidth. If the area above the threshold falls into separate pieces, the largest piece is the plume. The outline follows the cells of the density grid, 100 by 100 over the particles by default, so it looks a little blocky. Raise n for a finer outline.

Plot plume.density with the soundings and the outline before you trust it. A plume that misses the swath, or one that covers all of it, leaves no background to take.

The background#

Pass the soundings of the overpass and the plume to plume_background():

from stilt.observations import plume_background, read_tropomi_ch4

obs = read_tropomi_ch4(path, lon_range=(-113.5, -110.5), lat_range=(39.5, 42.0))
obs = obs[obs["good"]]

bg = plume_background(
    obs["longitude"], obs["latitude"], obs["value"], plume,
    uncertainties=obs["uncertainty"],
)
bg.value          # the background, ppb here
bg.uncertainty    # spread of the background soundings and their retrieval error, in quadrature
bg.sides          # the same for each side of the plume
obs["in_plume"] = bg.in_plume
obs["enhancement"] = obs["value"] - bg.value

The soundings inside the plume are the enhanced ones. PYSTILT draws the bounding box of those soundings and pads it on each side by pad times its size (0.1 by default). The background candidates are the soundings outside the plume within width degrees (0.5 by default) of that box, to the north, south, east, or west. The background is their median.

bg.sides has the count, mean, median, spread, and retrieval error of the soundings on each side. If the sides disagree, one of them is probably downwind of another source. Pass the upwind side, for example side="west", to use that side alone. The drift of the forward particles tells you which side is upwind.

Before any of that, trim drops the out-of-plume soundings above the 90th percentile of their values. This guards against enhanced air that the outline missed, at the cost of a slightly low background when the air is clean. trim=None keeps every sounding.

plume_background() raises an error when no sounding lies inside the plume, because then the overpass did not see the city. Pass only good-quality soundings. Missing values are ignored.

Choices made here#

  • The background comes from the swath, not from a model. It carries the instrument’s bias, which cancels when the enhancement is value - background from the same swath.

  • The outline is built from the density grid cells above the threshold rather than traced as a contour. It matches the grid exactly and never crosses itself. X-STILT traces the contour and repairs broken pieces. Keeping the largest piece does the same job here.

  • You choose the threshold. X-STILT raises it automatically when no low density falls in the release box. Here you look at the density and decide.

  • The side strips cover only the length of the padded box, so a northern strip does not run across the whole swath.

  • Each sounding is a point at its centre. A pixel that straddles the outline counts as inside or outside by where its centre falls.