From 4f7aea4837c2e6f3377bc04e67dbf4314882b7b9 Mon Sep 17 00:00:00 2001 From: DanRyanIrish Date: Wed, 25 Aug 2021 20:41:27 +0100 Subject: [PATCH 1/8] Define IRIS Meta object. --- sunraster/instr/iris.py | 204 ++++++++++++++++++++++++++++++++++++---- 1 file changed, 188 insertions(+), 16 deletions(-) diff --git a/sunraster/instr/iris.py b/sunraster/instr/iris.py index e3c215f7..4cc53e5f 100644 --- a/sunraster/instr/iris.py +++ b/sunraster/instr/iris.py @@ -1,14 +1,16 @@ import copy -from astropy.wcs import WCS -import numpy as np -from ndcube import NDCollection - 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'] @@ -136,15 +138,15 @@ def read_iris_spectrograph_level2_fits( # 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)] + ("raster position", 0, u.Quantity(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)] + 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): # Determine values of properties dependent on detector type. if "FUV" in hdulist[0].header["TDET{0}".format(window_fits_indices[i])]: @@ -206,10 +208,12 @@ 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=single_file_meta, mask=data_mask) + for ec in window_extra_coords: + cube.extra_coords.add(*ec) + 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, @@ -217,3 +221,171 @@ def read_iris_spectrograph_level2_fits( for window_name in spectral_windows_req] # Initialize an NDCollection object. return NDCollection(window_data_pairs, aligned_axes=(0, 1, 2), meta=top_meta) + + +class IRISSGMeta(Meta, metaclass=SlitSpectrographMetaABC): + def __init__(self, spectral_window, **kwargs): + spectral_windows = np.array([hdulist[0].header["TDESC{0}".format(i)] + for i in range(1, hdulist[0].header["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_windows.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} + Date:\t\t\t{self.date_reference} + OBS ID:\t\t\t{self.observing_mode_id} + """) + + 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_window_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 FOV_width_y(self): + """Width of the field of view in the Y (slit) direction.""" + return self.get("FOVY") * u.arcsec + + @property + def FOV_width_x(self): + """Width of the field of view 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: + return self.get("SUMSPTRF") + else: + return self.get("SUMSPTRN") From f721487f2bcaa9454da736616a27137066494704 Mon Sep 17 00:00:00 2001 From: DanRyanIrish Date: Thu, 26 Aug 2021 11:07:34 +0100 Subject: [PATCH 2/8] In IRIS SG reader, move all extra coords to metadata except time. Also move all metadata to the raster cube level. --- sunraster/instr/iris.py | 142 +++++++++------------------------------- 1 file changed, 31 insertions(+), 111 deletions(-) diff --git a/sunraster/instr/iris.py b/sunraster/instr/iris.py index 4cc53e5f..8886266c 100644 --- a/sunraster/instr/iris.py +++ b/sunraster/instr/iris.py @@ -1,4 +1,5 @@ import copy +import textwrap import astropy.units as u import numpy as np @@ -54,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)]) @@ -78,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, u.Quantity(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. @@ -169,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( @@ -210,29 +129,28 @@ def read_iris_spectrograph_level2_fits( # Appending NDCube instance to the corresponding window key in dictionary's list. cube = SpectrogramCube(hdulist[window_fits_indices[i]].data, wcs=wcs_, uncertainty=out_uncertainty, unit=DN_unit, - meta=single_file_meta, mask=data_mask) - for ec in window_extra_coords: - cube.extra_coords.add(*ec) + 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, spectral_window, **kwargs): - spectral_windows = np.array([hdulist[0].header["TDESC{0}".format(i)] - for i in range(1, hdulist[0].header["NWIN"] + 1)]) + 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_windows.sum() > 1: + elif window_mask.sum() > 1: raise ValueError( "Spectral window must be unique. " f"Input spectral window: {spectral_window}; " @@ -247,8 +165,10 @@ def __str__(self): 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): @@ -344,18 +264,18 @@ def number_raster_positions(self): self.get("NRASTERP") @property - def spectral_window_range(self): + 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 FOV_width_y(self): - """Width of the field of view in the Y (slit) direction.""" + 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 FOV_width_x(self): - """Width of the field of view in the X (rastering) direction.""" + 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 @@ -385,7 +305,7 @@ def spatial_summing_factor(self): @property def spectral_summing_factor(self): """Number of pixels summed together in the spectral direction.""" - if "fuv" in self.detector: + if "fuv" in self.detector.lower(): return self.get("SUMSPTRF") else: return self.get("SUMSPTRN") From 53d0b70075625e797f800fae8b331dc81d085d74 Mon Sep 17 00:00:00 2001 From: DanRyanIrish Date: Thu, 26 Aug 2021 11:11:26 +0100 Subject: [PATCH 3/8] Remove extra_coords kwarg from SpectrogramCube init. This is in accordance with API changes in ndcube 2.0. Create new celestial properties in SpectrogramCube and SpectrogramSequence to give the celestial coordinates of the pixels. --- sunraster/spectrogram.py | 106 ++++++++++++++++-------------- sunraster/spectrogram_sequence.py | 54 +++++++++++---- 2 files changed, 97 insertions(+), 63 deletions(-) diff --git a/sunraster/spectrogram.py b/sunraster/spectrogram.py index 9b2397d7..9f3480d1 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'] @@ -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 @@ -199,19 +205,20 @@ def __str__(self): 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 + if self._longitude_name or self._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([self.lon.min(), self.lon.max()]) + lon_range = u.Quantity([lon.min(), lon.max()]) + lat_range = u.Quantity([lat.min(), lat.max()]) else: lon_range = None - if self._latitude_name: - if self.lat.isscalar: - lat_range = self.lat - else: - lat_range = u.Quantity([self.lat.min(), self.lat.max()]) - else: lat_range = None if self._spectral_name: if self.spectral_axis.isscalar: @@ -237,15 +244,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,11 +263,8 @@ 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): @@ -283,35 +278,44 @@ 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 + 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) + return self._get_axis_coord_values(self._exposure_time_name, self._exposure_time_loc) @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 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("Celestial" + AXIS_NOT_FOUND_ERROR + + f"{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 +341,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 379226b4..5e1fdde7 100644 --- a/sunraster/spectrogram_sequence.py +++ b/sunraster/spectrogram_sequence.py @@ -1,10 +1,11 @@ import numbers import textwrap +import astropy.units as u import numpy as np +from astropy.coordinates import SkyCoord from ndcube import NDCubeSequence -import astropy.units as u from sunraster.spectrogram import SpectrogramABC @@ -64,6 +65,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 @@ -118,19 +127,24 @@ def __str__(self): time_period = start_time if start_time == stop_time else (start_time, stop_time) 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 +202,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) From b4e59a0a1a2030b30eca46c9e28ddc8bed5bf1c8 Mon Sep 17 00:00:00 2001 From: DanRyanIrish Date: Thu, 26 Aug 2021 11:27:18 +0100 Subject: [PATCH 4/8] Move calculation of coord names and storage locations from Spectrogram init to the properties which rely on them. --- sunraster/spectrogram.py | 47 ++++++++++++++++++++++--------- sunraster/spectrogram_sequence.py | 43 ++++++++++++++++++++++------ 2 files changed, 69 insertions(+), 21 deletions(-) diff --git a/sunraster/spectrogram.py b/sunraster/spectrogram.py index 9f3480d1..44ce4578 100644 --- a/sunraster/spectrogram.py +++ b/sunraster/spectrogram.py @@ -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", @@ -198,14 +198,17 @@ def __init__(self, data, wcs, unit=None, uncertainty=None, meta=None, 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 or self._latitude_name: + except ValueError as err: + if AXIS_NOT_FOUND_ERROR in err.args[0]: + time_period = None + else: + raise err + try: sc = self.celestial component_names = dict( [(item, key) for key, item in sc.representation_component_names.items()]) @@ -217,16 +220,22 @@ def __str__(self): else: lon_range = u.Quantity([lon.min(), lon.max()]) lat_range = u.Quantity([lat.min(), lat.max()]) - else: - lon_range = None - lat_range = None - if self._spectral_name: + 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__))} @@ -268,6 +277,9 @@ def __getitem__(self, item): @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}") @@ -276,8 +288,11 @@ def spectral_axis(self): @property def time(self): if not self._time_name: - raise ValueError("Time" + AXIS_NOT_FOUND_ERROR + - f"{SUPPORTED_TIME_NAMES}") + 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 @@ -305,6 +320,12 @@ def lat(self): @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 diff --git a/sunraster/spectrogram_sequence.py b/sunraster/spectrogram_sequence.py index 5e1fdde7..5e5d3e08 100644 --- a/sunraster/spectrogram_sequence.py +++ b/sunraster/spectrogram_sequence.py @@ -4,11 +4,13 @@ import astropy.units as u import numpy as np from astropy.coordinates import SkyCoord +from astropy.time import Time from ndcube import NDCubeSequence - -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" @@ -120,11 +122,28 @@ 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 or data0._latitude_name: @@ -229,8 +248,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] = \ From 9fec766832ac0fcd6e1a107233e26f67fdc8b9e8 Mon Sep 17 00:00:00 2001 From: DanRyanIrish Date: Thu, 26 Aug 2021 12:43:24 +0100 Subject: [PATCH 5/8] Re-implement SpectrogramCube.exposure_time assuming exposure time stored in meta. Also tidy up/bugfix repr strings. --- sunraster/spectrogram.py | 13 ++++++++----- sunraster/spectrogram_sequence.py | 9 +++++++-- 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/sunraster/spectrogram.py b/sunraster/spectrogram.py index 44ce4578..f4221161 100644 --- a/sunraster/spectrogram.py +++ b/sunraster/spectrogram.py @@ -202,7 +202,8 @@ def __str__(self): if self.time.isscalar: time_period = self.time else: - time_period = (str(self.time.min()), str(self.time.max())) + 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 @@ -297,10 +298,12 @@ def time(self): @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_values(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): diff --git a/sunraster/spectrogram_sequence.py b/sunraster/spectrogram_sequence.py index 5e5d3e08..4d4da06b 100644 --- a/sunraster/spectrogram_sequence.py +++ b/sunraster/spectrogram_sequence.py @@ -53,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): From f2d14a8c71938c7468daf2a7bf2cf3c4e38b899e Mon Sep 17 00:00:00 2001 From: DanRyanIrish Date: Thu, 26 Aug 2021 13:54:07 +0100 Subject: [PATCH 6/8] Change SPICEMeta to inherit from ndcube Meta. --- sunraster/instr/spice.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sunraster/instr/spice.py b/sunraster/instr/spice.py index 6dd745b9..ea42aa70 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 From 57d521709d26732c1965a14196635fdb18b20ca4 Mon Sep 17 00:00:00 2001 From: DanRyanIrish Date: Thu, 26 Aug 2021 13:55:38 +0100 Subject: [PATCH 7/8] Add some changelogs for updates required by ndcube 2.0. --- changelog/ndc2.1.breaking.rst | 1 + changelog/ndc2.1.feature.rst | 1 + changelog/ndc2.2.breaking.rst | 1 + changelog/ndc2.2.feature.rst | 1 + changelog/ndc2.3.breaking.rst | 1 + changelog/ndc2.3.feature.rst | 1 + 6 files changed, 6 insertions(+) create mode 100644 changelog/ndc2.1.breaking.rst create mode 100644 changelog/ndc2.1.feature.rst create mode 100644 changelog/ndc2.2.breaking.rst create mode 100644 changelog/ndc2.2.feature.rst create mode 100644 changelog/ndc2.3.breaking.rst create mode 100644 changelog/ndc2.3.feature.rst diff --git a/changelog/ndc2.1.breaking.rst b/changelog/ndc2.1.breaking.rst new file mode 100644 index 00000000..c455c282 --- /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 00000000..47612533 --- /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 00000000..95569698 --- /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 00000000..373249cc --- /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 00000000..53a79510 --- /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 00000000..3c25907f --- /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. From f32591dd8678b6d9a319c6383325b84b06aca017 Mon Sep 17 00:00:00 2001 From: DanRyanIrish Date: Fri, 27 Aug 2021 16:57:23 +0100 Subject: [PATCH 8/8] Fix bug in SpectrogramCube.__str__ --- sunraster/spectrogram.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/sunraster/spectrogram.py b/sunraster/spectrogram.py index f4221161..f75b34bb 100644 --- a/sunraster/spectrogram.py +++ b/sunraster/spectrogram.py @@ -336,8 +336,9 @@ def celestial(self): celestial_name = self._latitude_name celestial_loc = self._latitude_loc else: - raise ValueError("Celestial" + AXIS_NOT_FOUND_ERROR + - f"{SUPPORTED_LONGITUDE_NAMES + SUPPORTED_LATITUDE_NAMES}") + 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)