global ephemeris data cache

.gitignore

@@ -110,6 +110,7 @@ celerybeat.pid
 .env
 .venv
 env/
+venv/
 venv*/
 ENV/
 env.bak/

changelog/165.feature.rst

@@ -0,0 +1 @@
+add a internal ephemeris data cache. This allows to bypass hidden internal FIDO searches and downloads of the same data multiple times, improving performance when repeatedly querying for the same time ranges but also give the user a option to preload the cache with data.

stixpy/coordinates/tests/test_transforms.py

@@ -84,19 +84,32 @@ def test_stx_to_hpc_obstime_end():
     assert np.all(stix_coord.obstime_end.isclose(stix_coord_rt.obstime_end))
 
 
+def test_stx_to_hpc_obstime_end_x2():
+    # strange error that the call the second times crashes
+    test_stx_to_hpc_obstime_end()
+    test_stx_to_hpc_obstime_end()
+
+
 @pytest.mark.remote_data
 def test_get_aux_data():
     with pytest.raises(ValueError, match="No STIX pointing data found for time range"):
         _get_ephemeris_data(Time("2015-06-06"))  # Before the mission started
 
-    aux_data = _get_ephemeris_data(Time("2022-08-28T16:02:00"))
-    assert len(aux_data) == 1341
+    t1 = Time("2022-08-28T16:02:00")
+    aux_data = _get_ephemeris_data(t1)
+    assert len(aux_data) > 0
+    assert aux_data["time"].min() <= t1 <= aux_data["time_end"].max()
 
-    aux_data = _get_ephemeris_data(Time("2022-08-28T16:02:00"), end_time=Time("2022-08-28T16:04:00"))
-    assert len(aux_data) == 1341
+    t2 = Time("2022-08-28T16:04:00")
+    aux_data = _get_ephemeris_data(t1, end_time=t2)
+    assert len(aux_data) > 0
+    assert t1 <= aux_data["time"].min() < t2 <= aux_data["time_end"].max()
 
-    aux_data = _get_ephemeris_data(Time("2022-08-28T23:58:00"), end_time=Time("2022-08-29T00:02:00"))
-    assert len(aux_data) == 2691
+    t1 = Time("2022-08-28T23:58:00")
+    t2 = Time("2022-08-29T00:02:00")
+    aux_data = _get_ephemeris_data(t1, end_time=t2)
+    assert len(aux_data) > 0
+    assert aux_data["time"].min() <= t1 < t2 <= aux_data["time_end"].max()
 
 
 @pytest.mark.remote_data
@@ -135,9 +148,9 @@ def test_get_hpc_info_shapes():
     roll3, solo_heeq3, stix_pointing3 = get_hpc_info(t[5])
 
     assert_quantity_allclose(roll1[5], roll2)
-    assert_quantity_allclose(solo_heeq1[5, :], solo_heeq2[0, :])
-    assert_quantity_allclose(stix_pointing1[5, :], stix_pointing2[0, :])
+    assert_quantity_allclose(solo_heeq1[5, :], solo_heeq2)
+    assert_quantity_allclose(stix_pointing1[5, :], stix_pointing2)
 
     assert_quantity_allclose(roll3, roll2[0])
-    assert_quantity_allclose(solo_heeq3, solo_heeq2[0, :])
-    assert_quantity_allclose(stix_pointing3, stix_pointing2[0, :])
+    assert_quantity_allclose(solo_heeq3, solo_heeq2)
+    assert_quantity_allclose(stix_pointing3, stix_pointing2)

stixpy/coordinates/transforms.py

@@ -1,20 +1,21 @@
 import warnings
-from functools import lru_cache
 
 import astropy.coordinates as coord
 import astropy.units as u
 import numpy as np
 from astropy.coordinates import frame_transform_graph
 from astropy.coordinates.matrix_utilities import matrix_transpose, rotation_matrix
-from astropy.io import fits
-from astropy.table import QTable, vstack
+from astropy.table import vstack
 from astropy.time import Time
 from sunpy.coordinates import HeliographicStonyhurst, Helioprojective
 from sunpy.net import Fido
 from sunpy.net import attrs as a
 
 from stixpy.coordinates.frames import STIXImaging
+from stixpy.product.product_factory import Product
+from stixpy.product.sources.anc import Ephemeris
 from stixpy.utils.logging import get_logger
+from stixpy.utils.table_lru import TableLRUCache
 
 STIX_X_SHIFT = 26.1 * u.arcsec  # fall back to this when non sas solution available
 STIX_Y_SHIFT = 58.2 * u.arcsec  # fall back to this when non sas solution available
@@ -23,7 +24,23 @@
 
 logger = get_logger(__name__)
 
