API Reference

This page lists the principal Python modules and public classes used by PyStormTracker. The command-line orchestration function is available as pystormtracker.track.run_tracker.

Tracking orchestration

pystormtracker.track.run_tracker(infile: str, varname: str, outfile: str, start_time: str | datetime64 | None = None, end_time: str | datetime64 | None = None, mode: Literal['min', 'max'] = 'min', map_proj: Literal['global', 'nh_stereo', 'sh_stereo', 'healpix'] = 'global', resolution: float = 100.0, extent: tuple[float, float, float, float] | None = None, backend: Literal['serial', 'mpi', 'dask'] | None = None, n_workers: int | None = None, max_chunk_size: int | None = None, threshold: float | None = None, engine: str | None = None, algorithm: Literal['simple', 'hodges'] = 'simple', output_format: str = 'imilast', min_points: int = 1, w1: float | None = None, w2: float | None = None, dmax: float | None = None, phimax: float | None = None, n_iterations: int | None = None, min_lifetime: int | None = None, max_missing: int | None = None, zones: ndarray | None = None, adapt_params: ndarray | None = None, filter: bool | None = None, lmin: int = 5, lmax: int = 42, taper_points: int = 0, overlap: int = 3, subgrid_refine: bool | None = None) Tracks[source]

Orchestrates the storm tracking process from the CLI.

Trackers

class pystormtracker.simple.tracker.SimpleTracker[source]

Bases: object

A tracker implementing the PyStormTracker simple parallel algorithm.

preprocess_standard_track(data: xr.DataArray, lmin: int = 5, lmax: int = 42, taper_points: int = 0, map_proj: Literal['global', 'nh_stereo', 'sh_stereo', 'healpix'] = 'global', resolution: float = 100.0, extent: MapExtent | None = None, filter_type: Literal['sht', 'dct', 'auto'] = 'auto') xr.DataArray[source]

Applies standard spectral preprocessing using SHT or DCT. Optionally regrids to a Polar Stereographic or HEALPix projection.

track(infile: str | Path | xr.DataArray | xr.Dataset, varname: str, start_time: str | np.datetime64 | None = None, end_time: str | np.datetime64 | None = None, mode: Literal['min', 'max'] = 'min', map_proj: Literal['global', 'nh_stereo', 'sh_stereo', 'healpix'] = 'global', resolution: float = 100.0, extent: MapExtent | None = None, backend: Literal['serial', 'mpi', 'dask'] = 'serial', n_workers: int | None = None, max_chunk_size: int | None = None, threshold: float | None = None, engine: str | None = None, overlap: int = 3, min_points: int = 1, filter: bool = False, lmin: int = 5, lmax: int = 42, taper_points: int = 0, subgrid_refine: bool = False, **kwargs: float | int | str | None) Tracks[source]
class pystormtracker.hodges.tracker.HodgesTracker(w1: float = 0.2, w2: float = 0.8, dmax: float = 6.5, phimax: float = 0.5, n_iterations: int = 3, min_lifetime: int = 3, max_missing: int = 0, zones: ndarray[tuple[Any, ...], dtype[float64]] | None = None, adapt_params: ndarray[tuple[Any, ...], dtype[float64]] | None = None, use_standard_constraints: bool = True)[source]

Bases: Tracker

A tracker implementing the Hodges (TRACK) algorithm with adaptive constraints.

preprocess_standard_track(data: xr.DataArray, lmin: int = 5, lmax: int = 42, taper_points: int = 0, map_proj: Literal['global', 'nh_stereo', 'sh_stereo', 'healpix'] = 'global', resolution: float = 100.0, extent: MapExtent | None = None, filter_type: Literal['sht', 'dct', 'auto'] = 'auto') xr.DataArray[source]

Applies standard TRACK preprocessing: Tapering -> SHT or DCT Filter. Optionally regrids to a Polar Stereographic or HEALPix projection.

