Particle Weighting (Transforms)#

A footprint is the mean surface influence of a simulation’s particles. A transform changes how much each particle counts before that mean is taken, by rescaling the particle table’s foot column. Column instruments need this (an averaging kernel and a pressure weighting), and so does any species that decays during transport.

A transform is any object with one method:

def apply(self, particles: pd.DataFrame, context: TransformContext) -> pd.DataFrame

It receives the particle table and returns a new one. Transforms are listed per footprint in config.yaml or passed to stilt.Simulation.generate_footprint(), run once in order on the unweighted particles, and are recorded in the footprint’s netCDF so a stored footprint knows how it was weighted.

Built-in transforms#

footprints:
  column:
    grid: slv
    transforms:
      - kind: averaging_kernel
        levels: [0, 500, 1000, 2000]
        values: [1.0, 0.95, 0.8, 0.5]
      - kind: pressure_weighting
      - kind: first_order_lifetime
        lifetime_hours: 4.0
averaging_kernel — stilt.transforms.AveragingKernel

Multiplies each particle’s foot by the kernel interpolated at the particle’s release coordinate. levels are release heights AGL in metres by default; set coordinate: pres when the kernel is on pressure levels (hPa). Fold instrument-specific factors (for example TCCON’s wet-air scaling) into values. Adds an ak_weight column.

pressure_weighting — stilt.transforms.PressureWeighting

Weights each particle by the fraction of the column’s air mass it represents. A column instrument averages over air mass, but HYSPLIT releases column particles uniformly in height, so a plain particle mean over-weights the thin upper layers. Following X-STILT, the pressure weighting function is derived from the particles themselves: a hypsometric curve is fit to their first-step heights and pressures and evaluated at each release height, and each particle carries the slab of air centred on it. Nothing needs to be supplied; pass surface_pressure (hPa) to reference the profile to the retrieval’s surface pressure instead of the fitted value. Requires pres and zagl in varsiwant (both are defaults). Adds xpres (release pressure) and pwf columns.

Two things are worth knowing about the result. The weights sum to the fraction of the atmosphere’s mass the column covers (about 0.3 for a 0–3 km column), not to one: air above the column top cannot be reached by surface fluxes within the back-trajectory, so the footprint is complete and the remaining fraction belongs to the prior profile. And the weighted footprint’s magnitude does not depend on numpar. A column whose bottom sits above the ground still measures the air beneath it, so the lowest particle carries that whole sub-column; start the receptor at or near the surface unless you intend that.

first_order_lifetime — stilt.transforms.FirstOrderLifetime

Decays foot by exp(-age / lifetime), where age is the particle’s transport time from time_column (minutes by default).

The usual satellite or TCCON column weighting is averaging_kernel followed by pressure_weighting. Listing nothing leaves every particle with equal weight, which is standard STILT behaviour.

One kernel per receptor#

A satellite product gives every sounding its own averaging kernel, and an EM27 kernel changes with solar zenith angle. Instead of inline levels and values, name a table in the project:

footprints:
  column:
    grid: slv
    transforms:
      - kind: averaging_kernel
        table: kernels.parquet
        coordinate: pres
      - kind: pressure_weighting

The table has one receptor id per sounding and one level / value row per kernel point (Parquet or CSV). Build it with averaging_kernel_table() from the receptors you registered and the kernels from your product, and write it next to receptors.csv:

from stilt.transforms import averaging_kernel_table

table = averaging_kernel_table(receptors, levels=df.ak_pressure, values=df.ak)
table.to_parquet(model.project.directory / "kernels.parquet")

levels is one array per receptor, or a single array when every kernel shares a grid. When the footprint is generated the transform looks up the receptor’s id in the table, so it works the same in a notebook, under stilt run, and on Slurm and Kubernetes workers, which resolve the relative path against the project root. A receptor with no row is an error.

In Python#

The same classes work directly on a simulation:

from stilt.transforms import AveragingKernel, PressureWeighting

foot = sim.generate_footprint(
    "column",
    config,
    transforms=[AveragingKernel(levels=ak.z, values=ak.values), PressureWeighting()],
)

transforms= runs after config.transforms. To inspect what a transform did, apply it to the particle table yourself:

weighted = PressureWeighting().apply(sim.trajectories.data)
weighted.drop_duplicates("indx")[["xhgt", "xpres", "pwf"]]

Writing your own transform#

Write a pydantic model with the fields you want in YAML and an apply method. Return a new frame; never modify the one you are given.

# mypkg/transforms.py
import numpy as np
import pandas as pd
from pydantic import BaseModel

from stilt import TransformContext
from stilt.transforms import release_coordinate


class BoundaryLayerOnly(BaseModel):
    """Drop influence from particles released above ``max_height`` m AGL."""

    max_height: float = 1500.0

    def apply(self, particles: pd.DataFrame, context: TransformContext) -> pd.DataFrame:
        z = release_coordinate(particles, "xhgt")           # one value per particle
        keep = z.reindex(particles["indx"].to_numpy()) <= self.max_height
        out = particles.copy()
        out["foot"] = np.where(keep.to_numpy(), out["foot"], 0.0)
        return out

Reference it from config.yaml by its import path; the remaining keys are the model’s fields:

footprints:
  bl:
    grid: slv
    transforms:
      - kind: mypkg.transforms.BoundaryLayerOnly
        max_height: 1200

Three rules keep this predictable:

  • Return a copy. Trajectories.data must still be the unweighted table after a footprint is generated. The built-ins all particles.copy() first.

  • Be importable wherever it runs. Workers rebuild the model from config.yaml alone, so mypkg must be installed on the Slurm node or in the container. A footprint whose transform cannot be imported can still be read (the entry becomes an UnresolvedTransform), but a project config naming one fails validation with the import error.

  • Use the context for anything outside the table. context.receptor gives the receptor (its id keys any per-receptor input file), context.is_error says whether this is the error trajectory, and context.store is the project store, whose local_path(key) turns a path relative to the project root into a readable file wherever the worker runs. That is how the averaging_kernel table is found.

A plain class works too (kind imports it and calls cls(**keys)), but it cannot be written back to config or into a footprint’s netCDF, so prefer pydantic for anything declared in YAML. Any object with apply can be passed through the Python transforms= argument.

The science helpers the built-ins are made of are public in stilt.transforms: release_coordinate(), particle_pwf() and ak_weights(). Composing them is usually shorter than re-deriving the weighting.