diff --git a/specutils/io/default_loaders/tabular_fits.py b/specutils/io/default_loaders/tabular_fits.py index 7fd790ec2..74d0d4db8 100644 --- a/specutils/io/default_loaders/tabular_fits.py +++ b/specutils/io/default_loaders/tabular_fits.py @@ -1,7 +1,7 @@ import numpy as np from astropy.io import fits -from astropy.nddata import StdDevUncertainty +from astropy.nddata import StdDevUncertainty, Covariance from astropy.table import Table import astropy.units as u from astropy.wcs import WCS @@ -47,11 +47,11 @@ def tabular_fits_loader(file_obj, column_mapping=None, hdu=1, store_data_header= Parameters ---------- - file_obj : str, file-like, or :class:`~astropy.io.fits.HDUList` - FITS file name, object (provided from name by Astropy I/O Registry), - or HDU list (as resulting from `~astropy.io.fits.open`). + file_obj : str, file-like, or `~astropy.io.fits.HDUList` + FITS file name, object (provided from name by Astropy I/O Registry), or + `~astropy.io.fits.HDUList` (as resulting from `astropy.io.fits.open`). hdu : int - The HDU of the fits file (default: 1st extension) to read from + The HDU of the fits file (default: 1st extension) to read store_data_header : bool Defaults to ``False``, which stores the primary header in ``Spectrum.meta['header']``. Set to ``True`` to instead store the header from the specified data HDU. @@ -86,6 +86,9 @@ def tabular_fits_loader(file_obj, column_mapping=None, hdu=1, store_data_header= else: tab.meta = hdulist[0].header + # Determine if there is a covariance matrix + covar = Table.read(hdulist['COVAR']) if 'COVAR' in [h.name for h in hdulist] else None + # Minimal checks for wcs consistency with table data - # assume 1D spectral axis (having shape (0, NAXIS1), # or alternatively compare against shape of 1st column. @@ -96,9 +99,9 @@ def tabular_fits_loader(file_obj, column_mapping=None, hdu=1, store_data_header= # If no column mapping is given, attempt to parse the file using # unit information if column_mapping is None: - return generic_spectrum_from_table(tab, wcs=wcs) + return generic_spectrum_from_table(tab, wcs=wcs, covar=covar, **kwargs) - return spectrum_from_column_mapping(tab, column_mapping, wcs=wcs) + return spectrum_from_column_mapping(tab, column_mapping, wcs=wcs, covar=covar) @custom_writer("tabular-fits") @@ -131,6 +134,7 @@ def tabular_fits_writer(spectrum, file_name, hdu=1, update_header=False, store_d **kwargs Additional optional keywords passed to :func:`~astropy.io.fits.HDUList.writeto`. """ + # TODO: `hdu` is not used below. Is this necessary? if hdu < 1: raise ValueError(f'FITS does not support BINTABLE extension in HDU {hdu}.') @@ -144,7 +148,7 @@ def tabular_fits_writer(spectrum, file_name, hdu=1, update_header=False, store_d isinstance(keyword[1], hdr_types)]) # Strip header of FITS reserved keywords - for keyword in ['NAXIS', 'NAXIS1', 'NAXIS2']: + for keyword in ['EXTNAME', 'NAXIS', 'NAXIS1', 'NAXIS2']: header.remove(keyword, ignore_missing=True) # Add dispersion array and unit @@ -168,20 +172,26 @@ def tabular_fits_writer(spectrum, file_name, hdu=1, update_header=False, store_d colnames = [dispname, "flux"] # Include uncertainty - units to be inferred from spectrum.flux + covar = None if spectrum.uncertainty is not None: - try: - unc = ( - spectrum - .uncertainty - .represent_as(StdDevUncertainty) - .quantity - .to(funit, equivalencies=u.spectral_density(disp)) - ) - columns.append(unc.astype(ftype)) + if isinstance(spectrum.uncertainty, Covariance): + covar = spectrum.uncertainty.to_table() + columns.append(np.sqrt(spectrum.uncertainty.variance) * funit) colnames.append("uncertainty") - except RuntimeWarning: - raise ValueError("Could not convert uncertainty to StdDevUncertainty due" - " to divide-by-zero error.") + else: + try: + unc = ( + spectrum + .uncertainty + .represent_as(StdDevUncertainty) + .quantity + .to(funit, equivalencies=u.spectral_density(disp)) + ) + columns.append(unc.astype(ftype)) + colnames.append("uncertainty") + except RuntimeWarning: + raise ValueError("Could not convert uncertainty to StdDevUncertainty due" + " to divide-by-zero error.") # Add mask column if present if spectrum.mask is not None: @@ -197,6 +207,7 @@ def tabular_fits_writer(spectrum, file_name, hdu=1, update_header=False, store_d colnames.append('mask') # For > 1D data transpose from row-major format + # TODO: revisit this for c in range(1, len(columns)): if columns[c].ndim > 1: columns[c] = columns[c].T @@ -204,15 +215,14 @@ def tabular_fits_writer(spectrum, file_name, hdu=1, update_header=False, store_d tab = Table(columns, names=colnames) if store_data_header: hdu0 = fits.PrimaryHDU() - hdu1 = fits.BinTableHDU(data=tab, header=header) + hdu1 = fits.BinTableHDU(data=tab, header=header, name='DATA') else: hdu0 = fits.PrimaryHDU(header=header) - hdu1 = fits.BinTableHDU(data=tab) - - # This will overwrite any 'EXTNAME' previously read from a valid header; should it? - hdu1.header.update(EXTNAME='DATA') + hdu1 = fits.BinTableHDU(data=tab, name='DATA') hdulist = fits.HDUList([hdu0, hdu1]) + if covar is not None: + hdulist.append(fits.BinTableHDU(data=covar, name='COVAR')) # TODO: Use output_verify options to check for valid FITS hdulist.writeto(file_name, **kwargs) diff --git a/specutils/io/parsing_utils.py b/specutils/io/parsing_utils.py index 6d35028a7..e0420ee23 100644 --- a/specutils/io/parsing_utils.py +++ b/specutils/io/parsing_utils.py @@ -6,7 +6,7 @@ import contextlib from astropy.io import fits -from astropy.nddata import StdDevUncertainty +from astropy.nddata import StdDevUncertainty, Covariance from astropy.utils.exceptions import AstropyUserWarning import astropy.units as u import warnings @@ -52,7 +52,7 @@ def read_fileobj_or_hdulist(*args, **kwargs): hdulist.close() -def spectrum_from_column_mapping(table, column_mapping, wcs=None, verbose=False): +def spectrum_from_column_mapping(table, column_mapping, wcs=None, covar=None, verbose=False): """ Given a table and a mapping of the table column names to attributes on the Spectrum object, parse the information into a Spectrum. @@ -76,6 +76,10 @@ def spectrum_from_column_mapping(table, column_mapping, wcs=None, verbose=False) wcs : :class:`~astropy.wcs.WCS` or :class:`gwcs.WCS` WCS object passed to the Spectrum initializer. + covar : `astropy.table.Table`, optional + Table providing a covariance matrix for the uncertainties in coordinate + format; see `~astropy.nddata.Covariance`. + verbose : bool Print extra info. @@ -131,14 +135,27 @@ def spectrum_from_column_mapping(table, column_mapping, wcs=None, verbose=False) spec_kwargs.setdefault(kwarg_name, kwarg_val) # Ensure that the uncertainties are a subclass of NDUncertainty + if spec_kwargs.get('uncertainty') is None and covar is not None: + warnings.warn('Unable to parse uncertainty from provided table. Ignoring provided ' + 'covariance matrix data.') + covar = None if spec_kwargs.get('uncertainty') is not None: - spec_kwargs['uncertainty'] = StdDevUncertainty( - spec_kwargs.get('uncertainty')) + if covar is not None: + try: + spec_kwargs['uncertainty'] = Covariance.from_table(covar) + except (ValueError, TypeError): + warnings.warn('Unable to parse table with covariance data into into a Covariance ' + 'object. Ignoring covariance matrix data.') + covar = None + # NOTE: This is not an `else` or `elif` block in order to catch the + # change to covar=None when handling the exception above. + if covar is None: + spec_kwargs['uncertainty'] = StdDevUncertainty(spec_kwargs.get('uncertainty')) return Spectrum(**spec_kwargs, wcs=wcs, meta={'header': table.meta}) -def generic_spectrum_from_table(table, wcs=None): +def generic_spectrum_from_table(table, wcs=None, covar=None, **kwargs): """ Load spectrum from an Astropy table into a Spectrum object. Uses the following logic to figure out which column is which: @@ -158,9 +175,12 @@ def generic_spectrum_from_table(table, wcs=None): table : :class:`~astropy.table.Table` Table containing a column of ``flux``, and optionally ``spectral_axis`` and ``uncertainty`` as defined above. - wcs : :class:`~astropy.wcs.WCS` + wcs : :class:`~astropy.wcs.WCS`, optional A FITS WCS object. If this is present, the machinery will fall back - and default to using the ``wcs`` to find the dispersion information. + to using the ``wcs`` to find the dispersion information. + covar : `astropy.table.Table`, optional + Table providing a covariance matrix for the uncertainties in coordinate + format; see `~astropy.nddata.Covariance`. Returns ------- @@ -275,12 +295,27 @@ def _find_spectral_column(table, columns_to_search, spectral_axis): if table[err_column].ndim > 1: err = table[err_column].T elif flux.ndim > 1: # Repeat uncertainties over all flux columns + if covar is not None: + warnings.warn('When applying covariance, the dimensionality of the error ' + 'array must match the dimensionality of the flux array. Ignoring ' + 'covariance.') + covar = None err = np.tile(table[err_column], flux.shape[0], 1) else: err = table[err_column] - err = StdDevUncertainty(err.to(err.unit)) - if np.min(table[err_column]) <= 0.: - warnings.warn("Standard Deviation has values of 0 or less", AstropyUserWarning) + if covar is not None: + try: + err = Covariance.from_table(covar) + except (ValueError, TypeError): + warnings.warn('Unable to parse covariance table into a Covariance object. ' + 'Ignoring covariance matrix data.') + covar = None + # NOTE: This is not an `else` or `elif` block in order to catch the + # change to covar=None when handling the exception above. + if covar is None: + err = StdDevUncertainty(err.to(err.unit)) + if np.min(table[err_column]) <= 0.: + warnings.warn("Standard Deviation has values of 0 or less", AstropyUserWarning) else: err = None diff --git a/specutils/manipulation/extract_spectral_region.py b/specutils/manipulation/extract_spectral_region.py index a1a9aceac..536ee999e 100644 --- a/specutils/manipulation/extract_spectral_region.py +++ b/specutils/manipulation/extract_spectral_region.py @@ -193,6 +193,9 @@ def _get_joined_value(sps, key, unique_inds=None): uncert = sps[0].uncertainty if uncert is None: return None + if isinstance(uncert, Covariance): + raise NotImplementedError("Cannot yet combine spectral regions with " + "covariant uncertainties.") uncert._array = np.concatenate([sp.uncertainty._array for sp in sps]) return uncert[unique_inds] if unique_inds is not None else uncert elif key in concat_keys or key == 'spectral_axis': diff --git a/specutils/spectra/spectrum.py b/specutils/spectra/spectrum.py index 01d73dffd..bc849a19d 100644 --- a/specutils/spectra/spectrum.py +++ b/specutils/spectra/spectrum.py @@ -6,7 +6,7 @@ from astropy.coordinates import SpectralCoord from astropy.utils.decorators import lazyproperty from astropy.utils.decorators import deprecated -from astropy.nddata import NDUncertainty, NDIOMixin, NDArithmeticMixin, NDDataArray +from astropy.nddata import NDUncertainty, NDIOMixin, NDArithmeticMixin, NDDataArray, Covariance from gwcs.wcs import WCS as GWCS from .spectral_axis import SpectralAxis @@ -37,8 +37,10 @@ class Spectrum(OneDSpectrumMixin, NDCube, NDIOMixin, NDArithmeticMixin): Parameters ---------- flux : `~astropy.units.Quantity` - The flux data for this spectrum. This can be a simple `~astropy.units.Quantity`, - or an existing `~Spectrum` or `~ndcube.NDCube` object. + The flux data for this spectrum. This can be a simple + `~astropy.units.Quantity`, or an existing `~Spectrum` or + `~ndcube.NDCube` object. If an `~ndcube.NDCube` object, all other + arguments are ignored. spectral_axis : `~astropy.units.Quantity` or `~specutils.SpectralAxis` Dispersion information with the same shape as the last (or only) dimension of flux, or one greater than the last dimension of flux @@ -71,8 +73,10 @@ class Spectrum(OneDSpectrumMixin, NDCube, NDIOMixin, NDArithmeticMixin): values represent edges of the wavelength bin, or centers of the bin. uncertainty : `~astropy.nddata.NDUncertainty` Contains uncertainty information along with propagation rules for - spectrum arithmetic. Can take a unit, but if none is given, will use - the unit defined in the flux. + spectrum arithmetic. Can take a unit, but if none is given, will use the + unit defined in the flux. Note that functionality is limited for + `~astropy.nddata.Covariance` instances, particularly for + multidimensional data. mask : `~numpy.ndarray`-like Array where values in the flux to be masked are those that ``astype(bool)`` converts to True. (For example, integer arrays are not @@ -415,10 +419,13 @@ def __init__(self, flux=None, spectral_axis=None, spectral_axis_index=None, raise ValueError('Spectral axis must be strictly increasing or decreasing.') if hasattr(self, 'uncertainty') and self.uncertainty is not None: - if not flux.shape == self.uncertainty.array.shape: - raise ValueError( - "Flux axis ({}) and uncertainty ({}) shapes must be the " - "same.".format(flux.shape, self.uncertainty.array.shape)) + if isinstance(self.uncertainty, Covariance): + uncertainty_shape = self.uncertainty.data_shape + else: + uncertainty_shape = self.uncertainty.array.shape + if not flux.shape == uncertainty_shape: + raise ValueError(f"Flux axis ({flux.shape}) and uncertainty ({uncertainty_shape}) " + "shapes must be the same.") def __getitem__(self, item): """ @@ -498,11 +505,17 @@ def __getitem__(self, item): else: new_meta = deepcopy(self.meta) + if isinstance(self.uncertainty, Covariance): + new_unc = self.uncertainty.match_to_data_slice(item) + elif self.uncertainty is not None: + new_unc = self.uncertainty[item] + else: + new_unc = None + return self._copy( flux=self.flux[item], spectral_axis=self.spectral_axis[spec_item], - uncertainty=self.uncertainty[item] - if self.uncertainty is not None else None, + uncertainty=new_unc, mask=self.mask[item] if self.mask is not None else None, meta=new_meta, wcs=None, spectral_axis_index=new_spectral_axis_index) @@ -919,8 +932,13 @@ def __str__(self): # Add information about uncertainties if available if self.uncertainty: + _arr = ( + self.uncertainty.to_dense() + if isinstance(self.uncertainty, Covariance) + else self.uncertainty.array + ) result += (f'\nUncertainty={type(self.uncertainty).__name__} ' - f'({np.array2string(self.uncertainty.array, threshold=8)}' + f'({np.array2string(_arr, threshold=8)}' f' {self.uncertainty.unit})') return result diff --git a/specutils/tests/test_loaders.py b/specutils/tests/test_loaders.py index cb80200c8..934b373b7 100644 --- a/specutils/tests/test_loaders.py +++ b/specutils/tests/test_loaders.py @@ -6,6 +6,7 @@ import pytest import astropy.units as u import numpy as np +from scipy import sparse from astropy.io import fits from astropy.io.fits.verify import VerifyWarning from astropy.table import Table @@ -13,7 +14,7 @@ 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.nddata import StdDevUncertainty, InverseVariance, VarianceUncertainty, Covariance from astropy.tests.helper import quantity_allclose from astropy.utils.exceptions import AstropyUserWarning @@ -599,6 +600,48 @@ def test_tabular_fits_writer(tmp_path, spectral_axis): spectrum.uncertainty.quantity) +def test_tabular_fits_cov_io(tmp_path): + # Create a fake spectrum + wave = np.arange(4500., 5500., 1.) * u.AA + flux = np.full(len(wave), 10) * u.Jy + # And covariance + cov_diags = [ + np.ones(1000, dtype=float), + np.full(1000-1, 0.5, dtype=float), + np.full(1000-2, 0.2, dtype=float), + ] + cov = Covariance(array=sparse.diags(cov_diags, [0, 1, 2]), unit=u.Jy**2, assume_symmetric=True) + + # Create the Spectrum1D + spectrum = Spectrum1D(flux=flux, spectral_axis=wave, uncertainty=cov) + # Simple test of ingestion + assert np.array_equal(spectrum.uncertainty.to_dense(), cov.to_dense()), \ + 'Covariance not included correctly' + # Write it + tmpfile = str(tmp_path / '_tst.fits') + spectrum.write(tmpfile, format='tabular-fits') + + # Check the output file + with fits.open(tmpfile) as hdu: + assert len(hdu) == 3, 'Should contain 3 extensions, primary, DATA, COVAR' + assert np.array_equal([h.name for h in hdu], ['PRIMARY', 'DATA', 'COVAR']), \ + 'Extension names are wrong' + assert all([h.__class__.__name__ == 'BinTableHDU' for h in hdu[1:]]), \ + 'Data extensions should both be BinTableHDU' + assert len(hdu['COVAR'].data) == cov.stored_nnz, \ + 'Number of non-zero cov elements mismatch' + + # Read it + _spectrum = Spectrum1D.read(tmpfile) + assert spectrum.flux.unit == _spectrum.flux.unit + assert spectrum.spectral_axis.unit == _spectrum.spectral_axis.unit + assert quantity_allclose(spectrum.spectral_axis, _spectrum.spectral_axis) + assert quantity_allclose(spectrum.flux, _spectrum.flux) + assert np.array_equal(spectrum.uncertainty.to_dense(), _spectrum.uncertainty.to_dense()) + + # TODO: Add test for multidim + + @pytest.mark.parametrize("ndim", range(1, 4)) @pytest.mark.parametrize("spectral_axis", ['wavelength', 'frequency', 'energy', 'wavenumber'])