track(infile: str | Path | xr.DataArray | xr.Dataset, varname: str, start_time: str | np.datetime64 | None = None, end_time: str | np.datetime64 | None = None, mode: Literal['min', 'max'] = 'min', map_proj: Literal['global', 'nh_stereo', 'sh_stereo', 'healpix'] = 'global', resolution: float = 100.0, extent: MapExtent | None = None, backend: Literal['serial', 'mpi', 'dask'] = 'serial', n_workers: int | None = None, max_chunk_size: int | None = None, threshold: float | None = None, engine: str | None = None, overlap: int = 3, min_points: int = 1, filter: bool = True, lmin: int = 5, lmax: int = 42, taper_points: int = 0, subgrid_refine: bool = True, **kwargs: float | int | str | None) Tracks[source]

Runs the Hodges tracking algorithm. Supports chunked detection if max_chunk_size is provided. Detections are gathered before a single linking pass so chunk boundaries do not affect the result.

Parameters:
  • infile – Path to the input data file.

  • varname – Variable name to track.

  • start_time – Time range for tracking.

  • end_time – Time range for tracking.

  • mode – Search for ‘min’ or ‘max’ extrema.

  • backend – Processing backend (serial, mpi, dask).

  • n_workers – Number of parallel workers.

  • max_chunk_size – Number of steps per time chunk.

  • threshold – Intensity threshold for detection.

  • engine – Data loading engine (netcdf4, h5netcdf, etc).

  • overlap – Retained for cross-tracker API compatibility. Hodges gathers detections before linking and does not require overlap.

  • min_points – Minimum grid points per object.

  • filter – If True, apply spectral filtering.

  • lmin – Spectral truncation range (default T5-42).

  • lmax – Spectral truncation range (default T5-42).

  • taper_points – Boundary tapering points.

class pystormtracker.healpix.tracker.HealpixTracker(w1: float = 0.2, w2: float = 0.8, dmax: float = 6.5, phimax: float = 0.5, n_iterations: int = 3, min_lifetime: int = 3, max_missing: int = 0, zones: ndarray[tuple[Any, ...], dtype[float64]] | None = None, adapt_params: ndarray[tuple[Any, ...], dtype[float64]] | None = None, use_standard_constraints: bool = True)[source]

Bases: Tracker

A tracker specifically designed for 1D HEALPix grids.

preprocess_standard_track(data: DataArray, lmin: int = 5, lmax: int = 42, taper_points: int = 0) DataArray[source]

Apply spectral filtering and convert regular grids to HEALPix.

track(infile: str | Path | xr.DataArray | xr.Dataset, varname: str, start_time: str | np.datetime64 | None = None, end_time: str | np.datetime64 | None = None, mode: Literal['min', 'max'] = 'min', map_proj: Literal['global', 'nh_stereo', 'sh_stereo', 'healpix'] = 'global', resolution: float = 100.0, extent: MapExtent | None = None, backend: Literal['serial', 'mpi', 'dask'] = 'serial', n_workers: int | None = None, max_chunk_size: int | None = None, threshold: float | None = None, engine: str | None = None, overlap: int = 3, min_points: int = 1, filter: bool = True, lmin: int = 5, lmax: int = 42, taper_points: int = 0, subgrid_refine: bool = True, **kwargs: float | int | str | None) Tracks[source]

Detectors

class pystormtracker.simple.detector.SimpleDetector(pathname: str | Path | None, varname: str, time_range: TimeRange | None = None, global_start_idx: int = 0, global_total_steps: int | None = None, engine: str | None = None)[source]

Bases: object

A meteorological feature detector that treats fields as 2D images. Uses xarray for robust coordinate handling and lazy-loading.

detect(size: int = 5, threshold: float | None = None, minmaxmode: Literal['min', 'max'] = 'min', subgrid_refine: bool = False) list[tuple[datetime64, ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]], dict[str, ndarray[tuple[Any, ...], dtype[float64]]]]][source]
classmethod from_xarray(data: DataArray) SimpleDetector[source]

Creates a detector from an existing xarray DataArray.

get_time() ndarray[tuple[Any, ...], dtype[datetime64]] | None[source]
get_var(frame: int | tuple[int, int] | None = None) ndarray[tuple[Any, ...], dtype[float64]] | None[source]
get_xarray(start_time: str | datetime64 | None = None, end_time: str | datetime64 | None = None) DataArray[source]

Returns the requested data range as an xarray DataArray.

