diff --git a/changelog/ndc2.1.breaking.rst b/changelog/ndc2.1.breaking.rst new file mode 100644 index 0000000..c455c28 --- /dev/null +++ b/changelog/ndc2.1.breaking.rst @@ -0,0 +1 @@ +In IRIS spectrograph reader, all extra coords except time have been moved to the meta object. diff --git a/changelog/ndc2.1.feature.rst b/changelog/ndc2.1.feature.rst new file mode 100644 index 0000000..4761253 --- /dev/null +++ b/changelog/ndc2.1.feature.rst @@ -0,0 +1 @@ +Create a new `~sunraster.instr.iris.IRISSGMeta` metadata object. diff --git a/changelog/ndc2.2.breaking.rst b/changelog/ndc2.2.breaking.rst new file mode 100644 index 0000000..9556969 --- /dev/null +++ b/changelog/ndc2.2.breaking.rst @@ -0,0 +1 @@ +In IRIS spectrograph read, move all metadata to the meta objects of the raster cubes. diff --git a/changelog/ndc2.2.feature.rst b/changelog/ndc2.2.feature.rst new file mode 100644 index 0000000..373249c --- /dev/null +++ b/changelog/ndc2.2.feature.rst @@ -0,0 +1 @@ +Create new property `~sunraster.spectrogram.SpectrogramCube.celestial`, on `~sunraster.spectrogram.SpectrogramCube` to return a `~astropy.coordinates.SkyCoord` holding the celestial world coords of the pixels. diff --git a/changelog/ndc2.3.breaking.rst b/changelog/ndc2.3.breaking.rst new file mode 100644 index 0000000..53a7951 --- /dev/null +++ b/changelog/ndc2.3.breaking.rst @@ -0,0 +1 @@ +Remove extra_coords kwarg from `~sunraster.spectrogram.SpectrogramCube` __init__ in accordance with new ndcube 2.0 API. Extra coords can by added through the ndcube ExtraCoords.add API which is new in ndcube 2.0. diff --git a/changelog/ndc2.3.feature.rst b/changelog/ndc2.3.feature.rst new file mode 100644 index 0000000..3c25907 --- /dev/null +++ b/changelog/ndc2.3.feature.rst @@ -0,0 +1 @@ +Create new property `~sunraster.spectrogram.SpectrogramSequence.celestial`, on `~sunraster.spectrogram.SpectrogramSequence` to return a `~astropy.coordinates.SkyCoord` holding the celestial world coords of the pixels. diff --git a/sunraster/instr/iris.py b/sunraster/instr/iris.py index e3c215f..8886266 100644 --- a/sunraster/instr/iris.py +++ b/sunraster/instr/iris.py @@ -1,14 +1,17 @@ import copy - -from astropy.wcs import WCS -import numpy as np -from ndcube import NDCollection +import textwrap import astropy.units as u +import numpy as np +from astropy.coordinates import SkyCoord from astropy.io import fits from astropy.time import Time, TimeDelta +from astropy.wcs import WCS +from ndcube import Meta, NDCollection +from sunpy.coordinates import Helioprojective from sunraster import RasterSequence, SpectrogramCube +from sunraster.meta import SlitSpectrographMetaABC __all__ = ['read_iris_spectrograph_level2_fits'] @@ -52,8 +55,6 @@ def read_iris_spectrograph_level2_fits( hdulist = fits.open(filename, memmap=memmap, do_not_scale_image_data=memmap) hdulist.verify('fix') if f == 0: - # Determine number of raster positions in a scan - int(hdulist[0].header["NRASTERP"]) # Collecting the window observations. windows_in_obs = np.array([hdulist[0].header["TDESC{0}".format(i)] for i in range(1, hdulist[0].header["NWIN"] + 1)]) @@ -76,87 +77,37 @@ def read_iris_spectrograph_level2_fits( spectral_windows[missing_windows], filenames[0])) window_fits_indices = np.nonzero(np.in1d(windows_in_obs, spectral_windows))[0] + 1 - # Generate top level meta dictionary from first file - # main header. - top_meta = {"TELESCOP": hdulist[0].header["TELESCOP"], - "INSTRUME": hdulist[0].header["INSTRUME"], - "DATA_LEV": hdulist[0].header["DATA_LEV"], - "OBSID": hdulist[0].header["OBSID"], - "OBS_DESC": hdulist[0].header["OBS_DESC"], - "STARTOBS": Time(hdulist[0].header["STARTOBS"]), - "ENDOBS": Time(hdulist[0].header["ENDOBS"]), - "SAT_ROT": hdulist[0].header["SAT_ROT"] * u.deg, - "AECNOBS": int(hdulist[0].header["AECNOBS"]), - "FOVX": hdulist[0].header["FOVX"] * u.arcsec, - "FOVY": hdulist[0].header["FOVY"] * u.arcsec, - "SUMSPTRN": hdulist[0].header["SUMSPTRN"], - "SUMSPTRF": hdulist[0].header["SUMSPTRF"], - "SUMSPAT": hdulist[0].header["SUMSPAT"], - "NEXPOBS": hdulist[0].header["NEXPOBS"], - "NRASTERP": hdulist[0].header["NRASTERP"], - "KEYWDDOC": hdulist[0].header["KEYWDDOC"]} - # Initialize meta dictionary for each spectral_window - window_metas = {} - for i, window_name in enumerate(spectral_windows_req): - if "FUV" in hdulist[0].header["TDET{0}".format(window_fits_indices[i])]: - spectral_summing = hdulist[0].header["SUMSPTRF"] - else: - spectral_summing = hdulist[0].header["SUMSPTRN"] - window_metas[window_name] = { - "detector type": - hdulist[0].header["TDET{0}".format(window_fits_indices[i])], - "spectral window": - hdulist[0].header["TDESC{0}".format(window_fits_indices[i])], - "brightest wavelength": - hdulist[0].header["TWAVE{0}".format(window_fits_indices[i])], - "min wavelength": - hdulist[0].header["TWMIN{0}".format(window_fits_indices[i])], - "max wavelength": - hdulist[0].header["TWMAX{0}".format(window_fits_indices[i])], - "SAT_ROT": hdulist[0].header["SAT_ROT"], - "spatial summing": hdulist[0].header["SUMSPAT"], - "spectral summing": spectral_summing - } # Create a empty list for every spectral window and each # spectral window is a key for the dictionary. data_dict = dict([(window_name, list()) for window_name in spectral_windows_req]) - # Determine extra coords for this raster. + # Extract axis-aligned metadata. times = (Time(hdulist[0].header["STARTOBS"]) + TimeDelta(hdulist[-2].data[:, hdulist[-2].header["TIME"]], format='sec')) - np.arange(int(hdulist[0].header["NRASTERP"])) - pztx = hdulist[-2].data[:, hdulist[-2].header["PZTX"]] * u.arcsec - pzty = hdulist[-2].data[:, hdulist[-2].header["PZTY"]] * u.arcsec - xcenix = hdulist[-2].data[:, hdulist[-2].header["XCENIX"]] * u.arcsec - ycenix = hdulist[-2].data[:, hdulist[-2].header["YCENIX"]] * u.arcsec + fov_center = SkyCoord(Tx=hdulist[-2].data[:, hdulist[-2].header["XCENIX"]], + Ty=hdulist[-2].data[:, hdulist[-2].header["YCENIX"]], + unit=u.arcsec, frame=Helioprojective) obs_vrix = hdulist[-2].data[:, hdulist[-2].header["OBS_VRIX"]] * u.m / u.s ophaseix = hdulist[-2].data[:, hdulist[-2].header["OPHASEIX"]] exposure_times_fuv = hdulist[-2].data[:, hdulist[-2].header["EXPTIMEF"]] * u.s exposure_times_nuv = hdulist[-2].data[:, hdulist[-2].header["EXPTIMEN"]] * u.s - # If OBS is raster, include raster positions. Otherwise don't. - if top_meta["NRASTERP"] > 1: - general_extra_coords = [("time", 0, times), - ("raster position", 0, np.arange(top_meta["NRASTERP"])), - ("pztx", 0, pztx), ("pzty", 0, pzty), - ("xcenix", 0, xcenix), ("ycenix", 0, ycenix), - ("obs_vrix", 0, obs_vrix), ("ophaseix", 0, ophaseix)] - else: - general_extra_coords = [("time", 0, times), - ("pztx", 0, pztx), ("pzty", 0, pzty), - ("xcenix", 0, xcenix), ("ycenix", 0, ycenix), - ("obs_vrix", 0, obs_vrix), ("ophaseix", 0, ophaseix)] for i, window_name in enumerate(spectral_windows_req): + # Define metadata object for window. + meta = IRISSGMeta(hdulist[0].header, window_name, + data_shape=hdulist[window_fits_indices[i]].data.shape) # Determine values of properties dependent on detector type. - if "FUV" in hdulist[0].header["TDET{0}".format(window_fits_indices[i])]: + if "FUV" in meta.detector: exposure_times = exposure_times_fuv DN_unit = DN_UNIT["FUV"] readout_noise = READOUT_NOISE["FUV"] - elif "NUV" in hdulist[0].header["TDET{0}".format(window_fits_indices[i])]: + else: exposure_times = exposure_times_nuv DN_unit = DN_UNIT["NUV"] readout_noise = READOUT_NOISE["NUV"] - else: - raise ValueError("Detector type in FITS header not recognized.") + meta.add("exposure time", exposure_times, None, 0) + meta.add("exposure FOV center", fov_center, None, 0) + meta.add("observer radial velocity", obs_vrix, None, 0) + meta.add("orbital phase", ophaseix, None, 0) # Derive WCS, data and mask for NDCube from file. # Sit-and-stare have a CDELT of 0 which causes issues in astropy WCS. # In this case, set CDELT to a tiny non-zero number. @@ -167,36 +118,6 @@ def read_iris_spectrograph_level2_fits( data_mask = hdulist[window_fits_indices[i]].data == -200. else: data_mask = None - # Derive extra coords for this spectral window. - window_extra_coords = copy.deepcopy(general_extra_coords) - window_extra_coords.append(("exposure time", 0, exposure_times)) - # Collect metadata relevant to single files. - try: - date_obs = Time(hdulist[0].header["DATE_OBS"]) - except ValueError: - date_obs = None - try: - date_end = Time(hdulist[0].header["DATE_END"]) - except ValueError: - date_end = None - single_file_meta = {"SAT_ROT": hdulist[0].header["SAT_ROT"] * u.deg, - "DATE_OBS": date_obs, - "DATE_END": date_end, - "HLZ": bool(int(hdulist[0].header["HLZ"])), - "SAA": bool(int(hdulist[0].header["SAA"])), - "DSUN_OBS": hdulist[0].header["DSUN_OBS"] * u.m, - "IAECEVFL": hdulist[0].header["IAECEVFL"], - "IAECFLAG": hdulist[0].header["IAECFLAG"], - "IAECFLFL": hdulist[0].header["IAECFLFL"], - "KEYWDDOC": hdulist[0].header["KEYWDDOC"], - "detector type": - hdulist[0].header["TDET{0}".format(window_fits_indices[i])], - "spectral window": window_name, - "OBSID": hdulist[0].header["OBSID"], - "OBS_DESC": hdulist[0].header["OBS_DESC"], - "STARTOBS": Time(hdulist[0].header["STARTOBS"]), - "ENDOBS": Time(hdulist[0].header["ENDOBS"]) - } # Derive uncertainty of data if uncertainty: out_uncertainty = u.Quantity( @@ -206,14 +127,185 @@ def read_iris_spectrograph_level2_fits( else: out_uncertainty = None # Appending NDCube instance to the corresponding window key in dictionary's list. - data_dict[window_name].append( - SpectrogramCube(hdulist[window_fits_indices[i]].data, wcs=wcs_, - uncertainty=out_uncertainty, unit=DN_unit, meta=single_file_meta, - extra_coords=window_extra_coords, mask=data_mask)) + cube = SpectrogramCube(hdulist[window_fits_indices[i]].data, wcs=wcs_, + uncertainty=out_uncertainty, unit=DN_unit, + meta=meta, mask=data_mask) + cube.extra_coords.add("time", 0, times, physical_types="time") + data_dict[window_name].append(cube) hdulist.close() # Construct dictionary of SpectrogramSequences for spectral windows - window_data_pairs = [(window_name, RasterSequence(data_dict[window_name], common_axis=0, - meta=window_metas[window_name])) + window_data_pairs = [(window_name, RasterSequence(data_dict[window_name], common_axis=0)) for window_name in spectral_windows_req] # Initialize an NDCollection object. - return NDCollection(window_data_pairs, aligned_axes=(0, 1, 2), meta=top_meta) + return NDCollection(window_data_pairs, aligned_axes=(0, 1, 2)) + + +class IRISSGMeta(Meta, metaclass=SlitSpectrographMetaABC): + def __init__(self, header, spectral_window, **kwargs): + super().__init__(header, **kwargs) + spectral_windows = np.array([self["TDESC{0}".format(i)] + for i in range(1, self["NWIN"] + 1)]) + window_mask = np.array([spectral_window in window for window in spectral_windows]) + if window_mask.sum() < 1: + raise ValueError("Spectral window not found. " + f"Input spectral window: {spectral_window}; " + f"Spectral windows in header: {spectral_windows}") + elif window_mask.sum() > 1: + raise ValueError( + "Spectral window must be unique. " + f"Input spectral window: {spectral_window}; " + f"Ambiguous spectral windows in header: {spectral_windows[window_mask]}") + self._iwin = np.arange(len(spectral_windows))[window_mask][0] + 1 + + def __str__(self): + return textwrap.dedent(f"""\ + IRISMeta + -------- + Observatory:\t\t{self.observatory} + Instrument:\t\t{self.instrument} + Detector:\t\t{self.detector} + Spectral Window:\t{self.spectral_window} + Spectral Range:\t\t{self.spectral_range} + Date:\t\t\t{self.date_reference} + OBS ID:\t\t\t{self.observing_mode_id} + OBS Description:\t{self.observing_mode_description} + """) + + def __repr__(self): + return f"{object.__repr__(self)}\n{str(self)}" + + # ---------- IRIS-specific convenience methods ---------- + def _construct_time(self, key): + val = self.get(key, None) + if val is not None: + val = Time(val, format="fits", scale="utc") + return val + + # ---------- Inherited ABC properties ---------- + @property + def spectral_window(self): + return self.get(f"TDESC{self._iwin}") + + @property + def detector(self): + return self.get(f"TDET{self._iwin}") + + @property + def instrument(self): + return self.get("INSTRUME") + + @property + def observatory(self): + return self.get("TELESCOP") + + @property + def processing_level(self): + return self.get("DATA_LEV") + + @property + def distance_to_sun(self): + return self.get("DSUN_OBS") * u.m + + @property + def date_reference(self): + return self._construct_time("DATE_OBS") + + @property + def date_start(self): + return self.date_reference + + @property + def date_end(self): + return self._construct_time("DATE_END") + + @property + def observing_mode_id(self): + return int(self.get("OBSID")) + + # ---------- IRIS-specific metadata properties ---------- + + @property + def observing_mode_description(self): + return self.get("OBS_DESC") + + @property + def observing_campaign_start(self): + """Start time of observing campaign""" + return self._construct_time("STARTOBS") + + @property + def observing_campaign_end(self): + """End time of observing mode""" + return self._construct_time("ENDOBS") + + @property + def observation_includes_SAA(self): + """Whether IRIS passed through SAA during observations.""" + return bool(self.get("SAA")) + + @property + def satellite_rotation(self): + """Satellite roll from solar north""" + return self.get("SAT_ROT") * u.deg + + @property + def exposure_control_triggers_in_observation(self): + """Number of times automatic exposure control triggered during observing campaign.""" + return self.get("AECNOBS") + + @property + def exposure_control_triggers_in_raster(self): + """Number of times automatic exposure control was triggered during this raster.""" + return self.get("AECNRAS") + + @property + def number_raster_positions(self): + """ Number of positions in raster.""" + self.get("NRASTERP") + + @property + def spectral_range(self): + """The spectral range of the spectral window.""" + return [self.get(f"TWMIN{self._iwin}"), self.get(f"TWMAX{self._iwin}")] * u.AA + + @property + def raster_FOV_width_y(self): + """Width of the field of view of the raster in the Y (slit) direction.""" + return self.get("FOVY") * u.arcsec + + @property + def raster_FOV_width_x(self): + """Width of the field of view of the raster in the X (rastering) direction.""" + return self.get("FOVX") * u.arcsec + + @property + def FOV_center(self): + """Location of the center of the field of view.""" + return SkyCoord(Tx=self.get("XCEN"), Ty=self.get("YCEN"), unit=u.arcsec, + frame=Helioprojective) + + @property + def automatic_exposure_control_enabled(self): + return bool(self.get("IAECFLAG")) + + @property + def tracking_mode_enabled(self): + return bool(self.get("TR_MODE")) + + @property + def observatory_at_high_latitude(self): + """Whether IRIS passed through high Earth latitude during observations.""" + return bool(self.get("HLZ")) + + @property + def spatial_summing_factor(self): + """Number of pixels summed together in the spatial (Y/slit) direction.""" + return self.get("SUMSPAT") + + @property + def spectral_summing_factor(self): + """Number of pixels summed together in the spectral direction.""" + if "fuv" in self.detector.lower(): + return self.get("SUMSPTRF") + else: + return self.get("SUMSPTRN") diff --git a/sunraster/instr/spice.py b/sunraster/instr/spice.py index 6dd745b..ea42aa7 100644 --- a/sunraster/instr/spice.py +++ b/sunraster/instr/spice.py @@ -9,9 +9,9 @@ from astropy.time import Time, TimeDelta from astropy.wcs import WCS from sunpy.coordinates import HeliographicStonyhurst -from ndcube import NDCollection +from ndcube import NDCollection, Meta -from sunraster.meta import Meta, SlitSpectrographMetaABC +from sunraster.meta import SlitSpectrographMetaABC from sunraster import SpectrogramCube, SpectrogramSequence, RasterSequence diff --git a/sunraster/spectrogram.py b/sunraster/spectrogram.py index 9b2397d..f75b34b 100644 --- a/sunraster/spectrogram.py +++ b/sunraster/spectrogram.py @@ -1,12 +1,12 @@ import abc -import textwrap import numbers +import textwrap +import warnings -import numpy as np -from astropy.time import TimeDelta import astropy.units as u +import numpy as np +from astropy.time import Time from ndcube.ndcube import NDCube -from ndcube.utils.cube import convert_extra_coords_dict_to_input_format __all__ = ['SpectrogramCube'] @@ -20,7 +20,7 @@ "been undone since the unit does not include " "inverse time. To undo exposure time correction " "anyway, set 'force' kwarg to True.") -AXIS_NOT_FOUND_ERROR = " axis not found. If in extra_coords, axis name must be supported: " +AXIS_NOT_FOUND_ERROR = "axis not found. If in extra_coords, axis name must be supported:" # Define supported coordinate names for coordinate properties. SUPPORTED_LONGITUDE_NAMES = ["custom:pos.helioprojective.lon", "pos.helioprojective.lon", @@ -87,6 +87,12 @@ def lat(self): Return the latitude coordinates for each pixel. """ + @abc.abstractproperty + def celestial(self): + """ + Return the celestial coordinates for each pixel. + """ + @abc.abstractmethod def apply_exposure_time_correction(self, undo=False, force=False): """ @@ -163,11 +169,11 @@ class SpectrogramCube(NDCube, SpectrogramABC): Default is False. """ - def __init__(self, data, wcs, extra_coords=None, unit=None, uncertainty=None, meta=None, + def __init__(self, data, wcs, unit=None, uncertainty=None, meta=None, mask=None, instrument_axes=None, copy=False, **kwargs): # Initialize SpectrogramCube. super().__init__(data, wcs=wcs, uncertainty=uncertainty, mask=mask, meta=meta, - unit=unit, extra_coords=extra_coords, copy=copy, **kwargs) + unit=unit, copy=copy, **kwargs) # Determine labels and location of each key real world coordinate. self_extra_coords = self.extra_coords @@ -192,34 +198,45 @@ def __init__(self, data, wcs, extra_coords=None, unit=None, uncertainty=None, me self.instrument_axes = np.asarray(instrument_axes, dtype=str) def __str__(self): - if self._time_name: + try: if self.time.isscalar: time_period = self.time else: - time_period = (str(self.time.min()), str(self.time.max())) - else: - time_period = None - if self._longitude_name: - if self.lon.isscalar: - lon_range = self.lon + times = self.time + time_period = Time([times.min(), times.max()]).iso + except ValueError as err: + if AXIS_NOT_FOUND_ERROR in err.args[0]: + time_period = None else: - lon_range = u.Quantity([self.lon.min(), self.lon.max()]) - else: - lon_range = None - if self._latitude_name: - if self.lat.isscalar: - lat_range = self.lat + raise err + try: + sc = self.celestial + component_names = dict( + [(item, key) for key, item in sc.representation_component_names.items()]) + lon = getattr(sc, component_names["lon"]) + lat = getattr(sc, component_names["lat"]) + if sc.isscalar: + lon_range = lon + lat_range = lat else: - lat_range = u.Quantity([self.lat.min(), self.lat.max()]) - else: - lat_range = None - if self._spectral_name: + lon_range = u.Quantity([lon.min(), lon.max()]) + lat_range = u.Quantity([lat.min(), lat.max()]) + except ValueError as err: + if AXIS_NOT_FOUND_ERROR in err.args[0]: + lon_range = None + lat_range = None + else: + raise err + try: if self.spectral_axis.isscalar: spectral_range = self.spectral_axis else: spectral_range = u.Quantity([self.spectral_axis.min(), self.spectral_axis.max()]) - else: - spectral_range = None + except ValueError as err: + if AXIS_NOT_FOUND_ERROR in err.args[0]: + spectral_range = None + else: + raise err return (textwrap.dedent(f"""\ {self.__class__.__name__} {"".join(["-"] * len(self.__class__.__name__))} @@ -237,15 +254,6 @@ def __repr__(self): def __getitem__(self, item): # Slice SpectrogramCube using parent slicing. result = super().__getitem__(item) - - # As result is an NDCube, must put extra coord back into input format - # to create a SpectrogramCube. - if result.extra_coords is None: - extra_coords = None - else: - extra_coords = convert_extra_coords_dict_to_input_format(result.extra_coords, - result.missing_axes) - # Slice instrument_axes if it exists. # If item is a slice, cube dimensionality is not reduced # so instrument_axes need not be sliced. @@ -265,14 +273,14 @@ def __getitem__(self, item): # If slicing causes cube to be a scalar, set instrument_axes to None. if len(instrument_axes) == 0: instrument_axes = None - - return self.__class__(result.data, wcs=result.wcs, uncertainty=result.uncertainty, - mask=result.mask, meta=result.meta, unit=result.unit, - extra_coords=extra_coords, missing_axes=result.missing_axes, - instrument_axes=instrument_axes) + result.instrument_axes = instrument_axes + return result @property def spectral_axis(self): + if not self._spectral_name: + self._spectral_name, self._spectral_loc = _find_axis_name( + SUPPORTED_SPECTRAL_NAMES, self.wcs.world_axis_physical_types, self.extra_coords) if not self._spectral_name: raise ValueError("Spectral" + AXIS_NOT_FOUND_ERROR + f"{SUPPORTED_SPECTRAL_NAMES}") @@ -281,37 +289,58 @@ def spectral_axis(self): @property def time(self): if not self._time_name: - raise ValueError("Time" + AXIS_NOT_FOUND_ERROR + - f"{SUPPORTED_TIME_NAMES}") - times = self._get_axis_coord(self._time_name, self._time_loc) - if isinstance(times, (u.Quantity, TimeDelta)): - try: - if self.meta.date_reference: - times += self.meta.date_reference - except AttributeError: - pass - return times + self._time_name, self._time_loc = _find_axis_name(SUPPORTED_TIME_NAMES, + self.wcs.world_axis_physical_types, + self.extra_coords) + if not self._time_name: + raise ValueError(f"Time {AXIS_NOT_FOUND_ERROR} {SUPPORTED_TIME_NAMES}") + return Time(self._get_axis_coord(self._time_name, self._time_loc)) @property def exposure_time(self): - if not self._exposure_time_name: - raise ValueError("Exposure time" + AXIS_NOT_FOUND_ERROR + - f"{SUPPORTED_EXPOSURE_NAMES}") - return self._get_axis_coord(self._exposure_time_name, self._exposure_time_loc) + exposure_times = None + i = 0 + while exposure_times is None and i < len(SUPPORTED_EXPOSURE_NAMES): + exposure_times = self.meta.get(SUPPORTED_EXPOSURE_NAMES[i], None) + i += 1 + return exposure_times @property def lon(self): + warnings.warn("'.lon' is deprecated. Please use '.celestial'.") if not self._longitude_name: raise ValueError("Longitude" + AXIS_NOT_FOUND_ERROR + f"{SUPPORTED_LONGITUDE_NAMES}") - return self._get_axis_coord(self._longitude_name, self._longitude_loc) + return self._get_axis_coord_values(self._longitude_name, self._longitude_loc) @property def lat(self): + warnings.warn("'.lat' is deprecated. Please use '.celestial'.") if not self._latitude_name: raise ValueError("Latitude" + AXIS_NOT_FOUND_ERROR + f"{SUPPORTED_LATITUDE_NAMES}") - return self._get_axis_coord(self._latitude_name, self._latitude_loc) + return self._get_axis_coord_values(self._latitude_name, self._latitude_loc) + + @property + def celestial(self): + if not self._longitude_name: + self._longitude_name, self._longitude_loc = _find_axis_name( + SUPPORTED_LONGITUDE_NAMES, self.wcs.world_axis_physical_types, self.extra_coords) + if not self._latitude_name: + self._latitude_name, self._latitude_loc = _find_axis_name( + SUPPORTED_LATITUDE_NAMES, self.wcs.world_axis_physical_types, self.extra_coords) + if self._longitude_name: + celestial_name = self._longitude_name + celestial_loc = self._longitude_loc + elif self._latitude_name: + celestial_name = self._latitude_name + celestial_loc = self._latitude_loc + else: + raise ValueError( + f"Celestial {AXIS_NOT_FOUND_ERROR} " + f"{np.concatenate([SUPPORTED_LONGITUDE_NAMES, SUPPORTED_LATITUDE_NAMES])}") + return self._get_axis_coord(celestial_name, celestial_loc) + def apply_exposure_time_correction(self, undo=False, force=False): # Get exposure time in seconds. @@ -337,19 +366,23 @@ def apply_exposure_time_correction(self, undo=False, force=False): new_data, new_uncertainty, new_unit = _calculate_exposure_time_correction( self.data, self.uncertainty, self.unit, exposure_time_s, force=force) # Return new instance of SpectrogramCube with correction applied/undone. - - return self.__class__( - new_data, wcs=self.wcs, - extra_coords=convert_extra_coords_dict_to_input_format(self.extra_coords, - self.missing_axes), - unit=new_unit, uncertainty=new_uncertainty, meta=self.meta, mask=self.mask, - missing_axes=self.missing_axes) + result = copy.deepcopy(self) + result.data = new_data + result.uncertainty = new_uncertainty + result.unit = new_unit + return result def _get_axis_coord(self, axis_name, coord_loc): + if coord_loc == "wcs": + return self.axis_world_coords(axis_name)[0] + elif coord_loc == "extra_coords": + return self.axis_world_coords(wcs=self.extra_coords[axis_name])[0] + + def _get_axis_coord_values(self, axis_name, coord_loc): if coord_loc == "wcs": return self.axis_world_coords_values(axis_name)[0] elif coord_loc == "extra_coords": - return self.extra_coords[axis_name]["value"] + return self.axis_world_coords_values(wcs=self.extra_coords[axis_name])[0] def _find_axis_name(supported_names, world_axis_physical_types, extra_coords): diff --git a/sunraster/spectrogram_sequence.py b/sunraster/spectrogram_sequence.py index 379226b..4d4da06 100644 --- a/sunraster/spectrogram_sequence.py +++ b/sunraster/spectrogram_sequence.py @@ -1,13 +1,16 @@ import numbers import textwrap +import astropy.units as u import numpy as np +from astropy.coordinates import SkyCoord +from astropy.time import Time from ndcube import NDCubeSequence -import astropy.units as u - -from sunraster.spectrogram import SpectrogramABC - +from sunraster.spectrogram import (SpectrogramABC, _find_axis_name, + SUPPORTED_LATITUDE_NAMES, SUPPORTED_LONGITUDE_NAMES, + SUPPORTED_SPECTRAL_NAMES, SUPPORTED_TIME_NAMES) + __all__ = ['SpectrogramSequence', 'RasterSequence'] RASTER_AXIS_NAME = "raster scan" @@ -50,11 +53,16 @@ def spectral_axis(self): @property def time(self): - return np.concatenate([raster.time for raster in self.data]) + return Time(np.concatenate([raster.time for raster in self.data])) @property def exposure_time(self): - return np.concatenate([raster.exposure_time for raster in self.data]) + exposure_type = type(self.data[0].exposure_time) + exposure_time = np.concatenate([raster.exposure_time for raster in self.data]) + try: + return exposure_type(exposure_time) + except Exception: + return exposure_time @property def lon(self): @@ -64,6 +72,14 @@ def lon(self): def lat(self): return u.Quantity([raster.lat for raster in self.data]) + @property + def celestial(self): + sc = SkyCoord([raster.celestial for raster in self.data]) + sc_shape = list(sc.shape) + sc_shape.insert(0, len(self.data)) + sc_shape[1] = int(sc_shape[1] / sc_shape[0]) + return sc.reshape(sc_shape) + def apply_exposure_time_correction(self, undo=False, copy=False, force=False): """ Applies or undoes exposure time correction to data and uncertainty and @@ -111,26 +127,48 @@ def apply_exposure_time_correction(self, undo=False, copy=False, force=False): def __str__(self): data0 = self.data[0] + if not (data0._time_name and data0._longitude_name and data0._latitude_name + and data0._spectral_name): + for i, cube in enumerate(self): + self.data[i]._time_name, self.data[i]._time_loc = _find_axis_name( + SUPPORTED_TIME_NAMES, cube.wcs.world_axis_physical_types, + cube.extra_coords) + self.data[i]._longitude_name, self.data[i]._longitude_loc = _find_axis_name( + SUPPORTED_LONGITUDE_NAMES, cube.wcs.world_axis_physical_types, + cube.extra_coords) + self.data[i]._latitude_name, self.data[i]._latitude_loc = _find_axis_name( + SUPPORTED_LATITUDE_NAMES, cube.wcs.world_axis_physical_types, + cube.extra_coords) + self.data[i]._spectral_name, self.data[i]._spectral_loc = _find_axis_name( + SUPPORTED_SPECTRAL_NAMES, cube.wcs.world_axis_physical_types, + cube.extra_coords) + data0 = self.data[0] if data0._time_name: - start_time = data0.time.value if data0.time.isscalar else data0.time.value.squeeze()[0] + start_time = data0.time if data0.time.isscalar else data0.time.squeeze()[0] data_1 = self.data[-1] - stop_time = data_1.time.value if data_1.time.isscalar else data_1.time.value.squeeze()[-1] - time_period = start_time if start_time == stop_time else (start_time, stop_time) + stop_time = (data_1.time if data_1.time.isscalar else data_1.time.squeeze()[-1]) + time_period = (start_time if start_time == stop_time + else Time([start_time.iso, stop_time.iso])) else: time_period = None - if data0._longitude_name: - lons = self.lon - lon_min = lons.min() - lon_max = lons.max() - lon_range = lon_min if lon_min == lon_max else u.Quantity([lon_min, lon_max]) + if data0._longitude_name or data0._latitude_name: + sc = self.celestial + component_names = dict( + [(item, key) for key, item in sc.representation_component_names.items()]) + lon = getattr(sc, component_names["lon"]) + lat = getattr(sc, component_names["lat"]) + if sc.isscalar: + lon_range = lon + lat_range = lat + else: + lon_range = u.Quantity([lon.min(), lon.max()]) + lat_range = u.Quantity([lat.min(), lat.max()]) + if lon_range[0] == lon_range[1]: + lon_range = lon_range[0] + if lat_range[0] == lat_range[1]: + lat_range = lat_range[0] else: lon_range = None - if data0._latitude_name: - lats = self.lat - lat_min = lats.min() - lat_max = lats.max() - lat_range = lat_min if lat_min == lat_max else u.Quantity([lat_min, lat_max]) - else: lat_range = None if data0._spectral_name: spectral_vals = self.spectral_axis @@ -188,11 +226,23 @@ def __init__(self, data_list, common_axis, meta=None): SnS_dimensions = SpectrogramSequence.cube_like_dimensions raster_array_axis_physical_types = SpectrogramSequence.array_axis_physical_types SnS_array_axis_physical_types = SpectrogramSequence.cube_like_array_axis_physical_types - raster_axis_extra_coords = SpectrogramSequence.sequence_axis_extra_coords - SnS_axis_extra_coords = SpectrogramSequence.common_axis_extra_coords + raster_axis_coords = SpectrogramSequence.sequence_axis_coords + SnS_axis_coords = SpectrogramSequence.common_axis_coords plot_as_raster = SpectrogramSequence.plot plot_as_SnS = SpectrogramSequence.plot_as_cube + @property + def raster_axis_extra_coords(self): + warnings.warn( + "'.raster_axis_extra_coords is deprecated. Please use '.raster_axis_coords'.") + return self.raster_axis_coords + + @property + def SnS_axis_extra_coords(self): + warnings.warn( + "'.SnS_axis_extra_coords is deprecated. Please use '.SnS_axis_coords'.") + return self.SnS_axis_coords + def _set_single_scan_instrument_axes_types(self): if len(self.data) < 1: self._single_scan_instrument_axes_types = np.empty((0,), dtype=object) @@ -203,8 +253,16 @@ def _set_single_scan_instrument_axes_types(self): self._single_scan_instrument_axes_types[self._common_axis] = \ self._slit_step_axis_name # Spectral axis name. - spectral_raster_index = [physical_type == (self.data[0]._spectral_name,) - for physical_type in self.data[0].array_axis_physical_types] + # If spectral name not present in raster cube, try finding it. + if not self.data[0]._spectral_name: + for i, cube in enumerate(self): + self.data[i]._spectral_name, self.data[i]._spectral_name = _find_axis_name( + SUPPORTED_SPECTRAL_NAMES, cube.wcs.world_axis_physical_types, + cube.extra_coords) + spectral_name = self.data[0]._spectral_name + array_axis_physical_types = self.data[0].array_axis_physical_types + spectral_raster_index = [physical_type == (spectral_name,) + for physical_type in array_axis_physical_types] spectral_raster_index = np.arange(self.data[0].data.ndim)[spectral_raster_index] if len(spectral_raster_index) == 1: self._single_scan_instrument_axes_types[spectral_raster_index] = \