STILT-R Parity#

PYSTILT is tested against a fixed commit of STILT-R. This page says what “matches STILT-R” covers and how to move to a newer STILT-R commit.

Pinned upstream commit#

Field

Value

Upstream

uataq/stilt

Pinned SHA

0e290a68730e155bc2156e28af492a463228c88a

Pin location

.github/workflows/tests.yml (env STILT_R_SHA)

Last verified

2026-06-25

When the fidelity tests pass, PYSTILT matches STILT-R at this commit. It says nothing about STILT-R’s current main branch.

What matches#

PYSTILT’s footprints match STILT-R’s for the 20 fidelity scenarios in tests/fixtures/r_stilt_reference.py, within these tolerances:

Quantity

Tolerance

Footprint, per cell

rtol=1e-7, atol=1e-8

Footprint total

rtol=1e-6

Footprint peak

rtol=1e-7

Number of nonzero cells

identical

Intermediate tables (interpolation, rtime, gridding)

rtol=1e-12

Trajectory positions and foot after the HNF correction

rtol=1e-7

Grid coordinates, longitude/latitude

atol=1e-12 degrees

Grid coordinates, projected

atol=1e-3 m

The scenarios cover:

  • point, column, and three-location multipoint receptors

  • 6 h and 24 h backward runs, and a 6 h forward run

  • the hyper-near-field (HNF) plume correction on (12 scenarios) and off

  • smooth_factor of 0, 0.5, 1, and 2

  • 0.01° and 0.05° longitude/latitude grids and a UTM grid, all with square cells

  • a 2° box, a receptor in the domain’s corner, a receptor at the edge of the met grid, and a small 0.1° domain that particles leave early

  • receptor heights above ground, above sea level (kmsl=1), and at 0.5 m

  • winter (stable) and summer (convective) HRRR

  • hourly footprints and time_integrate=True

  • a wind-error run (siguverr=1 m/s)

The synthetic tests in tests/r_stilt/test_footprint_synth.py feed hand-made particle tables to both implementations to test single code paths: one Gaussian particle, cells on the grid boundary, the dateline, a global grid, and the latitude scaling of the kernel width. They need R but no meteorology or HYSPLIT, so they run anywhere R is installed:

STILT_R_DIR=$PWD/stilt-r-src uv run pytest tests/r_stilt -m r_only

