PyStormTracker's HodgesTracker implements the feature-identification and
trajectory-linking method developed by Hodges (1994, 1995, 1999), with selected
implementation details reconciled against
TRACK 1.5.4
(tag TRACK-1.5.4, commit 6ded301a5f5183d73e5b49c16019024b9a53eff7).
The papers are the scientific authority for the method. TRACK 1.5.4 is the implementation reference for source-specific behavior that is not fully specified in the papers. This page provides function-level TRACK source references for the supported workflow: source-dependent statements link directly to the tagged TRACK 1.5.4 implementation and relevant line ranges.
The implementation combines several layers that should not be conflated:
| Layer | Authority | PyStormTracker relationship |
|---|---|---|
| Feature identification and tracking method | Hodges (1994, 1995, 1999) | Scientific lineage |
| Source-specific tracking semantics | TRACK 1.5.4 | Implementation/parity reference |
| Rectangular and spherical B-spline construction | Dierckx/FITPACK through SciPy | Established numerical implementation |
| Spherical harmonics and HEALPix numerics | spharmgrid, ducc0 |
spharmgrid owns supported rectangular GL/CC operations; uses direct DUCC on special SHT paths |
| Spherical quadratic and intrinsic spherical B-spline optimization | PyStormTracker | Explicit extensions |
| Global one-to-one track assignment and exact timestamp-sequence identity | PyStormTracker | Comparison/validation extensions |
TRACK's spatial spectral-filter workflow is implemented by
spectral_filter().
The interactive wrapper is
spec_filt(),
while the Hoskins coefficient taper itself is implemented by
hoskins_filt().
PyStormTracker uses public spharmgrid filtering for supported rectangular
Gauss--Legendre and Clenshaw--Curtis fields. NumPy input reaches that operation
through a coordinate-aware adapter. Reduced-Gaussian filtering uses direct
DUCC; a non-triangular regular-grid spectral band is rejected because the
released public spharmgrid API exposes triangular selections.
Filtering is optional: both lmin and lmax must be supplied.
spectral_taper=1.0 retains the requested band
without an additional high-wavenumber coefficient taper. Spatial boundary
tapering through taper_points is a separate operation.
Global spherical-harmonic filtering is a spatial -> spectral -> spatial operation: spherical-harmonic coefficients are an intermediate representation, and feature identification operates on the synthesized spatial grid rather than on spectral coefficients.
TRACK's wind-derived vorticity route is implemented by
compute_vorticity().
PyStormTracker uses the public spharmgrid kinematics operation for default
rectangular xarray fields and for NumPy input through its coordinate-aware
adapter. Explicit-lmax calculations use direct ducc0 spin-1 vector
spherical harmonics because the released public same-grid vector composition
does not reproduce the established PST result; unsupported geometries are
rejected.
Hodges (1994) identifies coherent thresholded objects before selecting feature points. The following diagram is the conceptual sequence; the implementation path used to realize the object/candidate stages depends on the selected refinement workflow as described below.
flowchart LR
FIELD[/Input or derived spatial field/]
FIELD --> PREP["Optional preprocessing<br/>boundary taper / spectral filter"]
PREP --> GRID[/Prepared spatial grid/]
GRID --> THRESH[Apply object threshold]
GRID --> CONT[" "]
THRESH --> OBJ[Identify connected objects]
OBJ --> SIZE{Object retained?}
SIZE -->|no| DROP((discard))
SIZE -->|yes| EXT[Find object-local extrema]
EXT --> CAND((Candidate feature points))
classDef continuation fill:transparent,stroke:transparent,color:transparent;
class CONT continuation;
TRACK's frame-level workflow is implemented by
threshold().
Threshold membership is inclusive in
arrayd():
a normalized value equal to the threshold is retained. PyStormTracker therefore
uses >= for maxima and the sign-equivalent <= rule for minima.
TRACK segments the thresholded field with the hierarchical/quad-tree procedure
in
hierarc_segment().
For the default global periodic bspline workflow, PyStormTracker uses a
TRACK-shaped rectangular candidate path that preserves the explicit cyclic
endpoint, seam-object merging, source candidate ordering, and adjacent-extrema
grouping before SMOOPY/GDFP refinement. Other supported workflows use
PyStormTracker's iterative label-propagation CCL while preserving the supported
TRACK connectivity semantics. Global longitude wraps; projected and regional
grids do not.
min_object_grid_points is the minimum retained object size. TRACK removes an
object when point_num <= filt_pt_num in
object_filter(),
so min_object_grid_points=N corresponds to TRACK filt_pt_num=N-1.
TRACK's object-local feature-point search is implemented by
object_local_maxs().
It uses a 3x3 object-local neighborhood. exclude_boundary_extrema=True
corresponds to the source b_exc behavior; TRACK modes with tf >= 4 also
group adjacent tied extrema before later refinement. In PyStormTracker,
adjacent-extremum grouping is automatically part of the bspline,
spherical_bspline, and quadratic candidate workflows;
group_adjacent_extrema=True exposes optional grouping for the grid
refinement path.
Hodges (1995) extends feature-point location on the sphere using cubic
interpolation and local optimization. TRACK contains rectangular SMOOPY and
spherical SPHERY spline workflows. PyStormTracker's default
feature_refinement="bspline" follows the rectangular SMOOPY path used by
the reconciled TRACK workflow.
The diagram below shows the rectangular compatibility path. Periodic seam
restart and coordinate constraints are part of TRACK's non_lin_opt() wrapper
around GDFP rather than properties of GDFP itself.
flowchart LR
GRID[/Prepared spatial grid/]
CAND((Candidate feature points))
GRID --> SPLINE["FITPACK / SMOOPY<br/>B-spline surface"]
CAND --> OPT["TRACK non_lin_opt<br/>constraints + seam restart + GDFP"]
SPLINE --> OPT
OPT --> OK{Converged?}
OK -->|yes| REF[Refined position / value]
OK -->|no| ORIG[Original position / raw value]
REF --> DUP[Duplicate / DUFF handling]
ORIG --> DUP
DUP --> OUT["Refined feature points"]
TRACK's relevant source path is:
spline_smooth()provides the interactive B-spline workflow;surfit()selects SMOOPY or SPHERY and obtains the smoothing factor;smoopy_c()wraps the rectangular Dierckx spline routine;non_lin_opt()prepares coordinate constraints, performs periodic seam restart, handles optimizer failure, and suppresses duplicate refined extrema;gdfp_optimize()performs constrained Goldfarb-Davidon-Fletcher-Powell optimization; andupdate_h()updates the constraint/Hessian-like state used by that optimizer.
TRACK asks interactively for its B-spline smoothing factor and therefore has no
single non-interactive smoothing default. PyStormTracker uses
bspline_smoothing=0.0 by default, corresponding to interpolation. SciPy's
FITPACK implementation constructs the spline; PyStormTracker extracts knots and
coefficients and performs repeated evaluation and optimization in NumPy/Numba.
The mathematical spline lineage is Dierckx/FITPACK; the numerical interface is
provided by SciPy.
HodgesTracker(track_smoopy_optimization_scale=...) controls numerical scaling
in the rectangular GDFP optimizer. Its normal default is 1.0; the historical
TRACK-compatible validation setting is 0.01. It does not rescale stored fields
or detection thresholds. TRACK's line search is not invariant to this numerical
scaling.
The rectangular source-compatible path also preserves several non-obvious
non_lin_opt() semantics: a failed optimization retains the original grid
extremum and raw field value; periodic endpoint handling can trigger a seam
restart; and every rectangular candidate, including one whose optimization
failed, subsequently participates in source-order duplicate/DUFF handling.
These details are implementation compatibility, not new scientific criteria.
TRACK's spherical spline interface is exposed through
sphery_c().
PyStormTracker's spherical_bspline shares this Dierckx/FITPACK spline lineage,
but its optimizer is a PyStormTracker extension. It uses tangent-space
coordinates on
spherical_quadratic is likewise a PyStormTracker candidate-local tangent-space
extension using spherical logarithm/exponential maps. quadratic and grid
remain explicit alternatives.
PyStormTracker stores raw_value, object_gridcell_area_km2, and
object_moment_* diagnostics with final tracks. These are PyStormTracker
second-moment summaries, not equivalents of TRACK's optional object-shape
workflow. TRACK prepares that workflow in
boundary_find(),
shape_setup(),
and
anisotropy.c.
No direct TRACK-equivalence claim is made for the current PyStormTracker
morphology variables.
For three consecutive real feature points, Hodges (1999, Eq. 6) combines changes in tangent direction and displacement magnitude. In PyStormTracker's notation,
\psi =
0.5 w_1 \left(1 - \hat{\mathbf T}_1 \cdot \hat{\mathbf T}_2\right)
+ w_2 \left(1 - \frac{2\sqrt{d_1 d_2}}{d_1 + d_2}\right).
TRACK applies the directional 0.5 normalization when reading w1 in
mge_tracks(),
evaluates the spherical expression in
geod_dev(),
and dispatches real/phantom-point behavior through
devn().
The 0.5 factor normalizes the directional term, whose unscaled range is
0--2, so the directional contribution lies in 0--1 before weighting.
Before initialization, TRACK removes feature points that cannot connect within
the allowed displacement to the next frame, or to the previous frame if no
forward candidate qualifies. The displacement comparison is inclusive in
feature_pt_filter().
TRACK then creates a paired real/all-phantom workspace in
initialize_mge().
The nearest-candidate scan uses dist <= distm, so exact distance ties select
the later source-order candidate. PyStormTracker preserves this source-order
behavior because initialization can affect later equal-cost exchanges.
The following diagram shows the TRACK-shaped MGE algorithm within one temporal segment. PyStormTracker's segment planning and splicing are execution orchestration and are described separately below.
flowchart LR
DET["Refined feature points"]
DET --> PREF[Feature-point prefilter]
PREF --> INIT["Real / phantom workspace<br/>initialization"]
INIT --> FWD["Forward MGE sweeps<br/>until stable"]
FWD --> BWD{"Backward stage<br/>permitted?"}
BWD -->|yes| BACK["Backward MGE sweeps<br/>until stable"]
BACK --> NEXT{Another outer iteration?}
BWD -->|no| SPLIT[Split at phantom gaps]
NEXT -->|yes| FWD
NEXT -->|no| SPLIT
SPLIT --> SEG((Segment tracks))
TRACK identifies the implementation as a modified greedy exchange algorithm,
with the Sethi-Jain method and Salari-Sethi occlusion modification named in
mge_tracks.c.
The driver sets tot_term=3. Its outer loop
mge_tracks()
enters only for more than three frames. Within an active direction, complete
sweeps repeat until that direction makes no further exchange before control can
switch direction. A forward stage is permitted on each outer iteration, while
the backward stage is permitted only while tot_count < tot_term; the final
permitted outer iteration is therefore forward-only. PyStormTracker's
mge_max_iterations=3 reproduces this algorithmic bound; it is not a generic
timeout.
When a four-knot adaptive table is active, each active direction first applies
the source-shaped phantom-gap and directional constraint handling before its
MGE sweep. The directional exchange stages themselves are implemented by
fel_mge()
and
bel_mge().
PyStormTracker applies this TRACK-shaped linking procedure independently within overlapping temporal segments and then deterministically splices the resulting segment tracks. This segmentation/splicing is PyStormTracker execution architecture rather than part of the TRACK MGE algorithm shown above; see Architecture.
Hodges Dask execution has three independent controls:
| Control | Unit of concurrency | Default when omitted |
|---|---|---|
frame_workers |
concurrent frame tasks, including lazy source read, preprocessing, detection, and refinement | available process CPU concurrency |
sht_threads |
Threads per active spherical-harmonic transform; passed to spharmgrid or direct DUCC according to the path | one per active Dask/MPI transform |
mge_workers |
concurrent independent MGE segment-linking tasks | available process CPU concurrency |
segment_frames=62 and the two-frame overlap remain scientific segmentation
parameters; they are independent of mge_workers. MGE is not internally
parallelized.
The resolution helpers in pystormtracker.backends own these defaults.
Serial and MPI execution do not use Dask frame or MGE worker pools, so explicit
frame_workers and mge_workers values are rejected there. Explicit
sht_threads remains meaningful for serial SHT and for rank-local MPI SHT.
HodgesTracker no longer accepts workers; use frame_workers, sht_threads,
and mge_workers. SimpleTracker and HealpixTracker continue to accept
workers.
Supported rectangular SHT, triangular-band regridding, and default kinematics
calls use the public spharmgrid per-operation sht_threads argument. The
reduced-grid, HEALPix, polar, regional-DCT, and explicit-lmax vector paths
pass the resolved value to DUCC and use the existing direct pool
configuration. Native
environment values are logged at DEBUG level with the resolved execution
configuration.
Accepted exchanges are checked against displacement constraints. TRACK's
ub_disp()
provides the source upper-bound calculation. When a link fails, the directional
logic in
track_fail()
moves only the contiguous real section on the failing side into the first
compatible empty workspace interval. There is no separate generic bulk failure
cleanup after the bounded MGE loop.
After MGE, TRACK calls
track_split()
to separate real sections divided by phantom gaps. PyStormTracker performs the
same logical finalization before packing real trajectories into Tracks.
The following states are distinct:
- a phantom entry is an internal
feature_id=-1assignment at an existing workspace time; - an all-phantom workspace row is exchange workspace allocated alongside a real row;
- an existing input time can legitimately contain no detected feature; and
- a known missing input frame is a temporal jump between observed source frames, represented by temporal-gap metadata rather than by synthesizing an empty source frame.
TRACK records the number of known missing source frames on the preceding
observed frame. PyStormTracker derives that count from finalized source times
when time_step is known and never synthesizes the missing timestamps. For
missing_frame_parameters, row
min(nmiss, n_rows - 1) selects the TRACK-style (dmax, phimax) pair. Multiple
parameter rows therefore require a declared cadence; inferring cadence from the
shortest observed interval cannot distinguish missing frames when every
observed interval is already larger than the true source cadence.
The compact PyStormTracker zone and adaptive-table APIs describe one table.
TRACK's main MGE workflow can associate separate zone/adaptive tables with
missing-frame parameter rows; PyStormTracker does not silently reuse one table
for every row and rejects unsupported combinations. TRACK's separate legacy
post-link checker,
tr_miss_frame(),
is documented but not exposed as a separate public PyStormTracker workflow.
max_missing_steps is a separate PyStormTracker topology extension restricting
internal phantom runs during proposed exchanges. Leading and trailing phantoms
do not count toward it. The default None preserves TRACK MGE behavior with
respect to this extension.
Hodges (1999, Section 5) motivates spatially varying displacement limits and speed-dependent smoothness.
TRACK reads regional limits with
read_zones().
With a nonempty zone table, every used real feature endpoint must lie in a
zone; there is no silent fallback. Nonnegative longitude definitions are
interpreted in the 0..360 convention, boundaries are inclusive, and the
per-link displacement limit is the average of the two endpoint-zone limits.
TRACK also resets its global displacement value to the maximum zone value;
PyStormTracker mirrors that behavior. An empty table selects the global dmax
path.
TRACK requires four displacement cutoffs and four corresponding phimax values
in
read_adptp()
and precomputes three linear segments. The actual constraint is evaluated from
the mean of the two adjacent displacements by
phi().
PyStormTracker therefore accepts either a disabled table or four finite,
strictly increasing displacement knots. When active, static phimax is raised
to at least the maximum adaptive value. TRACK's directional zonal/adaptive
post-filter is implemented in
tr_zonal_filter().
PyStormTracker's filter_rsplice() implements the supported lifetime and
displacement semantics of TRACK's post-tracking splice workflow. This is
post-processing and should not be conflated with MGE trajectory construction.
TRACK's workflow driver is
splice_tracks().
For displacement filtering,
disp_filter()
can use cumulative travel distance or start-to-finish separation and removes a
track only when displacement is strictly less than the requested threshold.
A track exactly on the boundary is retained. TRACK's lower-level point-distance
helper is
measure().
Trajectory intercomparison is not part of the Hodges 1994/1995/1999 tracking algorithm. PyStormTracker uses TRACK's later ENSEMBLE utilities as the source reference for its default eligibility rule.
TRACK 1.5.4 defines
TOLMATCH = 2.0 and TOLNUM = 0.6.
For each candidate pair, compare_ensemble2.c computes
and accepts the pair when the overlap fraction is at least TOLNUM and the
selected separation is no greater than TOLMATCH; see
compare_ensemble2.c.
The source defaults are therefore 60% symmetric temporal overlap and 2 degrees
separation. TRACK permits mean or minimum separation; PyStormTracker uses
whole-overlap mean geodesic separation.
toverlap()
finds the common interval, while
trdist()
computes mean or minimum concurrent-point geodesic separation. TRACK then
selects the closest eligible candidate independently for each reference track,
so candidate reuse is possible.
The thresholds have published lineage but should not be attributed to the
original Hodges method papers. Hodges et al. (2003) used at least 60% temporal
overlap with a tighter 0.5-degree mean-separation condition. Wang, Swail, and
Zwiers (2006), comparing cyclone tracks on 2.5-degree unfiltered MSLP, retained
the 60% overlap criterion and used a 2.0-degree separation threshold. TRACK's
ENSEMBLE utility subsequently carries 0.6 and 2.0 as defaults.
PyStormTracker provides three pairing policies after applying eligibility:
nearestfollows TRACK's directed closest-eligible-candidate policy;mutual_nearestretains reciprocal nearest pairs, following the reciprocal-neighbor idea used by Blender and Schubert (2000), while retaining PyStormTracker's TRACK-style eligibility definition; andglobal_assignmentis a PyStormTracker deterministic one-to-one extension that maximizes matched-pair count, then total temporal overlap, then minimizes total mean separation.
topology_identical is a PyStormTracker validation diagnostic, not a Hodges
or TRACK matching criterion. For an already matched pair it is true only when
the complete timestamp arrays are exactly equal. It has no geographic or
intensity tolerance. same_time_range and same_point_count are separate
reported diagnostics.
Source-stage probes reproduce selected TRACK detection, workspace, MGE, constraint, failure, and splitting behavior directly. Those tests establish implementation correspondence for the stated configurations; they do not make TRACK independent scientific ground truth.
The strongest broad trajectory comparison currently available starts from
full-year 2024 six-hourly ERA5 MSLP on the F320 Gaussian source grid, retains
the TRACK T6--42 spectral band, reconstructs it onto the T42 Gaussian tracking
grid, and uses rectangular bspline refinement. Both implementations consume
the same TRACK-produced filtered spatial field, so this isolates the tracking
implementation and does not establish independent raw-ERA5
spectral-preprocessing identity.
| 2024 ERA5 MSLP, T6--42 on T42 grid | Raw | RSPLICE-filtered |
|---|---|---|
| TRACK tracks | 7,761 | 1,471 |
| PyStormTracker tracks | 7,859 | 1,470 |
| TRACK points | 60,654 | 30,998 |
| PyStormTracker points | 60,883 | 31,015 |
| Global-assignment F1 | 0.9921 | 0.9983 |
| Topology-identical matched pairs | 7,708 | 1,453 |
Directed nearest matching covers 7,750 of 7,761 TRACK raw trajectories
(99.86%). After TRACK-compatible RSPLICE filtering, global one-to-one
assignment identifies 1,468 common storms among approximately 1,470 tracks,
with F1 0.9983. topology_identical means only that a matched pair has exactly
the same complete timestamp sequence; it does not require identical center
coordinates or intensities.
| Area | TRACK/Hodges relationship | PyStormTracker status |
|---|---|---|
| Object identification | Hodges method; TRACK source correspondence | TRACK-shaped rectangular path for global bspline; iterative CCL for other supported paths |
| Rectangular B-spline centers | TRACK SMOOPY + coordinate-space GDFP | default bspline path |
| Spherical quadratic | not TRACK | explicit PyStormTracker extension |
| Intrinsic spherical B-spline optimization | not TRACK | explicit PyStormTracker extension |
| MGE workspace and exchange control | TRACK source correspondence | source-shaped Python/Numba implementation |
max_missing_steps |
not TRACK MGE behavior | optional extension; disabled by default |
| Directed overlap/separation comparison | TRACK ENSEMBLE utility lineage | nearest |
| Reciprocal nearest comparison | later intercomparison literature | mutual_nearest |
| Global one-to-one assignment | not TRACK | PyStormTracker extension |
| Exact timestamp-sequence identity | not TRACK | PyStormTracker validation diagnostic |
The narrative above links source where it affects scientific or software
behavior. This consolidated index replaces the former separate source-map page.
All links target the immutable TRACK-1.5.4 tag.
| Stage | TRACK 1.5.4 source | What it establishes | PyStormTracker relationship |
|---|---|---|---|
| Threshold workflow | threshold() |
Frame-level threshold/object driver | Workflow reference |
| Threshold membership | arrayd() |
Inclusive threshold test | Matching max/min semantics |
| Object segmentation | hierarc_segment() |
Hierarchical object segmentation | Rectangular bspline path preserves TRACK-shaped representation; other paths use PST CCL |
| Object construction | form_objects() |
Converts segmentation into object structures | Source behavior reference |
| Object-size filtering | object_filter() |
point_num <= filt_pt_num removal |
Maps to min_object_grid_points |
| Object-local extrema | object_local_maxs() |
3x3 extrema, boundary option, grouping | Detector source reference |
| Spectral filtering | spectral_filter() |
Spatial spectral-filter workflow | spharmgrid-backed implementation differs numerically |
| Spectral wrapper | spec_filt() |
Interactive filtering orchestration | Not copied by PST |
| Hoskins taper | hoskins_filt() |
Exponential coefficient taper | Correct source owner |
| Wind-derived vorticity | compute_vorticity() |
TRACK wind-to-vorticity workflow | PST uses spin-1 harmonics |
| Spline dispatch | surfit() |
SMOOPY/SPHERY selection and smoothing | Workflow reference |
| Rectangular spline | smoopy_c() |
Dierckx SMOOPY interface | bspline compatibility lineage |
| Spherical spline | sphery_c() |
Dierckx SPHERY interface | Spline lineage only for spherical_bspline |
| Nonlinear refinement driver | non_lin_opt() |
Constraints, seam restart, failure, duplicate handling | Rectangular compatibility behavior |
| GDFP optimizer | gdfp_optimize() |
Constrained variable-metric optimization | Native PST rectangular implementation |
| Constraint/Hessian update | update_h() |
Constraint-state update | Source-mapped implementation detail |
| Spline objective support | func.c |
Objective/spline evaluation support | Source-mapped implementation detail |
| Feature-point prefilter | feature_pt_filter() |
Inclusive adjacent-frame dmax eligibility |
Source-mapped |
| MGE initialization | initialize_mge() |
Greedy initialization and paired phantom rows | Source-mapped workspace |
| Spherical MGE cost | geod_dev(), devn() |
Real-point cost and phantom penalty | Source-mapped |
| MGE scheduler | mge_tracks() |
Three-stage outer control and final split | mge_max_iterations=3 source |
| Forward/backward exchange | fel_mge(), bel_mge() |
Directional MGE sweeps | hodges/mge.py lineage |
| Upper displacement bound | ub_disp() |
Exchange displacement bound | Source-mapped |
| Failure and final split | track_fail(), track_split() |
Failure relocation and phantom-gap splitting | Source-mapped |
Regional dmax |
read_zones() |
Regional displacement table | dmax_zones lineage |
| Adaptive smoothness | read_adptp(), phi() |
Four-knot piecewise-linear phimax |
adaptive_smoothness lineage |
| Directional constraints | tr_zonal_filter() |
Directional zonal/adaptive filtering | Source reference |
| Missing-frame checker | tr_miss_frame() |
Legacy post-link missing-frame workflow | Documented, not exposed as separate API |
| RSPLICE workflow | splice_tracks() |
TRACK postprocessing driver | Workflow lineage |
| RSPLICE displacement | disp_filter() |
Travel/end-to-end displacement and strict removal test | filter_rsplice() semantics |
| Distance helper | measure() |
Point separation helper | Source helper |
| Object boundary | boundary_find() |
Optional object-boundary representation | Not equivalent to PST moment diagnostics |
| Object shape setup | shape_setup() |
Optional shape workflow setup | Not equivalent to PST moment diagnostics |
| Anisotropy | anisotropy.c |
TRACK anisotropy/shape workflow | Not implemented as direct equivalent |
| Comparison defaults | compare_ensemble2.c |
TOLMATCH=2.0, TOLNUM=0.6 |
Default pair eligibility |
| Comparison eligibility | compare_ensemble2.c |
Overlap equation, thresholds, directed nearest selection | nearest lineage |
| Overlap/separation helpers | toverlap(), trdist() |
Common interval and mean/minimum separation | PST aligns exact common timestamps and uses mean separation |
- Hodges, K. I., 1994: A General Method for Tracking Analysis and Its Application to Meteorological Data. Mon. Wea. Rev., 122, 2573–2586. doi:10.1175/1520-0493(1994)122<2573:AGMFTA>2.0.CO;2.
- Hodges, K. I., 1995: Feature Tracking on the Unit Sphere. Mon. Wea. Rev., 123, 3458–3465. doi:10.1175/1520-0493(1995)123<3458:FTOTUS>2.0.CO;2.
- Hodges, K. I., 1999: Adaptive Constraints for Feature Tracking. Mon. Wea. Rev., 127, 1362–1373. doi:10.1175/1520-0493(1999)127<1362:ACFFT>2.0.CO;2.
- Dierckx, P., 1993: Curve and Surface Fitting with Splines. Oxford University Press.
- Smith, S. T., 1994: Optimization Techniques on Riemannian Manifolds. Fields Institute Communications, 3, 113–136.
- Edelman, A., T. A. Arias, and S. T. Smith, 1998: The Geometry of Algorithms with Orthogonality Constraints. SIAM J. Matrix Anal. Appl., 20(2), 303–353. doi:10.1137/S0895479895290954.
- Huang, W., K. A. Gallivan, and P.-A. Absil, 2015: A Broyden Class of Quasi-Newton Methods for Riemannian Optimization. SIAM J. Optim., 25(3), 1660–1685. doi:10.1137/140955483.
- Blender, R., and M. Schubert, 2000: Cyclone Tracking in Different Spatial and Temporal Resolutions. Mon. Wea. Rev., 128, 377–384.
- Hodges, K. I., B. J. Hoskins, J. Boyle, and C. Thorncroft, 2003: A Comparison of Recent Reanalysis Datasets Using Objective Feature Tracking: Storm Tracks and Tropical Easterly Waves. Mon. Wea. Rev., 131, 2012–2037. doi:10.1175/1520-0493(2003)131<2012:ACORRD>2.0.CO;2.
- Wang, X. L., V. R. Swail, and F. W. Zwiers, 2006: Climatology and Changes of Extratropical Cyclone Activity: Comparison of ERA-40 with NCEP-NCAR Reanalysis for 1958-2001. J. Climate, 19, 3145–3166. doi:10.1175/JCLI3781.1.