property lat: ndarray[tuple[Any, ...], dtype[float64]]
property lon: ndarray[tuple[Any, ...], dtype[float64]]
split(num: int) list[SimpleDetector][source]
class pystormtracker.hodges.detector.HodgesDetector(pathname: str | Path | None, varname: str, time_range: TimeRange | None = None, global_start_idx: int = 0, global_total_steps: int | None = None, engine: str | None = None)[source]

Bases: object

Feature detector based on the Hodges (TRACK) logic. Identifies local extrema (min/max) within thresholded objects.

detect(size: int = 5, threshold: float | None = None, minmaxmode: Literal['min', 'max'] = 'min', min_points: int = 1, subgrid_refine: bool = True) list[tuple[datetime64, ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]], dict[str, ndarray[tuple[Any, ...], dtype[float64]]]]][source]

Runs the feature detection on the selected time steps.

Parameters:
  • size – Diameter of local search window for extrema.

  • threshold – Intensity threshold for objects.

  • minmaxmode – Whether to search for local minima or maxima.

  • min_points – Minimum number of grid points in an object to be processed.

  • subgrid_refine – Whether to refine centers with a local quadratic fit.

classmethod from_xarray(data: DataArray) HodgesDetector[source]

Creates a detector from an existing xarray DataArray.

get_time() ndarray[tuple[Any, ...], dtype[datetime64]][source]
get_var(frame_idx: int | None = None) ndarray[tuple[Any, ...], dtype[float64]][source]
get_xarray(start_time: str | datetime64 | None = None, end_time: str | datetime64 | None = None) DataArray[source]

Returns the requested data range as an xarray DataArray.

property lat: ndarray[tuple[Any, ...], dtype[float64]]
property lon: ndarray[tuple[Any, ...], dtype[float64]]
class pystormtracker.healpix.detector.HealpixDetector(pathname: str | Path, varname: str, time_range: TimeRange | None = None, global_start_idx: int = 0, global_total_steps: int | None = None, engine: str | None = None)[source]

Bases: object

A meteorological feature detector that treats fields as 1D HEALPix maps. Uses xarray for lazy-loading and ducc0 for HEALPix grid math.

detect(threshold: float | None = None, minmaxmode: Literal['min', 'max'] = 'min', min_points: int = 1, subgrid_refine: bool = True) list[tuple[datetime64, ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]], dict[str, ndarray[tuple[Any, ...], dtype[float64]]]]][source]
classmethod from_xarray(data: DataArray, time_range: TimeRange | None = None, global_start_idx: int = 0, global_total_steps: int | None = None) HealpixDetector[source]
get_time() ndarray[tuple[Any, ...], dtype[datetime64]] | None[source]
get_var(frame: int | tuple[int, int] | None = None) ndarray[tuple[Any, ...], dtype[float64]] | None[source]
get_xarray() DataArray[source]
split(n: int) list[HealpixDetector][source]

Splits the detector into n smaller detectors with disjoint time ranges.

Core models

Tracks

class pystormtracker.models.tracks.TimeRange(start: datetime64, end: datetime64, step: timedelta64 | None = None)[source]

Bases: object

Metadata for the time range covered by a set of tracks.

end: datetime64
start: datetime64
step: timedelta64 | None
class pystormtracker.models.tracks.Track(track_id: int, tracks: Tracks)[source]

Bases: object

Represents a single storm track. In the array-backed architecture, it acts as a view into the parent Tracks object.

abs_dist(other: Track | Center) float[source]
append(center: Center) None[source]
extend(other: Track) None[source]
property indices: ndarray[tuple[Any, ...], dtype[int64]]
class pystormtracker.models.tracks.Tracks(track_ids: ndarray[tuple[Any, ...], dtype[int64]] | None = None, times: ndarray[tuple[Any, ...], dtype[datetime64]] | None = None, lats: ndarray[tuple[Any, ...], dtype[float64]] | None = None, lons: ndarray[tuple[Any, ...], dtype[float64]] | None = None, vars_dict: dict[str, ndarray[tuple[Any, ...], dtype[float64]]] | None = None, track_type: str = 'unknown')[source]

Bases: object

add_track(centers: list[Center]) Track[source]

Helper to append a new track from a list of Centers.

