SpiPy is a Python implementation of SPIRES (SPectral Inversion of REflectance from Snow), a spectral unmixing algorithm for analyzing snow reflectance data. The original MATLAB implementation is described in this IEEE paper.
Core Purpose: Invert satellite reflectance spectra to retrieve snow properties (grain size, dust concentration, fractional snow-covered area) using lookup tables generated from Mie-scattering theory.
This is a performance-critical scientific computing project with a hybrid architecture:
- Python layer: High-level API, data I/O, coordinate transformations
- C++ layer: Performance-critical bottlenecks (3000x speedup achieved)
- SWIG: Bindings between Python and C++
Critical files:
- spires/spires.cpp: Core C++ implementations (interpolation, optimization)
- spires/interpolator.py: LUT interpolation interface
- spires/invert.py: Main inversion algorithm
- spires/process.py: High-level processing pipeline
-
Lookup Table (LUT) Interpolator
- Interpolates pre-computed Mie-scattering reflectance spectra
- 4D interpolation: (band, solar_angle, dust, grain_size)
- C++ implementation using SWIG is ~3000x faster than pure Python
-
Spectral Inversion
- Nonlinear optimization to match observed reflectance to modeled reflectance
- Uses NLopt library (COBYLA algorithm) in C++
- Note: SLSQP doesn't work in C++ implementation (see issues in README)
-
Data Processing
- Handles MODIS, Sentinel-2, Landsat data
- Coordinate transformations, cloud masking, temporal analysis
- Uses xarray, rioxarray, gdal for geospatial operations
You MUST use conda-forge for dependencies, not apt:
nlopt: Optimization library (apt version missing C++ headers!)swig: Generate Python/C++ bindingsgcc,g++: C++ compilation- See README.md lines 32-48 for full installation
Core: numpy, h5py, scipy, xarray, netCDF4 Optional: gdal, geopandas, matplotlib, dask, rioxarray, pyproj
# Build SWIG extension (required after C++ changes)
python3 setup.py build_ext --inplace
# Install package
pip3 install .
# Build wheel
python -m build --wheel- Performance is paramount: This code has been heavily optimized. Check README.md lines 160-183 for performance benchmarks before/after changes.
- SWIG bindings: Changes to C++ files require rebuilding with
python3 setup.py build_ext --inplace - Test after C++ changes: Run
pytest --doctest-modulesand check examples notebooks
- Dimension ordering: C++ uses different array ordering than NumPy (see issues in README)
- Scaling in optimization: Different coordinate scales require problem scaling (see line 344-346)
- COBYLA vs SLSQP: SLSQP doesn't work in C++, using LN_COBYLA instead
- Git LFS: Test data stored in LFS, need
git lfs install
spires/
├── __init__.py # Package exports
├── interpolator.py # LUT interpolation wrapper
├── invert.py # Main inversion functions
├── process.py # High-level processing
├── core.py # Core utilities
├── spires.cpp # C++ optimized code
├── spires.h # C++ header
├── spires_wrap.cpp # SWIG-generated (don't edit!)
└── cobyla.cpp # COBYLA implementation
tests/
└── data/ # LUT files, test images (git-lfs)
examples/ # Jupyter notebooks with use cases
doc/ # Sphinx documentation
pytest --doctest-modules # Doctests in source filescd doc/
make htmlSee examples/ for Jupyter notebooks demonstrating:
- Basic inversion workflows
- Sentinel-2/MODIS data processing
- Cloud masking and postprocessing
- Temporal analysis
- Interpolator behavior differs between SWIG and scipy when coordinates aren't linspace
- SLSQP solver doesn't work in C++, using COBYLA
- NumPy data ownership issues may cause minor performance loss
- Potential optimization: Keep R_0 constant when inverting same location over time
- Switch from MATLAB .mat LUT files to HDF5 or similar
- Accept xarray inputs for interpolator and target/background spectra
SPIRES retrieves snow properties from multispectral satellite imagery by:
- Loading pre-computed reflectance lookup tables (Mie scattering theory)
- Iterating through each pixel's observed reflectance spectrum
- Optimizing model parameters (grain size, dust, fsca) to minimize difference from observed
- Using spatial constraints and background spectra for robustness
- fsca: Fractional Snow-Covered Area (0-1)
- grain_size: Effective snow grain radius (μm)
- dust: Dust concentration in snow
- R_0: Background (snow-free) reflectance spectrum
- solar_angle: Solar zenith angle at image acquisition
The C++ optimizations achieved dramatic speedups:
- Interpolation: 3000x faster (1.07 ms → 309 ns)
- Spectrum difference: 1000x faster (1.1 ms → 1 μs)
- Full optimization: 3000x faster (165 ms → 43 μs)
These speedups enable processing entire satellite images (millions of pixels) in reasonable time.
When working on this codebase, feel free to ask:
- "Should this change affect the C++ layer or stay in Python?"
- "Will this impact performance? Should I benchmark?"
- "Are there examples/tests that cover this use case?"
- "Does this need to work with MODIS, Sentinel-2, or both?"