Merge pull request #178 from DanRyanIrish/spice

Enable SPICE FITS reader to handle Dumbbell windows and multiple files

changelog/178.bugfix.1.rst

@@ -0,0 +1 @@
+Allow SPICE FITS reader to read handle dumbbell windows.

changelog/178.bugfix.2.rst

@@ -0,0 +1 @@
+Stop `~sunraster.spectrogram_sequence.SpectrogramSequence` crashing when time coord not 1-D.

changelog/178.feature.rst

@@ -0,0 +1 @@
+Enable SPICE FITS reader to handle multiple files.

changelog/178.trivial.rst

@@ -0,0 +1 @@
+Altered names of some SPICEMeta properties.

sunraster/instr/spice.py

@@ -1,5 +1,6 @@
-import textwrap
 import copy
+import numbers
+import textwrap
 
 import numpy as np
 from astropy.io import fits
@@ -11,39 +12,174 @@
 from ndcube import NDCollection
 
 from sunraster.meta import Meta, SlitSpectrographMetaABC
-from sunraster import SpectrogramCube
+from sunraster import SpectrogramCube, SpectrogramSequence, RasterSequence
 
 
 __all__ = ["read_spice_l2_fits", "SPICEMeta"]
 
 
-def read_spice_l2_fits(filename, windows=None, memmap=True):
+INCORRECT_OBSID_MESSAGE = "File has incorrect SPIOBSID."
+
+
+def read_spice_l2_fits(filenames, windows=None, memmap=True, read_dumbbells=False):
     """Read SPICE level 2 FITS file.
 
+    Parameters
+    ----------
+    filenames: iterable of `str`
+        The name(s), including path, of the SPICE FITS file(s) to read.
+
+    windows: iterable of `str`
+        The names of the windows to read.
+        All windows must of the same type: dumbbell and regular.
+        Default=None implies all narrow-slit or dumbbell windows read out
+        depending on value of read_dumbells kwarg.  See below.
+
+    memmap: `bool`
+        If True, FITS file is reading with memory mapping.
+
+    read_dumbbells: `bool`
+        Defines whether dumbbell or regular windows are returned.
+        If True, returns the dumbbell windows.
+        If False, returns regular windows.
+        Default=False
+        Ignored if windows kwarg is set.
+
+    Returns
+    -------
+    output: `ndcube.NDCollection` or `sunraster.SpectrogramCube`, `sunraster.RasterSequence`,
+        `sunraster.SpectrogramSequence`
+        A collection of spectrogram/raster cubes/sequences, one for each window.
+        If only one window present or requested, a single spectrogram cube
+        or sequence is returned.
+    """
+    # Sanitize inputs.
+    if isinstance(filenames, str):
+        filenames = [filenames]
+    # Read first file.
+    first_cubes = _read_single_spice_l2_fits(filenames[0], windows=windows, memmap=memmap,
+                                             read_dumbbells=read_dumbbells)
+    # Derive information for consistency checks between files and read subsequent files.
+    if len(filenames) > 1:
+        # Wrap windows from first file in lists
+        # so windows from other files can be appended.
+        cube_lists = dict([(key, [value]) for key, value in first_cubes.items()])
+        # Get info from first file for consistency checks between files.
+        first_meta = _get_meta_from_last_added(cube_lists)
+        first_obs_id = _get_obsid(first_meta)
+        if windows is None:
+            windows = list(cube_lists.keys())
+        # Read subsequent files and append output to relevant window in cube_lists.
+        for i, filename in enumerate(filenames[1:]):
+            try:
+                cube_lists = _read_single_spice_l2_fits(filename, windows=windows, memmap=memmap,
+                                                        read_dumbbells=read_dumbbells,
+                                                        output=cube_lists, spice_id=first_obs_id)
+            except ValueError as err:
+                err_message = err.args[0]
+                if INCORRECT_OBSID_MESSAGE in err_message:
+                    this_obs_id = err_message.split()[-1]
+                    raise ValueError(
+                        "All files must correspond to same observing campaign/SPICE OBS ID. "
+                         f"First file SPICE OBS ID: {first_obs_id}; "
+                         f"{i+1}th file SPICE OBS ID: {this_obs_id}")
+        # Depending on type of file, combine data from different files into
+        # SpectrogramSequences and RasterSequences.
+        is_raster = "ras" in first_meta.get("FILENAME") and \
+            not any([window[1].meta.contains_dumbbell for window in cube_lists.values()])
+        if is_raster:
+            sequence_class = RasterSequence
+        else:
+            sequence_class = SpectrogramSequence
+        window_sequences = [(key, sequence_class([v[0] for v in value], common_axis=-1))
+                            for key, value in cube_lists.items()]
+    else:
+        # If only one file being read, leave data in SpectrogramCube objects.
+        window_sequences = list(first_cubes.items())
+    if len(window_sequences) > 1:
+        # Data should be aligned along all axes except the spectral axis.
+        # But they should be aligned along all axes if they come from the
+        # same spectral window, e.g. because they are dumbbell windows.
+        first_sequence = window_sequences[0][1]
+        first_spectral_window = first_sequence[0].meta.spectral_window
+        if all([window[1][0].meta.spectral_window == first_spectral_window
+                for window in window_sequences]):
+            aligned_axes = tuple(range(len(first_sequence.dimensions)))
+        else:
+            aligned_axes = np.where(
+                np.asarray(first_sequence.world_axis_physical_types) != "em.wl")[0]
+            aligned_axes = tuple([int(i) for i in aligned_axes])
+    else:
+        aligned_axes = None
+    return NDCollection(window_sequences, aligned_axes=aligned_axes)
+
+
+def _get_meta_from_last_added(obj):
+    return list(obj.values())[0][-1].meta
+
+
+def _get_obsid(spice_meta):
+    return spice_meta.spice_observation_id
+
+
+def _read_single_spice_l2_fits(filename, windows=None, memmap=True, read_dumbbells=False,
+                               output=None, spice_id=None):
+    """Read SPICE level 2 FITS file(s).
+
     Parameters
     ----------
     filename: `str`
         The name, including path, of the SPICE FITS file to read.
 
     windows: iterable of `str`
         The names of the windows to read.
-        Default=None implies all windows read out.
+        All windows must of the same type: dumbbell and regular.
+        Default=None implies all narrow-slit or dumbbell windows read out
+        depending on value of read_dumbells kwarg.  See below.
 
     memmap: `bool`
         If True, FITS file is reading with memory mapping.
 
+    read_dumbbells: `bool`
+        Defines whether dumbbell or regular windows are returned.
+        If True, returns the dumbbell windows.
+        If False, returns regular windows.
+        Default=False
+        Ignored if windows kwarg is set.
+
+    output: `dict` of `list`s (optional)
+        A dictionary of lists with the same keys are the windows kwarg.
+        The output for each window will be appended to the list corresponding
+        the window's name.
+
+    spice_id: `int` (optional)
+        If not None, file must have a SPIOBSID equal to this value.
+        Otherwise an error is raised
+
     Returns
     -------
-    output: `ndcube.NDCollection` or `sunraster.SpectrogramCube`
+    output: `dict` of `sunraster.SpectrogramCube`
         A collection of spectrogram cubes, one for each window.
-        If only one window present or requested, a single spectrogram cube is returned.
     """
     window_cubes = []