append(obj: Track) None[source]
bulk_append(tids: ndarray[tuple[Any, ...], dtype[int64]], times: ndarray[tuple[Any, ...], dtype[datetime64]], lats: ndarray[tuple[Any, ...], dtype[float64]], lons: ndarray[tuple[Any, ...], dtype[float64]], vars_dict: dict[str, ndarray[tuple[Any, ...], dtype[float64]]]) None[source]

Appends multiple points to multiple tracks at once.

property head: list[Track]
sort() None[source]

Sorts tracks by their first point’s time, lat, then lon.

property tail: list[Track]
time_range: TimeRange | None
property unique_track_ids: list[int]
write(outfile: str | Path, format: str = 'imilast') None[source]

Exports tracks to a file in the specified format.

Parameters:
  • outfile (str | Path) – Output file path.

  • format (str) – Output format. Supports “imilast”, “hodges”, “tracks-json”.

Storm centers

class pystormtracker.models.center.Center(time: datetime64, lat: float, lon: float, vars: dict[str, float])[source]

Bases: object

Represents a detected storm center at a specific time and location.

abs_dist(center: Center) float[source]

Calculate clamped great-circle distance in kilometers.

lat: float
lat_dist(center: Center) float[source]

Calculates the latitudinal distance in km.

lon: float
lon_dist(center: Center) float[source]

Calculates the longitudinal distance in km, adjusted for latitude.

time: datetime64
vars: dict[str, float]

Data loader

class pystormtracker.io.data_loader.DataLoader(pathname: str | Path | DataArray | Dataset | None = None, engine: str | None = None)[source]

Handles optimized xarray loading for local and remote datasets. Supports NetCDF, GRIB, and Zarr formats, with thread-safe caching.

ensure_open() Dataset[source]

Ensures the xarray dataset is open and returns it.

get_coords() tuple[str, str, str][source]

Returns the mapped names for (time, lat, lon).

get_grid_metadata(varname: str | None = None) dict[str, ndarray][source]

Returns grid metadata (theta, nphi, phi0, ringstart) for SHT. Works for reduced Gaussian and HEALPix grids.

get_reduced_grid_pl(varname: str | None = None) ndarray | None[source]

Returns the ‘pl’ array (points per latitude) for a reduced grid.

is_global_longitude() bool[source]

Return whether a 1D longitude coordinate covers a periodic globe.

is_lat_reversed() bool[source]

Detects if the latitude coordinate is North-to-South (reversed). Returns True if lat[0] > lat[-1].

is_reduced_gaussian(varname: str | None = None) bool[source]

Detects if the dataset represents a reduced Gaussian grid.

Preprocessing

Kinematics (vorticity and divergence)

class pystormtracker.preprocessing.kinematics.Kinematics(R: float = 6371220.0, lmax: int | None = None, geometry: str = 'CC', lat_reverse: bool = False)[source]

Computes spatial derivatives and kinematic properties of the wind field.

compute(u: DataArray, v: DataArray, backend: Literal['serial', 'mpi', 'dask'] = 'serial', nthreads: int = 0) tuple[DataArray, DataArray][source]
compute(u: ndarray[tuple[Any, ...], dtype[float64]], v: ndarray[tuple[Any, ...], dtype[float64]], backend: Literal['serial', 'mpi', 'dask'] = 'serial', nthreads: int = 0) tuple[ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]]]

Computes vorticity and divergence from wind components.

Parameters:
  • u – Zonal wind component.

  • v – Meridional wind component.

  • backend – Parallelization backend (‘serial’, ‘mpi’, ‘dask’).

  • nthreads – Number of threads (for local computation).

Returns:

Divergence and relative vorticity.

Return type:

div, vort

class pystormtracker.preprocessing.kinematics.KinematicsKwargs[source]
pystormtracker.preprocessing.kinematics.apply_vort_div(u: DataArray, v: DataArray, R: float = 6371220.0, lmax: int | None = None, geometry: str = 'CC', nthreads: int = 0, backend: Literal['serial', 'mpi', 'dask'] = 'serial') tuple[DataArray, DataArray][source]

Xarray wrapper for computing relative vorticity and divergence.

