API Reference#

The API Reference provides an overview of all public objects, functions and methods implemented in SpiPy. All classes and function exposed in spires.* namespace plus those listed in the reference are public.

`LutInterpolator`([grid, reflectances, bands, ...])	A LutInterpolator provides a convenient interface to read LUT data from netCDF files and an interface to interpolate ideal snow reflectances.
`get_index`(coordinates, value)	Returns the interpolated index of a coordinate in a numpy array.
`get_index_linspace`(coordinates, value)	Returns the interpolated index of a coordinate in a numpy array assuming that coordinates are a linear space.
`is_linspace`(array)	Checks if a numpy array is linearly spaced.

`index_to_value`(index, coords)	Convert normalized index to coordinate value.
`snow_diff_3`(x, spectrum_target, solar_angle, ...)	Calculate spectral difference for 3-parameter snow model.
`snow_diff_4`(x, spectrum_target, ...)	Calculate spectral difference for 4-parameter snow model.
`snow_diff_softmax`(z, spectrum_target, ...)	Spectral-difference cost in unconstrained (softmax-reparameterized) coordinates.
`speedy_invert`(spectrum_target, ...[, ...])	Inverts the snow reflectance spectrum using nonlinear optimization.
`speedy_invert_array1d`(spectra_targets, ...)	Batch inversion of snow reflectance spectra for 1D arrays of observations.
`speedy_invert_array2d`(spectra_targets, ...)	Batch inversion of snow reflectance spectra for 2D spatial arrays.
`speedy_invert_scipy`(interpolator, ...[, ...])	Invert snow spectra using scipy.optimize.minimize.
`speedy_invert_scipy_normalized`(interpolator, ...)	Invert snow spectra with normalized parameter space.
`speedy_invert_scipy_softmax`(interpolator, ...)	Unconstrained scipy inversion via softmax reparameterization.
`speedy_invert_xarray`(spectra_targets, ...[, ...])	Batch inversion of snow reflectance spectra using xarray DataArrays.

`encode_results`(ds[, fill_value, fsca_scale, ...])	Encode an inversion result Dataset for compact storage.
`speedy_invert_array2d`(spectra_targets, ...)	Batch inversion of snow reflectance spectra for 2D spatial arrays.
`speedy_invert_dask`(spectra_targets, ...[, ...])	Parallel inversion of snow reflectance spectra using Dask and xarray.

`unique_elements_space`(spectra, ...)
`unique_elements_spacetime`(spectra, ...)
`uniquetol_1d`(array[, tol])	Finds unique vectors in an array across dimension 1

`load_lut`(lut_file)
`speedy_invert`(f, spectrum_target, ...[, ...])	Inverts a spectrum and determines f_sca, f_shade, dust_concentration and grain_size.

class spires.interpolator.LutInterpolator(grid=None, reflectances=None, bands=None, solar_angles=None, dust_concentrations=None, grain_sizes=None, lut_file=None)#

Bases: object

A LutInterpolator provides a convenient interface to read LUT data from netCDF files and an interface to interpolate ideal snow reflectances.

interpolate(band, solar_angle, dust_concentration, grain_size)#

Interpolate values using the c++ interpolator. Index lookup is performed in python.

Parameters:

band (int) – Band to interpolate the reflectance for. Has to be a value in self.bands
solar_angle (double) – solar angle to interpolate the reflectance for. Units of solar_angle has to match the units in self.solar_angles (e.g. degrees).
dust_concentration (double) – dust concentration to interpolate the reflectance for. Units of dust_concentration has to match the units in self.dust_concentrations (e.g. ppm)
grain_size (double) – grain size to interpolate the reflectance for. Units of grain_size has to match the units in self.grain_sizes (e.g. micro meters).

Return type:

Interpolated snow reflectance for given band and solar angle, dust concentration and grain size

Examples

>>> import spires
>>> interpolator = spires.LutInterpolator(lut_file='tests/data/lut_sentinel2b_b2to12_3um_dust.mat')
>>> reflectance = interpolator.interpolate(band=1, solar_angle=0, dust_concentration=0.1, grain_size=30)
>>> round(reflectance, 6)
0.992369

interpolate_all(solar_angle, dust_concentration, grain_size)#

Interpolate spectrum using the c++ interpolator for all bands in self.bands. Derive the index values using swig C++ functions (about 5x more efficient than self.interpolate_all_np_index).

Parameters:

solar_angle (double) – solar angle to interpolate the spectrum for. Units of solar_angle has to match the units in self.solar_angles (e.g. degrees).
dust_concentration (double) – dust concentration to interpolate the spectrum for. Units of dust_concentration has to match the units in self.dust_concentrations (e.g. ppm)
grain_size (double) – grain size to interpolate the spectrum for. Units of grain_size has to match the units in self.grain_sizes (e.g. micro meters).

