Merge pull request #43 from solo-spice/fov

Field-of-view plotting

README.rst

@@ -51,12 +51,18 @@ Documentation for this package is available on `Read the Docs <https://sospice.r
 -  Instrument modelling: ``instrument_modelling``
 
    -  ``Spice``: instrument calibration parameters, effective area,
-      quantum efficiency…
+      quantum efficiency...
    -  ``Study``: study parameters.
    -  ``Observation``: a SPICE observation with some study (including
       low-level functions used to compute the uncertainties on the
       data).
 
+- Other utilities: ``util``
+
+   - ``sigma_clipping``: sigma clipping (for cosmic rays removal).
+   - ``fov``: plot SPICE field-of-views on a background map.
+
+
 Package philosophy
 ------------------
 

docs/source/conf.py

@@ -4,20 +4,16 @@
 # list see the documentation:
 # https://www.sphinx-doc.org/en/master/usage/configuration.html
 
-# -- Path setup --------------------------------------------------------------
-
-# If extensions (or modules to document with autodoc) are in another directory,
-# add these directories to sys.path here. If the directory is relative to the
-# documentation root, use os.path.abspath to make it absolute, like shown here.
-#
-# import os
-# import sys
-# sys.path.insert(0, os.path.abspath('.'))
+import os
+import sys
 from datetime import datetime
 
 # from sospice import __version__
 from setuptools_scm import get_version
 
+# -- Path setup --------------------------------------------------------------
+
+sys.path.insert(0, os.path.abspath("../.."))
 
 # -- Project information -----------------------------------------------------
 

docs/source/utilities.rst

@@ -7,3 +7,35 @@ Local sigma clipping
 This provides a version of
 `astropy.stats.sigma_clip <https://docs.astropy.org/en/stable/api/astropy.stats.sigma_clip.html>`__
 working on a local neighbourhood.
+
+Plotting fields-of-view over a background map
+----------------------------------------------------
+
+Once observations are selected from a ``Catalog``, their Fields-Of-View (FOVs) can be plotted using ``plot_fovs_with_background()``. Different background maps can be selected: maps with some specific data (e.g. a HMI synoptic map or Solar Orbiter/EUI/FSI), a blank map with some projection (in development), or any map already plotted by the user.
+
+After
+
+.. code:: python
+
+    from sospice import Catalog, plot_fovs_with_background
+    cat = Catalog(release_tag="4.0")
+
+one can select for example all files for which ``DATE-BEG`` is on a given day
+
+.. code:: python
+
+    observations = cat.find_files(date_min="2022-03-08", date_max="2022-03-09", level="L2")
+
+and then plot them either with a HMI synoptic map as background
+
+.. code:: python
+
+   plot_fovs_with_background(observations, "HMI_synoptic")
+
+or with EUI/FSI data
+
+.. code:: python
+
+    plot_fovs_with_background(observations, "EUI/FSI")
+
+``plot_fovs_with_background()`` uses the lower-level methods of the ``FileMetadata`` and ``Catalog`` classes to compute and plot the FOVs. Advanced users can produce custom FOV plots using these methods and their options.

pyproject.toml

@@ -63,6 +63,10 @@ omit = [
     "__init__.py",
 ]
 
+[tool.pytest.ini_options]
+minversion = "6.0"
+addopts = "--ignore=local"
+
 [tool.towncrier]
 package = "sospice"
 filename = "CHANGELOG.rst"

sospice/__init__.py

@@ -4,3 +4,4 @@
 from .catalog.file_metadata import FileMetadata
 from .calibrate.uncertainties import spice_error
 from .util.sigma_clipping import sigma_clip
+from .util.fov import plot_fov_background, plot_fovs_with_background

sospice/calibrate/uncertainties.py

@@ -37,7 +37,6 @@ def spice_error(hdu=None, data=None, header=None, verbose=True):
     if header["LEVEL"] != "L2":
         raise RuntimeError("Level should be L2")
     data *= u.Unit(header["BUNIT"])
-    print(data.unit)
     study = Study()
     study.init_from_header(header)
     if verbose:

sospice/catalog/catalog.py

@@ -1,11 +1,15 @@
 from dataclasses import dataclass
+from pathlib import Path
+from itertools import cycle
 
+import matplotlib.colors as mcolors
 import pandas as pd
-from pathlib import Path
+import numpy as np
+
 from astropy.utils.data import download_file
 
 from .release import Release
-from .file_metadata import required_columns
+from .file_metadata import FileMetadata, required_columns
 
 
 @dataclass
