add option to fits reader to open L1 products also in UTC >L1 UTC per default (via keyword in header) (#447)

* change internal time system of L1 products to be SCETime or UTC

* fixes after review

* add warning when time out converting with spice

* rework to not auto convert times in dta table at processing time only on fits open

* cleanup after review

* pytz.UTC

Author: nicHoch

stixcore/io/product_processors/fits/processors.py

@@ -723,22 +723,23 @@ def generate_primary_header(cls, filename, product, *, version=0):
         # if not isinstance(product.obt_beg, SCETime):
         #     raise ValueError("Expected SCETime as time format")
 
+        scet_timerange = product.scet_timerange
         headers = FitsProcessor.generate_common_header(filename, product, version=version) + (
             # Name, Value, Comment
             # ('MJDREF', product.obs_beg.mjd),
             # ('DATEREF', product.obs_beg.fits),
-            ("OBT_BEG", product.scet_timerange.start.as_float().value, "Start acquisition time in OBT"),
-            ("OBT_END", product.scet_timerange.end.as_float().value, "End acquisition time in OBT"),
+            ("OBT_BEG", scet_timerange.start.as_float().value, "Start acquisition time in OBT"),
+            ("OBT_END", scet_timerange.end.as_float().value, "End acquisition time in OBT"),
             ("TIMESYS", "OBT", "System used for time keywords"),
             ("LEVEL", "L0", "Processing level of the data"),
-            ("DATE-OBS", product.scet_timerange.start.to_string(), "Depreciated, same as DATE-BEG"),
-            ("DATE-BEG", product.scet_timerange.start.to_string(), "Start time of observation"),
-            ("DATE-AVG", product.scet_timerange.avg.to_string(), "Average time of observation"),
-            ("DATE-END", product.scet_timerange.end.to_string(), "End time of observation"),
+            ("DATE-OBS", scet_timerange.start.to_string(), "Depreciated, same as DATE-BEG"),
+            ("DATE-BEG", scet_timerange.start.to_string(), "Start time of observation"),
+            ("DATE-AVG", scet_timerange.avg.to_string(), "Average time of observation"),
+            ("DATE-END", scet_timerange.end.to_string(), "End time of observation"),
             ("DATAMIN", product.dmin, "Minimum valid physical value"),
             ("DATAMAX", product.dmax, "Maximum valid physical value"),
             ("BUNIT", product.bunit, "Units of physical value, after application of BSCALE, BZERO"),
-            ("XPOSURE", product.exposure, "[s] shortest exposure time"),
+            ("XPOSURE", product.min_exposure, "[s] shortest exposure time"),
             ("XPOMAX", product.max_exposure, "[s] maximum exposure time"),
         )
 
@@ -782,7 +783,7 @@ def generate_primary_header(self, filename, product, *, version=0):
             ("DATAMIN", empty_if_nan(product.dmin), "Minimum valid physical value"),
             ("DATAMAX", empty_if_nan(product.dmax), "Maximum valid physical value"),
             ("BUNIT", product.bunit, "Units of physical value, after application of BSCALE, BZERO"),
-            ("XPOSURE", empty_if_nan(product.exposure), "[s] shortest exposure time"),
+            ("XPOSURE", empty_if_nan(product.min_exposure), "[s] shortest exposure time"),
             ("XPOMAX", empty_if_nan(product.max_exposure), "[s] maximum exposure time"),
         )
 
@@ -1000,7 +1001,7 @@ def generate_primary_header(self, filename, product, *, version=0):
             ("DATAMIN", empty_if_nan(product.dmin), "Minimum valid physical value"),
             ("DATAMAX", empty_if_nan(product.dmax), "Maximum valid physical value"),
             ("BUNIT", product.bunit, "Units of physical value, after application of BSCALE, BZERO"),
-            ("XPOSURE", empty_if_nan(product.exposure), "[s] shortest exposure time"),
+            ("XPOSURE", empty_if_nan(product.min_exposure), "[s] shortest exposure time"),
             ("XPOMAX", empty_if_nan(product.max_exposure), "[s] maximum exposure time"),
         )
 

stixcore/io/product_processors/tests/test_processors.py

@@ -191,13 +191,13 @@ def test_count_data_mixin(p_file):
     p = Product(p_file)
     assert p.dmin == p.data["counts"].min().value
     assert p.dmax == p.data["counts"].max().value