Return type:

The snow reflectance spectrum

Examples

>>> import spires
>>> interpolator = spires.LutInterpolator(lut_file='tests/data/lut_sentinel2b_b2to12_3um_dust.mat')
>>> interpolator.interpolate_all(solar_angle=0, dust_concentration=0.1, grain_size=30)
array([0.99236893, 0.987902  , 0.97464906, 0.96756319, 0.96042174,
       0.94498655, 0.1533866 , 0.18644477, 0.92160772])

interpolate_all_np_index(solar_angle, dust_concentration, grain_size)#

Interpolate spectrum using the c++ interpolator for all bands in self.bands. Derive the index values using numpy functions.

Parameters:

solar_angle (double) – solar angle to interpolate the spectrum for. Units of solar_angle has to match the units in self.solar_angles (e.g. degrees).
dust_concentration (double) – dust concentration to interpolate the spectrum for. Units of dust_concentration has to match the units in self.dust_concentrations (e.g. ppm)
grain_size (double) – grain size to interpolate the spectrum for. Units of grain_size has to match the units in self.grain_sizes (e.g. micro meter).

Return type:

The snow reflectance spectrum

Examples

>>> import spires
>>> interpolator = spires.LutInterpolator(lut_file='tests/data/lut_sentinel2b_b2to12_3um_dust.mat')
>>> interpolator.interpolate_all_np_index(solar_angle=0, dust_concentration=0.1, grain_size=30)
array([0.99236893, 0.987902  , 0.97464906, 0.96756319, 0.96042174,
       0.94498655, 0.1533866 , 0.18644477, 0.92160772])

interpolate_pts(pts)#

Compatability function; interpolate values using the c++ interpolator :param pts: :type pts: array_like with :param - pts[0]: :type - pts[0]: band :param - pts[1]: :type - pts[1]: solar_angle :param - pts[2]: :type - pts[2]: dust_concentration :param - pts[3]: :type - pts[3]: grain_size

Return type:: Interpolated snow reflectance for given band and solar angle, dust concentration and grain size

Examples

>>> import spires
>>> interpolator = spires.LutInterpolator(lut_file='tests/data/lut_sentinel2b_b2to12_3um_dust.mat')
>>> reflectance = interpolator.interpolate_pts([1, 0, 0.1, 30])
>>> round(reflectance, 6)
0.992369

interpolate_scipy(band, solar_angle, dust_concentration, grain_size)#

Interpolate values using the scipy RegularGridInterpolator interpolator

Parameters:

band (int)
solar_angle (double)
dust_concentration (double)
grain_size (double)

Return type:

Interpolated snow reflectance for given band and solar angle, dust concentration and grain size

interpolate_scipy_pts(pts)#

Interpolate values using the scipy RegularGridInterpolator interpolator

Parameters:: pts (array-like)
Return type:: Interpolated snow reflectance for given band and solar angle, dust concentration and grain size

load_mat(lut_file)#

Loads a LutInterpolator from a mat file. The file has to be structured in a fairly idiosyncratic format.

Parameters:: lut_file (string)

make_scipy_interpolator()#

Creates a scipy interpolator object using coordinates and reflectances

Return type:: None

make_scipy_interpolator_legacy()#

Creates a scipy interpolator object using coordinates and reflectances.

In the legacy implementation, the dimensions were (axis 2 and 3 swapped): - dim1 = bands - dim2 = solar_angles - dim3 = grain_sizes - dim4 = dust_concentrations

Return type:: None

to_xarray()#

verify_linspace()#

Verifies that all coordinates are linearly spaced.

Return type:: None

spires.interpolator.get_index(coordinates: array, value: float) → float#

Returns the interpolated index of a coordinate in a numpy array.

Parameters:

coordinates (array-like) – Array of coordinates.
value (float) – Value of the coordinate to find the coordinate index for

Returns:

interpolated_index – The interpolated index in the coordinates for value

Return type:

float

Examples

>>> import spires
>>> coordinates = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
>>> value = 1.5
>>> index = spires.get_index(coordinates, value)
>>> index == 1.5
True

spires.interpolator.get_index_linspace(coordinates: array, value: float) → float#

Returns the interpolated index of a coordinate in a numpy array assuming that coordinates are a linear space.

Parameters:

coordinates (array-like) – Array of coordinates.
value (float) – Value of the coordinate to find the coordinate index for

Returns:

interpolated_index – The interpolated index in the coordinates for value

Return type:

float

spires.interpolator.is_linspace(array)#

Checks if a numpy array is linearly spaced. :param array: :type array: numpy.ndarray

Returns:: is_linspace
Return type:: bool

spires.invert.index_to_value(index, coords)#

Convert normalized index to coordinate value.

