diff --git a/docs/mosviz/plugins.rst b/docs/mosviz/plugins.rst index 6317eef88c..7337a1a980 100644 --- a/docs/mosviz/plugins.rst +++ b/docs/mosviz/plugins.rst @@ -69,4 +69,4 @@ The :guilabel:`Remove` button can be used to remove a slit once it has been appl In order to plot a slit onto the image viewer, we need WCS information from an image and slit position from a 2D spectrum. WCS information is taken from the `meta` attribute of the :class:`~astropy.nddata.CCDData` object representing the data in the image viewer. The slit position is calculated using the `S_REGION` header extension value, located in the `meta` attribute of -the :class:`~spectral_cube.SpectralCube` object that is active in the 2D spectrum viewer. +the :class:`~specutils.Spectrum1D` object that is active in the 2D spectrum viewer. diff --git a/jdaviz/configs/mosviz/plugins/parsers.py b/jdaviz/configs/mosviz/plugins/parsers.py index fa7bff5186..fabf80bb4b 100644 --- a/jdaviz/configs/mosviz/plugins/parsers.py +++ b/jdaviz/configs/mosviz/plugins/parsers.py @@ -5,12 +5,12 @@ from pathlib import Path import numpy as np -from asdf.fits_embed import AsdfInFits +from astropy import units as u from astropy.io import fits +from astropy.io.registry import IORegistryError from astropy.wcs import WCS from glue.core.data import Data from glue.core.link_helpers import LinkSame -from spectral_cube import SpectralCube from specutils import Spectrum1D, SpectrumList, SpectrumCollection from jdaviz.configs.imviz.plugins.parsers import get_image_data_iterator @@ -110,14 +110,29 @@ def link_data_in_table(app, data_obj=None): # context manager, b) handling intra-row linkage of 1D and 2D spectra in a # loop, and c) handling inter-row linkage after that in one fell swoop. with app.data_collection.delay_link_manager_update(): - for index in range(len(mos_data.get_component('1D Spectra').data)): - spec_1d = mos_data.get_component('1D Spectra').data[index] - spec_2d = mos_data.get_component('2D Spectra').data[index] + + spectra_1d = mos_data.get_component('1D Spectra').data + spectra_2d = mos_data.get_component('2D Spectra').data + + # Link each 1D spectrum with its corresponding 2D spectrum + for index in range(len(spectra_1d)): + + spec_1d = spectra_1d[index] + spec_2d = spectra_2d[index] wc_spec_1d = app.session.data_collection[spec_1d].world_component_ids wc_spec_2d = app.session.data_collection[spec_2d].world_component_ids - wc_spec_ids.append(LinkSame(wc_spec_1d[0], wc_spec_2d[0])) + wc_spec_ids.append(LinkSame(wc_spec_1d[0], wc_spec_2d[1])) + + # Link each 1D spectrum to all other 1D spectra + first_spec_1d = spectra_1d[0] + wc_first_spec_1d = app.session.data_collection[first_spec_1d].world_component_ids + + for index in range(1, len(spectra_1d)): + spec_1d = spectra_1d[index] + wc_spec_1d = app.session.data_collection[spec_1d].world_component_ids + wc_spec_ids.append(LinkSame(wc_spec_1d[0], wc_first_spec_1d[0])) app.session.data_collection.add_link(wc_spec_ids) @@ -227,77 +242,59 @@ def mos_spec2d_parser(app, data_obj, data_labels=None, add_to_table=True, data_labels : str, optional The label applied to the glue data component. """ - # In the case where the data object is a string, attempt to parse it as - # a fits file. - # TODO: this current does not handle the case where the file in the path is - # anything but a fits file whose wcs can be extracted. - def _parse_as_cube(path): + def _parse_as_spectrum1d(path): + # Parse as a FITS file and assume the WCS is correct with fits.open(path) as hdulist: data = hdulist[1].data header = hdulist[1].header - if header['NAXIS'] == 2: - new_data = np.expand_dims(data, axis=1) - header['NAXIS'] = 3 - - header['NAXIS3'] = 1 - header['BUNIT'] = 'dN/s' - header['CUNIT3'] = 'um' - header['CTYPE3'] = 'WAVE' - - # Information not present in the SCI header has to be put there - # so spectral_cube won't choke. We cook up a simple linear wcs - # with the only intention of making the code run beyond the - # spectral_cube processing. There is no guarantee that this will - # result in the correct axis label values being displayed. - # - # This is a stopgap solution that will be replaced when specutils - # absorbs the functionality provided by spectral_cube. - - fa = AsdfInFits.open(path) - gwcs = fa.tree['meta']['wcs'] - - header['CTYPE1'] = 'RA---TAN' - header['CTYPE2'] = 'DEC--TAN' - header['CUNIT1'] = 'deg' - header['CUNIT2'] = 'deg' - - header['CRVAL1'] = gwcs.forward_transform.lon_4.value - header['CRVAL2'] = gwcs.forward_transform.lat_4.value - header['CRPIX1'] = gwcs.forward_transform.intercept_1.value - header['CRPIX2'] = gwcs.forward_transform.intercept_2.value - header['CDELT1'] = gwcs.forward_transform.slope_1.value - header['CDELT2'] = gwcs.forward_transform.slope_2.value - header['PC1_1'] = -1. - header['PC1_2'] = 0. - header['PC2_1'] = 0. - header['PC2_2'] = 1. - header['PC3_1'] = 1. - header['PC3_2'] = 0. - wcs = WCS(header) + return Spectrum1D(data, wcs=wcs) - meta = {'S_REGION': header['S_REGION'], 'INSTRUME': 'nirspec'} - - return SpectralCube(new_data, wcs=wcs, meta=meta) - - if isinstance(data_labels, str): - data_labels = [data_labels] - - # Coerce into list if needed + # Coerce into list-like object if not isinstance(data_obj, (list, tuple, SpectrumCollection)): data_obj = [data_obj] - data_obj = [_parse_as_cube(x) if _check_is_file(x) else x for x in data_obj] - - if data_labels is None: - data_labels = [f"2D Spectrum {i}" for i in range(len(data_obj))] - elif len(data_obj) != len(data_labels): - data_labels = [f"{data_labels} {i}" for i in range(len(data_obj))] + # If we're given a string, repeat it for each object + if isinstance(data_labels, str): + if len(data_obj) > 1: + data_labels = [f"{data_labels} {i}" for i in range(len(data_obj))] + else: + data_labels = [data_labels] + elif data_labels is None: + if len(data_obj) > 1: + data_labels = [f"2D Spectrum {i}" for i in range(len(data_obj))] + else: + data_labels = ['2D Spectrum'] with app.data_collection.delay_link_manager_update(): - for i in range(len(data_obj)): - app.add_data(data_obj[i], data_labels[i], notify_done=False) + for index, data in enumerate(data_obj): + # If we got a filepath, first try and parse using the Spectrum1D and + # SpectrumList parsers, and then fall back to parsing it as a generic + # FITS file. + if _check_is_file(data): + try: + data = Spectrum1D.read(data) + except IORegistryError: + try: + data = Spectrum1D.read(data) + except IORegistryError: + data = _parse_as_spectrum1d(data) + + # Copy (if present) region to top-level meta object + if ('header' in data.meta and + 'S_REGION' in data.meta['header'] and + 'S_REGION' not in data.meta): + data.meta['S_REGION'] = data.meta['header']['S_REGION'] + + # Set the instrument + # TODO: this should not be set to nirspec for all datasets + data.meta['INSTRUME'] = 'nirspec' + + # Get the corresponding label for this data product + label = data_labels[index] + + app.data_collection[label] = data if add_to_table: _add_to_table(app, data_labels, '2D Spectra') @@ -542,39 +539,27 @@ def mos_niriss_parser(app, data_dir, obs_label=""): with fits.open(fname, memmap=False) as temp: sci_hdus = [] + wav_hdus = {} for i in range(len(temp)): if "EXTNAME" in temp[i].header: if temp[i].header["EXTNAME"] == "SCI": sci_hdus.append(i) + wav_hdus[i] = ('WAVELENGTH', temp[i].header['EXTVER']) - # Now get a SpectralCube object for each SCI HDU + # Now get a Spectrum1D object for each SCI HDU for sci in sci_hdus: if temp[sci].header["SPORDER"] == 1: data = temp[sci].data meta = temp[sci].header - header = filter_wcs[filter_name] - - # Information that needs to be added to the header in order to load - # WCS information into SpectralCube. - # TODO: Use gwcs instead to avoid hardcoding information for 3rd wcs axis - new_data = np.expand_dims(data, axis=1) - header['NAXIS'] = 3 - - header['NAXIS3'] = 1 - header['BUNIT'] = 'dN/s' - header['CUNIT3'] = 'um' - header['WCSAXES'] = 3 - header['CRPIX3'] = 0.0 - header['CDELT3'] = 1E-06 - header['CUNIT3'] = 'm' - header['CTYPE3'] = 'WAVE' - header['CRVAL3'] = 0.0 - - wcs = WCS(header) - - spec2d = SpectralCube(new_data, wcs=wcs, meta=meta) + + # The wavelength is stored in a WAVELENGTH HDU. This is + # a 2D array, but in order to be able to use Spectrum1D + # we use the average wavelength for all image rows + wav = temp[wav_hdus[sci]].data.mean(axis=0) * u.micron + + spec2d = Spectrum1D(data * u.one, spectral_axis=wav, meta=meta) spec2d.meta['INSTRUME'] = 'NIRISS' diff --git a/jdaviz/configs/mosviz/plugins/viewers.py b/jdaviz/configs/mosviz/plugins/viewers.py index e1e5d09109..2dc8d61850 100644 --- a/jdaviz/configs/mosviz/plugins/viewers.py +++ b/jdaviz/configs/mosviz/plugins/viewers.py @@ -3,8 +3,6 @@ from glue_jupyter.bqplot.profile import BqplotProfileView from glue_jupyter.table import TableViewer from specutils import Spectrum1D -from spectral_cube import SpectralCube -from echo import delay_callback import astropy from astropy import units as u from astropy.utils.introspection import minversion @@ -78,7 +76,7 @@ def set_plot_axes(self): class MosvizProfile2DView(BqplotImageView): # Due to limitations in CCDData and 2D data that has spectral and spatial # axes, the default conversion class must handle cubes - default_class = SpectralCube + default_class = Spectrum1D tools = ['bqplot:panzoom_x'] @@ -86,7 +84,6 @@ def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) # Setup viewer option defaults self.state.aspect = 'auto' - self.state.add_callback('reference_data', self._update_world_axes, priority=100) def data(self, cls=None): return [layer_state.layer.get_object(cls=cls or self.default_class) @@ -94,13 +91,6 @@ def data(self, cls=None): if hasattr(layer_state, 'layer') and isinstance(layer_state.layer, BaseData)] - def _update_world_axes(self, data): - if data is not None: - with delay_callback(self.state, 'x_att_world', 'y_att_world'): - if 'Wave' in data.components: - self.state.x_att_world = data.id['Right Ascension'] - self.state.y_att_world = data.id['Wave'] - def set_plot_axes(self): self.figure.axes[0].visible = False diff --git a/jdaviz/configs/mosviz/tests/test_helper.py b/jdaviz/configs/mosviz/tests/test_helper.py index 39a76220c2..28232d761e 100644 --- a/jdaviz/configs/mosviz/tests/test_helper.py +++ b/jdaviz/configs/mosviz/tests/test_helper.py @@ -5,7 +5,6 @@ from astropy.nddata import CCDData from astropy.wcs import WCS from jdaviz.configs.mosviz.helper import Mosviz -from spectral_cube import SpectralCube from specutils import Spectrum1D, SpectrumCollection @@ -27,28 +26,17 @@ def spectrum1d(): @pytest.fixture def spectrum2d(): header = """ -WCSAXES = 3 / Number of coordinate axes -CRPIX1 = 1024.5 / Pixel coordinate of reference point +WCSAXES = 2 / Number of coordinate axes +CRPIX1 = 0.0 / Pixel coordinate of reference point CRPIX2 = 1024.5 / Pixel coordinate of reference point -CRPIX3 = 0.0 / Pixel coordinate of reference point -PC1_1 = -1.0 / Coordinate transformation matrix element -PC3_1 = 1.0 / Coordinate transformation matrix element -CDELT1 = 2.8685411111111E-05 / [deg] Coordinate increment at reference point +CDELT1 = 1E-06 / [m] Coordinate increment at reference point CDELT2 = 2.9256727777778E-05 / [deg] Coordinate increment at reference point -CDELT3 = 1E-06 / [m] Coordinate increment at reference point -CUNIT1 = 'deg' / Units of coordinate increment and value +CUNIT1 = 'm' / Units of coordinate increment and value CUNIT2 = 'deg' / Units of coordinate increment and value -CUNIT3 = 'm' / Units of coordinate increment and value -CTYPE1 = 'RA---TAN' / Right ascension, gnomonic projection -CTYPE2 = 'DEC--TAN' / Declination, gnomonic projection -CTYPE3 = 'WAVE' / Vacuum wavelength (linear) -CRVAL1 = 5.0 / [deg] Coordinate value at reference point +CTYPE1 = 'WAVE' / Vacuum wavelength (linear) +CTYPE2 = 'OFFSET' / Spatial offset +CRVAL1 = 0.0 / [m] Coordinate value at reference point CRVAL2 = 5.0 / [deg] Coordinate value at reference point -CRVAL3 = 0.0 / [m] Coordinate value at reference point -LONPOLE = 180.0 / [deg] Native longitude of celestial pole -LATPOLE = 5.0 / [deg] Native latitude of celestial pole -MJDREFI = 0.0 / [d] MJD of fiducial time, integer part -MJDREFF = 0.0 / [d] MJD of fiducial time, fractional part RADESYS = 'ICRS' / Equatorial coordinate system SPECSYS = 'BARYCENT' / Reference frame of spectral coordinates """ @@ -67,8 +55,8 @@ def spectrum2d(): new_hdr[key] = value wcs = WCS(new_hdr) - data = np.random.sample((15, 1, 1024)) - spectral_cube = SpectralCube(data, wcs=wcs) + data = np.random.sample((1024, 15)) * u.one + spectral_cube = Spectrum1D(data, wcs=wcs) return spectral_cube @@ -176,6 +164,7 @@ def test_load_list_of_spectrum1d(mosviz_app, spectrum1d): def test_load_spectrum2d(mosviz_app, spectrum2d): + label = "Test 2D Spectrum" mosviz_app.load_2d_spectra(spectrum2d, data_labels=label) @@ -187,7 +176,7 @@ def test_load_spectrum2d(mosviz_app, spectrum2d): data = mosviz_app.app.get_data_from_viewer('spectrum-2d-viewer') - assert list(data.values())[0].shape == (15, 1, 1024) + assert list(data.values())[0].shape == (1024, 15) assert list(data.keys())[0] == label diff --git a/setup.cfg b/setup.cfg index af5455a254..eb60363ad1 100644 --- a/setup.cfg +++ b/setup.cfg @@ -23,7 +23,7 @@ install_requires = bqplot>=0.12.29 bqplot-image-gl>=1.4.1 glue-core>=1.2.1 - glue-jupyter>=0.8 + glue-jupyter>=0.8.1 echo>=0.5.0 ipyvue>=1.4.1 ipyvuetify>=1.7.0 @@ -32,7 +32,7 @@ install_requires = voila>=0.2.4 pyyaml>=5.4.1 specutils>=1.4.0 - glue-astronomy>=0.2 + glue-astronomy>=0.3.2 click>=7.1.2 spectral-cube>=0.5 asteval>=0.9.23