Parameters:
  • u – Zonal wind DataArray.

  • v – Meridional wind DataArray.

  • R – Planetary radius in meters. Default is R_EARTH_METERS.

  • lmax – Maximum spherical harmonic degree.

  • geometry – Grid geometry (default ‘CC’).

  • nthreads – Number of threads.

  • backend – Parallelization backend. Options: ‘serial’, ‘mpi’, ‘dask’.

Returns:

Divergence and relative vorticity DataArrays.

Return type:

div, vort

pystormtracker.preprocessing.kinematics.compute_vort_div(u: ndarray[tuple[Any, ...], dtype[float64]], v: ndarray[tuple[Any, ...], dtype[float64]], R: float = 6371220.0, lmax: int | None = None, geometry: str = 'CC', nthreads: int = 0, lat_reverse: bool = False) tuple[ndarray[tuple[Any, ...], dtype[float64]], ndarray[tuple[Any, ...], dtype[float64]]][source]

Computes spatial divergence and relative vorticity from u and v wind components using ducc0.

Parameters:
  • u – Zonal wind (ntheta, nphi).

  • v – Meridional wind (ntheta, nphi).

  • R – Planetary radius in meters. Default is R_EARTH_METERS.

  • lmax – Maximum spherical harmonic degree. If None, derived from ntheta.

  • geometry – Grid geometry (for ducc0). Default ‘CC’.

  • nthreads – Number of threads (for ducc0).

  • lat_reverse – If True, assume latitude is North to South (reversed).

Returns:

Divergence (ntheta, nphi) vort: Relative vorticity (ntheta, nphi)

Return type:

div

Spectral filtering

class pystormtracker.preprocessing.spectral.DCTFilter(lmin: int = 5, lmax: int = 42, taper_val: float = 0.1)[source]

Spectral bandpass filter for regional grids using Discrete Cosine Transform.

class pystormtracker.preprocessing.spectral.FilterKwargs[source]
class pystormtracker.preprocessing.spectral.SHTFilter(lmin: int = 5, lmax: int = 42, lat_reverse: bool = False, taper_val: float = 0.1, geometry: Literal['CC', 'GL', 'DH', 'auto'] = 'auto')[source]

Spectral bandpass filter (truncation) for lat-lon grid data using ducc0.

filter(data: DataArray, backend: Literal['serial', 'mpi', 'dask'] = 'serial') DataArray[source]
filter(data: ndarray[tuple[Any, ...], dtype[float64]], backend: Literal['serial', 'mpi', 'dask'] = 'serial') ndarray[tuple[Any, ...], dtype[float64]]

Applies the filter to the input data.

Parameters:
  • data (xr.DataArray | np.ndarray) – Input data.

  • backend (str) – Parallelization backend. Options: ‘serial’, ‘mpi’, ‘dask’.

Returns:

The filtered data.

Return type:

xr.DataArray | np.ndarray

pystormtracker.preprocessing.spectral.apply_bandpass_mask_to_alm(alm: ndarray[tuple[Any, ...], dtype[complex128]], lmin: int, lmax: int, mmax: int | None = None, taper_val: float = 0.1) None[source]

Applies a tapered bandpass mask in-place to spherical harmonic coefficients. Uses the form from Hoskins and Sardeshmukh (1984), acting like a ∇⁴ smoother.

Parameters:
  • alm – Spherical harmonic coefficients.

  • lmin – Minimum wave number (hard cutoff).

  • lmax – Maximum wave number (tapered to taper_val).

  • mmax – Maximum m wave number.

  • taper_val – Value of the taper at lmax (default 0.1).

pystormtracker.preprocessing.spectral.apply_dct_filter(data: DataArray, lmin: int = 5, lmax: int = 42, backend: Literal['serial', 'mpi', 'dask'] = 'serial', taper_val: float = 0.1) DataArray[source]

Applies a DCT-based spectral bandpass filter to a regional DataArray. Uses ducc0.fft.dct for high-performance 2D DCT.

Parameters:
  • data – Input DataArray (regional).

  • lmin – Minimum total wave number equivalent to retain.

  • lmax – Maximum total wave number equivalent to retain.

  • backend – Parallelization backend.

  • taper_val – Value of the taper at lmax.

Returns:

The filtered regional data.

Return type:

xr.DataArray

