Load And Plot Results#

Each simulation is one receptor run under one variant (Configuration). It writes up to two files:

  • the trajectories, every particle’s path, as a Parquet file

  • the footprint, as a NetCDF file, when the variant has a grid

A variant declared with from: makes its footprint from another variant’s particles. It has no trajectory file of its own, and sim.trajectories returns the other variant’s.

This page shows how to plot these outputs, load them for analysis, and add footprints up over the areas you care about.

Quick look#

Open the project, pick a simulation, and plot it:

import stilt

model = stilt.Model(project="./my_project")
sim = next(iter(model.simulations))          # the first simulation

sim.footprint.plot.map()                     # footprint, summed over time
sim.trajectories.plot.map()                  # particle paths
sim.plot.map()                               # receptor, particles, and footprint together

To pick a particular simulation, index by receptor id and variant:

model.simulations.keys()                     # all (receptor, variant) ids
sim = model.simulations["202307151800_-111.848_40.766_10", "hrrr"]
sim = model.simulations["202307151800_-111.848_40.766_10/hrrr"]   # same thing

Plotting needs the visualization extra. If cartopy is installed, maps also show coastlines and state borders. There are a few other plots:

  • foot.plot.facet() draws one panel per hour.

  • receptor.plot.map() shows where the receptor is.

  • model.plot.availability() shows which receptor times and locations have results.

Footprints#

foot = sim.footprint                 # None if there is no footprint file
foot.data                            # an xarray.DataArray
foot.time_range                      # (start, end) of the footprint's hours
foot.receptor                        # the receptor it belongs to

foot.data has dimensions (time, lat, lon), or (time, y, x) on a projected grid. There is one map for each hour back from the receptor time, in units of ppm per (µmol m⁻² s⁻¹). With time_integrate: true in the footprint settings there is a single time step. To sum over time:

total = foot.integrate_over_time()

To open a footprint file without a model:

foot = stilt.Footprint.from_netcdf("path/to/..._foot.nc")

The file also records the receptor and the settings used to make it.

Many simulations at once#

model.simulations holds every receptor under every variant. Narrow it with sel, then load an output from the result:

sims = model.simulations.sel(
    variant="hrrr",
    time=slice("2023-07-01", "2023-07-31 23:00"),    # receptor times, both ends included
)
footprints = sims.footprint.load()                   # {simulation id: Footprint}
paths = sims.footprint.paths()                       # {simulation id: Path}
trajectories = sims.trajectories.load()

The results are dictionaries keyed by simulation id, so you always know which receptor a result belongs to:

for sid, foot in footprints.items():
    print(sid.receptor, float(foot.integrate_over_time().sum()))

sel takes these filters, each as one value or a list:

  • receptor, receptor ids

  • variant, variant names. The name of a realization group, such as hrrr-err, selects all of its realizations.

  • time, one receptor time or a slice of times

  • location, location ids

  • where, a function that takes a receptor and returns True to keep it

Each call narrows the one before. A receptor id or variant name that the project does not have raises KeyError. The other filters may select nothing.

Extra columns in receptors.csv are on each receptor as attrs. Use them with where to gather one satellite scene or one site:

scene = model.simulations.sel(variant="hrrr", where=lambda r: r.attrs["scene"] == "A")

To see what is left to do:

model.simulations.incomplete()   # a selection, like sel()
model.simulations.status()       # a DataFrame, one row per simulation

status() has a trajectory and a footprint column that say whether each output exists. They are blank where the variant does not make that output. The empty column marks footprints that are empty (Empty footprints), and the complete column says whether the simulation is done. model.status() returns the same table. It checks every simulation, so it is slow on a large project stored in the cloud. From the command line, stilt status prints the totals, per variant when there are several.

Trajectories#

traj = sim.trajectories        # None if there is no trajectory file
df = traj.data                 # pandas DataFrame, one row per particle per time step

The columns you are most likely to use:

Column

Meaning

long / lati

Particle longitude and latitude

zagl

Particle height above ground, m

time

Minutes from the receptor time (negative for a backward run)

