Skip to content

Commit

Permalink
WIP: "fix" incompatible WCS header from IRAF-style multispec formats
Browse files Browse the repository at this point in the history
  • Loading branch information
dhomeier committed Mar 26, 2024
1 parent 98dfcfd commit ac7b53b
Showing 1 changed file with 30 additions and 17 deletions.
47 changes: 30 additions & 17 deletions specutils/io/default_loaders/wcs_fits.py
Original file line number Diff line number Diff line change
Expand Up @@ -125,7 +125,6 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None,
raise ValueError('No HDU with spectral data found.')

header = hdulist[hdu].header
wcs = WCS(header)

if 'BUNIT' in header:
data = u.Quantity(hdulist[hdu].data, unit=header['BUNIT'])
Expand Down Expand Up @@ -188,23 +187,37 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None,
uunit = uunit**UNCERT_EXP[unc_type.lower()]
uncertainty = UNCERT_REF[unc_type](u.Quantity(uncertainty, unit=uunit))

# Have to translate to WCS-recognised keywords for IRAF-Multispec format:
if header.get('CTYPE1', 'WAVE') == 'MULTISPE':
header['SYSTEM'] = 'MULTISPE'
if 'WAT1_001' in header:

Check warning on line 193 in specutils/io/default_loaders/wcs_fits.py

View check run for this annotation

Codecov / codecov/patch

specutils/io/default_loaders/wcs_fits.py#L192-L193

Added lines #L192 - L193 were not covered by tests
# Try to extract from IRAF-style card or use Angstrom as default.
wat_dict = dict((rec.split('=') for rec in header['WAT1_001'].split()))
unit = wat_dict.get('units', 'Angstrom')
if hasattr(u, unit):
header['CUNIT1'] = unit

Check warning on line 198 in specutils/io/default_loaders/wcs_fits.py

View check run for this annotation

Codecov / codecov/patch

specutils/io/default_loaders/wcs_fits.py#L195-L198

Added lines #L195 - L198 were not covered by tests
else: # try with unit name stripped of excess plural 's'...
header['CUNIT1'] = unit.rstrip('s')
ctype_def = u.Unit(header['CUNIT1']).physical_type
ctype_def = 'WAVE' if 'length' in ctype_def else list(ctype_def)[0]
ctype1 = wat_dict.get('label', ctype_def)
header['CTYPE1'] = ctype1[:4].upper()
if verbose:
print(f"Extracted spectral axis '{header['CTYPE1']}' "

Check warning on line 206 in specutils/io/default_loaders/wcs_fits.py

View check run for this annotation

Codecov / codecov/patch

specutils/io/default_loaders/wcs_fits.py#L200-L206

Added lines #L200 - L206 were not covered by tests
f"with unit '{header['CUNIT1']}' from 'WAT1_001'")
elif header.get('CUNIT1', '') == '':
header['CUNIT1'] = 'Angstrom'
header['CTYPE1'] = 'WAVE'

Check warning on line 210 in specutils/io/default_loaders/wcs_fits.py

View check run for this annotation

Codecov / codecov/patch

specutils/io/default_loaders/wcs_fits.py#L208-L210

Added lines #L208 - L210 were not covered by tests

wcs = WCS(header)
if spectral_axis_unit is not None:
wcs.wcs.cunit[0] = str(spectral_axis_unit)
elif wcs.wcs.cunit[0] == '' and 'WAT1_001' in header:
# Try to extract from IRAF-style card or use Angstrom as default.
wat_dict = dict((rec.split('=') for rec in header['WAT1_001'].split()))
unit = wat_dict.get('units', 'Angstrom')
if hasattr(u, unit):
wcs.wcs.cunit[0] = unit
else: # try with unit name stripped of excess plural 's'...
wcs.wcs.cunit[0] = unit.rstrip('s')
if verbose:
print(f"Extracted spectral axis unit '{unit}' from 'WAT1_001'")
elif wcs.wcs.cunit[0] == '':
wcs.wcs.cunit[0] = 'Angstrom'

# Compatibility attribute for lookup_table (gwcs) WCS
# Compatibility attribute for lookup_table (gwcs) WCS; set physical_type for Spectrum1D.
wcs.unit = tuple(wcs.wcs.cunit)
if verbose:
print(f"WCS spectral axis unit '{wcs.unit}'")
print(f"WCS physical axes '{wcs.world_axis_physical_types}'")

Check warning on line 220 in specutils/io/default_loaders/wcs_fits.py

View check run for this annotation

Codecov / codecov/patch

specutils/io/default_loaders/wcs_fits.py#L219-L220

Added lines #L219 - L220 were not covered by tests

meta = {'header': header}

Expand Down Expand Up @@ -535,12 +548,12 @@ def _read_non_linear_iraf_fits(file_obj, spectral_axis_unit=None, flux_unit=None
wat_dict = dict((rec.split('=') for rec in header['WAT1_001'].split()))
unit = wat_dict.get('units', 'Angstrom')
if hasattr(u, unit):
spectral_axis_unit = unit
spectral_axis_unit = u.Unit(unit)

Check warning on line 551 in specutils/io/default_loaders/wcs_fits.py

View check run for this annotation

Codecov / codecov/patch

specutils/io/default_loaders/wcs_fits.py#L551

Added line #L551 was not covered by tests
else: # try with unit name stripped of excess plural 's'...
spectral_axis_unit = unit.rstrip('s')
spectral_axis_unit = u.Unit(unit.rstrip('s'))

Check warning on line 553 in specutils/io/default_loaders/wcs_fits.py

View check run for this annotation

Codecov / codecov/patch

specutils/io/default_loaders/wcs_fits.py#L553

Added line #L553 was not covered by tests
if verbose:
print(f"Extracted spectral axis unit '{spectral_axis_unit}' from 'WAT1_001'")
spectral_axis *= u.Unit(spectral_axis_unit)
spectral_axis *= spectral_axis_unit

Check warning on line 556 in specutils/io/default_loaders/wcs_fits.py

View check run for this annotation

Codecov / codecov/patch

specutils/io/default_loaders/wcs_fits.py#L556

Added line #L556 was not covered by tests

return spectral_axis, data, dict(header=header)

Expand Down

0 comments on commit ac7b53b

Please sign in to comment.