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, downloading meteorology from NOAA, and reading STILT-R output, install everything:

pip install "pystilt[complete]"

Quickstart

Describe a measurement, follow particles back from it through the meteorology, and calculate its footprint:

import stilt

receptor = stilt.PointReceptor(
    time="2023-07-15 18:00",   # UTC
    longitude=-111.848,
    latitude=40.766,
    altitude=10,               # meters above ground
)
met = {"directory": "/data/hrrr", "file_format": "%Y%m%d_%H", "file_tres": "6h"}
grid = stilt.Grid(xmin=-113.0, xmax=-110.5, ymin=40.0, ymax=42.0, xres=0.01, yres=0.01)

particles = stilt.run_trajectories(receptor, met, n_hours=-24, numpar=100)
foot = stilt.calc_footprint(particles, receptor, grid)

particles is a pandas DataFrame, one row per particle per step, and foot an xarray DataArray, one map per hour.

For many receptors, a project runs each one and keeps the results:

project = stilt.Project.init(
    "./my_project",
    receptors=[receptor],
    mets={"hrrr": met},
    variants={"hrrr": {}},
    n_hours=-24,
    numpar=100,
    grid=grid,
)
project.run()   # returns when the run is done

sim = project.simulation(receptor.id, "hrrr")
sim.particles, sim.footprint

Everything is saved in ./my_project. Later, open it again with stilt.Project("./my_project"). A second run() only runs what is missing. To change a setting, edit config.yaml.

Command line

stilt init ./my_project          # write a starter config.yaml and receptors.csv
stilt check ./my_project         # list what would stop a run, before any run
stilt run ./my_project           # run every receptor that is not done yet
stilt run ./my_project --cpus 8  # the same, eight receptors at a time
stilt status ./my_project        # what has finished, and the settings folders

stilt run exits with 0 when every simulation it ran is complete, 1 when some failed, and 3 when it was stopped, or a connection error stopped some, before they finished. An error it did not expect exits with 3 when it is a connection or I/O error (a bucket that cannot be reached), which a retry may get past, and with 4 for any other, which a retry would repeat. stilt check exits with 0 when it finds nothing and 1 when it lists anything. Every command exits with 2 when its command line is wrong; stilt run, stilt submit and stilt status also when a met directory does not exist or holds no met files.

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. The variants section declares them all; hrrr: {} runs the defaults.

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                         # realizations 0 to 3 of hrrr-err
    grid: null                              # particles only, no footprint
  hrrr-ak:                                  # another footprint from the hrrr particles
    transforms: [{kind: averaging_kernel, table: kernels}]

Each receptor under each variant is one simulation. Results go to the output directory (./output by default), in a folder named for the settings they were made with. Variants whose transport settings match share one HYSPLIT run and differ only in the footprint, as hrrr-ak does here. Adding a variant runs only what is new, and changing a setting writes to a new folder beside the old one.

On a cluster

Add an execution section to config.yaml and run the same command. PYSTILT splits the unfinished receptors among the tasks of one Slurm job array and submits it:

execution:
  backend: slurm
  n_workers: 200          # array tasks
  account: my-account
  partition: my-partition
  time: "02:00:00"
  mem: 4G

stilt submit ./my_project       # submits and returns; stilt run waits
stilt status ./my_project

Each task runs stilt run --task I/N, so a task is a command line you can read and rerun. A task that is preempted or nears its time limit requeues itself and skips the receptors it finished. Whether a simulation is finished is always decided by its files in the output directory. The same --task runs a project as a Kubernetes indexed Job.

Column and satellite soundings

Read the soundings into a table with one row per sounding, then make a receptor for each. Each sounding has its own averaging kernel, which goes into the project as a table:

import stilt
from stilt.observations import read_tropomi_ch4, receptors_from_soundings

df = read_tropomi_ch4(path)                      # or read_oco2, read_tccon, or your own reader
receptors, kernels = receptors_from_soundings(df[df.good], "slant", top=3000)

project = stilt.Project("./my_project")         # a project with a config.yaml
project.add_receptors(receptors)
project.add_table("kernels", kernels)
project.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
    coordinate: pres
  - kind: pressure_weighting

"column" makes vertical columns in place of slant lines of sight. The satellite tutorial goes from an orbit to modeled and observed enhancements.

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.BoundaryLayerOnly   # your own pydantic class with an apply() method
    max_height: 1200

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

Loading results

import pandas as pd

import stilt

sims = project.simulations                       # one row per receptor and variant
january = sims[(sims.variant == "hrrr") & sims.time.between("2023-01-01", "2023-01-31")]

project.status(january)                          # which results exist
footprints = project.footprints(january)         # one dataset: receptor, hour, lat, lon
mean = footprints.foot.sum("hour").mean("receptor")

windows = stilt.Mesh.from_windows([(-111.9, 40.7), (-111.8, 40.8)], size=0.05)
time_bins = pd.interval_range(
    start=pd.Timestamp("2022-12-31 00:00"),  # backward footprints reach a day back
    end=pd.Timestamp("2023-02-01 00:00"),
    freq="1h",
    closed="left",
)
H = project.jacobian(january, windows, time_bins)  # sparse: receptors by (hour, window)

stilt.sampling.sample_met adds met columns (mlht, dens, pres, zsfc) to a particle table, sampled from the met files at each particle's position and time:

from stilt.particles import particles_metadata
from stilt.sampling import sample_met

met_files = particles_metadata(sim.particles_path).met_files
sampled = sample_met(sim.particles.drop(columns="dens"), met_files, columns=["dens"])

If a simulation's particles never reach the grid, PYSTILT writes a footprint file with no cells and the reason. The simulation counts as finished, and project.footprints() lists it in attrs["empty"].

Runs made with STILT-R or X-STILT are read from their out/by-id folder, one simulation at a time, as a receptor, a particle table, and a footprint (pip install "pystilt[stilt-r]"):

from stilt.migration import read_stilt_r

for receptor, particles, foot in read_stilt_r("my_project/out/by-id"):
    print(receptor.id, len(particles), foot is not None)

stilt.migration.import_stilt_r writes them into a project's output directory instead, and prints the stilt-r variant that reads them, so project.status() and project.jacobian() work on runs STILT-R made.

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 page of the development docs has the details.

Roadmap

What is implemented and what is planned is on the roadmap page of the documentation.

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.

PYSTILT's wheels also carry a compiled hycs_std, the HYSPLIT build that STILT-R distributes. It is NOAA Air Resources Laboratory (ARL) software, not covered by the MIT License, and is subject to the HYSPLIT terms of use, which ask that publications using it cite NOAA ARL (Stein et al., 2015, doi:10.1175/BAMS-D-14-00110.1) and do not present results as official NOAA material.

Author

James Mineau - jmineau