+    dumbbell_label = "DUMBBELL"
     with fits.open(filename, memmap=memmap) as hdulist:
-        # Retrieve window names from FITS file.
+        if isinstance(spice_id, numbers.Integral) and hdulist[0].header["SPIOBSID"] != spice_id:
+            raise ValueError(f"{INCORRECT_OBSID_MESSAGE}  "
+                             f"Expected {spice_id}.  Got {hdulist[0].header['SPIOBSID']}.")
+        # Derive names of windows to be read.
         if windows is None:
-            windows = [hdu.header["EXTNAME"] for hdu in hdulist
-                       if hdu.header["EXTNAME"] != "VARIABLE_KEYWORDS"]
+            if read_dumbbells:
+                windows = [hdu.header["EXTNAME"] for hdu in hdulist
+                           if dumbbell_label in hdu.header["EXTNAME"]]
+            else:
+                windows = [hdu.header["EXTNAME"] for hdu in hdulist
+                           if (hdu.header["EXTNAME"] != "VARIABLE_KEYWORDS" and
+                               dumbbell_label not in hdu.header["EXTNAME"])]
+        dumbbells_requested = [dumbbell_label in window for window in windows]
+        if any(dumbbells_requested) and not all(dumbbells_requested):
+            raise ValueError("Cannot read dumbbell and other window types simultaneously.")
+        # Retrieve window names from FITS file.
         for i, hdu in enumerate(hdulist):
             if hdu.header["EXTNAME"] in windows:
                 # Define metadata object.
@@ -63,16 +199,15 @@ def read_spice_l2_fits(filename, windows=None, memmap=True):
                     data=data, wcs=wcs, mask=np.isnan(data), unit=u.adu,
                     extra_coords=[("exposure time", -1, exp_times)], meta=meta,
                     instrument_axes=("raster scan", "spectral", "slit", "slit step"))