pystormtracker.preprocessing.spectral.apply_sht_filter(data: DataArray, lmin: int = 5, lmax: int = 42, lat_reverse: bool = False, backend: Literal['serial', 'mpi', 'dask'] = 'serial', taper_val: float = 0.1, geometry: Literal['CC', 'GL', 'DH', 'auto'] = 'auto', out_geometry: str | None = None, out_ntheta: int | None = None, out_nphi: int | None = None) DataArray[source]

Applies a spectral bandpass filter to the input DataArray. Supports regular 2D grids and reduced Gaussian 1D grids.

Parameters:
  • data (xr.DataArray) – Input data.

  • lmin (int) – Minimum total wave number to retain. Defaults to 5.

  • lmax (int) – Maximum total wave number to retain. Defaults to 42.

  • lat_reverse (bool) – If True, assume latitude is South to North.

  • backend (str) – Parallelization backend. Options: ‘serial’, ‘mpi’, ‘dask’.

  • taper_val (float) – Value of the taper at lmax.

  • geometry (str) – Grid geometry (‘CC’, ‘GL’, ‘DH’, or ‘auto’).

  • out_geometry (str) – Output geometry (e.g., ‘CC’, ‘GL’, ‘HEALPix’).

  • out_ntheta (int) – Number of latitudes in output grid.

  • out_nphi (int) – Number of longitudes in output grid.

Returns:

The filtered (and possibly regridded) data.

Return type:

xr.DataArray

pystormtracker.preprocessing.spectral.is_global_grid(data: DataArray) bool[source]

Heuristic to determine if an xarray DataArray represents a global grid. Checks longitude range and presence of healpix attributes.

Regridding

class pystormtracker.preprocessing.regrid.SpectralRegridder(lmax: int | None = None, mmax: int | None = None)[source]

Spectral regridder for transforming data between Clenshaw-Curtis (CC), Gauss-Legendre (GL), and HEALPix grids using ducc0 SHT.

to_grid(data: DataArray, nlat: int, nlon: int, in_geometry: Literal['CC', 'GL'] = 'CC', out_geometry: Literal['CC', 'GL'] = 'CC', lat_reverse: bool = False, nthreads: int = 1, pl: ndarray[tuple[Any, ...], dtype[int32]] | None = None) DataArray[source]

Spectrally regrid to a regular 2D grid (CC or GL). Supports regular 2D and reduced Gaussian 1D inputs.

to_healpix(data: DataArray, nside: int, in_geometry: Literal['CC', 'GL'] = 'CC', lat_reverse: bool = False, nthreads: int = 1, pl: ndarray[tuple[Any, ...], dtype[int32]] | None = None) DataArray[source]

Spectrally regrid to a 1D HEALPix grid. Supports regular 2D and reduced Gaussian 1D inputs.

to_polar_stereo(data: xr.DataArray, hemisphere: Literal['nh', 'sh'] = 'nh', extent: MapExtent = (-13000.0, 13000.0, -13000.0, 13000.0), resolution: float = 100.0, lon_0: float = 0.0, filter_lmin: int | None = None, lmax: int | None = None, in_geometry: Literal['CC', 'GL'] = 'CC', lat_reverse: bool = False, nthreads: int = 1) xr.DataArray[source]

Spectrally regrid to a Polar Stereographic grid.

Parameters:
  • extent – Bounding box from pole in km (xmin, xmax, ymin, ymax).

  • resolution – Grid spacing in km.

  • lmax – Maximum total wave number. Overrides the constructor value.

Sub-grid refinement

pystormtracker.preprocessing.refinement.subgrid_refine(frame: ndarray[tuple[Any, ...], dtype[float64]], r: int, c: int, y: ndarray[tuple[Any, ...], dtype[float64]], x: ndarray[tuple[Any, ...], dtype[float64]], periodic_x: bool = True) tuple[float, float, float][source]

Refine an extremum with a quadratic fit over its 3x3 neighborhood.

Spatial tapering

class pystormtracker.preprocessing.taper.TaperFilter(n_points: int = 10)[source]

A filter that applies a cosine taper to the edges of the domain.

This is often used in TRACK before spherical harmonic filtering to minimize the Gibbs phenomenon (ringing artifacts) at the boundaries. Since spherical harmonic transforms effectively assume periodicity or symmetry, abrupt changes at the boundaries (e.g. between the North Pole and the first latitude row) can introduce high-frequency noise.

