diff --git a/stixcore/io/product_processors/fits/processors.py b/stixcore/io/product_processors/fits/processors.py index e8a1351a..90630eea 100644 --- a/stixcore/io/product_processors/fits/processors.py +++ b/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"), ) diff --git a/stixcore/io/product_processors/tests/test_processors.py b/stixcore/io/product_processors/tests/test_processors.py index bb0fe4ec..32c66cd6 100644 --- a/stixcore/io/product_processors/tests/test_processors.py +++ b/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 diff --git a/stixcore/processing/tests/test_end2end.py b/stixcore/processing/tests/test_end2end.py index 52a2cfeb..e3091873 100644 --- a/stixcore/processing/tests/test_end2end.py +++ b/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.end2end @pytest.mark.remote_data -@pytest.mark.end2end def test_identical(orig_fits, current_fits): error_c = 0 error_files = list() diff --git a/stixcore/processing/tests/test_publish.py b/stixcore/processing/tests/test_publish.py index fee75ec3..a7affcc5 100644 --- a/stixcore/processing/tests/test_publish.py +++ b/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 diff --git a/stixcore/products/level3/flarelist.py b/stixcore/products/level3/flarelist.py index aef66704..188f6193 100644 --- a/stixcore/products/level3/flarelist.py +++ b/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 diff --git a/stixcore/products/level3/flarelistproduct.py b/stixcore/products/level3/flarelistproduct.py index 4a3fbdaa..ae7cb0d4 100644 --- a/stixcore/products/level3/flarelistproduct.py +++ b/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) diff --git a/stixcore/products/product.py b/stixcore/products/product.py index 15c9d246..af08ff34 100644 --- a/stixcore/products/product.py +++ b/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): + 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): diff --git a/stixcore/products/tests/test_factory.py b/stixcore/products/tests/test_factory.py index 68a140ed..8e0d5074 100644 --- a/stixcore/products/tests/test_factory.py +++ b/stixcore/products/tests/test_factory.py @@ -2,12 +2,16 @@ import pytest +from astropy.time import Time +from astropy.units import Quantity + from stixcore.data.test import test_data from stixcore.products.level0.quicklookL0 import LightCurve as LCL0 from stixcore.products.level1.quicklookL1 import LightCurve as LCL1 from stixcore.products.levelb.binary import LevelB from stixcore.products.product import Product from stixcore.time import SCETime +from stixcore.time.datetime import SCETimeDelta def test_ql_lb(): @@ -22,6 +26,22 @@ def test_ql_lb(): assert lb_prod.obt_beg == SCETime(coarse=664148503, fine=10710) +def test_read_timeformat(): + lq_scet = Product(test_data.products.L1_LightCurve_fits[0]) + + assert type(lq_scet.data["time"][0]) is SCETime + assert type(lq_scet.data["timedel"][0]) is SCETimeDelta + + lq_utc = Product(test_data.products.L1_LightCurve_fits[0], get_timeformat_from_TIMESYS=True) + assert type(lq_utc.data["time"][0]) is Time + assert type(lq_utc.data["timedel"][0]) is Quantity + + assert abs((lq_scet.scet_timerange.start - lq_utc.scet_timerange.start).coarse) < 1 + assert abs((lq_scet.scet_timerange.end - lq_utc.scet_timerange.end).coarse) < 1 + assert abs((lq_scet.utc_timerange.start - lq_utc.utc_timerange.start).to("s").value) < 0.2 + assert abs((lq_scet.utc_timerange.end - lq_utc.utc_timerange.end).to("s").value) < 0.2 + + # The fits file times maybe off by onescet time bin need to regenerate and test @pytest.mark.xfail def test_ql_l0(): diff --git a/stixcore/util/scripts/end2end_testing.py b/stixcore/util/scripts/end2end_testing.py index cf8cfa49..eda001c7 100644 --- a/stixcore/util/scripts/end2end_testing.py +++ b/stixcore/util/scripts/end2end_testing.py @@ -169,6 +169,19 @@ def end2end_pipeline(indir, fitsdir): if __name__ == "__main__": + p_scet = Product(Path("/data/stix/out/test/e2e/orig/solo_L1_stix-ql-variance_20210626_V02U.fits")) + p_utc = Product( + Path("/data/stix/out/test/e2e/orig/solo_L1_stix-ql-variance_20210626_V02U.fits"), + get_timeformat_from_TIMESYS=True, + ) + + end2end_pipeline( + indir=Path("/data/stix/out/test/e2e/orig"), + fitsdir=Path("/data/stix/out/test/e2e/current"), + ) + + quit() + if len(sys.argv) > 2: zippath = Path(sys.argv[1]) datapath = Path(sys.argv[2])