From 421ae828fd7d4e5a89cfa092575eecfa74427ecf Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Thu, 20 Nov 2025 22:02:11 +0000 Subject: [PATCH 01/13] WIP fix spectrum object --- .pre-commit-config.yaml | 2 +- sunkit_spex/conftest.py | 44 ---- sunkit_spex/spectrum/spectrum.py | 213 +++++++++++++++++--- sunkit_spex/spectrum/tests/test_spectrum.py | 122 ++++++++++- 4 files changed, 299 insertions(+), 82 deletions(-) delete mode 100644 sunkit_spex/conftest.py diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 8e0601c7..64489eee 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -3,7 +3,7 @@ repos: - repo: https://github.com/astral-sh/ruff-pre-commit rev: "v0.14.10" hooks: - - id: ruff + - id: ruff-check args: ["--fix"] - id: ruff-format - repo: https://github.com/PyCQA/isort diff --git a/sunkit_spex/conftest.py b/sunkit_spex/conftest.py deleted file mode 100644 index 44d984ab..00000000 --- a/sunkit_spex/conftest.py +++ /dev/null @@ -1,44 +0,0 @@ -# This file is used to configure the behavior of pytest when using the Astropy -# test infrastructure. - - -# Uncomment the following line to treat all DeprecationWarnings as -# exceptions. For Astropy v2.0 or later, there are 2 additional keywords, -# as follow (although default should work for most cases). -# To ignore some packages that produce deprecation warnings on import -# (in addition to 'compiler', 'scipy', 'pygments', 'ipykernel', and -# 'setuptools'), add: -# modules_to_ignore_on_import=['module_1', 'module_2'] -# To ignore some specific deprecation warning messages for Python version -# MAJOR.MINOR or later, add: -# warnings_to_ignore_by_pyver={(MAJOR, MINOR): ['Message to ignore']} -# enable_deprecations_as_exceptions() - -# Uncomment and customize the following lines to add/remove entries from -# the list of packages for which version numbers are displayed when running -# the tests. Making it pass for KeyError is essential in some cases when -# the package uses other astropy affiliated packages. -# try: -# PYTEST_HEADER_MODULES['Astropy'] = 'astropy' -# PYTEST_HEADER_MODULES['scikit-image'] = 'skimage' -# del PYTEST_HEADER_MODULES['h5py'] -# except (NameError, KeyError): # NameError is needed to support Astropy < 1.0 -# pass - -# Uncomment the following lines to display the version number of the -# package rather than the version number of Astropy in the top line when -# running the tests. -# import os -# -# This is to figure out the package version, rather than -# using Astropy's -# try: -# from .version import version -# except ImportError: -# version = 'dev' -# -# try: -# packagename = os.path.basename(os.path.dirname(__file__)) -# TESTED_VERSIONS[packagename] = version -# except NameError: # Needed to support Astropy <= 1.0.0 -# pass diff --git a/sunkit_spex/spectrum/spectrum.py b/sunkit_spex/spectrum/spectrum.py index b1ca2d86..d86fdd67 100644 --- a/sunkit_spex/spectrum/spectrum.py +++ b/sunkit_spex/spectrum/spectrum.py @@ -117,6 +117,10 @@ def bin_edges(self): return self._bin_edges return self._edges_from_centers(self.value, self.unit) + def __array_finalize__(self, obj): + super().__array_finalize__(obj) + self._bin_edges = getattr(obj, "_bin_edges", None) + class Spectrum(NDCube): r""" @@ -134,15 +138,15 @@ class Spectrum(NDCube): ---------- data : `~astropy.units.Quantity` The data for this spectrum. This can be a simple `~astropy.units.Quantity`, - or an existing `~Spectrum1D` or `~ndcube.NDCube` object. + or an existing `~Spectrum` or `~ndcube.NDCube` object. uncertainty : `~astropy.nddata.NDUncertainty` Contains uncertainty information along with propagation rules for spectrum arithmetic. Can take a unit, but if none is given, will use the unit defined in the flux. spectral_axis : `~astropy.units.Quantity` or `~specutils.SpectralAxis` - Dispersion information with the same shape as the dimension specified by spectral_dimension + Dispersion information with the same shape as the dimension specified by spectral_axis_index of shape plus one if specifying bin edges. - spectral_dimension : `int` default 0 + spectral_axis_index : `int` default 0 The dimension of the data which represents the spectral information default to first dimension index 0. mask : `~numpy.ndarray`-like Array where values in the flux to be masked are those that @@ -157,7 +161,7 @@ class Spectrum(NDCube): >>> import numpy as np >>> import astropy.units as u >>> from sunkit_spex.spectrum import Spectrum - >>> spec = Spectrum(np.arange(1, 11)*u.watt, spectral_axis=np.arange(1, 12)*u.keV) + >>> spec = Spectrum(np.arange(1, 11)*u.watt,spectral_axis=np.arange(1, 12)*u.keV) >>> spec 0: + raise ValueError( + "Initializer contains unknown arguments(s): {}.".format(", ".join(map(str, unknown_kwargs))) + ) + + # Handle initializing from NDCube objects + if isinstance(data, NDCube): + if data.unit is None: + raise ValueError("Input NDCube missing unit parameter") + + # Change the data array from bare ndarray to a Quantity + q_data = data.data << u.Unit(data.unit) + + self.__init__(q_data, uncertainty=data.uncertainty, mask=data.mask, wcs=data.wcs) + return + + self._spectral_axis_index = spectral_axis_index + # If here data is should be an array or quantity + if spectral_axis_index is None and data is not None: + if data.ndim == 1: + self._spectral_axis_index = 0 + elif data is None: + self._spectral_axis_index = 0 + + # Ensure that the data argument is an astropy quantity + # if data is not None: + # if not isinstance(data, u.Quantity): + # raise ValueError("Data must be a `Quantity` object.") + # if data.isscalar: + # data = u.Quantity([data]) # Ensure that the unit information codified in the quantity object is # the One True Unit. @@ -191,41 +235,146 @@ def __init__( # If flux and spectral axis are both specified, check that their lengths # match or are off by one (implying the spectral axis stores bin edges) if data is not None and spectral_axis is not None: - if spectral_axis.shape[0] == data.shape[spectral_dimension]: + if spectral_axis.shape[0] == data.shape[self.spectral_axis_index]: bin_specification = "centers" - elif spectral_axis.shape[0] == data.shape[spectral_dimension] + 1: + elif spectral_axis.shape[0] == data.shape[self.spectral_axis_index] + 1: bin_specification = "edges" else: raise ValueError( - f"Spectral axis length ({spectral_axis.shape[0]}) must be the same size or one " - "greater (if specifying bin edges) than that of the spextral" - f"axis ({data.shape[spectral_dimension]})" + f"Spectral axis length ({spectral_axis.shape[0]}) must be the " + "same size or one greater (if specifying bin edges) than that " + f"of the corresponding flux axis ({data.shape[self.spectral_axis_index]})" ) + # If a WCS is provided, determine which axis is the spectral axis + if wcs is not None: + naxis = None + if hasattr(wcs, "naxis"): + naxis = wcs.naxis + # GWCS doesn't have naxis + elif hasattr(wcs, "world_n_dim"): + naxis = wcs.world_n_dim + + if naxis is not None and naxis > 1: + temp_axes = [] + phys_axes = wcs.world_axis_physical_types + if self._spectral_axis_index is None: + for i in range(len(phys_axes)): + if phys_axes[i] is None: + continue + if phys_axes[i][0:2] == "em" or phys_axes[i][0:5] == "spect" or phys_axes[i][7:12] == "Spect": + temp_axes.append(i) + if len(temp_axes) != 1: + raise ValueError( + f"Input WCS must have exactly one axis with spectral units, found {len(temp_axes)}" + ) + # Due to FITS conventions, the WCS axes are listed in opposite + # order compared to the data array. + self._spectral_axis_index = len(data.shape) - temp_axes[0] - 1 + + else: + if data is not None and data.ndim == 1: + self._spectral_axis_index = 0 + else: + if self.spectral_axis_index is None: + raise ValueError("WCS is 1D but flux is multi-dimensional. Please specify spectral_axis_index.") + # Attempt to parse the spectral axis. If none is given, try instead to # parse a given wcs. This is put into a GWCS object to # then be used behind-the-scenes for all operations. if spectral_axis is not None: # Ensure that the spectral axis is an astropy Quantity - if not isinstance(spectral_axis, u.Quantity): + if not isinstance(spectral_axis, (u.Quantity, SpectralAxis)): raise ValueError("Spectral axis must be a `Quantity` or `SpectralAxis` object.") - # If a spectral axis is provided as an astropy Quantity, convert it - # to a SpectralAxis object. + # If spectral axis is provided as an astropy Quantity, convert it + # to a specutils SpectralAxis object. if not isinstance(spectral_axis, SpectralAxis): - if spectral_axis.shape[0] == data.shape[spectral_dimension] + 1: - bin_specification = "edges" - else: - bin_specification = "centers" self._spectral_axis = SpectralAxis(spectral_axis, bin_specification=bin_specification) + # If a SpectralAxis object is provided, we assume it doesn't need + # information from other keywords added + else: + self._spectral_axis = spectral_axis - wcs = gwcs_from_array(self._spectral_axis) + if wcs is None: + wcs = gwcs_from_array(self._spectral_axis) - super().__init__( - data=data.value if isinstance(data, u.Quantity) else data, - wcs=wcs, - mask=mask, - meta=meta, - uncertainty=uncertainty, - **kwargs, + elif wcs is None: + # If no spectral axis or wcs information is provided, initialize + # with an empty gwcs based on the data. + if self.spectral_axis_index is None: + if data.ndim == 1: + self._spectral_axis_index = 0 + else: + raise ValueError("Must specify spectral_axis_index if no WCS or spectral axis is input.") + size = data.shape[self.spectral_axis_index] if not data.isscalar else 1 + wcs = gwcs_from_array( + np.arange(size) * u.Unit("pixel"), data.shape, spectral_axis_index=self.spectral_axis_index ) + + super().__init__(data=data.value if isinstance(data, u.Quantity) else data, wcs=wcs, **kwargs) + + # If no spectral_axis was provided, create a SpectralCoord based on + # the WCS + if spectral_axis is None: + # If the WCS doesn't have a spectral attribute, we assume it's the + # dummy GWCS we created or a solely spectral WCS + if hasattr(self.wcs, "spectral"): + # Handle generated 1D WCS that aren't set to spectral + if not self.wcs.is_spectral and self.wcs.naxis == 1: + spec_axis = self.wcs.pixel_to_world(np.arange(self.data.shape[self.spectral_axis_index])) + else: + spec_axis = self.wcs.spectral.pixel_to_world(np.arange(self.data.shape[self.spectral_axis_index])) + else: + # We now keep the entire GWCS, including spatial information, so we need to include + # all axes in the pixel_to_world call. Note that this assumes/requires that the + # dispersion is the same at all spatial locations. + wcs_args = [] + for i in range(len(self.data.shape)): + wcs_args.append(np.zeros(self.data.shape[self.spectral_axis_index])) + # Replace with arange for the spectral axis + wcs_args[self.spectral_axis_index] = np.arange(self.data.shape[self.spectral_axis_index]) + wcs_args.reverse() + temp_coords = self.wcs.pixel_to_world(*wcs_args) + # If there are spatial axes, temp_coords will have a SkyCoord and a SpectralCoord + if isinstance(temp_coords, list): + for coords in temp_coords: + if isinstance(coords, SpectralCoord): + spec_axis = coords + break + else: + # WCS axis ordering is reverse of numpy + spec_axis = temp_coords[len(temp_coords) - self.spectral_axis_index - 1] + else: + spec_axis = temp_coords + + try: + if spec_axis.unit.is_equivalent(u.one): + spec_axis = spec_axis * u.pixel + except AttributeError: + raise AttributeError(f"spec_axis does not have unit: {type(spec_axis)} {spec_axis}") + + self._spectral_axis = SpectralAxis(spec_axis) + + # make sure that spectral axis is strictly increasing + if not np.all(self._spectral_axis[1:] >= self._spectral_axis[:-1]): + raise ValueError("Spectral axis must be strictly increasing or decreasing.") + + if hasattr(self, "uncertainty") and self.uncertainty is not None: + if not data.shape == self.uncertainty.array.shape: + raise ValueError( + f"Data axis ({data.shape}) and uncertainty ({self.uncertainty.array.shape}) shapes must be the " + "same." + ) + + # @property + # def data(self): + # return u.Quantity(self._data, unit=self.unit, copy=False) + + @property + def spectral_axis(self): + return self._spectral_axis + + @property + def spectral_axis_index(self): + return self._spectral_axis_index diff --git a/sunkit_spex/spectrum/tests/test_spectrum.py b/sunkit_spex/spectrum/tests/test_spectrum.py index 1763bb05..3e967aa0 100644 --- a/sunkit_spex/spectrum/tests/test_spectrum.py +++ b/sunkit_spex/spectrum/tests/test_spectrum.py @@ -1,16 +1,128 @@ +from operator import add, mul, sub, truediv + import numpy as np +import pytest +from ndcube import NDCube +from ndcube.extra_coords import QuantityTableCoordinate, TimeTableCoordinate +from ndcube.wcs.wrappers import CompoundLowLevelWCS from numpy.testing import assert_array_equal import astropy.units as u +from astropy.time import Time +from astropy.wcs import WCS -from sunkit_spex.spectrum.spectrum import Spectrum +from sunkit_spex.spectrum.spectrum import SpectralAxis, Spectrum -def test_spectrum_bin_edges(): +def test_spectrum_quantity_bin_edges(): spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=np.arange(1, 12) * u.keV) assert_array_equal(spec._spectral_axis, [1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5] * u.keV) -def test_spectrum_bin_centers(): - spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=(np.arange(1, 11) - 0.5) * u.keV) - assert_array_equal(spec._spectral_axis, [0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5] * u.keV) +def test_spectrum_quantity_bin_centers(): + spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=(np.arange(1, 11) + 0.5) * u.keV) + assert_array_equal(spec._spectral_axis, [1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5] * u.keV) + + +def test_spectrum_spectral_axis_bin_edges(): + spec_axis = SpectralAxis(np.arange(1, 12) * u.keV, bin_specification="edges") + spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=spec_axis) + assert_array_equal(spec._spectral_axis, [1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5] * u.keV) + + +def test_spectrum_spectral_axis_bin_centers(): + spec_axis = SpectralAxis((np.arange(1, 11) + 0.5) * u.keV, bin_specification="centers") + spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=spec_axis) + assert_array_equal(spec._spectral_axis, [1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5] * u.keV) + + +def test_spectrum_from_ndcube_wcs(): + header = { + "CTYPE1": "TIME ", + "CUNIT1": "s", + "CDELT1": 10, + "CRPIX1": 0, + "CRVAL1": 0, + "CTYPE2": "ENER ", + "CUNIT2": "keV", + "CDELT2": 1, + "CRPIX2": 0, + "CRVAL2": 1, + "DATEREF": "2020-01-01T00:00:00", + } + wcs = WCS(header=header) + shape = (10, 5) + wcs.array_shape = shape + data = np.arange(np.prod(shape), dtype=int).reshape(shape) + cube = NDCube(data.T, wcs=wcs) + with pytest.raises(ValueError, match="Input NDCube missing unit.*"): + Spectrum(cube) + + cube = NDCube(data, wcs=wcs, unit=u.ph) + spec = Spectrum(cube) + assert isinstance(spec, Spectrum) + + +def test_spectrum_from_cube_wcs_tab(): + energy = np.arange(1, 11) * u.keV + energy_coord = QuantityTableCoordinate(energy, names="energy", physical_types="em.energy") + data = np.random.rand(len(energy)) + cube = NDCube(data=data, wcs=energy_coord.wcs, unit=u.ph) + spec = Spectrum(cube) + assert False + + +def test_spectrum_spectra_axis_detection(): + energy = np.arange(1, 11) * u.keV + energy_coord = QuantityTableCoordinate(energy, names="energy", physical_types="em.energy") + times = Time("2020-01-01") + 5 * np.arange(0, 11) * u.s + time_coord = TimeTableCoordinate(times, names="time", physical_types="time") + time_energy_wcs = (time_coord & energy_coord).wcs + data = np.arange(5 * 10).reshape(5, 10) + spec1 = Spectrum(data * u.ph, wcs=time_energy_wcs) + + energy_energy_wcs = (energy_coord & time_coord).wcs + data = np.arange(5 * 10).reshape(10, 5) + spec2 = Spectrum(data * u.ph, wcs=energy_energy_wcs) + + assert False + + +def test_spectrum_from_cubs_wcs_norm_tab(): + header = { + "CTYPE1": "TIME ", + "CUNIT1": "s", + "CDELT1": 10, + "CRPIX1": 0, + "CRVAL1": 0, + "DATEREF": "2020-01-01T00:00:00", + } + time_wcs = WCS(header=header) + energy = np.arange(10) * u.keV + energy_coord = QuantityTableCoordinate(energy, names="energy", physical_types="em.energy") + comp_wcs = CompoundLowLevelWCS(time_wcs, energy_coord.wcs) + cube = NDCube(np.arange(10 * 5).reshape(10, 5), unit=u.ph, wcs=comp_wcs) + spec = Spectrum(cube) + assert False + + +def test_slice(): + spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=(np.arange(1, 11) + 0.5) * u.keV) + sliced_spec = spec[5:] + assert sliced_spec.shape == (5,) + assert sliced_spec.spectral_axis.shape == (5,) + + +@pytest.mark.parametrize( + "op, value, res", + [ + (add, 2 * u.W, np.arange(1.0, 11) + 2), + (sub, 2 * u.W, np.arange(1.0, 11) - 2), + (mul, 2, np.arange(1.0, 11) * 2), + (truediv, 2 * u.W, np.arange(1.0, 11) / 2), + ], +) +def test_arithmetic_operators(op, value, res): + spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=(np.arange(1, 11) + 0.5) * u.keV) + res_spec = op(spec, value) + assert_array_equal(res_spec.data, res) From c944274905348e2b921e76d4eb3601e858ef9d6a Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Mon, 1 Dec 2025 17:00:30 +0000 Subject: [PATCH 02/13] WIP creation from various arrays, wcs, cube combos --- sunkit_spex/spectrum/conftest.py | 0 sunkit_spex/spectrum/spectrum.py | 111 +++++++++++--------- sunkit_spex/spectrum/tests/test_spectrum.py | 62 +++++++---- 3 files changed, 106 insertions(+), 67 deletions(-) create mode 100644 sunkit_spex/spectrum/conftest.py diff --git a/sunkit_spex/spectrum/conftest.py b/sunkit_spex/spectrum/conftest.py new file mode 100644 index 00000000..e69de29b diff --git a/sunkit_spex/spectrum/spectrum.py b/sunkit_spex/spectrum/spectrum.py index d86fdd67..99276501 100644 --- a/sunkit_spex/spectrum/spectrum.py +++ b/sunkit_spex/spectrum/spectrum.py @@ -7,9 +7,9 @@ from astropy.coordinates import SpectralCoord from astropy.modeling.tabular import Tabular1D from astropy.utils import lazyproperty +from astropy.wcs.wcsapi import sanitize_slices __all__ = ["SpectralAxis", "Spectrum", "gwcs_from_array"] -__doctest_requires__ = {"Spectrum": ["ndcube>=2.3"]} __doctest_requires__ = {"Spectrum": ["ndcube>=2.3"]} @@ -127,8 +127,8 @@ class Spectrum(NDCube): Spectrum container for data with one spectral axis. Note that "1D" in this case refers to the fact that there is only one - spectral axis. `Spectrum` can contain "vector 1D spectra" by having the - ``flux`` have a shape with dimension greater than 1. + spectral axis. `Spectrum` can contain ND data where + ``data`` have a shape with dimension greater than 1. Notes ----- @@ -176,12 +176,9 @@ def __init__( self, data, *, - uncertainty=None, spectral_axis=None, spectral_axis_index=None, wcs=None, - mask=None, - meta=None, **kwargs, ): # If the data argument is already a Spectrum (as it would @@ -207,10 +204,15 @@ def __init__( if data.unit is None: raise ValueError("Input NDCube missing unit parameter") + if spectral_axis is None: + raise ValueError("Spectral axis must be specified") + # Change the data array from bare ndarray to a Quantity q_data = data.data << u.Unit(data.unit) - self.__init__(q_data, uncertainty=data.uncertainty, mask=data.mask, wcs=data.wcs) + self.__init__( + q_data, wcs=data.wcs, mask=data.mask, uncertainty=data.uncertainty, spectral_axis=spectral_axis + ) return self._spectral_axis_index = spectral_axis_index @@ -221,33 +223,15 @@ def __init__( elif data is None: self._spectral_axis_index = 0 - # Ensure that the data argument is an astropy quantity - # if data is not None: - # if not isinstance(data, u.Quantity): - # raise ValueError("Data must be a `Quantity` object.") - # if data.isscalar: - # data = u.Quantity([data]) - # Ensure that the unit information codified in the quantity object is # the One True Unit. kwargs.setdefault("unit", data.unit if isinstance(data, u.Quantity) else kwargs.get("unit")) - # If flux and spectral axis are both specified, check that their lengths - # match or are off by one (implying the spectral axis stores bin edges) - if data is not None and spectral_axis is not None: - if spectral_axis.shape[0] == data.shape[self.spectral_axis_index]: - bin_specification = "centers" - elif spectral_axis.shape[0] == data.shape[self.spectral_axis_index] + 1: - bin_specification = "edges" - else: - raise ValueError( - f"Spectral axis length ({spectral_axis.shape[0]}) must be the " - "same size or one greater (if specifying bin edges) than that " - f"of the corresponding flux axis ({data.shape[self.spectral_axis_index]})" - ) - # If a WCS is provided, determine which axis is the spectral axis if wcs is not None: + if spectral_axis is None: + raise ValueError("Spectral axis must be specified") + naxis = None if hasattr(wcs, "naxis"): naxis = wcs.naxis @@ -279,27 +263,56 @@ def __init__( if self.spectral_axis_index is None: raise ValueError("WCS is 1D but flux is multi-dimensional. Please specify spectral_axis_index.") + # If data and spectral axis are both specified, check that their lengths + # match or are off by one (implying the spectral axis stores bin edges) + if data is not None and spectral_axis is not None: + if spectral_axis.shape[0] == data.shape[self.spectral_axis_index]: + bin_specification = "centers" + elif spectral_axis.shape[0] == data.shape[self.spectral_axis_index] + 1: + bin_specification = "edges" + else: + raise ValueError( + f"Spectral axis length ({spectral_axis.shape[0]}) must be the " + "same size or one greater (if specifying bin edges) than that " + f"of the corresponding data axis ({data.shape[self.spectral_axis_index]})" + ) + # Attempt to parse the spectral axis. If none is given, try instead to # parse a given wcs. This is put into a GWCS object to # then be used behind-the-scenes for all operations. - if spectral_axis is not None: - # Ensure that the spectral axis is an astropy Quantity - if not isinstance(spectral_axis, (u.Quantity, SpectralAxis)): - raise ValueError("Spectral axis must be a `Quantity` or `SpectralAxis` object.") - - # If spectral axis is provided as an astropy Quantity, convert it - # to a specutils SpectralAxis object. - if not isinstance(spectral_axis, SpectralAxis): - self._spectral_axis = SpectralAxis(spectral_axis, bin_specification=bin_specification) - # If a SpectralAxis object is provided, we assume it doesn't need - # information from other keywords added + + # Ensure that the spectral axis is an astropy Quantity or SpectralAxis + if not isinstance(spectral_axis, (u.Quantity, SpectralAxis)): + raise ValueError("Spectral axis must be a `Quantity` or `SpectralAxis` object.") + + # If spectral axis is provided as an astropy Quantity, convert it + # to a specutils SpectralAxis object. + if not isinstance(spectral_axis, SpectralAxis): + self._spectral_axis = SpectralAxis(spectral_axis, bin_specification=bin_specification) + # If a SpectralAxis object is provided, we assume it doesn't need + # information from other keywords added + else: + self._spectral_axis = spectral_axis + + # Check the spectral_axis matches the wcs + if wcs is not None: + if hasattr(wcs, "spectral"): + wcs_coords = wcs.spectral.pixel_to_world(np.arange(data.shape[self.spectral_axis_index])).to("keV") + elif wcs.pixel_n_dim == 1: + wcs_coords = wcs.pixel_to_world(np.arange(data.shape[self.spectral_axis_index])) else: - self._spectral_axis = spectral_axis + array_index = wcs.pixel_n_dim - self._spectral_axis_index - 1 + pixels = [0] * wcs.pixel_n_dim + pixels[array_index] = np.arange(data.shape[self.spectral_axis_index]) + wcs_coords = wcs.pixel_to_world(*pixels)[array_index] + + if not u.allclose(self._spectral_axis, wcs_coords): + raise ValueError(f"Spectral axis {self._spectral_axis} and wcs spectral axis {wcs_coords} must match.") - if wcs is None: - wcs = gwcs_from_array(self._spectral_axis) + if wcs is None: + wcs = gwcs_from_array(self._spectral_axis) - elif wcs is None: + if wcs is None: # If no spectral axis or wcs information is provided, initialize # with an empty gwcs based on the data. if self.spectral_axis_index is None: @@ -329,9 +342,7 @@ def __init__( # We now keep the entire GWCS, including spatial information, so we need to include # all axes in the pixel_to_world call. Note that this assumes/requires that the # dispersion is the same at all spatial locations. - wcs_args = [] - for i in range(len(self.data.shape)): - wcs_args.append(np.zeros(self.data.shape[self.spectral_axis_index])) + wcs_args = [np.zeros(self.data.shape[self.spectral_axis_index]) for i in range(len(self.data.shape))] # Replace with arange for the spectral axis wcs_args[self.spectral_axis_index] = np.arange(self.data.shape[self.spectral_axis_index]) wcs_args.reverse() @@ -367,9 +378,11 @@ def __init__( "same." ) - # @property - # def data(self): - # return u.Quantity(self._data, unit=self.unit, copy=False) + def __getitem__(self, item): + sliced_cube = super().__getitem__(item) + item = tuple(sanitize_slices(item, len(self.shape))) + sliced_spec_axis = self.spectral_axis[item[self.spectral_axis_index]] + return Spectrum(sliced_cube, spectral_axis=sliced_spec_axis) @property def spectral_axis(self): diff --git a/sunkit_spex/spectrum/tests/test_spectrum.py b/sunkit_spex/spectrum/tests/test_spectrum.py index 3e967aa0..6f462986 100644 --- a/sunkit_spex/spectrum/tests/test_spectrum.py +++ b/sunkit_spex/spectrum/tests/test_spectrum.py @@ -8,11 +8,14 @@ from numpy.testing import assert_array_equal import astropy.units as u +from astropy.tests.helper import assert_quantity_allclose from astropy.time import Time from astropy.wcs import WCS from sunkit_spex.spectrum.spectrum import SpectralAxis, Spectrum +rng = np.random.default_rng() + def test_spectrum_quantity_bin_edges(): spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=np.arange(1, 12) * u.keV) @@ -46,46 +49,67 @@ def test_spectrum_from_ndcube_wcs(): "CTYPE2": "ENER ", "CUNIT2": "keV", "CDELT2": 1, - "CRPIX2": 0, - "CRVAL2": 1, + "CRPIX2": 0.5, + "CRVAL2": 0.0, "DATEREF": "2020-01-01T00:00:00", } wcs = WCS(header=header) shape = (10, 5) wcs.array_shape = shape + + spec_axis_edges = np.arange(11) * u.keV + spec_axis_centers = spec_axis_edges[:-1] + np.diff(spec_axis_edges) * 0.5 + data = np.arange(np.prod(shape), dtype=int).reshape(shape) cube = NDCube(data.T, wcs=wcs) with pytest.raises(ValueError, match="Input NDCube missing unit.*"): Spectrum(cube) cube = NDCube(data, wcs=wcs, unit=u.ph) - spec = Spectrum(cube) + + with pytest.raises(ValueError, match="Spectral axis must be specified"): + Spectrum(cube) + + with pytest.raises(ValueError, match=r"Spectral axis"): + Spectrum(cube, spectral_axis=spec_axis_edges[:-1]) + + spec = Spectrum(cube, spectral_axis=spec_axis_edges) assert isinstance(spec, Spectrum) + assert spec.spectral_axis_index == 0 + assert spec.shape == (10, 5) + assert_quantity_allclose(spec_axis_centers, spec.wcs.pixel_to_world(0, np.arange(10))[1].to("keV")) + assert_quantity_allclose(spec_axis_edges, spec.spectral_axis.bin_edges) def test_spectrum_from_cube_wcs_tab(): - energy = np.arange(1, 11) * u.keV - energy_coord = QuantityTableCoordinate(energy, names="energy", physical_types="em.energy") - data = np.random.rand(len(energy)) + spec_axis_edges = np.arange(11) * u.keV + spec_axis_centers = spec_axis_edges[:-1] + np.diff(spec_axis_edges) * 0.5 + energy_coord = QuantityTableCoordinate(spec_axis_centers, names="energy", physical_types="em.energy") + data = rng.random(len(spec_axis_centers)) cube = NDCube(data=data, wcs=energy_coord.wcs, unit=u.ph) - spec = Spectrum(cube) - assert False + + spec = Spectrum(cube, spectral_axis=spec_axis_edges) + assert isinstance(spec, Spectrum) + + assert spec.spectral_axis_index == 0 + assert spec.shape == (10,) + assert_quantity_allclose(spec_axis_centers, spec.wcs.pixel_to_world(np.arange(10)).to("keV")) def test_spectrum_spectra_axis_detection(): - energy = np.arange(1, 11) * u.keV + energy = (np.arange(0, 10) + 0.5) * u.keV energy_coord = QuantityTableCoordinate(energy, names="energy", physical_types="em.energy") times = Time("2020-01-01") + 5 * np.arange(0, 11) * u.s time_coord = TimeTableCoordinate(times, names="time", physical_types="time") time_energy_wcs = (time_coord & energy_coord).wcs - data = np.arange(5 * 10).reshape(5, 10) - spec1 = Spectrum(data * u.ph, wcs=time_energy_wcs) - - energy_energy_wcs = (energy_coord & time_coord).wcs data = np.arange(5 * 10).reshape(10, 5) - spec2 = Spectrum(data * u.ph, wcs=energy_energy_wcs) + spec1 = Spectrum(data * u.ph, wcs=time_energy_wcs, spectral_axis=np.arange(11) * u.keV) + assert spec1.spectral_axis_index == 0 - assert False + energy_energy_wcs = (energy_coord & time_coord).wcs + data = np.arange(10 * 5).reshape(5, 10) + spec2 = Spectrum(data * u.ph, wcs=energy_energy_wcs, spectral_axis=np.arange(11) * u.keV) + assert spec2.spectral_axis_index == 1 def test_spectrum_from_cubs_wcs_norm_tab(): @@ -98,12 +122,14 @@ def test_spectrum_from_cubs_wcs_norm_tab(): "DATEREF": "2020-01-01T00:00:00", } time_wcs = WCS(header=header) - energy = np.arange(10) * u.keV + energy = (0.5 + np.arange(10)) * u.keV energy_coord = QuantityTableCoordinate(energy, names="energy", physical_types="em.energy") comp_wcs = CompoundLowLevelWCS(time_wcs, energy_coord.wcs) cube = NDCube(np.arange(10 * 5).reshape(10, 5), unit=u.ph, wcs=comp_wcs) - spec = Spectrum(cube) - assert False + spec = Spectrum(cube, spectral_axis=np.arange(11) * u.keV) + assert spec.shape == (10, 5) + assert spec.spectral_axis_index == 0 + assert_quantity_allclose(energy, spec.spectral_axis) def test_slice(): From 3c3fd47a5f35ab6be49a24dc1f6646ec7fb9436a Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Tue, 2 Dec 2025 12:34:30 +0000 Subject: [PATCH 03/13] Spectrum works and supports slicing and basic arithmatic --- sunkit_spex/spectrum/spectrum.py | 101 ++++++++++++-------- sunkit_spex/spectrum/tests/test_spectrum.py | 11 +++ 2 files changed, 73 insertions(+), 39 deletions(-) diff --git a/sunkit_spex/spectrum/spectrum.py b/sunkit_spex/spectrum/spectrum.py index 99276501..ad8cb950 100644 --- a/sunkit_spex/spectrum/spectrum.py +++ b/sunkit_spex/spectrum/spectrum.py @@ -1,3 +1,5 @@ +from copy import deepcopy + import numpy as np from gwcs import WCS as GWCS from gwcs import coordinate_frames as cf @@ -142,10 +144,10 @@ class Spectrum(NDCube): uncertainty : `~astropy.nddata.NDUncertainty` Contains uncertainty information along with propagation rules for spectrum arithmetic. Can take a unit, but if none is given, will use - the unit defined in the flux. + the unit defined in the data. spectral_axis : `~astropy.units.Quantity` or `~specutils.SpectralAxis` Dispersion information with the same shape as the dimension specified by spectral_axis_index - of shape plus one if specifying bin edges. + or shape plus one if specifying bin edges. spectral_axis_index : `int` default 0 The dimension of the data which represents the spectral information default to first dimension index 0. mask : `~numpy.ndarray`-like @@ -329,43 +331,43 @@ def __init__( # If no spectral_axis was provided, create a SpectralCoord based on # the WCS - if spectral_axis is None: - # If the WCS doesn't have a spectral attribute, we assume it's the - # dummy GWCS we created or a solely spectral WCS - if hasattr(self.wcs, "spectral"): - # Handle generated 1D WCS that aren't set to spectral - if not self.wcs.is_spectral and self.wcs.naxis == 1: - spec_axis = self.wcs.pixel_to_world(np.arange(self.data.shape[self.spectral_axis_index])) - else: - spec_axis = self.wcs.spectral.pixel_to_world(np.arange(self.data.shape[self.spectral_axis_index])) - else: - # We now keep the entire GWCS, including spatial information, so we need to include - # all axes in the pixel_to_world call. Note that this assumes/requires that the - # dispersion is the same at all spatial locations. - wcs_args = [np.zeros(self.data.shape[self.spectral_axis_index]) for i in range(len(self.data.shape))] - # Replace with arange for the spectral axis - wcs_args[self.spectral_axis_index] = np.arange(self.data.shape[self.spectral_axis_index]) - wcs_args.reverse() - temp_coords = self.wcs.pixel_to_world(*wcs_args) - # If there are spatial axes, temp_coords will have a SkyCoord and a SpectralCoord - if isinstance(temp_coords, list): - for coords in temp_coords: - if isinstance(coords, SpectralCoord): - spec_axis = coords - break - else: - # WCS axis ordering is reverse of numpy - spec_axis = temp_coords[len(temp_coords) - self.spectral_axis_index - 1] - else: - spec_axis = temp_coords - - try: - if spec_axis.unit.is_equivalent(u.one): - spec_axis = spec_axis * u.pixel - except AttributeError: - raise AttributeError(f"spec_axis does not have unit: {type(spec_axis)} {spec_axis}") - - self._spectral_axis = SpectralAxis(spec_axis) + # if spectral_axis is None: + # # If the WCS doesn't have a spectral attribute, we assume it's the + # # dummy GWCS we created or a solely spectral WCS + # if hasattr(self.wcs, "spectral"): + # # Handle generated 1D WCS that aren't set to spectral + # if not self.wcs.is_spectral and self.wcs.naxis == 1: + # spec_axis = self.wcs.pixel_to_world(np.arange(self.data.shape[self.spectral_axis_index])) + # else: + # spec_axis = self.wcs.spectral.pixel_to_world(np.arange(self.data.shape[self.spectral_axis_index])) + # else: + # # We now keep the entire GWCS, including spatial information, so we need to include + # # all axes in the pixel_to_world call. Note that this assumes/requires that the + # # dispersion is the same at all spatial locations. + # wcs_args = [np.zeros(self.data.shape[self.spectral_axis_index]) for i in range(len(self.data.shape))] + # # Replace with arange for the spectral axis + # wcs_args[self.spectral_axis_index] = np.arange(self.data.shape[self.spectral_axis_index]) + # wcs_args.reverse() + # temp_coords = self.wcs.pixel_to_world(*wcs_args) + # # If there are spatial axes, temp_coords will have a SkyCoord and a SpectralCoord + # if isinstance(temp_coords, list): + # for coords in temp_coords: + # if isinstance(coords, SpectralCoord): + # spec_axis = coords + # break + # else: + # # WCS axis ordering is reverse of numpy + # spec_axis = temp_coords[len(temp_coords) - self.spectral_axis_index - 1] + # else: + # spec_axis = temp_coords + # + # try: + # if spec_axis.unit.is_equivalent(u.one): + # spec_axis = spec_axis * u.pixel + # except AttributeError: + # raise AttributeError(f"spec_axis does not have unit: {type(spec_axis)} {spec_axis}") + # + # self._spectral_axis = SpectralAxis(spec_axis) # make sure that spectral axis is strictly increasing if not np.all(self._spectral_axis[1:] >= self._spectral_axis[:-1]): @@ -384,6 +386,27 @@ def __getitem__(self, item): sliced_spec_axis = self.spectral_axis[item[self.spectral_axis_index]] return Spectrum(sliced_cube, spectral_axis=sliced_spec_axis) + def _slice(self, item): + kwargs = super()._slice(item) + item = tuple(sanitize_slices(item, len(self.shape))) + + kwargs["spectral_axis_index"] = self.spectral_axis_index + kwargs["spectral_axis"] = self.spectral_axis[item[self.spectral_axis_index]] + return kwargs + + def _new_instance(self, **kwargs): + keys = ("unit", "wcs", "mask", "meta", "uncertainty", "psf", "spectral_axis") + full_kwargs = {k: deepcopy(getattr(self, k, None)) for k in keys} + # We Explicitly DO NOT deepcopy any data + full_kwargs["data"] = self.data + full_kwargs.update(kwargs) + new_spectrum = type(self)(**full_kwargs) + if self.extra_coords is not None: + new_spectrum._extra_coords = deepcopy(self.extra_coords) + if self.global_coords is not None: + new_spectrum._global_coords = deepcopy(self.global_coords) + return new_spectrum + @property def spectral_axis(self): return self._spectral_axis diff --git a/sunkit_spex/spectrum/tests/test_spectrum.py b/sunkit_spex/spectrum/tests/test_spectrum.py index 6f462986..79f1ba64 100644 --- a/sunkit_spex/spectrum/tests/test_spectrum.py +++ b/sunkit_spex/spectrum/tests/test_spectrum.py @@ -39,6 +39,17 @@ def test_spectrum_spectral_axis_bin_centers(): assert_array_equal(spec._spectral_axis, [1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5] * u.keV) +def test_spectrum_from_spectrum(): + spec_orig = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=np.arange(1, 12) * u.keV) + spec_new = Spectrum(spec_orig) + spec_orig == spec_new + + +def test_spectrum_unknow_keywords(): + with pytest.raises(ValueError, match=r"Initializer contains unknown arguments."): + Spectrum(np.arange(1, 11) * u.watt, spectral_axis=np.arange(1, 12) * u.keV, mykeyword="myvalue") + + def test_spectrum_from_ndcube_wcs(): header = { "CTYPE1": "TIME ", From 81589b38709ec3bbab856376d41f64ee33530f6b Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Fri, 30 Jan 2026 17:01:09 +0000 Subject: [PATCH 04/13] Tidy up and improve test coverage --- pytest.ini | 2 +- sunkit_spex/spectrum/spectrum.py | 92 ++------- sunkit_spex/spectrum/tests/test_spectrum.py | 207 +++++++++++++++++++- 3 files changed, 222 insertions(+), 79 deletions(-) diff --git a/pytest.ini b/pytest.ini index 7ec6ac41..433540ac 100644 --- a/pytest.ini +++ b/pytest.ini @@ -40,6 +40,6 @@ filterwarnings = # Oldestdeps issues ignore:`finfo.machar` is deprecated:DeprecationWarning ignore:Please use `convolve1d` from the `scipy.ndimage` namespace, the `scipy.ndimage.filters` namespace is deprecated.:DeprecationWarning - ignore::pyparsing.warnings.PyparsingDeprecationWarning ignore::FutureWarning:arviz.* ignore:The isiterable function.*:astropy.utils.exceptions.AstropyDeprecationWarning + ignore:'datfix' made the change:astropy.wcs.wcs.FITSFixedWarning diff --git a/sunkit_spex/spectrum/spectrum.py b/sunkit_spex/spectrum/spectrum.py index ad8cb950..28fddb5f 100644 --- a/sunkit_spex/spectrum/spectrum.py +++ b/sunkit_spex/spectrum/spectrum.py @@ -53,7 +53,7 @@ def pixel_to_world(self, *args, **kwargs): class SpectralAxis(SpectralCoord): - """ + r""" Coordinate object representing spectral values corresponding to a specific spectrum. Overloads SpectralCoord with additional information (currently only bin edges). @@ -77,6 +77,9 @@ def __new__(cls, value, *args, bin_specification="centers", **kwargs): ): raise ValueError("u.pix spectral axes should always be ascending") + if bin_specification == "edges" and value.size < 2: + raise ValueError('If bin_specification="centers" have at least two bin edges.') + # Convert to bin centers if bin edges were given, since SpectralCoord # only accepts centers if bin_specification == "edges": @@ -90,34 +93,22 @@ def __new__(cls, value, *args, bin_specification="centers", **kwargs): return obj - @staticmethod - def _edges_from_centers(centers, unit): - """ - Calculates interior bin edges based on the average of each pair of - centers, with the two outer edges based on extrapolated centers added - to the beginning and end of the spectral axis. - """ - a = np.insert(centers, 0, 2 * centers[0] - centers[1]) - b = np.append(centers, 2 * centers[-1] - centers[-2]) - edges = (a + b) / 2 - return edges * unit - @staticmethod def _centers_from_edges(edges): - """ + r""" Calculates the bin centers as the average of each pair of edges """ return (edges[1:] + edges[:-1]) / 2 @lazyproperty def bin_edges(self): - """ + r""" Calculates bin edges if the spectral axis was created with centers specified. """ - if hasattr(self, "_bin_edges"): + if hasattr(self, "_bin_edges") and self._bin_edges is not None: return self._bin_edges - return self._edges_from_centers(self.value, self.unit) + return None def __array_finalize__(self, obj): super().__array_finalize__(obj) @@ -186,8 +177,8 @@ def __init__( # If the data argument is already a Spectrum (as it would # be for internal arithmetic operations), avoid setup entirely. if isinstance(data, Spectrum): - self._spectral_axis = spectral_axis - self._spectral_dimension = spectral_axis_index + self._spectral_axis = data.spectral_axis + self._spectral_axis_index = data.spectral_axis_index super().__init__(data) return @@ -218,7 +209,7 @@ def __init__( return self._spectral_axis_index = spectral_axis_index - # If here data is should be an array or quantity + # If here data should be an array or quantity if spectral_axis_index is None and data is not None: if data.ndim == 1: self._spectral_axis_index = 0 @@ -267,6 +258,7 @@ def __init__( # If data and spectral axis are both specified, check that their lengths # match or are off by one (implying the spectral axis stores bin edges) + bin_specification = "centers" # default value if data is not None and spectral_axis is not None: if spectral_axis.shape[0] == data.shape[self.spectral_axis_index]: bin_specification = "centers" @@ -314,64 +306,12 @@ def __init__( if wcs is None: wcs = gwcs_from_array(self._spectral_axis) - if wcs is None: - # If no spectral axis or wcs information is provided, initialize - # with an empty gwcs based on the data. - if self.spectral_axis_index is None: - if data.ndim == 1: - self._spectral_axis_index = 0 - else: - raise ValueError("Must specify spectral_axis_index if no WCS or spectral axis is input.") - size = data.shape[self.spectral_axis_index] if not data.isscalar else 1 - wcs = gwcs_from_array( - np.arange(size) * u.Unit("pixel"), data.shape, spectral_axis_index=self.spectral_axis_index - ) - super().__init__(data=data.value if isinstance(data, u.Quantity) else data, wcs=wcs, **kwargs) - # If no spectral_axis was provided, create a SpectralCoord based on - # the WCS - # if spectral_axis is None: - # # If the WCS doesn't have a spectral attribute, we assume it's the - # # dummy GWCS we created or a solely spectral WCS - # if hasattr(self.wcs, "spectral"): - # # Handle generated 1D WCS that aren't set to spectral - # if not self.wcs.is_spectral and self.wcs.naxis == 1: - # spec_axis = self.wcs.pixel_to_world(np.arange(self.data.shape[self.spectral_axis_index])) - # else: - # spec_axis = self.wcs.spectral.pixel_to_world(np.arange(self.data.shape[self.spectral_axis_index])) - # else: - # # We now keep the entire GWCS, including spatial information, so we need to include - # # all axes in the pixel_to_world call. Note that this assumes/requires that the - # # dispersion is the same at all spatial locations. - # wcs_args = [np.zeros(self.data.shape[self.spectral_axis_index]) for i in range(len(self.data.shape))] - # # Replace with arange for the spectral axis - # wcs_args[self.spectral_axis_index] = np.arange(self.data.shape[self.spectral_axis_index]) - # wcs_args.reverse() - # temp_coords = self.wcs.pixel_to_world(*wcs_args) - # # If there are spatial axes, temp_coords will have a SkyCoord and a SpectralCoord - # if isinstance(temp_coords, list): - # for coords in temp_coords: - # if isinstance(coords, SpectralCoord): - # spec_axis = coords - # break - # else: - # # WCS axis ordering is reverse of numpy - # spec_axis = temp_coords[len(temp_coords) - self.spectral_axis_index - 1] - # else: - # spec_axis = temp_coords - # - # try: - # if spec_axis.unit.is_equivalent(u.one): - # spec_axis = spec_axis * u.pixel - # except AttributeError: - # raise AttributeError(f"spec_axis does not have unit: {type(spec_axis)} {spec_axis}") - # - # self._spectral_axis = SpectralAxis(spec_axis) - - # make sure that spectral axis is strictly increasing - if not np.all(self._spectral_axis[1:] >= self._spectral_axis[:-1]): - raise ValueError("Spectral axis must be strictly increasing or decreasing.") + # make sure that spectral axis is strictly increasing or strictly decreasing + is_strictly_increasing = np.all(self._spectral_axis[1:] > self._spectral_axis[:-1]) + if len(self._spectral_axis) > 1 and not (is_strictly_increasing): + raise ValueError("Spectral axis must be strictly increasing decreasing.") if hasattr(self, "uncertainty") and self.uncertainty is not None: if not data.shape == self.uncertainty.array.shape: diff --git a/sunkit_spex/spectrum/tests/test_spectrum.py b/sunkit_spex/spectrum/tests/test_spectrum.py index 79f1ba64..79f896fe 100644 --- a/sunkit_spex/spectrum/tests/test_spectrum.py +++ b/sunkit_spex/spectrum/tests/test_spectrum.py @@ -8,11 +8,13 @@ from numpy.testing import assert_array_equal import astropy.units as u +from astropy.coordinates import SpectralCoord +from astropy.nddata import StdDevUncertainty from astropy.tests.helper import assert_quantity_allclose from astropy.time import Time from astropy.wcs import WCS -from sunkit_spex.spectrum.spectrum import SpectralAxis, Spectrum +from sunkit_spex.spectrum.spectrum import SpectralAxis, Spectrum, gwcs_from_array rng = np.random.default_rng() @@ -73,7 +75,7 @@ def test_spectrum_from_ndcube_wcs(): data = np.arange(np.prod(shape), dtype=int).reshape(shape) cube = NDCube(data.T, wcs=wcs) - with pytest.raises(ValueError, match="Input NDCube missing unit.*"): + with pytest.raises(ValueError, match=r"Input NDCube missing unit.*"): Spectrum(cube) cube = NDCube(data, wcs=wcs, unit=u.ph) @@ -163,3 +165,204 @@ def test_arithmetic_operators(op, value, res): spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=(np.arange(1, 11) + 0.5) * u.keV) res_spec = op(spec, value) assert_array_equal(res_spec.data, res) + + +def test_spectral_axis_bin_edges_from_centers(): + """Test that bin_edges are correctly calculated when SpectralAxis is created with centers.""" + spec_axis = SpectralAxis(np.array([1.5, 2.5, 3.5, 4.5]) * u.keV, bin_specification="centers") + edges = spec_axis.bin_edges + assert edges is None + + +def test_spectral_axis_bin_edges_preserved(): + """Test that bin_edges are preserved when SpectralAxis is created with edges.""" + input_edges = np.array([1.0, 2.0, 3.0, 4.0, 5.0]) * u.keV + spec_axis = SpectralAxis(input_edges, bin_specification="edges") + assert_quantity_allclose(spec_axis.bin_edges, input_edges) + + +def test_spectral_axis_centers_from_edges(): + """Test that centers are correctly calculated from edges.""" + input_edges = np.array([1.0, 2.0, 3.0, 4.0, 5.0]) * u.keV + spec_axis = SpectralAxis(input_edges, bin_specification="edges") + assert_quantity_allclose(spec_axis, [1.5, 2.5, 3.5, 4.5] * u.keV) + + +def test_spectral_axis_single_center(): + """Test SpectralAxis handles single-element arrays.""" + spec_axis = SpectralAxis(np.array([5.0]) * u.keV, bin_specification="centers") + edges = spec_axis.bin_edges + assert edges is None + + +def test_spectral_axis_single_bin(): + """Test SpectralAxis handles single bins""" + with pytest.raises(ValueError, match="If bin_specification"): + SpectralAxis(np.array([5.0]) * u.keV, bin_specification="edges") + + spec_axis = SpectralAxis(np.array([5.0, 6.0]) * u.keV, bin_specification="edges") + edges = spec_axis.bin_edges + assert edges is not None + assert len(edges) == 2 + assert spec_axis == 5.5 * u.keV + + +def test_spectral_axis_empty_array(): + """Test SpectralAxis handles empty arrays.""" + edges = SpectralAxis(np.array([]), u.keV) + assert len(edges) == 0 + + +def test_spectral_axis_pixel_ascending(): + """Test that pixel spectral axes must be ascending.""" + with pytest.raises(ValueError, match=r"u\.pix spectral axes should always be ascending"): + SpectralAxis(np.array([5, 4, 3, 2, 1]) * u.pix) + + +def test_spectral_axis_pixel_ascending_valid(): + """Test that ascending pixel spectral axes are accepted.""" + spec_axis = SpectralAxis(np.array([1, 2, 3, 4, 5]) * u.pix) + assert len(spec_axis) == 5 + + +def test_spectrum_from_spectrum_inherits_attributes(): + """Test that Spectrum created from another Spectrum inherits spectral_axis and spectral_axis_index.""" + spec_orig = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=np.arange(1, 12) * u.keV) + spec_new = Spectrum(spec_orig) + + # Verify spectral_axis_index is inherited (Bug #1 fix) + assert spec_new.spectral_axis_index == spec_orig.spectral_axis_index + assert spec_new.spectral_axis_index == 0 + + # Verify spectral_axis is inherited + assert_quantity_allclose(spec_new.spectral_axis, spec_orig.spectral_axis) + + +def test_spectrum_from_spectrum_preserves_data(): + """Test that Spectrum created from another Spectrum preserves data.""" + data = np.arange(1, 11) * u.watt + spec_orig = Spectrum(data, spectral_axis=np.arange(1, 12) * u.keV) + spec_new = Spectrum(spec_orig) + + assert_array_equal(spec_new.data, spec_orig.data) + assert spec_new.unit == spec_orig.unit + + +def test_spectrum_strictly_increasing_spectral_axis(): + """Test that strictly increasing spectral axis is accepted.""" + spec = Spectrum(np.arange(1, 6) * u.watt, spectral_axis=np.array([1, 2, 3, 4, 5]) * u.keV) + assert_quantity_allclose(spec.spectral_axis, [1, 2, 3, 4, 5] * u.keV) + + +def test_spectrum_non_monotonic_spectral_axis_raises(): + """Test that non-monotonic spectral axis raises ValueError.""" + with pytest.raises(ValueError, match="strictly increasing"): + Spectrum(np.arange(1, 6) * u.watt, spectral_axis=np.array([1, 3, 2, 4, 5]) * u.keV) + + +def test_spectrum_duplicate_values_in_spectral_axis_raises(): + """Test that duplicate values in spectral axis raises ValueError.""" + with pytest.raises(ValueError, match="strictly increasing"): + Spectrum(np.arange(1, 5) * u.watt, spectral_axis=np.array([1, 2, 2, 3]) * u.keV) + + +def test_spectrum_single_element_spectral_axis(): + """Test that single-element spectral axis is accepted.""" + spec = Spectrum(np.array([5]) * u.watt, spectral_axis=np.array([10]) * u.keV) + assert spec.shape == (1,) + assert_quantity_allclose(spec.spectral_axis, [10] * u.keV) + + +def test_spectrum_spectral_axis_length_mismatch(): + """Test that mismatched spectral axis length raises ValueError.""" + with pytest.raises(ValueError, match="Spectral axis length"): + Spectrum(np.arange(1, 11) * u.watt, spectral_axis=np.arange(1, 5) * u.keV) + + +def test_spectrum_uncertainty_shape_mismatch(): + """Test that mismatched uncertainty shape raises ValueError.""" + data = np.arange(1, 11) * u.watt + uncertainty = StdDevUncertainty(np.arange(1, 6)) # Wrong shape + with pytest.raises(ValueError, match=r"Data axis .* and uncertainty .* shapes must be the same"): + Spectrum(data, spectral_axis=np.arange(1, 12) * u.keV, uncertainty=uncertainty) + + +def test_spectrum_with_valid_uncertainty(): + """Test Spectrum with correctly shaped uncertainty.""" + data = np.arange(1, 11) * u.watt + uncertainty = StdDevUncertainty(np.ones(10) * 0.1) + spec = Spectrum(data, spectral_axis=np.arange(1, 12) * u.keV, uncertainty=uncertainty) + assert spec.uncertainty is not None + assert spec.uncertainty.array.shape == data.shape + + +def test_gwcs_from_array_ascending(): + """Test gwcs_from_array with ascending array.""" + array = np.array([1, 2, 3, 4, 5]) * u.keV + wcs = gwcs_from_array(array) + + # Test pixel_to_world + result = wcs.pixel_to_world(0) + assert_quantity_allclose(result, 1 * u.keV) + + result = wcs.pixel_to_world(4) + assert_quantity_allclose(result, 5 * u.keV) + + +def test_gwcs_from_array_descending(): + """Test gwcs_from_array with descending array.""" + array = np.array([5, 4, 3, 2, 1]) * u.keV + wcs = gwcs_from_array(array) + + # Test pixel_to_world + result = wcs.pixel_to_world(0) + assert_quantity_allclose(result, 5 * u.keV) + + result = wcs.pixel_to_world(4) + assert_quantity_allclose(result, 1 * u.keV) + + +def test_gwcs_from_array_world_to_pixel(): + """Test gwcs_from_array world_to_pixel transformation.""" + array = np.array([1, 2, 3, 4, 5]) * u.keV + wcs = gwcs_from_array(array) + + # Test world_to_pixel - requires SpectralCoord + result = wcs.world_to_pixel(SpectralCoord(3 * u.keV)) + assert np.isclose(result, 2.0) + + +def test_slice_preserves_spectral_axis_index(): + """Test that slicing preserves spectral_axis_index.""" + spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=(np.arange(1, 11) + 0.5) * u.keV) + sliced = spec[2:7] + assert sliced.spectral_axis_index == spec.spectral_axis_index + + +def test_slice_updates_spectral_axis(): + """Test that slicing correctly updates spectral_axis.""" + spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=(np.arange(1, 11) + 0.5) * u.keV) + sliced = spec[2:5] + assert_quantity_allclose(sliced.spectral_axis, [3.5, 4.5, 5.5] * u.keV) + + +def test_slice_single_element(): + """Test slicing to a single element.""" + spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=(np.arange(1, 11) + 0.5) * u.keV) + sliced = spec[5:6] + assert sliced.shape == (1,) + assert_quantity_allclose(sliced.spectral_axis, [6.5] * u.keV) + + +def test_arithmetic_preserves_spectral_axis(): + """Test that arithmetic operations preserve spectral_axis.""" + spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=(np.arange(1, 11) + 0.5) * u.keV) + result = spec + 1 * u.watt + assert_quantity_allclose(result.spectral_axis, spec.spectral_axis) + + +def test_arithmetic_preserves_spectral_axis_index(): + """Test that arithmetic operations preserve spectral_axis_index.""" + spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=(np.arange(1, 11) + 0.5) * u.keV) + result = spec * 2 + assert result.spectral_axis_index == spec.spectral_axis_index From 0e2c88593a18106425b30dae1df4e07ee5fdb5e5 Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Fri, 6 Feb 2026 12:49:58 +0000 Subject: [PATCH 05/13] 1D Specrta working --- .gitignore | 1 + docs/reference/index.rst | 4 +- docs/reference/spectrum.rst | 7 + examples/spectrum.py | 220 ++++++++++++++++++++ pyproject.toml | 2 +- sunkit_spex/spectrum/spectrum.py | 204 +++++++++++++++--- sunkit_spex/spectrum/tests/test_spectrum.py | 135 ++++++++---- 7 files changed, 502 insertions(+), 71 deletions(-) create mode 100644 docs/reference/spectrum.rst create mode 100644 examples/spectrum.py diff --git a/.gitignore b/.gitignore index b52b152a..b2197db0 100644 --- a/.gitignore +++ b/.gitignore @@ -267,3 +267,4 @@ package.json # Log files generated by 'vagrant up' *.log +.worktree diff --git a/docs/reference/index.rst b/docs/reference/index.rst index 6c60dd2d..fc15b38d 100644 --- a/docs/reference/index.rst +++ b/docs/reference/index.rst @@ -9,7 +9,7 @@ Software and API. .. toctree:: :maxdepth: 2 - - fitting + spectrum models + fitting legacy diff --git a/docs/reference/spectrum.rst b/docs/reference/spectrum.rst new file mode 100644 index 00000000..5ce9a9e8 --- /dev/null +++ b/docs/reference/spectrum.rst @@ -0,0 +1,7 @@ +Models (`sunkit_spex.spectrum`) +******************************* + +``sunkit_spex.spectrum`` module contains objects for holding spectral data + +.. automodapi:: sunkit_spex.spectrum +.. automodapi:: sunkit_spex.spectrum.spectrum diff --git a/examples/spectrum.py b/examples/spectrum.py new file mode 100644 index 00000000..a6a68aad --- /dev/null +++ b/examples/spectrum.py @@ -0,0 +1,220 @@ +""" +======== +Spectrum +======== + +This example will demonstrate how to store spectral data in `~sunkit_spex.spectrum.Specutm` container +""" + +##################################################### +# +# Imports + +import numpy as np +from ndcube import NDMeta + +import astropy.units as u +from astropy.coordinates import SpectralCoord +from astropy.time import Time + +from sunkit_spex.spectrum import Spectrum + +rng = np.random.default_rng() +##################################################### +# +# 1D Spectrum +# ----------- +# Let's being with the simplest case a single spectrum that is a series of measurements as function of wavelength or +# energy. We will start of by creating some synthetic data and corresponding energy bins as well as some important metadata +# in this case the exposure time. + +data = rng.rand(50) * u.ct +energy = np.linspace(1, 50, 50) * u.keV +time = Time("2025-02-18T15:08") + +exposure_time = 5 * u.s + +##################################################### +# +# Once we have our synthetic data we can create our metadata container `NDMeta` and `Spectrum` object. + +meta = NDMeta() +meta.add("exposure_time", exposure_time) +meta.add("date-obs", time) + +spec_1d = Spectrum(data, spectral_axis=energy, meta=meta) +spec_1d + +##################################################### +# +# One of the key feature of the `Spectrum` object is the ability to slice, crop and perform other operations using +# standard sliceing methods: + +spec_1d_sliced = spec_1d[10:20] +spec_1d_sliced.shape +spec_1d_sliced.axis_world_coords_values() +spec_1d_sliced.meta +spec_1d_sliced.spectral_axis + +##################################################### +# +# High level coordinate objects such as SkyCoord and SpectralCoord + +spec_1d_crop = spec_1d.crop(SpectralCoord(10.5, unit=u.keV), SpectralCoord(20, unit=u.keV)) +spec_1d_crop.shape +spec_1d_crop.axis_world_coords_values() +spec_1d_crop.meta +spec_1d_crop.spectral_axis + +##################################################### +# +# And Quantities + +spec_1d_crop_value = spec_1d.crop_by_values((10.5 * u.keV), (20.5 * u.keV)) +spec_1d_crop_value.shape +spec_1d_crop_value.axis_world_coords_values() +spec_1d_crop_value.meta +spec_1d_crop_value.spectral_axis + +##################################################### +# +# 2D Spectrum (spectrogram or time v energy) +# ------------------------------------------ +# Let build on the previous example by increasing the dimensionality of the data in this case to a spectrogram or a +# series of spectra as a function of time. Here we will simulate a series of 10 spectra taken over 10 minutes. Again we +# begin by creating our synthetic data as before but additionally creating the time variable. + +# data = rng.random(10, 50) * u.ct +# energy = np.linspace(1, 50, 51) * u.keV +# times = Time("2025-02-18T15:08") + np.arange(10) * u.min +# exposure_time = np.arange(5, 15) * u.s + +##################################################### +# +# We are also going to demonstrate the power of the sliceable metadata, so in this example each of the individual +# spectra have different exposure times (this could be another important information regard the observation) + +# meta = NDMeta() +# meta.add("exposure_time", exposure_time, axes=(0,)) +# +# time_coord = TimeTableCoordinate(times, names="time", physical_types="time") +# energy_coord = QuantityTableCoordinate(energy, names="energy", physical_types="em.energy") +# wcs = (energy_coord | time_coord).wcs +# +# spec_2d_time_energy = Spectrum(data, spectral_axis=energy, wcs=wcs, spectral_axis_index=1, meta=meta) +# spec_2d_time_energy.axis_world_coords_values() + +###################################################### +# +# For this example we're assuming we don't have a full WCS for the data so we utilise `extra_coords` to store the +# temporal axes. + +# spec_2d_time_energy.extra_coords.add("time", (0,), times) +# +# spec_2d_time_energy +# spec_2d_time_energy.shape +# spec_2d_time_energy.axis_world_coords_values() +# spec_2d_time_energy.meta + +###################################################### +# +# Again all standard slicing works + +# spec_2d_time_energy[2:5] +# spec_2d_time_energy[:, 10:20] +# spec_2d_time_energy_sliced = spec_2d_time_energy[2:5, 10:20] + +###################################################### +# +# We can being to see the usefulness of the sliceable metadata notice how the exposure time entry has been sliced +# appropriately + +# spec_2d_time_energy_sliced +# spec_2d_time_energy_sliced.shape +# spec_2d_time_energy_sliced.axis_world_coords_values() +# spec_2d_time_energy_sliced.meta + +###################################################### +# +# The same can be archived using coordinate objects or values + +# spec_2d_time_energy_crop = spec_2d_time_energy.crop( +# [Time("2025-02-18T15:09"), SpectralCoord(10, unit=u.keV)], +# [Time("2025-02-18T15:014"), SpectralCoord(20, unit=u.keV)], +# ) +# spec_2d_time_energy_crop.crop_by_values() + +##################################################### +# +# 2D Spectrum ( e.g. detector v energy) +# + +# data = rng.rand(10, 50) * u.ct +# energy = np.linspace(1, 50, 50) * u.keV +# +# exposure_time = np.arange(10) * u.s +# labels = np.array([chr(97 + i) for i in range(10)]) +# +# meta = NDMeta() +# meta.add("exposure_time", exposure_time, axes=0) +# meta.add("detector", labels, axes=0) +# +# spec_2d_det_time = Spectrum(data, spectral_axis=energy, spectral_axis_index=1, meta=meta) +# spec_2d_det_time + +##################################################### +# +# 3D Spectrum ( e.g. detector v energy v time) +# + +# data = rng.random(10, 20, 30) * u.ct +# energy = np.linspace(1, 31, 31) * u.keV +# +# labels = np.array([chr(97 + i) for i in range(10)]) +# exposure_time = np.arange(10 * 20).reshape(10, 20) * u.s +# times = Time.now() + np.arange(20) * u.s +# +# meta = NDMeta() +# meta.add("exposure_time", exposure_time, axes=(0, 1)) +# meta.add("detector", labels, axes=(0,)) +# +# spec_3d_det_energy_time = Spectrum(data, spectral_axis=energy, spectral_axis_index=2, meta=meta) +# spec_3d_det_energy_time.extra_coords.add("time", (0,), times) +# +# spec_3d_det_energy_time[:, 10:15, :].meta +# spec_3d_det_energy_time[2:3, 10:15, :].meta + +##################################################### +# +# 4D Spectrum ( e.g. spatial v spatial v energy v time) +# + +# import numpy as np +# from ndcube import NDMeta +# +# import astropy.units as u +# from astropy.time import Time +# +# data = np.random.rand(10, 10, 20, 30) * u.ct +# energy = np.linspace(1, 31, 31) * u.keV +# exposure_time = np.arange(20) * u.s +# times = Time.now() + np.arange(20) * u.s +# +# meta = NDMeta() +# meta.add("exposure_time", exposure_time, axes=(2,)) +# +# wcs = astropy.wcs.WCS(naxis=2) +# wcs.wcs.ctype = "HPLT-TAN", "HPLN-TAN" +# wcs.wcs.cunit = "deg", "deg" +# wcs.wcs.cdelt = 0.5, 0.4 +# wcs.wcs.crpix = 5, 6 +# wcs.wcs.crval = 0.5, 1 +# wcs.wcs.cname = "HPC lat", "HPC lon" +# +# cube = NDCube(data=data, wcs=wcs, meta=meta) +# +# # Now instantiate the NDCube +# spec_4d_lon_lat_time_energy = Spectrum(data, wcs=wcs, spectral_axis=energy, spectral_axis_index=3, meta=meta) +# spec_4d_lon_lat_time_energy.extra_coords.add("time", (2,), times) +# +# spec_4d_lon_lat_time_energy diff --git a/pyproject.toml b/pyproject.toml index be250874..657b6020 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,7 +26,7 @@ dependencies = [ "scipy>=1.12", "sunpy>=7.0", "xarray>=2023.12", - "gwcs>=0.21.0", + "gwcs<=1.0.0", #until https://github.com/sunpy/ndcube/issues/913 is resolved "ndcube>=2.3", ] diff --git a/sunkit_spex/spectrum/spectrum.py b/sunkit_spex/spectrum/spectrum.py index 28fddb5f..6ed8eee2 100644 --- a/sunkit_spex/spectrum/spectrum.py +++ b/sunkit_spex/spectrum/spectrum.py @@ -1,3 +1,4 @@ +import copy from copy import deepcopy import numpy as np @@ -7,6 +8,7 @@ import astropy.units as u from astropy.coordinates import SpectralCoord +from astropy.modeling.mappings import Identity, Mapping from astropy.modeling.tabular import Tabular1D from astropy.utils import lazyproperty from astropy.wcs.wcsapi import sanitize_slices @@ -17,40 +19,183 @@ __doctest_requires__ = {"Spectrum": ["ndcube>=2.3"]} -def gwcs_from_array(array): +class SpectralGWCS(GWCS): + """ + This is a placeholder lookup-table GWCS created when a :class:`~specutils.Spectrum` is + instantiated with a ``spectral_axis`` and no WCS. + """ + + def __init__(self, *args, **kwargs): + self.original_unit = kwargs.pop("original_unit", "") + super().__init__(*args, **kwargs) + + def copy(self): + """ + Return a shallow copy of the object. + + Convenience method so user doesn't have to import the + :mod:`copy` stdlib module. + + .. warning:: + Use `deepcopy` instead of `copy` unless you know why you need a + shallow copy. + """ + return copy.copy(self) + + def deepcopy(self): + """ + Return a deep copy of the object. + + Convenience method so user doesn't have to import the + :mod:`copy` stdlib module. + """ + return copy.deepcopy(self) + + +def gwcs_from_array(array, flux_shape, spectral_axis_index=None): """ Create a new WCS from provided tabular data. This defaults to being - a GWCS object. + a GWCS object with a lookup table for the spectral axis and filler + pixel to pixel identity conversions for spatial axes, if they exist. """ orig_array = u.Quantity(array) + naxes = len(flux_shape) + + if naxes > 1: + if spectral_axis_index is None: + raise ValueError("spectral_axis_index must be set for multidimensional flux arrays") + # Axis order is reversed for WCS from numpy array + spectral_axis_index = naxes - spectral_axis_index - 1 + elif naxes == 1: + spectral_axis_index = 0 + + axes_order = list(np.arange(naxes)) + axes_type = [ + "SPATIAL", + ] * naxes + axes_type[spectral_axis_index] = "SPECTRAL" + + detector_frame = cf.CoordinateFrame( + naxes=naxes, + name="detector", + unit=[ + u.pix, + ] + * naxes, + axes_order=axes_order, + axes_type=axes_type, + ) + + if array.unit in ("", "pix", "pixel"): + # Spectrum was initialized without a wcs or spectral axis + spectral_frame = cf.CoordinateFrame( + naxes=1, + unit=[ + array.unit, + ], + axes_type=[ + "Spectral", + ], + axes_order=(spectral_axis_index,), + ) + else: + phys_types = None + # Note that some units have multiple physical types, so we can't just set the + # axis name to the physical type string. + if array.unit.physical_type == "length": + axes_names = [ + "wavelength", + ] + elif array.unit.physical_type == "frequency": + axes_names = [ + "frequency", + ] + elif array.unit.physical_type == "velocity": + axes_names = [ + "velocity", + ] + phys_types = [ + "spect.dopplerVeloc.optical", + ] + elif array.unit.physical_type == "wavenumber": + axes_names = [ + "wavenumber", + ] + elif array.unit.physical_type == "energy": + axes_names = [ + "energy", + ] + else: + raise ValueError("Spectral axis units must be one of length,frequency, velocity, energy, or wavenumber") - coord_frame = cf.CoordinateFrame(naxes=1, axes_type=("SPECTRAL",), axes_order=(0,)) - spec_frame = cf.SpectralFrame(unit=array.unit, axes_order=(0,)) + spectral_frame = cf.SpectralFrame( + unit=array.unit, axes_order=(spectral_axis_index,), axes_names=axes_names, axis_physical_types=phys_types + ) + + if naxes > 1: + axes_order.remove(spectral_axis_index) + spatial_frame = cf.CoordinateFrame( + naxes=naxes - 1, + unit=[ + "", + ] + * (naxes - 1), + axes_type=[ + "Spatial", + ] + * (naxes - 1), + axes_order=axes_order, + ) + output_frame = cf.CompositeFrame(frames=[spatial_frame, spectral_frame]) + else: + output_frame = spectral_frame # In order for the world_to_pixel transformation to automatically convert - # input units, the equivalencies in the lookup table have to be extended + # input units, the equivalencies in the look up table have to be extended # with spectral unit information. - SpectralTabular1D = type("SpectralTabular1D", (Tabular1D,), {"input_units_equivalencies": {"x0": u.spectral()}}) + SpectralTabular1D = type( + "SpectralTabular1D", (Tabular1D,), {"input_units_equivalencies": {"x0": u.spectral()}, "bounds_error": True} + ) + + # We pass through the pixel values of spatial axes with Identity and use a lookup + # table for the spectral axis values. We use Mapping to pipe the values to the correct + # model depending on which axis is the spectral axis + if naxes == 1: + forward_transform = SpectralTabular1D(np.arange(len(array)) * u.pix, lookup_table=array) + else: + axes_order.append(spectral_axis_index) + # WCS axis order is reverse of numpy array order + mapped_axes = axes_order + out_mapping = np.ones(len(mapped_axes)).astype(int) + for i in range(len(mapped_axes)): + out_mapping[mapped_axes[i]] = i + forward_transform = ( + Mapping(mapped_axes) + | Identity(naxes - 1) & SpectralTabular1D(np.arange(len(array)) * u.pix, lookup_table=array) + | Mapping(out_mapping) + ) - forward_transform = SpectralTabular1D(np.arange(len(array)), lookup_table=array) # If our spectral axis is in descending order, we have to flip the lookup # table to be ascending in order for world_to_pixel to work. if len(array) == 0 or array[-1] > array[0]: - forward_transform.inverse = SpectralTabular1D(array, lookup_table=np.arange(len(array))) + forward_transform.inverse = SpectralTabular1D(array, lookup_table=np.arange(len(array)) * u.pix) else: - forward_transform.inverse = SpectralTabular1D(array[::-1], lookup_table=np.arange(len(array))[::-1]) - - class SpectralGWCS(GWCS): - def pixel_to_world(self, *args, **kwargs): - if orig_array.unit == "": - return u.Quantity(super().pixel_to_world_values(*args, **kwargs)) - return super().pixel_to_world(*args, **kwargs).to(orig_array.unit, equivalencies=u.spectral()) + raise ValueError("Unsupported ") + # forward_transform.inverse = SpectralTabular1D( + # array[::-1], lookup_table=np.arange(len(array))[::-1]) - return SpectralGWCS(forward_transform=forward_transform, input_frame=coord_frame, output_frame=spec_frame) + tabular_gwcs = SpectralGWCS( + original_unit=orig_array.unit, + forward_transform=forward_transform, + input_frame=detector_frame, + output_frame=output_frame, + ) # Store the intended unit from the origin input array # tabular_gwcs._input_unit = orig_array.unit + return tabular_gwcs # noqa: RET504 + class SpectralAxis(SpectralCoord): r""" @@ -117,10 +262,10 @@ def __array_finalize__(self, obj): class Spectrum(NDCube): r""" - Spectrum container for data with one spectral axis. + Spectrum container for data which share a common spectral axis. Note that "1D" in this case refers to the fact that there is only one - spectral axis. `Spectrum` can contain ND data where + spectral axis. `Spectrum` can contain ND data where ``data`` have a shape with dimension greater than 1. Notes @@ -290,21 +435,24 @@ def __init__( # Check the spectral_axis matches the wcs if wcs is not None: - if hasattr(wcs, "spectral"): + wsc_coords = None + if hasattr(wcs, "spectral") and getattr(wcs, "is_spectral", False): wcs_coords = wcs.spectral.pixel_to_world(np.arange(data.shape[self.spectral_axis_index])).to("keV") elif wcs.pixel_n_dim == 1: wcs_coords = wcs.pixel_to_world(np.arange(data.shape[self.spectral_axis_index])) - else: - array_index = wcs.pixel_n_dim - self._spectral_axis_index - 1 - pixels = [0] * wcs.pixel_n_dim - pixels[array_index] = np.arange(data.shape[self.spectral_axis_index]) - wcs_coords = wcs.pixel_to_world(*pixels)[array_index] - - if not u.allclose(self._spectral_axis, wcs_coords): - raise ValueError(f"Spectral axis {self._spectral_axis} and wcs spectral axis {wcs_coords} must match.") + # else: + # array_index = wcs.pixel_n_dim - self._spectral_axis_index - 1 + # pixels = [0] * wcs.pixel_n_dim + # pixels[array_index] = np.arange(data.shape[self.spectral_axis_index]) + # wcs_coords = wcs.pixel_to_world(*pixels)[array_index] + if wsc_coords is not None: + if not u.allclose(self._spectral_axis, wcs_coords): + raise ValueError( + f"Spectral axis {self._spectral_axis} and wcs spectral axis {wcs_coords} must match." + ) if wcs is None: - wcs = gwcs_from_array(self._spectral_axis) + wcs = gwcs_from_array(self._spectral_axis, data.shape, spectral_axis_index=self.spectral_axis_index) super().__init__(data=data.value if isinstance(data, u.Quantity) else data, wcs=wcs, **kwargs) diff --git a/sunkit_spex/spectrum/tests/test_spectrum.py b/sunkit_spex/spectrum/tests/test_spectrum.py index 79f896fe..81d89100 100644 --- a/sunkit_spex/spectrum/tests/test_spectrum.py +++ b/sunkit_spex/spectrum/tests/test_spectrum.py @@ -2,23 +2,114 @@ import numpy as np import pytest +from gwcs import coordinate_frames as cf from ndcube import NDCube from ndcube.extra_coords import QuantityTableCoordinate, TimeTableCoordinate from ndcube.wcs.wrappers import CompoundLowLevelWCS from numpy.testing import assert_array_equal import astropy.units as u -from astropy.coordinates import SpectralCoord +from astropy.modeling import models from astropy.nddata import StdDevUncertainty from astropy.tests.helper import assert_quantity_allclose from astropy.time import Time from astropy.wcs import WCS -from sunkit_spex.spectrum.spectrum import SpectralAxis, Spectrum, gwcs_from_array +from sunkit_spex.spectrum.spectrum import SpectralAxis, SpectralGWCS, Spectrum, gwcs_from_array rng = np.random.default_rng() +def test_spectral_gwcs_init_and_copy(): + # Setup dummy transform and frames + trans = models.Identity(1) + # Create distinct frames with unique names + # Usually 'detector' or 'pixel' for input, and 'world' or 'spectral' for output + input_frame = cf.CoordinateFrame(naxes=1, axes_type=["SPECTRAL"], axes_order=(0,), name="pixel_frame") + + output_frame = cf.CoordinateFrame(naxes=1, axes_type=["SPECTRAL"], axes_order=(0,), name="world_frame") + sgwcs = SpectralGWCS( + forward_transform=trans, input_frame=input_frame, output_frame=output_frame, original_unit="Angstrom" + ) + + assert sgwcs.original_unit == "Angstrom" + + # Test shallow copy + sgwcs_copy = sgwcs.copy() + assert sgwcs_copy.original_unit == sgwcs.original_unit + assert sgwcs_copy is not sgwcs + + # Test deep copy + sgwcs_deepcopy = sgwcs.deepcopy() + assert sgwcs_deepcopy.original_unit == sgwcs.original_unit + assert sgwcs_deepcopy is not sgwcs + + +def test_gwcs_from_array_1d_wavelength(): + wavelengths = np.linspace(4000, 7000, 100) * u.AA + flux_shape = (100,) + + wcs = gwcs_from_array(wavelengths, flux_shape) + + assert isinstance(wcs, SpectralGWCS) + assert wcs.output_frame.unit[0] == u.AA + assert wcs.output_frame.axes_names[0] == "wavelength" + + # Test forward transform (pixel to world) + assert np.allclose(wcs(0), 4000) + assert np.allclose(wcs(99), 7000) + + # Test inverse transform (world to pixel) + assert np.allclose(wcs.invert(4000), 0) + + +def test_gwcs_from_array_3d_cube(): + # 3D cube: (Spatial, Spatial, Spectral) -> (y, x, lambda) + # In numpy: shape is (ny, nx, nlambda) + # We want spectral axis to be index 2 + n_lambda = 50 + flux_shape = (10, 20, n_lambda) + freqs = np.linspace(100, 200, n_lambda) * u.keV + + # Note: spectral_axis_index is relative to numpy shape + wcs = gwcs_from_array(freqs, flux_shape, spectral_axis_index=2) + + assert wcs.output_frame.naxes == 3 + assert wcs.forward_transform.n_inputs == 3 + + # Test mapping: (x, y, lambda_pix) -> (spatial, spatial, freq) + # GWCS/WCS usually expects (x, y, z) input order + world = wcs.pixel_to_world(0, 0, 0) # pixels for x, y, lambda + assert world[0] == 100 * u.keV + assert wcs.output_frame.frames[1].axes_names[0] == "energy" + + +def test_gwcs_from_array_descending(): + # Test descending spectral axis for inverse transform logic + waves = np.array([5000, 4000, 3000]) * u.AA + wcs = gwcs_from_array(waves, (3,)) + + assert np.allclose(wcs(0), 5000) + assert np.allclose(wcs(2), 3000) + + # Check that inverse works correctly on flipped table + assert np.allclose(wcs.invert(3000), 2) + assert np.allclose(wcs.invert(5000), 0) + + +def test_gwcs_from_array_invalid_units(): + data = np.arange(10) * u.Jy # Flux units are not valid for spectral axis + with pytest.raises(ValueError, match="Spectral axis units must be one of"): + gwcs_from_array(data, (10,)) + + +def test_gwcs_from_array_missing_index(): + data = np.linspace(1, 10, 10) * u.m + # 2D flux but no index provided + with pytest.raises(ValueError, match="spectral_axis_index must be set"): + gwcs_from_array(data, (10, 10)) + + def test_spectrum_quantity_bin_edges(): spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=np.arange(1, 12) * u.keV) assert_array_equal(spec._spectral_axis, [1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5] * u.keV) @@ -84,7 +175,7 @@ def test_spectrum_from_ndcube_wcs(): Spectrum(cube) with pytest.raises(ValueError, match=r"Spectral axis"): - Spectrum(cube, spectral_axis=spec_axis_edges[:-1]) + Spectrum(cube, spectral_axis=spec_axis_edges[1:-1]) spec = Spectrum(cube, spectral_axis=spec_axis_edges) assert isinstance(spec, Spectrum) @@ -296,42 +387,6 @@ def test_spectrum_with_valid_uncertainty(): assert spec.uncertainty.array.shape == data.shape -def test_gwcs_from_array_ascending(): - """Test gwcs_from_array with ascending array.""" - array = np.array([1, 2, 3, 4, 5]) * u.keV - wcs = gwcs_from_array(array) - - # Test pixel_to_world - result = wcs.pixel_to_world(0) - assert_quantity_allclose(result, 1 * u.keV) - - result = wcs.pixel_to_world(4) - assert_quantity_allclose(result, 5 * u.keV) - - -def test_gwcs_from_array_descending(): - """Test gwcs_from_array with descending array.""" - array = np.array([5, 4, 3, 2, 1]) * u.keV - wcs = gwcs_from_array(array) - - # Test pixel_to_world - result = wcs.pixel_to_world(0) - assert_quantity_allclose(result, 5 * u.keV) - - result = wcs.pixel_to_world(4) - assert_quantity_allclose(result, 1 * u.keV) - - -def test_gwcs_from_array_world_to_pixel(): - """Test gwcs_from_array world_to_pixel transformation.""" - array = np.array([1, 2, 3, 4, 5]) * u.keV - wcs = gwcs_from_array(array) - - # Test world_to_pixel - requires SpectralCoord - result = wcs.world_to_pixel(SpectralCoord(3 * u.keV)) - assert np.isclose(result, 2.0) - - def test_slice_preserves_spectral_axis_index(): """Test that slicing preserves spectral_axis_index.""" spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=(np.arange(1, 11) + 0.5) * u.keV) @@ -340,7 +395,7 @@ def test_slice_preserves_spectral_axis_index(): def test_slice_updates_spectral_axis(): - """Test that slicing correctly updates spectral_axis.""" + """Test that slicing correctly slices spectral_axis.""" spec = Spectrum(np.arange(1, 11) * u.watt, spectral_axis=(np.arange(1, 11) + 0.5) * u.keV) sliced = spec[2:5] assert_quantity_allclose(sliced.spectral_axis, [3.5, 4.5, 5.5] * u.keV) From 72ed25b346c7848a9e03ecbffd6ec1ad4719e68d Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Fri, 6 Feb 2026 13:17:33 +0000 Subject: [PATCH 06/13] really working for 1d spectra --- examples/spectrum.py | 2 +- sunkit_spex/spectrum/spectrum.py | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/spectrum.py b/examples/spectrum.py index a6a68aad..32bb4439 100644 --- a/examples/spectrum.py +++ b/examples/spectrum.py @@ -28,7 +28,7 @@ # energy. We will start of by creating some synthetic data and corresponding energy bins as well as some important metadata # in this case the exposure time. -data = rng.rand(50) * u.ct +data = rng.random(50) * u.ct energy = np.linspace(1, 50, 50) * u.keV time = Time("2025-02-18T15:08") diff --git a/sunkit_spex/spectrum/spectrum.py b/sunkit_spex/spectrum/spectrum.py index 6ed8eee2..bf0dada1 100644 --- a/sunkit_spex/spectrum/spectrum.py +++ b/sunkit_spex/spectrum/spectrum.py @@ -190,11 +190,12 @@ def gwcs_from_array(array, flux_shape, spectral_axis_index=None): input_frame=detector_frame, output_frame=output_frame, ) + tabular_gwcs.bounding_box = None # Store the intended unit from the origin input array # tabular_gwcs._input_unit = orig_array.unit - return tabular_gwcs # noqa: RET504 + return tabular_gwcs class SpectralAxis(SpectralCoord): From 80adc2557e6fa48850e3f7eb987769aa4827efb9 Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Fri, 6 Feb 2026 14:20:23 +0000 Subject: [PATCH 07/13] 2D spectrogram working --- changelog/239.feature.rst | 1 + examples/spectrum.py | 104 +++++++++++++++++++------------------- 2 files changed, 54 insertions(+), 51 deletions(-) create mode 100644 changelog/239.feature.rst diff --git a/changelog/239.feature.rst b/changelog/239.feature.rst new file mode 100644 index 00000000..2e11e9ce --- /dev/null +++ b/changelog/239.feature.rst @@ -0,0 +1 @@ +Add a new `sunkit_spex.spectrum.spectrum.Spectrum` object to hold spectral data. `~`~sunkit_spex.spectrum.spectrum.Spectrum` is based on `NDCube` and butils on it coordinate aware method and metadata handling. diff --git a/examples/spectrum.py b/examples/spectrum.py index 32bb4439..40d0ac14 100644 --- a/examples/spectrum.py +++ b/examples/spectrum.py @@ -12,6 +12,7 @@ import numpy as np from ndcube import NDMeta +from ndcube.extra_coords import QuantityTableCoordinate, TimeTableCoordinate import astropy.units as u from astropy.coordinates import SpectralCoord @@ -51,30 +52,30 @@ # standard sliceing methods: spec_1d_sliced = spec_1d[10:20] -spec_1d_sliced.shape -spec_1d_sliced.axis_world_coords_values() -spec_1d_sliced.meta -spec_1d_sliced.spectral_axis +print(spec_1d_sliced.shape) +print(spec_1d_sliced.axis_world_coords_values()) +print(spec_1d_sliced.meta) +print(spec_1d_sliced.spectral_axis) ##################################################### # # High level coordinate objects such as SkyCoord and SpectralCoord spec_1d_crop = spec_1d.crop(SpectralCoord(10.5, unit=u.keV), SpectralCoord(20, unit=u.keV)) -spec_1d_crop.shape -spec_1d_crop.axis_world_coords_values() -spec_1d_crop.meta -spec_1d_crop.spectral_axis +print(spec_1d_crop.shape) +print(spec_1d_crop.axis_world_coords_values()) +print(spec_1d_crop.meta) +print(spec_1d_crop.spectral_axis) ##################################################### # # And Quantities spec_1d_crop_value = spec_1d.crop_by_values((10.5 * u.keV), (20.5 * u.keV)) -spec_1d_crop_value.shape -spec_1d_crop_value.axis_world_coords_values() -spec_1d_crop_value.meta -spec_1d_crop_value.spectral_axis +print(spec_1d_crop_value.shape) +print(spec_1d_crop_value.axis_world_coords_values()) +print(spec_1d_crop_value.meta) +print(spec_1d_crop_value.spectral_axis) ##################################################### # @@ -84,70 +85,71 @@ # series of spectra as a function of time. Here we will simulate a series of 10 spectra taken over 10 minutes. Again we # begin by creating our synthetic data as before but additionally creating the time variable. -# data = rng.random(10, 50) * u.ct -# energy = np.linspace(1, 50, 51) * u.keV -# times = Time("2025-02-18T15:08") + np.arange(10) * u.min -# exposure_time = np.arange(5, 15) * u.s +data = rng.random((10, 50)) * u.ct +energy = np.linspace(1, 50, 51) * u.keV +times = Time("2025-02-18T15:08") + np.arange(10) * u.min +exposure_time = np.arange(5, 15) * u.s ##################################################### # # We are also going to demonstrate the power of the sliceable metadata, so in this example each of the individual # spectra have different exposure times (this could be another important information regard the observation) -# meta = NDMeta() -# meta.add("exposure_time", exposure_time, axes=(0,)) -# -# time_coord = TimeTableCoordinate(times, names="time", physical_types="time") -# energy_coord = QuantityTableCoordinate(energy, names="energy", physical_types="em.energy") -# wcs = (energy_coord | time_coord).wcs -# -# spec_2d_time_energy = Spectrum(data, spectral_axis=energy, wcs=wcs, spectral_axis_index=1, meta=meta) -# spec_2d_time_energy.axis_world_coords_values() +meta = NDMeta() +meta.add("exposure_time", exposure_time, axes=(0,)) -###################################################### -# -# For this example we're assuming we don't have a full WCS for the data so we utilise `extra_coords` to store the -# temporal axes. +time_coord = TimeTableCoordinate(times, names="time", physical_types="time") +energy_coord = QuantityTableCoordinate(energy, names="energy", physical_types="em.energy") +wcs = (energy_coord & time_coord).wcs -# spec_2d_time_energy.extra_coords.add("time", (0,), times) -# -# spec_2d_time_energy -# spec_2d_time_energy.shape -# spec_2d_time_energy.axis_world_coords_values() -# spec_2d_time_energy.meta +spec_2d_time_energy = Spectrum(data, spectral_axis=energy, wcs=wcs, spectral_axis_index=1, meta=meta) ###################################################### # # Again all standard slicing works -# spec_2d_time_energy[2:5] -# spec_2d_time_energy[:, 10:20] -# spec_2d_time_energy_sliced = spec_2d_time_energy[2:5, 10:20] +spec_2d_time_energy[2:5] +spec_2d_time_energy[:, 10:20] +spec_2d_time_energy_sliced = spec_2d_time_energy[2:5, 10:20] ###################################################### # # We can being to see the usefulness of the sliceable metadata notice how the exposure time entry has been sliced # appropriately -# spec_2d_time_energy_sliced -# spec_2d_time_energy_sliced.shape -# spec_2d_time_energy_sliced.axis_world_coords_values() -# spec_2d_time_energy_sliced.meta +print(spec_2d_time_energy_sliced.shape) +print(spec_2d_time_energy_sliced.axis_world_coords_values()) +print(spec_2d_time_energy_sliced.meta) +print(spec_2d_time_energy_sliced.spectral_axis) + +###################################################### +# +# The same can be archived using height level coordinate objects +# + +spec_2d_time_energy_crop = spec_2d_time_energy.crop( + [SpectralCoord(10, unit=u.keV), Time("2025-02-18T15:10")], [SpectralCoord(20, unit=u.keV), Time("2025-02-18T15:12")] +) + +print(spec_2d_time_energy_crop.shape) +print(spec_2d_time_energy_crop.axis_world_coords_values()) +print(spec_2d_time_energy_crop.meta) +print(spec_2d_time_energy_crop.spectral_axis) ###################################################### # -# The same can be archived using coordinate objects or values +# Or Quantities as before +spec_2d_time_energy_crop_values = spec_2d_time_energy.crop_by_values((10 * u.keV, 2 * u.min), (19.5 * u.keV, 4 * u.min)) -# spec_2d_time_energy_crop = spec_2d_time_energy.crop( -# [Time("2025-02-18T15:09"), SpectralCoord(10, unit=u.keV)], -# [Time("2025-02-18T15:014"), SpectralCoord(20, unit=u.keV)], -# ) -# spec_2d_time_energy_crop.crop_by_values() +print(spec_2d_time_energy_crop_values.shape) +print(spec_2d_time_energy_crop_values.axis_world_coords_values()) +print(spec_2d_time_energy_crop_values.meta) +print(spec_2d_time_energy_crop_values.spectral_axis) ##################################################### # # 2D Spectrum ( e.g. detector v energy) -# +# ------------------------------------- # data = rng.rand(10, 50) * u.ct # energy = np.linspace(1, 50, 50) * u.keV @@ -165,7 +167,7 @@ ##################################################### # # 3D Spectrum ( e.g. detector v energy v time) -# +# -------------------------------------------- # data = rng.random(10, 20, 30) * u.ct # energy = np.linspace(1, 31, 31) * u.keV @@ -187,7 +189,7 @@ ##################################################### # # 4D Spectrum ( e.g. spatial v spatial v energy v time) -# +# ----------------------------------------------------- # import numpy as np # from ndcube import NDMeta From b583a8bbe2d3a7244665a60f002c1fe71edc882a Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Fri, 6 Feb 2026 19:43:02 +0000 Subject: [PATCH 08/13] Fix GWCS version and start on 2d detector v energy example --- examples/spectrum.py | 31 +++++++++++++++++++++---------- pyproject.toml | 2 +- 2 files changed, 22 insertions(+), 11 deletions(-) diff --git a/examples/spectrum.py b/examples/spectrum.py index 40d0ac14..1c0eb5de 100644 --- a/examples/spectrum.py +++ b/examples/spectrum.py @@ -151,18 +151,29 @@ # 2D Spectrum ( e.g. detector v energy) # ------------------------------------- -# data = rng.rand(10, 50) * u.ct -# energy = np.linspace(1, 50, 50) * u.keV -# -# exposure_time = np.arange(10) * u.s -# labels = np.array([chr(97 + i) for i in range(10)]) +data = rng.random((10, 50)) * u.ct +energy = np.linspace(1, 50, 50) * u.keV + +exposure_time = np.arange(10) * u.s +labels = np.array([f"det_+{chr(97 + i)}" for i in range(10)]) + +meta = NDMeta() +meta.add("exposure_time", exposure_time, axes=0) +meta.add("detector", labels, axes=0) + +spec_2d_det_time = Spectrum(data, spectral_axis=energy, spectral_axis_index=1, meta=meta) +spec_2d_det_time + + +##################################################### # -# meta = NDMeta() -# meta.add("exposure_time", exposure_time, axes=0) -# meta.add("detector", labels, axes=0) + +spec_2d_det_time.crop((SpectralCoord(10 * u.keV), None), (SpectralCoord(20 * u.keV), None)) + +##################################################### # -# spec_2d_det_time = Spectrum(data, spectral_axis=energy, spectral_axis_index=1, meta=meta) -# spec_2d_det_time + +spec_2d_det_time.crop_by_values((10 * u.keV, 0), (20 * u.keV, 2)) ##################################################### # diff --git a/pyproject.toml b/pyproject.toml index 657b6020..9ce6a48a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,7 +26,7 @@ dependencies = [ "scipy>=1.12", "sunpy>=7.0", "xarray>=2023.12", - "gwcs<=1.0.0", #until https://github.com/sunpy/ndcube/issues/913 is resolved + "gwcs<1.0.0", #until https://github.com/sunpy/ndcube/issues/913 is resolved "ndcube>=2.3", ] From c8d929a5785c43df936964fe4aac4eba34355c75 Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Fri, 13 Feb 2026 17:13:11 +0000 Subject: [PATCH 09/13] Fix tests --- sunkit_spex/spectrum/spectrum.py | 7 +------ sunkit_spex/spectrum/tests/test_spectrum.py | 19 +++---------------- 2 files changed, 4 insertions(+), 22 deletions(-) diff --git a/sunkit_spex/spectrum/spectrum.py b/sunkit_spex/spectrum/spectrum.py index bf0dada1..9ee36569 100644 --- a/sunkit_spex/spectrum/spectrum.py +++ b/sunkit_spex/spectrum/spectrum.py @@ -177,12 +177,7 @@ def gwcs_from_array(array, flux_shape, spectral_axis_index=None): # If our spectral axis is in descending order, we have to flip the lookup # table to be ascending in order for world_to_pixel to work. - if len(array) == 0 or array[-1] > array[0]: - forward_transform.inverse = SpectralTabular1D(array, lookup_table=np.arange(len(array)) * u.pix) - else: - raise ValueError("Unsupported ") - # forward_transform.inverse = SpectralTabular1D( - # array[::-1], lookup_table=np.arange(len(array))[::-1]) + forward_transform.inverse = SpectralTabular1D(array, lookup_table=np.arange(len(array)) * u.pix) tabular_gwcs = SpectralGWCS( original_unit=orig_array.unit, diff --git a/sunkit_spex/spectrum/tests/test_spectrum.py b/sunkit_spex/spectrum/tests/test_spectrum.py index 81d89100..3fec06a1 100644 --- a/sunkit_spex/spectrum/tests/test_spectrum.py +++ b/sunkit_spex/spectrum/tests/test_spectrum.py @@ -56,11 +56,11 @@ def test_gwcs_from_array_1d_wavelength(): assert wcs.output_frame.axes_names[0] == "wavelength" # Test forward transform (pixel to world) - assert np.allclose(wcs(0), 4000) - assert np.allclose(wcs(99), 7000) + assert np.allclose(wcs(0), 4000 << u.AA) + assert np.allclose(wcs(99), 7000 << u.AA) # Test inverse transform (world to pixel) - assert np.allclose(wcs.invert(4000), 0) + assert np.allclose(wcs.invert(4000).value, 0) def test_gwcs_from_array_3d_cube(): @@ -84,19 +84,6 @@ def test_gwcs_from_array_3d_cube(): assert wcs.output_frame.frames[1].axes_names[0] == "energy" -def test_gwcs_from_array_descending(): - # Test descending spectral axis for inverse transform logic - waves = np.array([5000, 4000, 3000]) * u.AA - wcs = gwcs_from_array(waves, (3,)) - - assert np.allclose(wcs(0), 5000) - assert np.allclose(wcs(2), 3000) - - # Check that inverse works correctly on flipped table - assert np.allclose(wcs.invert(3000), 2) - assert np.allclose(wcs.invert(5000), 0) - - def test_gwcs_from_array_invalid_units(): data = np.arange(10) * u.Jy # Flux units are not valid for spectral axis with pytest.raises(ValueError, match="Spectral axis units must be one of"): From deb49715f4dae907f079142e6b830de2541930ad Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Tue, 17 Feb 2026 14:00:21 +0000 Subject: [PATCH 10/13] Fix changelog --- changelog/239.feature.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/changelog/239.feature.rst b/changelog/239.feature.rst index 2e11e9ce..72a338c9 100644 --- a/changelog/239.feature.rst +++ b/changelog/239.feature.rst @@ -1 +1 @@ -Add a new `sunkit_spex.spectrum.spectrum.Spectrum` object to hold spectral data. `~`~sunkit_spex.spectrum.spectrum.Spectrum` is based on `NDCube` and butils on it coordinate aware method and metadata handling. +Add a new `sunkit_spex.spectrum.spectrum.Spectrum` object to hold spectral data. `~sunkit_spex.spectrum.spectrum.Spectrum` is based on `NDCube` and butils on it coordinate aware methods and metadata handling. From bab89258925122d98ac5a30745850c3af7ef6704 Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Tue, 17 Feb 2026 14:25:03 +0000 Subject: [PATCH 11/13] Futz with gwcs version pins --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 9ce6a48a..ce3dacd1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,7 +26,7 @@ dependencies = [ "scipy>=1.12", "sunpy>=7.0", "xarray>=2023.12", - "gwcs<1.0.0", #until https://github.com/sunpy/ndcube/issues/913 is resolved + "gwcs>=0.25.0,<1.0.0", #until https://github.com/sunpy/ndcube/issues/913 is resolved "ndcube>=2.3", ] From 562c16547e07b8033277f3de3dc4b49bd2ac9876 Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Tue, 17 Feb 2026 15:05:01 +0000 Subject: [PATCH 12/13] More futzing no idea --- pyproject.toml | 4 ++-- pytest.ini | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index ce3dacd1..e8f44c46 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,10 +23,10 @@ dependencies = [ "numdifftools>=0.9.42", "numpy>=1.26", # Note: keeping support for numpy 1.x for now "parfive>=2.1", - "scipy>=1.12", + "scipy>=1.14.1", "sunpy>=7.0", "xarray>=2023.12", - "gwcs>=0.25.0,<1.0.0", #until https://github.com/sunpy/ndcube/issues/913 is resolved + "gwcs>=0.26.0,<1.0.0", #until https://github.com/sunpy/ndcube/issues/913 is resolved "ndcube>=2.3", ] diff --git a/pytest.ini b/pytest.ini index 433540ac..4f886703 100644 --- a/pytest.ini +++ b/pytest.ini @@ -43,3 +43,4 @@ filterwarnings = ignore::FutureWarning:arviz.* ignore:The isiterable function.*:astropy.utils.exceptions.AstropyDeprecationWarning ignore:'datfix' made the change:astropy.wcs.wcs.FITSFixedWarning + ignore::pyparsing.warnings.PyparsingDeprecationWarning From 6b12c4ed2b8d4ce45bb64ab9556d8b6064c9509c Mon Sep 17 00:00:00 2001 From: Shane Maloney Date: Tue, 17 Feb 2026 15:56:35 +0000 Subject: [PATCH 13/13] Example doc build --- examples/spectrum.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/spectrum.py b/examples/spectrum.py index 1c0eb5de..413a0970 100644 --- a/examples/spectrum.py +++ b/examples/spectrum.py @@ -168,12 +168,12 @@ ##################################################### # -spec_2d_det_time.crop((SpectralCoord(10 * u.keV), None), (SpectralCoord(20 * u.keV), None)) +# spec_2d_det_time.crop((SpectralCoord(10 * u.keV), None), (SpectralCoord(20 * u.keV), None)) ##################################################### # -spec_2d_det_time.crop_by_values((10 * u.keV, 0), (20 * u.keV, 2)) +# spec_2d_det_time.crop_by_values((10 * u.keV, 0), (20 * u.keV, 2)) ##################################################### #