Allow spectral axis to be anywhere, instead of forcing it to be last (a…

…stropy#1033)

* Starting to work on flexible spectral axis location

Debugging initial spectrum creation

Set private attribute here

Working on debugging failing tests

More things are temporarily broken, but I don't want to lose this work so I'm committing here

Set spectral axis index to 0 if flux is None

Working through test failures

Fix codestyle

Allow passing spectral_axis_index to wcs_fits loader

Require specification of spectral_axis_index if WCS is 1D and flux is multi-D

Decrement spectral_axis_index when slicing with integers

Propagate spectral_axis_index through resampling

Fix last test to account for spectral axis staying first

Fix codestyle

Specify spectral_axis_index in SDSS plate loader

Greatly simply extract_bounding_spectral_region

Account for variable spectral axis location in moment calculation, fix doc example

Working on SpectrumCollection moment handling...not sure this is the way

Need to add one to the axis index here

Update narrative docs to reflect updates

* Add back in the option to move the spectral axis to last, for back-compatibility

Work around pixel unit slicing failure

Change order on crop example

Fix spectral slice handling in tuple input case (e.g. crop)

Update output of crop example

* Apply suggestions from code review

Co-authored-by: Adam Ginsburg <[email protected]>

Apply suggestion from code review

Add helpful comment

* Address review comment about move_spectral_axis, more docs

* Add suggested line to docstring

Co-authored-by: Erik Tollerud <[email protected]>

* Add convenience method

Make this a docstring

* Add v2.0.0 changelog section

---------

Co-authored-by: Erik Tollerud <[email protected]>

CHANGES.rst

@@ -1,4 +1,16 @@
-1.12.0 (unreleased)
+2.0.0 (unreleased)
+------------------
+
+New Features
+^^^^^^^^^^^^
+
+- Spectral axis can now be any axis, rather than being forced to be last. See docs
+  for more details. [#1033]
+
+Other Changes and Additions
+^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+1.11.0 (unreleased)
 -------------------
 
 New Features

docs/index.rst

@@ -23,10 +23,22 @@ details about the underlying principles, see
 `APE13 <https://github.com/astropy/astropy-APEs/blob/main/APE13.rst>`_, the
 guiding document for spectroscopic development in the Astropy Project.
 
-.. note::
-    While specutils is available for general use, the API is in an early enough
-    development stage that some interfaces may change if user feedback and
-    experience warrants it.
+
+Changes in version 2
+====================
+
+Specutils version 2 implemented a major change in that `~specutils.Spectrum1D`
+no longer forces the spectral axis to be last for multi-dimensional data. This
+was motivated by the desire for greater flexibility to allow for interoperability
+with other packages that may wish to use ``specutils`` classes as the basis for
+their own, and by the desire for consistency with the axis order that results
+from a simple ``astropy.io.fits.read`` of a file. The legacy behavior can be
+replicated by setting ``move_spectral_axis='last'`` when creating a new
+`~specutils.Spectrum1D` object.
+
+For a summary of other changes in version 2, please see the
+`release notes <https://github.com/astropy/specutils/releases>`_.
+
 
 Getting started with :ref:`specutils <specutils>`
 =================================================

docs/spectral_cube.rst

@@ -41,14 +41,14 @@ The cube has  74x74 spaxels with 4563 spectral axis points in each one:
 .. code-block:: python
 
     >>> sc.shape # doctest: +REMOTE_DATA
-    (74, 74, 4563)
+    (4563, 74, 74)
 
 
 Print the contents of 3 spectral axis points in a 3x3 spaxel array:
 
 .. code-block:: python
 
-    >>> sc[30:33,30:33,2000:2003] # doctest: +REMOTE_DATA
+    >>> sc[2000:2003,30:33,30:33] # doctest: +REMOTE_DATA
     <Spectrum1D(flux=<Quantity [[[0.48920232, 0.4987253 , 0.5098349 ],
                 [0.493365  , 0.4964812 , 0.5223962 ],
                 [0.49446177, 0.4909543 , 0.5304416 ]],
@@ -144,9 +144,8 @@ Moments
 =======
 
 The `~specutils.analysis.moment` function can be used to compute moments of any order
-along one of the cube's axes. By default, ``axis=-1``, which computes moments
-along the spectral axis (remember that the spectral axis is always last in a
-:class:`~specutils.Spectrum1D`).
+along one of the cube's axes. By default, ``axis='spectral'``, in which case the moment
+is computed along the spectral axis.
 
 .. code-block:: python
 
@@ -193,8 +192,8 @@ cube, using `~specutils.manipulation.spectral_slab` and
 
     # Convert flux density to microJy and correct negative flux offset for
     # this particular dataset
-    ha_flux = (np.sum(subspec.flux.value, axis=(0,1)) + 0.0093) * 1.0E-6*u.Jy
-    ha_flux_wide = (np.sum(subspec_wide.flux.value, axis=(0,1)) + 0.0093) * 1.0E-6*u.Jy
+    ha_flux = (np.sum(subspec.flux.value, axis=(1,2)) + 0.0093) * 1.0E-6*u.Jy
+    ha_flux_wide = (np.sum(subspec_wide.flux.value, axis=(1,2)) + 0.0093) * 1.0E-6*u.Jy
 
     # Compute moment maps for H-alpha line
     moment0_halpha = moment(subspec, order=0)

docs/spectrum1d.rst

@@ -231,19 +231,32 @@ Providing a FITS-style WCS
     >>> spec.wcs.pixel_to_world(np.arange(3))  # doctest: +FLOAT_CMP
     <SpectralCoord [6.5388e-07, 6.5398e-07, 6.5408e-07] m>
 
+When creating a `~specutils.Spectrum1D` using a WCS, you can also use the
+``move_spectral_axis`` argument to force the spectral axis to a certain dimension
+of a multi-dimenasional flux array. Prior to ``specutils`` version 2.0, the flux
+array was always reordered such that the spectral axis corresponded to the last
+flux axis - this behavior can be reproduced by setting ``move_spectral_axis=-1``
+or ``move_spectral_axis='last'``. Note that the relevant axes in the flux, mask,
+and uncertainty arrays are simply swapped, and the swap is also reflected in the
+resulting WCS. No check is currently done to ensure that the resulting array has
+the spatial axes (most often RA and Dec) in any particular order.
 
 Multi-dimensional Data Sets
 ---------------------------
 
 `~specutils.Spectrum1D` also supports the multidimensional case where you
-have, say, an ``(n_spectra, n_pix)``
+have, for example, an ``(n_spectra, n_pix)``
 shaped data set where each ``n_spectra`` element provides a different flux
-data array and so ``flux`` and ``uncertainty`` may be multidimensional as
-long as the last dimension matches the shape of spectral_axis This is meant
+data array. ``flux`` and ``uncertainty`` may be multidimensional as
+long as one dimension matches the shape of the spectral_axis. This is meant
 to allow fast operations on collections of spectra that share the same
 ``spectral_axis``. While it may seem to conflict with the “1D” in the class
 name, this name scheme is meant to communicate the presence of a single
-common spectral axis.
+common spectral axis. In cases where the flux axis corresponding to the spectral
+axis cannot be determined automatically (for example, if multiple flux axes
+have the same length as the spectral axis), the spectral axis must be specified
+with the ``spectral_axis_index`` argument when initializing the
+`~specutils.Spectrum1D`.
 
 .. note:: The case where each flux data array is related to a *different* spectral
           axis is encapsulated in the :class:`~specutils.SpectrumCollection`
@@ -263,8 +276,7 @@ common spectral axis.
                0.33281393, 0.59830875, 0.18673419, 0.67275604, 0.94180287] Jy>
 
 While the above example only shows two dimensions, this concept generalizes to
-any number of dimensions for `~specutils.Spectrum1D`, as long as the spectral
-axis is always the last.
+any number of dimensions for `~specutils.Spectrum1D`.
 
 
 Slicing
@@ -323,72 +335,23 @@ value will apply to the lower bound input.
     ...          'SPECSYS': 'BARYCENT', 'RADESYS': 'ICRS', 'EQUINOX': 2000.0,
     ...          'LONPOLE': 180.0, 'LATPOLE': 27.004754})
     >>> spec = Spectrum1D(flux=np.random.default_rng(12345).random((20, 5, 10)) * u.Jy, wcs=w)  # doctest: +IGNORE_WARNINGS
-    >>> lower = [SpectralCoord(4.9, unit=u.um), SkyCoord(ra=205, dec=26, unit=u.deg)]
-    >>> upper = [SpectralCoord(4.9, unit=u.um), SkyCoord(ra=205.5, dec=27.5, unit=u.deg)]
+    >>> lower = [SkyCoord(ra=205, dec=26, unit=u.deg), SpectralCoord(4.9, unit=u.um)]
+    >>> upper = [SkyCoord(ra=205.5, dec=27.5, unit=u.deg), SpectralCoord(4.9, unit=u.um)]
     >>> spec.crop(lower, upper)  # doctest: +IGNORE_WARNINGS +FLOAT_CMP
-    <Spectrum1D(flux=<Quantity [[[0.70861236],
-                [0.5663815 ],
-                [0.0606386 ],
-                [0.13811995],
-                [0.8974065 ]],
-    <BLANKLINE>
-               [[0.97618597],
-                [0.870499  ],
-                 [0.01522275],
-                 [0.59180312],
-                 [0.29160346]],
-    <BLANKLINE>
-               [[0.13274479],
-                [0.99381789],
-                [0.767089  ],
-                [0.73765335],
-                [0.36185756]],
-    <BLANKLINE>
-               [[0.09635759],
-                [0.38060021],
-                [0.63263223],
-                [0.60785285],
-                [0.4914904 ]],
-    <BLANKLINE>
-               [[0.83173506],
-                [0.1679683 ],
-                [0.39415721],
-                [0.25540459],
-                [0.39779061]],
-    <BLANKLINE>
-               [[0.27625437],
-                [0.90381118],
-                [0.65453332],
-                [0.64141404],
-                [0.36941242]],
-    <BLANKLINE>
-               [[0.90620901],
-                [0.41437874],
-                [0.71279457],
-                [0.41105785],
-                [0.7459662 ]],
-    <BLANKLINE>
-               [[0.57160214],
-                [0.45227951],
-                [0.29112881],
-                [0.44654224],
-                [0.65356417]],
-    <BLANKLINE>
-               [[0.60171961],
-                [0.46916617],
-                [0.3919754 ],
-                [0.01304226],
-                [0.32085301]],
-    <BLANKLINE>
-               [[0.00470581],
-                [0.54743546],
-                [0.24740782],
-                [0.77903598],
-                [0.63457146]]] Jy>, spectral_axis=<SpectralAxis
-       (observer to target:
-          radial_velocity=0.0 km / s
-          redshift=0.0)
-      [4.90049987e-06] m>)>
+    <Spectrum1D(flux=<Quantity [[[0.70861236, 0.97618597, 0.13274479, 0.09635759, 0.83173506,
+             0.27625437, 0.90620901, 0.57160214, 0.60171961, 0.00470581],
+            [0.5663815 , 0.870499  , 0.99381789, 0.38060021, 0.1679683 ,
+             0.90381118, 0.41437874, 0.45227951, 0.46916617, 0.54743546],
+            [0.0606386 , 0.01522275, 0.767089  , 0.63263223, 0.39415721,
+             0.65453332, 0.71279457, 0.29112881, 0.3919754 , 0.24740782],
+            [0.13811995, 0.59180312, 0.73765335, 0.60785285, 0.25540459,
+             0.64141404, 0.41105785, 0.44654224, 0.01304226, 0.77903598],
+            [0.8974065 , 0.29160346, 0.36185756, 0.4914904 , 0.39779061,
+             0.36941242, 0.7459662 , 0.65356417, 0.32085301, 0.63457146]]] Jy>, spectral_axis=<SpectralAxis
+    (observer to target:
+      radial_velocity=0.0 km / s
+      redshift=0.0)
+    [4.90049987e-06] m>)>
 
 Collapsing
 ----------

specutils/analysis/moment.py

@@ -5,13 +5,14 @@
 
 import numpy as np
 from ..manipulation import extract_region
+from ..spectra import SpectrumCollection
 from .utils import computation_wrapper
 
 
 __all__ = ['moment']
 
 
-def moment(spectrum, regions=None, order=0, axis=-1):
+def moment(spectrum, regions=None, order=0, axis='spectral'):
     """
     Estimate the moment of the spectrum.
 
@@ -28,9 +29,9 @@ def moment(spectrum, regions=None, order=0, axis=-1):
     order : int
         The order of the moment to be calculated. Default=0
 
-    axis : int
-        Axis along which a moment is calculated. Default=-1, computes along
-        the last axis (spectral axis).
+    axis : int, str
+        Axis along which a moment is calculated. Default='spectral', which computes
+        along the spectral axis.
 
 
     Returns
@@ -43,10 +44,22 @@ def moment(spectrum, regions=None, order=0, axis=-1):
                                order=order, axis=axis)
 
 
-def _compute_moment(spectrum, regions=None, order=0, axis=-1):
+def _compute_moment(spectrum, regions=None, order=0, axis='spectral'):
     """
     This is a helper function for the above `moment()` method.
     """
+    if axis == "spectral":
+        if isinstance(spectrum, SpectrumCollection):
+            axes = [spec.spectral_axis_index for spec in spectrum]
+            if not np.all([x==axes[0] for x in axes]):
+                raise ValueError("All spectra in SpectrumCollection must have the same "
+                                 "spectral_axis_index for simultaneous moment calculation.")
+            # SpectumCollection adds a leading axis when it stacks the spectra.
+            axis = axes[0]+1
+
+        else:
+            axis = spectrum.spectral_axis_index
+
     if regions is not None:
         calc_spectrum = extract_region(spectrum, regions)
     else:
@@ -64,9 +77,12 @@ def _compute_moment(spectrum, regions=None, order=0, axis=-1):
         return np.sum(flux, axis=axis)
 
     dispersion = spectral_axis
+    # We now have to account for the spectral axis being anywhere, not always last
     if len(flux.shape) > len(spectral_axis.shape):
-        _shape = flux.shape[:-1] + (1,)
-        dispersion = np.tile(spectral_axis, _shape)
+        for i in range(flux.ndim):
+            if i != calc_spectrum.spectral_axis_index:
+                dispersion = np.expand_dims(dispersion, i)
+                dispersion = np.repeat(dispersion, flux.shape[i], i)
 
     if order == 1:
         return np.sum(flux * dispersion, axis=axis) / np.sum(flux, axis=axis)

specutils/fitting/fitmodels.py

@@ -276,7 +276,7 @@ def fit_lines(spectrum, model, fitter=fitting.LevMarLSQFitter(calc_uncertainties
     exclude_regions : list of `~specutils.SpectralRegion`
         List of regions to exclude in the fitting.
     weights : array-like or 'unc', optional
-        If 'unc', the unceratinties from the spectrum object are used to
+        If 'unc', the uncertainties from the spectrum object are used to
         to calculate the weights. If array-like, represents the weights to
         use in the fitting.  Note that if a mask is present on the spectrum, it
         will be applied to the ``weights`` as it would be to the spectrum

specutils/io/default_loaders/sdss.py

@@ -239,4 +239,5 @@ def spPlate_loader(file_obj, limit=None, **kwargs):
         meta['plugmap'] = Table.read(hdulist[5])[0:limit]
 
     return Spectrum1D(flux=flux, wcs=fixed_wcs,
-                      uncertainty=uncertainty, meta=meta, mask=mask)
+                      uncertainty=uncertainty, meta=meta, mask=mask,
+                      spectral_axis_index=-1)

specutils/io/default_loaders/wcs_fits.py

@@ -105,6 +105,8 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None,
     if verbose:
         print("Spectrum file looks like wcs1d-fits")
 
+    spectral_axis_index = kwargs.get("spectral_axis_index")
+
     with read_fileobj_or_hdulist(file_obj, **kwargs) as hdulist:
         if hdu is None:
             for ext in ('FLUX', 'SCI', 'DATA', 'PRIMARY'):
@@ -212,7 +214,8 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None,
     if wcs.naxis > 4:
         raise ValueError('FITS file input to wcs1d_fits_loader is > 4D')
 
-    return Spectrum1D(flux=data, wcs=wcs, mask=mask, uncertainty=uncertainty, meta=meta)
+    return Spectrum1D(flux=data, wcs=wcs, mask=mask, uncertainty=uncertainty,
+                      meta=meta, spectral_axis_index=spectral_axis_index)
 
 
 @custom_writer("wcs1d-fits")

specutils/manipulation/extract_spectral_region.py

@@ -1,5 +1,3 @@
-import sys
-
 from math import floor, ceil  # faster than int(np.floor/ceil(float))
 
 import numpy as np
@@ -162,7 +160,13 @@ def extract_region(spectrum, region, return_single_spectrum=False):
                                         flux=[]*spectrum.flux.unit)
             extracted_spectrum.append(empty_spectrum)
         else:
-            extracted_spectrum.append(spectrum[..., left_index:right_index])
+            slices = [slice(None),] * len(spectrum.shape)
+            slices[spectrum.spectral_axis_index] = slice(left_index, right_index)
+            if len(slices) == 1:
+                slices = slices[0]
+            else:
+                slices = tuple(slices)
+            extracted_spectrum.append(spectrum[slices])
 
     # If there is only one subregion in the region then we will
     # just return a spectrum.
@@ -270,33 +274,10 @@ def extract_bounding_spectral_region(spectrum, region):
     if len(region) == 1:
         return extract_region(spectrum, region)
 
-    min_left = sys.maxsize
-    max_right = -sys.maxsize - 1
-
     # Look for indices that bound the entire set of sub-regions.
-    index_list = [_subregion_to_edge_pixels(sr, spectrum) for sr in region._subregions]
-
-    for left_index, right_index in index_list:
-        if left_index is not None:
-            min_left = min(left_index, min_left)
-        if right_index is not None:
-            max_right = max(right_index, max_right)
-
-    # If both indices are out of bounds then return an empty spectrum
-    if min_left is None and max_right is None:
-        empty_spectrum = Spectrum1D(spectral_axis=[]*spectrum.spectral_axis.unit,
-                                    flux=[]*spectrum.flux.unit)
-        return empty_spectrum
-    else:
-        # If only one index is out of bounds then set it to
-        # the lower or upper extent
-        if min_left is None:
-            min_left = 0
-
-        if max_right is None:
-            max_right = len(spectrum.spectral_axis)
+    min_list = [min(sr) for sr in region._subregions]
+    max_list = [max(sr) for sr in region._subregions]
 
-        if min_left > max_right:
-            min_left, max_right = max_right, min_left
+    single_region = SpectralRegion(min(min_list), max(max_list))
 
-        return spectrum[..., min_left:max_right]
+    return extract_region(spectrum, single_region)