From f2154395ab0f95153f8163d7b7bdd902d99ad45e Mon Sep 17 00:00:00 2001 From: Clare Shanahan Date: Tue, 14 Nov 2023 16:31:32 -0500 Subject: [PATCH 1/5] . --- specutils/io/default_loaders/tabular_fits.py | 1 + 1 file changed, 1 insertion(+) diff --git a/specutils/io/default_loaders/tabular_fits.py b/specutils/io/default_loaders/tabular_fits.py index 9dc7ef73d..221a857bf 100644 --- a/specutils/io/default_loaders/tabular_fits.py +++ b/specutils/io/default_loaders/tabular_fits.py @@ -116,6 +116,7 @@ def tabular_fits_writer(spectrum, file_name, hdu=1, update_header=False, **kwarg raise ValueError(f'FITS does not support BINTABLE extension in HDU {hdu}.') header = spectrum.meta.get('header', fits.header.Header()).copy() + print('header', header) if update_header: hdr_types = (str, int, float, complex, bool, From 993d35b11f34a88c186b81ddd5448dd8a4439622 Mon Sep 17 00:00:00 2001 From: Clare Shanahan Date: Wed, 15 Nov 2023 17:42:09 -0500 Subject: [PATCH 2/5] very basic first commit --- .../io/default_loaders/mef_tabular_fits.py | 89 +++++++++++++++++++ 1 file changed, 89 insertions(+) create mode 100644 specutils/io/default_loaders/mef_tabular_fits.py diff --git a/specutils/io/default_loaders/mef_tabular_fits.py b/specutils/io/default_loaders/mef_tabular_fits.py new file mode 100644 index 000000000..0805bdc94 --- /dev/null +++ b/specutils/io/default_loaders/mef_tabular_fits.py @@ -0,0 +1,89 @@ +# a hack around tabular-fits that supports putting header info in the primary hdu + +from astropy.wcs import WCS +from specutils import Spectrum1D +import astropy.units as u +from astropy.io import fits +import numpy as np +from astropy.table import Table + +def mef_tabular_fits_writer(spectrum, file_name, hdu=0, update_header=False, **kwargs): + + # no matter what the dimensionality of 'flux', there will always be a primary header and + # a single BinTableHDU extension. the `hdu` arg controls where `meta.header` is dumped. + hdulist = [fits.PrimaryHDU(), None] + + if hdu > 1: + raise ValueError('`hdu`, which controls which extension `meta.header` is written to, must either be 0 or 1') + + # we expect anything that should be written out to the file header to be in `meta.header` + # This should be a `fits.header.Header` object, so convert to if meta.header is a dictionary + header = spectrum.meta.get('header', fits.header.Header()) # if no meta.header, create empty `fits.header.Header` + if not isinstance(header, fits.header.Header): + if isinstance(header, dict): + header = fits.header.Header(header) + else: + raise ValueError('`Spectrum1d.meta.header must be `fits.header.Header` or dictionary.') + + if update_header: + hdr_types = (str, int, float, complex, bool, + np.floating, np.integer, np.complexfloating, np.bool_) + header.update([keyword for keyword in spectrum.meta.items() if + isinstance(keyword[1], hdr_types)]) + + # Strip header of FITS reserved keywords + for keyword in ['NAXIS', 'NAXIS1', 'NAXIS2']: + header.remove(keyword, ignore_missing=True) + + # Add dispersion array and unit + wtype = kwargs.pop('wtype', spectrum.spectral_axis.dtype) + wunit = u.Unit(kwargs.pop('wunit', spectrum.spectral_axis.unit)) + + disp = spectrum.spectral_axis.to(wunit, equivalencies=u.spectral()) + + # Mapping of spectral_axis types to header TTYPE1 (no "torque/work" types!) + dispname = str(wunit.physical_type) + if dispname == "length": + dispname = "wavelength" + elif "energy" in dispname: + dispname = "energy" + + # Add flux array and unit + ftype = kwargs.pop('ftype', spectrum.flux.dtype) + funit = u.Unit(kwargs.pop('funit', spectrum.flux.unit)) + flux = spectrum.flux.to(funit, equivalencies=u.spectral_density(disp)) + + columns = [disp.astype(wtype), flux.astype(ftype)] + colnames = [dispname, "flux"] + + # Include uncertainty - units to be inferred from spectrum.flux + 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)) + colnames.append("uncertainty") + except RuntimeWarning: + raise ValueError("Could not convert uncertainty to StdDevUncertainty due" + " to divide-by-zero error.") + + # For > 1D data transpose from row-major format + for c in range(1, len(columns)): + if columns[c].ndim > 1: + columns[c] = columns[c].T + + tab = Table(columns, names=colnames, meta=header) + hdulist[1] = fits.BinTableHDU(tab) + + # now figure out where meta.header (and if specified, addl meta keys) should go + if hdu==0: + hdulist[0].header = header + else: # must otherwise be 1, we already checked this + hdulist[1].header = header + + hdulist = fits.HDUList(hdulist) \ No newline at end of file From c08d2704976ce20781b688c83d4c58504778387b Mon Sep 17 00:00:00 2001 From: Clare Shanahan Date: Wed, 15 Nov 2023 17:52:21 -0500 Subject: [PATCH 3/5] . --- specutils/io/default_loaders/mef_tabular_fits.py | 4 +++- specutils/io/default_loaders/tabular_fits.py | 1 - 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/specutils/io/default_loaders/mef_tabular_fits.py b/specutils/io/default_loaders/mef_tabular_fits.py index 0805bdc94..0f9957556 100644 --- a/specutils/io/default_loaders/mef_tabular_fits.py +++ b/specutils/io/default_loaders/mef_tabular_fits.py @@ -86,4 +86,6 @@ def mef_tabular_fits_writer(spectrum, file_name, hdu=0, update_header=False, **k else: # must otherwise be 1, we already checked this hdulist[1].header = header - hdulist = fits.HDUList(hdulist) \ No newline at end of file + hdulist = fits.HDUList(hdulist) + + hdulist.writeto(file_name) \ No newline at end of file diff --git a/specutils/io/default_loaders/tabular_fits.py b/specutils/io/default_loaders/tabular_fits.py index 221a857bf..9dc7ef73d 100644 --- a/specutils/io/default_loaders/tabular_fits.py +++ b/specutils/io/default_loaders/tabular_fits.py @@ -116,7 +116,6 @@ def tabular_fits_writer(spectrum, file_name, hdu=1, update_header=False, **kwarg raise ValueError(f'FITS does not support BINTABLE extension in HDU {hdu}.') header = spectrum.meta.get('header', fits.header.Header()).copy() - print('header', header) if update_header: hdr_types = (str, int, float, complex, bool, From a980b55e0c939f2aa7ba6f39217efe5b43cd1a7a Mon Sep 17 00:00:00 2001 From: Clare Shanahan Date: Wed, 15 Nov 2023 18:11:29 -0500 Subject: [PATCH 4/5] . --- .../io/default_loaders/mef_tabular_fits.py | 169 +++++++++--------- 1 file changed, 87 insertions(+), 82 deletions(-) diff --git a/specutils/io/default_loaders/mef_tabular_fits.py b/specutils/io/default_loaders/mef_tabular_fits.py index 0f9957556..33dca8afe 100644 --- a/specutils/io/default_loaders/mef_tabular_fits.py +++ b/specutils/io/default_loaders/mef_tabular_fits.py @@ -7,85 +7,90 @@ import numpy as np from astropy.table import Table -def mef_tabular_fits_writer(spectrum, file_name, hdu=0, update_header=False, **kwargs): - - # no matter what the dimensionality of 'flux', there will always be a primary header and - # a single BinTableHDU extension. the `hdu` arg controls where `meta.header` is dumped. - hdulist = [fits.PrimaryHDU(), None] - - if hdu > 1: - raise ValueError('`hdu`, which controls which extension `meta.header` is written to, must either be 0 or 1') - - # we expect anything that should be written out to the file header to be in `meta.header` - # This should be a `fits.header.Header` object, so convert to if meta.header is a dictionary - header = spectrum.meta.get('header', fits.header.Header()) # if no meta.header, create empty `fits.header.Header` - if not isinstance(header, fits.header.Header): - if isinstance(header, dict): - header = fits.header.Header(header) - else: - raise ValueError('`Spectrum1d.meta.header must be `fits.header.Header` or dictionary.') - - if update_header: - hdr_types = (str, int, float, complex, bool, - np.floating, np.integer, np.complexfloating, np.bool_) - header.update([keyword for keyword in spectrum.meta.items() if - isinstance(keyword[1], hdr_types)]) - - # Strip header of FITS reserved keywords - for keyword in ['NAXIS', 'NAXIS1', 'NAXIS2']: - header.remove(keyword, ignore_missing=True) - - # Add dispersion array and unit - wtype = kwargs.pop('wtype', spectrum.spectral_axis.dtype) - wunit = u.Unit(kwargs.pop('wunit', spectrum.spectral_axis.unit)) - - disp = spectrum.spectral_axis.to(wunit, equivalencies=u.spectral()) - - # Mapping of spectral_axis types to header TTYPE1 (no "torque/work" types!) - dispname = str(wunit.physical_type) - if dispname == "length": - dispname = "wavelength" - elif "energy" in dispname: - dispname = "energy" - - # Add flux array and unit - ftype = kwargs.pop('ftype', spectrum.flux.dtype) - funit = u.Unit(kwargs.pop('funit', spectrum.flux.unit)) - flux = spectrum.flux.to(funit, equivalencies=u.spectral_density(disp)) - - columns = [disp.astype(wtype), flux.astype(ftype)] - colnames = [dispname, "flux"] - - # Include uncertainty - units to be inferred from spectrum.flux - 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)) - colnames.append("uncertainty") - except RuntimeWarning: - raise ValueError("Could not convert uncertainty to StdDevUncertainty due" - " to divide-by-zero error.") - - # For > 1D data transpose from row-major format - for c in range(1, len(columns)): - if columns[c].ndim > 1: - columns[c] = columns[c].T - - tab = Table(columns, names=colnames, meta=header) - hdulist[1] = fits.BinTableHDU(tab) - - # now figure out where meta.header (and if specified, addl meta keys) should go - if hdu==0: - hdulist[0].header = header - else: # must otherwise be 1, we already checked this - hdulist[1].header = header - - hdulist = fits.HDUList(hdulist) - - hdulist.writeto(file_name) \ No newline at end of file +def mef_tabular_fits_writer(spectrum, file_name='my_func_output.fits', hdu=0, update_header=False, **kwargs): + + # no matter what the dimensionality of 'flux', there will always be a primary header and + # a single BinTableHDU extension. the `hdu` arg controls where `meta.header` is dumped. + hdulist = [fits.PrimaryHDU(), None] + pri_header_initial = hdulist[0].header + + if hdu > 1: + raise ValueError('`hdu`, which controls which extension `meta.header` is written to, must either be 0 or 1') + + # we expect anything that should be written out to the file header to be in `meta.header` + # This should be a `fits.header.Header` object, so convert to if meta.header is a dictionary + header = spectrum.meta.get('header', fits.header.Header()) # if no meta.header, create empty `fits.header.Header` + + if hdu == 0: + header.update(pri_header_initial) # update with initial values created from primary header + + if not isinstance(header, fits.header.Header): + if isinstance(header, dict): + header = fits.header.Header(header) + else: + raise ValueError('`Spectrum1d.meta.header must be `fits.header.Header` or dictionary.') + + if update_header: + hdr_types = (str, int, float, complex, bool, + np.floating, np.integer, np.complexfloating, np.bool_) + header.update([keyword for keyword in spectrum.meta.items() if + isinstance(keyword[1], hdr_types)]) + + # Strip header of FITS reserved keywords + for keyword in ['NAXIS', 'NAXIS1', 'NAXIS2']: + header.remove(keyword, ignore_missing=True) + + # Add dispersion array and unit + wtype = kwargs.pop('wtype', spectrum.spectral_axis.dtype) + wunit = u.Unit(kwargs.pop('wunit', spectrum.spectral_axis.unit)) + + disp = spectrum.spectral_axis.to(wunit, equivalencies=u.spectral()) + + # Mapping of spectral_axis types to header TTYPE1 (no "torque/work" types!) + dispname = str(wunit.physical_type) + if dispname == "length": + dispname = "wavelength" + elif "energy" in dispname: + dispname = "energy" + + # Add flux array and unit + ftype = kwargs.pop('ftype', spectrum.flux.dtype) + funit = u.Unit(kwargs.pop('funit', spectrum.flux.unit)) + flux = spectrum.flux.to(funit, equivalencies=u.spectral_density(disp)) + + columns = [disp.astype(wtype), flux.astype(ftype)] + colnames = [dispname, "flux"] + + # Include uncertainty - units to be inferred from spectrum.flux + 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)) + colnames.append("uncertainty") + except RuntimeWarning: + raise ValueError("Could not convert uncertainty to StdDevUncertainty due" + " to divide-by-zero error.") + + # For > 1D data transpose from row-major format + for c in range(1, len(columns)): + if columns[c].ndim > 1: + columns[c] = columns[c].T + + tab = Table(columns, names=colnames, meta=header) + hdulist[1] = fits.BinTableHDU(tab) + + # now figure out where meta.header (and if specified, addl meta keys) should go + if hdu==0: + hdulist[0].header = header + else: # must otherwise be 1, we already checked this + hdulist[1].header = header + + hdulist = fits.HDUList(hdulist) + + hdulist.writeto(file_name) \ No newline at end of file From 1b9ccf5c34e679c5a3fca575c2e55aa7b14bf1ac Mon Sep 17 00:00:00 2001 From: Clare Shanahan Date: Wed, 15 Nov 2023 18:30:25 -0500 Subject: [PATCH 5/5] . --- .../io/default_loaders/mef_tabular_fits.py | 129 +++++++++--------- 1 file changed, 65 insertions(+), 64 deletions(-) diff --git a/specutils/io/default_loaders/mef_tabular_fits.py b/specutils/io/default_loaders/mef_tabular_fits.py index 33dca8afe..8efe52661 100644 --- a/specutils/io/default_loaders/mef_tabular_fits.py +++ b/specutils/io/default_loaders/mef_tabular_fits.py @@ -30,67 +30,68 @@ def mef_tabular_fits_writer(spectrum, file_name='my_func_output.fits', hdu=0, up else: raise ValueError('`Spectrum1d.meta.header must be `fits.header.Header` or dictionary.') - if update_header: - hdr_types = (str, int, float, complex, bool, - np.floating, np.integer, np.complexfloating, np.bool_) - header.update([keyword for keyword in spectrum.meta.items() if - isinstance(keyword[1], hdr_types)]) - - # Strip header of FITS reserved keywords - for keyword in ['NAXIS', 'NAXIS1', 'NAXIS2']: - header.remove(keyword, ignore_missing=True) - - # Add dispersion array and unit - wtype = kwargs.pop('wtype', spectrum.spectral_axis.dtype) - wunit = u.Unit(kwargs.pop('wunit', spectrum.spectral_axis.unit)) - - disp = spectrum.spectral_axis.to(wunit, equivalencies=u.spectral()) - - # Mapping of spectral_axis types to header TTYPE1 (no "torque/work" types!) - dispname = str(wunit.physical_type) - if dispname == "length": - dispname = "wavelength" - elif "energy" in dispname: - dispname = "energy" - - # Add flux array and unit - ftype = kwargs.pop('ftype', spectrum.flux.dtype) - funit = u.Unit(kwargs.pop('funit', spectrum.flux.unit)) - flux = spectrum.flux.to(funit, equivalencies=u.spectral_density(disp)) - - columns = [disp.astype(wtype), flux.astype(ftype)] - colnames = [dispname, "flux"] - - # Include uncertainty - units to be inferred from spectrum.flux - 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)) - colnames.append("uncertainty") - except RuntimeWarning: - raise ValueError("Could not convert uncertainty to StdDevUncertainty due" - " to divide-by-zero error.") - - # For > 1D data transpose from row-major format - for c in range(1, len(columns)): - if columns[c].ndim > 1: - columns[c] = columns[c].T - - tab = Table(columns, names=colnames, meta=header) - hdulist[1] = fits.BinTableHDU(tab) - - # now figure out where meta.header (and if specified, addl meta keys) should go - if hdu==0: - hdulist[0].header = header - else: # must otherwise be 1, we already checked this - hdulist[1].header = header - - hdulist = fits.HDUList(hdulist) - - hdulist.writeto(file_name) \ No newline at end of file + if update_header: + hdr_types = (str, int, float, complex, bool, + np.floating, np.integer, np.complexfloating, np.bool_) + header.update([keyword for keyword in spectrum.meta.items() if + isinstance(keyword[1], hdr_types)]) + + # Strip header of FITS reserved keywords + for keyword in ['NAXIS', 'NAXIS1', 'NAXIS2']: + header.remove(keyword, ignore_missing=True) + + # Add dispersion array and unit + wtype = kwargs.pop('wtype', spectrum.spectral_axis.dtype) + wunit = u.Unit(kwargs.pop('wunit', spectrum.spectral_axis.unit)) + + disp = spectrum.spectral_axis.to(wunit, equivalencies=u.spectral()) + + # Mapping of spectral_axis types to header TTYPE1 (no "torque/work" types!) + dispname = str(wunit.physical_type) + if dispname == "length": + dispname = "wavelength" + elif "energy" in dispname: + dispname = "energy" + + # Add flux array and unit + ftype = kwargs.pop('ftype', spectrum.flux.dtype) + funit = u.Unit(kwargs.pop('funit', spectrum.flux.unit)) + flux = spectrum.flux.to(funit, equivalencies=u.spectral_density(disp)) + + columns = [disp.astype(wtype), flux.astype(ftype)] + colnames = [dispname, "flux"] + + # Include uncertainty - units to be inferred from spectrum.flux + 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)) + colnames.append("uncertainty") + except RuntimeWarning: + raise ValueError("Could not convert uncertainty to StdDevUncertainty due" + " to divide-by-zero error.") + + # For > 1D data transpose from row-major format + for c in range(1, len(columns)): + if columns[c].ndim > 1: + columns[c] = columns[c].T + + tab = Table(columns, names=colnames) + hdulist[1] = fits.BinTableHDU(tab) + + # now figure out where meta.header (and if specified, addl meta keys) should go + if hdu==0: + print('here') + hdulist[0].header.update(header) + else: # must otherwise be 1, we already checked this + hdulist[1].header.update(header) + + hdulist = fits.HDUList(hdulist) + + hdulist.writeto(file_name, overwrite=True) \ No newline at end of file