What does not match#

  • The NetCDF layout. STILT-R writes dimensions (x, y, time) with a fill value of -1. PYSTILT writes CF-1.8 NetCDF with dimensions (time, lat, lon). Read PYSTILT files as ordinary CF NetCDF. Tools that expect STILT-R’s layout will not read them.

  • Smoothing on cells that are not square. When xres and yres differ, a smoothing kernel wider than one cell has a different number of cells along longitude and along latitude. STILT-R’s permute.f90 reads the kernel as if the two numbers were equal. With xres larger than yres, its footprint total comes out 5 to 13 % high on 24 h runs. With yres larger than xres, it reads memory past the end of the kernel. The footprint then loses mass (15 % in one case), holds values near 1e300, or R crashes (uataq/stilt#116; the fix is uataq/stilt#148). PYSTILT adds each kernel along the axis each cell count belongs to, so smoothing keeps the footprint total. Its footprint equals STILT-R’s with permute.f90 reading the kernel that way, to within 2e-15 of the peak. On square cells the two agree.

  • Smoothing on a projected grid. The kernel width, 0.06 * (var(x) + var(y))**0.25 * sqrt(|rtime| / 1440), is made for coordinates in degrees. STILT-R applies it to the projected coordinates, so on a grid in meters the kernel is a few meters wide and the footprint is hardly smoothed. On a 24 h run on 1 km UTM cells, STILT-R’s footprint has a third of the nonzero cells of the same run on a 0.01° grid (9,771 and 28,706), and PYSTILT’s has 26,385. PYSTILT works the width out in degrees from the particles’ longitudes and latitudes, as on a longitude/latitude grid, then multiplies it by the length of one degree of longitude at the particles’ mean position, measured in the grid’s projection and units. Each step’s longitudes are first taken within 180° of that step’s circular mean, so a cloud is measured as one wherever it is, on both sides of 180° included. The kernel is then as wide from east to west, on the ground, as on a longitude/latitude grid. On the projected grid it is round. On a longitude/latitude grid it is 1 / cos(lat) times wider from north to south, because it has one width in degrees on both axes. The fidelity suite’s UTM scenario compares unsmoothed footprints, and test_projected_smoothing_intentionally_differs_from_r guards the difference. Longitude/latitude grids are not affected.

  • The early-time interpolation on a projected grid. In the first 100 minutes, when the particles’ median step is longer than a cell along either axis, both add positions between the output steps. STILT-R measures the step in degrees and compares it with xres and yres in the grid’s units, so on a grid in meters it never interpolates. On a 24 h run on 100 m UTM cells, unsmoothed, STILT-R’s rule gives a footprint with 12,341 nonzero cells, against 15,139 for the same run on a 0.001° grid, which does interpolate. PYSTILT measures the step on the projected grid, in its units and along its axes, and its footprint has 14,563. It then interpolates as STILT-R does on longitude/latitude cells of the same size. The positions added are interpolated in degrees on any grid. test_utm_early_interpolation_intentionally_differs_from_r guards the difference. Longitude/latitude grids are not affected.

  • Smoothing a cloud near 180° on a global grid. On a grid that goes all the way around (xmin=-180, xmax=180), STILT-R leaves the particles’ longitudes as HYSPLIT writes them, from -180 to 180. A cloud on both sides of 180° then has a variance of thousands of square degrees and gets a kernel many times too wide, and what the kernel spreads past ±180° is lost. For a cloud of 300 particles (standard deviation 0.3°, a day after release) on a 0.25° grid, a cloud at 175° covers 36 cells and keeps its whole footprint in both programs; at 179.9° STILT-R covers 1,242 cells and loses 42 % of the footprint, and PYSTILT covers 34 cells and keeps all of it. PYSTILT takes the longitudes of each time step within 180° of that step’s mean longitude (the circular mean, so the seam is on the far side of the globe from the cloud) before it measures the spread. The cloud is measured as one wherever it is, and a real spread of 10° stays 10°. PYSTILT also adds the smoothing that passes one edge of such a grid to the cells at the other edge, so the footprint of a cloud near 180° is the footprint of the same cloud anywhere else, shifted. The longitudes are unwrapped this way on every longitude/latitude grid, which also fixes a cloud around 0° on a wide grid that crosses 180°, where STILT-R’s shift to 0° to 360° cuts the cloud in two. A grid whose cells do not fit 360° (0.7° cells leave 0.2°) has an edge, and smoothing is cut there as before. A cloud that is not cut by the seam is measured as before, bit for bit. test_global_grid_across_180_intentionally_differs_from_r guards the difference.

  • Forward runs with the HNF correction. STILT-R accumulates the plume from the far end of a forward trajectory back toward the release point, so its plume is widest at release. PYSTILT accumulates outward from the release point. The two foot values differ on purpose, and test_forward_hnf_foot_intentionally_differs_from_r guards the difference.

  • The HNF correction with ``veght`` above 1. HYSPLIT reads such a veght as meters above ground, and so does PYSTILT’s correction. STILT-R multiplies it by the mixed-layer height, which removes most of the footprint after the first hour (uataq/stilt#142). No fidelity scenario sets veght above 1.

  • The HNF correction of a receptor above sea level. PYSTILT grows the plume from the release height above ground, so it needs zsfc in varsiwant for an MSL receptor. STILT-R grows it from the receptor height as given, which with kmsl=1 is above sea level, so over high terrain its correction does nothing. The msl_altitude scenario compares the uncorrected foot_no_hnf_dilution, and test_msl_hnf_foot_intentionally_differs_from_r guards the difference.

  • Untested inputs. Runs longer than 24 h backward, latitudes above 80°, and grids finer than 0.001° have not been compared.

  • Other STILT-R commits. Nothing detects upstream changes to calc_footprint.r, permute.f90, or calc_trajectory.r. Bump the pinned SHA by hand (below).

Large-scale comparison#

Outside CI, PYSTILT was compared with STILT-R on 200 receptors at the WBB tower in the Salt Lake Valley, 2016 to 2024 (35 m above ground, 1000 particles, 24 h backward, HRRR, krand=2 and seed=42, the same v5.1.0 hycs_std on both sides). For each receptor, the PYSTILT trajectory was compared with a separate STILT-R calc_trajectory run, and the PYSTILT footprint with STILT-R’s calc_footprint of the same particles, at 0.01°, 0.05°, and 0.1°. All 200 receptors matched.

Quantity

Agreement (worst of 200 receptors)

Trajectory, per particle

absolute 1.0e-17, relative 5.6e-16 (float64 round-off)

Footprint, per-cell relative difference

median ~2e-8, max ~6e-8

Footprint nonzero cells

identical (for example 250,644 cells at 0.01°)

Footprint total

relative ~2e-8

The trajectories are identical, so the footprint differences come only from the order in which floating-point sums are added. The relative difference is the same for small cells (under 0.01 % of the peak) as for large ones. Footprints are stored as float32, whose precision is about 1.2e-7, so the two outputs agree to the precision the file can hold. The mass-weighted relative difference, which is what an inversion sums over, is also about 2e-8.

Reading the same met files on both sides over the full date range needs the find_met_files fix in the pinned STILT-R commit. HRRR archives can hold the same file both at the top level (as a symlink) and in a YYYY/MM/ subdirectory, and older STILT-R listed it twice.

Moving to a newer STILT-R commit#

  1. Change STILT_R_SHA in .github/workflows/tests.yml.

  2. Run the fidelity tests locally:

    git clone https://github.com/uataq/stilt stilt-r-src
git -C stilt-r-src checkout NEW_SHA   # the commit to pin
STILT_R_DIR=$PWD/stilt-r-src STILT_TEST_MET_DIR=tests/met_cache \
  uv run pytest tests/r_stilt/ -v -m fidelity

  3. If every scenario passes at the current tolerances, commit the change and update “Last verified” above.

  4. If a scenario fails, diff calc_footprint.r, permute.f90, and calc_trajectory.r between the two commits to find the upstream change. Then do one of the following:

    • make the same change in PYSTILT,

    • loosen that scenario’s tolerance, with a comment saying why, or

    • stay on the old commit.

Files to diff on every bump#

These STILT-R files determine the numbers PYSTILT is compared against:

  • r/src/calc_footprint.r: footprint kernel and gridding

  • r/src/calc_trajectory.r: the HYSPLIT wrapper, random seed, and krand

  • r/src/permute.f90: the Fortran loop that adds each particle’s kernel to the grid

The R helper scripts in tests/fixtures/r_helpers/ call these files for the synthetic tests. They may need updating if a function signature changes upstream.