Linearly interpolates between coordinate values based on a normalized index in the range [0, 1].

Parameters:

index (float) – Normalized index value between 0 and 1.
coords (numpy.ndarray) – Array of coordinate values to interpolate between.

Returns:

Interpolated coordinate value.

Return type:

float

Notes

Used internally by speedy_invert_scipy_normalized to convert normalized optimization parameters back to physical units.

spires.invert.snow_diff_3(x, spectrum_target, solar_angle, interpolator, shade)#

Calculate spectral difference for 3-parameter snow model.

Computes the Euclidean distance between observed and modeled spectra using a simplified 3-parameter model where shade fills the non-snow fraction.

\[\begin{split}\begin{align} R_{model} & = R_{pure snow}( \phi_{sun}, c_{dust}, s_{grain}) * f_{sca} \\ & + R_{shade} * (1-f_{sca}) \end{align}\end{split}\]

Parameters:

x (array-like) – Model parameters (note: only first 3 are used): - x[0] : f_sca - Fractional snow-covered area (0-1) - x[1] : dust_concentration - Dust concentration in snow (ppm) - x[2] : grain_size - Effective snow grain radius (μm)
spectrum_target (numpy.ndarray) – The observed mixed spectrum to match.
solar_angle (float) – Solar zenith angle of the observation (degrees).
interpolator (spires.interpolator.LutInterpolator) – Callable object that returns modeled snow spectrum given solar_angle, dust_concentration, and grain_size.
shade (numpy.ndarray) – Ideal shade endmember spectrum.

Returns:

Euclidean distance between modeled and target spectra.

Return type:

float

Notes

This 3-parameter model assumes the non-snow fraction is entirely shade (no background component). Use when f_sca is near 100% to avoid overfitting.

Examples

>>> import spires
>>> import numpy as np
>>> interpolator = spires.interpolator.LutInterpolator(lut_file='tests/data/lut_sentinel2b_b2to12_3um_dust.mat')
>>> f_sca = 0.482
>>> dust_concentration = 1000  # ppm
>>> grain_size = 220  # μm
>>> solar_angle = 55.73733298
>>> x = [f_sca, dust_concentration, grain_size]
>>> spectrum_target = np.array([0.3424,0.366,0.3624,0.38932347,0.41624767,0.39567757,0.07043362,0.06267947, 0.3792])
>>> shade = np.array([0,0,0,0,0,0,0,0,0])
>>> spires.snow_diff_3(x=x, spectrum_target=spectrum_target,
...                    solar_angle=solar_angle, interpolator=interpolator, shade=shade)
0.06984199561833446

spires.invert.snow_diff_4(x, spectrum_target, spectrum_background, solar_angle, interpolator, shade)#

Calculate spectral difference for 4-parameter snow model.

Computes the Euclidean distance between observed and modeled spectra using a 4-parameter linear mixing model with snow, shade, and background components.

\[\begin{split}\begin{align} R_{model} & = R_{pure snow}( \phi_{sun}, c_{dust}, s_{grain}) * f_{sca} \\ & + R_{shade} * f_{shade} \\ & + R_{0} * (1 - f_{sca} - f_{shade}) \end{align}\end{split}\]

Parameters:

x (array-like) – Model parameters: - x[0] : f_sca - Fractional snow-covered area (0-1) - x[1] : f_shade - Fractional shaded area (0-1) - x[2] : dust_concentration - Dust concentration in snow (ppm) - x[3] : grain_size - Effective snow grain radius (μm)
spectrum_target (numpy.ndarray) – The observed mixed spectrum to match.
spectrum_background (numpy.ndarray) – The background (snow-free, R_0) spectrum.
solar_angle (float) – Solar zenith angle of the observation (degrees).
interpolator (spires.interpolator.LutInterpolator) – Callable object that returns modeled snow spectrum given solar_angle, dust_concentration, and grain_size.
shade (numpy.ndarray) – Ideal shade endmember spectrum.

Returns:

Euclidean distance between modeled and target spectra.

Return type:

float

Notes

If f_sca is within 2%, consider using 3-parameter solution (snow_diff_3) to avoid overfitting.

Examples

>>> import spires
>>> import numpy as np
>>> interpolator = spires.interpolator.LutInterpolator(lut_file='tests/data/lut_sentinel2b_b2to12_3um_dust.mat')
>>> f_sca = 0.482
>>> f_shade = 0.065
>>> dust_concentration = 1000  # ppm
>>> grain_size = 220  # μm
>>> solar_angle = 55.73733298
>>> x = [f_sca, f_shade, dust_concentration, grain_size]
>>> spectrum_target = np.array([0.3424,0.366,0.3624,0.38932347,0.41624767,0.39567757,0.07043362,0.06267947, 0.3792])
>>> spectrum_background = np.array([0.0182,0.0265,0.0283,0.056067,0.095432,0.12036866,0.12491679,0.07888655,0.1406])
>>> shade = np.array([0,0,0,0,0,0,0,0,0])
>>> diff = spires.snow_diff_4(x=x, spectrum_target=spectrum_target, spectrum_background=spectrum_background,
...                    solar_angle=solar_angle, interpolator=interpolator, shade=shade)
>>> diff
0.08870043573321955

