Skip to content

Commit

Permalink
Use exact domain limits for polynomial wavelength solutions in _read_…
Browse files Browse the repository at this point in the history 
…non_linear_iraf_wcs (#1199)

* Enable non-integer domain limits for non_linear_iraf_wcs; improve logging

* Update tests for Chebychev + Legendre wavelength solutions

* Update changelog
  • Loading branch information
@dhomeier
dhomeier authored Dec 9, 2024
1 parent f8a8d6a commit 9ea3ad1
Show file tree
Hide file tree
Showing 3 changed files with 67 additions and 23 deletions.
3 changes: 3 additions & 0 deletions CHANGES.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,9 @@ New Features
Bug Fixes
^^^^^^^^^

- Fixed loaders for IRAF-encoded non-linear Chebychev + Legendre wavelength
solutions that were imposing integer-valued domain limits. [#1199]

Other Changes and Additions
^^^^^^^^^^^^^^^^^^^^^^^^^^^

Expand Down
28 changes: 18 additions & 10 deletions specutils/io/default_loaders/wcs_fits.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -620,11 +620,11 @@ def _read_non_linear_iraf_wcs(header, wcsdim, verbose=False):
'alow': lambda x: int(float(x)), 'ahigh': lambda x: int(float(x)),
'weight': float, 'zeropoint': float,
'ftype': int, 'order': lambda x: int(float(x)),
'pmin': lambda x: int(float(x)), 'pmax': lambda x: int(float(x))}
'pmin': lambda x: float(x), 'pmax': lambda x: float(x)}

ctypen = header['CTYPE{:d}'.format(wcsdim)]
if verbose:
print(f'Attempting to read CTYPE{wcsdim}: {ctypen}')
print(f'Attempting to read CTYPE{wcsdim}: {ctypen} from WAT{wcsdim}*')
if ctypen == 'MULTISPE':
# This is extracting all header cards for f'WAT{wcsdim}_*' into a list
wat_head = header['WAT{:d}*'.format(wcsdim)]
Expand All @@ -640,6 +640,8 @@ def _read_non_linear_iraf_wcs(header, wcsdim, verbose=False):
for key in wat_head:
wat_string += f'{header[key]:68s}' # Keep header from stripping trailing blanks!
wat_array = shlex.split(wat_string.replace('=', ' '))
if verbose:
print(f'Extracted WAT{wcsdim} array[{len(wat_array)}]')
if len(wat_array) % 2 == 0:
for i in range(0, len(wat_array), 2):
# if wat_array[i] not in wcs_dict.keys():
Expand All @@ -655,13 +657,13 @@ def _read_non_linear_iraf_wcs(header, wcsdim, verbose=False):
for n, sp in enumerate(specn):
spec = wat_wcs_dict[sp].split()
wcs_dict = dict((k, wcs_parser[k](spec[i])) for i, k in enumerate(wcs_parser.keys()))
wcs_dict['fpar'] = [float(i) for i in spec[15:]]
wcs_dict['fpar'] = [float(i) for i in spec[15:15+wcs_dict['order']]]

if verbose:
print(f'Retrieving model for {sp}: {wcs_dict["dtype"]} {wcs_dict["ftype"]}')
math_model = _set_math_model(wcs_dict=wcs_dict)
math_model = _set_math_model(wcs_dict=wcs_dict, verbose=verbose)

spectral_axis[n] = math_model(range(1, wcs_dict['pnum'] + 1)) / (1. + wcs_dict['z'])
spectral_axis[n] = math_model(np.arange(wcs_dict['pnum']) + 1) / (1. + wcs_dict['z'])

