From 12a643d030b71f194c5f85f92acf6a295b374145 Mon Sep 17 00:00:00 2001 From: Nicky Hochmuth Date: Wed, 17 Dec 2025 11:44:43 +0100 Subject: [PATCH 1/6] mixings in columns --- stixcore/io/FlareListManager.py | 45 ++-- .../io/product_processors/fits/processors.py | 20 +- stixcore/processing/FLtoFL.py | 15 +- stixcore/processing/FlareListL3.py | 2 +- stixcore/processing/pipeline_daily.py | 106 +++++---- stixcore/products/__init__.py | 1 + stixcore/products/level3/flarelist.py | 122 ++++++----- stixcore/products/level3/flarelistproduct.py | 14 +- stixcore/products/level3/processing.py | 203 ++++++++++++++++++ stixcore/products/product.py | 2 +- 10 files changed, 392 insertions(+), 138 deletions(-) create mode 100644 stixcore/products/level3/processing.py diff --git a/stixcore/io/FlareListManager.py b/stixcore/io/FlareListManager.py index 38359be6..e9d32e23 100644 --- a/stixcore/io/FlareListManager.py +++ b/stixcore/io/FlareListManager.py @@ -12,6 +12,7 @@ from astropy.time import Time from stixcore.config.config import CONFIG +from stixcore.io.product_processors.fits.processors import CreateUtcColumn from stixcore.products.level3.flarelist import FlarelistSC, FlarelistSDC from stixcore.products.product import Product from stixcore.util.logging import get_logger @@ -188,12 +189,17 @@ def get_data(self, *, start, end, fido_client): data["flare_id"] = Column( mt["flare_id"].astype(int), description=f"unique flare id for flarelist {self.flarelistname}" ) - data["start_UTC"] = Column(0, description="start time of flare") - data["start_UTC"] = [Time(d, format="isot", scale="utc") for d in mt["start_UTC"]] + CreateUtcColumn( + data, + [Time(d, format="isot", scale="utc") for d in mt["start_UTC"]], + "start_UTC", + description="start time of flare", + ) + data["duration"] = Column(mt["duration"].astype(float) * u.s, description="duration of flare") - data["end_UTC"] = Column(0, description="end time of flare") + data["end_UTC"] = CreateUtcColumn(description="end time of flare") data["end_UTC"] = [Time(d, format="isot", scale="utc") for d in mt["end_UTC"]] - data["peak_UTC"] = Column(0, description="flare peak time") + data["peak_UTC"] = CreateUtcColumn(description="flare peak time") data["peak_UTC"] = [Time(d, format="isot", scale="utc") for d in mt["peak_UTC"]] data["att_in"] = Column(mt["att_in"].astype(bool), description="was attenuator in during flare") data["bkg_baseline"] = Column(mt["LC0_BKG"] * u.ct, description="background baseline at 4-10 keV") @@ -277,7 +283,7 @@ def get_data(self, *, start, end, fido_client): data.add_index("flare_id") - # add energy axis for the lightcurve peek time data for each flare + # add energy axis for the lightcurve peak time data for each flare # the energy bins are taken from the daily ql-lightcurve products # as the definition of the lc energy chanel's are will change only very seldom # the ql-lightcurve products assume a constant definition for an entire day. @@ -439,13 +445,28 @@ def get_data(self, *, start, end, fido_client): data["flare_id"] = Column( mt["flare_id"].astype(int), description=f"unique flare id for flarelist {self.flarelistname}" ) - data["start_UTC"] = Column(0, description="start time of flare") - data["start_UTC"] = [Time(d, format="isot", scale="utc") for d in mt["start_UTC"]] + + CreateUtcColumn( + data, + [Time(d, format="isot", scale="utc") for d in mt["start_UTC"]], + "start_UTC", + description="start time of flare", + ) data["duration"] = Column(mt["duration"].astype(float) * u.s, description="duration of flare") - data["end_UTC"] = Column(0, description="end time of flare") - data["end_UTC"] = [Time(d, format="isot", scale="utc") for d in mt["end_UTC"]] - data["peak_UTC"] = Column(0, description="flare peak time") - data["peak_UTC"] = [Time(d, format="isot", scale="utc") for d in mt["peak_UTC"]] + + CreateUtcColumn( + data, + [Time(d, format="isot", scale="utc") for d in mt["end_UTC"]], + "end_UTC", + description="end time of flare", + ) + CreateUtcColumn( + data, + [Time(d, format="isot", scale="utc") for d in mt["peak_UTC"]], + "peak_UTC", + description="flare peak time", + ) + data["att_in"] = Column(mt["att_in"].astype(bool), description="was attenuator in during flare") data["bkg_baseline"] = Column(mt["LC0_BKG"] * u.ct, description="background baseline at 4-10 keV") data["GOES_class"] = Column( @@ -528,7 +549,7 @@ def get_data(self, *, start, end, fido_client): data.add_index("flare_id") - # add energy axis for the lightcurve peek time data for each flare + # add energy axis for the lightcurve peak time data for each flare # the energy bins are taken from the daily ql-lightcurve products # as the definition of the lc energy chanel's are will change only very seldom # the ql-lightcurve products assume a constant definition for an entire day. diff --git a/stixcore/io/product_processors/fits/processors.py b/stixcore/io/product_processors/fits/processors.py index 90630eea..7ca925b2 100644 --- a/stixcore/io/product_processors/fits/processors.py +++ b/stixcore/io/product_processors/fits/processors.py @@ -58,6 +58,24 @@ def set_bscale_unsigned(table_hdu): return table_hdu +def CreateUtcColumn(table, data, colname, description="UTC Time"): + """ + Create UTC time column for FITS tables. + + Parameters + ---------- + description : `str` + Description for the column + + Returns + ------- + `astropy.table.Column` + Column representing UTC time + """ + table[colname] = data + table[colname].info.description = description + + def add_default_tuint(table_hdu): """ Add a default empty string tunit if not already defined @@ -1139,7 +1157,7 @@ def write_fits(self, prod, *, version=0): # Add comment and history [primary_hdu.header.add_comment(com) for com in prod.comment] [primary_hdu.header.add_history(com) for com in prod.history] - primary_hdu.header.update({"HISTORY": "Processed by STIXCore ANC"}) + primary_hdu.header.update({"HISTORY": "Processed by STIXCore L3"}) if hasattr(prod, "maps") and len(prod.maps) > 0: # fig = plt.figure(figsize=(12, 6)) diff --git a/stixcore/processing/FLtoFL.py b/stixcore/processing/FLtoFL.py index 087f0d36..7f08acf1 100644 --- a/stixcore/processing/FLtoFL.py +++ b/stixcore/processing/FLtoFL.py @@ -15,7 +15,7 @@ ) from stixcore.products.level3.flarelist import ( FlareList, - FlarePeekPreviewMixin, + FlarePeakPreviewMixin, FlarePositionMixin, FlareSOOPMixin, ) @@ -107,7 +107,7 @@ def test_for_processing( """ try: c_header = fits.getheader(candidate) - f_data_end = datetime.fromisoformat(c_header["DATE-END"]) + # f_data_end = datetime.fromisoformat(c_header["DATE-END"]) f_create_date = datetime.fromisoformat(c_header["DATE"]) cfn = get_complete_file_name_and_path(candidate) @@ -128,8 +128,9 @@ def test_for_processing( # safety margin of 1day until we process higher products with position and pointing # only use flown spice kernels not predicted once as pointing information # can be "very off" - if f_data_end > (Spice.instance.get_mk_date(meta_kernel_type="flown") - timedelta(hours=24)): - return TestForProcessingResult.NotSuitable + # TODO redo + # if f_data_end > (Spice.instance.get_mk_date(meta_kernel_type="flown") - timedelta(hours=24)): + # return TestForProcessingResult.NotSuitable # safety margin of x until we start with processing the list files if f_create_date >= (datetime.now() - self.cadence): @@ -187,9 +188,9 @@ def process_fits_files( if issubclass(out_product, FlareSOOPMixin) and not issubclass(in_product, FlareSOOPMixin): out_product.add_soop(data) - # add peek preview images if not already present - if issubclass(out_product, FlarePeekPreviewMixin) and not issubclass(in_product, FlarePeekPreviewMixin): - out_product.add_peek_preview(data, energy, file_path.name, fido_client, img_processor, month=month) + # add peak preview images if not already present + if issubclass(out_product, FlarePeakPreviewMixin) and not issubclass(in_product, FlarePeakPreviewMixin): + out_product.add_peak_preview(data, energy, file_path.name, fido_client, img_processor, month=month) out_prod = out_product(control=control, data=data, month=month, energy=energy) out_prod.parent = file_path.name diff --git a/stixcore/processing/FlareListL3.py b/stixcore/processing/FlareListL3.py index 89e6d5c8..855bdde2 100644 --- a/stixcore/processing/FlareListL3.py +++ b/stixcore/processing/FlareListL3.py @@ -24,7 +24,7 @@ class FlareListL3(SingleProductProcessingStepMixin): """Processing step from a FlareListManager to monthly solo_L3_stix-flarelist-*.fits file.""" - STARTDATE = date(2024, 1, 1) + STARTDATE = date(2025, 1, 1) def __init__(self, flm: FlareListManager, output_dir: Path): """Crates a new Processor. diff --git a/stixcore/processing/pipeline_daily.py b/stixcore/processing/pipeline_daily.py index 09b19181..0b1a5a20 100644 --- a/stixcore/processing/pipeline_daily.py +++ b/stixcore/processing/pipeline_daily.py @@ -9,34 +9,27 @@ from stixcore.config.config import CONFIG from stixcore.ephemeris.manager import Spice, SpiceKernelManager +from stixcore.io.FlareListManager import SDCFlareListManager from stixcore.io.ProcessingHistoryStorage import ProcessingHistoryStorage from stixcore.io.product_processors.fits.processors import ( # FitsANCProcessor,; FitsL3Processor, + FitsANCProcessor, FitsL2Processor, + FitsL3Processor, ) from stixcore.io.product_processors.plots.processors import PlotProcessor from stixcore.io.RidLutManager import RidLutManager from stixcore.processing.AspectANC import AspectANC +from stixcore.processing.FlareListL3 import FlareListL3 +from stixcore.processing.FLtoFL import FLtoFL from stixcore.processing.LL import LL03QL from stixcore.processing.pipeline import PipelineStatus from stixcore.processing.SingleStep import SingleProcessingStepResult from stixcore.products.level1.quicklookL1 import LightCurve +from stixcore.products.level3.flarelist import FlarelistSDC, FlarelistSDCLoc from stixcore.products.lowlatency.quicklookLL import LightCurveL3 from stixcore.soop.manager import SOOPManager from stixcore.util.logging import STX_LOGGER_DATE_FORMAT, STX_LOGGER_FORMAT, get_logger -# from stixpy.net.client import STIXClient -# from stixcore.io.FlareListManager import SCFlareListManager, SDCFlareListManager -# from stixcore.processing.FlareListL3 import FlareListL3 -# from stixcore.processing.FLtoFL import FLtoFL -# from stixcore.products.level3.flarelist import ( -# FlarelistSC, -# FlarelistSCLoc, -# FlarelistSCLocImg, -# FlarelistSDC, -# FlarelistSDCLoc, -# FlarelistSDCLocImg, -# ) - logger = get_logger(__name__) @@ -215,8 +208,8 @@ def run_daily_pipeline(args): # SCFlareListManager.instance = SCFlareListManager(flare_lut_file, fido_client, update=True) # TODO reactivate once flarelist processing is finalized - # flare_lut_file = Path(CONFIG.get("Pipeline", "flareid_sdc_lut_file")) - # SDCFlareListManager.instance = SDCFlareListManager(flare_lut_file, update=False) + flare_lut_file = Path(CONFIG.get("Pipeline", "flareid_sdc_lut_file")) + SDCFlareListManager.instance = SDCFlareListManager(flare_lut_file, update=False) RidLutManager.instance = RidLutManager(Path(CONFIG.get("Publish", "rid_lut_file")), update=False) @@ -249,19 +242,19 @@ def run_daily_pipeline(args): aspect_anc_processor = AspectANC(fits_in_dir, fits_out_dir) # TODO reactivate once flarelist processing is finalized - # flarelist_sdc = FlareListL3(SDCFlareListManager.instance, fits_out_dir) + flarelist_sdc = FlareListL3(SDCFlareListManager.instance, fits_out_dir) # flarelist_sc = FlareListL3(SCFlareListManager.instance, fits_out_dir) - # fl_to_fl = FLtoFL( - # fits_in_dir, - # fits_out_dir, - # products_in_out=[ - # (FlarelistSDC, FlarelistSDCLoc), - # (FlarelistSDCLoc, FlarelistSDCLocImg), - # (FlarelistSC, FlarelistSCLoc), - # (FlarelistSCLoc, FlarelistSCLocImg), - # ], - # cadence=timedelta(seconds=1), - # ) + fl_to_fl = FLtoFL( + fits_in_dir, + fits_out_dir, + products_in_out=[ + (FlarelistSDC, FlarelistSDCLoc), + # (FlarelistSDCLoc, FlarelistSDCLocImg), + # (FlarelistSC, FlarelistSCLoc), + # (FlarelistSCLoc, FlarelistSCLocImg), + ], + cadence=timedelta(seconds=1), + ) ll03ql = LL03QL( fits_in_dir, fits_out_dir, in_product=LightCurve, out_product=LightCurveL3, cadence=timedelta(seconds=1) @@ -270,24 +263,25 @@ def run_daily_pipeline(args): plot_writer = PlotProcessor(fits_out_dir) l2_fits_writer = FitsL2Processor(fits_out_dir) # TODO reactivate once flarelist processing is finalized - # l3_fits_writer = FitsL3Processor(fits_out_dir) - # anc_fits_writer = FitsANCProcessor(fits_out_dir) + l3_fits_writer = FitsL3Processor(fits_out_dir) + anc_fits_writer = FitsANCProcessor(fits_out_dir) - hk_in_files = aspect_anc_processor.get_processing_files(phs) + # hk_in_files = aspect_anc_processor.get_processing_files(phs) + hk_in_files = [] - ll_candidates = ll03ql.get_processing_files(phs) - # ll_candidates = [] + # ll_candidates = ll03ql.get_processing_files(phs) + ll_candidates = [] # TODO reactivate once flarelist processing is finalized - # fl_sdc_months = flarelist_sdc.find_processing_months(phs) - # fl_sdc_months = [] + fl_sdc_months = flarelist_sdc.find_processing_months(phs) + fl_sdc_months = [] # TODO reactivate once flarelist processing is finalized # fl_sc_months = flarelist_sc.find_processing_months(phs) # fl_sc_months = [] # TODO reactivate once flarelist processing is finalized - # fl_to_fl_files = fl_to_fl.get_processing_files(phs) + fl_to_fl_files = fl_to_fl.get_processing_files(phs) # fl_to_fl_files = [] # all processing files should be terminated before the next step as the different @@ -306,16 +300,16 @@ def run_daily_pipeline(args): ) ) # TODO reactivate once flarelist processing is finalized - # jobs.append( - # executor.submit( - # flarelist_sdc.process_fits_files, - # fl_sdc_months, - # soopmanager=SOOPManager.instance, - # spice_kernel_path=Spice.instance.meta_kernel_path, - # processor=l3_fits_writer, - # config=CONFIG, - # ) - # ) + jobs.append( + executor.submit( + flarelist_sdc.process_fits_files, + fl_sdc_months, + soopmanager=SOOPManager.instance, + spice_kernel_path=Spice.instance.meta_kernel_path, + processor=l3_fits_writer, + config=CONFIG, + ) + ) # jobs.append( # executor.submit( @@ -330,17 +324,17 @@ def run_daily_pipeline(args): # # TODO a owen processing step for each flarelist file? # # for fl_to_fl_file in fl_to_fl_files: - # jobs.append( - # executor.submit( - # fl_to_fl.process_fits_files, - # fl_to_fl_files, - # soopmanager=SOOPManager.instance, - # spice_kernel_path=Spice.instance.meta_kernel_path, - # fl_processor=anc_fits_writer, - # img_processor=l3_fits_writer, - # config=CONFIG, - # ) - # ) + jobs.append( + executor.submit( + fl_to_fl.process_fits_files, + fl_to_fl_files, + soopmanager=SOOPManager.instance, + spice_kernel_path=Spice.instance.meta_kernel_path, + fl_processor=anc_fits_writer, + img_processor=l3_fits_writer, + config=CONFIG, + ) + ) jobs.append( executor.submit( diff --git a/stixcore/products/__init__.py b/stixcore/products/__init__.py index 7863b1bf..974838ed 100644 --- a/stixcore/products/__init__.py +++ b/stixcore/products/__init__.py @@ -10,5 +10,6 @@ from stixcore.products.level2.housekeepingL2 import * from stixcore.products.level2.quicklookL2 import * from stixcore.products.level2.scienceL2 import * +from stixcore.products.level3.flarelist import * from stixcore.products.levelb.binary import LevelB from stixcore.products.lowlatency.quicklookLL import * diff --git a/stixcore/products/level3/flarelist.py b/stixcore/products/level3/flarelist.py index 188f6193..80e30601 100644 --- a/stixcore/products/level3/flarelist.py +++ b/stixcore/products/level3/flarelist.py @@ -24,7 +24,8 @@ from stixcore.config.config import CONFIG from stixcore.ephemeris.manager import Spice -from stixcore.products.level3.flarelistproduct import PeekPreviewImage +from stixcore.products.level3.flarelistproduct import PeakPreviewImage +from stixcore.products.level3.processing import estimate_stix_flare_location from stixcore.products.product import CountDataMixin, GenericProduct, L2Mixin, read_qtable from stixcore.soop.manager import SOOPManager from stixcore.time import SCETime, SCETimeRange @@ -39,7 +40,7 @@ "FlareSOOPMixin", "FlareList", "FlarelistSDCLocImg", - "FlarePeekPreviewMixin", + "FlarePeakPreviewMixin", "FlarelistSC", "FlarelistSCLoc", "FlarelistSCLocImg", @@ -85,13 +86,16 @@ def add_flare_position( fido_client: STIXClient, *, filter_function=lambda x: True, - peek_time_colname="peak_UTC", + peak_time_colname="peak_UTC", start_time_colname="start_UTC", end_time_colname="end_UTC", keep_all_flares=True, month=None, ): - data["flare_position"] = [SkyCoord(0, 0, frame="icrs", unit="deg") for i in range(0, len(data))] + # helio_frame = Helioprojective(observer="earth") + # SkyCoord(HeliographicStonyhurst(0 * u.deg, 0 * u.deg)) + # SkyCoord(0 * u.deg, 0 * u.deg, frame=helio_frame) + data["flare_position"] = [SkyCoord(HeliographicStonyhurst(0 * u.deg, 0 * u.deg)) for i in range(0, len(data))] data["anc_ephemeris_path"] = Column(" " * 500, dtype=str, description="TDB") data["cpd_path"] = Column(" " * 500, dtype=str, description="TDB") @@ -106,11 +110,11 @@ def add_flare_position( total_flares = len(data) day_asp_ephemeris_cache = dict() - + flare_positions = [] for i, row in enumerate(data): - if filter_function(row): + if filter_function(row) and i < 200: pass_filter += 1 - peak_time = row[peek_time_colname] + peak_time = row[peak_time_colname] start_time = row[start_time_colname] end_time = row[end_time_colname] @@ -187,15 +191,26 @@ def add_flare_position( best_cpd_idx = 0 data[i]["cpd_path"] = cpd_res["path"][best_cpd_idx] - # do the calculations with stixpy + try: + stixpy_cpd = STIXPYProduct(Path(data[i]["cpd_path"])) + coord, map = estimate_stix_flare_location(stixpy_cpd) + + roll, solo_xyz, pointing = get_hpc_info(start_time, end_time) + solo = HeliographicStonyhurst(*solo_xyz, obstime=peak_time, representation_type="cartesian") - data[i]["flare_position"] = SkyCoord(1, 1, frame="icrs", unit="deg") - data[i]["_position_status"] = True - data[i]["_position_message"] = "OK" + # data[i]["flare_position"] = coord.transform_to(Helioprojective(observer=solo)) + flare_positions.append(coord.transform_to(Helioprojective(observer=solo))) + data[i]["_position_status"] = True + data[i]["_position_message"] = "OK" + except Exception as e: + flare_positions.append(None) + data[i]["_position_message"] = f"Error: {type(e)}" else: to_remove.append(i) + flare_positions.append(None) + data["flare_position"] = flare_positions if not keep_all_flares: data.remove_rows(to_remove) @@ -203,7 +218,8 @@ def add_flare_position( f"Flare position calculated for month {month} with {total_flares} flares, " f"passed filter: {pass_filter} no ephemeris data found for {no_ephemeris} " f"flares, no CPD data found for {no_cpd} flares, many CPD data found for " - f"{many_cpd} flares, one CPD data found for {one_cpd} flares" + f"{many_cpd} flares, one CPD data found for {one_cpd} flares." + f"finally {len(data)} flares remaining" ) @@ -212,14 +228,14 @@ class FlareSOOPMixin: @classmethod def add_soop( - self, data, *, peek_time_colname="peak_UTC", start_time_colname="start_UTC", end_time_colname="end_UTC" + self, data, *, peak_time_colname="peak_UTC", start_time_colname="start_UTC", end_time_colname="end_UTC" ): soop_encoded_type = list() soop_id = list() soop_type = list() for row in data: - soops = SOOPManager.instance.find_soops(start=row[peek_time_colname]) + soops = SOOPManager.instance.find_soops(start=row[peak_time_colname]) if soops: soop = soops[0] soop_encoded_type.append(soop.encodedSoopType) @@ -235,13 +251,13 @@ def add_soop( data["soop_type"] = Column(soop_type, dtype=str, description="name of the SOOP campaign") -class FlarePeekPreviewMixin: - """Mixin class to add peek preview images to flare list products. - This class provides a method to generate and add peek preview images +class FlarePeakPreviewMixin: + """Mixin class to add peak preview images to flare list products. + This class provides a method to generate and add peak preview images to the flare list data. The images are generated based on the flare's peak time, start time, and end time, using the STIXPy library for visibility calculations and image reconstruction. - The generated images are stored in the 'peek_preview_path' column of the data. + The generated images are stored in the 'peak_preview_path' column of the data. The method also updates the status and message columns to indicate the success or failure of the image generation process. @@ -249,7 +265,7 @@ class FlarePeekPreviewMixin: """ @classmethod - def add_peek_preview( + def add_peak_preview( cls, data, energies, @@ -257,7 +273,7 @@ def add_peek_preview( fido_client: STIXClient, img_processor, *, - peek_time_colname="peak_UTC", + peak_time_colname="peak_UTC", start_time_colname="start_UTC", end_time_colname="end_UTC", anc_ephemeris_path_colname="anc_ephemeris_path", @@ -266,17 +282,17 @@ def add_peek_preview( keep_all_flares=True, month=None, ): - data["peek_preview_path"] = Column(" " * 500, dtype=str, description="TDB") - data["preview_start_UTC"] = [Time(d, format="isot", scale="utc") for d in data[peek_time_colname]] - data["preview_end_UTC"] = [Time(d, format="isot", scale="utc") for d in data[peek_time_colname]] - data["_peek_preview_status"] = Column(False, dtype=bool, description="TDB") - data["_peek_preview_message"] = Column(" " * 500, dtype=str, description="TDB") + data["peak_preview_path"] = Column(" " * 500, dtype=str, description="TDB") + data["preview_start_UTC"] = [Time(d, format="isot", scale="utc") for d in data[peak_time_colname]] + data["preview_end_UTC"] = [Time(d, format="isot", scale="utc") for d in data[peak_time_colname]] + data["_peak_preview_status"] = Column(False, dtype=bool, description="TDB") + data["_peak_preview_message"] = Column(" " * 500, dtype=str, description="TDB") to_remove = [] products = [] images = 0 for i, row in enumerate(data): - peak_time = row[peek_time_colname] + peak_time = row[peak_time_colname] row[start_time_colname] row[end_time_colname] @@ -286,8 +302,8 @@ def add_peek_preview( status = False message = "" - peek_preview_start = row[peek_time_colname] - peek_preview_end = row[peek_time_colname] + peak_preview_start = row[peak_time_colname] + peak_preview_end = row[peak_time_colname] if anc_ephemeris_path.exists() and cpd_path.exists(): try: @@ -295,18 +311,18 @@ def add_peek_preview( # do the imaging with stixpy preview_data = data[i : i + 1] - del preview_data["peek_preview_path"] - del preview_data["_peek_preview_status"] - del preview_data["_peek_preview_message"] + del preview_data["peak_preview_path"] + del preview_data["_peak_preview_status"] + del preview_data["_peak_preview_message"] - peek_preview_start = row[peek_time_colname] - 10 * u.s - peek_preview_end = row[peek_time_colname] + 10 * u.s + peak_preview_start = row[peak_time_colname] - 10 * u.s + peak_preview_end = row[peak_time_colname] + 10 * u.s - preview_data["preview_start_UTC"] = peek_preview_start - preview_data["preview_end_UTC"] = peek_preview_end + preview_data["preview_start_UTC"] = peak_preview_start + preview_data["preview_end_UTC"] = peak_preview_end cpd_sci = STIXPYProduct(cpd_path) - time_range_sci = [peek_preview_start, peek_preview_end] + time_range_sci = [peak_preview_start, peak_preview_end] maps = [] for energy_range in [[4, 20], [20, 120]] * u.keV: # flare_position = preview_data['flare_position'][0] @@ -388,7 +404,7 @@ def add_peek_preview( maps.append((map_with_erange, header)) - ppi = PeekPreviewImage( + ppi = PeakPreviewImage( control=QTable(), data=preview_data, month=month, @@ -407,18 +423,18 @@ def add_peek_preview( status = False message = str(e) - data[i]["preview_start_UTC"] = peek_preview_start - data[i]["preview_end_UTC"] = peek_preview_end - data[i]["peek_preview_path"] = "test" - data[i]["_peek_preview_status"] = status - data[i]["_peek_preview_message"] = message + data[i]["preview_start_UTC"] = peak_preview_start + data[i]["preview_end_UTC"] = peak_preview_end + data[i]["peak_preview_path"] = "test" + data[i]["_peak_preview_status"] = status + data[i]["_peak_preview_message"] = message if not keep_all_flares: data.remove_rows(to_remove) logger.info( f"Flare images created for month {month} with {len(data)} flares, " - f"{len(products)} peek previews created, with total {images} images" + f"{len(products)} peak previews created, with total {images} images" ) return products @@ -558,7 +574,7 @@ def add_flare_position(cls, data, fido_client: STIXClient, *, month=None): data, fido_client, filter_function=cls.filter_flare_function, - peek_time_colname="peak_UTC", + peak_time_colname="peak_UTC", start_time_colname="start_UTC", end_time_colname="end_UTC", keep_all_flares=False, @@ -570,7 +586,7 @@ def is_datasource_for(cls, *, service_type, service_subtype, ssid, **kwargs): return kwargs["level"] == "L3" and service_type == 0 and service_subtype == 0 and ssid == 3 -class FlarelistSDCLocImg(FlarelistSDCLoc, FlarePeekPreviewMixin): +class FlarelistSDCLocImg(FlarelistSDCLoc, FlarePeakPreviewMixin): """Flarelist product class for StixDataCenter flares. In ANC product format. @@ -589,14 +605,14 @@ def enhance_from_product(self, in_prod: GenericProduct): pass @classmethod - def add_peek_preview(cls, data, energies, parent, fido_client: STIXClient, img_processor, *, month=None): - super().add_peek_preview( + def add_peak_preview(cls, data, energies, parent, fido_client: STIXClient, img_processor, *, month=None): + super().add_peak_preview( data, energies, parent, fido_client, img_processor, - peek_time_colname="peak_UTC", + peak_time_colname="peak_UTC", start_time_colname="start_UTC", end_time_colname="end_UTC", anc_ephemeris_path_colname="anc_ephemeris_path", @@ -695,7 +711,7 @@ def add_flare_position(cls, data, fido_client: STIXClient, *, month=None): data, fido_client, filter_function=cls.filter_flare_function, - peek_time_colname="peak_UTC", + peak_time_colname="peak_UTC", start_time_colname="start_UTC", end_time_colname="end_UTC", keep_all_flares=False, @@ -707,7 +723,7 @@ def is_datasource_for(cls, *, service_type, service_subtype, ssid, **kwargs): return kwargs["level"] == "L3" and service_type == 0 and service_subtype == 0 and ssid == 7 -class FlarelistSCLocImg(FlarelistSCLoc, FlarePeekPreviewMixin): +class FlarelistSCLocImg(FlarelistSCLoc, FlarePeakPreviewMixin): """Flarelist product class for StixCore flares. In ANC product format. @@ -726,14 +742,14 @@ def enhance_from_product(self, in_prod: GenericProduct): pass @classmethod - def add_peek_preview(cls, data, energies, parent, fido_client: STIXClient, img_processor, *, month=None): - super().add_peek_preview( + def add_peak_preview(cls, data, energies, parent, fido_client: STIXClient, img_processor, *, month=None): + super().add_peak_preview( data, energies, parent, fido_client, img_processor, - peek_time_colname="peak_UTC", + peak_time_colname="peak_UTC", start_time_colname="start_UTC", end_time_colname="end_UTC", anc_ephemeris_path_colname="anc_ephemeris_path", diff --git a/stixcore/products/level3/flarelistproduct.py b/stixcore/products/level3/flarelistproduct.py index ae7cb0d4..032864e1 100644 --- a/stixcore/products/level3/flarelistproduct.py +++ b/stixcore/products/level3/flarelistproduct.py @@ -6,7 +6,7 @@ from stixcore.time.datetime import SCETime, SCETimeRange from stixcore.util.logging import get_logger -__all__ = ["FlareListProduct", "PeekPreviewImage"] +__all__ = ["FlareListProduct", "PeakPreviewImage"] logger = get_logger(__name__) @@ -22,17 +22,17 @@ def from_timerange(cls, timerange: SCETimeRange, *, flarelistparent: str = ""): pass -class PeekPreviewImage(FlareListProduct): +class PeakPreviewImage(FlareListProduct): PRODUCT_PROCESSING_VERSION = 1 Level = "L3" Type = "sci" - Name = "peekpreviewimg" + Name = "peakpreviewimg" def __init__(self, control, data, energy, maps, parents, *, product_name_suffix="", **kwargs): super().__init__(service_type=0, service_subtype=0, ssid=5, control=control, data=data, energy=energy, **kwargs) - self.name = f"{PeekPreviewImage.Name}-{product_name_suffix}" - self.level = PeekPreviewImage.Level - self.type = PeekPreviewImage.Type + self.name = f"{PeakPreviewImage.Name}-{product_name_suffix}" + self.level = PeakPreviewImage.Level + self.type = PeakPreviewImage.Type self.energy = energy self.maps = maps self.parents = parents @@ -62,4 +62,4 @@ def split_to_files(self): @classmethod def is_datasource_for(cls, *, service_type, service_subtype, ssid, **kwargs): - return kwargs["level"] == PeekPreviewImage.Level and service_type == 0 and service_subtype == 0 and ssid == 5 + return kwargs["level"] == PeakPreviewImage.Level and service_type == 0 and service_subtype == 0 and ssid == 5 diff --git a/stixcore/products/level3/processing.py b/stixcore/products/level3/processing.py new file mode 100644 index 00000000..e6225a1e --- /dev/null +++ b/stixcore/products/level3/processing.py @@ -0,0 +1,203 @@ +######################################################### +### Temporary code should be come from stixpy finally ### +######################################################### + +import numpy as np +import stixpy.calibration.visibility +import stixpy.coordinates.transforms +import sunpy.map +import sunpy.time +import xrayvision.imaging +from stixpy.coordinates.transforms import STIXImaging +from sunpy.coordinates import HeliographicStonyhurst + +import astropy.units as u +from astropy.coordinates import SkyCoord +from astropy.time import Time + + +def construct_stix_calibrated_visibilities( + cpd_sci, + flare_location, + time_range=None, + energy_range=None, + subcollimators=None, + cpd_bkg=None, + time_range_bkg=None, + **kwargs, +): + """ + Constructs calibrated STIX visibilities from STIX compressed pixel data. + + Extra kwargs are passed to `stixpy.calibration.visibility.create_meta_pixels` + + Parameters + ---------- + cpd_sci: `stixpy.product.Product` + The STIX pixel data. Assumed to be already background subtracted. + flare_location: `astropy.coordinates.SkyCoord` + The flare location. Frame must be convertible to `stixpy.coordinates.transdforms.STIXImaging`. + time_range: `sunpy.time.TimeRange` (optional) + The time range over which to estimate the flare location. + Default is all times in cpd_sci. + energy_range: `astropy.units.Quantity` in spectral units (optional) + Length-2 quantity giving the lower and upper bounds of the energy range to use for imaging. + Default is all finite energies in cpd_sci. + subcollimators: `iterable` of `str` + The labels of the subcollimators to use in estimating the flare locations, e.g. + ``['10a', '10b', '10c',...]`` + Default is all subcollimators in cpd_sci. + cpd_bkg: stixpy.product.Product` (optional) + The background to subtract from the pixel data before determining the flare location. + If None and time_range_bkg is also None, no background is subtracted. + time_range_bkg: `sunpy.time.TimeRange` + The time range within cpd_bkg to use for the background subtraction. + If None, entire time range of cpd_bkg is used to determine background. + If not None, and cpd_bkg is None, the background is determined from this time range + applied to cpd_sci. + + Returns + ------- + vis: `xrayvision.visibility.Visibilities` + The calibrated STIX visibilities. + """ + # Sanitze inputs. + if time_range is None: + time_range = cpd_sci.time_range + times = Time([time_range.start, time_range.end]) + if energy_range is None: + energy_range = u.Quantity( + [ + cpd_sci.energies["e_low"][np.isfinite(cpd_sci.energies["e_low"])][0], + cpd_sci.energies["e_high"][np.isfinite(cpd_sci.energies["e_high"])][-1], + ] + ) + no_shadowing = kwargs.pop("no_shadowing", True) + # Generate meta_pixels from pixel_data. + meta_pixels = stixpy.calibration.visibility.create_meta_pixels( + cpd_sci, + time_range=times, + energy_range=energy_range, + flare_location=flare_location, + no_shadowing=no_shadowing, + **kwargs, + ) + # Subtract background if background pixel data provided. + if cpd_bkg is None and time_range_bkg is not None: + cpd_bkg = cpd_sci + if cpd_bkg is not None: + if time_range_bkg is None: + time_range_bkg = cpd_bkg.time_range + times_bkg = Time([time_range_bkg.start, time_range_bkg.end]) + meta_pixels_bkg = stixpy.calibration.visibility.create_meta_pixels( + cpd_bkg, + time_range=times_bkg, + energy_range=energy_range, + flare_location=[0, 0] * u.arcsec, + no_shadowing=no_shadowing, + **kwargs, + ) + meta_pixels = _subtract_background_from_stix_meta_pixels(meta_pixels, meta_pixels_bkg) + # Generate and calibrate visibilities. + vis = stixpy.calibration.visibility.create_visibility(meta_pixels) + vis = stixpy.calibration.visibility.calibrate_visibility(vis, flare_location=SkyCoord(flare_location)) + if subcollimators is not None: + idx_subcol = np.argwhere(np.isin(vis.meta["vis_labels"], subcollimators)).ravel() + vis = vis[idx_subcol] + return vis + + +def estimate_stix_flare_location( + cpd_sci, time_range=None, energy_range=None, subcollimators=None, cpd_bkg=None, time_range_bkg=None +): + """ + Estimates flare location from STIX compressed pixel data. + + FLare location is assumed to be the location of the brightest pixel in the backprojection + map calculated from the input STIX pixel data. + + Parameters + ---------- + cpd_sci: `stixpy.product.Product` + The STIX pixel data. Assumed to be already background subtracted. + time_range: `sunpy.time.TimeRange` (optional) + The time range over which to estimate the flare location. + Default is all times in cpd_sci. + energy_range: `astropy.units.Quantity` in spectral units (optional) + Length-2 quantity giving the lower and upper bounds of the energy range to use for imaging. + Default is all finite energies in cpd_sci. + subcollimators: `iterable` of `str` (optional) + The labels of the subcollimators to include in the output Visibilities object, e.g. + ``['10a', '10b', '10c',...]`` + Default is all subcollimators in cpd_sci + cpd_bkg: stixpy.product.Product` (optional) + The background to subtract from the pixel data before determining the flare location. + Passed to construct_stix_calibrated_visibilities(). + time_range_bkg: `sunpy.time.TimeRange` + The time range within cpd_bkg to use for the background subtraction. + Passed to construct_stix_calibrated_visibilities(). + + Returns + ------- + flare_loc: `astropy.coordinates.SkyCoord` + The estimated flare location. + map_bp: `sunpy.map.Map` + The backprojection map from whose brightest pixel the flare location was estimated. + """ + if time_range is None: + time_range = cpd_sci.time_range + if subcollimators is None: + subcollimators = ["10a", "10b", "10c", "9a", "9b", "9c", "8a", "8b", "8c", "7a", "7b", "7c"] + # Construct STIX location and centre of FOV. + roll, solo_xyz, pointing = stixpy.coordinates.transforms.get_hpc_info(time_range.start, time_range.end) + solo = HeliographicStonyhurst(*solo_xyz, obstime=time_range.center, representation_type="cartesian") + fov_centre = STIXImaging( + 0 * u.arcsec, 0 * u.arcsec, obstime=time_range.start, obstime_end=time_range.end, observer=solo + ) + # Generate calibrated visibilities using coarser subcollimators. + vis = construct_stix_calibrated_visibilities( + cpd_sci, fov_centre, time_range=time_range, energy_range=energy_range, subcollimators=subcollimators + ) + # Produce backprojected image and find brightest pixel. Use this for flare location. + imsize = [512, 512] * u.pixel # number of pixels of the map to reconstruct + plate_scale = [10, 10] * u.arcsec / u.pixel # pixel size in arcsec + bp_image = xrayvision.imaging.vis_to_image(vis, imsize, pixel_size=plate_scale) + max_idx = np.argwhere(bp_image == bp_image.max()).ravel() + # Calculate WCS for backprojected image in STIXImaging and HPC frames. + # Recalculate STIX HPC info as slightly different times will have been used + # than input times due to onboard STIX time binning. + vis_tr = sunpy.time.TimeRange(vis.meta["time_range"]) + roll, solo_xyz, pointing = stixpy.coordinates.transforms.get_hpc_info(vis_tr.start, vis_tr.end) + solo = HeliographicStonyhurst(*solo_xyz, obstime=vis_tr.center, representation_type="cartesian") + coord = STIXImaging(0 * u.arcsec, 0 * u.arcsec, obstime=vis_tr.start, obstime_end=vis_tr.end, observer=solo) + header_bp = sunpy.map.make_fitswcs_header( + bp_image, coord, telescope="STIX", observatory="Solar Orbiter", scale=plate_scale + ) + map_bp = sunpy.map.Map(bp_image, header_bp) + wcs_bp = map_bp.wcs + # Estimate flare location from brightest pixel in backprojection image + flare_loc = wcs_bp.array_index_to_world(*max_idx) + return flare_loc, map_bp + + +def _subtract_background_from_stix_meta_pixels(meta_pixels_sci, meta_pixels_bkg): + """ + Estimates flare location from STIX pixel data. + + Parameters + ---------- + meta_pixels_sci: `dict` + The STIX meta pixel representing the observations. Format must be + same as output from `stixpy.calibration.visibility.create_meta_pixels`. + meta_pixels_bkg: `dict` + The STIX meta pixels representing the background. Format must be + same as output from `stixpy.calibration.visibility.create_meta_pixels`. + """ + meta_pixels_bkg_sub = { + **meta_pixels_sci, + "abcd_rate_kev_cm": meta_pixels_sci["abcd_rate_kev_cm"] - meta_pixels_bkg["abcd_rate_kev_cm"], + "abcd_rate_error_kev_cm": np.sqrt( + meta_pixels_sci["abcd_rate_error_kev_cm"] ** 2 + meta_pixels_bkg["abcd_rate_error_kev_cm"] ** 2 + ), + } + return meta_pixels_bkg_sub diff --git a/stixcore/products/product.py b/stixcore/products/product.py index 02712266..27415375 100644 --- a/stixcore/products/product.py +++ b/stixcore/products/product.py @@ -74,7 +74,7 @@ def read_qtable(file, hdu, hdul=None): `astropy.table.QTable` The corrected QTable with correct data types """ - qtable = QTable.read(file, hdu) + qtable = QTable.read(file, hdu, astropy_native=True) if hdul is None: hdul = fits.open(file) From 0b5f6c4a1c9904d434d95e5e190c6b783cabc0c2 Mon Sep 17 00:00:00 2001 From: Nicky Hochmuth Date: Mon, 26 Jan 2026 17:57:32 +0100 Subject: [PATCH 2/6] skyCoords components as single columns --- stixcore/processing/FlareListL3.py | 2 +- stixcore/products/level3/flarelist.py | 119 +++++++++++++++++++++++--- stixcore/soop/manager.py | 2 + 3 files changed, 110 insertions(+), 13 deletions(-) diff --git a/stixcore/processing/FlareListL3.py b/stixcore/processing/FlareListL3.py index 855bdde2..5a6e8b3b 100644 --- a/stixcore/processing/FlareListL3.py +++ b/stixcore/processing/FlareListL3.py @@ -24,7 +24,7 @@ class FlareListL3(SingleProductProcessingStepMixin): """Processing step from a FlareListManager to monthly solo_L3_stix-flarelist-*.fits file.""" - STARTDATE = date(2025, 1, 1) + STARTDATE = date(2025, 7, 1) def __init__(self, flm: FlareListManager, output_dir: Path): """Crates a new Processor. diff --git a/stixcore/products/level3/flarelist.py b/stixcore/products/level3/flarelist.py index 80e30601..9379f02f 100644 --- a/stixcore/products/level3/flarelist.py +++ b/stixcore/products/level3/flarelist.py @@ -1,5 +1,6 @@ from pathlib import Path from datetime import datetime +from itertools import groupby import numpy as np from stixpy.calibration.visibility import ( @@ -10,7 +11,7 @@ from stixpy.coordinates.transforms import get_hpc_info from stixpy.net.client import STIXClient from stixpy.product import Product as STIXPYProduct -from sunpy.coordinates import HeliographicStonyhurst, Helioprojective +from sunpy.coordinates import HeliographicStonyhurst, Helioprojective, SphericalScreen from sunpy.map import make_fitswcs_header from sunpy.net import attrs as a from sunpy.time import TimeRange @@ -95,12 +96,26 @@ def add_flare_position( # helio_frame = Helioprojective(observer="earth") # SkyCoord(HeliographicStonyhurst(0 * u.deg, 0 * u.deg)) # SkyCoord(0 * u.deg, 0 * u.deg, frame=helio_frame) - data["flare_position"] = [SkyCoord(HeliographicStonyhurst(0 * u.deg, 0 * u.deg)) for i in range(0, len(data))] + + n = len(data) + + data["flareposition_obs_hgs_x"] = Column( + np.zeros(n, dtype=float) * u.km, description="HeliographicStonyhurst X of observer" + ) + data["flareposition_obs_hgs_y"] = Column( + np.zeros(n, dtype=float) * u.km, description="HeliographicStonyhurst Y of observer" + ) + data["flareposition_obs_hgs_z"] = Column( + np.zeros(n, dtype=float) * u.km, description="HeliographicStonyhurst Z of observer" + ) + data["flareposition_hp_tx"] = Column(np.zeros(n, dtype=float) * u.arcsec, description="Helioprojective Tx") + data["flareposition_hp_ty"] = Column(np.zeros(n, dtype=float) * u.arcsec, description="Helioprojective Ty") data["anc_ephemeris_path"] = Column(" " * 500, dtype=str, description="TDB") data["cpd_path"] = Column(" " * 500, dtype=str, description="TDB") data["_position_status"] = Column(False, dtype=bool, description="TDB") data["_position_message"] = Column(" " * 500, dtype=str, description="TDB") + to_remove = [] pass_filter = 0 no_ephemeris = 0 @@ -110,9 +125,9 @@ def add_flare_position( total_flares = len(data) day_asp_ephemeris_cache = dict() - flare_positions = [] + for i, row in enumerate(data): - if filter_function(row) and i < 200: + if filter_function(row): # and i < 200: pass_filter += 1 peak_time = row[peak_time_colname] start_time = row[start_time_colname] @@ -183,7 +198,7 @@ def add_flare_position( cpd_res["duration"][i] = header["OBT_END"] - header["OBT_BEG"] # TODO: add more criteria to select the best CPD file - cpd_res.sort(["tbins", "duration"]) + cpd_res.sort(["inc_peak", "tbins", "duration"], reverse=True) # cpd_res.pprint() best_cpd_idx = 0 else: @@ -193,24 +208,63 @@ def add_flare_position( try: stixpy_cpd = STIXPYProduct(Path(data[i]["cpd_path"])) - coord, map = estimate_stix_flare_location(stixpy_cpd) + time_range = TimeRange(max(peak_time - 20 * u.s, start_time), min(peak_time + 20 * u.s, end_time)) + overlaps = calculate_overlap(stixpy_cpd.time_range, time_range) + if overlaps is None: + logger.warning( + f"CPD data does not cover time range around peak time {time_range.start} to {time_range.end}" + ) + time_range = stixpy_cpd.time_range + contains_peak_time = False + else: + contains_peak_time = True + time_range = overlaps + + mask = (stixpy_cpd.data["time"] >= time_range.start) & (stixpy_cpd.data["time"] <= time_range.end) + data_at_peak = stixpy_cpd.data[mask] + if len(np.unique(data_at_peak["rcr"])) > 1: + logger.warning( + f"Multiple rcr values found for flare at time {time_range.start} : {time_range.end}" + ) + # allow a larger time range for finding a constant rcr sequence + if contains_peak_time: + time_range = TimeRange( + max(peak_time - 40 * u.s, start_time), min(peak_time + 40 * u.s, end_time) + ) + mask = (stixpy_cpd.data["time"] >= time_range.start) & ( + stixpy_cpd.data["time"] <= time_range.end + ) + data_at_peak = stixpy_cpd.data[mask] + length, start_idx, rcr = longest_constant_sequence(data_at_peak["rcr"].value) + time_range = TimeRange( + data_at_peak["time"][start_idx], data_at_peak["time"][start_idx + length - 1] + ) + logger.info( + f"Using time range {time_range.start} to {time_range.end} for flare at around {peak_time} with constant rcr={rcr}" + ) + + coord, map = estimate_stix_flare_location(stixpy_cpd, time_range=time_range) roll, solo_xyz, pointing = get_hpc_info(start_time, end_time) solo = HeliographicStonyhurst(*solo_xyz, obstime=peak_time, representation_type="cartesian") + with SphericalScreen(solo, only_off_disk=True): + center_hpc = coord.transform_to(Helioprojective(observer=solo)) + + data[i]["flareposition_obs_hgs_x"] = solo_xyz[0].to(u.km) + data[i]["flareposition_obs_hgs_y"] = solo_xyz[1].to(u.km) + data[i]["flareposition_obs_hgs_z"] = solo_xyz[2].to(u.km) + data[i]["flareposition_hp_tx"] = center_hpc.Tx.to(u.arcsec) + data[i]["flareposition_hp_ty"] = center_hpc.Ty.to(u.arcsec) - # data[i]["flare_position"] = coord.transform_to(Helioprojective(observer=solo)) - flare_positions.append(coord.transform_to(Helioprojective(observer=solo))) data[i]["_position_status"] = True data[i]["_position_message"] = "OK" except Exception as e: - flare_positions.append(None) + data[i]["_position_status"] = False data[i]["_position_message"] = f"Error: {type(e)}" - + logger.warn(f"Error calculating flare position for flare at time {start_time} : {end_time}: {e}") else: to_remove.append(i) - flare_positions.append(None) - data["flare_position"] = flare_positions if not keep_all_flares: data.remove_rows(to_remove) @@ -762,3 +816,44 @@ def add_peak_preview(cls, data, energies, parent, fido_client: STIXClient, img_p @classmethod def is_datasource_for(cls, *, service_type, service_subtype, ssid, **kwargs): return kwargs["level"] == "L3" and service_type == 0 and service_subtype == 0 and ssid == 8 + + +def longest_constant_sequence(state_array): + """Find the longest sequence where state is constant. + In case of equal length, prefer the one with the lower state value.""" + if len(state_array) == 0: + return 0, None, None + + max_length = 0 + max_state = None + max_start_idx = None + current_idx = 0 + + for state, group in groupby(state_array): + length = len(list(group)) + # Update if longer, OR if equal length but lower state value + if length > max_length or (length == max_length and (max_state is None or state < max_state)): + max_length = length + max_state = state + max_start_idx = current_idx + current_idx += length + + return max_length, max_start_idx, max_state + + +def calculate_overlap(range1, range2): + """Calculate the overlap between two TimeRanges. + Returns the overlap duration and the overlapping TimeRange, or None if no overlap.""" + + # Check if they intersect first + if not range1.intersects(range2): + return None + + # Calculate intersection boundaries + overlap_start = max(range1.start, range2.start) + overlap_end = min(range1.end, range2.end) + + # Create the overlapping TimeRange + overlap_range = TimeRange(overlap_start, overlap_end) + + return overlap_range diff --git a/stixcore/soop/manager.py b/stixcore/soop/manager.py index b3dd242e..142bfa42 100644 --- a/stixcore/soop/manager.py +++ b/stixcore/soop/manager.py @@ -529,6 +529,8 @@ def add_soop_file_to_index(self, path, *, rebuild_index=True, **args): all_soop_file = Path(CONFIG.get("SOOP", "soop_files_download")) / f"{plan}.{version}.all.json" if not all_soop_file.exists(): + # TODO remove ove SOOP API is working reliably + return self.download_all_soops_from_api(plan, version, all_soop_file) with open(all_soop_file) as f_all: From 1e3f10ca0f8abd7625fb011d63e3f3c17bc7e054 Mon Sep 17 00:00:00 2001 From: Nicky Hochmuth Date: Mon, 23 Mar 2026 10:29:33 +0100 Subject: [PATCH 3/6] SkyCoord for flarePosition Fits de/serialize hook for converting into fits serializeable icrs frame --- .../io/product_processors/fits/processors.py | 3 +- stixcore/products/level3/flarelist.py | 104 +++++++++++++----- stixcore/products/product.py | 15 +++ stixcore/products/tests/test_flarelist.py | 86 +++++++++++++++ 4 files changed, 180 insertions(+), 28 deletions(-) create mode 100644 stixcore/products/tests/test_flarelist.py diff --git a/stixcore/io/product_processors/fits/processors.py b/stixcore/io/product_processors/fits/processors.py index 7ca925b2..a9887344 100644 --- a/stixcore/io/product_processors/fits/processors.py +++ b/stixcore/io/product_processors/fits/processors.py @@ -1144,7 +1144,8 @@ def write_fits(self, prod, *, version=0): elif fitspath_complete.exists(): logger.warning("Complete Fits file %s exists will be overridden", fitspath.name) - data = prod.data + data = prod.data.copy() + prod.on_serialize(data) primary_header, header_override = self.generate_primary_header(filename, prod, version=version) diff --git a/stixcore/products/level3/flarelist.py b/stixcore/products/level3/flarelist.py index 9379f02f..ee4e750b 100644 --- a/stixcore/products/level3/flarelist.py +++ b/stixcore/products/level3/flarelist.py @@ -97,24 +97,12 @@ def add_flare_position( # SkyCoord(HeliographicStonyhurst(0 * u.deg, 0 * u.deg)) # SkyCoord(0 * u.deg, 0 * u.deg, frame=helio_frame) - n = len(data) - - data["flareposition_obs_hgs_x"] = Column( - np.zeros(n, dtype=float) * u.km, description="HeliographicStonyhurst X of observer" - ) - data["flareposition_obs_hgs_y"] = Column( - np.zeros(n, dtype=float) * u.km, description="HeliographicStonyhurst Y of observer" - ) - data["flareposition_obs_hgs_z"] = Column( - np.zeros(n, dtype=float) * u.km, description="HeliographicStonyhurst Z of observer" - ) - data["flareposition_hp_tx"] = Column(np.zeros(n, dtype=float) * u.arcsec, description="Helioprojective Tx") - data["flareposition_hp_ty"] = Column(np.zeros(n, dtype=float) * u.arcsec, description="Helioprojective Ty") - data["anc_ephemeris_path"] = Column(" " * 500, dtype=str, description="TDB") data["cpd_path"] = Column(" " * 500, dtype=str, description="TDB") data["_position_status"] = Column(False, dtype=bool, description="TDB") data["_position_message"] = Column(" " * 500, dtype=str, description="TDB") + tx_list, ty_list = [], [] + solo_x_list, solo_y_list, solo_z_list, peak_time_list = [], [], [], [] to_remove = [] pass_filter = 0 @@ -127,12 +115,11 @@ def add_flare_position( day_asp_ephemeris_cache = dict() for i, row in enumerate(data): - if filter_function(row): # and i < 200: + peak_time = row[peak_time_colname] + start_time = row[start_time_colname] + end_time = row[end_time_colname] + if filter_function(row): # and i < 60: pass_filter += 1 - peak_time = row[peak_time_colname] - start_time = row[start_time_colname] - end_time = row[end_time_colname] - day = peak_time.to_datetime().date() if day in day_asp_ephemeris_cache: @@ -249,21 +236,54 @@ def add_flare_position( solo = HeliographicStonyhurst(*solo_xyz, obstime=peak_time, representation_type="cartesian") with SphericalScreen(solo, only_off_disk=True): center_hpc = coord.transform_to(Helioprojective(observer=solo)) - - data[i]["flareposition_obs_hgs_x"] = solo_xyz[0].to(u.km) - data[i]["flareposition_obs_hgs_y"] = solo_xyz[1].to(u.km) - data[i]["flareposition_obs_hgs_z"] = solo_xyz[2].to(u.km) - data[i]["flareposition_hp_tx"] = center_hpc.Tx.to(u.arcsec) - data[i]["flareposition_hp_ty"] = center_hpc.Ty.to(u.arcsec) + tx_list.append(center_hpc.Tx) + ty_list.append(center_hpc.Ty) + solo_x_list.append(solo.cartesian.x) + solo_y_list.append(solo.cartesian.y) + solo_z_list.append(solo.cartesian.z) + peak_time_list.append(peak_time) data[i]["_position_status"] = True data[i]["_position_message"] = "OK" except Exception as e: data[i]["_position_status"] = False data[i]["_position_message"] = f"Error: {type(e)}" - logger.warn(f"Error calculating flare position for flare at time {start_time} : {end_time}: {e}") + logger.warning(f"Error calculating flare position for flare at time {start_time} : {end_time}: {e}") + tx_list.append(np.nan * u.arcsec) + ty_list.append(np.nan * u.arcsec) + solo_x_list.append(np.nan * u.km) + solo_y_list.append(np.nan * u.km) + solo_z_list.append(np.nan * u.km) + peak_time_list.append(peak_time) else: to_remove.append(i) + tx_list.append(np.nan * u.arcsec) + ty_list.append(np.nan * u.arcsec) + solo_x_list.append(np.nan * u.km) + solo_y_list.append(np.nan * u.km) + solo_z_list.append(np.nan * u.km) + peak_time_list.append(peak_time) + data[i]["_position_status"] = False + data[i]["_position_message"] = "flare did not pass the filter function" + + solo_times = Time(peak_time_list) + hgs_coords = SkyCoord( + u.Quantity(solo_x_list), + u.Quantity(solo_y_list), + u.Quantity(solo_z_list), + frame=HeliographicStonyhurst(obstime=solo_times), + representation_type="cartesian", + ) + + hp_coords = SkyCoord( + u.Quantity(tx_list), u.Quantity(ty_list), frame=Helioprojective(obstime=solo_times, observer=hgs_coords) + ) + + data["location_hgs"] = hgs_coords + # description="Flare location in Heliographic Stonyhurst coordinates" + + data["location_hp"] = hp_coords + # description="Flare location in Helioprojective coordinates" if not keep_all_flares: data.remove_rows(to_remove) @@ -276,6 +296,34 @@ def add_flare_position( f"finally {len(data)} flares remaining" ) + def on_serialize(self, data): + for col_name in ("location_hgs", "location_hp"): + if col_name in data.colnames: + icrs = data[col_name].icrs + icrs_coord = SkyCoord(icrs.ra, icrs.dec, icrs.distance, frame="icrs") + col_idx = data.colnames.index(col_name) + data.remove_column(col_name) + data.add_column(icrs_coord, name=col_name, index=col_idx) + s = super() + if hasattr(s, "on_serialize"): + s.on_serialize(data) + + def on_deserialize(self, data, *, peak_time_colname=None): + peak_col = peak_time_colname or self.peak_time_colname + if peak_col not in data.colnames: + logger.warning(f"on_deserialize: column '{peak_col}' not found, skipping location transform") + else: + obstime = Time(data[peak_col]) + if "location_hgs" in data.colnames: + data["location_hgs"] = data["location_hgs"].transform_to(HeliographicStonyhurst(obstime=obstime)) + if "location_hp" in data.colnames: + data["location_hp"] = data["location_hp"].transform_to( + Helioprojective(obstime=obstime, observer=data["location_hgs"]) + ) + s = super() + if hasattr(s, "on_deserialize"): + s.on_deserialize(data) + class FlareSOOPMixin: """_summary_""" @@ -614,6 +662,7 @@ def __init__(self, *, service_type=0, service_subtype=0, ssid=3, data, month, ** self.name = FlarelistSDCLoc.NAME self.ssid = 3 + self.peak_time_colname = "peak_UTC" def enhance_from_product(self, in_prod: GenericProduct): pass @@ -631,7 +680,7 @@ def add_flare_position(cls, data, fido_client: STIXClient, *, month=None): peak_time_colname="peak_UTC", start_time_colname="start_UTC", end_time_colname="end_UTC", - keep_all_flares=False, + keep_all_flares=True, month=month, ) @@ -751,6 +800,7 @@ def __init__(self, *, service_type=0, service_subtype=0, ssid=7, data, month, ** self.name = FlarelistSCLoc.NAME self.ssid = 7 + self.peak_time_colname = "peak_UTC" def enhance_from_product(self, in_prod: GenericProduct): pass diff --git a/stixcore/products/product.py b/stixcore/products/product.py index 27415375..7df9c3e8 100644 --- a/stixcore/products/product.py +++ b/stixcore/products/product.py @@ -336,6 +336,9 @@ def __call__(self, *args, **kwargs): month=month, ) + if hasattr(p, "on_deserialize") and callable(getattr(p, "on_deserialize")): + p.on_deserialize(p.data) + if hasattr(p, "get_additional_extensions") and data is not None: for _, name in p.get_additional_extensions(): # read the additional extension data @@ -609,6 +612,18 @@ def max_exposure(self): # default for FITS HEADER return 0.0 + def on_serialize(self, data): + """Hook called before writing data to FITS. Mixins override and chain via super().""" + s = super() + if hasattr(s, "on_serialize"): + s.on_serialize(data) + + def on_deserialize(self, data): + """Hook called after reading data from FITS. Mixins override and chain via super().""" + s = super() + if hasattr(s, "on_deserialize"): + s.on_deserialize(data) + def find_parent_products(self, root): """ Convenient way to get access to the parent products. diff --git a/stixcore/products/tests/test_flarelist.py b/stixcore/products/tests/test_flarelist.py new file mode 100644 index 00000000..ff847b6c --- /dev/null +++ b/stixcore/products/tests/test_flarelist.py @@ -0,0 +1,86 @@ +from datetime import date + +import numpy as np +import pytest +from sunpy.coordinates import HeliographicStonyhurst, Helioprojective + +import astropy.units as u +from astropy.coordinates import SkyCoord +from astropy.io import fits +from astropy.table import QTable +from astropy.time import Time + +from stixcore.io.product_processors.fits.processors import FitsL3Processor +from stixcore.products.level3.flarelist import FlarelistSDCLoc +from stixcore.products.product import Product + +N = 10 + + +@pytest.fixture +def flare_data(): + peak_times = Time("2022-01-01T12:00:00") + np.arange(N) * 600 * u.s + + hgs_coords = SkyCoord( + lon=np.linspace(0, 30, N) * u.deg, + lat=np.linspace(-5, 5, N) * u.deg, + radius=np.ones(N) * 1.0 * u.AU, + frame=HeliographicStonyhurst(obstime=peak_times), + ) + + hp_coords = SkyCoord( + Tx=np.linspace(-300, 300, N) * u.arcsec, + Ty=np.linspace(-200, 200, N) * u.arcsec, + frame=Helioprojective(obstime=peak_times, observer=hgs_coords), + ) + + data = QTable() + data["peak_UTC"] = peak_times + data["start_UTC"] = peak_times - 60 * u.s + data["end_UTC"] = peak_times + 60 * u.s + data["duration"] = np.ones(N) * 120 * u.s + data["lc_peak"] = np.ones((N, 5)) * u.ct / u.s + data["location_hgs"] = hgs_coords + data["location_hp"] = hp_coords + + return data + + +def test_flarelist_sdcloc_location_roundtrip(flare_data, tmp_path): + prod = FlarelistSDCLoc( + data=flare_data, + month=date(2022, 1, 1), + control=QTable(), + ) + + # minimal header bypasses the Spice-dependent header generation chain + header = fits.Header() + header["LEVEL"] = "L3" + header["STYPE"] = 0 + header["SSTYPE"] = 0 + header["SSID"] = 3 + header["DATE-BEG"] = "2022-01-01T00:00:00" + prod.fits_header = header + + # energy/additional_header_keywords are not set for freshly created products + prod.energy = None + prod._additional_header_keywords = [] + + orig_hgs_lon = prod.data["location_hgs"].lon.copy() + orig_hgs_lat = prod.data["location_hgs"].lat.copy() + orig_hp_tx = prod.data["location_hp"].Tx.copy() + orig_hp_ty = prod.data["location_hp"].Ty.copy() + + # write via FitsL3Processor — calls on_serialize internally, prod.data unchanged + writer = FitsL3Processor(tmp_path) + written_file_name = writer.write_fits(prod) + assert len(written_file_name) == 1 + + # read back via Product factory — calls on_deserialize internally + recovered = Product(written_file_name[0]) + + assert isinstance(recovered, FlarelistSDCLoc) + assert u.allclose(recovered.data["location_hgs"].lon, orig_hgs_lon, atol=1e-6 * u.deg) + assert u.allclose(recovered.data["location_hgs"].lat, orig_hgs_lat, atol=1e-6 * u.deg) + assert u.allclose(recovered.data["location_hp"].Tx, orig_hp_tx, atol=1e-3 * u.arcsec) + assert u.allclose(recovered.data["location_hp"].Ty, orig_hp_ty, atol=1e-3 * u.arcsec) From a888dcb38a47be708dfa3491e2a9f47fe1249171 Mon Sep 17 00:00:00 2001 From: Nicky Hochmuth Date: Thu, 26 Mar 2026 09:45:09 +0100 Subject: [PATCH 4/6] change on_deserialize flow --- .../io/product_processors/fits/processors.py | 5 +- stixcore/products/level3/flarelist.py | 124 ++++++++++-------- stixcore/products/product.py | 21 +-- stixcore/products/tests/test_flarelist.py | 67 ++++++---- 4 files changed, 128 insertions(+), 89 deletions(-) diff --git a/stixcore/io/product_processors/fits/processors.py b/stixcore/io/product_processors/fits/processors.py index a9887344..fda79ca1 100644 --- a/stixcore/io/product_processors/fits/processors.py +++ b/stixcore/io/product_processors/fits/processors.py @@ -823,10 +823,11 @@ def generate_primary_header(self, filename, product, *, version=0): if default[0] not in soop_key_names: soop_headers += tuple([default]) + scet_range = product.scet_timerange time_headers = ( # Name, Value, Comment - ("OBT_BEG", product.scet_timerange.start.as_float().value, "Start of acquisition time in OBT"), - ("OBT_END", product.scet_timerange.end.as_float().value, "End of acquisition time in OBT"), + ("OBT_BEG", scet_range.start.as_float().value, "Start of acquisition time in OBT"), + ("OBT_END", scet_range.end.as_float().value, "End of acquisition time in OBT"), ("TIMESYS", "UTC", "System used for time keywords"), ("LEVEL", "L1", "Processing level of the data"), ("DATE-OBS", product.utc_timerange.start.fits, "Start of acquisition time in UTC"), diff --git a/stixcore/products/level3/flarelist.py b/stixcore/products/level3/flarelist.py index ee4e750b..5abcc611 100644 --- a/stixcore/products/level3/flarelist.py +++ b/stixcore/products/level3/flarelist.py @@ -11,7 +11,7 @@ from stixpy.coordinates.transforms import get_hpc_info from stixpy.net.client import STIXClient from stixpy.product import Product as STIXPYProduct -from sunpy.coordinates import HeliographicStonyhurst, Helioprojective, SphericalScreen +from sunpy.coordinates import HeliographicStonyhurst, Helioprojective from sunpy.map import make_fitswcs_header from sunpy.net import attrs as a from sunpy.time import TimeRange @@ -77,7 +77,21 @@ def make_stix_fitswcs_header(data, flare_position, *, scale, exposure, rotation_ return header -class FlarePositionMixin: +class _SerializeMixin: + """No-op chain terminator for on_serialize/on_deserialize. + + Functional mixins inherit from this so super() calls always land safely + instead of hitting object and raising AttributeError. + """ + + def on_serialize(self, data): + pass + + def on_deserialize(self, data, **kwargs): + pass + + +class FlarePositionMixin(_SerializeMixin): """_summary_""" @classmethod @@ -101,8 +115,8 @@ def add_flare_position( data["cpd_path"] = Column(" " * 500, dtype=str, description="TDB") data["_position_status"] = Column(False, dtype=bool, description="TDB") data["_position_message"] = Column(" " * 500, dtype=str, description="TDB") - tx_list, ty_list = [], [] - solo_x_list, solo_y_list, solo_z_list, peak_time_list = [], [], [], [] + # tx_list, ty_list = [], [] + solo_cartesian_list = [] to_remove = [] pass_filter = 0 @@ -118,6 +132,7 @@ def add_flare_position( peak_time = row[peak_time_colname] start_time = row[start_time_colname] end_time = row[end_time_colname] + logger.info(f"Processing flare {i}/{len(data)} at time {start_time} : {end_time} (peak at {peak_time})") if filter_function(row): # and i < 60: pass_filter += 1 day = peak_time.to_datetime().date() @@ -137,6 +152,8 @@ def add_flare_position( logger.warning(f"No ephemeris data found for flare at time {start_time} : {end_time}") data[i]["_position_message"] = "no ephemeris data found" no_ephemeris += 1 + solo_cartesian_list.append((np.nan * u.km, np.nan * u.km, np.nan * u.km)) + continue data[i]["anc_ephemeris_path"] = anc_res["path"][0] @@ -153,6 +170,8 @@ def add_flare_position( logger.warning(f"No CPD data found for flare at time {start_time} : {end_time}") data[i]["_position_message"] = "no CPD data found" no_cpd += 1 + solo_cartesian_list.append((np.nan * u.km, np.nan * u.km, np.nan * u.km)) + continue if len(cpd_res) > 1: logger.debug(f"Many CPD data found for flare at time {start_time} : {end_time}") @@ -234,14 +253,11 @@ def add_flare_position( roll, solo_xyz, pointing = get_hpc_info(start_time, end_time) solo = HeliographicStonyhurst(*solo_xyz, obstime=peak_time, representation_type="cartesian") - with SphericalScreen(solo, only_off_disk=True): - center_hpc = coord.transform_to(Helioprojective(observer=solo)) - tx_list.append(center_hpc.Tx) - ty_list.append(center_hpc.Ty) - solo_x_list.append(solo.cartesian.x) - solo_y_list.append(solo.cartesian.y) - solo_z_list.append(solo.cartesian.z) - peak_time_list.append(peak_time) + # with SphericalScreen(solo, only_off_disk=True): + # center_hpc = coord.transform_to(Helioprojective(observer=solo)) + # tx_list.append(center_hpc.Tx) + # ty_list.append(center_hpc.Ty) + solo_cartesian_list.append((solo.cartesian.x, solo.cartesian.y, solo.cartesian.z)) data[i]["_position_status"] = True data[i]["_position_message"] = "OK" @@ -249,40 +265,36 @@ def add_flare_position( data[i]["_position_status"] = False data[i]["_position_message"] = f"Error: {type(e)}" logger.warning(f"Error calculating flare position for flare at time {start_time} : {end_time}: {e}") - tx_list.append(np.nan * u.arcsec) - ty_list.append(np.nan * u.arcsec) - solo_x_list.append(np.nan * u.km) - solo_y_list.append(np.nan * u.km) - solo_z_list.append(np.nan * u.km) - peak_time_list.append(peak_time) + # tx_list.append(np.nan * u.arcsec) + # ty_list.append(np.nan * u.arcsec) + solo_cartesian_list.append((np.nan * u.km, np.nan * u.km, np.nan * u.km)) + else: to_remove.append(i) - tx_list.append(np.nan * u.arcsec) - ty_list.append(np.nan * u.arcsec) - solo_x_list.append(np.nan * u.km) - solo_y_list.append(np.nan * u.km) - solo_z_list.append(np.nan * u.km) - peak_time_list.append(peak_time) + # tx_list.append(np.nan * u.arcsec) + # ty_list.append(np.nan * u.arcsec) + solo_cartesian_list.append((np.nan * u.km, np.nan * u.km, np.nan * u.km)) data[i]["_position_status"] = False data[i]["_position_message"] = "flare did not pass the filter function" - solo_times = Time(peak_time_list) + solo_times = Time(data[peak_time_colname]) + solo_x, solo_y, solo_z = zip(*solo_cartesian_list) hgs_coords = SkyCoord( - u.Quantity(solo_x_list), - u.Quantity(solo_y_list), - u.Quantity(solo_z_list), + u.Quantity(solo_x), + u.Quantity(solo_y), + u.Quantity(solo_z), frame=HeliographicStonyhurst(obstime=solo_times), representation_type="cartesian", ) - hp_coords = SkyCoord( - u.Quantity(tx_list), u.Quantity(ty_list), frame=Helioprojective(obstime=solo_times, observer=hgs_coords) - ) + # hp_coords = SkyCoord( + # u.Quantity(tx_list), u.Quantity(ty_list), frame=Helioprojective(obstime=solo_times, observer=hgs_coords) + # ) data["location_hgs"] = hgs_coords # description="Flare location in Heliographic Stonyhurst coordinates" - data["location_hp"] = hp_coords + # data["location_hp"] = hp_coords # description="Flare location in Helioprojective coordinates" if not keep_all_flares: @@ -297,35 +309,31 @@ def add_flare_position( ) def on_serialize(self, data): - for col_name in ("location_hgs", "location_hp"): - if col_name in data.colnames: - icrs = data[col_name].icrs - icrs_coord = SkyCoord(icrs.ra, icrs.dec, icrs.distance, frame="icrs") - col_idx = data.colnames.index(col_name) - data.remove_column(col_name) - data.add_column(icrs_coord, name=col_name, index=col_idx) - s = super() - if hasattr(s, "on_serialize"): - s.on_serialize(data) - - def on_deserialize(self, data, *, peak_time_colname=None): + logger.info("FlarePositionMixin on_serialize called, transforming location columns to ICRS for serialization") + + if "location_hgs" in data.colnames: + icrs = data["location_hgs"].icrs + icrs_coord = SkyCoord(icrs.ra, icrs.dec, icrs.distance, frame="icrs") + col_idx = data.colnames.index("location_hgs") + data.remove_column("location_hgs") + data.add_column(icrs_coord, name="location_icrs", index=col_idx) + super().on_serialize(data) + + def on_deserialize(self, data, *, peak_time_colname=None, **kwargs): + logger.info( + "FlarePositionMixin on_deserialize called, transforming location columns back to heliographic coordinates" + ) peak_col = peak_time_colname or self.peak_time_colname if peak_col not in data.colnames: logger.warning(f"on_deserialize: column '{peak_col}' not found, skipping location transform") else: obstime = Time(data[peak_col]) - if "location_hgs" in data.colnames: - data["location_hgs"] = data["location_hgs"].transform_to(HeliographicStonyhurst(obstime=obstime)) - if "location_hp" in data.colnames: - data["location_hp"] = data["location_hp"].transform_to( - Helioprojective(obstime=obstime, observer=data["location_hgs"]) - ) - s = super() - if hasattr(s, "on_deserialize"): - s.on_deserialize(data) + if "location_icrs" in data.colnames: + data["location_hgs"] = data["location_icrs"].transform_to(HeliographicStonyhurst(obstime=obstime)) + super().on_deserialize(data, **kwargs) -class FlareSOOPMixin: +class FlareSOOPMixin(_SerializeMixin): """_summary_""" @classmethod @@ -352,6 +360,14 @@ def add_soop( data["soop_id"] = Column(soop_id, dtype=str, description="SOOP ID") data["soop_type"] = Column(soop_type, dtype=str, description="name of the SOOP campaign") + # def on_serialize(self, data): + # logger.info("FlareSOOPMixin on_serialize called, but no special handling implemented for SOOP data") + # super().on_serialize(data) + + # def on_deserialize(self, data, **kwargs): + # logger.info("FlareSOOPMixin on_deserialize called, but no special handling implemented for SOOP data") + # super().on_deserialize(data, **kwargs) + class FlarePeakPreviewMixin: """Mixin class to add peak preview images to flare list products. diff --git a/stixcore/products/product.py b/stixcore/products/product.py index 7df9c3e8..3dce6632 100644 --- a/stixcore/products/product.py +++ b/stixcore/products/product.py @@ -613,16 +613,21 @@ def max_exposure(self): return 0.0 def on_serialize(self, data): - """Hook called before writing data to FITS. Mixins override and chain via super().""" - s = super() - if hasattr(s, "on_serialize"): - s.on_serialize(data) + """Hook called before writing data to FITS. Mixins override and chain via super(). - def on_deserialize(self, data): + Uses getattr so plain products without functional mixins are safe — + the functional mixins (FlarePositionMixin etc.) appear after GenericProduct + in the MRO, so pass would stop the chain before reaching them. + """ + serialize = getattr(super(), "on_serialize", None) + if serialize is not None: + serialize(data) + + def on_deserialize(self, data, **kwargs): """Hook called after reading data from FITS. Mixins override and chain via super().""" - s = super() - if hasattr(s, "on_deserialize"): - s.on_deserialize(data) + deserialize = getattr(super(), "on_deserialize", None) + if deserialize is not None: + deserialize(data, **kwargs) def find_parent_products(self, root): """ diff --git a/stixcore/products/tests/test_flarelist.py b/stixcore/products/tests/test_flarelist.py index ff847b6c..842e3960 100644 --- a/stixcore/products/tests/test_flarelist.py +++ b/stixcore/products/tests/test_flarelist.py @@ -2,12 +2,13 @@ import numpy as np import pytest -from sunpy.coordinates import HeliographicStonyhurst, Helioprojective +from sunpy.coordinates import HeliographicStonyhurst import astropy.units as u from astropy.coordinates import SkyCoord from astropy.io import fits from astropy.table import QTable +from astropy.tests.helper import assert_quantity_allclose from astropy.time import Time from stixcore.io.product_processors.fits.processors import FitsL3Processor @@ -21,19 +22,19 @@ def flare_data(): peak_times = Time("2022-01-01T12:00:00") + np.arange(N) * 600 * u.s + lon = np.linspace(0, 30, N) + lat = np.linspace(-5, 5, N) + # mark same rows fully NaN so the entire SkyCoord row is invalid + lon[2] = lat[2] = np.nan + lon[7] = lat[7] = np.nan + hgs_coords = SkyCoord( - lon=np.linspace(0, 30, N) * u.deg, - lat=np.linspace(-5, 5, N) * u.deg, + lon=lon * u.deg, + lat=lat * u.deg, radius=np.ones(N) * 1.0 * u.AU, frame=HeliographicStonyhurst(obstime=peak_times), ) - hp_coords = SkyCoord( - Tx=np.linspace(-300, 300, N) * u.arcsec, - Ty=np.linspace(-200, 200, N) * u.arcsec, - frame=Helioprojective(obstime=peak_times, observer=hgs_coords), - ) - data = QTable() data["peak_UTC"] = peak_times data["start_UTC"] = peak_times - 60 * u.s @@ -41,12 +42,12 @@ def flare_data(): data["duration"] = np.ones(N) * 120 * u.s data["lc_peak"] = np.ones((N, 5)) * u.ct / u.s data["location_hgs"] = hgs_coords - data["location_hp"] = hp_coords return data -def test_flarelist_sdcloc_location_roundtrip(flare_data, tmp_path): +@pytest.fixture +def written_fits(flare_data, tmp_path): prod = FlarelistSDCLoc( data=flare_data, month=date(2022, 1, 1), @@ -61,26 +62,42 @@ def test_flarelist_sdcloc_location_roundtrip(flare_data, tmp_path): header["SSID"] = 3 header["DATE-BEG"] = "2022-01-01T00:00:00" prod.fits_header = header - - # energy/additional_header_keywords are not set for freshly created products prod.energy = None prod._additional_header_keywords = [] + writer = FitsL3Processor(tmp_path) + written = writer.write_fits(prod) + assert len(written) == 1 + + return prod, written[0] + + +def test_flarelist_sdcloc_location_roundtrip(written_fits): + prod, fits_path = written_fits orig_hgs_lon = prod.data["location_hgs"].lon.copy() orig_hgs_lat = prod.data["location_hgs"].lat.copy() - orig_hp_tx = prod.data["location_hp"].Tx.copy() - orig_hp_ty = prod.data["location_hp"].Ty.copy() - - # write via FitsL3Processor — calls on_serialize internally, prod.data unchanged - writer = FitsL3Processor(tmp_path) - written_file_name = writer.write_fits(prod) - assert len(written_file_name) == 1 # read back via Product factory — calls on_deserialize internally - recovered = Product(written_file_name[0]) + recovered = Product(fits_path) assert isinstance(recovered, FlarelistSDCLoc) - assert u.allclose(recovered.data["location_hgs"].lon, orig_hgs_lon, atol=1e-6 * u.deg) - assert u.allclose(recovered.data["location_hgs"].lat, orig_hgs_lat, atol=1e-6 * u.deg) - assert u.allclose(recovered.data["location_hp"].Tx, orig_hp_tx, atol=1e-3 * u.arcsec) - assert u.allclose(recovered.data["location_hp"].Ty, orig_hp_ty, atol=1e-3 * u.arcsec) + assert_quantity_allclose(recovered.data["location_hgs"].lon, orig_hgs_lon, atol=1e-6 * u.deg, equal_nan=True) + assert_quantity_allclose(recovered.data["location_hgs"].lat, orig_hgs_lat, atol=1e-6 * u.deg, equal_nan=True) + + +def test_flarelist_sdcloc_fits_stores_icrs(written_fits): + prod, fits_path = written_fits + orig_hgs_lon = prod.data["location_hgs"].lon.copy() + orig_hgs_lat = prod.data["location_hgs"].lat.copy() + + # read the DATA extension directly — no on_deserialize, raw FITS content + raw = QTable.read(fits_path, hdu="DATA", astropy_native=True) + + assert "location_hgs" not in raw.colnames, "HGS column should not be stored in FITS" + assert "location_icrs" in raw.colnames, "ICRS column should be present in FITS" + + # manually transform ICRS back to HGS and compare with original + obstime = Time(raw["peak_UTC"]) + hgs = raw["location_icrs"].transform_to(HeliographicStonyhurst(obstime=obstime)) + assert_quantity_allclose(hgs.lon, orig_hgs_lon, atol=1e-6 * u.deg, equal_nan=True) + assert_quantity_allclose(hgs.lat, orig_hgs_lat, atol=1e-6 * u.deg, equal_nan=True) From 4bbc2e85d8fe45d6945f797451114905d32bebbd Mon Sep 17 00:00:00 2001 From: Nicky Hochmuth Date: Thu, 26 Mar 2026 09:55:48 +0100 Subject: [PATCH 5/6] enable soop update --- stixcore/soop/manager.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/stixcore/soop/manager.py b/stixcore/soop/manager.py index 142bfa42..b3dd242e 100644 --- a/stixcore/soop/manager.py +++ b/stixcore/soop/manager.py @@ -529,8 +529,6 @@ def add_soop_file_to_index(self, path, *, rebuild_index=True, **args): all_soop_file = Path(CONFIG.get("SOOP", "soop_files_download")) / f"{plan}.{version}.all.json" if not all_soop_file.exists(): - # TODO remove ove SOOP API is working reliably - return self.download_all_soops_from_api(plan, version, all_soop_file) with open(all_soop_file) as f_all: From 776b72132a3cf393037dd461ff3bee92c258cae8 Mon Sep 17 00:00:00 2001 From: Nicky Hochmuth Date: Fri, 27 Mar 2026 14:36:33 +0100 Subject: [PATCH 6/6] experimental fix for reading aspect burst data --- stixcore/products/product.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/stixcore/products/product.py b/stixcore/products/product.py index 3dce6632..01c29478 100644 --- a/stixcore/products/product.py +++ b/stixcore/products/product.py @@ -74,7 +74,10 @@ def read_qtable(file, hdu, hdul=None): `astropy.table.QTable` The corrected QTable with correct data types """ - qtable = QTable.read(file, hdu, astropy_native=True) + astropy_native = True + if (hdu.upper() == "DATA") and (file.name.startswith("solo_L0_stix-sci-aspect-burst")): + astropy_native = False + qtable = QTable.read(file, hdu, astropy_native=astropy_native) if hdul is None: hdul = fits.open(file) @@ -91,11 +94,12 @@ def read_qtable(file, hdu, hdul=None): if hasattr(dtype, "subdtype"): dtype = dtype.base - + # qtable[col.name] = qtable[col.name].astype(dtype) if col.coord_type != "UTC": qtable[col.name] = qtable[col.name].astype(dtype) else: - qtable[col.name].format = "isot" + # qtable[col.name].format = "isot" + pass return qtable