spires.invert.snow_diff_softmax(z, spectrum_target, spectrum_background, solar_angle, interpolator, shade, dust_min, dust_max, grain_min, grain_max)#

Spectral-difference cost in unconstrained (softmax-reparameterized) coordinates.

Maps an unconstrained vector z = [z_snow, z_shade, z_dust, z_grain] to physical parameters such that the simplex constraints (f_sca, f_shade, f_bg ≥ 0, f_sca + f_shade + f_bg = 1) and the box bounds on dust/grain are satisfied by construction. This lets unconstrained solvers (Nelder-Mead, L-BFGS-B, BFGS) replace COBYLA / SLSQP on the fractional sub-problem.

Reparameterization#

Fractions via softmax with z_bg pinned to 0 (gauge fix):

\[(f_{sca}, f_{shade}, f_{bg}) = \mathrm{softmax}(z_{snow}, z_{shade}, 0)\]

Dust and grain via sigmoid-scaled-to-bounds:

\[d = d_{min} + (d_{max} - d_{min})\,\sigma(z_{dust}),\quad g = g_{min} + (g_{max} - g_{min})\,\sigma(z_{grain})\]

The cost itself is the same Euclidean distance as snow_diff_4().

spires.invert.speedy_invert(spectrum_target, spectrum_background, solar_angle, spectrum_shade=None, bands=None, solar_angles=None, dust_concentrations=None, grain_sizes=None, reflectances=None, interpolator=None, lut_dataarray=None, max_eval=100, x0=array([5.0e-01, 5.0e-02, 1.0e+01, 2.5e+02]), algorithm=2)#

Inverts the snow reflectance spectrum using nonlinear optimization.

Parameters:

spectrum_target (numpy.ndarray) – The mixed spectrum to invert. Must be same length as spectrum_background. Must have same band order as spectrum_background and bands.
spectrum_background (numpy.ndarray) – The background (snow-free, R_0) spectrum.
solar_angle (float) – The solar zenith angle of the spectrum target (degrees).
spectrum_shade (numpy.ndarray, optional) – The ideal shaded spectrum. Must be same length as spectrum_target. If None, uses zeros (default: None).
bands (numpy.ndarray, optional) – Band wavelength coordinates of reflectances. Required if interpolator not provided.
solar_angles (numpy.ndarray, optional) – Solar angle coordinates of reflectances. Required if interpolator not provided.
dust_concentrations (numpy.ndarray, optional) – Dust concentration coordinates of reflectances (ppm). Required if interpolator not provided.
grain_sizes (numpy.ndarray, optional) – Grain size coordinates of reflectances (μm). Required if interpolator not provided.
reflectances (numpy.ndarray, optional) – 4D snow reflectance lookup table with dimensions (bands, solar_angles, dust_concentrations, grain_sizes). Required if interpolator not provided.
interpolator (spires.interpolator.LutInterpolator, optional) – Pre-configured interpolator. If provided, overrides individual LUT parameters.
lut_dataarray (xarray.DataArray, optional) – Not currently used. Reserved for future xarray support.
max_eval (int, optional) – Maximum number of optimization iterations. Default is 100.
x0 (array-like, optional) – Initial guess for [fsca, fshade, dust_conc, grain_size]. Default is [0.5, 0.05, 10, 250].
algorithm (int, optional) – Optimization algorithm to use (default: 2). 1 = LN_COBYLA (Constrained Optimization BY Linear Approximations), 2 = LN_NELDERMEAD (Nelder-Mead simplex), 3 = LD_SLSQP (Sequential Least Squares Programming, not working in C++).

Returns:

Optimization results as (fsca, fshade, dust_concentration, grain_size) where:

fsca : float - Fractional snow-covered area (0-1)
fshade : float - Fractional shaded area (0-1)
dust_concentration : float - Dust concentration in snow (ppm)
grain_size : float - Effective snow grain radius (μm)

Return type:

tuple

Examples