-__all__ = ["get_hpc_info", "stixim_to_hpc", "hpc_to_stixim"]
+__all__ = ["get_hpc_info", "stixim_to_hpc", "hpc_to_stixim", "STIX_EPHEMERIS_CACHE", "load_ephemeris_fits_to_cache"]
+
+# Create a global cache for STIX ephemeris data
+STIX_EPHEMERIS_CACHE = TableLRUCache("STIX_EPHEMERIS_CACHE", maxsize=300000, default_bin_duration=64 * u.s)
+
+
+def load_ephemeris_fits_to_cache(anc_file):
+    """
+    Load ephemeris data from a fits file into the ephemeris cache.
+    """
+    logger.info(f"Loading STIX ephemeris data from {anc_file}. into cache")
+    try:
+        anc = Product(anc_file, data_only=True)
+        if isinstance(anc, Ephemeris):
+            STIX_EPHEMERIS_CACHE.put(anc.data)
+    except (OSError, ValueError) as e:
+        logger.error(f"Error loading STIX ephemeris data from {anc_file}: {e}")
 
 
 def _get_rotation_matrix_and_position(obstime, obstime_end=None):
@@ -57,7 +74,7 @@ def _get_rotation_matrix_and_position(obstime, obstime_end=None):
 
 def get_hpc_info(times, end_time=None):
     r"""
-    Get STIX pointing and SO location from L2 aspect files.
+    Get STIX pointing and SO location from ANC aspect files.
 
     Parameters
     ----------
@@ -68,11 +85,20 @@ def get_hpc_info(times, end_time=None):
     -------
 
     """
-    aux = _get_ephemeris_data(times.min(), end_time or times.max())
+    start_time = times.min()
+    end_time = end_time or times.max()
+
+    aux = _get_ephemeris_data(start_time, end_time)
+
+    # indices = np.argwhere(
+    #     ((aux["time"] >= start_time) | (aux["time_end"] >= start_time))
+    #     & ((aux["time_end"] <= end_time) | (aux["time"] <= end_time))
+    # )
 
     indices = np.argwhere((aux["time"] >= times.min()) & (aux["time"] <= times.max()))
     if end_time is not None:
         indices = np.argwhere((aux["time"] >= times.min()) & (aux["time"] <= end_time))
+
     indices = indices.flatten()
 
     if end_time is not None and times.size == 1 and indices.size >= 2:
@@ -120,8 +146,8 @@ def get_hpc_info(times, end_time=None):
             << aux["solo_loc_heeq_zxy"].unit
         )
 
-        sas_x = np.interp(x, xp, aux["y_srf"])
-        sas_y = np.interp(x, xp, aux["z_srf"])
+        sas_x = np.interp(x, xp, aux["y_srf"].value) << aux["y_srf"].unit
+        sas_y = np.interp(x, xp, aux["z_srf"].value) << aux["z_srf"].unit
         if x.size == 1:
             good_sas = [True] if np.interp(x, xp, aux["sas_ok"]).astype(bool) else []
         else:
@@ -154,10 +180,9 @@ def get_hpc_info(times, end_time=None):
     return roll, solo_heeq, stix_pointing
 
 
-@lru_cache
 def _get_ephemeris_data(start_time, end_time=None):
     r"""
-    Search, download and read L2 pointing data.
+    Search, in cache or download and read ANC pointing data.
 
     Parameters
     ----------
@@ -170,39 +195,44 @@ def _get_ephemeris_data(start_time, end_time=None):
     """
     if end_time is None:
         end_time = start_time
-    # Find, download, read aux file with pointing, sas and position information
-    logger.debug(f"Searching for AUX data: {start_time} - {end_time}")
-    query = Fido.search(
-        a.Time(start_time, end_time),
-        a.Instrument.stix,
-        a.Level.anc,
-        a.stix.DataType.asp,
-        a.stix.DataProduct.asp_ephemeris,
-    )
-    if len(query["stix"]) == 0:
-        raise ValueError(f"No STIX pointing data found for time range {start_time} to {end_time}.")
-
-    logger.debug(f"Downloading {len(query['stix'])} AUX files")
-    aux_files = Fido.fetch(query["stix"])
-    if len(aux_files.errors) > 0:
-        raise ValueError("There were errors downloading the data.")
-    # Read and extract data
-    logger.debug("Loading and extracting AUX data")
-
-    aux_data = []
-    for aux_file in aux_files:
-        hdu = fits.getheader(aux_file, ext=0)
-        aux = QTable.read(aux_file, hdu=2)
-        date_beg = Time(hdu.get("DATE-BEG"))
-        aux["time"] = (
-            date_beg + aux["time"] - 32 * u.s
-        )  # Shift AUX data by half a time bin (starting time vs. bin centre)
-        aux_data.append(aux)
-
-    aux = vstack(aux_data)
-    aux.sort(keys=["time"])
-
-    return aux
+
+    logger.warning(f"Getting STIX ephemeris data for {start_time} to {end_time}")
+    from_cache = STIX_EPHEMERIS_CACHE.get(start_time, end_time)
+    if from_cache is None:
+        # Find, download, read aux file with pointing, sas and position information
+        logger.warning(f"FIDO searching for AUX data: {start_time} - {end_time}")
+        query = Fido.search(
+            a.Time(start_time, end_time),
+            a.Instrument.stix,
+            a.Level.anc,
+            a.stix.DataType.asp,
+            a.stix.DataProduct.asp_ephemeris,
+        )
+        if len(query["stix"]) == 0:
+            raise ValueError(f"No STIX pointing data found for time range {start_time} to {end_time}.")
+        else:
+            query["stix"].filter_for_latest_version()
+        logger.debug(f"Downloading {len(query['stix'])} AUX files")
+        aux_files = Fido.fetch(query["stix"])
+        if len(aux_files.errors) > 0:
+            raise ValueError("There were errors downloading the data.")
+        # Read and extract data
+        logger.debug("Loading and extracting ANC data")
+
+        anc_data = []
+        for aux_file in aux_files:
+            anc = Product(aux_file, data_only=True)
+            anc_data.append(anc.data)
+
+        ephemeris = vstack(anc_data)
+        ephemeris.sort(keys=["time"])
+
+        STIX_EPHEMERIS_CACHE.put(ephemeris)
+
+        return ephemeris
+    else:
+        logger.debug(f"Using cached ANC data: {start_time} - {end_time}")
+        return from_cache
 
 
 @frame_transform_graph.transform(coord.FunctionTransform, STIXImaging, Helioprojective)

