From 8d788d053d18046317dc4ff171a08b5850a07d1e Mon Sep 17 00:00:00 2001 From: Derek Homeier Date: Tue, 14 Feb 2023 03:01:17 +0100 Subject: [PATCH 1/5] WIP: fix wcs1d-fits loader for celestial WCS --- specutils/io/default_loaders/wcs_fits.py | 6 ++- specutils/tests/test_loaders.py | 54 +++++++++++++++++++++++- 2 files changed, 58 insertions(+), 2 deletions(-) diff --git a/specutils/io/default_loaders/wcs_fits.py b/specutils/io/default_loaders/wcs_fits.py index 492eae4c3..af901d42b 100644 --- a/specutils/io/default_loaders/wcs_fits.py +++ b/specutils/io/default_loaders/wcs_fits.py @@ -158,7 +158,11 @@ def wcs1d_fits_writer(spectrum, file_name, hdu=0, update_header=False, **kwargs) # Verify spectral axis constructed from WCS wl = spectrum.spectral_axis - dwl = (wcs.all_pix2world(np.arange(len(wl)), 0) - wl.value) / wl.value + # Not sure why the extra check is necessary for FITS WCS + if hasattr(wcs, 'celestial') and wcs.celestial.naxis > 0: + dwl = (wcs.spectral.all_pix2world(np.arange(len(wl)), 0) - wl.value) / wl.value + else: + dwl = (wcs.all_pix2world(np.arange(len(wl)), 0) - wl.value) / wl.value if np.abs(dwl).max() > 1.e-10: m = np.abs(dwl).argmax() raise ValueError('Relative difference between WCS spectral axis and' diff --git a/specutils/tests/test_loaders.py b/specutils/tests/test_loaders.py index 9901e71b3..7b23574b3 100644 --- a/specutils/tests/test_loaders.py +++ b/specutils/tests/test_loaders.py @@ -10,7 +10,7 @@ from astropy.io.fits.verify import VerifyWarning from astropy.table import Table from astropy.units import UnitsWarning -from astropy.wcs import FITSFixedWarning, WCS +from astropy.wcs import FITSFixedWarning, InconsistentAxisTypesError, WCS from astropy.io.registry import IORegistryError from astropy.modeling import models from astropy.tests.helper import quantity_allclose @@ -740,6 +740,58 @@ def test_wcs1d_fits_writer(tmp_path, spectral_axis): assert quantity_allclose(spec.flux, spectrum.flux) +@pytest.mark.parametrize("spectral_axis", + ['WAVE', 'FREQ', 'ENER', 'WAVN']) +def test_wcs1d_fits_cube(tmpdir, spectral_axis): + """Test write/read for Spectrum1D spectral cube with WCS spectral_axis.""" + wlunits = {'WAVE': 'Angstrom', 'FREQ': 'GHz', 'ENER': 'eV', 'WAVN': 'cm**-1'} + # Header dictionary for constructing WCS + hdr = {'CTYPE1': spectral_axis, 'CUNIT1': wlunits[spectral_axis], + 'CRPIX1': 1, 'CRVAL1': 1, 'CDELT1': 0.01, + 'WCSAXES': 3, 'DISPAXIS': 0, + 'CRPIX2': 38.0, 'CRPIX3': 38.0, 'CRVAL2': 205.4384, 'CRVAL3': 27.004754, + 'CTYPE2': 'RA---TAN', 'CTYPE3': 'DEC--TAN', 'CUNIT2': 'deg', 'CUNIT3': 'deg', + 'CDELT2': 3.61111097865634E-05, 'CDELT3': 3.61111097865634E-05, + 'PC1_1': 1.0, 'PC12 ': 0, 'PC1_3': 0, + 'PC2_1': 0, 'PC2_2': -1.0, 'PC2_3': 0, + 'PC3_1': 0, 'PC3_2': 0, 'PC3_3': 1.0} + # Create a small data set + flux = np.arange(1, 121).reshape((10, 4, 3))**2 * 1.e-14 * u.Jy + wlu = u.Unit(hdr['CUNIT1']) + wl0 = hdr['CRVAL1'] + dwl = hdr['CDELT1'] + disp = np.arange(wl0, wl0 + (flux.shape[0] - 0.5) * dwl, dwl) * wlu + + spectrum = Spectrum1D(flux=flux, wcs=WCS(hdr)) + tmpfile = str(tmpdir.join('_tst.fits')) + spectrum.write(tmpfile, hdu=0) + + # Broken reader! + with pytest.raises(InconsistentAxisTypesError, match='Unmatched celestial axes'): + # Read it in and check against the original + spec = Spectrum1D.read(tmpfile, format='wcs1d-fits') + assert spec.flux.unit == spectrum.flux.unit + assert spec.flux.shape == spectrum.flux.shape + assert spec.spectral_axis.unit == spectrum.spectral_axis.unit + assert quantity_allclose(spec.spectral_axis, spectrum.spectral_axis) + assert quantity_allclose(spec.spectral_axis, disp) + assert quantity_allclose(spec.flux, spectrum.flux) + + # Read from HDUList + hdulist = fits.open(tmpfile) + w = WCS(hdulist[0].header) + assert w.naxis == 3 + assert w.axis_type_names == [spectral_axis, 'RA', 'DEC'] + assert w.array_shape == spectrum.flux.shape + + with pytest.raises(InconsistentAxisTypesError, match='Unmatched celestial axes'): + spec = Spectrum1D.read(hdulist, format='wcs1d-fits') + assert isinstance(spec, Spectrum1D) + assert spec.flux.shape == spectrum.flux.shape + assert quantity_allclose(spec.spectral_axis, spectrum.spectral_axis) + assert quantity_allclose(spec.flux, spectrum.flux) + + @pytest.mark.filterwarnings('ignore:Card is too long') @pytest.mark.parametrize("hdu", range(3)) def test_wcs1d_fits_hdus(tmp_path, hdu): From d98d48abdcd0861a3fb1b77f2f327193e91badeb Mon Sep 17 00:00:00 2001 From: Derek Homeier Date: Tue, 14 Feb 2023 03:26:57 +0100 Subject: [PATCH 2/5] DNM: read in full WCS from wcs1d-fits --- specutils/io/default_loaders/wcs_fits.py | 6 - specutils/tests/test_loaders.py | 200 +++++++++++++++++------ 2 files changed, 150 insertions(+), 56 deletions(-) diff --git a/specutils/io/default_loaders/wcs_fits.py b/specutils/io/default_loaders/wcs_fits.py index af901d42b..d1ad7e59c 100644 --- a/specutils/io/default_loaders/wcs_fits.py +++ b/specutils/io/default_loaders/wcs_fits.py @@ -118,12 +118,6 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None, if wcs.naxis > 4: raise ValueError('FITS file input to wcs1d_fits_loader is > 4D') - elif wcs.naxis > 1: - for i in range(wcs.naxis - 1, 0, -1): - try: - wcs = wcs.dropaxis(i) - except _wcs.NonseparableSubimageCoordinateSystemError as e: - raise ValueError(f'WCS cannot be reduced to 1D: {e!r} {wcs}') return Spectrum1D(flux=data, wcs=wcs, meta=meta) diff --git a/specutils/tests/test_loaders.py b/specutils/tests/test_loaders.py index 7b23574b3..d667371bb 100644 --- a/specutils/tests/test_loaders.py +++ b/specutils/tests/test_loaders.py @@ -705,11 +705,13 @@ def test_tabular_fits_compressed(compress, tmp_path): @pytest.mark.parametrize("spectral_axis", - ['wavelength', 'frequency', 'energy', 'wavenumber']) -def test_wcs1d_fits_writer(tmp_path, spectral_axis): + ['WAVE', 'FREQ', 'ENER', 'WAVN']) +@pytest.mark.parametrize("with_mask", [False, True]) +@pytest.mark.parametrize("uncertainty", + [None, StdDevUncertainty, VarianceUncertainty, InverseVariance]) +def test_wcs1d_fits_writer(tmp_path, spectral_axis, with_mask, uncertainty): """Test write/read for Spectrum1D with WCS-constructed spectral_axis.""" - wlunits = {'wavelength': 'Angstrom', 'frequency': 'GHz', 'energy': 'eV', - 'wavenumber': 'cm**-1'} + wlunits = {'WAVE': 'Angstrom', 'FREQ': 'GHz', 'ENER': 'eV', 'WAVN': 'cm**-1'} # Header dictionary for constructing WCS hdr = {'CTYPE1': spectral_axis, 'CUNIT1': wlunits[spectral_axis], 'CRPIX1': 1, 'CRVAL1': 1, 'CDELT1': 0.01} @@ -719,9 +721,20 @@ def test_wcs1d_fits_writer(tmp_path, spectral_axis): wl0 = hdr['CRVAL1'] dwl = hdr['CDELT1'] disp = np.arange(wl0, wl0 + (len(flux) - 0.5) * dwl, dwl) * wlu + tmpfile = tmp_path / 'wcs_tst.fits' - spectrum = Spectrum1D(flux=flux, wcs=WCS(hdr)) - tmpfile = str(tmp_path / '_tst.fits') + if with_mask: + mask = np.array([0, 0, 0, 0, 1, 0, 0, 0, 1, 0], dtype=np.uint16) + else: + mask = None + + # ToDo: test with explicit (and different from flux) units. + if uncertainty is None: + spectrum = Spectrum1D(flux=flux, wcs=WCS(hdr), mask=mask) + assert spectrum.uncertainty is None + else: + unc = uncertainty(0.1 * np.sqrt(np.abs(flux.value))) + spectrum = Spectrum1D(flux=flux, wcs=WCS(hdr), mask=mask, uncertainty=unc) spectrum.write(tmpfile, hdu=0) # Read it in and check against the original @@ -731,18 +744,31 @@ def test_wcs1d_fits_writer(tmp_path, spectral_axis): assert quantity_allclose(spec.spectral_axis, spectrum.spectral_axis) assert quantity_allclose(spec.spectral_axis, disp) assert quantity_allclose(spec.flux, spectrum.flux) + assert np.all(spec.mask == spectrum.mask) + if uncertainty is None: + assert spec.uncertainty is None + else: + assert quantity_allclose(spec.uncertainty.quantity, spectrum.uncertainty.quantity) # Read from HDUList - with fits.open(tmpfile) as hdulist: - spec = Spectrum1D.read(hdulist, format='wcs1d-fits') + hdulist = fits.open(tmpfile) + spec = Spectrum1D.read(hdulist, format='wcs1d-fits') assert isinstance(spec, Spectrum1D) assert quantity_allclose(spec.spectral_axis, spectrum.spectral_axis) assert quantity_allclose(spec.flux, spectrum.flux) + assert np.all(spec.mask == spectrum.mask) + if uncertainty is None: + assert spec.uncertainty is None + else: + assert quantity_allclose(spec.uncertainty.quantity, spectrum.uncertainty.quantity) @pytest.mark.parametrize("spectral_axis", ['WAVE', 'FREQ', 'ENER', 'WAVN']) -def test_wcs1d_fits_cube(tmpdir, spectral_axis): +@pytest.mark.parametrize("with_mask", [False, True]) +@pytest.mark.parametrize("uncertainty", + [None, StdDevUncertainty, VarianceUncertainty, InverseVariance]) +def test_wcs1d_fits_cube(tmp_path, spectral_axis, with_mask, uncertainty): """Test write/read for Spectrum1D spectral cube with WCS spectral_axis.""" wlunits = {'WAVE': 'Angstrom', 'FREQ': 'GHz', 'ENER': 'eV', 'WAVN': 'cm**-1'} # Header dictionary for constructing WCS @@ -760,36 +786,110 @@ def test_wcs1d_fits_cube(tmpdir, spectral_axis): wlu = u.Unit(hdr['CUNIT1']) wl0 = hdr['CRVAL1'] dwl = hdr['CDELT1'] - disp = np.arange(wl0, wl0 + (flux.shape[0] - 0.5) * dwl, dwl) * wlu + disp = np.arange(wl0, wl0 + (flux.shape[2] - 0.5) * dwl, dwl) * wlu + tmpfile = tmp_path / 'wcs_tst.fits' - spectrum = Spectrum1D(flux=flux, wcs=WCS(hdr)) - tmpfile = str(tmpdir.join('_tst.fits')) - spectrum.write(tmpfile, hdu=0) + if with_mask: + die = np.random.Generator(np.random.MT19937(23)) + mask = die.choice([0, 0, 0, 0, 0, 1], size=flux.shape).astype(np.uint16) + else: + mask = None + + if uncertainty is None: + spectrum = Spectrum1D(flux=flux, wcs=WCS(hdr), mask=mask) + assert spectrum.uncertainty is None + with pytest.warns(AstropyUserWarning, match='No uncertainty array found'): + spectrum.write(tmpfile, hdu=0, uncertainty_name='STD') + else: + unc = uncertainty(0.1 * np.sqrt(np.abs(flux.value))) + spectrum = Spectrum1D(flux=flux, wcs=WCS(hdr), mask=mask, uncertainty=unc) + spectrum.write(tmpfile, hdu=0) # Broken reader! - with pytest.raises(InconsistentAxisTypesError, match='Unmatched celestial axes'): - # Read it in and check against the original - spec = Spectrum1D.read(tmpfile, format='wcs1d-fits') - assert spec.flux.unit == spectrum.flux.unit - assert spec.flux.shape == spectrum.flux.shape - assert spec.spectral_axis.unit == spectrum.spectral_axis.unit - assert quantity_allclose(spec.spectral_axis, spectrum.spectral_axis) - assert quantity_allclose(spec.spectral_axis, disp) - assert quantity_allclose(spec.flux, spectrum.flux) + # Read it in and check against the original + spec = Spectrum1D.read(tmpfile, format='wcs1d-fits') + assert spec.flux.unit == spectrum.flux.unit + assert spec.flux.shape == spectrum.flux.shape + assert spec.spectral_axis.unit == spectrum.spectral_axis.unit + assert quantity_allclose(spec.spectral_axis, spectrum.spectral_axis) + assert quantity_allclose(spec.spectral_axis, disp) + assert quantity_allclose(spec.flux, spectrum.flux) + assert np.all(spec.mask == spectrum.mask) + if uncertainty is None: + assert spec.uncertainty is None + else: + assert quantity_allclose(spec.uncertainty.quantity, spectrum.uncertainty.quantity) # Read from HDUList - hdulist = fits.open(tmpfile) - w = WCS(hdulist[0].header) + with fits.open(tmpfile) as hdulist: + w = WCS(hdulist[0].header) + spec = Spectrum1D.read(hdulist, format='wcs1d-fits') + assert w.naxis == 3 assert w.axis_type_names == [spectral_axis, 'RA', 'DEC'] - assert w.array_shape == spectrum.flux.shape + assert isinstance(spec, Spectrum1D) + assert spec.flux.shape == spectrum.flux.shape + assert quantity_allclose(spec.spectral_axis, spectrum.spectral_axis) + assert quantity_allclose(spec.flux, spectrum.flux) + assert np.all(spec.mask == spectrum.mask) + if uncertainty is None: + assert spec.uncertainty is None + else: + assert quantity_allclose(spec.uncertainty.quantity, spectrum.uncertainty.quantity) - with pytest.raises(InconsistentAxisTypesError, match='Unmatched celestial axes'): - spec = Spectrum1D.read(hdulist, format='wcs1d-fits') - assert isinstance(spec, Spectrum1D) - assert spec.flux.shape == spectrum.flux.shape - assert quantity_allclose(spec.spectral_axis, spectrum.spectral_axis) - assert quantity_allclose(spec.flux, spectrum.flux) + +@pytest.mark.parametrize("uncertainty_rsv", ['STD', 'ERR', 'UNCERT', 'VAR', 'IVAR']) +def test_wcs1d_fits_uncertainty(tmp_path, uncertainty_rsv): + """ + Test Spectrum1D.write with custom `uncertainty` names, + ensure it raises on illegal (reserved) names. + """ + # Header dictionary for constructing WCS + hdr = {'CTYPE1': 'WAVE', 'CUNIT1': 'Angstrom', + 'CRPIX1': 1, 'CRVAL1': 1, 'CDELT1': 0.01} + # Reserved EXTNAMEs for uncertainty types + UNCERT_REF = {'STD': StdDevUncertainty, 'ERR': StdDevUncertainty, 'UNCERT': StdDevUncertainty, + 'VAR': VarianceUncertainty, 'IVAR': InverseVariance} + # Alternative EXTNAMEs for uncertainty types + UNCERT_ALT = {'std': 'StdErr', 'var': 'VARIAN', 'ivar': 'InvVar'} + + # Create a small data set + flux = np.arange(1, 11)**2 * 1.e-14 * u.Jy + wlu = u.Unit(hdr['CUNIT1']) + wl0 = hdr['CRVAL1'] + dwl = hdr['CDELT1'] + disp = np.arange(wl0, wl0 + (len(flux) - 0.5) * dwl, dwl) * wlu + tmpfile = tmp_path / 'wcs_tst.fits' + + mask = np.array([0, 0, 0, 0, 1, 0, 0, 0, 1, 0], dtype=np.uint16) + + # Set uncertainty to mismatched type + uncertainty = [u for n, u in UNCERT_REF.items() if u != UNCERT_REF[uncertainty_rsv]][0] + unc = uncertainty(0.1 * np.sqrt(np.abs(flux.value))) + spectrum = Spectrum1D(flux=flux, wcs=WCS(hdr), mask=mask, uncertainty=unc) + + with pytest.raises(ValueError, match=f"Illegal label for uncertainty: '{uncertainty_rsv}' " + f"is reserved for {UNCERT_REF[uncertainty_rsv]}, not {uncertainty}."): + spectrum.write(tmpfile, format='wcs1d-fits', uncertainty_name=uncertainty_rsv) + + uncertainty_type = spectrum.uncertainty.uncertainty_type + uncertainty_alt = UNCERT_ALT[uncertainty_type] + spectrum.write(tmpfile, format='wcs1d-fits', uncertainty_name=uncertainty_alt) + + # Check EXTNAME (uncertainty is in last HDU) + with fits.open(tmpfile) as hdulist: + assert hdulist[-1].name == uncertainty_alt.upper() + + # Read it in and check against the original + with pytest.raises(ValueError, match=f"Invalid uncertainty type: '{uncertainty_alt}'; should"): + spec = Spectrum1D.read(tmpfile, uncertainty_hdu=2, uncertainty_type=uncertainty_alt) + spec = Spectrum1D.read(tmpfile, uncertainty_hdu=2, uncertainty_type=uncertainty_type) + assert spec.flux.unit == spectrum.flux.unit + assert spec.spectral_axis.unit == spectrum.spectral_axis.unit + assert quantity_allclose(spec.uncertainty.quantity, spectrum.uncertainty.quantity) + spec = Spectrum1D.read(tmpfile, uncertainty_hdu=uncertainty_alt, + uncertainty_type=uncertainty_type) + assert quantity_allclose(spec.uncertainty.quantity, spectrum.uncertainty.quantity) @pytest.mark.filterwarnings('ignore:Card is too long') @@ -804,7 +904,7 @@ def test_wcs1d_fits_hdus(tmp_path, hdu): flux = np.arange(1, 11)**2 * 1.e-14 * flu spectrum = Spectrum1D(flux=flux, wcs=WCS(hdr)) - tmpfile = str(tmp_path / '_tst.fits') + tmpfile = tmp_path / 'tst.fits' spectrum.write(tmpfile, hdu=hdu, format='wcs1d-fits') # Read it in and check against the original @@ -816,7 +916,7 @@ def test_wcs1d_fits_hdus(tmp_path, hdu): assert quantity_allclose(hdulist[hdu].data * flu, flux) # Test again with automatic format selection by filename pattern - tmpfile = str(tmp_path / '_wcs.fits') + tmpfile = tmp_path / 'wcs.fits' spectrum.write(tmpfile, hdu=hdu) with fits.open(tmpfile) as hdulist: assert hdulist[hdu].is_image @@ -824,11 +924,10 @@ def test_wcs1d_fits_hdus(tmp_path, hdu): @pytest.mark.parametrize("spectral_axis", - ['wavelength', 'frequency', 'energy', 'wavenumber']) + ['WAVE', 'FREQ', 'ENER', 'WAVN']) def test_wcs1d_fits_multid(tmp_path, spectral_axis): """Test spectrum with WCS-1D spectral_axis and higher dimension in flux.""" - wlunits = {'wavelength': 'Angstrom', 'frequency': 'GHz', 'energy': 'eV', - 'wavenumber': 'cm**-1'} + wlunits = {'WAVE': 'Angstrom', 'FREQ': 'GHz', 'ENER': 'eV', 'WAVN': 'cm**-1'} # Header dictionary for constructing WCS hdr = {'CTYPE1': spectral_axis, 'CUNIT1': wlunits[spectral_axis], 'CRPIX1': 1, 'CRVAL1': 1, 'CDELT1': 0.01} @@ -839,12 +938,12 @@ def test_wcs1d_fits_multid(tmp_path, spectral_axis): dwl = hdr['CDELT1'] disp = np.arange(wl0, wl0 + len(flux[1:]) * dwl, dwl) * wlu - # Construct up to 4D flux array, write and read (no auto-identify) + # Construct 2D to 4D flux array, write and read (no auto-identify) shape = [-1, 1] for i in range(2, 5): flux = flux * np.arange(i, i+5).reshape(*shape) spectrum = Spectrum1D(flux=flux, wcs=WCS(hdr)) - tmpfile = str(tmp_path / f'_{i}d.fits') + tmpfile = tmp_path / f'wcs_{i}d.fits' spectrum.write(tmpfile, format='wcs1d-fits') spec = Spectrum1D.read(tmpfile, format='wcs1d-fits') @@ -857,34 +956,35 @@ def test_wcs1d_fits_multid(tmp_path, spectral_axis): # Test exception for NAXIS > 4 flux = flux * np.arange(i+1, i+6).reshape(*shape) spectrum = Spectrum1D(flux=flux, wcs=WCS(hdr)) - tmpfile = str(tmp_path / f'_{i+1}d.fits') + tmpfile = tmp_path / f'wcs_{i+1}d.fits' spectrum.write(tmpfile, format='wcs1d-fits') with pytest.raises(ValueError, match='input to wcs1d_fits_loader is > 4D'): spec = Spectrum1D.read(tmpfile, format='wcs1d-fits') -@pytest.mark.parametrize("spectral_axis", ['wavelength', 'frequency']) +@pytest.mark.parametrize("spectral_axis", ['WAVE', 'FREQ']) def test_wcs1d_fits_non1d(tmp_path, spectral_axis): """Test exception on trying to load FITS with 2D flux and irreducible WCS spectral_axis. """ - wlunits = {'wavelength': 'Angstrom', 'frequency': 'GHz', 'energy': 'eV', - 'wavenumber': 'cm**-1'} + wlunits = {'WAVE': 'Angstrom', 'FREQ': 'GHz', 'ENER': 'eV', 'WAVN': 'cm**-1'} # Header dictionary for constructing WCS hdr = {'CTYPE1': spectral_axis, 'CUNIT1': wlunits[spectral_axis], 'CRPIX1': 1, 'CRVAL1': 1, 'CDELT1': 0.01} # Create a small 2D data set flux = np.arange(1, 11)**2 * np.arange(4).reshape(-1, 1) * 1.e-14 * u.Jy spectrum = Spectrum1D(flux=flux, wcs=WCS(hdr)) - tmpfile = str(tmp_path / f'_{2}d.fits') + tmpfile = tmp_path / f'wcs_{2}d.fits' spectrum.write(tmpfile, format='wcs1d-fits') # Reopen file and update header with off-diagonal element - with fits.open(tmpfile, mode='update') as hdulist: - hdulist[0].header.update([('PC1_2', 0.2)]) + hdulist = fits.open(tmpfile, mode='update') + hdulist[0].header.update([('PC1_2', 0.2)]) + hdulist.close() - with pytest.raises(ValueError, match='WCS cannot be reduced to 1D'): + with pytest.raises(ValueError, + match='Non-zero off-diagonal matrix elements excluded from the subimage.'): Spectrum1D.read(tmpfile, format='wcs1d-fits') @@ -894,7 +994,7 @@ def test_wcs1d_fits_non1d(tmp_path, spectral_axis): def test_wcs1d_fits_compressed(compress, tmp_path): """Test automatic recognition of supported compression formats for IMAGE/WCS. """ - ext = {'gzip': '.gz', 'bzip2': '.bz2', 'xz': '.xz'} + ext = {'gzip': 'gz', 'bzip2': 'bz2', 'xz': 'xz'} if compress == 'bzip2' and not HAS_BZ2: pytest.xfail("Python installation has no bzip2 support") if compress == 'xz' and not HAS_LZMA: @@ -911,14 +1011,14 @@ def test_wcs1d_fits_compressed(compress, tmp_path): disp = np.arange(wl0, wl0 + (len(flux) - 0.5) * dwl, dwl) * wlu spectrum = Spectrum1D(flux=flux, wcs=WCS(hdr)) - tmpfile = str(tmp_path / '_tst.fits') + tmpfile = tmp_path / 'wcs_tst.fits' spectrum.write(tmpfile, hdu=0) # Deliberately not using standard filename pattern to test header info. with warnings.catch_warnings(): warnings.simplefilter('ignore', FITSFixedWarning) os.system(f'{compress} {tmpfile}') - spec = Spectrum1D.read(tmpfile + ext[compress]) + spec = Spectrum1D.read(tmpfile.with_suffix(f'.fits.{ext[compress]}')) assert isinstance(spec, Spectrum1D) assert quantity_allclose(spec.spectral_axis, disp) @@ -927,7 +1027,7 @@ def test_wcs1d_fits_compressed(compress, tmp_path): # Try again without compression suffix: with warnings.catch_warnings(): warnings.simplefilter('ignore', FITSFixedWarning) - os.system(f'mv {tmpfile}{ext[compress]} {tmpfile}') + os.system(f'mv {tmpfile}.{ext[compress]} {tmpfile}') spec = Spectrum1D.read(tmpfile) assert isinstance(spec, Spectrum1D) From 6f4908d58fba9bed1182bfc169c9749081ee399f Mon Sep 17 00:00:00 2001 From: Derek Homeier Date: Fri, 17 Feb 2023 00:44:36 +0100 Subject: [PATCH 3/5] Support HDU naming, add mask HDU r/w for wcs1d-fits --- specutils/io/default_loaders/wcs_fits.py | 59 ++++++++++++++++++------ 1 file changed, 44 insertions(+), 15 deletions(-) diff --git a/specutils/io/default_loaders/wcs_fits.py b/specutils/io/default_loaders/wcs_fits.py index d1ad7e59c..6dbf32397 100644 --- a/specutils/io/default_loaders/wcs_fits.py +++ b/specutils/io/default_loaders/wcs_fits.py @@ -3,7 +3,7 @@ from astropy import units as u from astropy.io import fits -from astropy.wcs import WCS, _wcs +from astropy.wcs import WCS from astropy.modeling import models from astropy.utils.exceptions import AstropyUserWarning @@ -36,7 +36,8 @@ def identify_wcs1d_fits(origin, *args, **kwargs): with read_fileobj_or_hdulist(*args, **kwargs) as hdulist: return (hdulist[hdu].header.get('WCSDIM', 1) == 1 and (hdulist[hdu].header['NAXIS'] == 1 or - hdulist[hdu].header.get('WCSAXES', 0) == 1 )and not + hdulist[hdu].header.get('WCSAXES', 0) == 1 or + hdulist[hdu].header.get('DISPAXIS', -1) >= 0 ) and not hdulist[hdu].header.get('MSTITLE', 'n').startswith('2dF-SDSS LRG') and not # Check in CTYPE1 key for linear solution (single spectral axis) hdulist[hdu].header.get('CTYPE1', 'w').upper().startswith('MULTISPE')) @@ -45,7 +46,7 @@ def identify_wcs1d_fits(origin, *args, **kwargs): @data_loader("wcs1d-fits", identifier=identify_wcs1d_fits, dtype=Spectrum1D, extensions=['fits', 'fit'], priority=5) def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None, - hdu=0, verbose=False, **kwargs): + hdu=None, verbose=False, **kwargs): """ Loader for single spectrum-per-HDU spectra in FITS files, with the spectral axis stored in the header as FITS-WCS. The flux unit of the spectrum is @@ -67,9 +68,9 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None, Units of the flux for this spectrum. If not given (or None), the unit will be inferred from the BUNIT keyword in the header. Note that this unit will attempt to convert from BUNIT if BUNIT is present. - hdu : int - The index of the HDU to load into this spectrum. - verbose : bool + hdu : int, str or None. optional + The index or name of the HDU to load into this spectrum (default: search 1st). + verbose : bool. optional Print extra info. **kwargs Extra keywords for :func:`~specutils.io.parsing_utils.read_fileobj_or_hdulist`. @@ -86,6 +87,17 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None, print("Spectrum file looks like wcs1d-fits") with read_fileobj_or_hdulist(file_obj, **kwargs) as hdulist: + if hdu is None: + for ext in ('FLUX', 'SCI', 'PRIMARY'): + # For now rely on extension containing spectral data. + if ext in hdulist and ( + isinstance(hdulist[ext], (fits.ImageHDU, fits.PrimaryHDU)) and + hdulist[ext].data is not None): + hdu = ext + break + else: + raise ValueError('No HDU with spectral data found.') + header = hdulist[hdu].header wcs = WCS(header) @@ -96,6 +108,12 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None, else: data = u.Quantity(hdulist[hdu].data, unit=flux_unit) + # TODO: add additional HDU data like 'IVAR' + if 'MASK' in hdulist: + mask = hdulist['MASK'].data + else: + mask = None + 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: @@ -119,11 +137,11 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None, if wcs.naxis > 4: raise ValueError('FITS file input to wcs1d_fits_loader is > 4D') - return Spectrum1D(flux=data, wcs=wcs, meta=meta) + return Spectrum1D(flux=data, wcs=wcs, mask=mask, meta=meta) @custom_writer("wcs1d-fits") -def wcs1d_fits_writer(spectrum, file_name, hdu=0, update_header=False, **kwargs): +def wcs1d_fits_writer(spectrum, file_name, hdu=0, update_header=False, flux_name='FLUX', **kwargs): """ Write spectrum with spectral axis defined by its WCS to (primary) IMAGE_HDU of a FITS file. @@ -133,14 +151,16 @@ def wcs1d_fits_writer(spectrum, file_name, hdu=0, update_header=False, **kwargs) spectrum : :class:`~specutils.Spectrum1D` file_name : str The path to the FITS file - hdu : int + hdu : int, optional Header Data Unit in FITS file to write to (base 0; default primary HDU) - update_header : bool + update_header : bool, optional Update FITS header with all compatible entries in `spectrum.meta` - unit : str or :class:`~astropy.units.Unit` - Unit for the flux (and associated uncertainty) - dtype : str or :class:`~numpy.dtype` - Floating point type for storing flux array + flux_name : str, optional + HDU name to store flux spectrum under (default 'FLUX') + unit : str or :class:`~astropy.units.Unit`, optional + Unit for the flux (and associated uncertainty; defaults to `spectrum.flux.unit`) + dtype : str or :class:`~numpy.dtype`, optional + Floating point type for storing flux array (defaults to `spectrum.flux.dtype`) """ # Create HDU list from WCS try: @@ -160,7 +180,7 @@ def wcs1d_fits_writer(spectrum, file_name, hdu=0, update_header=False, **kwargs) if np.abs(dwl).max() > 1.e-10: m = np.abs(dwl).argmax() raise ValueError('Relative difference between WCS spectral axis and' - f'spectral_axis at {m:}: dwl[m]') + f'spectral_axis at {m:}: {dwl[m]}') if update_header: hdr_types = (str, int, float, complex, bool, @@ -180,6 +200,15 @@ def wcs1d_fits_writer(spectrum, file_name, hdu=0, update_header=False, **kwargs) funit = u.Unit(kwargs.pop('unit', spectrum.flux.unit)) flux = spectrum.flux.to(funit, equivalencies=u.spectral_density(wl)) hdulist[0].data = flux.value.astype(ftype) + if flux_name is not None: + hdulist[0].name = flux_name + + # Append mask array (duplicate WCS for that extension)? + if spectrum.mask is not None: + hdulist.append(wcs.to_fits()[0]) + hdulist[1].data = spectrum.mask + hdulist[1].name = 'MASK' + hdu += 1 if hasattr(funit, 'long_names') and len(funit.long_names) > 0: comment = f'[{funit.long_names[0]}] {funit.physical_type}' From abb8c1eda65bfd6133846871da8b86720787533a Mon Sep 17 00:00:00 2001 From: Derek Homeier Date: Wed, 15 Mar 2023 19:24:26 +0100 Subject: [PATCH 4/5] Handle custom uncertainty names and test them --- specutils/io/default_loaders/wcs_fits.py | 182 ++++++++++++++++++----- specutils/tests/test_loaders.py | 23 ++- 2 files changed, 152 insertions(+), 53 deletions(-) diff --git a/specutils/io/default_loaders/wcs_fits.py b/specutils/io/default_loaders/wcs_fits.py index 6dbf32397..7c72855a9 100644 --- a/specutils/io/default_loaders/wcs_fits.py +++ b/specutils/io/default_loaders/wcs_fits.py @@ -5,6 +5,7 @@ from astropy.io import fits from astropy.wcs import WCS from astropy.modeling import models +from astropy.nddata import StdDevUncertainty, InverseVariance, VarianceUncertainty from astropy.utils.exceptions import AstropyUserWarning import numpy as np @@ -16,6 +17,14 @@ __all__ = ['wcs1d_fits_loader', 'non_linear_wcs1d_fits', 'non_linear_multispec_fits'] +UNCERT_REF = {'STD': StdDevUncertainty, + 'ERR': StdDevUncertainty, + 'UNCERT': StdDevUncertainty, + 'VAR': VarianceUncertainty, + 'IVAR': InverseVariance} + +UNCERT_EXP = {'std': 1, 'var': 2, 'ivar': -2} + def identify_wcs1d_fits(origin, *args, **kwargs): """ @@ -27,17 +36,16 @@ def identify_wcs1d_fits(origin, *args, **kwargs): # Default FITS format is BINTABLE in 1st extension HDU, unless IMAGE is # indicated via naming pattern or (explicitly) selecting primary HDU. if origin == 'write': - return ((args[0].endswith(('wcs.fits', 'wcs1d.fits', 'wcs.fit')) or - (args[0].endswith(('.fits', '.fit')) and whdu == 0)) and not - hasattr(args[2], 'uncertainty')) + return (args[0].endswith(('wcs.fits', 'wcs1d.fits', 'wcs.fit')) or + (args[0].endswith(('.fits', '.fit')) and whdu == 0)) hdu = kwargs.get('hdu', 0) # Check if number of axes is one and dimension of WCS is one with read_fileobj_or_hdulist(*args, **kwargs) as hdulist: return (hdulist[hdu].header.get('WCSDIM', 1) == 1 and (hdulist[hdu].header['NAXIS'] == 1 or - hdulist[hdu].header.get('WCSAXES', 0) == 1 or - hdulist[hdu].header.get('DISPAXIS', -1) >= 0 ) and not + hdulist[hdu].header.get('WCSAXES', 0) == 1 or + hdulist[hdu].header.get('DISPAXIS', -1) >= 0 ) and not hdulist[hdu].header.get('MSTITLE', 'n').startswith('2dF-SDSS LRG') and not # Check in CTYPE1 key for linear solution (single spectral axis) hdulist[hdu].header.get('CTYPE1', 'w').upper().startswith('MULTISPE')) @@ -46,7 +54,8 @@ def identify_wcs1d_fits(origin, *args, **kwargs): @data_loader("wcs1d-fits", identifier=identify_wcs1d_fits, dtype=Spectrum1D, extensions=['fits', 'fit'], priority=5) def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None, - hdu=None, verbose=False, **kwargs): + hdu=None, verbose=False, mask_hdu=True, + uncertainty_hdu=True, uncertainty_type=None, **kwargs): """ Loader for single spectrum-per-HDU spectra in FITS files, with the spectral axis stored in the header as FITS-WCS. The flux unit of the spectrum is @@ -68,10 +77,20 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None, Units of the flux for this spectrum. If not given (or None), the unit will be inferred from the BUNIT keyword in the header. Note that this unit will attempt to convert from BUNIT if BUNIT is present. - hdu : int, str or None. optional - The index or name of the HDU to load into this spectrum (default: search 1st). + hdu : int, str or None, optional + The index or name of the HDU to load into this spectrum + (default: find 1st applicable `ImageHDU`). verbose : bool. optional Print extra info. + mask_hdu : int, str, bool or None, optional + The index or name of the HDU to read mask from + (default: try to autodetect; `False`|`None`: do not read in). + uncertainty_hdu : int, str, bool or None, optional + The index or name of the HDU to read uncertainy from + (default: try to autodetect; `False`|`None`: do not read in). + uncertainty_type : str or None, optional + The ``uncertainty_type`` of `~astropy.nddata.NDUncertainty` + (one of 'std', 'var', 'ivar'; default: try to infer from HDU EXTNAME). **kwargs Extra keywords for :func:`~specutils.io.parsing_utils.read_fileobj_or_hdulist`. @@ -88,15 +107,22 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None, with read_fileobj_or_hdulist(file_obj, **kwargs) as hdulist: if hdu is None: - for ext in ('FLUX', 'SCI', 'PRIMARY'): + for ext in ('FLUX', 'SCI', 'DATA', 'PRIMARY'): # For now rely on extension containing spectral data. if ext in hdulist and ( isinstance(hdulist[ext], (fits.ImageHDU, fits.PrimaryHDU)) and hdulist[ext].data is not None): - hdu = ext - break - else: - raise ValueError('No HDU with spectral data found.') + if hdu is None or hdulist[ext] == hdulist[hdu]: + hdu = ext + continue + else: + warnings.warn(f"Found multiple data HDUs '{hdu}' and '{ext}', " + f"will read '{hdu}'! Please use `hdu=` " + "to select a specific one.", AstropyUserWarning) + break + + if hdu is None: + raise ValueError('No HDU with spectral data found.') header = hdulist[hdu].header wcs = WCS(header) @@ -108,11 +134,54 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None, else: data = u.Quantity(hdulist[hdu].data, unit=flux_unit) - # TODO: add additional HDU data like 'IVAR' - if 'MASK' in hdulist: - mask = hdulist['MASK'].data - else: - mask = None + # Read mask from specified HDU or try to auto-detect it. + mask = None + if mask_hdu is True: + for ext in ('MASK', 'DQ', 'QUALITY'): + if ext in hdulist and (isinstance(hdulist[ext], fits.ImageHDU) and + hdulist[ext].data is not None): + mask = hdulist[ext].data + break + elif isinstance(mask_hdu, (int, str)): + if mask_hdu in hdulist: + mask = hdulist[mask_hdu].data + else: + warnings.warn(f"No HDU '{mask_hdu}' for mask found in file.", + AstropyUserWarning) + + uncertainty = None + if uncertainty_hdu is True: + for ext in UNCERT_REF: + if ext in hdulist and (isinstance(hdulist[ext], fits.ImageHDU) and + hdulist[ext].data is not None): + uncertainty = hdulist[ext].data + break + elif isinstance(uncertainty_hdu, (int, str)): + if uncertainty_hdu in hdulist: + ext = uncertainty_hdu + uncertainty = hdulist[ext].data + else: + warnings.warn(f"No HDU '{uncertainty_hdu}' for uncertainty found in file.", + AstropyUserWarning) + + if uncertainty is not None: + if uncertainty_type is None: + if hdulist[ext].name in UNCERT_REF: + unc_type = hdulist[ext].name + else: + warnings.warn(f"Could not determine uncertainty type for HDU '{ext}' " + f"('{hdulist[ext].name}'), assuming 'StdDev'.", + AstropyUserWarning) + unc_type = 'STD' + elif str.upper(uncertainty_type) in UNCERT_REF: + unc_type = str.upper(uncertainty_type) + else: + raise ValueError(f"Invalid uncertainty type: '{uncertainty_type}'; " + "should be one of 'std', 'var', 'ivar'.") + uunit = u.Unit(header.get('BUNIT', flux_unit)) + if unc_type != 'STD': + uunit = uunit**UNCERT_EXP[unc_type.lower()] + uncertainty = UNCERT_REF[unc_type](u.Quantity(uncertainty, unit=uunit)) if spectral_axis_unit is not None: wcs.wcs.cunit[0] = str(spectral_axis_unit) @@ -134,14 +203,16 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None, meta = {'header': header} + # Is this restriction still appropriate? if wcs.naxis > 4: raise ValueError('FITS file input to wcs1d_fits_loader is > 4D') - return Spectrum1D(flux=data, wcs=wcs, mask=mask, meta=meta) + return Spectrum1D(flux=data, wcs=wcs, mask=mask, uncertainty=uncertainty, meta=meta) @custom_writer("wcs1d-fits") -def wcs1d_fits_writer(spectrum, file_name, hdu=0, update_header=False, flux_name='FLUX', **kwargs): +def wcs1d_fits_writer(spectrum, file_name, hdu=0, update_header=False, + flux_name='FLUX', mask_name='', uncertainty_name='', **kwargs): """ Write spectrum with spectral axis defined by its WCS to (primary) IMAGE_HDU of a FITS file. @@ -157,10 +228,14 @@ def wcs1d_fits_writer(spectrum, file_name, hdu=0, update_header=False, flux_name Update FITS header with all compatible entries in `spectrum.meta` flux_name : str, optional HDU name to store flux spectrum under (default 'FLUX') + mask_name : str or `None`, optional + HDU name to store mask under (default 'MASK'; `None`: do not save) + uncertainty_name : str or `None`, optional + HDU name to store uncertainty under (default set from type; `None`: do not save) unit : str or :class:`~astropy.units.Unit`, optional - Unit for the flux (and associated uncertainty; defaults to `spectrum.flux.unit`) + Unit for the flux (and associated uncertainty; defaults to ``spectrum.flux.unit``) dtype : str or :class:`~numpy.dtype`, optional - Floating point type for storing flux array (defaults to `spectrum.flux.dtype`) + Floating point type for storing flux array (defaults to ``spectrum.flux.dtype``) """ # Create HDU list from WCS try: @@ -171,16 +246,16 @@ def wcs1d_fits_writer(spectrum, file_name, hdu=0, update_header=False, flux_name raise ValueError(f'Only Spectrum1D objects with valid WCS can be written as wcs1d: {err}') # Verify spectral axis constructed from WCS - wl = spectrum.spectral_axis + disp = spectrum.spectral_axis # Not sure why the extra check is necessary for FITS WCS if hasattr(wcs, 'celestial') and wcs.celestial.naxis > 0: - dwl = (wcs.spectral.all_pix2world(np.arange(len(wl)), 0) - wl.value) / wl.value + ddisp = (wcs.spectral.all_pix2world(np.arange(len(disp)), 0) - disp.value) / disp.value else: - dwl = (wcs.all_pix2world(np.arange(len(wl)), 0) - wl.value) / wl.value - if np.abs(dwl).max() > 1.e-10: - m = np.abs(dwl).argmax() + ddisp = (wcs.all_pix2world(np.arange(len(disp)), 0) - disp.value) / disp.value + if np.abs(ddisp).max() > 1.e-10: + m = np.abs(ddisp).argmax() raise ValueError('Relative difference between WCS spectral axis and' - f'spectral_axis at {m:}: {dwl[m]}') + f'spectral_axis at {m:}: {ddisp[m]}') if update_header: hdr_types = (str, int, float, complex, bool, @@ -189,26 +264,51 @@ def wcs1d_fits_writer(spectrum, file_name, hdu=0, update_header=False, flux_name (isinstance(keyword[1], hdr_types) and keyword[0] not in ('NAXIS', 'NAXIS1', 'NAXIS2'))]) - # Cannot include uncertainty in IMAGE_HDU - maybe provide option to - # separately write this to BINARY_TBL extension later. - if spectrum.uncertainty is not None: - warnings.warn("Saving uncertainties in wcs1d format is not yet supported!", - AstropyUserWarning) - # Add flux array and unit ftype = kwargs.pop('dtype', spectrum.flux.dtype) funit = u.Unit(kwargs.pop('unit', spectrum.flux.unit)) - flux = spectrum.flux.to(funit, equivalencies=u.spectral_density(wl)) + flux = spectrum.flux.to(funit, equivalencies=u.spectral_density(disp)) hdulist[0].data = flux.value.astype(ftype) if flux_name is not None: hdulist[0].name = flux_name - # Append mask array (duplicate WCS for that extension)? - if spectrum.mask is not None: - hdulist.append(wcs.to_fits()[0]) - hdulist[1].data = spectrum.mask - hdulist[1].name = 'MASK' + # Append mask array + if spectrum.mask is not None and mask_name is not None: + hdulist.append(fits.ImageHDU()) + hdulist[-1].data = spectrum.mask + if mask_name == '': + hdulist[-1].name = 'MASK' + else: + hdulist[-1].name = mask_name hdu += 1 + # Warn if saving was requested (per explicitly choosing extension name). + elif mask_name is not None and mask_name != '': + warnings.warn("No mask found in this Spectrum1D, none saved.", AstropyUserWarning) + + # Append uncertainty array + if spectrum.uncertainty is not None and uncertainty_name is not None: + hdulist.append(fits.ImageHDU()) + uncertainty_name = str.upper(uncertainty_name) + uncertainty_type = type(spectrum.uncertainty) + # uncertainty - units to be inferred from and consistent with spectrum.flux. + if uncertainty_name in UNCERT_REF and uncertainty_type != UNCERT_REF[uncertainty_name]: + raise ValueError(f"Illegal label for uncertainty: '{uncertainty_name}' is reserved " + f"for {UNCERT_REF[uncertainty_name]}, not {uncertainty_type}.") + + if spectrum.uncertainty.uncertainty_type in ('var', 'ivar'): + uunit = funit**UNCERT_EXP[spectrum.uncertainty.uncertainty_type] + else: + uunit = funit + if uncertainty_name == '': + uncertainty_name = spectrum.uncertainty.uncertainty_type.upper() + sig = spectrum.uncertainty.quantity.to_value(uunit, equivalencies=u.spectral_density(disp)) + hdulist[-1].data = sig.astype(ftype) + hdulist[-1].name = uncertainty_name + hdu += 1 + # Warn if saving was requested (per explicitly choosing extension name). + elif uncertainty_name is not None and uncertainty_name != '': + warnings.warn("No uncertainty array found in this Spectrum1D, none saved.", + AstropyUserWarning) if hasattr(funit, 'long_names') and len(funit.long_names) > 0: comment = f'[{funit.long_names[0]}] {funit.physical_type}' @@ -383,7 +483,7 @@ def _read_non_linear_iraf_fits(file_obj, spectral_axis_unit=None, flux_unit=None Returns ------- - Tuple of data to pass to SpectrumCollection() or Spectrum1D(): + Tuple of data to pass to `SpectrumCollection`() or `Spectrum1D`(): spectral_axis : :class:`~astropy.units.Quantity` The spectral axis or axes as constructed from WCS(hdulist[0].header). diff --git a/specutils/tests/test_loaders.py b/specutils/tests/test_loaders.py index d667371bb..761c4decd 100644 --- a/specutils/tests/test_loaders.py +++ b/specutils/tests/test_loaders.py @@ -10,11 +10,13 @@ from astropy.io.fits.verify import VerifyWarning from astropy.table import Table from astropy.units import UnitsWarning -from astropy.wcs import FITSFixedWarning, InconsistentAxisTypesError, WCS +from astropy.wcs import FITSFixedWarning, WCS from astropy.io.registry import IORegistryError from astropy.modeling import models +from astropy.nddata import StdDevUncertainty, InverseVariance, VarianceUncertainty from astropy.tests.helper import quantity_allclose -from astropy.nddata import StdDevUncertainty +from astropy.utils.exceptions import AstropyUserWarning + from numpy.testing import assert_allclose from .conftest import remote_access @@ -751,8 +753,9 @@ def test_wcs1d_fits_writer(tmp_path, spectral_axis, with_mask, uncertainty): assert quantity_allclose(spec.uncertainty.quantity, spectrum.uncertainty.quantity) # Read from HDUList - hdulist = fits.open(tmpfile) - spec = Spectrum1D.read(hdulist, format='wcs1d-fits') + with fits.open(tmpfile) as hdulist: + spec = Spectrum1D.read(hdulist, format='wcs1d-fits') + assert isinstance(spec, Spectrum1D) assert quantity_allclose(spec.spectral_axis, spectrum.spectral_axis) assert quantity_allclose(spec.flux, spectrum.flux) @@ -845,22 +848,17 @@ def test_wcs1d_fits_uncertainty(tmp_path, uncertainty_rsv): ensure it raises on illegal (reserved) names. """ # Header dictionary for constructing WCS - hdr = {'CTYPE1': 'WAVE', 'CUNIT1': 'Angstrom', - 'CRPIX1': 1, 'CRVAL1': 1, 'CDELT1': 0.01} + hdr = {'CTYPE1': 'WAVE', 'CUNIT1': 'um', 'CRPIX1': 1, 'CRVAL1': 1, 'CDELT1': 0.001} # Reserved EXTNAMEs for uncertainty types UNCERT_REF = {'STD': StdDevUncertainty, 'ERR': StdDevUncertainty, 'UNCERT': StdDevUncertainty, 'VAR': VarianceUncertainty, 'IVAR': InverseVariance} # Alternative EXTNAMEs for uncertainty types UNCERT_ALT = {'std': 'StdErr', 'var': 'VARIAN', 'ivar': 'InvVar'} - # Create a small data set - flux = np.arange(1, 11)**2 * 1.e-14 * u.Jy - wlu = u.Unit(hdr['CUNIT1']) - wl0 = hdr['CRVAL1'] - dwl = hdr['CDELT1'] - disp = np.arange(wl0, wl0 + (len(flux) - 0.5) * dwl, dwl) * wlu tmpfile = tmp_path / 'wcs_tst.fits' + # Create a small data set + flux = np.arange(1, 11)**2 * 1.e-14 * u.Jy mask = np.array([0, 0, 0, 0, 1, 0, 0, 0, 1, 0], dtype=np.uint16) # Set uncertainty to mismatched type @@ -872,6 +870,7 @@ def test_wcs1d_fits_uncertainty(tmp_path, uncertainty_rsv): f"is reserved for {UNCERT_REF[uncertainty_rsv]}, not {uncertainty}."): spectrum.write(tmpfile, format='wcs1d-fits', uncertainty_name=uncertainty_rsv) + # Set permitted custom name uncertainty_type = spectrum.uncertainty.uncertainty_type uncertainty_alt = UNCERT_ALT[uncertainty_type] spectrum.write(tmpfile, format='wcs1d-fits', uncertainty_name=uncertainty_alt) From 1d5a59fda82001ced92438857baa2976f15d4e01 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen <39831871+rosteen@users.noreply.github.com> Date: Fri, 17 Mar 2023 14:11:01 -0400 Subject: [PATCH 5/5] Apply suggestions from code review Co-authored-by: P. L. Lim <2090236+pllim@users.noreply.github.com> --- specutils/io/default_loaders/wcs_fits.py | 6 +++--- specutils/tests/test_loaders.py | 5 ++--- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/specutils/io/default_loaders/wcs_fits.py b/specutils/io/default_loaders/wcs_fits.py index 7c72855a9..446c2e63b 100644 --- a/specutils/io/default_loaders/wcs_fits.py +++ b/specutils/io/default_loaders/wcs_fits.py @@ -117,8 +117,8 @@ def wcs1d_fits_loader(file_obj, spectral_axis_unit=None, flux_unit=None, continue else: warnings.warn(f"Found multiple data HDUs '{hdu}' and '{ext}', " - f"will read '{hdu}'! Please use `hdu=` " - "to select a specific one.", AstropyUserWarning) + f"will read '{hdu}'. Please use `hdu=` " + "to select a specific extension.", AstropyUserWarning) break if hdu is None: @@ -483,7 +483,7 @@ def _read_non_linear_iraf_fits(file_obj, spectral_axis_unit=None, flux_unit=None Returns ------- - Tuple of data to pass to `SpectrumCollection`() or `Spectrum1D`(): + Tuple of data to pass to `~specutils.SpectrumCollection` or `~specutils.Spectrum1D`: spectral_axis : :class:`~astropy.units.Quantity` The spectral axis or axes as constructed from WCS(hdulist[0].header). diff --git a/specutils/tests/test_loaders.py b/specutils/tests/test_loaders.py index 761c4decd..02dd46cd1 100644 --- a/specutils/tests/test_loaders.py +++ b/specutils/tests/test_loaders.py @@ -978,9 +978,8 @@ def test_wcs1d_fits_non1d(tmp_path, spectral_axis): spectrum.write(tmpfile, format='wcs1d-fits') # Reopen file and update header with off-diagonal element - hdulist = fits.open(tmpfile, mode='update') - hdulist[0].header.update([('PC1_2', 0.2)]) - hdulist.close() + with fits.open(tmpfile, mode='update') as hdulist: + hdulist[0].header.update([('PC1_2', 0.2)]) with pytest.raises(ValueError, match='Non-zero off-diagonal matrix elements excluded from the subimage.'):