>>> import spires
>>> import numpy as np
>>> spectrum_target = np.array([0.3424,0.366,0.3624,0.38932347,0.41624767,0.39567757,0.07043362,0.06267947, 0.3792])
>>> spectrum_background = np.array([0.0182,0.0265,0.0283,0.056067,0.095432,0.12036866,0.12491679,0.07888655,0.1406])
>>> solar_angle = 55.73733298
>>> interpolator = spires.interpolator.LutInterpolator(lut_file='tests/data/lut_sentinel2b_b2to12_3um_dust.mat')
>>> spires.speedy_invert(spectrum_target=spectrum_target, spectrum_background=spectrum_background,
...                      solar_angle=solar_angle, interpolator=interpolator, algorithm=1)
(0.4089303296055291, 0.155201675059351, 138.79357872804923, 364.58404302094834)

spires.invert.speedy_invert_array1d(spectra_targets, spectra_backgrounds, obs_solar_angles, spectrum_shade=None, bands=None, solar_angles=None, dust_concentrations=None, grain_sizes=None, reflectances=None, interpolator=None, lut_dataarray=None, max_eval=100, x0=array([5.0e-01, 5.0e-02, 1.0e+01, 2.5e+02]), algorithm=2)#

Batch inversion of snow reflectance spectra for 1D arrays of observations.

Efficiently processes multiple pixels/observations sequentially using optimized C++ implementations for improved performance.

Parameters:

spectra_targets (numpy.ndarray) – 2D array of mixed spectra to invert with shape (n_observations, n_bands). Must have same length as spectra_backgrounds along first dimension.
spectra_backgrounds (numpy.ndarray) – 2D array of background (snow-free, R_0) spectra with shape (n_observations, n_bands). Must have same length as spectra_targets along first dimension.
obs_solar_angles (numpy.ndarray) – 1D array of solar zenith angles (degrees) for each observation. Must have same length as first dimension of spectra_targets.
spectrum_shade (numpy.ndarray, optional) – 1D array representing the ideal shaded spectrum for all observations. Must have same length as number of bands. If None, uses zeros (default: None).
bands (numpy.ndarray, optional) – Band wavelength coordinates of reflectances. Required if interpolator not provided.
solar_angles (numpy.ndarray, optional) – Solar angle coordinates of reflectances. Required if interpolator not provided.
dust_concentrations (numpy.ndarray, optional) – Dust concentration coordinates of reflectances (ppm). Required if interpolator not provided.
grain_sizes (numpy.ndarray, optional) – Grain size coordinates of reflectances (μm). Required if interpolator not provided.
reflectances (numpy.ndarray, optional) – 4D snow reflectance lookup table with dimensions (bands, solar_angles, dust_concentrations, grain_sizes). Required if interpolator not provided.
interpolator (spires.interpolator.LutInterpolator, optional) – Pre-configured interpolator. If provided, overrides individual LUT parameters.
lut_dataarray (xarray.DataArray, optional) – Not currently used. Reserved for future xarray support.
max_eval (int, optional) – Maximum number of optimization iterations per observation (default: 100).
x0 (array-like, optional) – Initial guess for [fsca, fshade, dust_conc, grain_size]. Default is [0.5, 0.05, 10, 250].
algorithm (int, optional) – Optimization algorithm to use (default: 2). 1 = LN_COBYLA (Constrained Optimization BY Linear Approximations), 2 = LN_NELDERMEAD (Nelder-Mead simplex), 3 = LD_SLSQP (Sequential Least Squares Programming, not working in C++).

Returns:

2D array of shape (n_observations, 4) containing inversion results: - results[:, 0] : Fractional snow-covered area (0-1) - results[:, 1] : Fractional shaded area (0-1) - results[:, 2] : Dust concentration in snow (ppm) - results[:, 3] : Effective snow grain radius (μm)

Return type:

numpy.ndarray

Examples

>>> import spires
>>> import numpy as np
>>> spectra_targets = np.array([[0.3424,0.366,0.3624,0.38932347,0.41624767,0.39567757,0.0704336,0.06267947,0.3792],
...                            [0.2866,0.3046,0.324,0.34468558,0.35373732,0.35651454,0.1807259,0.16601688,0.3488]])
>>> spectra_backgrounds = np.array([[0.0182,0.0265,0.0283,0.0560674,0.0954323,0.1203686,0.1249167,0.0788865,0.1406],
...                                [0.1002,0.1492,0.2088,0.2179780,0.2314920,0.2514020,0.3103066,0.2875081,0.2546]])
>>> obs_solar_angles = np.array([55.73733298, 55.83733298])
>>> interpolator = spires.interpolator.LutInterpolator(lut_file='tests/data/lut_sentinel2b_b2to12_3um_dust.mat')
>>> spires.speedy_invert_array1d(spectra_targets=spectra_targets, spectra_backgrounds=spectra_backgrounds,
...                            obs_solar_angles=obs_solar_angles, interpolator=interpolator, algorithm=1)
array([[4.06627881e-01, 1.45134251e-01, 1.37503982e+02, 3.61158500e+02],
       [2.63873228e-01, 1.83226478e-01, 1.94343159e+02, 3.80170927e+02]])