datetime

Time of the row, UTC

foot

The particle’s influence from the surface at this step, in ppm per (µmol m⁻² s⁻¹)

indx

Particle number

xhgt

Release height, for column and multipoint receptors

mlht, sigw, tlgr, pres

Mixed-layer height, vertical velocity spread, Lagrangian time scale, and pressure

To open a trajectory file without a model:

traj = stilt.Trajectories.from_parquet("path/to/..._traj.parquet")

Empty footprints#

Sometimes a simulation runs fine but no particle ever reaches the footprint grid. Usually the grid is too small or is not upwind. PYSTILT then writes a small <receptor id>_foot.empty file instead of a NetCDF. The simulation counts as finished, so reruns skip it. sim.footprint is None, sim.empty_reason says why, and load() and paths() leave the simulation out because there is nothing to load. The empty column of model.simulations.status() lists them. If you see many, make your footprint grid bigger.

An empty footprint is not a footprint of zeros. It means the transport never connected the receptor to your grid, so treating it as “the model says zero” in a comparison or an inversion would be wrong. Drop those observations, or find out why the particles never arrived.

Adding footprints up over areas#

Footprints are always calculated on the regular grid in your config, as in STILT-R. To get the influence of other areas, such as counties, hexagons, or small windows around point sources, use aggregate(). It adds up the footprint cells in each area, and the hours in each time bin:

import pandas as pd
import stilt

bins = pd.interval_range(
    start=foot.time_range[0], end=foot.time_range[1], freq="1h"
)
state = stilt.Grid(xmin=-112.3, xmax=-111.6, ymin=40.4, ymax=41.0,
                   xres=0.02, yres=0.02)
by_cell = foot.aggregate(state, time_bins=bins)        # index == state.index

The result is a DataFrame with one row per area and one column per time bin, labelled by the start of the bin. A footprint cell that straddles two areas is split between them by area, so the total influence is kept. This is what you want before multiplying by emissions. Influence that falls outside every area is dropped.

The target can be:

  • a stilt.Grid. Rows follow grid.index, one (lon, lat) pair per cell.

  • a stilt.Mesh of polygons with ids: a shapefile (Mesh.from_file), H3 hexagons (Mesh.from_h3), or windows around points (Mesh.from_windows). Rows are the polygon ids.

  • a stilt.Zones, which merges the cells of a grid or mesh into larger groups by label.

For cells given some other way, such as an xarray grid or a list of cell centres, build the stilt.Grid they lie on and select the rows you need from the result.

sources = stilt.Mesh.from_windows(
    [(-111.97, 40.515), (-112.015, 40.779)], 0.01, ids=["landfill", "wwtp"]
)
by_source = foot.aggregate(sources, time_bins=bins)    # index == ["landfill", "wwtp"]

counties = stilt.Mesh.from_file("counties.shp", ids="NAME")
by_county = foot.aggregate(counties, time_bins=bins)

sectors = stilt.Zones.from_labels(state, labels)       # one label per cell of state
by_sector = foot.aggregate(sectors, time_bins=bins)

Polygons in another coordinate system are reprojected onto the footprint grid. The overlaps between the footprint grid and your areas are worked out once and reused, so adding up thousands of footprints is fast. Polygon overlaps use shapely. If the optional exactextract package is installed (pip install exactextract), PYSTILT uses it instead. It gives the same result and is about a hundred times faster on large grids.

The footprint grid must be fine enough to resolve your areas. aggregate warns when the smallest area spans fewer than two footprint cells. In that case, calculate the footprint again from its particles on a finer grid:

hexes = stilt.Mesh.from_h3(8, bounds=state)
grid = stilt.Grid.from_geometry(hexes, cells_per_target=4)
fine = sim.generate_footprint(sim.footprint_config.model_copy(update={"grid": grid}))
by_hex = fine.aggregate(hexes, time_bins=bins)

Grid.from_geometry picks a grid that covers the areas with at least four cells across the smallest one. generate_footprint applies the variant’s particle transforms, as the stored footprint did, and does not overwrite the stored file.