Skip to content

Commit

Permalink
Manga cube loaded using Spectrum1D, still some things to iron out
Browse files Browse the repository at this point in the history
  • Loading branch information
javerbukh committed Apr 27, 2021
1 parent a3f28ee commit 1572a3b
Show file tree
Hide file tree
Showing 3 changed files with 29 additions and 542 deletions.
60 changes: 21 additions & 39 deletions jdaviz/configs/cubeviz/plugins/parsers.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,11 @@
import numpy as np
from astropy.io import fits
from spectral_cube import SpectralCube
from spectral_cube.io.fits import FITSReadError
from specutils import Spectrum1D

from astropy.wcs import WCS
from astropy import units as u

from jdaviz.core.registries import data_parser_registry

__all__ = ['parse_data']
Expand Down Expand Up @@ -44,13 +46,12 @@ def parse_data(app, file_obj, data_type=None, data_label=None):
# generic enough to work with other file types (e.g. ASDF). For now, this
# supports MaNGA and JWST data.
if isinstance(file_obj, fits.hdu.hdulist.HDUList):
_new_parse_hdu(app, file_obj, file_name=data_label)
_parse_hdu(app, file_obj, file_name=data_label)
elif isinstance(file_obj, str) and os.path.exists(file_obj):
file_name = os.path.basename(file_obj)

with fits.open(file_obj) as hdulist:
# hdulist = fits.open(file_obj)
_new_parse_hdu(app, hdulist, file_name=data_label or file_name)
_parse_hdu(app, hdulist, file_name=data_label or file_name)

# If the data types are custom data objects, use explicit parsers. Note
# that this relies on the glue-astronomy machinery to turn the data object
Expand All @@ -63,60 +64,41 @@ def parse_data(app, file_obj, data_type=None, data_label=None):
_parse_spectrum1d(app, file_obj)


def _new_parse_hdu(app, hdulist, fule_name=None):
pass


def _parse_hdu(app, hdulist, file_name=None):
if file_name is None:
if hasattr(hdulist, 'file_name'):
file_name = hdulist.file_name

file_name = file_name or "Unknown HDU object"

# WCS may only exist in a single extension (in this case, just the flux
# flux extention), so first find and store then wcs information.
wcs = None

for hdu in hdulist:
if hdu.data is None or not hdu.is_image:
continue

try:
sc = Spectrum1D.read(hdu)
except (ValueError, FITSReadError):
continue
else:
wcs = sc.wcs

# Now loop through and attempt to parse the fits extensions as spectral
# cube object. If the wcs fails to parse in any case, use the wcs
# information we scraped above.
for hdu in hdulist:
data_label = f"{file_name}[{hdu.name}]"

if hdu.data is None or not hdu.is_image:
if hdu.data is None or not hdu.is_image or len(hdu.data.shape) != 3:
continue

# This is supposed to fail on attempting to load anything that
# isn't cube-shaped. But it's not terribly reliable
wcs = WCS(hdu.header)

if 'BUNIT' in hdu.header:
flux = hdu.data * u.Unit(hdu.header['BUNIT'])
else:
# TODO Add warning that actual units were not found and using fake ones
logging.warn("No flux units found in hdu, using Jansky as a stand-in")
flux = hdu.data * u.Jy

flux = np.moveaxis(flux, 1, 2)
flux = np.moveaxis(flux, 0, 1)

try:
sc = Spectrum1D.read(hdu, format="MaNGA cube")
except (ValueError, OSError):
# This will fail if the parsing of the wcs does not provide
# proper celestial axes
try:
hdu.header.update(wcs.to_header())
sc = Spectrum1D.read(hdu, format="MaNGA cube")
except (ValueError, AttributeError) as e:
logging.warn(e)
continue
except FITSReadError as e:
sc = Spectrum1D(flux=flux, wcs=wcs)
app.data_collection[data_label] = sc
except Exception as e:
logging.warn(e)
continue

app.data_collection[data_label] = sc

# If the data type is some kind of integer, assume it's the mask/dq
if hdu.data.dtype in (np.int, np.uint, np.uint32) or \
any(x in hdu.name.lower() for x in EXT_TYPES['mask']):
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -45,7 +45,7 @@ def _on_linked_changed(self, event):
def _on_data_added(self, msg):
if len(msg.data.shape) == 3 and \
isinstance(msg.viewer, BqplotImageView):
self.max_value = msg.data.shape[2] - 1
self.max_value = msg.data.shape[0] - 1

if msg.viewer not in self._watched_viewers:
self._watched_viewers.append(msg.viewer)
Expand All @@ -55,7 +55,7 @@ def _on_data_added(self, msg):

def _slider_value_updated(self, value):
if len(value) > 0:
self.slider = float(value[2])
self.slider = float(value[0])

@observe('slider')
def _on_slider_updated(self, event):
Expand All @@ -66,4 +66,4 @@ def _on_slider_updated(self, event):

if self.linked:
for viewer in self._watched_viewers:
viewer.state.slices = (0, 0, value)
viewer.state.slices = (value, 0, 0)
Loading

0 comments on commit 1572a3b

Please sign in to comment.