Spectral Filtering Accuracy

PyStormTracker supports multiple backends for spherical harmonic transforms (SHT) used in spectral filtering and kinematic derivative calculations. The choice of engine can be controlled via the sht_engine parameter in the Python API.

Methodology

Accuracy is evaluated using Root Mean Square Error (RMSE), Relative Error, and Pearson Correlation Coefficient for ERA5 Mean Sea Level Pressure (MSL) and Wind (U/V) data.

Two spatial resolutions are tested:

  • 2.5°x2.5°: A low-resolution grid (73x144) where T42 truncation (requiring ~85 latitudes) results in aliasing.

  • 0.25°x0.25°: A high-resolution grid (721x1440) where T42 truncation is alias-free.

SHT Engines

1. ducc0 (Default)

ducc0 is the standard backend for PyStormTracker. It is a high-performance C++ library that provides excellent multi-frame performance and bit-wise parity with modern SHT implementations.

  • Pros: Extremely fast, robust, and handles aliased coarse grids gracefully. No external C dependencies (self-contained).

  • Cons: Slightly less accurate than specialized legacy libraries like SHTns on specific high-resolution scalar fields.

2. JAX (Experimental)

The JAX backend provides a JAX-native implementation of the SHT algorithms. It is designed for researchers looking to leverage GPU acceleration or automatic differentiation in their preprocessing pipelines.

  • Pros: GPU/TPU support via JAX. Machine-precision parity with ducc0 for resolved harmonics.

  • Cons: Requires the jax extra (pip install pystormtracker[jax]). Subject to significant aliasing errors on coarse grids (e.g., 2.5°) where lmax is close to the grid’s Nyquist frequency.

3. SHTns (Legacy Reference)

SHTns was the primary backend in earlier versions. While no longer the default, its performance metrics serve as a high-precision benchmark for scalar filtering.

  • Pros: Exceptional accuracy on high-resolution alias-free grids for scalar fields.

  • Cons: Slower than ducc0 for multi-frame workloads. Requires complex external C compilation. Significant numerical discrepancy (~10⁻⁶ RMSE) in kinematic derivative calculations compared to Spherepack (NCL).

Accuracy & Performance Metrics

Spectral Filtering (MSL)

Parity: JAX vs ducc0

Resolution

Target lmax

Mean Relative Error

Max Relative Error

RMSE

0.25°x0.25°

42

2.77e-15

5.49e-09

4.03e-16

1.0°x1.0°

42

2.67e-14

1.15e-10

5.93e-16

2.5°x2.5°

42

5.81e-03

2.71e+00

7.24e-04

Note: The high relative error at 2.5° is primarily due to differences in integration quadrature and aliasing treatment between the JAX-native matrix-vector approach and ducc0’s optimized kernels.

Cross-Engine Accuracy (vs NCL Reference)

Resolution: 2.5°x2.5° (ERA5) Note: Higher error in this resolution is primarily due to aliasing on the coarse grid.

Engine

Truncation

RMSE (Pa)

Rel. Error

Correlation

Time (s)

SHTns

T5-42

0.45711949

6.40e-04

0.999999897272

0.0350

ducc0

T5-42

0.05266872

7.37e-05

0.999999998663

0.0030

SHTns

T0-42

0.45729944

4.53e-06

0.999999918845

0.0009

ducc0

T0-42

0.05369308

5.31e-07

0.999999998880

0.0019

Resolution: 0.25°x0.25° (ERA5) Note: This resolution satisfies the sampling theorem for T42, resulting in near-perfect parity.

Engine

Truncation

RMSE (Pa)

Rel. Error

Correlation

Time (s)

SHTns

T5-42

0.00003643

5.16e-08

1.000000000000

0.0665

ducc0

T5-42

0.01276583

1.81e-05

0.999999999928

0.0041

SHTns

T0-42

0.00486497

4.81e-08

0.999999999994

0.0141

ducc0

T0-42

0.02114745

2.09e-07

0.999999999831

0.0038

Kinematic Derivatives (Vorticity & Divergence)

Evaluated on 0.25°x0.25° ERA5 grid. Comparison against NCL uv2vrdvF reference.

Engine

Variable

RMSE (s⁻¹)

Correlation

Time (s)

SHTns

Vorticity

9.05e-06

0.9899

0.1017

ducc0

Vorticity

1.74e-14

1.0000

0.0576

jax

Vorticity

1.12e-10

1.0000

3.49

JAX timings are measured on CPU; significantly faster performance is expected on CUDA-enabled devices.

Impact of Polar Optimization (SHTns)

For high-resolution grids (0.25°), disabling the polar optimization threshold (polar_opt=0.0) ensures maximum precision. A comparison with the default setting (1e-10) shows:

  • Accuracy: RMSE difference is negligible (~10⁻¹⁵ Pa).

  • Performance: The default setting is approximately 1.25% faster.

Disabling the optimization is recommended for cases where strict bit-wise parity with double-precision ground truth is required.

Summary

For standard storm tracking applications, all engines are scientifically equivalent. ducc0 is the recommended default for its balance of speed, robustness, and ease of installation. It also provides the highest parity with legacy Spherepack references for kinematics. The JAX backend is an excellent choice for high-resolution datasets where alias-free transforms ensure parity while enabling GPU acceleration. SHTns (with polar_opt=0.0) remains a high-precision benchmark for scalar fields, though it was passed over as the primary engine due to its numerical discrepancies in derivative calculations compared to standard meteorological tools.