PYSTILT

PYSTILT is a Python version of the STILT atmospheric transport model. It follows air backward in time from a measurement with HYSPLIT and computes a footprint, a map of which upwind surface areas influenced the measurement and by how much.

Status

PYSTILT is alpha software. There are no backward compatibility guarantees before v1.0, and the public API may still change.

The transport core is stable. The test suite covers HYSPLIT runs, trajectories and footprints, agreement with STILT-R, and runs on a local machine and on Slurm.

Installation

pip install pystilt

For plotting, projected grids, aggregation over shapefiles, and cloud features (NOAA downloads, s3:// and gs:// projects, the Postgres work queue, Kubernetes), install everything:

pip install "pystilt[complete]"

Quickstart

Describe a measurement, point PYSTILT at meteorology and a footprint grid, and run:

import stilt

receptor = stilt.PointReceptor(
    time="2023-07-15 18:00",   # UTC
    longitude=-111.848,
    latitude=40.766,
    altitude=10,               # metres above ground
)

model = stilt.Model(
    project="./my_project",
    receptors=[receptor],
    n_hours=-24,
    numpar=100,
    mets={
        "hrrr": {
            "directory": "/data/hrrr",
            "file_format": "%Y%m%d_%H",
            "file_tres": "6h",
        }
    },
    grid={
        "xmin": -113.0, "xmax": -110.5,
        "ymin": 40.0, "ymax": 42.0,
        "xres": 0.01, "yres": 0.01,
    },
)

model.run()   # returns when the run is done

sim = model.simulations[receptor.id, "hrrr"]
traj = sim.trajectories
foot = sim.footprint

Everything is saved in ./my_project. Later you can open it again with stilt.Model(project="./my_project"), and a second run() only runs what is missing.

Command line

stilt init ./my_project          # write a starter config.yaml and receptors.csv
stilt run ./my_project           # run every receptor that is not done yet
stilt run ./my_project --n-workers 8   # the same, in 8 parallel processes
stilt status ./my_project        # what has finished

To run on a Slurm cluster, add an execution section to config.yaml. stilt run then submits a job array and returns. No database is needed.

execution:
  backend: slurm
  n_workers: 200          # array tasks
  account: my-account     # other keys become #SBATCH options
  partition: my-partition
  time: "02:00:00"

Variants

A variant is a named set of settings. Every receptor runs once under each variant. The top-level settings in config.yaml are the defaults, and a variant lists only what it changes. With no variants section there is one variant per met, named after it.

variants:
  hrrr: {}                                  # the defaults
  hrrr-zi08: {ziscale: 0.8}                 # mixed-layer height scaled by 0.8
  hrrr-err:                                 # wind-error ensemble for transport_error
    siguverr: 2.6
    tluverr: 260
    zcoruverr: 450
    horcoruverr: 14
    realizations: 4                         # hrrr-err-0 .. hrrr-err-3
    grid: null                              # particles only, no footprint
  hrrr-ak:                                  # another footprint from the hrrr particles
    from: hrrr
    transforms: [{kind: averaging_kernel, table: kernels.parquet}]

Each receptor under each variant is one simulation, stored in simulations/by-id/<receptor>/<variant>/. Adding a variant runs only the new simulations. PYSTILT refuses to change the settings of a variant that has already run. Give the new settings a new name, or delete the old outputs with stilt rm --variant NAME and run again.

Queue workers

For large batches in the cloud, PYSTILT can hand out receptors from a PostgreSQL work queue. Point PYSTILT_DB_URL at the database, add the receptors to the queue, and start workers:

export PYSTILT_DB_URL=postgresql://user:pass@host:5432/pystilt

stilt register ./my_project      # save the inputs and queue every receptor
stilt pull-worker ./my_project   # run queued receptors until the queue is empty
stilt serve ./my_project         # or keep waiting for new work

The same from Python:

import stilt
from stilt.execution import pull_receptors

model = stilt.Model(project="./my_project")
model.register()
pull_receptors(model)
print(model.status())

The queue only tracks which receptors are pending, done, or failed. Whether a simulation is finished is always decided by its output files in the project.

Column and satellite soundings

Read the soundings into a table with one row per sounding. Group and thin the rows with stilt.observations, then make one receptor per row. Each sounding has its own averaging kernel, so write the kernels to a table in the project:

import stilt
from stilt.observations import group_by_overpass, read_tropomi_ch4
from stilt.transforms import averaging_kernel_table

df = read_tropomi_ch4(path)                      # or read_oco2, read_tccon, or your own reader
df["overpass"] = group_by_overpass(df["time"])   # label rows by overpass

model = stilt.Model(project="./my_project")      # a project with a config.yaml
receptors = [stilt.ColumnReceptor(r.time, r.longitude, r.latitude, 0, 3000) for r in df.itertuples()]
model.register(receptors=receptors)

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

Then name the table in config.yaml:

grid: {xmin: -114.0, xmax: -111.0, ymin: 39.0, ymax: 42.0, xres: 0.01, yres: 0.01}
transforms:
  - kind: averaging_kernel
    table: kernels.parquet
    coordinate: pres
  - kind: pressure_weighting

For slanted lines of sight, use slant_points and Receptor.from_points. pressure_altitudes converts a retrieval's pressure levels to altitudes. See the observations guide and the slant columns guide.

Particle transforms

A transform changes how much each particle counts before the footprint is calculated. List them in config.yaml, as a default or for one variant:

transforms:
  - kind: averaging_kernel
    levels: [0.0, 1000.0, 2000.0]
    values: [1.0, 0.8, 0.5]
  - kind: pressure_weighting
  - kind: first_order_lifetime
    lifetime_hours: 4.0
  - kind: mypkg.transforms.MyWeighting   # your own pydantic class with an apply() method
    some_field: 3

A transform is any object with an apply(particles, context) method. The transforms guide covers the column-weighting science and how to write your own.

Loading results

import pandas as pd

for sim in model.simulations.sel(variant="hrrr"):
    traj = sim.trajectories
    foot = sim.footprint

# Load the footprints of a selection of simulations
footprints = model.simulations.sel(
    variant="hrrr", time=slice("2023-01-01", "2023-01-31")
).footprint.load()

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

for sim_id, footprint in footprints.items():       # keyed by (receptor, variant)
    hourly = footprint.aggregate(target=coords, time_bins=time_bins)

If a simulation's particles never reach the grid, PYSTILT writes a small .empty file instead of a NetCDF. The simulation counts as finished, and load() leaves it out.

STILT-R parity

PYSTILT footprints match the STILT-R footprints to a relative tolerance of 1e-7 in every grid cell. The test suite checks this against a pinned STILT-R commit. The NetCDF files are laid out differently from STILT-R's, so read them as standard CF-1.8 NetCDF. The STILT-R parity section of the development docs has the details.

Roadmap

PYSTILT borrows from two sister projects. X-STILT is the source of its column and satellite science, and stiltctl of its queue-based execution. The tables below show what PYSTILT has taken over so far. The full roadmap has more detail.

Execution and orchestration (from stiltctl)

Column and satellite science (from X-STILT)

PYSTILT does not try to match every X-STILT feature. It takes over X-STILT's approach to observations and column weighting, and leaves the rest.

Use of AI coding agents

This project is developed with the help of AI coding agents, directed and reviewed by the maintainer, who owns the design and the science.

Documentation

The full documentation is at https://jmineau.github.io/PYSTILT/.

Contributing

Contributions are welcome. See CONTRIBUTING.md for how to set up and submit changes.

License

PYSTILT is released under the MIT License. See the LICENSE file.

Author

James Mineau - jmineau