PYSTILT

PYSTILT is a Python implementation of the STILT Lagrangian atmospheric transport model. It runs backward trajectories with HYSPLIT and computes receptor footprints that map where upwind surface fluxes influence a measurement.

The project is in alpha and focused on a unified execution model that works for one-off runs, large batch runs, and streaming queue workers.

Status

PYSTILT is in alpha development. No backward compatibility guarantees before v1.0.

The core transport is stable: HYSPLIT execution, trajectory and footprint generation, numerical R-STILT parity, and the local and SLURM execution paths are all exercised by the test suite. The public API may change while the package settles.

Choose a workflow

• One-off transport runs for local analysis and notebooks: use Model.run() or stilt run.
• Queue-backed batch or service runs for cloud execution: use Model.register(), stilt register, stilt pull-worker, and stilt serve with a PostgreSQL work queue configured via PYSTILT_DB_URL. Slurm needs no database: stilt run --backend slurm.
• Observation-driven workflows for science-facing code: use stilt.observations to turn normalized observations into Receptor objects before feeding them into the same runtime.

Roadmap

PYSTILT draws design and science inspiration from two sister projects: X-STILT for column and satellite science workflows, and stiltctl for cloud-native execution patterns. The tables below track what has been absorbed and what remains in scope. See the full roadmap for more details.

Execution and orchestration (from stiltctl)

Column and satellite science (from X-STILT)

Full X-STILT feature parity is not a goal. PYSTILT absorbs X-STILT's observation-layer design and column-weighting concepts without trying to replicate every script.

Installation

pip install pystilt

For Slurm, Kubernetes, projections, plotting, and cloud object stores:

pip install "pystilt[complete]"

Quickstart: one-off run

Define a receptor, configure meteorology and footprint grid, then run:

import pandas as pd
import stilt

receptor = stilt.Receptor(
    time=pd.Timestamp("2023-07-15 18:00", tz="UTC"),
    latitude=40.766,
    longitude=-111.848,
    altitude=10,
)

model = stilt.Model(
    project="./my_project",
    receptors=[receptor],
    config=stilt.ModelConfig(
        n_hours=-24,
        numpar=100,
        mets={
            "hrrr": stilt.MetConfig(
                directory="/data/hrrr",
                file_format="hrrr_%Y%m%d.arl",
                file_tres="1h",
            )
        },
        footprints={
            "default": stilt.FootprintConfig(
                grid=stilt.Grid(
                    xmin=-113.0,
                    xmax=-110.5,
                    ymin=40.0,
                    ymax=42.0,
                    xres=0.01,
                    yres=0.01,
                )
            )
        },
    ),
)

handle = model.run()
handle.wait()

sim = list(model.simulations.values())[0]
traj = sim.trajectories
foot = sim.get_footprint("default")

Quickstart: queue/service runtime

# Queue workers require a PostgreSQL work queue.
export PYSTILT_DB_URL=postgresql://user:pass@host:5432/pystilt

# Initialize project files (config.yaml and receptors.csv)
stilt init ./my_project

# Run with local workers (blocks until complete)
stilt run ./my_project --backend local --n-workers 8

# Persist inputs and enqueue every simulation (receptors x mets)
stilt register ./my_project

# Drain queue from worker processes (batch mode)
stilt pull-worker ./my_project

# Long-lived queue workers (streaming mode)
stilt serve ./my_project

# Check project status
stilt status ./my_project

The same queue model is available in Python:

import stilt
from stilt.execution import pull_simulations

model = stilt.Model(project="./my_project")
model.register()
pull_simulations(model, follow=False)  # batch mode
print(model.status())

Workers claim simulations from the queue and record done/failed there; whether outputs exist is always read from the project itself.

Quickstart: observation layer

stilt.observations sits above Receptor for measurements and retrievals that are not already receptors. Your reader produces Observation objects; PYSTILT groups, selects, and turns them into receptors:

import stilt
from stilt.observations import Observation, build_point_receptor, group_by_overpass

observations = [
    Observation(sensor="tower", species="co2", time="2023-01-01 12:00:00",
                latitude=40.77, longitude=-111.85, altitude=30.0, observation_id="tower-001"),
    Observation(sensor="tower", species="co2", time="2023-01-01 12:05:00",
                latitude=40.78, longitude=-111.84, altitude=30.0, observation_id="tower-002"),
]

model = stilt.Model(project="./my_project")  # existing project config on disk
for scene in group_by_overpass(observations):
    model.register(receptors=scene.receptors(build_point_receptor))

scene.receptors() takes any observation-to-receptor callable, so a custom instrument is a reader plus, when needed, your own builder and transform. See the observations guide.

This layer currently focuses on:

• normalized Observation and Scene objects
• horizontal, viewing, and line-of-sight geometry
• overpass scenes, sounding selection, jitter
• observation-to-receptor conversion

See docs/advanced/observations.rst for the intended workflow boundary.

Particle transforms

Per-footprint transforms rescale each particle's influence before the footprint is rasterized. Declare them in config, or pass them in Python:

footprints:
  column:
    grid: slv
    transforms:
      - kind: averaging_kernel
        levels: [0.0, 1000.0, 2000.0]
        values: [1.0, 0.8, 0.5]
      - kind: pressure_weighting      # derived from the particles, X-STILT style
      - kind: first_order_lifetime
        lifetime_hours: 4.0
      - kind: mypkg.transforms.MyWeighting   # your own pydantic class with apply()
        some_field: 3

A transform is any object with apply(particles, context). See the transforms guide for the column-weighting science and for writing your own.

Accessing results

import pandas as pd

for sim in model.simulations.values():
    traj = sim.trajectories
    foot = sim.get_footprint("default")

# Load footprints across all matching simulations
footprints = model.footprints["default"].load(
    time_range=("2023-01-01", "2023-01-31")
)

coords = [(-111.9, 40.7), (-111.8, 40.8)]
time_bins = pd.interval_range(
    start=pd.Timestamp("2023-01-01 00:00", tz="UTC"),
    end=pd.Timestamp("2023-01-02 00:00", tz="UTC"),
    freq="1h",
)

for footprint in footprints:
    hourly = footprint.aggregate(target=coords, time_bins=time_bins)

If a footprint is tracked as complete-empty, no NetCDF file is expected for that footprint. The model APIs treat it as a successful terminal outcome while skipping missing file loads.

R-STILT parity

PYSTILT footprints match the uataq/stilt R implementation on numerical values at rtol=1e-7 per cell, validated by end-to-end fidelity scenarios against a pinned upstream commit. NetCDF output is not byte-compatible with R-STILT; it should be read as generic CF-1.8 NetCDF. See STILT-R.md for more details.

Documentation

Full documentation is available at https://jmineau.github.io/PYSTILT/

Contributing

Contributions are welcome! Please see CONTRIBUTING.md for guidelines.

License

This project is licensed under the MIT License - see the LICENSE file for details.

Author

James Mineau - jmineau