-    assert p.exposure == p.data["timedel"].min().as_float().to_value()
+    assert p.min_exposure == p.data["timedel"].min().as_float().to_value()
     assert p.max_exposure == p.data["timedel"].max().as_float().to_value()
 
     test_data = {
         "DATAMAX": p.dmax,
         "DATAMIN": p.dmin,
-        "XPOSURE": p.exposure,
+        "XPOSURE": p.min_exposure,
         "XPOMAX": p.max_exposure,
         "BUNIT": "counts",
     }
@@ -257,7 +257,7 @@ def test_level1_processor_generate_primary_header(product, soop_manager):
     product.dmax = 1
     product.dunit = ""
     product.max_exposure = 1
-    product.exposure = 1
+    product.min_exposure = 1
     product.service_type = 1
     product.service_subtype = 2
     product.ssid = 3

stixcore/processing/tests/test_end2end.py

@@ -73,8 +73,8 @@ def test_complete(orig_fits, current_fits):
         raise ValueError(f"{error_c} errors out of {len(orig_fits)}\nnumber of fits files differ")
 
 
-@pytest.mark.remote_data
 @pytest.mark.end2end
+@pytest.mark.remote_data
 def test_identical(orig_fits, current_fits):
     error_c = 0
     error_files = list()

stixcore/processing/tests/test_publish.py

@@ -138,7 +138,7 @@ def test_publish_fits_to_esa_incomplete(product, out_dir):
         product.date_obs = beg
         product.date_beg = beg
         product.date_end = end
-        product.exposure = 2
+        product.min_exposure = 2
         product.max_exposure = 3
         product.dmin = 2
         product.dmax = 3
@@ -247,7 +247,7 @@ def test_fits_incomplete_switch_over(out_dir):
             product.date_obs = beg
             product.date_beg = beg
             product.date_end = end
-            product.exposure = 2
+            product.min_exposure = 2
             product.max_exposure = 3
             product.dmin = 2
             product.dmax = 3
@@ -385,7 +385,7 @@ def test_publish_fits_to_esa(product, out_dir):
     product.date_obs = beg
     product.date_beg = beg
     product.date_end = end
-    product.exposure = 2
+    product.min_exposure = 2
     product.max_exposure = 3
     product.dmin = 2
     product.dmax = 3

stixcore/products/CAL/energy.py

@@ -81,7 +81,7 @@ def bunit(self):
         return "keV"
 
     @property
-    def exposure(self):
+    def min_exposure(self):
         # default for FITS HEADER
         return self.control["integration_time"].min().to_value(u.s)
 

stixcore/products/level3/flarelist.py

@@ -499,6 +499,9 @@ def utc_timerange(self):
     @property
     def scet_timerange(self):
         tr = self.utc_timerange
+        logger.warning(
+            "scet_timerange will be approximated using Spice. Better to work with utc_timerange property to avoid automatic time conversion"
+        )
         start = SCETime.from_string(Spice.instance.datetime_to_scet(tr.start)[2:])
         end = SCETime.from_string(Spice.instance.datetime_to_scet(tr.end)[2:])
         return SCETimeRange(start=start, end=end)
@@ -515,7 +518,7 @@ def dmax(self):
         return (self.data["lc_peak"].sum(axis=1)).max().value if len(self.data) > 0 else np.nan
 
     @property
-    def exposure(self):
+    def min_exposure(self):
         return self.data["duration"].min().to_value("s") if len(self.data) > 0 else np.nan
 
     @property
@@ -633,6 +636,9 @@ def utc_timerange(self):
     @property
     def scet_timerange(self):
         tr = self.utc_timerange
+        logger.warning(
+            "scet_timerange will be approximated using Spice. Better to work with utc_timerange property to avoid automatic time conversion"
+        )
         start = SCETime.from_string(Spice.instance.datetime_to_scet(tr.start)[2:])
         end = SCETime.from_string(Spice.instance.datetime_to_scet(tr.end)[2:])
         return SCETimeRange(start=start, end=end)
@@ -649,7 +655,7 @@ def dmax(self):
         return (self.data["lc_peak"].sum(axis=1)).max().value if len(self.data) > 0 else np.nan
 
     @property
-    def exposure(self):
+    def min_exposure(self):
         return self.data["duration"].min().to_value("s") if len(self.data) > 0 else np.nan
 
     @property

stixcore/products/level3/flarelistproduct.py

@@ -4,9 +4,12 @@
 from stixcore.ephemeris.manager import Spice
 from stixcore.products.product import GenericProduct, L3Mixin
 from stixcore.time.datetime import SCETime, SCETimeRange
+from stixcore.util.logging import get_logger
 
 __all__ = ["FlareListProduct", "PeekPreviewImage"]
 