if verbose:
print(f'Constructed spectral axis of shape {spectral_axis.shape}')
Expand Down Expand Up @@ -710,19 +712,25 @@ def _set_math_model(wcs_dict, verbose=False):
return _none()
elif wcs_dict['dtype'] == 0:
if verbose:
print('Setting model for DTYPE={dtype:d}'.format(**wcs_dict))
print('Setting linear model for DTYPE={dtype:d}'.format(**wcs_dict))
return _linear_solution(wcs_dict=wcs_dict)
elif wcs_dict['dtype'] == 1:
if verbose:
print('Setting model for DTYPE={dtype:d}'.format(**wcs_dict))
print('Setting log-linear model for DTYPE={dtype:d}'.format(**wcs_dict))
return _log_linear(wcs_dict=wcs_dict)
elif wcs_dict['dtype'] == 2:
if verbose:
print('Setting model for DTYPE={dtype:d} FTYPE={ftype:d}'.format(**wcs_dict))
print('Setting non-linear model for DTYPE={dtype:d} FTYPE={ftype:d}'.format(**wcs_dict))
if wcs_dict['ftype'] == 1:
return _chebyshev(wcs_dict=wcs_dict)
model = _chebyshev(wcs_dict=wcs_dict)
if verbose:
print(f'Chebyshev series of degree {model.degree} on domain {model.domain}: {model.parameters}')
return model
elif wcs_dict['ftype'] == 2:
return _non_linear_legendre(wcs_dict=wcs_dict)
model = _non_linear_legendre(wcs_dict=wcs_dict)
if verbose:
print(f'Legendre series of degree {model.degree} on domain {model.domain}: {model.parameters}')
return model
elif wcs_dict['ftype'] == 3:
return _non_linear_cspline(wcs_dict=wcs_dict)
elif wcs_dict['ftype'] == 4:
Expand Down
59 changes: 46 additions & 13 deletions specutils/tests/test_loaders.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -398,7 +398,8 @@ def test_iraf_linear(remote_data_path):
@pytest.mark.filterwarnings('ignore:non-ASCII characters are present in the FITS file header')
@remote_access([{'id': '3359180', 'filename': 'log-linear_fits_solution.fits'}])
def test_iraf_log_linear(remote_data_path):

"""Non-linear wavelength solution for DTYPE=1 (log-linear) encoded IRAF-style (not implemented).
"""
with pytest.raises(NotImplementedError):
assert Spectrum1D.read(remote_data_path, format='iraf')

Expand All @@ -407,43 +408,75 @@ def test_iraf_log_linear(remote_data_path):
@pytest.mark.filterwarnings('ignore:Read spectral axis of shape')
@remote_access([{'id': '3359190', 'filename': 'non-linear_fits_solution_cheb.fits'}])
def test_iraf_non_linear_chebyshev(remote_data_path):
chebyshev_model = models.Chebyshev1D(degree=2, domain=[1616, 3259])
"""Read non-linear wavelength solution for FTYPE=1 (Chebyshev series) encoded IRAF-style.
"""
chebyshev_model = models.Chebyshev1D(degree=2, domain=[1616.37, 3259.98])
chebyshev_model.c0.value = 5115.64008186
chebyshev_model.c1.value = 535.515983712
chebyshev_model.c2.value = -0.779265625182

wavelength_axis = chebyshev_model(range(1, 4097)) * u.angstrom

spectrum_1d = Spectrum1D.read(remote_data_path, format='iraf')

assert isinstance(spectrum_1d, Spectrum1D)
assert_allclose(wavelength_axis, spectrum_1d.wavelength)
assert_allclose(spectrum_1d.wavelength, wavelength_axis, rtol=1e-10)

# Read from HDUList
with fits.open(remote_data_path) as hdulist:
spectrum_1d = Spectrum1D.read(hdulist, format='iraf')
assert isinstance(spectrum_1d, Spectrum1D)
assert_allclose(wavelength_axis, spectrum_1d.wavelength)
assert isinstance(spectrum_1d, Spectrum1D)
assert_allclose(spectrum_1d.wavelength, wavelength_axis, rtol=1e-10)
assert_allclose(spectrum_1d.wavelength[[0, 1, -1]],
[3514.56625403, 3515.2291341, 6190.37186578] * u.angstrom, rtol=1e-10)

# Read pmin, pmax as integer values (standard interpretation pre #1196)
hdulist[0].header['WAT2_002'] = '.39 1. 0. 1 3 1616 3259 5115.64008185559 535.515983711607 -0.7'
chebyshev_model.domain = [1616, 3259]
wavelength_axis = chebyshev_model(range(1, 4097)) * u.angstrom

