Flexible way to fill DATAMAX/MIN BUNIT EXPOSURE in FITS header keywor…

…ds (#393)

* bump all product versions to 2

* add properties for dedicated fits header keywords with class inhe. / override

* add test and minor fixes after review

* fix XPOSURE/XPOMAX header error

* fix tests for dmax...

* relax end2end fitsdiff

* less tolerance in end2end test

stixcore/io/fits/processors.py

@@ -695,13 +695,6 @@ def generate_primary_header(cls, filename, product, *, version=0):
         #     raise ValueError(f"Try to crate FITS file L0 for {product.level} data product")
         # if not isinstance(product.obt_beg, SCETime):
         #     raise ValueError("Expected SCETime as time format")
-        dmin = 0.0
-        dmax = 0.0
-        bunit = ' '
-        if 'counts' in product.data.colnames:
-            dmax = product.data['counts'].max().value
-            dmin = product.data['counts'].min().value
-            bunit = 'counts'
 
         headers = FitsProcessor.generate_common_header(filename, product, version=version) + (
             # Name, Value, Comment
@@ -716,9 +709,11 @@ def generate_primary_header(cls, filename, product, *, version=0):
             ('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'),
-            ('DATAMIN', dmin, 'Minimum valid physical value'),
-            ('DATAMAX', dmax, 'Maximum valid physical value'),
-            ('BUNIT', bunit, 'Units of physical value, after application of BSCALE, BZERO')
+            ('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'),
+            ('XPOMAX', product.max_exposure, '[s] maximum exposure time')
         )
 
         return headers
@@ -745,20 +740,12 @@ def generate_primary_header(self, filename, product, *, version=0):
 
         headers = FitsProcessor.generate_common_header(filename, product, version=version)
 
-        dmin = 0.0
-        dmax = 0.0
-        bunit = ' '
-        exposure = 0.0
-        if 'counts' in product.data.colnames:
-            dmax = product.data['counts'].max().value
-            dmin = product.data['counts'].min().value
-            bunit = 'counts'
-            exposure = product.data['timedel'].as_float().min().to_value('s')
         data_headers = (
-            ('DATAMIN', dmin, 'Minimum valid physical value'),
-            ('DATAMAX', dmax, 'Maximum valid physical value'),
-            ('BUNIT', bunit, 'Units of physical value, after application of BSCALE, BZERO'),
-            ('XPOSURE', exposure, '[s] shortest exposure time')
+            ('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'),
+            ('XPOMAX', product.max_exposure, '[s] maximum exposure time')
         )
 
         soop_keywords = SOOPManager.instance.get_keywords(start=product.utc_timerange.start,
@@ -972,6 +959,11 @@ def generate_primary_header(self, filename, product, *, version=0):
             ('VERS_CFG', str(stixcore.__version_conf__),
              'Version of the common instrument configuration package'),
             ('HISTORY', 'Processed by STIXCore L2'),
+            ('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'),
+            ('XPOMAX', product.max_exposure, '[s] maximum exposure time')
         )
 
         return L1headers, L2headers

stixcore/io/fits/tests/test_processors.py

@@ -7,6 +7,7 @@
 
 from stixcore.data.test import test_data
 from stixcore.io.fits.processors import FitsL0Processor, FitsL1Processor, FitsLBProcessor
+from stixcore.products.product import Product
 from stixcore.soop.manager import SOOPManager
 from stixcore.time import SCETime, SCETimeRange
 
@@ -145,6 +146,7 @@ def test_level0_processor_generate_primary_header(datetime, product):
     product.obs_beg = SCETime(coarse=0, fine=0)
     product.obs_avg = SCETime(coarse=0, fine=2 ** 15)
     product.obs_end = SCETime(coarse=1, fine=2 ** 15)
+
     product.scet_timerange = SCETimeRange(start=product.obs_beg, end=product.obs_end)
     product.raw = ['packet1.xml', 'packet2.xml']
     product.parent = ['lb1.fits', 'lb2.fts']
@@ -177,6 +179,31 @@ def test_level0_processor_generate_primary_header(datetime, product):
             assert value == test_data[name]
 
 
+@pytest.mark.parametrize('p_file', [test_data.products.L0_LightCurve_fits[0],
+                                    test_data.products.L1_LightCurve_fits[0]],
+                         ids=["ql_l0", "ql_l1"])
+def test_count_data_mixin(p_file):
+    processor = FitsL0Processor('some/path')
+    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.max_exposure == p.data["timedel"].max().as_float().to_value()
+
+    test_data = {
+        "DATAMAX": p.dmax,
+        "DATAMIN": p.dmin,
+        "XPOSURE": p.exposure,
+        "XPOMAX": p.max_exposure,
+        "BUNIT": "counts"
+    }
+
+    header = processor.generate_primary_header('a_filename.fits', p)
+    for name, value, *comment in header:
+        if name in test_data.keys():
+            assert value == test_data[name]
+
+
 def test_level1_processor_init():
     pro = FitsL1Processor('some/path')
     assert pro.archive_path == 'some/path'

stixcore/processing/find.py

@@ -171,7 +171,7 @@ def find_fits(args):
     logger.info(f'#candidates: {n_candidates}')
     logger.info(f'#found: {len(found)}')
 
-    for f in found:
+    for f in sorted(found):
         print(str(f))
         if args.copy_dest:
             fits_parent_path = f.relative_to(fits_dir)

stixcore/processing/tests/test_end2end.py

@@ -79,8 +79,9 @@ def test_identical(orig_fits, current_fits):
             error_c += 1
             warnings.warn(f"no corresponding file found for {cfits} in the original fits files")
             continue
-        diff = FITSDiff(ofits, cfits,
-                        ignore_keywords=['CHECKSUM', 'DATASUM', 'DATE', 'VERS_SW', 'HISTORY'])
+        diff = FITSDiff(ofits, cfits, atol=0.00001, rtol=0.00001,
+                        ignore_keywords=['CHECKSUM', 'DATASUM', 'DATE',
+                                         'VERS_SW', 'VERS_CFG', 'HISTORY'])
         if not diff.identical:
             error_c += 1
             warnings.warn(diff.report())

stixcore/processing/tests/test_publish.py

@@ -6,6 +6,7 @@
 import numpy as np
 import pytest
 
+import astropy.units as u
 from astropy.io import fits
 from astropy.table import QTable
 
@@ -138,6 +139,11 @@ 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.max_exposure = 3
+        product.dmin = 2
+        product.dmax = 3
+        product.bunit = "count"
         product.split_to_files.return_value = [product]
         product.get_energies = False
         product.get_processing_version.return_value = 2
@@ -203,6 +209,7 @@ def test_fits_incomplete_switch_over(out_dir):
             product.data = QTable({"time": t,
                                    "timedel": t-beg,
                                    "fcounts": np.array([1, 2]),
+                                   "counts": np.array([1, 2]) * u.deg_C,
                                    "control_index": [1, 1]})
             product.raw = ['packet1.xml', 'packet2.xml']
             product.parent = ['packet1.xml', 'packet2.xml']
@@ -218,6 +225,11 @@ def test_fits_incomplete_switch_over(out_dir):
             product.date_obs = beg
             product.date_beg = beg
             product.date_end = end
+            product.exposure = 2
+            product.max_exposure = 3
+            product.dmin = 2
+            product.dmax = 3
+            product.bunit = "counts"
             product.split_to_files.return_value = [product]
             product.get_energies = False
             product.get_processing_version.return_value = 2
@@ -256,6 +268,7 @@ def test_fits_incomplete_switch_over(out_dir):
             p = Product(get_complete_file_name_and_path(f))
             # old data.fcounts = [t1: 1, t2: 2]
             p.data['fcounts'][0] = 3
+            p.data['counts'] = 3 * u.deg_C
             # remove the second time stamp
             p.data.remove_row(1)
             # new data.fcounts = [t1: 3]
@@ -339,6 +352,11 @@ 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.max_exposure = 3
+    product.dmin = 2
+    product.dmax = 3
+    product.bunit = "count"
     product.split_to_files.return_value = [product]
     product.get_processing_version.return_value = 2
     product.get_energies = False

stixcore/products/level0/quicklookL0.py

@@ -18,7 +18,13 @@
     rebin_proportional,
     unscale_triggers,
 )