stixpy/product/product.py

@@ -23,7 +23,7 @@ def __init_subclass__(cls, **kwargs):
 
 
 class GenericProduct(BaseProduct):
-    def __init__(self, *, meta, control, data, idb_versions=None, energies=None):
+    def __init__(self, *, meta, control, data, idb_versions=None, energies=None, **kwargs):
         """
         Generic product composed of meta, control, data and optionally idb, and energy information
 

stixpy/product/product_factory.py

@@ -116,7 +116,16 @@ def _read_file(self, fname, **kwargs):
             raise FileError(f"File '{fname}' is not a STIX fits file.")
 
         data = {"meta": hdul[0].header}
-        for name in ["CONTROL", "DATA", "IDB_VERSIONS", "ENERGIES"]:
+
+        # determine extensions to load
+        # allow for just data extension if requested
+        if kwargs.get("data_only", False):
+            extensions = ["CONTROL", "DATA"]
+        # normally read all extensions
+        else:
+            extensions = ["CONTROL", "DATA", "IDB_VERSIONS", "ENERGIES"]
+
+        for name in extensions:
             try:
                 data[name.lower()] = read_qtable(fname, hdu=name)
             except KeyError as e:

stixpy/product/sources/__init__.py

@@ -1,3 +1,4 @@
+from stixpy.product.sources.anc import *
 from stixpy.product.sources.housekeeping import *
 from stixpy.product.sources.quicklook import *
 from stixpy.product.sources.science import *

stixpy/product/sources/anc.py

@@ -0,0 +1,55 @@
+import astropy.units as u
+from astropy.time import Time
+from astropy.units import Quantity
+from sunpy.time import TimeRange
+
+from stixpy.product.product import L1Product
+
+__all__ = ["ANCProduct", "Ephemeris"]
+
+
+class ANCProduct(L1Product):
+    """
+    Basic ANC
+    """
+
+    @property
+    def time(self) -> Time:
+        return self.data["time"]
+
+    @property
+    def exposure_time(self) -> Quantity[u.s]:
+        return self.data["timedel"].to(u.s)
+
+    @property
+    def time_range(self) -> TimeRange:
+        """
+        A `sunpy.time.TimeRange` for the data.
+        """
+        return TimeRange(self.time[0] - self.exposure_time[0] / 2, self.time[-1] + self.exposure_time[-1] / 2)
+
+
+class Ephemeris(ANCProduct):
+    """
+    Ephemeris data in daily files normal with 64s time resolution
+    """
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        # TODO remove this when we have a better solution
+        # Shift ANC data by half a time bin (starting time vs. bin centre)
+        self.data["time_end"] = self.data["time"] + 32 * u.s
+        self.data["time"] = self.data["time"] - 32 * u.s
+        self.data["timedel"] = 64 * u.s
+        self.data["time_utc"] = [Time(t, format="isot", scale="utc") for t in self.data["time_utc"]]
+
+    @classmethod
+    def is_datasource_for(cls, *, meta, **kwargs):
+        """Determines if meta data meach Raw Pixel Data"""
+        service_subservice_ssid = tuple(meta[name] for name in ["STYPE", "SSTYPE", "SSID"])
+        level = meta["level"]
+        if service_subservice_ssid == (0, 0, 1) and level == "ANC":
+            return True
+
+    def __repr__(self):
+        return f"{self.__class__.__name__}\n    {self.time_range}"