spires.invert.speedy_invert_array2d(spectra_targets, spectra_backgrounds, obs_solar_angles, max_eval=100, x0=array([5.0e-01, 5.0e-02, 1.0e+01, 2.5e+02]), algorithm=2, bands=None, solar_angles=None, dust_concentrations=None, grain_sizes=None, reflectances=None, interpolator=None)#

Batch inversion of snow reflectance spectra for 2D spatial arrays.

Processes entire images or 2D grids of observations efficiently using optimized C++ implementations. Ideal for processing satellite imagery or gridded data.

Parameters:

spectra_targets (numpy.ndarray) – 3D array of mixed spectra to invert with shape (ny, nx, n_bands): - dim 0: y spatial dimension - dim 1: x spatial dimension - dim 2: spectral bands (must match order of spectra_backgrounds)
spectra_backgrounds (numpy.ndarray) – 3D array of background (snow-free, R_0) spectra with shape (ny, nx, n_bands): - dim 0: y spatial dimension (must match spectra_targets) - dim 1: x spatial dimension (must match spectra_targets) - dim 2: spectral bands (must match order of spectra_targets)
obs_solar_angles (numpy.ndarray) – 2D array of solar zenith angles (degrees) with shape (ny, nx). One angle per spatial location.
max_eval (int, optional) – Maximum number of optimization iterations per pixel (default: 100).
x0 (array-like, optional) – Initial guess for [fsca, fshade, dust_conc, grain_size]. Default is [0.5, 0.05, 10, 250].
algorithm (int, optional) – Optimization algorithm to use (default: 2). 1 = LN_COBYLA (Constrained Optimization BY Linear Approximations), 2 = LN_NELDERMEAD (Nelder-Mead simplex), 3 = LD_SLSQP (Sequential Least Squares Programming, not working in C++).
bands (numpy.ndarray, optional) – Band wavelength coordinates of reflectances. Required if interpolator not provided.
solar_angles (numpy.ndarray, optional) – Solar angle coordinates of reflectances. Required if interpolator not provided.
dust_concentrations (numpy.ndarray, optional) – Dust concentration coordinates of reflectances (ppm). Required if interpolator not provided.
grain_sizes (numpy.ndarray, optional) – Grain size coordinates of reflectances (μm). Required if interpolator not provided.
reflectances (numpy.ndarray, optional) – 4D snow reflectance lookup table with dimensions (bands, solar_angles, dust_concentrations, grain_sizes). Required if interpolator not provided.
interpolator (spires.interpolator.LutInterpolator, optional) – Pre-configured interpolator. If provided, overrides individual LUT parameters.

Returns:

3D array of shape (ny, nx, 4) containing inversion results: - results[:, :, 0] : Fractional snow-covered area (0-1) - results[:, :, 1] : Fractional shaded area (0-1) - results[:, :, 2] : Dust concentration in snow (ppm) - results[:, :, 3] : Effective snow grain radius (μm)

Return type:

numpy.ndarray

Notes

The shade spectrum is automatically set to zeros for all pixels. Future versions may support spatially-varying shade spectra.

spires.invert.speedy_invert_scipy(interpolator: LutInterpolator, spectrum_target, spectrum_background, solar_angle, shade=None, scipy_options=None, mode=3, method='SLSQP')#

Invert snow spectra using scipy.optimize.minimize.

Alternative implementation using SciPy’s optimization routines instead of NLopt. Provides compatibility with legacy code and additional solver options.

Parameters:

interpolator (spires.interpolator.LutInterpolator) – Interpolator object with: - Attributes: bands, solar_angles, dust_concentrations, grain_sizes - Method: interpolate_all(solar_angle, dust_concentration, grain_size)
spectrum_target (numpy.ndarray) – Target spectrum to be inverted. Must be same shape as spectrum_background.
spectrum_background (numpy.ndarray) – Background (snow-free, R_0) spectrum. Must be same shape as spectrum_target.
solar_angle (float) – Solar zenith angle of observation (degrees). Must use same units as interpolator coordinates.
shade (numpy.ndarray, optional) – Ideal shade endmember spectrum. Must be same shape as spectrum_target. If None, uses zeros (default: None).
scipy_options (dict, optional) – SciPy solver options. Default: {‘disp’: False, ‘iprint’: 100, ‘maxiter’: 1000, ‘ftol’: 1e-9}
mode (int, optional) – Number of parameters in model (default: 3). 3 = Simplified model (f_sca, dust, grain_size). 4 = Full model (f_sca, f_shade, dust, grain_size). Use mode=3 when f_sca is near 100% to avoid overfitting.
method (str, optional) – SciPy optimization method (default: ‘SLSQP’). Common options: ‘SLSQP’, ‘L-BFGS-B’, ‘TNC’.

