diff --git a/CHANGES.rst b/CHANGES.rst index 0e44bf9..32afec1 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,6 +1,16 @@ -1.3.1 (unreleased) +1.4.0 (unreleased) ================== +- ``read_fits_spec()`` now uses ``astropy.table.QTable.read`` for parsing to + ensure that the correct ``TUNITn`` is read. As a result, ``wave_unit`` and + ``flux_unit`` keywords are deprecated and no longer used in that function. + Additionally, if any ``TUNITn`` in the table is invalid, regardless whether + the column is used or not, an exception will now be raised. The inputs + for ``wave_col`` and ``flux_col`` are now case-sensitive. [#384] + +- ``read_spec()`` now detects whether given filename is FITS more consistently + w.r.t. ``astropy``. [#384] + 1.3.0 (2023-11-28) ================== diff --git a/synphot/reddening.py b/synphot/reddening.py index 6ee9b9f..584cd51 100644 --- a/synphot/reddening.py +++ b/synphot/reddening.py @@ -7,6 +7,7 @@ # THIRD-PARTY import numpy as np from astropy import units as u +from astropy.io.fits.connect import is_fits # LOCAL from synphot import exceptions, specio, units @@ -137,8 +138,8 @@ def to_fits(self, filename, wavelengths=None, **kwargs): def from_file(cls, filename, **kwargs): """Create a reddening law from file. - If filename has 'fits' or 'fit' suffix, it is read as FITS. - Otherwise, it is read as ASCII. + If filename is recognized by ``astropy.io.fits`` as FITS, + it is read as such. Otherwise, it is read as ASCII. Parameters ---------- @@ -156,13 +157,12 @@ def from_file(cls, filename, **kwargs): Empirical reddening law. """ - if 'flux_unit' not in kwargs: + if is_fits("", filename, None): + if 'flux_col' not in kwargs: + kwargs['flux_col'] = 'Av/E(B-V)' + elif 'flux_unit' not in kwargs: kwargs['flux_unit'] = cls._internal_flux_unit - if ((filename.endswith('fits') or filename.endswith('fit')) and - 'flux_col' not in kwargs): - kwargs['flux_col'] = 'Av/E(B-V)' - header, wavelengths, rvs = specio.read_spec(filename, **kwargs) return cls(Empirical1D, points=wavelengths, lookup_table=rvs, @@ -217,13 +217,12 @@ def from_extinction_model(cls, modelname, **kwargs): filename = cfgitem() - if 'flux_unit' not in kwargs: + if is_fits("", filename, None): + if 'flux_col' not in kwargs: + kwargs['flux_col'] = 'Av/E(B-V)' + elif 'flux_unit' not in kwargs: kwargs['flux_unit'] = cls._internal_flux_unit - if ((filename.endswith('fits') or filename.endswith('fit')) and - 'flux_col' not in kwargs): - kwargs['flux_col'] = 'Av/E(B-V)' - header, wavelengths, rvs = specio.read_remote_spec(filename, **kwargs) header['filename'] = filename header['descrip'] = cfgitem.description diff --git a/synphot/specio.py b/synphot/specio.py index 91e2a45..99ac459 100644 --- a/synphot/specio.py +++ b/synphot/specio.py @@ -12,7 +12,10 @@ from astropy import log from astropy import units as u from astropy.io import ascii, fits +from astropy.io.fits.connect import is_fits +from astropy.table import QTable from astropy.utils.data import get_readable_fileobj +from astropy.utils.decorators import deprecated_renamed_argument from astropy.utils.exceptions import AstropyUserWarning # LOCAL @@ -88,7 +91,7 @@ def read_spec(filename, fname='', **kwargs): elif not fname: # pragma: no cover raise exceptions.SynphotError('Cannot determine filename.') - if fname.endswith('fits') or fname.endswith('fit'): + if is_fits("", fname, None): read_func = read_fits_spec else: read_func = read_ascii_spec @@ -143,12 +146,15 @@ def read_ascii_spec(filename, wave_unit=u.AA, flux_unit=units.FLAM, **kwargs): return header, wavelengths, fluxes +@deprecated_renamed_argument( + ["wave_unit", "flux_unit"], [None, None], ["1.4", "1.4"], + alternative='TUNITn as per FITS standards') def read_fits_spec(filename, ext=1, wave_col='WAVELENGTH', flux_col='FLUX', wave_unit=u.AA, flux_unit=units.FLAM): """Read FITS spectrum. - Wavelength and flux units are extracted from ``TUNIT1`` and ``TUNIT2`` - keywords, respectively, from data table (not primary) header. + Wavelength and flux units are extracted from respective ``TUNITn`` + keywords, from data table (not primary) header. If these keywords are not present, units are taken from ``wave_unit`` and ``flux_unit`` instead. @@ -161,12 +167,14 @@ def read_fits_spec(filename, ext=1, wave_col='WAVELENGTH', flux_col='FLUX', FITS extension with table data. Default is 1. wave_col, flux_col : str - Wavelength and flux column names (case-insensitive). + Wavelength and flux column names (case-sensitive). wave_unit, flux_unit : str or `~astropy.units.Unit` - Wavelength and flux units, which default to Angstrom and FLAM, - respectively. These are *only* used if ``TUNIT1`` and ``TUNIT2`` - keywords are not present in table (not primary) header. + Wavelength and flux units. These are *no longer used*. + Define your units in the respective ``TUNITn`` + keywords in table (not primary) header. + + .. deprecated:: 1.4 Returns ------- @@ -179,35 +187,23 @@ def read_fits_spec(filename, ext=1, wave_col='WAVELENGTH', flux_col='FLUX', """ try: fs = fits.open(filename) - header = dict(fs[str('PRIMARY')].header) - wave_dat = fs[ext].data.field(wave_col).copy() - flux_dat = fs[ext].data.field(flux_col).copy() - fits_wave_unit = fs[ext].header.get('TUNIT1') - fits_flux_unit = fs[ext].header.get('TUNIT2') - - if fits_wave_unit is not None: - try: - wave_unit = units.validate_unit(fits_wave_unit) - except (exceptions.SynphotError, ValueError) as e: # pragma: no cover # noqa: E501 - warnings.warn( - '{0} from FITS header is not valid wavelength unit, using ' - '{1}: {2}'.format(fits_wave_unit, wave_unit, e), - AstropyUserWarning) - - if fits_flux_unit is not None: - try: - flux_unit = units.validate_unit(fits_flux_unit) - except (exceptions.SynphotError, ValueError) as e: # pragma: no cover # noqa: E501 - warnings.warn( - '{0} from FITS header is not valid flux unit, using ' - '{1}: {2}'.format(fits_flux_unit, flux_unit, e), - AstropyUserWarning) - - wave_unit = units.validate_unit(wave_unit) - flux_unit = units.validate_unit(flux_unit) - - wavelengths = wave_dat * wave_unit - fluxes = flux_dat * flux_unit + subhdu = fs[ext] + + # Need to fix table units + for key in subhdu.header["TUNIT*"]: + val = subhdu.header[key] + if not val: + continue + newval = units.validate_unit(val) + subhdu.header[key] = newval.to_string("fits") + + t = QTable.read(subhdu) + header = dict(fs["PRIMARY"].header) + t_col_wave = t[wave_col] + wavelengths = t_col_wave.value * (t_col_wave.unit or u.dimensionless_unscaled) # noqa: E501 + t_col_flux = t[flux_col] + fluxes = t_col_flux.value * (t_col_flux.unit or u.dimensionless_unscaled) # noqa: E501 + finally: if isinstance(filename, str): fs.close() diff --git a/synphot/spectrum.py b/synphot/spectrum.py index f1a6173..f730e2d 100644 --- a/synphot/spectrum.py +++ b/synphot/spectrum.py @@ -14,6 +14,7 @@ # ASTROPY from astropy import log from astropy import units as u +from astropy.io.fits.connect import is_fits from astropy.modeling import Model from astropy.modeling.core import CompoundModel from astropy.modeling.models import RedshiftScaleFactor, Scale @@ -1921,8 +1922,8 @@ def to_fits(self, filename, wavelengths=None, **kwargs): def from_file(cls, filename, **kwargs): """Creates a bandpass from file. - If filename has 'fits' or 'fit' suffix, it is read as FITS. - Otherwise, it is read as ASCII. + If filename is recognized by ``astropy.io.fits`` as FITS, + it is read as such. Otherwise, it is read as ASCII. Parameters ---------- @@ -1940,13 +1941,12 @@ def from_file(cls, filename, **kwargs): Empirical bandpass. """ - if 'flux_unit' not in kwargs: + if is_fits("", filename, None): + if 'flux_col' not in kwargs: + kwargs['flux_col'] = 'THROUGHPUT' + elif 'flux_unit' not in kwargs: kwargs['flux_unit'] = cls._internal_flux_unit - if ((filename.endswith('fits') or filename.endswith('fit')) and - 'flux_col' not in kwargs): - kwargs['flux_col'] = 'THROUGHPUT' - header, wavelengths, throughput = specio.read_spec(filename, **kwargs) return cls(Empirical1D, points=wavelengths, lookup_table=throughput, keep_neg=True, meta={'header': header}) @@ -2009,13 +2009,12 @@ def from_filter(cls, filtername, **kwargs): filename = cfgitem() - if 'flux_unit' not in kwargs: + if is_fits("", filename, None): + if 'flux_col' not in kwargs: + kwargs['flux_col'] = 'THROUGHPUT' + elif 'flux_unit' not in kwargs: kwargs['flux_unit'] = cls._internal_flux_unit - if ((filename.endswith('fits') or filename.endswith('fit')) and - 'flux_col' not in kwargs): - kwargs['flux_col'] = 'THROUGHPUT' - header, wavelengths, throughput = specio.read_remote_spec( filename, **kwargs) header['filename'] = filename diff --git a/synphot/tests/test_binning.py b/synphot/tests/test_binning.py index 79dd482..1815342 100644 --- a/synphot/tests/test_binning.py +++ b/synphot/tests/test_binning.py @@ -198,7 +198,7 @@ def setup_class(self): get_pkg_data_filename( os.path.join('data', 'hst_acs_hrc_f555w.fits'), package='synphot.tests'), - flux_col='THROUGHPUT', flux_unit=u.dimensionless_unscaled) + flux_col='THROUGHPUT') # Binned data. bins = generate_wavelengths( diff --git a/synphot/tests/test_specio.py b/synphot/tests/test_specio.py index 7a4e56c..ff3bc4a 100644 --- a/synphot/tests/test_specio.py +++ b/synphot/tests/test_specio.py @@ -3,8 +3,6 @@ # STDLIB import os -import shutil -import tempfile # THIRD-PARTY import numpy as np @@ -15,10 +13,12 @@ from astropy.io import fits from astropy.tests.helper import assert_quantity_allclose from astropy.utils.data import get_pkg_data_filename -from astropy.utils.exceptions import AstropyUserWarning +from astropy.utils.exceptions import ( + AstropyUserWarning, AstropyDeprecationWarning) # LOCAL from synphot import exceptions, specio, units +from synphot.spectrum import SpectralElement @pytest.mark.remote_data @@ -55,7 +55,6 @@ class TestReadWriteFITS: """Test read/write FITS spectrum.""" def setup_class(self): self.epsilon = 0.00031 - self.outdir = tempfile.mkdtemp() self.wave = np.array([1000.0, 2000.0, 2000.0 + self.epsilon, 3000.0, 4000.0, 5000.0], dtype=np.float64) * u.AA self.flux = np.array([0.1, 100.2, 10.0, 0.0, 6.5, 1.2], @@ -63,9 +62,9 @@ def setup_class(self): self.prihdr = {'PEDIGREE': 'DUMMY'} self.scihdr = {'SPEC_SRC': 'RANDOM'} - def test_array_data(self): + def test_array_data(self, tmp_path): """Data as Numpy array.""" - outfile = os.path.join(self.outdir, 'outspec1.fits') + outfile = str(tmp_path / 'outspec1.fits') # Write it out with pytest.warns(AstropyUserWarning, match=r'rows are thrown out'): @@ -76,7 +75,8 @@ def test_array_data(self): wave_unit=self.wave.unit, flux_unit=self.flux.unit) # Read it back in and check values (flux_unit should be ignored) - hdr, wave, flux = specio.read_spec(outfile, flux_unit='foo') + with pytest.warns(AstropyDeprecationWarning, match=r"\"flux_unit\" was deprecated"): # noqa: E501 + hdr, wave, flux = specio.read_spec(outfile, flux_unit='foo') # Compare data np.testing.assert_allclose( @@ -95,9 +95,9 @@ def test_array_data(self): assert sci_hdr['SPEC_SRC'] == 'RANDOM' assert sci_hdr['TFORM2'].lower() == 'e' - def test_quantity_data(self): + def test_quantity_data(self, tmp_path): """Data as Quantity.""" - outfile = os.path.join(self.outdir, 'outspec2.fits') + outfile = str(tmp_path / 'outspec2.fits') # Write it out (flux_unit should be ignored) specio.write_fits_spec( @@ -105,7 +105,8 @@ def test_quantity_data(self): ext_header=self.scihdr, precision='double', flux_unit='foo') # Read it back in and check values (flux_unit should be ignored) - hdr, wave, flux = specio.read_spec(outfile, flux_unit='foo') + with pytest.warns(AstropyDeprecationWarning, match=r"\"flux_unit\" was deprecated"): # noqa: E501 + hdr, wave, flux = specio.read_spec(outfile, flux_unit='foo') # Compare data (trim_zero=True, pad_zero_ends=True) np.testing.assert_allclose( @@ -125,9 +126,9 @@ def test_quantity_data(self): assert sci_hdr['SPEC_SRC'] == 'RANDOM' assert sci_hdr['TFORM2'].lower() == 'd' - def test_exceptions(self): + def test_exceptions(self, tmp_path): """Test for appropriate exceptions.""" - outfile = os.path.join(self.outdir, 'outspec3.fits') + outfile = str(tmp_path / 'outspec3.fits') # Shape mismatch with pytest.raises(exceptions.SynphotError): @@ -149,5 +150,55 @@ def test_exceptions(self): specio.write_fits_spec( outfile, self.wave, np.arange(6), overwrite=True) - def teardown_class(self): - shutil.rmtree(self.outdir) + +def test_read_nonstandard_fits_cols_01(tmp_path): + """See https://github.com/spacetelescope/synphot_refactor/issues/372""" + pix = np.arange(5, dtype=int) + 1 + wav = (pix * 0.1) * u.micron + trace = np.array([0, 0.5, 1, 0.9, 0]) + coldefs = fits.ColDefs([ + fits.Column(name="X", format="I", array=pix), + fits.Column(name="WAVELENGTH", format="E", + unit=wav.unit.to_string(format="fits"), array=wav.value), + fits.Column(name="TRACE", format="E", array=trace)]) + hdulist = fits.HDUList([ + fits.PrimaryHDU(), + fits.BinTableHDU.from_columns(coldefs)]) + outfile = str(tmp_path / "jwst_niriss_soss_trace.fits") + hdulist.writeto(outfile, overwrite=True) + + tr = SpectralElement.from_file(outfile, flux_col="TRACE") + assert_quantity_allclose(tr.waveset, wav) + assert_quantity_allclose(tr(wav), trace, atol=1e-7) + + +def test_read_nonstandard_fits_cols_02(tmp_path): + """See https://github.com/spacetelescope/synphot_refactor/issues/372""" + + wav = (np.arange(5) + 1) * u.nm + flux_unit_str = "ph/s/m2/micron/arcsec2" # Invalid but should not matter. + flux = np.ones(5) + thru = np.array([0, 0.5, 1, 0.9, 0]) + coldefs = fits.ColDefs([ + fits.Column(name="lam", format="E", + unit=wav.unit.to_string(format="fits"), array=wav.value), + fits.Column(name="flux", format="E", + unit=flux_unit_str, array=flux), + fits.Column(name="dflux1", format="E", + unit=flux_unit_str, array=flux), + fits.Column(name="dflux2", format="E", + unit=flux_unit_str, array=flux), + fits.Column(name="trans", format="E", unit="1", array=thru)]) + hdulist = fits.HDUList([ + fits.PrimaryHDU(), + fits.BinTableHDU.from_columns(coldefs)]) + outfile = str(tmp_path / "skytable.fits") + hdulist.writeto(outfile, overwrite=True) + + with pytest.warns(u.UnitsWarning, match="'ph/s/m2/micron/arcsec2'"): # noqa: E501 + header, wavelengths, transmission = specio.read_spec( + outfile, wave_col="lam", flux_col="trans") + + assert header["SIMPLE"] + assert_quantity_allclose(wavelengths, wav) + assert_quantity_allclose(transmission, thru) diff --git a/synphot/thermal.py b/synphot/thermal.py index 093d412..8ad4ff7 100644 --- a/synphot/thermal.py +++ b/synphot/thermal.py @@ -4,6 +4,7 @@ # ASTROPY from astropy import units as u from astropy.io import fits +from astropy.io.fits.connect import is_fits # LOCAL from synphot import exceptions, specio, units @@ -124,7 +125,7 @@ def from_file(cls, filename, temperature_key='DEFT', Invalid inputs. """ - if not (filename.endswith('fits') or filename.endswith('fit')): + if not is_fits("", filename, None): raise exceptions.SynphotError('Only FITS format is supported.') # Extra info from table header @@ -138,9 +139,6 @@ def from_file(cls, filename, temperature_key='DEFT', beam_fill_factor = tab_hdr.get('BEAMFILL', 1) - if 'flux_unit' not in kwargs: - kwargs['flux_unit'] = cls._internal_flux_unit - if 'flux_col' not in kwargs: kwargs['flux_col'] = 'EMISSIVITY' diff --git a/synphot/units.py b/synphot/units.py index e3c4f9f..c1b7ead 100644 --- a/synphot/units.py +++ b/synphot/units.py @@ -295,6 +295,18 @@ def validate_unit(input_unit): output_unit = THROUGHPUT elif input_unit_lowcase == 'jy': output_unit = u.Jy + elif input_unit_lowcase == "flam": + output_unit = FLAM + elif input_unit_lowcase == "fnu": + output_unit = FNU + elif input_unit_lowcase == "photlam": + output_unit = PHOTLAM + elif input_unit_lowcase == "photnu": + output_unit = PHOTNU + elif input_unit_lowcase == "none": + output_unit = u.dimensionless_unscaled + elif input_unit_lowcase == "sec": + output_unit = u.s # Work around mag unit limitations elif input_unit_lowcase in ('stmag', 'mag(st)'):