spectrum_1d = Spectrum1D.read(hdulist, format='iraf')
assert isinstance(spectrum_1d, Spectrum1D)
assert_allclose(spectrum_1d.wavelength, wavelength_axis, rtol=1e-10)
assert_allclose(spectrum_1d.wavelength[[0, 1, -1]],
[3514.41405321, 3515.07718045, 6191.20308524] * u.angstrom, rtol=1e-10)


@pytest.mark.filterwarnings('ignore:Flux unit was not provided')
@pytest.mark.filterwarnings('ignore:non-ASCII characters are present in the FITS file header')
@pytest.mark.filterwarnings('ignore:Read spectral axis of shape')
@remote_access([{'id': '3359194', 'filename': 'non-linear_fits_solution_legendre.fits'}])
def test_iraf_non_linear_legendre(remote_data_path):

legendre_model = models.Legendre1D(degree=3, domain=[21, 4048])
"""Read non-linear wavelength solution for FTYPE=2 (Legendre series) encoded IRAF-style.
"""
legendre_model = models.Legendre1D(degree=3, domain=[21.64, 4048.55])
legendre_model.c0.value = 5468.67555891
legendre_model.c1.value = 835.332144466
legendre_model.c2.value = -6.02202094803
legendre_model.c3.value = -1.13142953897

wavelength_axis = legendre_model(range(1, 4143)) * u.angstrom

spectrum_1d = Spectrum1D.read(remote_data_path, format='iraf')

assert isinstance(spectrum_1d, Spectrum1D)
assert_allclose(wavelength_axis, spectrum_1d.wavelength)
assert_allclose(spectrum_1d.wavelength, wavelength_axis, rtol=1e-10)

# Read from HDUList
with fits.open(remote_data_path) as hdulist:
spectrum_1d = Spectrum1D.read(hdulist, format='iraf')
assert isinstance(spectrum_1d, Spectrum1D)
assert_allclose(spectrum_1d.wavelength, wavelength_axis, rtol=1e-10)
assert_allclose(spectrum_1d.wavelength[[0, 1, -1]],
[4619.77414264, 4620.19462372, 6334.43272858] * u.angstrom, rtol=1e-10)

# Read pmin, pmax as integer values (standard interpretation pre #1196)
hdulist[0].header['WAT2_002'] = '.00 1. 0. 2 4 21 4048.00000000000 5468.67555890614 835.3321444656'
legendre_model.domain = [21, 4048]
wavelength_axis = legendre_model(range(1, 4143)) * u.angstrom

spectrum_1d = Spectrum1D.read(hdulist, format='iraf')
assert isinstance(spectrum_1d, Spectrum1D)
assert_allclose(spectrum_1d.wavelength, wavelength_axis, rtol=1e-10)
assert_allclose(spectrum_1d.wavelength[[0, 1, -1]],
[4620.04343851, 4620.46391004 , 6334.65282602] * u.angstrom, rtol=1e-10)


@pytest.mark.filterwarnings('ignore:non-ASCII characters are present in the FITS file header')
Expand All @@ -466,7 +499,7 @@ def test_iraf_non_linear_cubic_spline(remote_data_path):
@pytest.mark.remote_data
def test_iraf_multispec_chebyshev():
"""Test loading of SpectrumCollection from IRAF MULTISPEC format FITS file -
nonlinear 2D WCS with Chebyshev solution (FTYPE=1).
nonlinear 2D WCS with Chebyshev solution (DTYPE=2, FTYPE=1).
"""
iraf_url = 'https://github.com/astropy/specutils/raw/legacy-specutils/specutils/io/tests/files'
# Read reference ASCII spectrum from remote file.
Expand All @@ -486,7 +519,7 @@ def test_iraf_multispec_chebyshev():
@pytest.mark.remote_data
def test_iraf_multispec_legendre():
"""Test loading of SpectrumCollection from IRAF MULTISPEC format FITS file -
nonlinear 2D WCS with Legendre solution (FTYPE=2).
nonlinear 2D WCS with Legendre solution (DTYPE=2, FTYPE=2).
"""
iraf_url = 'https://github.com/astropy/specutils/raw/legacy-specutils/specutils/io/tests/files'
# Read reference ASCII spectrum from remote file.
Expand Down

0 comments on commit 9ea3ad1

Please sign in to comment.