Returns:

(res, model_refl) where:

res : scipy.optimize.OptimizeResult Optimization result object. res.x contains: [f_sca, f_shade, dust_concentration, grain_size]
model_refl : numpy.ndarray The optimized modeled reflectance spectrum.

Return type:

tuple

See also

scipy.optimize.OptimizeResult: Documentation of result object
speedy_invert: NLopt-based implementation (faster)

Examples

>>> import spires
>>> import numpy as np
>>> interpolator = spires.interpolator.LutInterpolator(lut_file='tests/data/lut_sentinel2b_b2to12_3um_dust.mat')
>>> interpolator.make_scipy_interpolator_legacy()
>>> spectrum_target = np.array([0.3424,0.366,0.3624,0.38932347,0.41624767,0.39567757,0.07043362,0.06267947, 0.3792])
>>> spectrum_background = np.array([0.0182,0.0265,0.0283,0.056067,0.095432,0.12036866,0.12491679,0.07888655,0.1406])
>>> solar_angle = 24.0
>>> res, model_refl = spires.speedy_invert_scipy(interpolator=interpolator,
...                                              spectrum_target=spectrum_target,
...                                              spectrum_background=spectrum_background,
...                                              solar_angle=solar_angle,
...                                              mode=3, method='SLSQP')
>>> res.x
array([4.36429085e-01, 5.63570915e-01, 9.91000000e+02, 4.12331162e+01])

spires.invert.speedy_invert_scipy_normalized(interpolator: LutInterpolator, spectrum_target, spectrum_background, solar_angle, spectrum_shade=None, method='COBYLA')#

Invert snow spectra with normalized parameter space.

Performs optimization with all parameters scaled to [0, 1] range to improve convergence for solvers like COBYLA that don’t support parameter-specific step sizes.

Parameters:

interpolator (spires.interpolator.LutInterpolator) – Interpolator object with lookup table and coordinate arrays.
spectrum_target (numpy.ndarray) – Target spectrum to be inverted.
spectrum_background (numpy.ndarray) – Background (snow-free, R_0) spectrum. Must be same shape as spectrum_target.
solar_angle (float) – Solar zenith angle of observation (degrees).
spectrum_shade (numpy.ndarray, optional) – Ideal shade endmember spectrum. Must be same shape as spectrum_target. If None, uses zeros (default: None).
method (str, optional) – SciPy optimization method (default: ‘COBYLA’). COBYLA is recommended as it handles the normalized space well.

Returns:

(res, model_refl) where:

res : scipy.optimize.OptimizeResult Optimization result with res.x containing: [f_sca, f_shade, dust_concentration, grain_size] (dust and grain_size are converted back to physical units)
model_refl : numpy.ndarray The optimized modeled reflectance spectrum.

Return type:

tuple

Notes

This function internally normalizes dust_concentration and grain_size to [0, 1] for optimization, then converts back to physical units. This improves convergence for algorithms that assume similar scales across parameters.

spires.invert.speedy_invert_scipy_softmax(interpolator: LutInterpolator, spectrum_target, spectrum_background, solar_angle, shade=None, scipy_options=None, method='Nelder-Mead', z0=None)#

Unconstrained scipy inversion via softmax reparameterization.

Drops the inequality constraint 1 - f_sca - f_shade ≥ 0 and the box bounds by optimizing in an unconstrained space (see snow_diff_softmax()). Returns the same (res, model_refl) shape as speedy_invert_scipy(), with res.x = [f_sca, f_shade, dust, grain] in physical units.

Parameters:

method (str, optional) – Any unconstrained scipy method (‘Nelder-Mead’, ‘BFGS’, ‘L-BFGS-B’, ‘Powell’). Default ‘Nelder-Mead’.
z0 (array-like, optional) – Initial guess in z-space (length 4). If None, defaults to zeros, which corresponds to f = (1/3, 1/3, 1/3) and dust/grain at the bounds midpoint.

spires.invert.speedy_invert_xarray(spectra_targets, spectra_backgrounds, obs_solar_angles, lut_dataarray, spectrum_shade=None, max_eval=100, x0=array([5.0e-01, 5.0e-02, 1.0e+01, 2.5e+02]), algorithm=2)#

Batch inversion of snow reflectance spectra using xarray DataArrays.

Provides a high-level interface for processing geospatial data with coordinate information preserved. Automatically handles dimension ordering and metadata.

Parameters:

