From b6851a6dc604e4db380e6a3af105d4ecb11ad4c7 Mon Sep 17 00:00:00 2001 From: Riley Thai Date: Sat, 12 Oct 2024 09:44:52 +1100 Subject: [PATCH 01/16] fix: mwmVisit BOSS HDU default loader fail - no longer checks for "date_obs"; calculate that yourself - also adds "sdss_id" to metadata now --- specutils/io/default_loaders/sdss_v.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/specutils/io/default_loaders/sdss_v.py b/specutils/io/default_loaders/sdss_v.py index 514432978..dc5f4a0e6 100644 --- a/specutils/io/default_loaders/sdss_v.py +++ b/specutils/io/default_loaders/sdss_v.py @@ -312,6 +312,7 @@ def load_sdss_apVisit_list(file_obj, **kwargs): # BOSS REDUX products (specLite, specFull, custom coadd files, etc) + @data_loader( "SDSS-V spec", identifier=spec_sdss5_identify, @@ -578,16 +579,16 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): meta["snr"] = np.array(hdulist[hdu].data["snr"]) # Add identifiers (obj, telescope, mjd, datatype) - # TODO: need to see what metadata we're interested in for the MWM files. meta["telescope"] = hdulist[hdu].data["telescope"] meta["instrument"] = hdulist[hdu].header.get("INSTRMNT") - try: + try: # get obj if exists meta["obj"] = hdulist[hdu].data["obj"] except KeyError: pass + + # choose between mwmVisit/Star via KeyError except try: - meta["date"] = hdulist[hdu].data["date_obs"] - meta["mjd"] = hdulist[hdu].data["mjd"] + meta['mjd'] = hdulist[hdu].data['mjd'] meta["datatype"] = "mwmVisit" except KeyError: meta["mjd"] = (str(hdulist[hdu].data["min_mjd"][0]) + "->" + @@ -595,6 +596,7 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): meta["datatype"] = "mwmStar" finally: meta["name"] = hdulist[hdu].name + meta["sdss_id"] = hdulist[hdu].data['sdss_id'] return Spectrum1D( spectral_axis=spectral_axis, From 51e8a2e9e3eb3c9900eee834ddbed3510024f603 Mon Sep 17 00:00:00 2001 From: Riley Thai Date: Sat, 12 Oct 2024 11:39:51 +1100 Subject: [PATCH 02/16] fix: SDSS-V SpectrumList loader ambiguity + add: BOSS-only mwm test cases - added new test cases for BOSS-only mwmVisit and mwmStar files - added new checks to SpectrumList mwmVisit/mwmStar test to check verified filetype is correct - forced override on default SpectrumList loaders -- now SpectrumList is no longer ambiguous and doesn't require a format specification - relevant areas in tests are updated accordingly - added print warnings to when HDU is not specified on Spectrum1D loaders for files with multiple spectra. - ensured tests now remove tempfiles with os.remove - arguably, this could work better with tmpfile, but i don't know how tests are deployed on the server-side --- specutils/io/default_loaders/sdss_v.py | 15 ++- .../io/default_loaders/tests/test_sdss_v.py | 106 +++++++++--------- 2 files changed, 62 insertions(+), 59 deletions(-) diff --git a/specutils/io/default_loaders/sdss_v.py b/specutils/io/default_loaders/sdss_v.py index dc5f4a0e6..e3195d31f 100644 --- a/specutils/io/default_loaders/sdss_v.py +++ b/specutils/io/default_loaders/sdss_v.py @@ -183,9 +183,10 @@ def load_sdss_apStar_1D(file_obj, idx: int = 0, **kwargs): @data_loader( - "SDSS-V apStar multi", + "SDSS-V apStar", identifier=apStar_identify, dtype=SpectrumList, + force=True, priority=10, extensions=["fits"], ) @@ -259,9 +260,10 @@ def load_sdss_apVisit_1D(file_obj, **kwargs): @data_loader( - "SDSS-V apVisit multi", + "SDSS-V apVisit", identifier=apVisit_identify, dtype=SpectrumList, + force=True, priority=10, extensions=["fits"], ) @@ -338,6 +340,7 @@ def load_sdss_spec_1D(file_obj, *args, hdu: Optional[int] = None, **kwargs): """ if hdu is None: # TODO: how should we handle this -- multiple things in file, but the user cannot choose. + print('HDU not specified. Loading coadd spectrum (HDU1)') hdu = 1 # defaulting to coadd # raise ValueError("HDU not specified! Please specify a HDU to load.") elif hdu in [2, 3, 4]: @@ -348,9 +351,10 @@ def load_sdss_spec_1D(file_obj, *args, hdu: Optional[int] = None, **kwargs): @data_loader( - "SDSS-V spec multi", + "SDSS-V spec", identifier=spec_sdss5_identify, dtype=SpectrumList, + force=True, priority=5, extensions=["fits"], ) @@ -463,14 +467,17 @@ def load_sdss_mwm_1d(file_obj, hdu: Optional[int] = None, **kwargs): for i in range(len(hdulist)): if hdulist[i].header.get("DATASUM") != "0": hdu = i + print('HDU not specified. Loading spectrum at (HDU{})'. + format(i)) break return _load_mwmVisit_or_mwmStar_hdu(hdulist, hdu, **kwargs) @data_loader( - "SDSS-V mwm multi", + "SDSS-V mwm", identifier=mwm_identify, + force=True, dtype=SpectrumList, priority=20, extensions=["fits"], diff --git a/specutils/io/default_loaders/tests/test_sdss_v.py b/specutils/io/default_loaders/tests/test_sdss_v.py index f70490482..15f9916cb 100644 --- a/specutils/io/default_loaders/tests/test_sdss_v.py +++ b/specutils/io/default_loaders/tests/test_sdss_v.py @@ -1,8 +1,9 @@ +import os + import numpy as np import pytest - from astropy.io import fits -from astropy.units import Unit, Angstrom +from astropy.units import Angstrom, Unit from specutils import Spectrum1D, SpectrumList @@ -25,34 +26,7 @@ def generate_apogee_hdu(observatory="APO", with_wl=True, datasum="0"): fits.Column(name="apred", array=[b"1.2"], format="3A"), fits.Column(name="obj", array=[b"2M19534321+6705175"], format="18A"), fits.Column(name="telescope", array=[b"apo25m"], format="6A"), - fits.Column(name="min_mjd", array=[59804], format="J"), - fits.Column(name="max_mjd", array=[59866], format="J"), - fits.Column(name="n_entries", array=[-1], format="J"), - fits.Column(name="n_visits", array=[5], format="J"), - fits.Column(name="n_good_visits", array=[5], format="J"), - fits.Column(name="n_good_rvs", array=[5], format="J"), - fits.Column(name="snr", array=[46.56802], format="E"), - fits.Column(name="mean_fiber", array=[256.0], format="E"), - fits.Column(name="std_fiber", array=[0.0], format="E"), - fits.Column(name="spectrum_flags", array=[1048576], format="J"), - fits.Column(name="v_rad", array=[-56.7284381], format="E"), - fits.Column(name="e_v_rad", array=[5.35407624], format="E"), - fits.Column(name="std_v_rad", array=[10.79173857], format="E"), - fits.Column(name="median_e_v_rad", array=[16.19418386], format="E"), - fits.Column(name="doppler_teff", array=[7169.0107], format="E"), - fits.Column(name="doppler_e_teff", array=[9.405238], format="E"), - fits.Column(name="doppler_logg", array=[2.981389], format="E"), - fits.Column(name="doppler_e_logg", array=[0.01916536], format="E"), - fits.Column(name="doppler_fe_h", array=[-1.20532212], format="E"), - fits.Column(name="doppler_e_fe_h", array=[0.0093738], format="E"), - fits.Column(name="doppler_rchi2", array=[1.1424173], format="E"), - fits.Column(name="doppler_flags", array=[0], format="J"), - fits.Column(name="xcorr_v_rad", array=[np.nan], format="E"), - fits.Column(name="xcorr_v_rel", array=[np.nan], format="E"), - fits.Column(name="xcorr_e_v_rel", array=[np.nan], format="E"), - fits.Column(name="ccfwhm", array=[np.nan], format="E"), - fits.Column(name="autofwhm", array=[np.nan], format="E"), - fits.Column(name="n_components", array=[1], format="J"), + fits.Column(name="snr", array=[50], format="E"), ] if with_wl: columns.append( @@ -60,6 +34,14 @@ def generate_apogee_hdu(observatory="APO", with_wl=True, datasum="0"): array=wl, format="8575E", dim="(8575)")) + columns += [ + fits.Column(name="min_mjd", array=[59804], format="J"), + fits.Column(name="max_mjd", array=[59866], format="J"), + ] + else: + columns += [ + fits.Column(name="mjd", array=[59804], format="J"), + ] columns += [ fits.Column(name="flux", array=flux, format="8575E", dim="(8575)"), fits.Column(name="ivar", array=ivar, format="8575E", dim="(8575)"), @@ -110,32 +92,23 @@ def generate_boss_hdu(observatory="APO", with_wl=True, datasum="0"): fits.Column(name="sdss_id", array=[42], format="K"), fits.Column(name="run2d", array=["6_1_2"], format="6A"), fits.Column(name="telescope", array=["apo25m"], format="6A"), - fits.Column(name="min_mjd", array=[54], format="J"), - fits.Column(name="max_mjd", array=[488], format="J"), - fits.Column(name="n_visits", array=[1], format="J"), - fits.Column(name="n_good_visits", array=[1], format="J"), - fits.Column(name="n_good_rvs", array=[1], format="J"), - fits.Column(name="v_rad", array=[0], format="E"), - fits.Column(name="e_v_rad", array=[1], format="E"), - fits.Column(name="std_v_rad", array=[1], format="E"), - fits.Column(name="median_e_v_rad", array=[3], format="E"), - fits.Column(name="xcsao_teff", array=[5000], format="E"), - fits.Column(name="xcsao_e_teff", array=[10], format="E"), - fits.Column(name="xcsao_logg", array=[4], format="E"), - fits.Column(name="xcsao_e_logg", array=[3], format="E"), - fits.Column(name="xcsao_fe_h", array=[0], format="E"), - fits.Column(name="xcsao_e_fe_h", array=[5], format="E"), - fits.Column(name="xcsao_meanrxc", array=[0], format="E"), fits.Column(name="snr", array=[50], format="E"), - fits.Column(name="gri_gaia_transform_flags", array=[1], format="J"), - fits.Column(name="zwarning_flags", array=[0], format="J"), ] + if with_wl: columns.append( fits.Column(name="wavelength", array=wl, format="4648E", dim="(4648)")) + columns += [ + fits.Column(name="min_mjd", array=[54], format="J"), + fits.Column(name="max_mjd", array=[488], format="J"), + ] + else: + columns += [ + fits.Column(name="mjd", array=[59804], format="J"), + ] columns += [ fits.Column(name="flux", array=flux, format="4648E", dim="(4648)"), fits.Column(name="ivar", array=ivar, format="4648E", dim="(4648)"), @@ -483,13 +456,20 @@ def test_mwm_1d(file_obj, hdu, with_wl, hduflags): assert len(data.flux.value) == length assert data.spectral_axis.unit == Angstrom assert data.flux.unit == Unit("1e-17 erg / (s cm2 Angstrom)") + os.remove(tmpfile) @pytest.mark.parametrize( "file_obj, with_wl, hduflags", [ + ("mwm-temp", False, [0, 0, 1, 1]), + ("mwm-temp", True, [0, 0, 1, 1]), + ("mwm-temp", False, [0, 1, 1, 0]), + ("mwm-temp", True, [0, 1, 1, 0]), + ("mwm-temp", False, [1, 1, 0, 0]), + ("mwm-temp", True, [1, 1, 0, 0]), ("mwm-temp", False, [1, 1, 1, 1]), - ("mwm-temp", True, [0, 1, 0, 1]), + ("mwm-temp", True, [1, 1, 1, 1]), ], ) def test_mwm_list(file_obj, with_wl, hduflags): @@ -497,7 +477,7 @@ def test_mwm_list(file_obj, with_wl, hduflags): tmpfile = str(file_obj) + ".fits" mwm_HDUList(hduflags, with_wl).writeto(tmpfile, overwrite=True) - data = SpectrumList.read(tmpfile, format="SDSS-V mwm multi") + data = SpectrumList.read(tmpfile) assert isinstance(data, SpectrumList) for i in range(len(data)): assert isinstance(data[i], Spectrum1D) @@ -509,10 +489,15 @@ def test_mwm_list(file_obj, with_wl, hduflags): else: raise ValueError( "INSTRMNT tag in test HDU header is not set properly.") + if with_wl: + assert data[i].meta['datatype'].lower() == 'mwmstar' + else: + assert data[i].meta['datatype'].lower() == 'mwmvisit' assert len(data[i].spectral_axis.value) == length assert len(data[i].flux.value) == length assert data[i].spectral_axis.unit == Angstrom assert data[i].flux.unit == Unit("1e-17 erg / (s cm2 Angstrom)") + os.remove(tmpfile) @pytest.mark.parametrize( @@ -530,6 +515,7 @@ def test_mwm_1d_fail_spec(file_obj, hdu, hduflags): mwm_HDUList(hduflags, True).writeto(tmpfile, overwrite=True) with pytest.raises(IndexError): Spectrum1D.read(tmpfile, hdu=hdu) + os.remove(tmpfile) @pytest.mark.parametrize( @@ -546,6 +532,7 @@ def test_mwm_1d_fail(file_obj, with_wl): with pytest.raises(ValueError): Spectrum1D.read(tmpfile) + os.remove(tmpfile) @pytest.mark.parametrize( @@ -561,7 +548,8 @@ def test_mwm_list_fail(file_obj, with_wl): mwm_HDUList([0, 0, 0, 0], with_wl).writeto(tmpfile, overwrite=True) with pytest.raises(ValueError): - SpectrumList.read(tmpfile, format="SDSS-V mwm multi") + SpectrumList.read(tmpfile) + os.remove(tmpfile) @pytest.mark.parametrize( @@ -589,6 +577,7 @@ def test_spec_1d(file_obj, n_spectra): assert len(data.mask) == 10 assert data[i].meta["header"].get("foobar") == "barfoo" + os.remove(tmpfile) @pytest.mark.parametrize( @@ -603,7 +592,7 @@ def test_spec_list(file_obj, n_spectra): tmpfile = str(file_obj) + ".fits" spec_HDUList(n_spectra).writeto(tmpfile, overwrite=True) - data = SpectrumList.read(tmpfile, format="SDSS-V spec multi") + data = SpectrumList.read(tmpfile) assert isinstance(data, SpectrumList) assert len(data) == n_spectra + 1 for i in range(n_spectra): @@ -613,6 +602,7 @@ def test_spec_list(file_obj, n_spectra): assert data[i].spectral_axis.unit == Angstrom assert len(data[i].mask) == 10 assert data[i].meta["header"].get("foobar") == "barfoo" + os.remove(tmpfile) @pytest.mark.parametrize( @@ -631,6 +621,7 @@ def test_spec_1d_fail_hdu(file_obj, hdu): with pytest.raises(ValueError): Spectrum1D.read(tmpfile, hdu=hdu) + os.remove(tmpfile) @pytest.mark.parametrize( @@ -653,6 +644,7 @@ def test_apStar_1D(file_obj, idx): assert data.spectral_axis.unit == Angstrom assert data.meta["header"].get("foobar") == "barfoo" + os.remove(tmpfile) @pytest.mark.parametrize( @@ -666,7 +658,7 @@ def test_apStar_list(file_obj, n_spectra): tmpfile = str(file_obj) + ".fits" apStar_HDUList(n_spectra).writeto(tmpfile, overwrite=True) - data = SpectrumList.read(tmpfile, format="SDSS-V apStar multi") + data = SpectrumList.read(tmpfile) assert isinstance(data, SpectrumList) assert len(data) == n_spectra for i in range(len(data)): @@ -675,6 +667,7 @@ def test_apStar_list(file_obj, n_spectra): assert data[i].flux.unit == Unit("1e-17 erg / (s cm2 Angstrom)") assert len(data[i].spectral_axis.value) == 10 assert data[i].meta["header"].get("foobar") == "barfoo" + os.remove(tmpfile) @pytest.mark.parametrize( @@ -689,7 +682,8 @@ def test_apStar_fail_list(file_obj): apStar_HDUList(1).writeto(tmpfile, overwrite=True) with pytest.raises(ValueError): - SpectrumList.read(tmpfile, format="SDSS-V apStar multi") + SpectrumList.read(tmpfile) + os.remove(tmpfile) @pytest.mark.parametrize( @@ -707,6 +701,7 @@ def test_apVisit_1D(file_obj): assert np.array_equal(data.spectral_axis.value, np.arange(1, 31, 1)) assert len(data.flux.value) == 30 assert data.meta["header"].get("foobar") == "barfoo" + os.remove(tmpfile) @pytest.mark.parametrize( @@ -719,6 +714,7 @@ def test_apVisit_list(file_obj): tmpfile = str(file_obj) + ".fits" apVisit_HDUList().writeto(tmpfile, overwrite=True) - data = SpectrumList.read(tmpfile, format="SDSS-V apVisit multi") + data = SpectrumList.read(tmpfile) assert isinstance(data, SpectrumList) assert len(data) == 3 + os.remove(tmpfile) From 4bee13694d22e1a5da5d6e224a0b1c776849b60a Mon Sep 17 00:00:00 2001 From: Riley Thai Date: Sun, 13 Oct 2024 09:28:59 +1100 Subject: [PATCH 03/16] add: changelog -> CHANGES.rst - three points outlining changes listed in PR as per #1185 --- CHANGES.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index 7b0dd3146..af3ca4290 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -10,12 +10,18 @@ Bug Fixes - Fixed ``Spectrum1D.with_flux_unit()`` not converting uncertainty along with flux unit. [#1181] +- Fixed ``mwmVisit`` SDSS-V ``Spectrum1D`` and ``SpectrumList`` default loader being unable to load files containing only BOSS instrument spectra. [#1185] + +- Fixed automatic format detection for SDSS-V ``SpectrumList`` default loaders. [#1185] + Other Changes and Additions ^^^^^^^^^^^^^^^^^^^^^^^^^^^ - Replaced ``LevMarLSQFitter`` with ``TRFLSQFitter`` as the former is no longer recommended by ``astropy``. [#1180] +- "Multi" loaders have been removed from SDSS-V ``SpectrumList`` default loaders. [#1185] + 1.17.0 (2024-10-04) ------------------- From a859dd6c9b540b6f651a3b0d882be6b418aecde7 Mon Sep 17 00:00:00 2001 From: "P. L. Lim" <2090236+pllim@users.noreply.github.com> Date: Mon, 14 Oct 2024 10:28:49 -0400 Subject: [PATCH 04/16] FEAT: with_spectral_axis_and_flux_units (#1184) --- CHANGES.rst | 3 + specutils/spectra/spectrum_mixin.py | 71 ++++++++------ specutils/tests/test_spectrum1d.py | 137 +++++++++++++++------------- 3 files changed, 121 insertions(+), 90 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 7b0dd3146..d2a0ebb60 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -4,6 +4,9 @@ New Features ^^^^^^^^^^^^ +- New ``Spectrum1D.with_spectral_axis_and_flux_units`` method to convert both + spectral axis and flux units at the same time. [#1184] + Bug Fixes ^^^^^^^^^ diff --git a/specutils/spectra/spectrum_mixin.py b/specutils/spectra/spectrum_mixin.py index 6479c0df5..06fa32123 100644 --- a/specutils/spectra/spectrum_mixin.py +++ b/specutils/spectra/spectrum_mixin.py @@ -86,6 +86,26 @@ def new_flux_unit(self, unit, equivalencies=None, suppress_conversion=False): return self.with_flux_unit(unit, equivalencies=equivalencies, suppress_conversion=suppress_conversion) + def _convert_flux(self, unit, equivalencies=None, suppress_conversion=False): + """This is always done in-place. + Also see :meth:`with_flux_unit`.""" + + if not suppress_conversion: + if equivalencies is None: + equivalencies = eq.spectral_density(self.spectral_axis) + + new_data = self.flux.to(unit, equivalencies=equivalencies) + + self._data = new_data.value + self._unit = new_data.unit + else: + self._unit = u.Unit(unit) + + if self.uncertainty is not None: + self.uncertainty = StdDevUncertainty( + self.uncertainty.represent_as(StdDevUncertainty).quantity.to( + unit, equivalencies=equivalencies)) + def with_flux_unit(self, unit, equivalencies=None, suppress_conversion=False): """Returns a new spectrum with a different flux unit. If uncertainty is defined, it will be converted to @@ -107,30 +127,14 @@ def with_flux_unit(self, unit, equivalencies=None, suppress_conversion=False): Returns ------- - `~specutils.Spectrum1D` + new_spec : `~specutils.Spectrum1D` A new spectrum with the converted flux array (and uncertainty, if applicable). """ new_spec = deepcopy(self) - - if not suppress_conversion: - if equivalencies is None: - equivalencies = eq.spectral_density(self.spectral_axis) - - new_data = self.flux.to( - unit, equivalencies=equivalencies) - - new_spec._data = new_data.value - new_spec._unit = new_data.unit - else: - new_spec._unit = u.Unit(unit) - - if self.uncertainty is not None: - new_spec.uncertainty = StdDevUncertainty( - self.uncertainty.represent_as(StdDevUncertainty).quantity.to( - unit, equivalencies=equivalencies)) - + new_spec._convert_flux( + unit, equivalencies=equivalencies, suppress_conversion=suppress_conversion) return new_spec @property @@ -175,7 +179,7 @@ def velocity(self): Returns ------- - ~`astropy.units.Quantity` + new_data : `~astropy.units.Quantity` The converted dispersion array in the new dispersion space. """ if self.rest_value is None: @@ -202,8 +206,7 @@ def with_spectral_unit(self, unit, velocity_convention=None, self.with_spectral_axis_unit(unit, velocity_convention=velocity_convention, rest_value=rest_value) - def with_spectral_axis_unit(self, unit, velocity_convention=None, - rest_value=None): + def with_spectral_axis_unit(self, unit, velocity_convention=None, rest_value=None): """ Returns a new spectrum with a different spectral axis unit. Note that this creates a new object using the converted spectral axis and thus drops the original WCS, if it existed, @@ -230,11 +233,9 @@ def with_spectral_axis_unit(self, unit, velocity_convention=None, even if your spectrum has air wavelength units """ - - velocity_convention = velocity_convention if velocity_convention is not None else self.velocity_convention # noqa + velocity_convention = velocity_convention if velocity_convention is not None else self.velocity_convention # noqa rest_value = rest_value if rest_value is not None else self.rest_value - unit = self._new_wcs_argument_validation(unit, velocity_convention, - rest_value) + unit = self._new_wcs_argument_validation(unit, velocity_convention, rest_value) # Store the original unit information and WCS for posterity meta = deepcopy(self._meta) @@ -252,6 +253,24 @@ def with_spectral_axis_unit(self, unit, velocity_convention=None, return self.__class__(flux=self.flux, spectral_axis=new_spectral_axis, meta=meta, uncertainty=self.uncertainty, mask=self.mask) + def with_spectral_axis_and_flux_units(self, spectral_axis_unit, flux_unit, + velocity_convention=None, rest_value=None, + flux_equivalencies=None, suppress_flux_conversion=False): + """Perform :meth:`with_spectral_axis_unit` and :meth:`with_flux_unit` together. + See the respective methods for input and output definitions. + + Returns + ------- + new_spec : `~specutils.Spectrum1D` + Spectrum in requested units. + + """ + new_spec = self.with_spectral_axis_unit( + spectral_axis_unit, velocity_convention=velocity_convention, rest_value=rest_value) + new_spec._convert_flux( + flux_unit, equivalencies=flux_equivalencies, suppress_conversion=suppress_flux_conversion) + return new_spec + def _new_wcs_argument_validation(self, unit, velocity_convention, rest_value): # Allow string specification of units, for example diff --git a/specutils/tests/test_spectrum1d.py b/specutils/tests/test_spectrum1d.py index 7b8e1aac4..39ea23278 100644 --- a/specutils/tests/test_spectrum1d.py +++ b/specutils/tests/test_spectrum1d.py @@ -25,7 +25,7 @@ def test_empty_spectrum(): def test_create_from_arrays(): spec = Spectrum1D(spectral_axis=np.arange(50) * u.AA, - flux=np.random.randn(50) * u.Jy) + flux=np.ones(50) * u.Jy) assert isinstance(spec.spectral_axis, SpectralCoord) assert spec.spectral_axis.size == 50 @@ -36,7 +36,7 @@ def test_create_from_arrays(): # Test creating spectrum with unknown arguments with pytest.raises(ValueError): spec = Spectrum1D(wavelength=np.arange(1, 50) * u.nm, - flux=np.random.randn(48) * u.Jy) + flux=np.ones(48) * u.Jy) def test_create_from_multidimensional_arrays(): @@ -47,7 +47,7 @@ def test_create_from_multidimensional_arrays(): """ freqs = np.arange(50) * u.GHz - flux = np.random.random((5, len(freqs))) * u.Jy + flux = np.ones((5, len(freqs))) * u.Jy spec = Spectrum1D(spectral_axis=freqs, flux=flux) assert (spec.frequency == freqs).all() @@ -55,15 +55,15 @@ def test_create_from_multidimensional_arrays(): # Mis-matched lengths should raise an exception (unless freqs is one longer # than flux, in which case it's interpreted as bin edges) - freqs = np.arange(50) * u.GHz - flux = np.random.random((5, len(freqs)-10)) * u.Jy + flux = np.ones((5, len(freqs) - 10)) * u.Jy with pytest.raises(ValueError): spec = Spectrum1D(spectral_axis=freqs, flux=flux) def test_create_from_quantities(): - spec = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, - flux=np.random.randn(49) * u.Jy) + wav = np.arange(1, 50) * u.nm + flux = np.ones(49) * u.Jy + spec = Spectrum1D(spectral_axis=wav, flux=flux) assert isinstance(spec.spectral_axis, SpectralCoord) assert spec.spectral_axis.unit == u.nm @@ -72,13 +72,12 @@ def test_create_from_quantities(): # Mis-matched lengths should raise an exception (unless freqs is one longer # than flux, in which case it's interpreted as bin edges) with pytest.raises(ValueError): - spec = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, - flux=np.random.randn(47) * u.Jy) + spec = Spectrum1D(spectral_axis=wav, flux=np.ones(47) * u.Jy) def test_create_implicit_wcs(): spec = Spectrum1D(spectral_axis=np.arange(50) * u.AA, - flux=np.random.randn(50) * u.Jy) + flux=np.ones(50) * u.Jy) assert isinstance(spec.wcs, gwcs.wcs.WCS) @@ -90,7 +89,7 @@ def test_create_implicit_wcs(): def test_create_implicit_wcs_with_spectral_unit(): spec = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, - flux=np.random.randn(49) * u.Jy) + flux=np.ones(49) * u.Jy) assert isinstance(spec.wcs, gwcs.wcs.WCS) @@ -102,8 +101,8 @@ def test_create_implicit_wcs_with_spectral_unit(): def test_create_with_spectral_coord(): - spectral_coord = SpectralCoord(np.arange(5100, 5150)*u.AA, radial_velocity=u.Quantity(1000.0, "km/s")) - flux = np.random.randn(50)*u.Jy + spectral_coord = SpectralCoord(np.arange(5100, 5150) * u.AA, radial_velocity=u.Quantity(1000.0, "km/s")) + flux = np.ones(50) * u.Jy spec = Spectrum1D(spectral_axis=spectral_coord, flux=flux) assert spec.radial_velocity == u.Quantity(1000.0, "km/s") @@ -137,26 +136,24 @@ def test_spectral_axis_conversions(): assert np.all(spec.spectral_axis == np.array([400, 500]) * u.angstrom) assert spec.spectral_axis.unit == u.angstrom - spec = Spectrum1D(spectral_axis=np.arange(50) * u.AA, - flux=np.random.randn(50) * u.Jy) + flux = np.ones(49) * u.Jy + spec = Spectrum1D(spectral_axis=np.arange(flux.size) * u.AA, flux=flux) assert spec.wavelength.unit == u.AA - spec = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, - flux=np.random.randn(49) * u.Jy) + spec = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, flux=flux) assert spec.frequency.unit == u.GHz with pytest.raises(ValueError): spec.velocity - spec = Spectrum1D(spectral_axis=np.arange(100, 150) * u.nm, - flux=np.random.randn(49) * u.Jy) + spec = Spectrum1D(spectral_axis=np.arange(100, 150) * u.nm, flux=flux) - new_spec = spec.with_spectral_axis_unit(u.km/u.s, rest_value=125*u.um, - velocity_convention="relativistic") + new_spec = spec.with_spectral_axis_unit(u.km / u.s, rest_value=125 * u.um, + velocity_convention="relativistic") - assert new_spec.spectral_axis.unit == u.km/u.s + assert new_spec.spectral_axis.unit == u.km / u.s assert new_spec.wcs.world_axis_units[0] == "km.s**-1" # Make sure meta stored the old WCS correctly assert new_spec.meta["original_wcs"].world_axis_units[0] == "nm" @@ -164,10 +161,10 @@ def test_spectral_axis_conversions(): wcs_dict = {"CTYPE1": "WAVE", "CRVAL1": 3.622e3, "CDELT1": 8e-2, "CRPIX1": 0, "CUNIT1": "Angstrom"} - wcs_spec = Spectrum1D(flux=np.random.randn(49) * u.Jy, wcs=WCS(wcs_dict), + wcs_spec = Spectrum1D(flux=flux, wcs=WCS(wcs_dict), meta={'header': wcs_dict.copy()}) - new_spec = wcs_spec.with_spectral_axis_unit(u.km/u.s, rest_value=125*u.um, - velocity_convention="relativistic") + new_spec = wcs_spec.with_spectral_axis_unit(u.km / u.s, rest_value=125 * u.um, + velocity_convention="relativistic") new_spec.meta['original_wcs'].wcs.crval = [3.777e-7] new_spec.meta['header']['CRVAL1'] = 3777.0 @@ -175,9 +172,26 @@ def test_spectral_axis_conversions(): assert wcs_spec.meta['header']['CRVAL1'] == 3622. +def test_spectral_axis_and_flux_conversions(): + """A little bit from both sets of tests.""" + spec = Spectrum1D(spectral_axis=np.arange(100, 150) * u.nm, + flux=np.ones(49) * u.Jy) + + new_spec = spec.with_spectral_axis_and_flux_units( + u.km / u.s, u.uJy, rest_value=125 * u.um, velocity_convention="relativistic") + + assert new_spec.spectral_axis.unit == u.km/u.s + assert new_spec.wcs.world_axis_units[0] == "km.s**-1" + # Make sure meta stored the old WCS correctly + assert new_spec.meta["original_wcs"].world_axis_units[0] == "nm" + assert new_spec.meta["original_spectral_axis_unit"] == "nm" + assert new_spec.flux.unit == u.uJy + assert_allclose(new_spec.flux.value, 1000000) + + def test_spectral_slice(): spec = Spectrum1D(spectral_axis=np.linspace(100, 1000, 10) * u.nm, - flux=np.random.random(10) * u.Jy) + flux=np.ones(10) * u.Jy) sliced_spec = spec[300*u.nm:600*u.nm] assert np.all(sliced_spec.spectral_axis == [300, 400, 500] * u.nm) @@ -192,7 +206,7 @@ def test_spectral_slice(): # Test higher dimensional slicing spec = Spectrum1D(spectral_axis=np.linspace(100, 1000, 10) * u.nm, - flux=np.random.random((10, 10)) * u.Jy) + flux=np.ones((10, 10)) * u.Jy) sliced_spec = spec[300*u.nm:600*u.nm] assert np.all(sliced_spec.spectral_axis == [300, 400, 500] * u.nm) @@ -302,7 +316,7 @@ def test_flux_unit_conversion(): def test_wcs_transformations(): # Test with a GWCS spec = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, - flux=np.random.randn(49) * u.Jy) + flux=np.ones(49) * u.Jy) pix_axis = spec.wcs.world_to_pixel(np.arange(20, 30) * u.nm) disp_axis = spec.wcs.pixel_to_world(np.arange(20, 30)) @@ -359,24 +373,19 @@ def test_create_explicit_fitswcs(): def test_create_with_uncertainty(): spec = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, - flux=np.random.sample(49) * u.Jy, - uncertainty=StdDevUncertainty(np.random.sample(49) * 0.1)) + flux=np.ones(49) * u.Jy, + uncertainty=StdDevUncertainty(np.ones(49) * 0.1)) assert isinstance(spec.uncertainty, StdDevUncertainty) - - spec = Spectrum1D(spectral_axis=np.arange(1, 50) * u.nm, - flux=np.random.sample(49) * u.Jy, - uncertainty=StdDevUncertainty(np.random.sample(49) * 0.1)) - assert spec.flux.unit == spec.uncertainty.unit # If flux and uncertainty are different sizes then raise exception - wavelengths = np.arange(0, 10) - flux=100*np.abs(np.random.randn(3, 4, 10))*u.Jy - uncertainty = StdDevUncertainty(np.abs(np.random.randn(3, 2, 10))*u.Jy) + wavelengths = np.arange(10) * u.um + flux= np.ones((3, 4, 10)) * u.Jy + uncertainty = StdDevUncertainty(np.ones((3, 2, 10)) * u.Jy) with pytest.raises(ValueError): - Spectrum1D(spectral_axis=wavelengths*u.um, flux=flux, uncertainty=uncertainty) + Spectrum1D(spectral_axis=wavelengths, flux=flux, uncertainty=uncertainty) @pytest.mark.parametrize("flux_unit", ["adu", "ct/s", "count"]) @@ -407,7 +416,7 @@ def test_read_linear_solution(remote_data_path): def test_energy_photon_flux(): spec = Spectrum1D(spectral_axis=np.linspace(100, 1000, 10) * u.nm, - flux=np.random.randn(10)*u.Jy) + flux=np.ones(10) * u.Jy) assert spec.energy.size == 10 assert spec.photon_flux.size == 10 assert spec.photon_flux.unit == u.photon * u.cm**-2 * u.s**-1 * u.nm**-1 @@ -415,7 +424,7 @@ def test_energy_photon_flux(): def test_flux_nans_propagate_to_mask(): """Check that indices in input flux with NaNs get propagated to the mask""" - flux = np.random.randn(10) + flux = np.ones(10) nan_idx = [0, 3, 5] flux[nan_idx] = np.nan spec = Spectrum1D(spectral_axis=np.linspace(100, 1000, 10) * u.nm, @@ -424,16 +433,15 @@ def test_flux_nans_propagate_to_mask(): def test_repr(): - spec_with_wcs = Spectrum1D(spectral_axis=np.linspace(100, 1000, 10) * u.nm, - flux=np.random.random(10) * u.Jy) + wav = np.linspace(100, 1000, 10) * u.nm + flux = np.ones(10) * u.Jy + spec_with_wcs = Spectrum1D(spectral_axis=wav, flux=flux) result = repr(spec_with_wcs) assert result.startswith(' Date: Tue, 15 Oct 2024 15:47:08 -0400 Subject: [PATCH 05/16] Fix mixed unit spectral region extraction (#1187) * Fix extraction when only one of region/spectrum are in a reversed spectral axis unit * Changelog * Add test * Get rid of more pointless random calls * Update specutils/tests/test_region_extract.py Co-authored-by: P. L. Lim <2090236+pllim@users.noreply.github.com> * Better test improvement --------- Co-authored-by: P. L. Lim <2090236+pllim@users.noreply.github.com> --- CHANGES.rst | 3 ++ .../manipulation/extract_spectral_region.py | 6 ++- specutils/tests/test_region_extract.py | 37 +++++++++---------- 3 files changed, 26 insertions(+), 20 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index d2a0ebb60..f7fafc58c 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -13,6 +13,9 @@ Bug Fixes - Fixed ``Spectrum1D.with_flux_unit()`` not converting uncertainty along with flux unit. [#1181] +- Fixed extracting a spectral region when one of spectrum/region is in wavelength + and the other is in frequency units. [#1187] + Other Changes and Additions ^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/specutils/manipulation/extract_spectral_region.py b/specutils/manipulation/extract_spectral_region.py index 720754130..a9e5e2963 100644 --- a/specutils/manipulation/extract_spectral_region.py +++ b/specutils/manipulation/extract_spectral_region.py @@ -103,7 +103,11 @@ def _subregion_to_edge_pixels(subregion, spectrum): right_index = _edge_value_to_pixel(right_reg_in_spec_unit, spectrum, order, "right") - return left_index, right_index + # If the spectrum is in wavelength and region is in Hz (for example), these still might be reversed + if left_index < right_index: + return left_index, right_index + else: + return right_index, left_index def extract_region(spectrum, region, return_single_spectrum=False): diff --git a/specutils/tests/test_region_extract.py b/specutils/tests/test_region_extract.py index f40eb5b9f..2c46f6f12 100644 --- a/specutils/tests/test_region_extract.py +++ b/specutils/tests/test_region_extract.py @@ -18,10 +18,9 @@ def test_region_simple(simulated_spectra): - np.random.seed(42) spectrum = simulated_spectra.s1_um_mJy_e1 - uncertainty = StdDevUncertainty(0.1*np.random.random(len(spectrum.flux))*u.mJy) + uncertainty = StdDevUncertainty(np.full(spectrum.flux.shape, 0.1) * u.mJy) spectrum.uncertainty = uncertainty region = SpectralRegion(0.6*u.um, 0.8*u.um) @@ -112,10 +111,8 @@ def test_pixel_spectralaxis_extraction(): def test_slab_simple(simulated_spectra): - np.random.seed(42) - spectrum = simulated_spectra.s1_um_mJy_e1 - uncertainty = StdDevUncertainty(0.1*np.random.random(len(spectrum.flux))*u.mJy) + uncertainty = StdDevUncertainty(np.full(spectrum.flux.shape, 0.1) * u.mJy) spectrum.uncertainty = uncertainty sub_spectrum = spectral_slab(spectrum, 0.6*u.um, 0.8*u.um) @@ -148,6 +145,16 @@ def test_region_ghz(simulated_spectra): assert_quantity_allclose(sub_spectrum.flux, sub_spectrum_flux_expected) +def test_region_ghz_spectrum_wave(simulated_spectra): + spectrum = simulated_spectra.s1_um_mJy_e1 + region = SpectralRegion(499654.09666667*u.GHz, 374740.5725*u.GHz) + + sub_spectrum = extract_region(spectrum, region) + + sub_spectrum_flux_expected = FLUX_ARRAY * u.mJy + assert_quantity_allclose(sub_spectrum.flux, sub_spectrum_flux_expected) + + def test_region_simple_check_ends(simulated_spectra): np.random.seed(42) @@ -168,13 +175,11 @@ def test_region_simple_check_ends(simulated_spectra): assert sub_spectrum.spectral_axis.value[-1] == 25 -def test_region_empty(simulated_spectra): - np.random.seed(42) - +def test_region_empty(): empty_spectrum = Spectrum1D(spectral_axis=[]*u.um, flux=[]*u.Jy) # Region past upper range of spectrum - spectrum = Spectrum1D(spectral_axis=np.linspace(1, 25, 25)*u.um, flux=np.random.random(25)*u.Jy) + spectrum = Spectrum1D(spectral_axis=np.linspace(1, 25, 25)*u.um, flux=np.ones(25)*u.Jy) region = SpectralRegion(28*u.um, 30*u.um) sub_spectrum = extract_region(spectrum, region) @@ -185,7 +190,7 @@ def test_region_empty(simulated_spectra): assert sub_spectrum.flux.unit == empty_spectrum.flux.unit # Region below lower range of spectrum - spectrum = Spectrum1D(spectral_axis=np.linspace(1, 25, 25)*u.um, flux=np.random.random(25)*u.Jy) + spectrum = Spectrum1D(spectral_axis=np.linspace(1, 25, 25)*u.um, flux=np.ones(25)*u.Jy) region = SpectralRegion(0.1*u.um, 0.3*u.um) sub_spectrum = extract_region(spectrum, region) @@ -212,10 +217,8 @@ def test_region_empty(simulated_spectra): def test_region_descending(simulated_spectra): - np.random.seed(42) - spectrum = simulated_spectra.s1_um_mJy_e1 - uncertainty = StdDevUncertainty(0.1*np.random.random(len(spectrum.flux))*u.mJy) + uncertainty = StdDevUncertainty(np.full(spectrum.flux.shape, 0.1) * u.mJy) spectrum.uncertainty = uncertainty region = SpectralRegion(0.8*u.um, 0.6*u.um) @@ -244,10 +247,8 @@ def test_descending_spectral_axis(simulated_spectra): def test_region_two_sub(simulated_spectra): - np.random.seed(42) - spectrum = simulated_spectra.s1_um_mJy_e1 - uncertainty = StdDevUncertainty(0.1*np.random.random(len(spectrum.flux))*u.mJy) + uncertainty = StdDevUncertainty(np.full(spectrum.flux.shape, 0.1) * u.mJy) spectrum.uncertainty = uncertainty region = SpectralRegion([(0.6*u.um, 0.8*u.um), (0.86*u.um, 0.89*u.um)]) @@ -284,10 +285,8 @@ def test_region_two_sub(simulated_spectra): def test_bounding_region(simulated_spectra): - np.random.seed(42) - spectrum = simulated_spectra.s1_um_mJy_e1 - uncertainty = StdDevUncertainty(0.1*np.random.random(len(spectrum.flux))*u.mJy) + uncertainty = StdDevUncertainty(np.full(spectrum.flux.shape, 0.1) * u.mJy) spectrum.uncertainty = uncertainty region = SpectralRegion([(0.6*u.um, 0.8*u.um), (0.86*u.um, 0.89*u.um)]) From b70c393bcba1e7a2c310d940560336750e2444d3 Mon Sep 17 00:00:00 2001 From: riley <127359621+rileythai@users.noreply.github.com> Date: Wed, 16 Oct 2024 21:59:48 +1100 Subject: [PATCH 06/16] fix: HDU unspecified print message Co-authored-by: Brian Cherinka --- specutils/io/default_loaders/sdss_v.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/specutils/io/default_loaders/sdss_v.py b/specutils/io/default_loaders/sdss_v.py index e3195d31f..cf9a297b3 100644 --- a/specutils/io/default_loaders/sdss_v.py +++ b/specutils/io/default_loaders/sdss_v.py @@ -467,8 +467,7 @@ def load_sdss_mwm_1d(file_obj, hdu: Optional[int] = None, **kwargs): for i in range(len(hdulist)): if hdulist[i].header.get("DATASUM") != "0": hdu = i - print('HDU not specified. Loading spectrum at (HDU{})'. - format(i)) + print(f'HDU not specified. Loading spectrum at (HDU{i}).') break return _load_mwmVisit_or_mwmStar_hdu(hdulist, hdu, **kwargs) From 06d31b474ec48c7386c11fa2aae8b5a690a8be31 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Wed, 16 Oct 2024 12:19:02 -0400 Subject: [PATCH 07/16] Update changelog for 1.18.0 release --- CHANGES.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index f7fafc58c..02d5aa057 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,4 +1,4 @@ -1.18.0 (unreleased) +1.18.0 (2024-10-16) ------------------- New Features From e362fdd1b64e48183a0028f852ccdb16c1f33822 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Wed, 16 Oct 2024 15:18:18 -0400 Subject: [PATCH 08/16] Get min and max instead of specific indices --- specutils/manipulation/resample.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/specutils/manipulation/resample.py b/specutils/manipulation/resample.py index e208d2d14..b02596c30 100644 --- a/specutils/manipulation/resample.py +++ b/specutils/manipulation/resample.py @@ -461,8 +461,8 @@ def resample1d(self, orig_spectrum, fin_spec_axis): if self.extrapolation_treatment == 'zero_fill': fill_val = 0 - origedges = orig_spectrum.spectral_axis.bin_edges - off_edges = (fin_spec_axis < origedges[0]) | (origedges[-1] < fin_spec_axis) + orig_edges = orig_spectrum.spectral_axis.bin_edges + off_edges = (fin_spec_axis < np.min(orig_edges)) | (np.max(orig_edges) < fin_spec_axis) out_flux_val[off_edges] = fill_val if new_unc is not None: new_unc.array[off_edges] = fill_val From a34ebb8d5552c3cfbea5843de9d7a02cdd29e898 Mon Sep 17 00:00:00 2001 From: Ricky O'Steen Date: Wed, 16 Oct 2024 15:21:25 -0400 Subject: [PATCH 09/16] Changelog back to development --- CHANGES.rst | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/CHANGES.rst b/CHANGES.rst index 02d5aa057..78b724bed 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,3 +1,15 @@ +1.19.0 (unreleased) +------------------- + +New Features +^^^^^^^^^^^^ + +Bug Fixes +^^^^^^^^^ + +Other Changes and Additions +^^^^^^^^^^^^^^^^^^^^^^^^^^^ + 1.18.0 (2024-10-16) ------------------- From 62747c4bf89b71ed7389b82a25d9c10387f20b61 Mon Sep 17 00:00:00 2001 From: Riley Thai Date: Thu, 17 Oct 2024 09:17:44 +1100 Subject: [PATCH 10/16] feat: condense loaders into only SpectrumList of 1D flux - all loaders now only load for a single datatype, avoiding prior knowledge of SDSS datatypes - updated to only load as SpectrumList - updated to load all visits in mwmVisit files as individual Spectrum1D objects in the SpectrumList - relevant tests removed - relevant import __all__ adjusted --- specutils/io/default_loaders/sdss_v.py | 249 +++++------------- .../io/default_loaders/tests/test_sdss_v.py | 175 +----------- 2 files changed, 74 insertions(+), 350 deletions(-) diff --git a/specutils/io/default_loaders/sdss_v.py b/specutils/io/default_loaders/sdss_v.py index e3195d31f..804a9f375 100644 --- a/specutils/io/default_loaders/sdss_v.py +++ b/specutils/io/default_loaders/sdss_v.py @@ -1,23 +1,21 @@ """Register reader functions for various spectral formats.""" +import warnings from typing import Optional import numpy as np -from astropy.units import Unit, Quantity, Angstrom -from astropy.nddata import StdDevUncertainty, InverseVariance -from astropy.io.fits import HDUList, BinTableHDU, ImageHDU +from astropy.io.fits import BinTableHDU, HDUList, ImageHDU +from astropy.nddata import InverseVariance, StdDevUncertainty +from astropy.units import Angstrom, Quantity, Unit +from astropy.utils.exceptions import AstropyUserWarning from ...spectra import Spectrum1D, SpectrumList -from ..registers import data_loader from ..parsing_utils import read_fileobj_or_hdulist +from ..registers import data_loader __all__ = [ - "load_sdss_apVisit_1D", "load_sdss_apVisit_list", - "load_sdss_apStar_1D", "load_sdss_apStar_list", - "load_sdss_spec_1D", "load_sdss_spec_list", - "load_sdss_mwm_1d", "load_sdss_mwm_list", ] @@ -119,16 +117,9 @@ def _fetch_flux_unit(hdu): # APOGEE files -@data_loader( - "SDSS-V apStar", - identifier=apStar_identify, - dtype=Spectrum1D, - priority=10, - extensions=["fits"], -) -def load_sdss_apStar_1D(file_obj, idx: int = 0, **kwargs): +def _load_sdss_apStar_1D(file_obj, idx: int = 0, **kwargs): """ - Load an apStar file as a Spectrum1D. + Subloader to load an apStar file as a Spectrum1D. Parameters ---------- @@ -206,59 +197,15 @@ def load_sdss_apStar_list(file_obj, **kwargs): """ with read_fileobj_or_hdulist(file_obj, **kwargs) as hdulist: nvisits = hdulist[0].header.get("NVISITS") - if nvisits <= 1: - raise ValueError( - "Only 1 visit in this file. Use Spectrum1D.read() instead.") + if nvisits < 1: + raise RuntimeError( + "No visits in this file. Use Spectrum1D.read() instead.") return SpectrumList([ - load_sdss_apStar_1D(file_obj, idx=i, **kwargs) + _load_sdss_apStar_1D(file_obj, idx=i, **kwargs) for i in range(nvisits) ]) -@data_loader( - "SDSS-V apVisit", - identifier=apVisit_identify, - dtype=Spectrum1D, - priority=10, - extensions=["fits"], -) -def load_sdss_apVisit_1D(file_obj, **kwargs): - """ - Load an apVisit file as a Spectrum1D. - - Parameters - ---------- - file_obj : str, file-like, or HDUList - FITS file name, file object, or HDUList. - - Returns - ------- - Spectrum1D - The spectrum contained in the file - """ - flux_unit = Unit("1e-17 erg / (Angstrom cm2 s)") - - with read_fileobj_or_hdulist(file_obj, memmap=False, **kwargs) as hdulist: - # Fetch and flatten spectrum parameters - spectral_axis = Quantity(hdulist[4].data.flatten(), unit=Angstrom) - flux_unit = _fetch_flux_unit(hdulist[1]) - flux = Quantity(hdulist[1].data.flatten(), unit=flux_unit) - e_flux = StdDevUncertainty(hdulist[2].data.flatten()) - - # Get metadata and attach bitmask and MJD - meta = dict() - meta["header"] = hdulist[0].header - mask = hdulist[3].data.flatten() - # NOTE: specutils considers 0/False as valid values, simlar to numpy convention - mask = mask != 0 - - return Spectrum1D(spectral_axis=spectral_axis, - flux=flux, - mask=mask, - uncertainty=e_flux, - meta=meta) - - @data_loader( "SDSS-V apVisit", identifier=apVisit_identify, @@ -313,49 +260,12 @@ def load_sdss_apVisit_list(file_obj, **kwargs): # BOSS REDUX products (specLite, specFull, custom coadd files, etc) - - -@data_loader( - "SDSS-V spec", - identifier=spec_sdss5_identify, - dtype=Spectrum1D, - priority=5, - extensions=["fits"], -) -def load_sdss_spec_1D(file_obj, *args, hdu: Optional[int] = None, **kwargs): - """ - Load a given BOSS spec file as a Spectrum1D object. - - Parameters - ---------- - file_obj : str, file-like, or HDUList - FITS file name, file object, or HDUList.. - hdu : int - The specified HDU to load a given spectra from. - - Returns - ------- - Spectrum1D - The spectrum contained in the file at the specified HDU. - """ - if hdu is None: - # TODO: how should we handle this -- multiple things in file, but the user cannot choose. - print('HDU not specified. Loading coadd spectrum (HDU1)') - hdu = 1 # defaulting to coadd - # raise ValueError("HDU not specified! Please specify a HDU to load.") - elif hdu in [2, 3, 4]: - raise ValueError("Invalid HDU! HDU{} is not spectra.".format(hdu)) - with read_fileobj_or_hdulist(file_obj, memmap=False, **kwargs) as hdulist: - # directly load the coadd at HDU1 - return _load_BOSS_HDU(hdulist, hdu, **kwargs) - - @data_loader( "SDSS-V spec", identifier=spec_sdss5_identify, dtype=SpectrumList, force=True, - priority=5, + priority=10, extensions=["fits"], ) def load_sdss_spec_list(file_obj, **kwargs): @@ -431,55 +341,12 @@ def _load_BOSS_HDU(hdulist: HDUList, hdu: int, **kwargs): # MWM LOADERS -@data_loader( - "SDSS-V mwm", - identifier=mwm_identify, - dtype=Spectrum1D, - priority=20, - extensions=["fits"], -) -def load_sdss_mwm_1d(file_obj, hdu: Optional[int] = None, **kwargs): - """ - Load an unspecified spec file as a Spectrum1D. - - Parameters - ---------- - file_obj : str, file-like, or HDUList - FITS file name, file object, or HDUList.. - hdu : int - Specified HDU to load. - - Returns - ------- - Spectrum1D - The spectrum contained in the file - """ - with read_fileobj_or_hdulist(file_obj, memmap=False, **kwargs) as hdulist: - # Check if file is empty first - datasums = [] - for i in range(1, len(hdulist)): - datasums.append(int(hdulist[i].header.get("DATASUM"))) - if (np.array(datasums) == 0).all(): - raise ValueError("Specified file is empty.") - - # TODO: how should we handle this -- multiple things in file, but the user cannot choose. - if hdu is None: - for i in range(len(hdulist)): - if hdulist[i].header.get("DATASUM") != "0": - hdu = i - print('HDU not specified. Loading spectrum at (HDU{})'. - format(i)) - break - - return _load_mwmVisit_or_mwmStar_hdu(hdulist, hdu, **kwargs) - - @data_loader( "SDSS-V mwm", identifier=mwm_identify, force=True, dtype=SpectrumList, - priority=20, + priority=10, extensions=["fits"], ) def load_sdss_mwm_list(file_obj, **kwargs): @@ -494,11 +361,7 @@ def load_sdss_mwm_list(file_obj, **kwargs): Returns ------- SpectrumList - The spectra contained in the file, where: - Spectrum1D - A given spectra of nD flux - None - If there are no spectra for that spectrograph/observatory + The spectra contained in the file, each of 1D flux. """ spectra = SpectrumList() with read_fileobj_or_hdulist(file_obj, memmap=False, **kwargs) as hdulist: @@ -507,18 +370,18 @@ def load_sdss_mwm_list(file_obj, **kwargs): for hdu in range(1, len(hdulist)): datasums.append(int(hdulist[hdu].header.get("DATASUM"))) if (np.array(datasums) == 0).all(): - raise ValueError("Specified file is empty.") + raise RuntimeError("Specified file is empty.") # Now load file for hdu in range(1, len(hdulist)): if hdulist[hdu].header.get("DATASUM") == "0": # Skip zero data HDU's - # TODO: validate if we want this printed warning or not. - # it might get annoying & fill logs with useless alerts. - print("WARNING: HDU{} ({}) is empty.".format( - hdu, hdulist[hdu].name)) + warnings.warn( + "WARNING: HDU{} ({}) is empty.".format( + hdu, hdulist[hdu].name), AstropyUserWarning) continue - spectra.append(_load_mwmVisit_or_mwmStar_hdu(hdulist, hdu)) + spectra_list = _load_mwmVisit_or_mwmStar_hdu(hdulist, hdu) + [spectra.append(spectra_list[i]) for i in range(len(spectra_list))] return spectra @@ -560,7 +423,7 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): spectral_axis = Quantity(wavelength, unit=Angstrom) # Fetch flux, e_flux - # NOTE:: flux info is not written into mwm files + # NOTE: flux info is not written into mwm files flux_unit = Unit("1e-17 erg / (Angstrom cm2 s)") # NOTE: hardcoded unit flux = Quantity(hdulist[hdu].data["flux"], unit=flux_unit) e_flux = InverseVariance(array=hdulist[hdu].data["ivar"]) @@ -570,14 +433,6 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): # NOTE: specutils considers 0/False as valid values, simlar to numpy convention mask = mask != 0 - # collapse shape if 1D spectra in 2D array - # NOTE: this fixes a jdaviz handling bug for 2D of shape 1, - # it could be that it's expected to be parsed this way. - if flux.shape[0] == 1: - flux = flux[0] - e_flux = e_flux[0] - mask = mask[0] - # Create metadata meta = dict() meta["header"] = hdulist[0].header @@ -598,17 +453,59 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): meta['mjd'] = hdulist[hdu].data['mjd'] meta["datatype"] = "mwmVisit" except KeyError: - meta["mjd"] = (str(hdulist[hdu].data["min_mjd"][0]) + "->" + - str(hdulist[hdu].data["max_mjd"][0])) + meta["min_mjd"] = str(hdulist[hdu].data["min_mjd"][0]) + meta["max_mjd"] = str(hdulist[hdu].data["max_mjd"][0]) meta["datatype"] = "mwmStar" finally: meta["name"] = hdulist[hdu].name meta["sdss_id"] = hdulist[hdu].data['sdss_id'] - return Spectrum1D( - spectral_axis=spectral_axis, - flux=flux, - uncertainty=e_flux, - mask=mask, - meta=meta, - ) + # drop back a list of Spectrum1Ds to unpack + metadicts = _split_mwm_meta_dict(meta) + return [ + Spectrum1D( + spectral_axis=spectral_axis, + flux=flux[i], + uncertainty=e_flux[i], + mask=mask[i], + meta=metadicts[i], + ) for i in range(flux.shape[0]) + ] + + +def _split_mwm_meta_dict(d): + """ + Metadata sub-loader subfunction for MWM files. + + Parameters + ---------- + d: dict + Initial metadata dictionary. + + Returns + ------- + dicts: list[dict] + List of dicts with unpacked metadata for length > 1 array objects. + + """ + # find the length of entries + N = max(len(v) if isinstance(v, np.ndarray) else 1 for v in d.values()) + + # create N dictionaries to hold the split results + dicts = [{} for _ in range(N)] + + for key, value in d.items(): + if isinstance(value, np.ndarray): + # Ensure that the length of the list matches N + if len(value) != N: + # an error case we ignore + continue + # distribute each element to the corresponding metadict + for i in range(N): + dicts[i][key] = value[i] + else: + # if it's a single object, copy it to each metadict + for i in range(N): + dicts[i][key] = value + + return dicts diff --git a/specutils/io/default_loaders/tests/test_sdss_v.py b/specutils/io/default_loaders/tests/test_sdss_v.py index 15f9916cb..9993efc2d 100644 --- a/specutils/io/default_loaders/tests/test_sdss_v.py +++ b/specutils/io/default_loaders/tests/test_sdss_v.py @@ -425,40 +425,6 @@ def spec_HDUList(n_spectra): # TEST MWM loaders -@pytest.mark.parametrize( - "file_obj, hdu, with_wl, hduflags", - [ - ("mwm-temp", None, False, [0, 0, 1, 0]), - ("mwm-temp", 3, False, [0, 0, 1, 0]), - ("mwm-temp", None, True, [0, 1, 1, 0]), - ("mwm-temp", 2, True, [0, 1, 1, 0]), - ], -) -def test_mwm_1d(file_obj, hdu, with_wl, hduflags): - """Test mwm Spectrum1D loader""" - tmpfile = str(file_obj) + ".fits" - mwm_HDUList(hduflags, with_wl).writeto(tmpfile, overwrite=True) - - if hdu is None: - data = Spectrum1D.read(tmpfile) - else: - data = Spectrum1D.read(tmpfile, hdu=hdu) - assert isinstance(data, Spectrum1D) - assert isinstance(data.meta["header"], fits.Header) - if data.meta["instrument"].lower() == "apogee": - length = 8575 - elif data.meta["instrument"].lower() == "boss": - length = 4648 - else: - raise ValueError( - "INSTRMNT tag in test HDU header is not set properly.") - assert len(data.spectral_axis.value) == length - assert len(data.flux.value) == length - assert data.spectral_axis.unit == Angstrom - assert data.flux.unit == Unit("1e-17 erg / (s cm2 Angstrom)") - os.remove(tmpfile) - - @pytest.mark.parametrize( "file_obj, with_wl, hduflags", [ @@ -500,41 +466,6 @@ def test_mwm_list(file_obj, with_wl, hduflags): os.remove(tmpfile) -@pytest.mark.parametrize( - "file_obj, hdu, hduflags", - [ - ("mwm-temp", 1, [0, 1, 0, 1]), - ("mwm-temp", 2, [1, 0, 1, 1]), - ("mwm-temp", 3, [0, 0, 0, 1]), - ("mwm-temp", 4, [0, 1, 1, 0]), - ], -) -def test_mwm_1d_fail_spec(file_obj, hdu, hduflags): - """Test mwm Spectrum1D loader fail on bad spec""" - tmpfile = str(file_obj) + ".fits" - mwm_HDUList(hduflags, True).writeto(tmpfile, overwrite=True) - with pytest.raises(IndexError): - Spectrum1D.read(tmpfile, hdu=hdu) - os.remove(tmpfile) - - -@pytest.mark.parametrize( - "file_obj, with_wl", - [ - ("mwm-temp", False), - ("mwm-temp", True), - ], -) -def test_mwm_1d_fail(file_obj, with_wl): - """Test mwm Spectrum1D loader fail on empty""" - tmpfile = str(file_obj) + ".fits" - mwm_HDUList([0, 0, 0, 0], with_wl).writeto(tmpfile, overwrite=True) - - with pytest.raises(ValueError): - Spectrum1D.read(tmpfile) - os.remove(tmpfile) - - @pytest.mark.parametrize( "file_obj, with_wl", [ @@ -547,39 +478,11 @@ def test_mwm_list_fail(file_obj, with_wl): tmpfile = str(file_obj) + ".fits" mwm_HDUList([0, 0, 0, 0], with_wl).writeto(tmpfile, overwrite=True) - with pytest.raises(ValueError): + with pytest.raises(RuntimeError): SpectrumList.read(tmpfile) os.remove(tmpfile) -@pytest.mark.parametrize( - "file_obj,n_spectra", - [ - ("spec-temp", 1), - ("spec-temp", 5), - ], -) -def test_spec_1d(file_obj, n_spectra): - """Test BOSS spec loader""" - tmpfile = str(file_obj) + ".fits" - spec_HDUList(n_spectra).writeto(tmpfile, overwrite=True) - - idxs = [1] + list(np.arange(5, 5 + n_spectra, 1)) - - for i in idxs: - data = Spectrum1D.read(tmpfile, hdu=i) - assert isinstance(data, Spectrum1D) - assert len(data.flux.value) == 10 - assert data.flux.unit == Unit("1e-17 erg / (s cm2 Angstrom)") - - assert len(data.spectral_axis.value) == 10 - assert data.spectral_axis.unit == Angstrom - assert len(data.mask) == 10 - - assert data[i].meta["header"].get("foobar") == "barfoo" - os.remove(tmpfile) - - @pytest.mark.parametrize( "file_obj,n_spectra", [ @@ -605,48 +508,6 @@ def test_spec_list(file_obj, n_spectra): os.remove(tmpfile) -@pytest.mark.parametrize( - "file_obj,hdu", - [ - ("spec-temp", 2), - ("spec-temp", 3), - ("spec-temp", 4), - ], - ids=["Fail on HDU2", "Fail on HDU3", "Fail on HDU4"], -) -def test_spec_1d_fail_hdu(file_obj, hdu): - """Test if fail on reading HDU2, HDU3, or HDU4 (non-spectra)""" - tmpfile = str(file_obj) + ".fits" - spec_HDUList(5).writeto(tmpfile, overwrite=True) - - with pytest.raises(ValueError): - Spectrum1D.read(tmpfile, hdu=hdu) - os.remove(tmpfile) - - -@pytest.mark.parametrize( - "file_obj,idx", - [ - ("apStar-temp", 1), - ("apStar-temp", 4), - ], -) -def test_apStar_1D(file_obj, idx): - tmpfile = str(file_obj) + ".fits" - apStar_HDUList(6).writeto(tmpfile, overwrite=True) - - data = Spectrum1D.read(tmpfile, idx=idx) - assert isinstance(data, Spectrum1D) - assert len(data.flux.value) == 10 - assert data.flux.unit == Unit("1e-17 erg / (s cm2 Angstrom)") - - assert len(data.spectral_axis.value) == 10 - assert data.spectral_axis.unit == Angstrom - - assert data.meta["header"].get("foobar") == "barfoo" - os.remove(tmpfile) - - @pytest.mark.parametrize( "file_obj,n_spectra", [ @@ -670,40 +531,6 @@ def test_apStar_list(file_obj, n_spectra): os.remove(tmpfile) -@pytest.mark.parametrize( - "file_obj", - [ - ("apStar-temp"), - ], -) -def test_apStar_fail_list(file_obj): - """Test if fail on reading 1D apStar as list""" - tmpfile = str(file_obj) + ".fits" - apStar_HDUList(1).writeto(tmpfile, overwrite=True) - - with pytest.raises(ValueError): - SpectrumList.read(tmpfile) - os.remove(tmpfile) - - -@pytest.mark.parametrize( - "file_obj", - [ - ("apVisit-temp"), - ], -) -def test_apVisit_1D(file_obj): - tmpfile = str(file_obj) + ".fits" - apVisit_HDUList().writeto(tmpfile, overwrite=True) - - data = Spectrum1D.read(tmpfile) - assert isinstance(data, Spectrum1D) - assert np.array_equal(data.spectral_axis.value, np.arange(1, 31, 1)) - assert len(data.flux.value) == 30 - assert data.meta["header"].get("foobar") == "barfoo" - os.remove(tmpfile) - - @pytest.mark.parametrize( "file_obj", [ From bafede99efaa463fec4adca23d356ee058ccbb78 Mon Sep 17 00:00:00 2001 From: Riley Thai Date: Thu, 17 Oct 2024 09:21:10 +1100 Subject: [PATCH 11/16] chore: changelog update - describes changes shown previously --- CHANGES.rst | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index af3ca4290..13ed32551 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -20,7 +20,9 @@ Other Changes and Additions - Replaced ``LevMarLSQFitter`` with ``TRFLSQFitter`` as the former is no longer recommended by ``astropy``. [#1180] -- "Multi" loaders have been removed from SDSS-V ``SpectrumList`` default loaders. [#1185] +- All "Multi" loaders have been merged into the standard SDSS-V ``SpectrumList`` default loaders. [#1185] + +- All ``Spectrum1D`` loaders have been removed from SDSS-V default loaders. [#1185] 1.17.0 (2024-10-04) ------------------- From 122b3684f44ec9c44a1f0fc1a6c3296865a961cf Mon Sep 17 00:00:00 2001 From: Riley Thai Date: Fri, 18 Oct 2024 17:17:46 +1100 Subject: [PATCH 12/16] revert: readd all Spectrum1D loaders and tests revert back to 4bee1369 --- CHANGES.rst | 4 +- specutils/io/default_loaders/sdss_v.py | 249 +++++++++++++----- .../io/default_loaders/tests/test_sdss_v.py | 175 +++++++++++- 3 files changed, 351 insertions(+), 77 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 13ed32551..af3ca4290 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -20,9 +20,7 @@ Other Changes and Additions - Replaced ``LevMarLSQFitter`` with ``TRFLSQFitter`` as the former is no longer recommended by ``astropy``. [#1180] -- All "Multi" loaders have been merged into the standard SDSS-V ``SpectrumList`` default loaders. [#1185] - -- All ``Spectrum1D`` loaders have been removed from SDSS-V default loaders. [#1185] +- "Multi" loaders have been removed from SDSS-V ``SpectrumList`` default loaders. [#1185] 1.17.0 (2024-10-04) ------------------- diff --git a/specutils/io/default_loaders/sdss_v.py b/specutils/io/default_loaders/sdss_v.py index 804a9f375..e3195d31f 100644 --- a/specutils/io/default_loaders/sdss_v.py +++ b/specutils/io/default_loaders/sdss_v.py @@ -1,21 +1,23 @@ """Register reader functions for various spectral formats.""" -import warnings from typing import Optional import numpy as np -from astropy.io.fits import BinTableHDU, HDUList, ImageHDU -from astropy.nddata import InverseVariance, StdDevUncertainty -from astropy.units import Angstrom, Quantity, Unit -from astropy.utils.exceptions import AstropyUserWarning +from astropy.units import Unit, Quantity, Angstrom +from astropy.nddata import StdDevUncertainty, InverseVariance +from astropy.io.fits import HDUList, BinTableHDU, ImageHDU from ...spectra import Spectrum1D, SpectrumList -from ..parsing_utils import read_fileobj_or_hdulist from ..registers import data_loader +from ..parsing_utils import read_fileobj_or_hdulist __all__ = [ + "load_sdss_apVisit_1D", "load_sdss_apVisit_list", + "load_sdss_apStar_1D", "load_sdss_apStar_list", + "load_sdss_spec_1D", "load_sdss_spec_list", + "load_sdss_mwm_1d", "load_sdss_mwm_list", ] @@ -117,9 +119,16 @@ def _fetch_flux_unit(hdu): # APOGEE files -def _load_sdss_apStar_1D(file_obj, idx: int = 0, **kwargs): +@data_loader( + "SDSS-V apStar", + identifier=apStar_identify, + dtype=Spectrum1D, + priority=10, + extensions=["fits"], +) +def load_sdss_apStar_1D(file_obj, idx: int = 0, **kwargs): """ - Subloader to load an apStar file as a Spectrum1D. + Load an apStar file as a Spectrum1D. Parameters ---------- @@ -197,15 +206,59 @@ def load_sdss_apStar_list(file_obj, **kwargs): """ with read_fileobj_or_hdulist(file_obj, **kwargs) as hdulist: nvisits = hdulist[0].header.get("NVISITS") - if nvisits < 1: - raise RuntimeError( - "No visits in this file. Use Spectrum1D.read() instead.") + if nvisits <= 1: + raise ValueError( + "Only 1 visit in this file. Use Spectrum1D.read() instead.") return SpectrumList([ - _load_sdss_apStar_1D(file_obj, idx=i, **kwargs) + load_sdss_apStar_1D(file_obj, idx=i, **kwargs) for i in range(nvisits) ]) +@data_loader( + "SDSS-V apVisit", + identifier=apVisit_identify, + dtype=Spectrum1D, + priority=10, + extensions=["fits"], +) +def load_sdss_apVisit_1D(file_obj, **kwargs): + """ + Load an apVisit file as a Spectrum1D. + + Parameters + ---------- + file_obj : str, file-like, or HDUList + FITS file name, file object, or HDUList. + + Returns + ------- + Spectrum1D + The spectrum contained in the file + """ + flux_unit = Unit("1e-17 erg / (Angstrom cm2 s)") + + with read_fileobj_or_hdulist(file_obj, memmap=False, **kwargs) as hdulist: + # Fetch and flatten spectrum parameters + spectral_axis = Quantity(hdulist[4].data.flatten(), unit=Angstrom) + flux_unit = _fetch_flux_unit(hdulist[1]) + flux = Quantity(hdulist[1].data.flatten(), unit=flux_unit) + e_flux = StdDevUncertainty(hdulist[2].data.flatten()) + + # Get metadata and attach bitmask and MJD + meta = dict() + meta["header"] = hdulist[0].header + mask = hdulist[3].data.flatten() + # NOTE: specutils considers 0/False as valid values, simlar to numpy convention + mask = mask != 0 + + return Spectrum1D(spectral_axis=spectral_axis, + flux=flux, + mask=mask, + uncertainty=e_flux, + meta=meta) + + @data_loader( "SDSS-V apVisit", identifier=apVisit_identify, @@ -260,12 +313,49 @@ def load_sdss_apVisit_list(file_obj, **kwargs): # BOSS REDUX products (specLite, specFull, custom coadd files, etc) + + +@data_loader( + "SDSS-V spec", + identifier=spec_sdss5_identify, + dtype=Spectrum1D, + priority=5, + extensions=["fits"], +) +def load_sdss_spec_1D(file_obj, *args, hdu: Optional[int] = None, **kwargs): + """ + Load a given BOSS spec file as a Spectrum1D object. + + Parameters + ---------- + file_obj : str, file-like, or HDUList + FITS file name, file object, or HDUList.. + hdu : int + The specified HDU to load a given spectra from. + + Returns + ------- + Spectrum1D + The spectrum contained in the file at the specified HDU. + """ + if hdu is None: + # TODO: how should we handle this -- multiple things in file, but the user cannot choose. + print('HDU not specified. Loading coadd spectrum (HDU1)') + hdu = 1 # defaulting to coadd + # raise ValueError("HDU not specified! Please specify a HDU to load.") + elif hdu in [2, 3, 4]: + raise ValueError("Invalid HDU! HDU{} is not spectra.".format(hdu)) + with read_fileobj_or_hdulist(file_obj, memmap=False, **kwargs) as hdulist: + # directly load the coadd at HDU1 + return _load_BOSS_HDU(hdulist, hdu, **kwargs) + + @data_loader( "SDSS-V spec", identifier=spec_sdss5_identify, dtype=SpectrumList, force=True, - priority=10, + priority=5, extensions=["fits"], ) def load_sdss_spec_list(file_obj, **kwargs): @@ -341,12 +431,55 @@ def _load_BOSS_HDU(hdulist: HDUList, hdu: int, **kwargs): # MWM LOADERS +@data_loader( + "SDSS-V mwm", + identifier=mwm_identify, + dtype=Spectrum1D, + priority=20, + extensions=["fits"], +) +def load_sdss_mwm_1d(file_obj, hdu: Optional[int] = None, **kwargs): + """ + Load an unspecified spec file as a Spectrum1D. + + Parameters + ---------- + file_obj : str, file-like, or HDUList + FITS file name, file object, or HDUList.. + hdu : int + Specified HDU to load. + + Returns + ------- + Spectrum1D + The spectrum contained in the file + """ + with read_fileobj_or_hdulist(file_obj, memmap=False, **kwargs) as hdulist: + # Check if file is empty first + datasums = [] + for i in range(1, len(hdulist)): + datasums.append(int(hdulist[i].header.get("DATASUM"))) + if (np.array(datasums) == 0).all(): + raise ValueError("Specified file is empty.") + + # TODO: how should we handle this -- multiple things in file, but the user cannot choose. + if hdu is None: + for i in range(len(hdulist)): + if hdulist[i].header.get("DATASUM") != "0": + hdu = i + print('HDU not specified. Loading spectrum at (HDU{})'. + format(i)) + break + + return _load_mwmVisit_or_mwmStar_hdu(hdulist, hdu, **kwargs) + + @data_loader( "SDSS-V mwm", identifier=mwm_identify, force=True, dtype=SpectrumList, - priority=10, + priority=20, extensions=["fits"], ) def load_sdss_mwm_list(file_obj, **kwargs): @@ -361,7 +494,11 @@ def load_sdss_mwm_list(file_obj, **kwargs): Returns ------- SpectrumList - The spectra contained in the file, each of 1D flux. + The spectra contained in the file, where: + Spectrum1D + A given spectra of nD flux + None + If there are no spectra for that spectrograph/observatory """ spectra = SpectrumList() with read_fileobj_or_hdulist(file_obj, memmap=False, **kwargs) as hdulist: @@ -370,18 +507,18 @@ def load_sdss_mwm_list(file_obj, **kwargs): for hdu in range(1, len(hdulist)): datasums.append(int(hdulist[hdu].header.get("DATASUM"))) if (np.array(datasums) == 0).all(): - raise RuntimeError("Specified file is empty.") + raise ValueError("Specified file is empty.") # Now load file for hdu in range(1, len(hdulist)): if hdulist[hdu].header.get("DATASUM") == "0": # Skip zero data HDU's - warnings.warn( - "WARNING: HDU{} ({}) is empty.".format( - hdu, hdulist[hdu].name), AstropyUserWarning) + # TODO: validate if we want this printed warning or not. + # it might get annoying & fill logs with useless alerts. + print("WARNING: HDU{} ({}) is empty.".format( + hdu, hdulist[hdu].name)) continue - spectra_list = _load_mwmVisit_or_mwmStar_hdu(hdulist, hdu) - [spectra.append(spectra_list[i]) for i in range(len(spectra_list))] + spectra.append(_load_mwmVisit_or_mwmStar_hdu(hdulist, hdu)) return spectra @@ -423,7 +560,7 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): spectral_axis = Quantity(wavelength, unit=Angstrom) # Fetch flux, e_flux - # NOTE: flux info is not written into mwm files + # NOTE:: flux info is not written into mwm files flux_unit = Unit("1e-17 erg / (Angstrom cm2 s)") # NOTE: hardcoded unit flux = Quantity(hdulist[hdu].data["flux"], unit=flux_unit) e_flux = InverseVariance(array=hdulist[hdu].data["ivar"]) @@ -433,6 +570,14 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): # NOTE: specutils considers 0/False as valid values, simlar to numpy convention mask = mask != 0 + # collapse shape if 1D spectra in 2D array + # NOTE: this fixes a jdaviz handling bug for 2D of shape 1, + # it could be that it's expected to be parsed this way. + if flux.shape[0] == 1: + flux = flux[0] + e_flux = e_flux[0] + mask = mask[0] + # Create metadata meta = dict() meta["header"] = hdulist[0].header @@ -453,59 +598,17 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): meta['mjd'] = hdulist[hdu].data['mjd'] meta["datatype"] = "mwmVisit" except KeyError: - meta["min_mjd"] = str(hdulist[hdu].data["min_mjd"][0]) - meta["max_mjd"] = str(hdulist[hdu].data["max_mjd"][0]) + meta["mjd"] = (str(hdulist[hdu].data["min_mjd"][0]) + "->" + + str(hdulist[hdu].data["max_mjd"][0])) meta["datatype"] = "mwmStar" finally: meta["name"] = hdulist[hdu].name meta["sdss_id"] = hdulist[hdu].data['sdss_id'] - # drop back a list of Spectrum1Ds to unpack - metadicts = _split_mwm_meta_dict(meta) - return [ - Spectrum1D( - spectral_axis=spectral_axis, - flux=flux[i], - uncertainty=e_flux[i], - mask=mask[i], - meta=metadicts[i], - ) for i in range(flux.shape[0]) - ] - - -def _split_mwm_meta_dict(d): - """ - Metadata sub-loader subfunction for MWM files. - - Parameters - ---------- - d: dict - Initial metadata dictionary. - - Returns - ------- - dicts: list[dict] - List of dicts with unpacked metadata for length > 1 array objects. - - """ - # find the length of entries - N = max(len(v) if isinstance(v, np.ndarray) else 1 for v in d.values()) - - # create N dictionaries to hold the split results - dicts = [{} for _ in range(N)] - - for key, value in d.items(): - if isinstance(value, np.ndarray): - # Ensure that the length of the list matches N - if len(value) != N: - # an error case we ignore - continue - # distribute each element to the corresponding metadict - for i in range(N): - dicts[i][key] = value[i] - else: - # if it's a single object, copy it to each metadict - for i in range(N): - dicts[i][key] = value - - return dicts + return Spectrum1D( + spectral_axis=spectral_axis, + flux=flux, + uncertainty=e_flux, + mask=mask, + meta=meta, + ) diff --git a/specutils/io/default_loaders/tests/test_sdss_v.py b/specutils/io/default_loaders/tests/test_sdss_v.py index 9993efc2d..15f9916cb 100644 --- a/specutils/io/default_loaders/tests/test_sdss_v.py +++ b/specutils/io/default_loaders/tests/test_sdss_v.py @@ -425,6 +425,40 @@ def spec_HDUList(n_spectra): # TEST MWM loaders +@pytest.mark.parametrize( + "file_obj, hdu, with_wl, hduflags", + [ + ("mwm-temp", None, False, [0, 0, 1, 0]), + ("mwm-temp", 3, False, [0, 0, 1, 0]), + ("mwm-temp", None, True, [0, 1, 1, 0]), + ("mwm-temp", 2, True, [0, 1, 1, 0]), + ], +) +def test_mwm_1d(file_obj, hdu, with_wl, hduflags): + """Test mwm Spectrum1D loader""" + tmpfile = str(file_obj) + ".fits" + mwm_HDUList(hduflags, with_wl).writeto(tmpfile, overwrite=True) + + if hdu is None: + data = Spectrum1D.read(tmpfile) + else: + data = Spectrum1D.read(tmpfile, hdu=hdu) + assert isinstance(data, Spectrum1D) + assert isinstance(data.meta["header"], fits.Header) + if data.meta["instrument"].lower() == "apogee": + length = 8575 + elif data.meta["instrument"].lower() == "boss": + length = 4648 + else: + raise ValueError( + "INSTRMNT tag in test HDU header is not set properly.") + assert len(data.spectral_axis.value) == length + assert len(data.flux.value) == length + assert data.spectral_axis.unit == Angstrom + assert data.flux.unit == Unit("1e-17 erg / (s cm2 Angstrom)") + os.remove(tmpfile) + + @pytest.mark.parametrize( "file_obj, with_wl, hduflags", [ @@ -466,6 +500,41 @@ def test_mwm_list(file_obj, with_wl, hduflags): os.remove(tmpfile) +@pytest.mark.parametrize( + "file_obj, hdu, hduflags", + [ + ("mwm-temp", 1, [0, 1, 0, 1]), + ("mwm-temp", 2, [1, 0, 1, 1]), + ("mwm-temp", 3, [0, 0, 0, 1]), + ("mwm-temp", 4, [0, 1, 1, 0]), + ], +) +def test_mwm_1d_fail_spec(file_obj, hdu, hduflags): + """Test mwm Spectrum1D loader fail on bad spec""" + tmpfile = str(file_obj) + ".fits" + mwm_HDUList(hduflags, True).writeto(tmpfile, overwrite=True) + with pytest.raises(IndexError): + Spectrum1D.read(tmpfile, hdu=hdu) + os.remove(tmpfile) + + +@pytest.mark.parametrize( + "file_obj, with_wl", + [ + ("mwm-temp", False), + ("mwm-temp", True), + ], +) +def test_mwm_1d_fail(file_obj, with_wl): + """Test mwm Spectrum1D loader fail on empty""" + tmpfile = str(file_obj) + ".fits" + mwm_HDUList([0, 0, 0, 0], with_wl).writeto(tmpfile, overwrite=True) + + with pytest.raises(ValueError): + Spectrum1D.read(tmpfile) + os.remove(tmpfile) + + @pytest.mark.parametrize( "file_obj, with_wl", [ @@ -478,11 +547,39 @@ def test_mwm_list_fail(file_obj, with_wl): tmpfile = str(file_obj) + ".fits" mwm_HDUList([0, 0, 0, 0], with_wl).writeto(tmpfile, overwrite=True) - with pytest.raises(RuntimeError): + with pytest.raises(ValueError): SpectrumList.read(tmpfile) os.remove(tmpfile) +@pytest.mark.parametrize( + "file_obj,n_spectra", + [ + ("spec-temp", 1), + ("spec-temp", 5), + ], +) +def test_spec_1d(file_obj, n_spectra): + """Test BOSS spec loader""" + tmpfile = str(file_obj) + ".fits" + spec_HDUList(n_spectra).writeto(tmpfile, overwrite=True) + + idxs = [1] + list(np.arange(5, 5 + n_spectra, 1)) + + for i in idxs: + data = Spectrum1D.read(tmpfile, hdu=i) + assert isinstance(data, Spectrum1D) + assert len(data.flux.value) == 10 + assert data.flux.unit == Unit("1e-17 erg / (s cm2 Angstrom)") + + assert len(data.spectral_axis.value) == 10 + assert data.spectral_axis.unit == Angstrom + assert len(data.mask) == 10 + + assert data[i].meta["header"].get("foobar") == "barfoo" + os.remove(tmpfile) + + @pytest.mark.parametrize( "file_obj,n_spectra", [ @@ -508,6 +605,48 @@ def test_spec_list(file_obj, n_spectra): os.remove(tmpfile) +@pytest.mark.parametrize( + "file_obj,hdu", + [ + ("spec-temp", 2), + ("spec-temp", 3), + ("spec-temp", 4), + ], + ids=["Fail on HDU2", "Fail on HDU3", "Fail on HDU4"], +) +def test_spec_1d_fail_hdu(file_obj, hdu): + """Test if fail on reading HDU2, HDU3, or HDU4 (non-spectra)""" + tmpfile = str(file_obj) + ".fits" + spec_HDUList(5).writeto(tmpfile, overwrite=True) + + with pytest.raises(ValueError): + Spectrum1D.read(tmpfile, hdu=hdu) + os.remove(tmpfile) + + +@pytest.mark.parametrize( + "file_obj,idx", + [ + ("apStar-temp", 1), + ("apStar-temp", 4), + ], +) +def test_apStar_1D(file_obj, idx): + tmpfile = str(file_obj) + ".fits" + apStar_HDUList(6).writeto(tmpfile, overwrite=True) + + data = Spectrum1D.read(tmpfile, idx=idx) + assert isinstance(data, Spectrum1D) + assert len(data.flux.value) == 10 + assert data.flux.unit == Unit("1e-17 erg / (s cm2 Angstrom)") + + assert len(data.spectral_axis.value) == 10 + assert data.spectral_axis.unit == Angstrom + + assert data.meta["header"].get("foobar") == "barfoo" + os.remove(tmpfile) + + @pytest.mark.parametrize( "file_obj,n_spectra", [ @@ -531,6 +670,40 @@ def test_apStar_list(file_obj, n_spectra): os.remove(tmpfile) +@pytest.mark.parametrize( + "file_obj", + [ + ("apStar-temp"), + ], +) +def test_apStar_fail_list(file_obj): + """Test if fail on reading 1D apStar as list""" + tmpfile = str(file_obj) + ".fits" + apStar_HDUList(1).writeto(tmpfile, overwrite=True) + + with pytest.raises(ValueError): + SpectrumList.read(tmpfile) + os.remove(tmpfile) + + +@pytest.mark.parametrize( + "file_obj", + [ + ("apVisit-temp"), + ], +) +def test_apVisit_1D(file_obj): + tmpfile = str(file_obj) + ".fits" + apVisit_HDUList().writeto(tmpfile, overwrite=True) + + data = Spectrum1D.read(tmpfile) + assert isinstance(data, Spectrum1D) + assert np.array_equal(data.spectral_axis.value, np.arange(1, 31, 1)) + assert len(data.flux.value) == 30 + assert data.meta["header"].get("foobar") == "barfoo" + os.remove(tmpfile) + + @pytest.mark.parametrize( "file_obj", [ From 0c5fe8f38c4053911e00289344d3a10c5dbc8422 Mon Sep 17 00:00:00 2001 From: Riley Thai Date: Fri, 18 Oct 2024 17:28:10 +1100 Subject: [PATCH 13/16] feat: visit specification on mwmVisit load - readded print -> warnings conversion - can specify the visit to load on mwmVisit load. - added relevant tests for the new mwmVisit case --- specutils/io/default_loaders/sdss_v.py | 121 +++++++++++++----- .../io/default_loaders/tests/test_sdss_v.py | 65 ++++++---- 2 files changed, 124 insertions(+), 62 deletions(-) diff --git a/specutils/io/default_loaders/sdss_v.py b/specutils/io/default_loaders/sdss_v.py index e3195d31f..9ddeaed81 100644 --- a/specutils/io/default_loaders/sdss_v.py +++ b/specutils/io/default_loaders/sdss_v.py @@ -1,10 +1,12 @@ """Register reader functions for various spectral formats.""" +import warnings from typing import Optional import numpy as np from astropy.units import Unit, Quantity, Angstrom from astropy.nddata import StdDevUncertainty, InverseVariance from astropy.io.fits import HDUList, BinTableHDU, ImageHDU +from astropy.utils.exceptions import AstropyUserWarning from ...spectra import Spectrum1D, SpectrumList from ..registers import data_loader @@ -313,8 +315,6 @@ def load_sdss_apVisit_list(file_obj, **kwargs): # BOSS REDUX products (specLite, specFull, custom coadd files, etc) - - @data_loader( "SDSS-V spec", identifier=spec_sdss5_identify, @@ -340,7 +340,8 @@ def load_sdss_spec_1D(file_obj, *args, hdu: Optional[int] = None, **kwargs): """ if hdu is None: # TODO: how should we handle this -- multiple things in file, but the user cannot choose. - print('HDU not specified. Loading coadd spectrum (HDU1)') + warnings.warn('HDU not specified. Loading coadd spectrum (HDU1)', + AstropyUserWarning) hdu = 1 # defaulting to coadd # raise ValueError("HDU not specified! Please specify a HDU to load.") elif hdu in [2, 3, 4]: @@ -438,7 +439,10 @@ def _load_BOSS_HDU(hdulist: HDUList, hdu: int, **kwargs): priority=20, extensions=["fits"], ) -def load_sdss_mwm_1d(file_obj, hdu: Optional[int] = None, **kwargs): +def load_sdss_mwm_1d(file_obj, + hdu: Optional[int] = None, + visit: Optional[int] = None, + **kwargs): """ Load an unspecified spec file as a Spectrum1D. @@ -446,8 +450,10 @@ def load_sdss_mwm_1d(file_obj, hdu: Optional[int] = None, **kwargs): ---------- file_obj : str, file-like, or HDUList FITS file name, file object, or HDUList.. - hdu : int + hdu : Optional[int] Specified HDU to load. + visit : Optional[int] + Specified visit index to load. Returns ------- @@ -467,11 +473,22 @@ def load_sdss_mwm_1d(file_obj, hdu: Optional[int] = None, **kwargs): for i in range(len(hdulist)): if hdulist[i].header.get("DATASUM") != "0": hdu = i - print('HDU not specified. Loading spectrum at (HDU{})'. - format(i)) + warnings.warn( + 'HDU not specified. Loading spectrum at (HDU{})'. + format(i), AstropyUserWarning) break - return _load_mwmVisit_or_mwmStar_hdu(hdulist, hdu, **kwargs) + # load spectra + spectra_list = _load_mwmVisit_or_mwmStar_hdu(hdulist, hdu) + + # if mwmVisit AND visit not specified, load first visit in the HDU + if visit is None: + if spectra_list[0].meta['datatype'].lower() == 'mwmvisit': + warnings.warn( + 'Visit to load not specified. Loading spectrum at index 0. Specify with (visit=...)', + AstropyUserWarning) + visit = 0 + return spectra_list[visit] @data_loader( @@ -493,12 +510,8 @@ def load_sdss_mwm_list(file_obj, **kwargs): Returns ------- - SpectrumList - The spectra contained in the file, where: - Spectrum1D - A given spectra of nD flux - None - If there are no spectra for that spectrograph/observatory + SpectrumList[Spectrum1D] + A list spectra from each visit with each instrument at each observatory (mwmVisit), or the coadd from each instrument/observatory (mwmStar). """ spectra = SpectrumList() with read_fileobj_or_hdulist(file_obj, memmap=False, **kwargs) as hdulist: @@ -513,12 +526,12 @@ def load_sdss_mwm_list(file_obj, **kwargs): for hdu in range(1, len(hdulist)): if hdulist[hdu].header.get("DATASUM") == "0": # Skip zero data HDU's - # TODO: validate if we want this printed warning or not. - # it might get annoying & fill logs with useless alerts. - print("WARNING: HDU{} ({}) is empty.".format( - hdu, hdulist[hdu].name)) + warnings.warn( + "WARNING: HDU{} ({}) is empty.".format( + hdu, hdulist[hdu].name), AstropyUserWarning) continue - spectra.append(_load_mwmVisit_or_mwmStar_hdu(hdulist, hdu)) + spectra_list = _load_mwmVisit_or_mwmStar_hdu(hdulist, hdu) + [spectra.append(spectra_list[i]) for i in range(len(spectra_list))] return spectra @@ -535,8 +548,8 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): Returns ------- - Spectrum1D - The spectrum with nD flux contained in the HDU. + list[Spectrum1D] + List of spectrum with 1D flux contained in the HDU. """ if hdulist[hdu].header.get("DATASUM") == "0": @@ -571,12 +584,7 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): mask = mask != 0 # collapse shape if 1D spectra in 2D array - # NOTE: this fixes a jdaviz handling bug for 2D of shape 1, # it could be that it's expected to be parsed this way. - if flux.shape[0] == 1: - flux = flux[0] - e_flux = e_flux[0] - mask = mask[0] # Create metadata meta = dict() @@ -597,18 +605,61 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): try: meta['mjd'] = hdulist[hdu].data['mjd'] meta["datatype"] = "mwmVisit" + except KeyError: - meta["mjd"] = (str(hdulist[hdu].data["min_mjd"][0]) + "->" + - str(hdulist[hdu].data["max_mjd"][0])) + meta["min_mjd"] = str(hdulist[hdu].data["min_mjd"][0]) + meta["max_mjd"] = str(hdulist[hdu].data["max_mjd"][0]) meta["datatype"] = "mwmStar" finally: meta["name"] = hdulist[hdu].name meta["sdss_id"] = hdulist[hdu].data['sdss_id'] - return Spectrum1D( - spectral_axis=spectral_axis, - flux=flux, - uncertainty=e_flux, - mask=mask, - meta=meta, - ) + # drop back a list of Spectrum1Ds to unpack + metadicts = _split_mwm_meta_dict(meta) + return [ + Spectrum1D( + spectral_axis=spectral_axis, + flux=flux[i], + uncertainty=e_flux[i], + mask=mask[i], + meta=metadicts[i], + ) for i in range(flux.shape[0]) + ] + + +def _split_mwm_meta_dict(d): + """ + Metadata sub-loader subfunction for MWM files. + + Parameters + ---------- + d: dict + Initial metadata dictionary. + + Returns + ------- + dicts: list[dict] + List of dicts with unpacked metadata for length > 1 array objects. + + """ + # find the length of entries + N = max(len(v) if isinstance(v, np.ndarray) else 1 for v in d.values()) + + # create N dictionaries to hold the split results + dicts = [{} for _ in range(N)] + + for key, value in d.items(): + if isinstance(value, np.ndarray): + # Ensure that the length of the list matches N + if len(value) != N: + # an error case we ignore + continue + # distribute each element to the corresponding metadict + for i in range(N): + dicts[i][key] = value[i] + else: + # if it's a single object, copy it to each metadict + for i in range(N): + dicts[i][key] = value + + return dicts diff --git a/specutils/io/default_loaders/tests/test_sdss_v.py b/specutils/io/default_loaders/tests/test_sdss_v.py index 15f9916cb..56f76ba9a 100644 --- a/specutils/io/default_loaders/tests/test_sdss_v.py +++ b/specutils/io/default_loaders/tests/test_sdss_v.py @@ -8,13 +8,16 @@ from specutils import Spectrum1D, SpectrumList -def generate_apogee_hdu(observatory="APO", with_wl=True, datasum="0"): +def generate_apogee_hdu(observatory="APO", + with_wl=True, + datasum="0", + nvisits=1): wl = (10**(4.179 + 6e-6 * np.arange(8575))).reshape((1, -1)) - flux = np.zeros_like(wl) - ivar = np.zeros_like(wl) - pixel_flags = np.zeros_like(wl) - continuum = np.zeros_like(wl) - nmf_rectified_model_flux = np.zeros_like(wl) + flux = np.array([np.zeros_like(wl)] * nvisits) + ivar = np.array([np.zeros_like(wl)] * nvisits) + pixel_flags = np.array([np.zeros_like(wl)] * nvisits) + continuum = np.array([np.zeros_like(wl)] * nvisits) + nmf_rectified_model_flux = np.array([np.zeros_like(wl)] * nvisits) columns = [ fits.Column(name="spectrum_pk_id", array=[159783564], format="K"), @@ -79,10 +82,13 @@ def generate_apogee_hdu(observatory="APO", with_wl=True, datasum="0"): return fits.BinTableHDU.from_columns(columns, header=header) -def generate_boss_hdu(observatory="APO", with_wl=True, datasum="0"): +def generate_boss_hdu(observatory="APO", with_wl=True, datasum="0", nvisits=1): wl = (10**(3.5523 + 1e-4 * np.arange(4648))).reshape((1, -1)) - flux = ivar = continuum = pixel_flags = nmf_rectified_model_flux = np.zeros_like( - wl) + flux = np.array([np.zeros_like(wl)] * nvisits) + ivar = np.array([np.zeros_like(wl)] * nvisits) + pixel_flags = np.array([np.zeros_like(wl)] * nvisits) + continuum = np.array([np.zeros_like(wl)] * nvisits) + nmf_rectified_model_flux = np.array([np.zeros_like(wl)] * nvisits) columns = [ fits.Column(name="spectrum_pk_id", array=[0], format="K"), fits.Column(name="release", array=["sdss5"], format="5A"), @@ -285,7 +291,7 @@ def fake_primary_hdu(): ])) -def mwm_HDUList(hduflags, with_wl): +def mwm_HDUList(hduflags, with_wl, **kwargs): hdulist = [fake_primary_hdu()] for i, flag in enumerate(hduflags): obs = ["APO", "LCO"] @@ -293,12 +299,14 @@ def mwm_HDUList(hduflags, with_wl): hdulist.append( generate_boss_hdu(obs[i % 2], with_wl=with_wl, - datasum=str(flag))) + datasum=str(flag), + **kwargs)) else: hdulist.append( generate_apogee_hdu(obs[i % 2], with_wl=with_wl, - datasum=str(flag))) + datasum=str(flag), + **kwargs)) print(hdulist) return fits.HDUList(hdulist) @@ -426,23 +434,25 @@ def spec_HDUList(n_spectra): # TEST MWM loaders @pytest.mark.parametrize( - "file_obj, hdu, with_wl, hduflags", + "file_obj, hdu, visit, with_wl, hduflags, nvisits", [ - ("mwm-temp", None, False, [0, 0, 1, 0]), - ("mwm-temp", 3, False, [0, 0, 1, 0]), - ("mwm-temp", None, True, [0, 1, 1, 0]), - ("mwm-temp", 2, True, [0, 1, 1, 0]), + ("mwm-temp", None, None, False, [0, 0, 1, 0], 1), # visit + ("mwm-temp", 3, None, False, [0, 0, 1, 0], 1), + ("mwm-temp", None, 2, False, [0, 1, 1, 0], 3), # multi-ext visits + ("mwm-temp", 2, 2, False, [0, 1, 1, 0], 3), + ("mwm-temp", None, None, True, [0, 0, 1, 0], 1), # star + ("mwm-temp", 3, None, True, [0, 0, 1, 0], 1), + ("mwm-temp", None, None, True, [0, 1, 1, 0], 1), + ("mwm-temp", 2, None, True, [0, 1, 1, 0], 1), ], ) -def test_mwm_1d(file_obj, hdu, with_wl, hduflags): +def test_mwm_1d(file_obj, hdu, visit, with_wl, hduflags, nvisits): """Test mwm Spectrum1D loader""" tmpfile = str(file_obj) + ".fits" - mwm_HDUList(hduflags, with_wl).writeto(tmpfile, overwrite=True) + mwm_HDUList(hduflags, with_wl, nvisits=nvisits).writeto(tmpfile, + overwrite=True) - if hdu is None: - data = Spectrum1D.read(tmpfile) - else: - data = Spectrum1D.read(tmpfile, hdu=hdu) + data = Spectrum1D.read(tmpfile, hdu=hdu, visit=visit) assert isinstance(data, Spectrum1D) assert isinstance(data.meta["header"], fits.Header) if data.meta["instrument"].lower() == "apogee": @@ -463,19 +473,20 @@ def test_mwm_1d(file_obj, hdu, with_wl, hduflags): "file_obj, with_wl, hduflags", [ ("mwm-temp", False, [0, 0, 1, 1]), - ("mwm-temp", True, [0, 0, 1, 1]), ("mwm-temp", False, [0, 1, 1, 0]), - ("mwm-temp", True, [0, 1, 1, 0]), ("mwm-temp", False, [1, 1, 0, 0]), - ("mwm-temp", True, [1, 1, 0, 0]), ("mwm-temp", False, [1, 1, 1, 1]), + ("mwm-temp", True, [0, 0, 1, 1]), + ("mwm-temp", True, [0, 1, 1, 0]), + ("mwm-temp", True, [1, 1, 0, 0]), ("mwm-temp", True, [1, 1, 1, 1]), ], ) def test_mwm_list(file_obj, with_wl, hduflags): """Test mwm SpectrumList loader""" tmpfile = str(file_obj) + ".fits" - mwm_HDUList(hduflags, with_wl).writeto(tmpfile, overwrite=True) + mwm_HDUList(hduflags, with_wl, + nvisits=1 if with_wl else 3).writeto(tmpfile, overwrite=True) data = SpectrumList.read(tmpfile) assert isinstance(data, SpectrumList) From ca4d8c3394f189a8197c225f131b6efa6c5a8570 Mon Sep 17 00:00:00 2001 From: Riley Thai Date: Sat, 19 Oct 2024 13:19:23 +1100 Subject: [PATCH 14/16] revert: completely rollback to initial case - no longer unpacks nD spectrum in all cases for mwm - will still unpack flux from 2D to 1D in the event its single visit or coadd --- specutils/io/default_loaders/sdss_v.py | 78 +++++--------------------- 1 file changed, 15 insertions(+), 63 deletions(-) diff --git a/specutils/io/default_loaders/sdss_v.py b/specutils/io/default_loaders/sdss_v.py index 9ddeaed81..014e3b1f2 100644 --- a/specutils/io/default_loaders/sdss_v.py +++ b/specutils/io/default_loaders/sdss_v.py @@ -470,7 +470,7 @@ def load_sdss_mwm_1d(file_obj, # TODO: how should we handle this -- multiple things in file, but the user cannot choose. if hdu is None: - for i in range(len(hdulist)): + for i in range(1, len(hdulist)): if hdulist[i].header.get("DATASUM") != "0": hdu = i warnings.warn( @@ -478,17 +478,8 @@ def load_sdss_mwm_1d(file_obj, format(i), AstropyUserWarning) break - # load spectra - spectra_list = _load_mwmVisit_or_mwmStar_hdu(hdulist, hdu) - - # if mwmVisit AND visit not specified, load first visit in the HDU - if visit is None: - if spectra_list[0].meta['datatype'].lower() == 'mwmvisit': - warnings.warn( - 'Visit to load not specified. Loading spectrum at index 0. Specify with (visit=...)', - AstropyUserWarning) - visit = 0 - return spectra_list[visit] + # load spectra and return + return _load_mwmVisit_or_mwmStar_hdu(hdulist, hdu) @data_loader( @@ -530,8 +521,7 @@ def load_sdss_mwm_list(file_obj, **kwargs): "WARNING: HDU{} ({}) is empty.".format( hdu, hdulist[hdu].name), AstropyUserWarning) continue - spectra_list = _load_mwmVisit_or_mwmStar_hdu(hdulist, hdu) - [spectra.append(spectra_list[i]) for i in range(len(spectra_list))] + spectra.append(_load_mwmVisit_or_mwmStar_hdu(hdulist, hdu)) return spectra @@ -585,6 +575,10 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): # collapse shape if 1D spectra in 2D array # it could be that it's expected to be parsed this way. + if flux.shape[0] == 1: + flux = np.ravel(flux) + e_flux = e_flux[0] # different class + mask = np.ravel(mask) # Create metadata meta = dict() @@ -605,7 +599,6 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): try: meta['mjd'] = hdulist[hdu].data['mjd'] meta["datatype"] = "mwmVisit" - except KeyError: meta["min_mjd"] = str(hdulist[hdu].data["min_mjd"][0]) meta["max_mjd"] = str(hdulist[hdu].data["max_mjd"][0]) @@ -615,51 +608,10 @@ def _load_mwmVisit_or_mwmStar_hdu(hdulist: HDUList, hdu: int, **kwargs): meta["sdss_id"] = hdulist[hdu].data['sdss_id'] # drop back a list of Spectrum1Ds to unpack - metadicts = _split_mwm_meta_dict(meta) - return [ - Spectrum1D( - spectral_axis=spectral_axis, - flux=flux[i], - uncertainty=e_flux[i], - mask=mask[i], - meta=metadicts[i], - ) for i in range(flux.shape[0]) - ] - - -def _split_mwm_meta_dict(d): - """ - Metadata sub-loader subfunction for MWM files. - - Parameters - ---------- - d: dict - Initial metadata dictionary. - - Returns - ------- - dicts: list[dict] - List of dicts with unpacked metadata for length > 1 array objects. - - """ - # find the length of entries - N = max(len(v) if isinstance(v, np.ndarray) else 1 for v in d.values()) - - # create N dictionaries to hold the split results - dicts = [{} for _ in range(N)] - - for key, value in d.items(): - if isinstance(value, np.ndarray): - # Ensure that the length of the list matches N - if len(value) != N: - # an error case we ignore - continue - # distribute each element to the corresponding metadict - for i in range(N): - dicts[i][key] = value[i] - else: - # if it's a single object, copy it to each metadict - for i in range(N): - dicts[i][key] = value - - return dicts + return Spectrum1D( + spectral_axis=spectral_axis, + flux=flux, + uncertainty=e_flux, + mask=mask, + meta=meta, + ) From 9960f980022a308237a91b4b3a129d4d8655df1f Mon Sep 17 00:00:00 2001 From: Riley Thai Date: Sat, 26 Oct 2024 14:10:12 +1100 Subject: [PATCH 15/16] fix: test cases - fixed flux array length check so that the 2D mwm visits are checked properly --- .../io/default_loaders/tests/test_sdss_v.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/specutils/io/default_loaders/tests/test_sdss_v.py b/specutils/io/default_loaders/tests/test_sdss_v.py index 56f76ba9a..5be49cc86 100644 --- a/specutils/io/default_loaders/tests/test_sdss_v.py +++ b/specutils/io/default_loaders/tests/test_sdss_v.py @@ -438,7 +438,7 @@ def spec_HDUList(n_spectra): [ ("mwm-temp", None, None, False, [0, 0, 1, 0], 1), # visit ("mwm-temp", 3, None, False, [0, 0, 1, 0], 1), - ("mwm-temp", None, 2, False, [0, 1, 1, 0], 3), # multi-ext visits + ("mwm-temp", None, 2, False, [0, 1, 1, 0], 2), # multi-ext visits ("mwm-temp", 2, 2, False, [0, 1, 1, 0], 3), ("mwm-temp", None, None, True, [0, 0, 1, 0], 1), # star ("mwm-temp", 3, None, True, [0, 0, 1, 0], 1), @@ -463,7 +463,10 @@ def test_mwm_1d(file_obj, hdu, visit, with_wl, hduflags, nvisits): raise ValueError( "INSTRMNT tag in test HDU header is not set properly.") assert len(data.spectral_axis.value) == length - assert len(data.flux.value) == length + assert data.flux.value.shape[-1] == length + if nvisits > 1: + assert data.flux.value.shape[0] == nvisits + assert data.spectral_axis.unit == Angstrom assert data.flux.unit == Unit("1e-17 erg / (s cm2 Angstrom)") os.remove(tmpfile) @@ -485,8 +488,9 @@ def test_mwm_1d(file_obj, hdu, visit, with_wl, hduflags, nvisits): def test_mwm_list(file_obj, with_wl, hduflags): """Test mwm SpectrumList loader""" tmpfile = str(file_obj) + ".fits" - mwm_HDUList(hduflags, with_wl, - nvisits=1 if with_wl else 3).writeto(tmpfile, overwrite=True) + nvisits = 1 if with_wl else 3 + mwm_HDUList(hduflags, with_wl, nvisits=nvisits).writeto(tmpfile, + overwrite=True) data = SpectrumList.read(tmpfile) assert isinstance(data, SpectrumList) @@ -505,7 +509,9 @@ def test_mwm_list(file_obj, with_wl, hduflags): else: assert data[i].meta['datatype'].lower() == 'mwmvisit' assert len(data[i].spectral_axis.value) == length - assert len(data[i].flux.value) == length + assert data[i].flux.value.shape[-1] == length + if nvisits > 1: + assert data[i].flux.value.shape[0] == nvisits assert data[i].spectral_axis.unit == Angstrom assert data[i].flux.unit == Unit("1e-17 erg / (s cm2 Angstrom)") os.remove(tmpfile) From 87fa13ffa63677d887b2b8b821d53c585676f29b Mon Sep 17 00:00:00 2001 From: Riley Thai Date: Sat, 26 Oct 2024 14:28:00 +1100 Subject: [PATCH 16/16] fix: split test cases for user warning + remove useless warning - removed HDU is empty warning on mwm SpecList, not really worth a warning - moved test cases which are designed to throw exceptions into their own test function --- specutils/io/default_loaders/sdss_v.py | 3 - .../io/default_loaders/tests/test_sdss_v.py | 56 +++++++++++++++---- 2 files changed, 45 insertions(+), 14 deletions(-) diff --git a/specutils/io/default_loaders/sdss_v.py b/specutils/io/default_loaders/sdss_v.py index 014e3b1f2..ab5f92d45 100644 --- a/specutils/io/default_loaders/sdss_v.py +++ b/specutils/io/default_loaders/sdss_v.py @@ -517,9 +517,6 @@ def load_sdss_mwm_list(file_obj, **kwargs): for hdu in range(1, len(hdulist)): if hdulist[hdu].header.get("DATASUM") == "0": # Skip zero data HDU's - warnings.warn( - "WARNING: HDU{} ({}) is empty.".format( - hdu, hdulist[hdu].name), AstropyUserWarning) continue spectra.append(_load_mwmVisit_or_mwmStar_hdu(hdulist, hdu)) return spectra diff --git a/specutils/io/default_loaders/tests/test_sdss_v.py b/specutils/io/default_loaders/tests/test_sdss_v.py index 5be49cc86..fae2bf1ff 100644 --- a/specutils/io/default_loaders/tests/test_sdss_v.py +++ b/specutils/io/default_loaders/tests/test_sdss_v.py @@ -1,9 +1,11 @@ import os +import warnings # noqa ; required for pytest import numpy as np import pytest from astropy.io import fits from astropy.units import Angstrom, Unit +from astropy.utils.exceptions import AstropyUserWarning from specutils import Spectrum1D, SpectrumList @@ -432,27 +434,59 @@ def spec_HDUList(n_spectra): return hdulist +@pytest.mark.parametrize( + "file_obj, hdu, with_wl, hduflags, nvisits", + [ + ("mwm-temp", None, False, [0, 0, 1, 0], 1), # visit + ("mwm-temp", None, False, [0, 1, 1, 0], 3), # multi-ext visits + ("mwm-temp", None, True, [0, 0, 1, 0], 1), # star + ("mwm-temp", None, True, [0, 1, 1, 0], 1), + ], +) +def test_mwm_1d_nohdu(file_obj, hdu, with_wl, hduflags, nvisits): + """Test mwm Spectrum1D loader when HDU isn't specified""" + tmpfile = str(file_obj) + ".fits" + mwm_HDUList(hduflags, with_wl, nvisits=nvisits).writeto(tmpfile, + overwrite=True) + + with pytest.warns(AstropyUserWarning): + data = Spectrum1D.read(tmpfile, hdu=hdu) + assert isinstance(data, Spectrum1D) + assert isinstance(data.meta["header"], fits.Header) + if data.meta["instrument"].lower() == "apogee": + length = 8575 + elif data.meta["instrument"].lower() == "boss": + length = 4648 + else: + raise ValueError( + "INSTRMNT tag in test HDU header is not set properly.") + assert len(data.spectral_axis.value) == length + assert data.flux.value.shape[-1] == length + if nvisits > 1: + assert data.flux.value.shape[0] == nvisits + + assert data.spectral_axis.unit == Angstrom + assert data.flux.unit == Unit("1e-17 erg / (s cm2 Angstrom)") + os.remove(tmpfile) + + # TEST MWM loaders @pytest.mark.parametrize( - "file_obj, hdu, visit, with_wl, hduflags, nvisits", + "file_obj, hdu, with_wl, hduflags, nvisits", [ - ("mwm-temp", None, None, False, [0, 0, 1, 0], 1), # visit - ("mwm-temp", 3, None, False, [0, 0, 1, 0], 1), - ("mwm-temp", None, 2, False, [0, 1, 1, 0], 2), # multi-ext visits - ("mwm-temp", 2, 2, False, [0, 1, 1, 0], 3), - ("mwm-temp", None, None, True, [0, 0, 1, 0], 1), # star - ("mwm-temp", 3, None, True, [0, 0, 1, 0], 1), - ("mwm-temp", None, None, True, [0, 1, 1, 0], 1), - ("mwm-temp", 2, None, True, [0, 1, 1, 0], 1), + ("mwm-temp", 3, False, [0, 0, 1, 0], 1), + ("mwm-temp", 3, False, [0, 1, 1, 0], 5), + ("mwm-temp", 3, True, [0, 0, 1, 0], 1), + ("mwm-temp", 2, True, [0, 1, 1, 0], 1), ], ) -def test_mwm_1d(file_obj, hdu, visit, with_wl, hduflags, nvisits): +def test_mwm_1d(file_obj, hdu, with_wl, hduflags, nvisits): """Test mwm Spectrum1D loader""" tmpfile = str(file_obj) + ".fits" mwm_HDUList(hduflags, with_wl, nvisits=nvisits).writeto(tmpfile, overwrite=True) - data = Spectrum1D.read(tmpfile, hdu=hdu, visit=visit) + data = Spectrum1D.read(tmpfile, hdu=hdu) assert isinstance(data, Spectrum1D) assert isinstance(data.meta["header"], fits.Header) if data.meta["instrument"].lower() == "apogee":