@@ -160,7 +164,7 @@ def find_files_by_keywords(self, **kwargs):
 
     def find_files_by_date_range(self, date_min=None, date_max=None):
         """
-        Find files in some date range.
+        Find files with DATE-BEG in some date range.
 
         Parameters
         ----------
@@ -260,3 +264,133 @@ def find_files(
                 .T
             )
         return df
+
+    def mid_time(self, method=None):
+        """
+        Find "middle time" for observations in catalog
+
+        Parameters
+        ----------
+        method: str
+            Method for determining middle time. Can be
+
+            * "midrange" (default): middle of time range, from beginning of first observation to end of last observation
+            * "mean": mean of observation times (not weighted by observations durations)
+
+        """
+        if method is None or method == "midrange":
+            begin_min = self["DATE-BEG"].min()
+            begin_max = self["DATE-BEG"].max()
+            last_telapse = self[self["DATE-BEG"] == begin_max].TELAPSE.max()
+            end_max = begin_max + pd.Timedelta(seconds=last_telapse)
+            return begin_min + (end_max - begin_min) / 2
+        elif method == "mean":
+            begin_mean = self["DATE-BEG"].mean()
+            telapse_mean = pd.Timedelta(seconds=self.TELAPSE.mean())
+            return begin_mean + telapse_mean
+        elif method == "barycenter":
+            mid_observation = self["DATE-BEG"] + self.apply(
+                lambda row: pd.Timedelta(seconds=row.TELAPSE / 2), axis=1
+            )
+            weight = self.TELAPSE
+            t0 = mid_observation.iloc[0]
+            return t0 + ((mid_observation - t0) * weight).sum() / weight.sum()
+        else:
+            raise RuntimeError("Invalid method")
+
+    @classmethod
+    def _format_time_range(cls, row, timespec="minutes"):
+        """
+        Format time range for observation
+
+        Parameters
+        ----------
+        row: pd.Series
+            Catalog row
+        timespec: str
+            Time terms specification for pandas.Timestamp.isoformat()
+
+        Return
+        ------
+        str
+            Formatted time range
+
+        The end of the time range is known from the single observation in `row`
+        thanks to an additional element `last_DATE-BEG' in the Series.
+        All dates are DATE-BEG, we don't compute a DATE-END.
+        """
+        t = [row["DATE-BEG"]]
+        is_range = ("last_DATE-BEG" in row.index) and (row["last_DATE-BEG"] != t[0])
+        if is_range:
+            t.append(row["last_DATE-BEG"])
+        t_str = [tt.isoformat(timespec=timespec) for tt in t]
+        if is_range and t_str[0][:10] == t_str[1][:10]:
+            t_str[1] = t_str[1][10:]
+        return " - ".join(t_str)
+
+    def plot_fov(self, ax, **kwargs):
+        """
+        Plot SPICE FOVs on a background map
+
+        Parameters
+        ----------
+        ax: matplotlib.axes.Axes
+            Axes (with relevant projection)
+        color: str or list
+            Color(s) cycle for drawing the FOVs
+        kwargs: dict
+            Keyword arguments, passed to FileMetadata.plot_fov()
+        """
+        time_range_length = self["DATE-BEG"].max() - self["DATE-BEG"].min()
+        if time_range_length > pd.Timedelta(days=60):
+            print(
+                f"Time range length is {time_range_length}, this is long, and probably not what you want; aborting"
+            )
+            return
+        merge_by_spiobsid = True
+        if merge_by_spiobsid:
+            groups = self.groupby("SPIOBSID")
+            fovs = groups.first()
+            fovs_last = groups.last()
+            fovs["last_DATE-BEG"] = fovs_last["DATE-BEG"]
+            fovs.reset_index(inplace=True)
+        else:
+            fovs = Catalog(data_frame=self[list(required_columns)])
+        # label at the position of the plot
+        fovs["fov_text"] = fovs.apply(Catalog._format_time_range, axis=1)
+        # label at the level of the plot (will be de-duplicated afterwards)
+        fovs["fov_label"] = fovs.apply(
+            lambda row: f"{row.STUDY} ({row.MISOSTUD})", axis=1
+        )
+        # color(s)
+        color = kwargs.pop("color", None)
+        studies = sorted(list(self.STUDY.unique()))
+        colors = (
+            mcolors.TABLEAU_COLORS
+            if color is None
+            else color
+            if type(color) is list
+            else [color]
+        )
+        study_color = dict(zip(studies, cycle(colors)))
+        fovs["fov_color"] = fovs.apply(lambda row: study_color[row.STUDY], axis=1)
+        fovs.apply(
+            lambda row: FileMetadata(row).plot_fov(ax, **kwargs),
+            axis=1,
+        )
+        if merge_by_spiobsid:
+            # also plot last FOV, with dashes
+            fovs_last.reset_index(inplace=True)
+            fovs_last["fov_color"] = fovs.fov_color
+            fovs_last["fov_linestyle"] = ":"
+            fovs_last = fovs_last[fovs_last.RASTERNO != 0]
+            fovs_last.apply(
+                lambda row: FileMetadata(row).plot_fov(ax, **kwargs),
+                axis=1,
+            )
+        # De-duplicate labels for legend (an alternative would be
+        # to provide labels only to the first instance of each study)
+        handles, labels = ax.get_legend_handles_labels()
+        unique_indices = [labels.index(x) for x in sorted(set(labels))]
+        handles = list(np.array(handles)[unique_indices])
+        ax.legend(handles=handles)