mwyau · mwyau · Mar 25, 2026 · Mar 25, 2026 · Mar 25, 2026 · Mar 25, 2026
diff --git a/README.md b/README.md
@@ -59,10 +59,9 @@ Full documentation, including API references and advanced usage examples, is ava
 - Python 3.11+
 - (Optional) OpenMPI for MPI support.
 - **SHT Backends**: Supported engines for filtering and derivatives:
-  - `ducc0`: **Core dependency**. High-precision C++ backend. Default for derivatives.
-  - `pyshtools`: **Core dependency**. Includes Fortran backend.
-  - `shtns`: (Optional) High-performance C library (`pip install PyStormTracker[pre]`). Recommended for large datasets.
-- **Windows**: GRIB support is experimental. SHTns is not supported on Windows; `ducc0` or `shtools` will be used automatically.
+  - `ducc0`: **Core dependency**. High-precision C++ library (Distinctly Useful Code Collection) providing high performance spherical harmonic transforms.
+  - `shtns`: (Optional) High-performance C library (`pip install PyStormTracker[shtns]`). Recommended for large datasets.
+- **Windows**: GRIB support is experimental. SHTns is not supported on Windows; `ducc0` will be used automatically.
 
 ### From PyPI
 You can install the latest stable version of PyStormTracker directly from PyPI:
@@ -135,7 +134,7 @@ stormtracker -i data.nc -v msl -o my_tracks.txt
 | `--threshold` | `-t` | Intensity threshold for feature detection. |
 | `--filter-range` | | Spectral filter range (min-max). Default '5-42'. |
 | `--no-filter` | | Disable default T5-42 spectral filtering. |
-| `--engine` | | SHT backend: `auto`, `shtns`, `ducc0`, `shtools`. |
+| `--sht-engine` | | SHT backend: `auto`, `shtns`, `ducc0`. |
 | `--num` | `-n` | Number of time steps to process. |
 | **Performance** | | |
 | `--backend` | `-b` | `serial`, `dask`, or `mpi`. Auto-detected by default. |

diff --git a/docs/changelog.md b/docs/changelog.md
@@ -2,7 +2,7 @@
 
 ## v0.5.0.dev
 ### Features
-- **Spectral Backends**: Integrated support for `shtns`, `ducc0`, and `shtools` for all spherical harmonic operations.
+- **Spectral Backends**: Integrated support for `shtns` and `ducc0` for all spherical harmonic operations.
 - **High-Precision Derivatives**: New `vodv` module for computing relative vorticity and divergence using spin-1 vector harmonics.
 - **Planetary Constants**: Standardized Earth radius to 6,371,220 m across derivatives and tracking geometry.
 - **pst-convert & JSON Support**: New utility for format conversion and a GPU-optimized JSON format for track data.
@@ -15,7 +15,7 @@
 
 ### Performance
 - **SHTns Optimization**: Optimized thread handling for spherical harmonic transforms.
-- **Backend Fallback**: Implemented `pyshtools` as the default backend for spherical harmonic filtering, with `shtns` as an optional, high-performance alternative.
+- **Backend Fallback**: Implemented native `ducc0` usage for spherical harmonic filtering and transformations, removing the previous wrapper dependency. `shtns` remains an optional alternative.
 - **ARM64 Support**: Restored full support for **ARM64** architectures in Docker and CI by making `shtns` an optional dependency.
 
 ### CI/CD & Testing

diff --git a/docs/hodges.md b/docs/hodges.md
@@ -11,12 +11,12 @@ This document details the architecture, mathematical implementation, and design
 ## 1. Feature Identification Design
 
 ### 1.1 Preprocessing (Spectral Filtering & Derivatives)
-**Design Choice**: Integrated native spherical harmonic filtering (e.g., T42 truncation) and high-precision derivative calculation using the `shtns`, `ducc0`, or `shtools` backends.
+**Design Choice**: Integrated native spherical harmonic filtering (e.g., T42 truncation) and high-precision derivative calculation using the `shtns` or `ducc0` backends.
 - **References**: `spec_filt.c`, `uv2vr.c`.
-- **CLI Options**: `--filter-range` (e.g., `5-42`), `--no-filter`, `--engine` (auto/shtns/ducc0/shtools), and `--taper`.
+- **CLI Options**: `--filter-range` (e.g., `5-42`), `--no-filter`, `--engine` (auto/shtns/ducc0), and `--taper`.
 
 **Reasoning**: 
-Original TRACK workflows typically require offline spectral filtering to remove the planetary background and high-frequency noise. `PyStormTracker` incorporates this directly into its preprocessing module for on-the-fly execution. By default, the Hodges algorithm applies a T5-42 band-pass filter unless `--no-filter` is specified. The system also supports high-precision **Relative Vorticity** and **Divergence** calculation from wind components using spin-1 vector harmonics, ensuring bit-wise parity with NCL/Spherepack when using the `ducc0` or `shtools` backends.
+Original TRACK workflows typically require offline spectral filtering to remove the planetary background and high-frequency noise. `PyStormTracker` incorporates this directly into its preprocessing module for on-the-fly execution. By default, the Hodges algorithm applies a T5-42 band-pass filter unless `--no-filter` is specified. The system also supports high-precision **Relative Vorticity** and **Divergence** calculation from wind components using spin-1 vector harmonics, ensuring bit-wise parity with NCL/Spherepack when using the `ducc0` backend.
 
 ### 1.2 Object-Based Detection
 **Design Choice**: Feature detection is implemented as a multi-stage pipeline: `Thresholding -> Connected Component Labeling (CCL) -> Object Filtering -> Local Extrema`.