-from stixcore.products.product import Control, Data, EnergyChannelsMixin, GenericProduct
+from stixcore.products.product import (
+    Control,
+    CountDataMixin,
+    Data,
+    EnergyChannelsMixin,
+    GenericProduct,
+)
 from stixcore.time import SCETime, SCETimeDelta, SCETimeRange
 from stixcore.util.logging import get_logger
 
@@ -32,7 +38,7 @@
 QLNIX00405_off = 0.1
 
 
-class QLProduct(GenericProduct, EnergyChannelsMixin):
+class QLProduct(CountDataMixin, GenericProduct, EnergyChannelsMixin):
     """Generic QL product class composed of control and data."""
     def __init__(self, *, service_type, service_subtype, ssid, control, data,
                  idb_versions=defaultdict(SCETimeRange), **kwargs):
@@ -400,6 +406,14 @@ def from_levelb(cls, levelb, parent=''):
                    idb_versions=idb_versions,
                    packets=packets)
 
+    @property
+    def dmin(self):
+        return self.data['spectra'].min().value
+
+    @property
+    def dmax(self):
+        return self.data['spectra'].max().value
+
     @classmethod
     def _get_time(cls, control, num_energies, packets, pad_before, pad_after):
         times = []
@@ -494,6 +508,19 @@ def from_levelb(cls, levelb, parent=''):
                    idb_versions=idb_versions,
                    packets=packets)
 
