specutils icon indicating copy to clipboard operation
specutils copied to clipboard

Address wcs1d loader failures on various IRAF multispec formats

Open dhomeier opened this issue 1 year ago • 1 comments

Opening this both as a WIP and issue for tracking incompatibilities of the wcs1d-fits loader with the IRAF 'equispec' and 'multispec' (Echelle?) formats. These are largely from reports on the list or even off-list or in the Facebook group etc., so I have yet to find out if I can link to publicly available example files.

The "fix" implemented here so far basically just works around a WCS error, loading a MULTISPEC file without a valid spectral WCS; this seems to have been ignored for some time, but with recent specutils versions results in a

File "~/opt/mambaforge/envs/py312/lib/python3.12/site-packages/specutils/spectra/spectrum1d.py", line 197, in __init__
    raise ValueError("Input WCS must have exactly one axis with "
ValueError: Input WCS must have exactly one axis with spectral units, found 0

The underlying problem is that the files are using a WCS defined as

>>> from astropy.wcs import WCS
>>> wcs_mult = {"CTYPE1": "MULTISPEC", "CTYPE2": "MULTISPEC", "CTYPE3": "LINEAR",
                "CD1_1": 1, "CD2_2": 1, "CD2_3": 1, "LTM1_1": 1, "LTM2_2": 1, "LTM3_3": 1, "WCSDIM": 3,
                "WAT0_001": 'system=multispec', "WAT1_001": 'wtype=multispec label=Wavelength units=angstroms',
                "WAT2_001": 'wtype=multispec spec1 = "1 1 2 1. 0.97876143129469 4062 0. -11.77 10',
                "WAT2_002": '.75 1. 0. 2 5 218.715530395508 3952.29809570313 5358.835101501 -1345',
                "WAT2_003": '.90493401415 22.7260636873979 10.6581837361319 0.218203523536446"',
                "BANDID1": "Optimally extracted spectrum", "BANDID2": "Straight sum of spectrum; no CR cleaning",
                "BANDID3": "Background fit (per cross-dispersed pixel)", "BANDID4": "Sigma per pixel.",
                "DCLOG1": 'REFSPEC1 = master_blue_comp'}
>>> WCS(wcs_mult)
WARNING: FITSFixedWarning: 'celfix' made the change 'Linear transformation matrix is singular'. [astropy.wcs.wcs]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "~/opt/mambaforge/envs/py312/lib/python3.12/site-packages/astropy/wcs/wcs.py", line 600, in __init__
    self.wcs.set()
astropy.wcs._wcs.SingularMatrixError: ERROR 3 in wcsset() at line 2808 of file cextern/wcslib/C/wcs.c:
Linear transformation matrix is singular.
ERROR 3 in linset() at line 718 of file cextern/wcslib/C/lin.c:
PCi_ja matrix is singular.

which to the best of my knowledge is simply not supported by WCSLIB (I think spatial WCS of similar construction are, but I don't see a way to define a valid header for a spectral WCS). The kludge here is catching that error and working around the WCS failure by just constructing a spectral_axis from the pixel scale, but of course without having solved the WCS provides no wavelength calibration – it merely allows to read in the spectrum at all, similar to older specutils versions.

As it does not seem that Astropy WCS will be able to do this, a real fix will probably require solving for the spectral axis in specutils, as is e.g. implemented in nonlinearwave.

dhomeier avatar Mar 26 '24 21:03 dhomeier

Codecov Report

Attention: Patch coverage is 12.50000% with 21 lines in your changes are missing coverage. Please review.

Project coverage is 72.04%. Comparing base (98dfcfd) to head (ac7b53b).

Files Patch % Lines
specutils/io/default_loaders/wcs_fits.py 12.50% 21 Missing :warning:
Additional details and impacted files
@@            Coverage Diff             @@
##             main    #1127      +/-   ##
==========================================
- Coverage   72.21%   72.04%   -0.17%     
==========================================
  Files          62       62              
  Lines        4430     4440      +10     
==========================================
  Hits         3199     3199              
- Misses       1231     1241      +10     

:umbrella: View full report in Codecov by Sentry.
:loudspeaker: Have feedback on the report? Share it here.

codecov[bot] avatar Mar 26 '24 21:03 codecov[bot]