Merge pull request #384 from pllim/fits-io-table

BUG/API: read_fits_spec now more flexible but at a cost

Author: pllim

CHANGES.rst

@@ -1,6 +1,15 @@
-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. [#384]
+
+- ``read_spec()`` now detects whether given filename is FITS more consistently
+  w.r.t. ``astropy``. [#384]
+
 1.3.0 (2023-11-28)
 ==================
 

docs/synphot/overview.rst

@@ -396,9 +396,7 @@ keywords (unless you overwrite them with non-default values in
 * ``TTYPE2`` set to "FLUX" (for source spectrum) or "THROUGHPUT"
   (for bandpass).
 
-While these were required in ASTROLIB PYSYNPHOT, they are optional here in that
-default units would be applied, where applicable, if they are missing from
-the header. Regardless, setting them is highly recommended:
+These were required in ASTROLIB PYSYNPHOT, as well as this package:
 
 * ``TUNIT1`` set to :ref:`supported wavelength unit name <synphot_units>`.
 * ``TUNIT2`` set to :ref:`supported flux unit name <synphot_units>`

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:  # pragma: no cover
             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:  # pragma: no cover
             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

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.
 
@@ -164,9 +170,11 @@ def read_fits_spec(filename, ext=1, wave_col='WAVELENGTH', flux_col='FLUX',
         Wavelength and flux column names (case-insensitive).
 
     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
     -------
@@ -177,37 +185,42 @@ def read_fits_spec(filename, ext=1, wave_col='WAVELENGTH', flux_col='FLUX',
         Wavelength and flux of the spectrum.
 
     """
+    wave_col = wave_col.lower()
+    flux_col = flux_col.lower()
+
     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()  # Must be generic to handle mag  # noqa: E501
+
+        with warnings.catch_warnings():
+            warnings.filterwarnings("ignore", category=u.UnitsWarning,
+                                    message=".* did not parse as fits unit")
+            t = QTable.read(subhdu)
+        header = dict(fs["PRIMARY"].header)
+
+        # https://github.com/astropy/astropy/issues/16221
+        lower_colnames = [c.lower() for c in t.colnames]
+        t_col_wave = t.columns[lower_colnames.index(wave_col)]
+        if t_col_wave.unit:
+            t_col_wave_unit = units.validate_unit(t_col_wave.unit.to_string())
+        else:
+            t_col_wave_unit = u.dimensionless_unscaled
+        t_col_flux = t.columns[lower_colnames.index(flux_col)]
+        if t_col_flux.unit:
+            t_col_flux_unit = units.validate_unit(t_col_flux.unit.to_string())
+        else:
+            t_col_flux_unit = u.dimensionless_unscaled
+        wavelengths = t_col_wave.value * t_col_wave_unit
+        fluxes = t_col_flux.value * t_col_flux_unit
+
     finally:
         if isinstance(filename, str):
             fs.close()

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:  # pragma: no cover
             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:  # pragma: no cover
             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

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(

synphot/tests/test_reddening.py

@@ -3,8 +3,6 @@
 
 # STDLIB
 import os
-import shutil
-import tempfile
 
 # THIRD-PARTY
 import numpy as np
@@ -156,32 +154,26 @@ def test_redlaw_from_model_exception():
         ReddeningLaw.from_extinction_model('foo')
 
 
-class TestWriteReddeningLaw:
+@pytest.mark.parametrize('ext_hdr', [None, {'foo': 'foo'}])
+def test_write_reddening_law(tmp_path, ext_hdr):
     """Test ReddeningLaw ``to_fits()`` method."""
-    def setup_class(self):
-        self.outdir = tempfile.mkdtemp()
-        self.x = np.linspace(1000, 5000, 5)
-        self.y = np.linspace(1, 5, 5) * 0.1
-        self.redlaw = ReddeningLaw(
-            Empirical1D, points=self.x, lookup_table=self.y)
-
-    @pytest.mark.parametrize('ext_hdr', [None, {'foo': 'foo'}])
-    def test_write(self, ext_hdr):
-        outfile = os.path.join(self.outdir, 'outredlaw.fits')
-
-        if ext_hdr is None:
-            self.redlaw.to_fits(outfile, overwrite=True)
-        else:
-            self.redlaw.to_fits(outfile, overwrite=True, ext_header=ext_hdr)
-
-        # Read it back in and check
-        redlaw2 = ReddeningLaw.from_file(outfile)
-        np.testing.assert_allclose(redlaw2.waveset.value, self.x)
-        np.testing.assert_allclose(redlaw2(self.x).value, self.y)
-
-        if ext_hdr is not None:
-            hdr = fits.getheader(outfile, 1)
-            assert 'foo' in hdr
-
-    def teardown_class(self):
-        shutil.rmtree(self.outdir)
+    x = np.linspace(1000, 5000, 5)
+    y = np.linspace(1, 5, 5) * 0.1
+    redlaw = ReddeningLaw(
+        Empirical1D, points=x, lookup_table=y, meta={"expr": "ebv(test)"})
+
+    outfile = str(tmp_path / 'outredlaw.fits')
+
+    if ext_hdr is None:
+        redlaw.to_fits(outfile, overwrite=True)
+    else:
+        redlaw.to_fits(outfile, overwrite=True, ext_header=ext_hdr)
+
+    # Read it back in and check
+    redlaw2 = ReddeningLaw.from_file(outfile)
+    np.testing.assert_allclose(redlaw2.waveset.value, x)
+    np.testing.assert_allclose(redlaw2(x).value, y)
+
+    if ext_hdr is not None:
+        hdr = fits.getheader(outfile, 1)
+        assert 'foo' in hdr