+    @property
+    def dmin(self):
+        return self.data['variance'].min()
+
+    @property
+    def dmax(self):
+        return self.data['variance'].max()
+
+    @property
+    def bunit(self):
+        # TODO define
+        return ' '
+
     @classmethod
     def is_datasource_for(cls, *, service_type, service_subtype, ssid, **kwargs):
         return (kwargs['level'] == 'L0' and service_type == 21
@@ -554,6 +581,19 @@ def from_levelb(cls, levelb, parent=''):
                    idb_versions=idb_versions,
                    packets=packets)
 
+    @property
+    def dmin(self):
+        return min([self.data['loc_y'].min(), self.data['loc_z'].min()])
+
+    @property
+    def dmax(self):
+        return max([self.data['loc_y'].max(), self.data['loc_z'].max()])
+
+    @property
+    def bunit(self):
+        # TODO define
+        return ' '
+
     @classmethod
     def is_datasource_for(cls, *, service_type, service_subtype, ssid, **kwargs):
         return (kwargs['level'] == 'L0' and service_type == 21
@@ -784,6 +824,21 @@ def from_levelb(cls, levelb, parent=''):
                    idb_versions=idb_versions,
                    packets=packets)
 
+    @property
+    def dmin(self):
+        # TODO define
+        return 0.0
+
+    @property
+    def dmax(self):
+        # TODO define
+        return 0.0
+
+    @property
+    def bunit(self):
+        # TODO define
+        return ''
+
     @classmethod
     def is_datasource_for(cls, *, service_type, service_subtype, ssid, **kwargs):
         return (kwargs['level'] == 'L0' and service_type == 21

stixcore/products/level0/scienceL0.py

@@ -16,6 +16,7 @@
 )
 from stixcore.products.product import (
     ControlSci,
+    CountDataMixin,
     Data,
     EnergyChannelsMixin,
     FitsHeaderMixin,
@@ -71,7 +72,7 @@ class NotCombineException(Exception):
     pass
 
 
-class ScienceProduct(GenericProduct, EnergyChannelsMixin, FitsHeaderMixin):
+class ScienceProduct(CountDataMixin, GenericProduct, EnergyChannelsMixin, FitsHeaderMixin):
     """Generic science data product class composed of control and data."""
     def __init__(self, *, service_type, service_subtype, ssid, control, data, **kwargs):
         """Create a generic science data product composed of control and data.
@@ -682,6 +683,21 @@ def from_levelb(cls, levelb, parent=''):
         prod.add_additional_header_keywords(additional_header_keywords)
         return prod
 
+    @property
+    def dmin(self):
+        # TODO define columns for dmin/max
+        return 0.0
+
+    @property
+    def dmax(self):
+        # TODO define columns for dmin/max
+        return 0.0
+
+    @property
+    def bunit(self):
+        # TODO define columns for dmin/max
+        return ' '
+
     @classmethod
     def is_datasource_for(cls, *, service_type, service_subtype, ssid, **kwargs):
         return (kwargs['level'] == 'L0' and service_type == 21
@@ -938,6 +954,24 @@ def from_levelb(cls, levelb, parent=''):
                    idb_versions=idb_versions,
                    packets=packets)
 
+    @property
+    def dmin(self):
+        return min([self.data['cha_diode0'].min(),
+                    self.data['cha_diode1'].min(),
+                    self.data['chb_diode0'].min(),
+                    self.data['chb_diode1'].min()])
+
+    @property
+    def dmax(self):
+        return max([self.data['cha_diode0'].max(),
+                    self.data['cha_diode1'].max(),
+                    self.data['chb_diode0'].max(),
+                    self.data['chb_diode1'].max()])
+
+    @property
+    def bunit(self):
+        return ' '
+
     @classmethod
     def is_datasource_for(cls, *, service_type, service_subtype, ssid, **kwargs):
         return (kwargs['level'] == 'L0' and service_type == 21