astropy / specutils

An Astropy coordinated package for astronomical spectroscopy. Maintainers: @rosteen @keflavich @eteq
http://specutils.readthedocs.io/en/latest/
166 stars 127 forks source link

Address wcs1d loader failures on various IRAF multispec formats #1127

Open dhomeier opened 6 months ago

dhomeier commented 6 months ago

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.

codecov[bot] commented 6 months ago

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 ```diff @@ 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.