filter(data: DataArray | ndarray[tuple[Any, ...], dtype[float64]]) DataArray | ndarray[tuple[Any, ...], dtype[float64]][source]

Applies a cosine taper to the edges of the input data.

Parameters:

data (xr.DataArray | np.ndarray) – Input data.

Returns:

The tapered data.

Return type:

xr.DataArray | np.ndarray

Post-processing and analysis

Secondary-variable sampling

pystormtracker.sample.sample_tracks(tracks: Tracks, ds: Dataset, varname: str, method: Literal['nearest', 'bilinear', 'mean', 'max', 'min'] = 'nearest', radius_km: float = 0.0, output_varname: str | None = None) Tracks[source]

Samples a variable from a NetCDF dataset along storm tracks.

Parameters:
  • tracks – The Tracks object to update.

  • ds – The xarray Dataset containing the variable to sample.

  • varname – The name of the variable in the dataset.

  • method – The sampling method (‘nearest’, ‘bilinear’, ‘mean’, ‘max’, ‘min’).

  • radius_km – The radius in km for spatial operations (mean, max, min).

  • output_varname – The name to store in the track’s ‘vars’ dictionary. Defaults to varname.

Returns:

The updated Tracks object.

Track matching

pystormtracker.metrics.compare.match_tracks(tracks_ref: Tracks, tracks_comp: Tracks, max_dist_km: float = 440.0, min_overlap_fraction: float = 0.1) dict[int, int][source]

Matches tracks from a comparison set to a reference set based on spatial proximity and temporal overlap, mirroring the algorithm in TRACK’s trdist_eps.c.

Parameters:
  • tracks_ref – The reference Tracks object.

  • tracks_comp – The comparison Tracks object.

  • max_dist_km – Maximum mean geodetic distance (in km) allowed between tracks. Default is 440.0 km (~4 degrees).

  • min_overlap_fraction – Minimum overlap ratio required. Ratio = (2.0 * overlap_len) / (len(A) + len(B)). Default is 0.1 (10%).

Returns:

A dictionary mapping comparison track IDs to reference track IDs.

Spherical weighting kernels

class pystormtracker.metrics.weighting.WeightType(*values)[source]

Supported weight types for spherical kernel estimation.

pystormtracker.metrics.weighting.calculate_spherical_weight(dist_km: float, radius_km: float, weight_type: int, kappa: float) float[source]

Computes weight based on kernel type.

Parameters:
  • dist_km – Geodesic distance in kilometers.

  • radius_km – Radius of influence in kilometers.

  • weight_type – Integer ID from WeightType enum.

  • kappa – Smoothing parameter for Fisher kernel.

Returns:

Computed weight.

Return type:

float

Eulerian track metrics

pystormtracker.metrics.eulerian.compute_eke(u: DataArray, v: DataArray, freq: str = 'MS', outfile: str | Path | None = None) DataArray[source]

Computes Eddy Kinetic Energy (EKE) using a 24-hour difference filter. EKE = 1/2 * [Var(u) + Var(v)]

Parameters:
  • u – Zonal wind DataArray.

  • v – Meridional wind DataArray.

  • freq – Resampling frequency (default “MS” for monthly start).

  • outfile – Optional path to save the result directly to disk.

Returns:

The computed EKE.

Return type:

xr.DataArray

pystormtracker.metrics.eulerian.compute_high_wind_index(ds: Dataset, u_var: str, v_var: str, freq: str = 'MS', percentile: float = 0.95, outfile: str | Path | None = None) DataArray[source]

Computes the high-wind index (e.g., 95th percentile of wind speed). Used as a weather impact metric (Yau and Chang 2020).

Parameters:
  • ds – Dataset containing u and v wind components.

  • u_var – Variable name for zonal wind.

  • v_var – Variable name for meridional wind.

  • freq – Resampling frequency (default “MS” for monthly start).

  • percentile – Quantile to calculate (default 0.95).

  • outfile – Optional path to save the result directly to disk.

Returns:

The computed high-wind index.

Return type:

xr.DataArray

pystormtracker.metrics.eulerian.compute_variance_metric(da: DataArray, freq: str = 'MS', outfile: str | Path | None = None) DataArray[source]

