Cannot read IRAF MULTISPEC with a third dimension

[hamogu]

At this point, specutils reads the MULTISPEC format just fine for many cases. Thanks for the great work everyone!

However, I have "in the wild" (i.e. I don't have access to the pipeline that created it and know next to nothing about that) as spectrum that has 3 dimensions.
Axis 1 is the wavelength, axis 2 the order (as is usual in MULTISPEC, but there is an additional linear axis 3 with just two entries (0 for the science data and 1 for the background). This third dimension has CTYPE3=LINEAR so that the data format is

>>>hdus.info()
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     251   (2036, 83, 2)   float32   

It turns out that the third, linear, dimension is passed to specutils.io.default_loaders.wcs_fits._read_linear_iraf_wcs. That routine makes two crucial assumptions that break down here:

1. I assumes to find a header keyword "DC-FLAG", which distinguishes between linear and log-linear scales. Since we don't have a 1D WCS though, but a 3D one, that information is instead part of the WAT3_001 keyword for the third dimension of the MULTISPEC. When I patch that by changing the header and putting a DC-FLAG keyword, then I get further, but still fail in _read_linear_iraf_wcs because
2. _read_linear_iraf_wcs is hardcoded to use the first entry of all WCS information (under the assumptions that it sees a 1D dataset), e.g. 'pnum': wcs._naxis[0] and it does not account for the fact that wcs.wcs.cd is a matrix (in this case [3,3]) and not a simple vector thus 'cdelt': wcs.wcs.cd[0] gives a vector, not a number, and that then confuses the math model that's supposed to create the right function to generate the axis values.
3. I also get a warning ("UserWarning: linear Solution: Try using format='wcs1d-fits' instead") which is inappropriate here, since only one out of three axes is linear and this is clearly not a wcs1d-fits format.

Error message and an example file are attached below.
As a solution to this we could make _read_linear_iraf_wcs take a keyword for the dimension it's looking at, but it might be easier to make a separate function _read_linear_iraf_wcs_as_part_of_a_multispec since the function is only 8 lines long anyway.

However, I'm not sure if we should fix this at all or if this format is just too rare to encounter to worry about. Some thoughts:

• I think I can work around this problem by calling the appropriate _read_*_iraf_wcs by hand and making my own wavelength scale. However, very few users are as deep in the IRAF MULTISPEC (I am the author of Implement reader for IRAF WCS = multispec files, superceeds #14 #18 and Bug in wcs1d-fits reader and flux_unit #1245 after all), so I assume most users would just give up if it doesn't work, especially since the error message is to obscure.
• Currently, SpectrumCollection is essentially a 1D list of 1D spectra. What would I return here? Extend SpectrumCollection to allow more dimensions? Return a list of SpectrumCollection objects? Or allows the spectrum in SpectrumCollection to be a 1D list of 2D spectra? Or maybe it's better to have the user select which index in the third dimension to use with a new keyword in the reader (in my use case, I really only need the extracted spectra ([:,:, 0]), not the backgrounds ([:, :, 1])?

I we can't decide the last question than we can't fix it...

Error message:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[129], [line 1](vscode-notebook-cell:?execution_count=129&line=1)
----> [1](vscode-notebook-cell:?execution_count=129&line=1) ham2011 = SpectrumCollection.read(hdus)

File ~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/nddata/mixins/ndio.py:59, in NDDataRead.__call__(self, *args, **kwargs)
     58 def __call__(self, *args, **kwargs):
---> [59](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/nddata/mixins/ndio.py:59)     return self.registry.read(self._cls, *args, **kwargs)

File ~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/io/registry/core.py:221, in UnifiedInputRegistry.read(self, cls, format, cache, *args, **kwargs)
    218         kwargs.update({"filename": path})
    220 reader = self.get_reader(format, cls)
--> [221](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/io/registry/core.py:221) data = reader(*args, **kwargs)
    223 if not isinstance(data, cls):
    224     # User has read with a subclass where only the parent class is
    225     # registered.  This returns the parent class, so try coercing
    226     # to desired subclass.
    227     try:

File ~/code/specutils/specutils/io/default_loaders/wcs_fits.py:456, in non_linear_multispec_fits(file_obj, **kwargs)
    427 @data_loader('iraf', identifier=identify_iraf_wcs, dtype=SpectrumCollection, extensions=['fits'])
    428 def non_linear_multispec_fits(file_obj, **kwargs):
    429     """Load SpectrumCollection with WCS spectral axes from FITS files written by IRAF
    430 
    431     Loader for files containing 2D spectra, i.e. flux rows of equal length but
   (...)    454     :class:`~specutils.SpectrumCollection`
    455     """
--> [456](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/code/specutils/specutils/io/default_loaders/wcs_fits.py:456)     spectral_axis, flux, meta = _read_non_linear_iraf_fits(file_obj, **kwargs)
    458     if spectral_axis.ndim == 1:
    459         warnings.warn(f'Read 1D spectral axis of length {spectral_axis.shape[0]} - '
    460                       'consider loading into Spectrum.')

File ~/code/specutils/specutils/io/default_loaders/wcs_fits.py:521, in _read_non_linear_iraf_fits(file_obj, spectral_axis_unit, flux_unit, verbose, **kwargs)
    519     warnings.warn("linear Solution: Try using `format='wcs1d-fits'` instead")
    520     wcs = WCS(header)
--> [521](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/code/specutils/specutils/io/default_loaders/wcs_fits.py:521)     spectral_axis = _read_linear_iraf_wcs(wcs=wcs, dc_flag=header['DC-FLAG'])
    522 elif ctypen == 'MULTISPE':
    523     if verbose:

File ~/code/specutils/specutils/io/default_loaders/wcs_fits.py:585, in _read_linear_iraf_wcs(wcs, dc_flag)
    577 wcs_dict = {'crval': wcs.wcs.crval[0],
    578             'crpix': wcs.wcs.crpix[0],
    579             'cdelt': wcs.wcs.cd[0],
    580             'dtype': dc_flag,
    581             'pnum': wcs._naxis[0]}
    583 math_model = _set_math_model(wcs_dict=wcs_dict)
--> [585](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/code/specutils/specutils/io/default_loaders/wcs_fits.py:585) spectral_axis = math_model(range(wcs_dict['pnum']))
    587 return spectral_axis

File ~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:386, in _ModelMeta._handle_special_methods.<locals>.__call__(self, model_set_axis, with_bounding_box, fill_value, equivalencies, inputs_map, *inputs, **new_inputs)
    385 def __call__(self, *inputs, **kwargs):
--> [386](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:386)     """Evaluate this model on the supplied inputs."""
    387     return super(cls, self).__call__(*inputs, **kwargs)

File ~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:387, in _ModelMeta._handle_special_methods.<locals>.__call__(self, *inputs, **kwargs)
    385 def __call__(self, *inputs, **kwargs):
    386     """Evaluate this model on the supplied inputs."""
--> [387](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:387)     return super(cls, self).__call__(*inputs, **kwargs)

File ~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:1092, in Model.__call__(self, *args, **kwargs)
   1089 fill_value = kwargs.pop("fill_value", np.nan)
   1091 # prepare for model evaluation (overridden in CompoundModel)
-> [1092](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:1092) evaluate, inputs, broadcasted_shapes, kwargs = self._pre_evaluate(
   1093     *args, **kwargs
   1094 )
   1096 outputs = self._generic_evaluate(evaluate, inputs, fill_value, with_bbox)
   1098 # post-process evaluation results (overridden in CompoundModel)

File ~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:938, in Model._pre_evaluate(self, *args, **kwargs)
    934 """
    935 Model specific input setup that needs to occur prior to model evaluation.
    936 """
    937 # Broadcast inputs into common size
--> [938](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:938) inputs, broadcasted_shapes = self.prepare_inputs(*args, **kwargs)
    940 # Setup actual model evaluation method
    941 parameters = self._param_sets(raw=True, units=self._has_units)

File ~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:2122, in Model.prepare_inputs(self, model_set_axis, equivalencies, *inputs, **kwargs)
   2118 # The input formatting required for single models versus a multiple
   2119 # model set are different enough that they've been split into separate
   2120 # subroutines
   2121 if self._n_models == 1:
-> [2122](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:2122)     return self._prepare_inputs_single_model(params, inputs, **kwargs)
   2123 else:
   2124     return self._prepare_inputs_model_set(
   2125         params, inputs, model_set_axis, **kwargs
   2126     )

File ~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:1982, in Model._prepare_inputs_single_model(self, params, inputs, **kwargs)
   1976 except ValueError as exc:
   1977     exc.add_note(
   1978         f"self input argument {self.inputs[idx]!r} of shape"
   1979         f" {input_shape!r} cannot be broadcast with parameter"
   1980         f" {param.name!r} of shape {param.shape!r}.",
   1981     )
-> [1982](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:1982)     raise exc
   1984 if len(broadcast) > len(max_broadcast):
   1985     max_broadcast = broadcast

File ~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:1973, in Model._prepare_inputs_single_model(self, params, inputs, **kwargs)
   1969 try:
   1970     # bypass the broadcast_shapes call for performance reasons
   1971     # if parameter is a scalar
   1972     if self.standard_broadcasting and param.shape:
-> [1973](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/astropy/modeling/core.py:1973)         broadcast = np.broadcast_shapes(input_shape, param.shape)
   1974     else:
   1975         broadcast = input_shape

File ~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/numpy/lib/stride_tricks.py:473, in broadcast_shapes(*args)
    435 """
    436 Broadcast the input shapes into a single shape.
    437 
   (...)    470 (5, 6, 7)
    471 """
    472 arrays = [np.empty(x, dtype=[]) for x in args]
--> [473](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/numpy/lib/stride_tricks.py:473) return _broadcast_shape(*arrays)

File ~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/numpy/lib/stride_tricks.py:422, in _broadcast_shape(*args)
    417 """Returns the shape of the arrays that would result from broadcasting the
    418 supplied arrays against each other.
    419 """
    420 # use the old-iterator because np.nditer does not handle size 0 arrays
    421 # consistently
--> [422](https://file+.vscode-resource.vscode-cdn.net/Users/guenther/projects/Chandraprojects/TYC-4144-329-2/~/miniforge3/envs/ciao-4.17/lib/python3.11/site-packages/numpy/lib/stride_tricks.py:422) b = np.broadcast(*args[:32])
    423 # unfortunately, it cannot handle 32 or more arguments directly
    424 for pos in range(32, len(args), 31):
    425     # ironically, np.broadcast does not properly handle np.broadcast
    426     # objects (it treats them as scalars)
    427     # use broadcasting to avoid allocating the full array

ValueError: shape mismatch: objects cannot be broadcast to a single shape.  Mismatch is between arg 0 with shape (2036,) and arg 1 with shape (3,).
self input argument 'x' of shape (2036,) cannot be broadcast with parameter 'slope' of shape (3,).

edited to add more detail on warnings and error messages

[hamogu]

d155.wav.fits.gz

[hamogu]

Here is a full example to reproduce with the unzipped version of the fits file that I attached to the previous comment:

from specutils import Spectrum, SpectrumCollection
from astropy.io import fits

hdus = fits.open('data/spectra/Hamilton/d155.wav.fits')

# Work around the fact that _read_linear_iraf_wcs expects a DC-FLAR keyword and not a WAT3_001 keyword
hdus[0].header['DC-FLAG'] = 0

ham2011 = SpectrumCollection.read(hdus)

[rosteen]

• Currently, SpectrumCollection is essentially a 1D list of 1D spectra. What would I return here? Extend SpectrumCollection to allow more dimensions?

I haven't had time to look into/think about this yet, but I wanted to at least point out that the spectra in a SpectrumCollection can already be multi-dimensional.