Merge develop

.github/workflows/python-package.yml

@@ -16,12 +16,12 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        python-version: ["3.8", "3.11"]
+        python-version: ["3.9", "3.12"]
 
     steps:
-    - uses: actions/checkout@v3
+    - uses: actions/checkout@v4
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v3
+      uses: actions/setup-python@v5
       with:
         python-version: ${{ matrix.python-version }}
     - name: Install dependencies
@@ -37,7 +37,10 @@ jobs:
         flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics
     - name: Test with pytest
       run: |
-        python -m setuptools_scm
+        python -c "import setuptools_scm; version = setuptools_scm.get_version(); setuptools_scm.dump_version('sospice', version, '_version.py')"
+        ls -l _version.py || echo "_version.py not found"
+        ls -l sospice/_version.py ||  echo "sospice/_version.py not found"
         pytest --cov --cov-report=xml
+# Codecov should be updated to @v4 but then there is a key not certified / token not found error
     - name: Upload coverage reports to Codecov with GitHub Action
       uses: codecov/codecov-action@v3

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
 ------------------
 

changelog/35.feature.rst

@@ -0,0 +1 @@
+Add local sigma clipping.

changelog/43.feature

@@ -0,0 +1 @@
+Add field of view plotting on a background (including a SDO/HMI synoptic map or a Solar Orbiter/EUI/FSI image).

changelog/README.rst

@@ -12,7 +12,7 @@ The ``CHANGELOG`` will be read by users, so this description should be aimed at
 
 Make sure to use full sentences with correct case and punctuation, for example::
 
-    Add support for PSF correction`.
+    Add support for PSF correction.
 
 Please try to use Sphinx intersphinx using backticks.
 

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/index.rst

@@ -31,6 +31,7 @@ Documentation for this package is available on `Read the Docs <https://sospice.r
    getting_started
    calibration
    catalog
+   utilities
    citation
    api
    develop

docs/source/utilities.rst

@@ -0,0 +1,41 @@
+Utilities
+=========
+
+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"

requirements.txt

@@ -1,6 +1,8 @@
 astropy>=4.0
+beautifulsoup4>=4.0
 black>=23.0
 build>=0.9.0
+drms>=0.7.0
 flake8>=6.0
 matplotlib>=3.7.0
 numpy>=1.24.0
@@ -13,10 +15,12 @@ pytest>=7.0
 pytest-cov>=4.0
 requests>=2.30
 scipy>=1.10
-setuptools>=45
-setuptools_scm[toml]>=7.0
+setuptools>=61
+setuptools_scm[toml]>=8.0
 Sphinx>=6.2,<7.0
 sphinx-rtd-theme>=1.2.0
 sunpy>=4.0
+sunpy-soar>=1.10
 towncrier>=22.0.0
 twine>=4.0
+zeep>=4.0

sospice/__init__.py

@@ -3,3 +3,5 @@
 from .catalog.release import Release
 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
@@ -18,7 +22,7 @@
     filename: str
         A file name (or URL) for the catalog
     release_tag: str
-        A release tag. The catalog is fetched online and dowloaded to the astropy cache.
+        A release tag. The catalog is fetched online and downloaded to the astropy cache.
     data_frame: pandas.DataFrame
         A pandas DataFrame to be used as SPICE catalog. Some basic checks are made to ensure
         that is can be used as a SPICE catalog.
@@ -160,7 +164,7 @@
 
     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 @@
                 .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)