-                window_cubes.append((meta.get("EXTNAME"), spectrogram))
-
-    if len(windows) > 1:
-        aligned_axes = np.where(np.asarray(spectrogram.world_axis_physical_types) != "em.wl")[0]
-        aligned_axes = tuple([int(i) for i in aligned_axes])
-        output = NDCollection(window_cubes, aligned_axes=aligned_axes)
+                window_name = meta.get("EXTNAME")
+                if output is None:
+                    window_cubes.append((window_name, spectrogram))
+                else:
+                    output[window_name].append(spectrogram)
+    if output is None:
+        return dict(window_cubes)
     else:
-        output = spectrogram
-
-    return output
+        return output
 
 
 def _convert_fits_comments_to_key_value_pairs(fits_header):
@@ -124,11 +259,20 @@ def __repr__(self):
     @property
     def spectral_window(self):
         spectral_window = self.get("EXTNAME")
+        # Remove redundant text associated with dumbbells.
+        joiner = "_"
+        if self.contains_dumbbell:
+            dummy_txt = ""
+            spectral_window = spectral_window.replace("DUMBBELL", dummy_txt)
+            spectral_window = spectral_window.replace("UPPER", dummy_txt)
+            spectral_window = spectral_window.replace("LOWER", dummy_txt)
+            spectral_window = joiner.join(list(filter((dummy_txt).__ne__,
+                                                      spectral_window.split(joiner))))
+        # Remove other redundant text from window name.
         redundant_txt = "WINDOW"
         if redundant_txt in spectral_window:
-            spectral_window = np.asanyarray(spectral_window.split("_"))
-            idx = np.array([redundant_txt not in window_chunk for window_chunk in spectral_window])
-            spectral_window = "_".join(spectral_window[idx])
+            spectral_window = joiner.join([comp for comp in spectral_window.split(joiner)
+                                           if "WINDOW" not in comp])
         return spectral_window
 
     @property
@@ -156,11 +300,11 @@ def rsun_angular(self):
         return self._construct_quantity("RSUN_ARC")
 
     @property
-    def observing_mode_id(self):
+    def spice_observation_id(self):
         return self.get("SPIOBSID")
 
     @property
-    def observatory_radial_velocity(self):
+    def observer_radial_velocity(self):
         return self._construct_quantity("OBS_VR")
 
     @property
@@ -180,7 +324,7 @@ def date_end(self):
         return self._construct_time("DATE-END")
 
     @property
-    def observer_coordinate(self):
+    def observer_location(self):
         lon_unit = u.deg
         lat_unit = u.deg
         radius_unit = u.m
@@ -214,10 +358,6 @@ def bias_frame_subtracted_onboard(self):
     def window_type(self):
         return self.get("WIN_TYPE")
 
-    @property
-    def window_table_id(self):
-        return self.get("WINTABID")
-
     @property
     def slit_id(self):
         return self.get("SLIT_ID")
@@ -227,10 +367,14 @@ def slit_width(self):
         return self._construct_quantity("SLIT_WID")
 
     @property
-    def dumbbell(self):
-        dumbell_types = ["none", "lower", "upper"]
-        dumbell_idx = self.get("DUMBBELL")
-        return dumbell_types[dumbell_idx]
+    def contains_dumbbell(self):
+        return self.get("DUMBBELL") in [1, 2]
+
+    @property
+    def dumbbell_type(self):
+        dumbbell_types = [None, "lower", "upper"]
+        dumbbell_idx = self.get("DUMBBELL")
+        return dumbbell_types[dumbbell_idx]
 
     @property
     def solar_B0(self):

sunraster/meta.py

@@ -24,7 +24,7 @@ def processing_level(self):
         pass
 
     @abc.abstractproperty
-    def observer_coordinate(self):
+    def observer_location(self):
         """Coordinate of observatory location based on header info."""
         pass
 
@@ -78,7 +78,7 @@ def observing_mode_id(self):
         pass
 
     @abc.abstractproperty
-    def observatory_radial_velocity(self):
+    def observer_radial_velocity(self):
         """Velocity of observatory in direction of source."""
         pass
 

sunraster/spectrogram_sequence.py

@@ -112,9 +112,9 @@ def apply_exposure_time_correction(self, undo=False, copy=False, force=False):
     def __str__(self):
         data0 = self.data[0]
         if data0._time_name:
-            start_time = data0.time.value if data0.time.isscalar else data0.time[0].value
+            start_time = data0.time.value if data0.time.isscalar else data0.time.value.squeeze()[0]
             data_1 = self.data[-1]
-            stop_time = data_1.time.value if data_1.time.isscalar else data_1.time[-1].value
+            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)
         else:
             time_period = None

sunraster/tests/__init__.py

@@ -0,0 +1,7 @@
+import os.path
+
+import sunraster
+
+__all__ = ["test_data_path"]
+
+test_data_dir = os.path.join(os.path.dirname(sunraster.__file__), "tests", "data")