Computes the Eulerian variance metric using a 24-hour difference filter. e.g., Var(SLP) = [SLP(t + 24h) - SLP(t)]^2 averaged over the period.

Parameters:
  • da – Input DataArray (e.g., SLP or z500).

  • freq – Resampling frequency (default “MS” for monthly start).

  • outfile – Optional path to save the result directly to disk.

Returns:

The computed variance metric.

Return type:

xr.DataArray

Lagrangian track metrics

pystormtracker.metrics.tracks.compute_track_metrics(tracks: Tracks, grid_lat: ndarray[tuple[Any, ...], dtype[float64]], grid_lon: ndarray[tuple[Any, ...], dtype[float64]], radius_km: float = 500.0, kernel: Literal['constant', 'fisher', 'cressman', 'linear', 'quadratic'] = 'constant', kappa: float = 20.0, varname: str | None = None, is_min: bool = False, monthly: bool = True) Dataset[source]

Computes storm track metrics on a 2D spatial grid using weighted estimators. Supports 5 Lagrangian metrics (Yau and Chang 2020, Hodges 1999, Simmonds 2026): - cyclone_amplitude - cyclone_frequency (weighted) - track_frequency (weighted) - aca (Accumulated Cyclone Activity) - ata (Accumulated Track Activity)

Parameters:
  • tracks – Tracks object containing the storm tracks.

  • grid_lat – 1D array of latitude coordinates.

  • grid_lon – 1D array of longitude coordinates.

  • radius_km – Radius of influence in km. Default 500km (Yau & Chang).

  • kernel – Kernel type: ‘constant’, ‘fisher’, ‘cressman’, ‘linear’, ‘quadratic’.

  • kappa – Smoothing parameter for Fisher kernel (default 20.0).

  • varname – Variable in tracks.vars to use as amplitude.

  • is_min – If True, tracks are defined by minima (e.g., SLP).

  • monthly – If True (default), metrics are aggregated into monthly values.

Returns:

Dataset containing the computed metrics.

Return type:

xr.Dataset

CORMAX and CCA/PCA cross-validation

class pystormtracker.metrics.cross_validation.CCAModel(*args, **kwargs)[source]

Operations used from a fitted xeofs CCA model.

pystormtracker.metrics.cross_validation.compute_cormax(impact_da: DataArray, metric_da: DataArray, search_lon: float = 60.0, search_lat: float = 20.0) DataArray[source]

Computes the CORMAX score: maximum one-point correlation within a local region. For each grid point in impact_da, find the max correlation with metric_da within a search_lon x search_lat window (Yau and Chang 2020).

Parameters:
  • impact_da – Weather impact anomalies (time, lat, lon).

  • metric_da – Storm track metric anomalies (time, lat, lon).

  • search_lon – Longitude window width in degrees (default 60).

  • search_lat – Latitude window width in degrees (default 20).

Returns:

CORMAX scores (lat, lon).

Return type:

xr.DataArray

pystormtracker.metrics.cross_validation.find_best_cca_truncation(X: DataArray, Y: DataArray, max_modes: int = 15, leave_n_out: int = 3, pca: bool = True) Dataset[source]

Evaluates CCA truncation numbers (M) to find the best hyperparameter.

Parameters:
  • X – Predictor field, dimensions (time, lat, lon).

  • Y – Predictand field, dimensions (time, lat, lon).

  • max_modes – Maximum number of PCA/CCA modes to test.

  • leave_n_out – Samples to hold out (default 3 for seasonal).

  • pca – If True, performs PCA pre-filtering before CCA.

Returns:

Dataset containing ‘acc’ and ‘fve’ scores.

Return type:

xr.Dataset

pystormtracker.metrics.cross_validation.train_cca_model(X: DataArray, Y: DataArray, n_modes: int) CCAModel[source]

Trains a CCA model with the given number of modes on the full dataset.

Parameters:
  • X – Predictor field, dimensions (time, lat, lon).

  • Y – Predictand field, dimensions (time, lat, lon).

  • n_modes – Number of PCA/CCA modes to use (usually the best M).

  • pca – If True, performs PCA pre-filtering before CCA.

Returns:

The trained xeofs CCA model.