diff --git a/pyproject.toml b/pyproject.toml
@@ -28,7 +28,7 @@ dependencies = [
     "h5netcdf[h5py]>=1.0.0",
     "numpy>=1.24.0",
     "numba>=0.60.0",
-    "pyshtools[ducc]>=4.10.4",
+    "ducc0>=0.33.0",
     "xarray>=2024.9.0",
 ]
 
@@ -123,7 +123,7 @@ disallow_untyped_calls = false
 exclude = ["worktrees", ".venv", "notebooks"]
 
 [[tool.mypy.overrides]]
-module = ["dask.*", "numba.*"]
+module = ["dask.*", "numba.*", "ducc0", "ducc0.*"]
 ignore_missing_imports = true
 
 [tool.pytest.ini_options]
@@ -157,4 +157,4 @@ exclude_lines = [
 
 [tool.uv]
 link-mode = "copy"
-no-binary-package = ["shtns"]
+no-binary-package = ["ducc0", "shtns"]
diff --git a/src/pystormtracker/cli.py b/src/pystormtracker/cli.py
@@ -57,6 +57,7 @@ def run_tracker(
     lmax: int = constants.LMAX_DEFAULT,
     taper_points: int = constants.TAPER_DEFAULT,
     overlap: int = model_constants.OVERLAP_DEFAULT,
+    sht_engine: Literal["auto", "shtns", "ducc0"] = "auto",
 ) -> None:
     """Orchestrates the storm tracking process from the CLI."""
     timer: dict[str, float] = {}
@@ -152,6 +153,7 @@ def run_tracker(
         lmax=lmax,
         taper_points=taper_points,
         overlap=overlap,
+        sht_engine=sht_engine,
     )
 
     # Export Phase
@@ -238,6 +240,13 @@ def parse_args() -> Namespace:
     )
     # Default is determined in main() based on algorithm
 
+    general.add_argument(
+        "--sht-engine",
+        choices=["auto", "shtns", "ducc0"],
+        default="auto",
+        help="SHT backend for filtering and derivatives. Default 'auto'.",
+    )
+
     general.add_argument(
         "-n", "--num", type=int, help="Number of time steps to process."
     )
@@ -455,6 +464,7 @@ def main() -> None:
         lmax=lmax,
         taper_points=args.taper,
         overlap=args.overlap,
+        sht_engine=args.sht_engine,
     )
 
 

diff --git a/src/pystormtracker/hodges/tracker.py b/src/pystormtracker/hodges/tracker.py
@@ -80,6 +80,7 @@ def preprocess_standard_track(
         lmin: int = constants.LMIN_DEFAULT,
         lmax: int = constants.LMAX_DEFAULT,
         taper_points: int = constants.TAPER_DEFAULT,
+        sht_engine: Literal["auto", "shtns", "ducc0"] = "auto",
     ) -> xr.DataArray:
         """
         Applies standard TRACK preprocessing: Tapering -> Spherical Harmonic Filter.
@@ -94,7 +95,7 @@ def preprocess_standard_track(
             data = cast(xr.DataArray, taper.filter(data))
 
         # 2. Spectral Filtering
-        spectral_filter = SpectralFilter(lmin=lmin, lmax=lmax)
+        spectral_filter = SpectralFilter(lmin=lmin, lmax=lmax, sht_engine=sht_engine)
         data = spectral_filter.filter(data)
 
         return data
@@ -161,6 +162,7 @@ def track(
         lmin: int = constants.LMIN_DEFAULT,
         lmax: int = constants.LMAX_DEFAULT,
         taper_points: int = constants.TAPER_DEFAULT,
+        sht_engine: Literal["auto", "shtns", "ducc0"] = "auto",
         **kwargs: float | int | str | None,
     ) -> Tracks:
         """
@@ -181,6 +183,8 @@ def track(
             min_points: Minimum grid points per object.
             filter: If True, apply spectral filtering.
             lmin, lmax: Spectral truncation range (default T5-42).
+            taper_points: Boundary tapering points.
+            sht_engine: SHT backend engine.
         """
         import timeit
 
@@ -199,7 +203,13 @@ def track(
         data_xr = detector_peek.get_xarray(start_time, end_time)
 
         if filter:
-            data_xr = self.preprocess_standard_track(data_xr, lmin=lmin, lmax=lmax)
+            data_xr = self.preprocess_standard_track(
+                data_xr,
+                lmin=lmin,
+                lmax=lmax,
+                taper_points=taper_points,
+                sht_engine=sht_engine,
+            )
         t1 = timeit.default_timer()
         print(f"    [Serial] Preprocessing time: {t1 - t0:.4f}s")
 

diff --git a/src/pystormtracker/models/tracker.py b/src/pystormtracker/models/tracker.py
@@ -41,5 +41,6 @@ def track(
         lmin: int = 5,
         lmax: int = 42,
         taper_points: int = 0,
+        sht_engine: Literal["auto", "shtns", "ducc0"] = "auto",
         **kwargs: float | int | str | None,
     ) -> Tracks: ...