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

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

@@ -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=None``, 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

@@ -234,14 +234,18 @@ 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`
@@ -261,8 +265,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
@@ -321,72 +324,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=None):
     """
     Estimate the moment of the spectrum.
 
@@ -43,10 +44,21 @@ 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=None):
     """
     This is a helper function for the above `moment()` method.
     """
+    if axis is None:
+        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.")
+            axis = axes[0] + 1
+
+        else:
+            axis = spectrum.spectral_axis_index
+
     if regions is not None:
         calc_spectrum = extract_region(spectrum, regions)
     else:
@@ -64,9 +76,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(len(flux.shape)):
+            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

@@ -231,4 +231,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

@@ -85,6 +85,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:
         header = hdulist[hdu].header
         wcs = WCS(header)
@@ -125,7 +127,7 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None,
             except _wcs.NonseparableSubimageCoordinateSystemError as e:
                 raise ValueError(f'WCS cannot be reduced to 1D: {e!r} {wcs}')
 
-    return Spectrum1D(flux=data, wcs=wcs, meta=meta)
+    return Spectrum1D(flux=data, wcs=wcs, 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
@@ -176,7 +174,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.
@@ -284,33 +288,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)

specutils/manipulation/resample.py

@@ -214,7 +214,8 @@ def resample1d(self, orig_spectrum, fin_spec_axis):
         # geometry based, so won't work to just let quantity math handle it.
         resampled_spectrum = Spectrum1D(flux=out_flux,
                                         spectral_axis=np.array(fin_spec_axis) * orig_spectrum.spectral_axis.unit,
-                                        uncertainty=out_uncertainty)
+                                        uncertainty=out_uncertainty,
+                                        spectral_axis_index = orig_spectrum.spectral_axis_index)
 
         return resampled_spectrum
 
@@ -290,7 +291,8 @@ def resample1d(self, orig_spectrum, fin_spec_axis):
 
         return Spectrum1D(spectral_axis=fin_spec_axis,
                           flux=out_flux,
-                          uncertainty=new_unc)
+                          uncertainty=new_unc,
+                          spectral_axis_index = orig_spectrum.spectral_axis_index)
 
 
 class SplineInterpolatedResampler(ResamplerBase):
@@ -367,4 +369,5 @@ def resample1d(self, orig_spectrum, fin_spec_axis):
 
         return Spectrum1D(spectral_axis=fin_spec_axis,
                           flux=out_flux_val*orig_spectrum.flux.unit,
-                          uncertainty=new_unc)
+                          uncertainty=new_unc,
+                          spectral_axis_index = orig_spectrum.spectral_axis_index)