spectra_targets (xarray.DataArray) – Mixed spectra to invert with dimensions (y, x, band). Will be automatically transposed if needed.
spectra_backgrounds (xarray.DataArray) – Background (snow-free, R_0) spectra with dimensions (y, x, band). Must have same spatial dimensions as spectra_targets.
obs_solar_angles (xarray.DataArray) – Solar zenith angles (degrees) with dimensions (y, x). One angle per spatial location.
lut_dataarray (xarray.DataArray) – Lookup table with dimensions (band, solar_angle, dust_concentration, grain_size). Coordinates are extracted and used for interpolation.
spectrum_shade (numpy.ndarray, optional) – 1D array representing the ideal shaded spectrum. Must have same length as number of bands. If None, uses zeros (default: None).
max_eval (int, optional) – Maximum number of optimization iterations per pixel (default: 100).
x0 (array-like, optional) – Initial guess for [fsca, fshade, dust_conc, grain_size]. Default is [0.5, 0.05, 10, 250].
algorithm (int, optional) – Optimization algorithm to use (default: 2). 1 = LN_COBYLA (Constrained Optimization BY Linear Approximations), 2 = LN_NELDERMEAD (Nelder-Mead simplex), 3 = LD_SLSQP (Sequential Least Squares Programming, not working in C++).

Returns:

Return type:

numpy.ndarray

Notes

Currently returns a numpy array. Future versions will return an xarray.DataArray with appropriate coordinates and metadata (see TODO comment in code).

Dask-parallel inversion of snow reflectance spectra.

spires.parallel.encode_results(ds, fill_value=-1, fsca_scale=100, fshade_scale=100, fraction_dtype=<class 'numpy.int8'>, concentration_dtype=<class 'numpy.int16'>)#

Encode an inversion result Dataset for compact storage.

Replaces NaNs with fill_value, scales the fractional variables, and casts to integer dtypes. The original Dataset is unchanged; a new one is returned with _FillValue, scale_factor, and add_offset attrs set so CF-aware readers (xarray, GDAL) decode back to physical units.

Parameters:

ds (xarray.Dataset) – Output of speedy_invert_dask with variables fsca, fshade, dust_concentration, grain_size in physical units (NaN for nodata).
fill_value (int, optional) – Sentinel for NaN pixels (default: -1).
fsca_scale (float, optional) – Multiplier applied before integer cast (default: 100, i.e. percent).
fshade_scale (float, optional) – Multiplier applied before integer cast (default: 100, i.e. percent).
fraction_dtype (numpy dtype, optional) – Integer type for fsca/fshade (default: np.int8).
concentration_dtype (numpy dtype, optional) – Integer type for dust_concentration/grain_size (default: np.int16).

Returns:

Dataset with the same variable names, encoded for storage.

Return type:

xarray.Dataset

spires.parallel.speedy_invert_dask(spectra_targets, spectra_backgrounds, obs_solar_angles, interpolator, spectrum_shade=None, max_eval=100, x0=array([5.0e-01, 5.0e-02, 1.0e+01, 2.5e+02]), algorithm=2, client=None, scatter_lut=True)#

Parallel inversion of snow reflectance spectra using Dask and xarray.

This method enables distributed processing of large satellite imagery datasets by leveraging Dask’s parallel computation capabilities through xarray.apply_ufunc. It’s particularly useful for processing time series of satellite imagery where the data is too large to fit in memory.

Parameters:

spectra_targets (xarray.DataArray) – Mixed spectra to invert. Must have a ‘band’ dimension and can have any combination of spatial (x, y) and temporal (time) dimensions. Shape: (time, y, x, band) or (y, x, band).
spectra_backgrounds (xarray.DataArray) – Background (snow-free, R_0) spectra with same dimensions as targets except potentially missing the time dimension if using static backgrounds. Shape: (y, x, band).
obs_solar_angles (xarray.DataArray) – Solar zenith angles (degrees) for each observation. Shape: (time, y, x) or (y, x).
interpolator (spires.interpolator.LutInterpolator) – Lookup table interpolator object.
spectrum_shade (numpy.ndarray, optional) – Currently unused: speedy_invert_array2d hardcodes a zero shade spectrum. Accepted for forward compatibility.
max_eval (int, optional) – Maximum optimization iterations per pixel (default: 100).
x0 (array-like, optional) – Initial guess: [fsca, fshade, dust_conc (ppm), grain_size (μm)].
algorithm (int, optional) – NLopt algorithm code (see speedy_invert_array2d).
client (dask.distributed.Client, optional) – Dask client. If None, uses the default client when one is active.
scatter_lut (bool, optional) – Broadcast the LUT to all workers (default: True).

Notes

The C++ inversion releases the Python GIL, so Dask clients with threads_per_worker > 1 get real parallel speedup and share a single LUT copy per worker process — usually preferable to many single-threaded workers when memory or IPC is a concern.

Returns:: Dataset with variables fsca, fshade, dust_concentration, grain_size, preserving input coordinates and dimensions.
Return type:: xarray.Dataset