+logger = get_logger(__name__)
+
 
 class FlareListProduct(GenericProduct, L3Mixin):
     """Product not based on direct TM data but on time ranges defined in flare lists.
@@ -47,6 +50,9 @@ def utc_timerange(self):
     @property
     def scet_timerange(self):
         tr = self.utc_timerange
+        logger.warning(
+            "scet_timerange will be approximated using Spice. Better to work with utc_timerange property to avoid automatic time conversion"
+        )
         start = SCETime.from_string(Spice.instance.datetime_to_scet(tr.start)[2:])
         end = SCETime.from_string(Spice.instance.datetime_to_scet(tr.end)[2:])
         return SCETimeRange(start=start, end=end)

stixcore/products/product.py

@@ -3,6 +3,8 @@
 from itertools import chain
 
 import numpy as np
+import pytz
+from sunpy.time.timerange import TimeRange
 from sunpy.util.datatype_factory_base import (
     BasicRegistrationFactory,
     MultipleMatchError,
@@ -17,6 +19,7 @@
 
 import stixcore.processing.decompression as decompression
 import stixcore.processing.engineering as engineering
+from stixcore.ephemeris.manager import Spice
 from stixcore.idb.manager import IDBManager
 from stixcore.time import SCETime, SCETimeDelta, SCETimeRange
 from stixcore.tmtc.packet_factory import Packet
@@ -47,7 +50,7 @@
 
 # date when the min integration time was changed from 1.0s to 0.5s needed to fix count and time
 # offset issue
-MIN_INT_TIME_CHANGE = datetime(2021, 9, 6, 13)
+MIN_INT_TIME_CHANGE = datetime(2021, 9, 6, 13, tzinfo=pytz.UTC)
 
 
 def read_qtable(file, hdu, hdul=None):
@@ -211,7 +214,9 @@ def get_cls_processing_version(cls):
 
 class ProductFactory(BasicRegistrationFactory):
     def __call__(self, *args, **kwargs):
-        if len(args) == 1 and len(kwargs) == 0:
+        get_timeformat_from_TIMESYS = kwargs.get("get_timeformat_from_TIMESYS", False)
+
+        if len(args) == 1:
             if isinstance(args[0], (str, Path)):
                 file_path = Path(args[0])
                 pri_header = fits.getheader(file_path)
@@ -258,8 +263,24 @@ def __call__(self, *args, **kwargs):
                         ssid = 34
 
                 if level not in ["LB", "LL01"] and "timedel" in data.colnames and "time" in data.colnames:
-                    data["timedel"] = SCETimeDelta(data["timedel"])
-                    offset = SCETime.from_float(pri_header["OBT_BEG"] * u.s)
+                    if level in ["L0", "L1"] and not get_timeformat_from_TIMESYS:
+                        # L0 and L1 date are open by default in SCETime format so we can directly apply the timedelta
+                        data["timedel"] = SCETimeDelta(data["timedel"])
+                        offset = SCETime.from_float(pri_header["OBT_BEG"] * u.s)
+                    else:
+                        # in L2 and higher the time format should not be in SCETime format
+                        # select the time format based on available header keywords
+                        offset = None
+                        if pri_header.get("TIMESYS", "") == "UTC":
+                            try:
+                                offset = Time(pri_header["DATE-OBS"])
+                            except ValueError:
+                                offset = None
+
+                        # fallback to OBT_BEG if no TIMESYS=UTC or DATE-OBS is present or can not be parsed
+                        if offset is None:
+                            offset = SCETime.from_float(pri_header["OBT_BEG"] * u.s)
+                            data["timedel"] = SCETimeDelta(data["timedel"])
 
                     try:
                         control["time_stamp"] = SCETime.from_float(control["time_stamp"])
@@ -535,10 +556,25 @@ def __init__(
 
     @property
     def scet_timerange(self):
-        return SCETimeRange(
-            start=self.data["time"][0] - self.data["timedel"][0] / 2,
-            end=self.data["time"][-1] + self.data["timedel"][-1] / 2,
-        )
+        if isinstance(self.data["time"], SCETime):
+            return SCETimeRange(
+                start=self.data["time"][0] - self.data["timedel"][0] / 2,
+                end=self.data["time"][-1] + self.data["timedel"][-1] / 2,
+            )
+        else:
+            logger.warning(
+                "internal time format is not in SCETime format, scet_timerange will be approximated using Spice. Better to work with utc_timerange property to avoid automatic time conversion"
+            )
+            start_str = Spice.instance.datetime_to_scet((self.data["time"][0] - self.data["timedel"][0] / 2).datetime)
+            end_str = Spice.instance.datetime_to_scet((self.data["time"][-1] + self.data["timedel"][-1] / 2).datetime)
+            if "/" in start_str:
+                start_str = start_str.split("/")[-1]
+            if "/" in end_str:
+                end_str = end_str.split("/")[-1]
+            return SCETimeRange(
+                start=SCETime.from_string(start_str),
+                end=SCETime.from_string(end_str),
+            )
 
     @property
     def raw(self):
@@ -564,7 +600,7 @@ def bunit(self):
         return " "
 
     @property
-    def exposure(self):
+    def min_exposure(self):
         # default for FITS HEADER
         return 0.0
 
@@ -644,6 +680,12 @@ def __add__(self, other):
         if not isinstance(other, type(self)):
             raise TypeError(f"Products must of same type not {type(self)} and {type(other)}")
 
+        if "time" in self.data.colnames and "time" in other.data.colnames:
+            if type(self.data["time"]) is not type(other.data["time"]):
+                raise TypeError(
+                    f"Products must have the same time format not {type(self.data['time'])} and {type(other.data['time'])}"
+                )
+
         # make a deep copy of the data and control
         other_control = other.control[:]
         other_data = other.data[:]
@@ -667,8 +709,10 @@ def __add__(self, other):
 
         # Fits write we do np.around(time - start_time).as_float().to(u.cs)).astype("uint32"))
         # So need to do something similar here to avoid comparing un-rounded value to rounded values
-        data["time_float"] = np.around((data["time"] - data["time"].min()).as_float().to("cs"))
-
+        if isinstance(data["time"], SCETime):
+            data["time_float"] = np.around((data["time"] - data["time"].min()).as_float().to("cs"))
+        else:  # datetime or Time
+            data["time_float"] = np.around((data["time"] - data["time"].min()).to("cs"))
         # remove duplicate data based on time bin and sort the data
         data = unique(data, keys=["time_float"])
         # data.sort(["time_float"])
@@ -863,12 +907,18 @@ def bunit(self):
         return "counts"
 
     @property
-    def exposure(self):
-        return self.data["timedel"].as_float().min().to_value("s")
+    def min_exposure(self):
+        if isinstance(self.data["timedel"], SCETimeDelta):
+            return self.data["timedel"].as_float().min().to_value("s")
+        else:
+            return self.data["timedel"].min().to_value("s")
 
     @property
     def max_exposure(self):
-        return self.data["timedel"].as_float().max().to_value("s")
+        if isinstance(self.data["timedel"], SCETimeDelta):
+            return self.data["timedel"].as_float().max().to_value("s")
+        else:
+            return self.data["timedel"].max().to_value("s")
 
 
 class EnergyChannelsMixin:
@@ -925,7 +975,13 @@ class L1Mixin(FitsHeaderMixin):
 
     @property
     def utc_timerange(self):
-        return self.scet_timerange.to_timerange()
+        if isinstance(self.data["time"], SCETime):
+            return self.scet_timerange.to_timerange()
+        else:
+            return TimeRange(
+                (self.data["time"][0] - self.data["timedel"][0] / 2),
+                (self.data["time"][-1] + self.data["timedel"][-1] / 2),
+            )
 
     @classmethod
     def from_level0(cls, l0product, parent=""):
@@ -951,10 +1007,10 @@ def from_level0(cls, l0product, parent=""):
             if idbs[0] < (2, 26, 36) and len(l1.data) > 1:
                 # Check if request was at min configured time resolution
                 if (
-                    l1.utc_timerange.start.datetime < MIN_INT_TIME_CHANGE
+                    l0product.scet_timerange.start.to_datetime() < MIN_INT_TIME_CHANGE
                     and l1.data["timedel"].as_float().min() == 1 * u.s
                 ) or (
-                    l1.utc_timerange.start.datetime >= MIN_INT_TIME_CHANGE
+                    l0product.scet_timerange.start.to_datetime() >= MIN_INT_TIME_CHANGE
                     and l1.data["timedel"].as_float().min() == 0.5 * u.s
                 ):
                     l1.data["timedel"][1:-1] = l1.data["timedel"][:-2]
@@ -972,7 +1028,13 @@ def from_level0(cls, l0product, parent=""):
 class L2Mixin(FitsHeaderMixin):
     @property
     def utc_timerange(self):
-        return self.scet_timerange.to_timerange()
+        if isinstance(self.data["time"], SCETime):
+            return self.scet_timerange.to_timerange()
+        else:
+            return TimeRange(
+                (self.data["time"][0] - self.data["timedel"][0] / 2).datetime,
+                (self.data["time"][-1] + self.data["timedel"][-1] / 2).datetime,
+            )
 
     @classmethod
     def get_additional_extensions(cls):