diff --git a/.gitignore b/.gitignore index c942678b..d4d6d599 100644 --- a/.gitignore +++ b/.gitignore @@ -166,3 +166,4 @@ stixcore/data/test/idb/v2.26.38/idb.sqlite .python-version sunpy/_version.py stixcore/_version.py +.claude/* diff --git a/stixcore/data/stixcore.ini b/stixcore/data/stixcore.ini index cdc8b347..3d7afdd7 100644 --- a/stixcore/data/stixcore.ini +++ b/stixcore/data/stixcore.ini @@ -36,3 +36,15 @@ soop_files_download = ./stixcore/data/soop ecc_path = /opt/stix_det_cal/bin/ [Processing] flarelist_sdc_min_count = 1000 +# background-data-file lookup (find_background_file_for_time) +flarelist_bkg_window_past_days = 21 +flarelist_bkg_window_future_days = 21 +flarelist_bkg_min_duration_s = 500 +# require the background file to share the flare's on-board ELUT configuration (ELUTManager lookup) +flarelist_bkg_require_same_elut = True +# optional stricter selection (default OFF - all three move common cases, not just edge cases): +# purpose preference: prefer purpose=="Background" unless a subject-only match is this many days closer +# (0.0 disables); exclude "elevated" backgrounds; drop requests whose comment references a flare id +flarelist_bkg_purpose_penalty_days = 0.0 +flarelist_bkg_exclude_elevated = False +flarelist_bkg_exclude_flare_comment = True diff --git a/stixcore/data/test/publish/rid_lut.csv b/stixcore/data/test/publish/rid_lut.csv index a6f03173..1eae888e 100644 --- a/stixcore/data/test/publish/rid_lut.csv +++ b/stixcore/data/test/publish/rid_lut.csv @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:2ba33aa45a7d045dcbc08d1ddddf8bdf8d01dde93e01fdf4db0c12e1dfb0b4bf -size 1005 +oid sha256:577485ea766335f15fb1ac4e40f1803d96e82919fd77164c429029a11f7afcac +size 1383 diff --git a/stixcore/ephemeris/manager.py b/stixcore/ephemeris/manager.py index 787e9490..adee19ff 100644 --- a/stixcore/ephemeris/manager.py +++ b/stixcore/ephemeris/manager.py @@ -388,6 +388,20 @@ def get_sun_disc_size(self, *, date): return rsun_arc + def get_earth_solo_time_shift(self, *, date): + """gets Time(Sun to Earth) - Time(Sun to S/C) + + Returns + ------- + `astropy.units.Quantity` + Time difference between Sun to Earth and Sun to S/C in seconds + """ + et = spiceypy.scs2e(SOLAR_ORBITER_ID, str(date)) + solo_sun_hg, sun_solo_lt = spiceypy.spkezr("SOLO", et, "SUN_EARTH_CEQU", "None", "Sun") + sun_earth_hee, sun_earth_lt = spiceypy.spkezr("Earth", et, "SOLO_HEE", "None", "Sun") + + return (sun_earth_lt - sun_solo_lt) * u.s + def get_position(self, *, date, frame): """ Get the position of SolarOrbiter at the given date in the given coordinate frame. diff --git a/stixcore/io/FlareListManager.py b/stixcore/io/FlareListManager.py index 38359be6..86f492e4 100644 --- a/stixcore/io/FlareListManager.py +++ b/stixcore/io/FlareListManager.py @@ -1,28 +1,357 @@ import sys import time +from pathlib import Path from datetime import datetime, timedelta +from collections import namedtuple import numpy as np import pandas as pd from stixdcpy.net import Request as stixdcpy_req +from stixpy.calibration.livetime import get_livetime_fraction +from stixpy.product import Product as STIXPYProduct from sunpy.net import attrs as a import astropy.units as u from astropy.table import Column, QTable, vstack from astropy.time import Time +from stixcore.calibration.elut_manager import ELUTManager from stixcore.config.config import CONFIG +from stixcore.io.product_processors.fits.processors import CreateUtcColumn +from stixcore.io.RidLutManager import ( + DEFAULT_BKG_EXCLUDE_KEYWORDS, + DEFAULT_BKG_KEYWORDS, + RidLutManager, + search_background_candidates, +) from stixcore.products.level3.flarelist import FlarelistSC, FlarelistSDC from stixcore.products.product import Product from stixcore.util.logging import get_logger from stixcore.util.singleton import Singleton from stixcore.util.util import url_to_path -__all__ = ["FlareListManager", "SDCFlareListManager", "SCFlareListManager"] +__all__ = [ + "FlareListManager", + "SDCFlareListManager", + "SCFlareListManager", + "compute_ql_count_rate", + "build_month_timeline", + "nearest_bin_index", + "max_rcr_in_window", + "find_background_file_for_time", + "BackgroundSelection", +] logger = get_logger(__name__) +def compute_ql_count_rate(counts, timedel, triggers, energy_delta, *, n_detectors): + """Reproduce stixpy's QL count-rate normalization -> ``ct / (s * keV)``. + + Mirrors ``stixpy.timeseries.quicklook`` (lightcurve uses ``n_detectors=16``, + background uses ``n_detectors=1``). Pure, no I/O. + + Parameters + ---------- + counts : `~astropy.units.Quantity` + Raw counts, shape ``(N, 5)`` in ``ct``. + timedel : `~astropy.units.Quantity` + Bin durations, shape ``(N,)``. + triggers : array-like + Trigger counts, shape ``(N,)`` or ``(N, 1)``. + energy_delta : `~astropy.units.Quantity` + Channel widths, shape ``(5,)`` in ``keV``. + n_detectors : int + 16 for the lightcurve, 1 for the background detector. + + Returns + ------- + `~astropy.units.Quantity` + Count rate, shape ``(N, 5)`` in ``ct / (s * keV)``. + """ + timedel = timedel.to(u.s) + trig = np.asarray(triggers).reshape(-1) + live_frac, *_ = get_livetime_fraction(trig / (n_detectors * timedel)) + return counts / ((timedel * live_frac).reshape(-1, 1) * energy_delta) + + +def build_month_timeline(daily_data_tables): + """Stack per-day QTables into one time-sorted timeline with unique timestamps. + + ``None``/empty inputs are ignored; an empty list yields an empty ``QTable``. + After sorting by ``time`` duplicate timestamps (day-boundary overlaps) are + dropped so the timeline is strictly increasing. + """ + tables = [t for t in daily_data_tables if t is not None and len(t) > 0] + if not tables: + return QTable() + timeline = vstack(tables, metadata_conflicts="silent") + timeline.sort("time") + if len(timeline) > 1: + keep = np.ones(len(timeline), dtype=bool) + keep[1:] = np.diff(timeline["time"].jd) > 0 + timeline = timeline[keep] + return timeline + + +def nearest_bin_index(times, target, tol): + """Index of the bin in ``times`` nearest ``target``. + + Returns ``None`` if ``times`` is empty or the nearest gap exceeds ``tol`` + (both ``target`` and ``times`` are `~astropy.time.Time`, ``tol`` a duration). + """ + if times is None or len(times) == 0: + return None + dt = np.abs((times - target).to_value(u.s)) + j = int(np.argmin(dt)) + if dt[j] > tol.to_value(u.s): + return None + return j + + +def max_rcr_in_window(times, rcr, start, end, *, fallback): + """Highest ``rcr`` for bins with ``start <= time <= end``; ``fallback`` if none.""" + if times is None or len(times) == 0: + return fallback + mask = (times >= start) & (times <= end) + if not np.any(mask): + return fallback + return int(np.asarray(rcr)[mask].max()) + + +#: Result of :func:`find_background_file_for_time`. ``path`` is a `~pathlib.Path` +#: (or ``None`` when nothing qualifies), ``rid`` the selected BSD request id +#: (``-1`` when none), and ``valid_from``/``valid_to`` the `~astropy.time.Time` +#: interval over which this selection stays valid for a time-ordered caller. +BackgroundSelection = namedtuple("BackgroundSelection", ["path", "rid", "valid_from", "valid_to"]) + + +def _rid_from_filename(path): + """Parse the BSD request id embedded in a science FITS filename. + + Mirrors ``stixcore.processing.publish`` (the 6th ``_``-separated segment is + ``-``). Returns ``None`` when the name doesn't carry one. + """ + parts = Path(path).name.split("_") + if len(parts) <= 5: + return None + try: + return int(parts[5].replace(".fits", "").split("-")[0]) + except ValueError: + return None + + +def _elut_id(time): + """Identifier of the ELUT active at ``time`` (the resolved ELUT filename), or + ``None`` when the ELUT index has no unambiguous entry for it. Metadata-only + lookup via `~stixcore.calibration.elut_manager.ELUTManager` — reads no FITS.""" + try: + return ELUTManager.instance._find_elut_file(time.to_datetime()) + except Exception as e: + logger.debug(f"no ELUT resolved for {getattr(time, 'isot', time)}: {e}") + return None + + +def _effective_crossover(t0, sa, pa, sb, pb): + """Earliest ``t >= t0`` (all in JD days) at which candidate ``b``'s effective + distance ``|t - sb| + pb`` drops to/below ``a``'s ``|t - sa| + pa``; ``None`` if + it never does. Bounds the nearest-in-time validity interval under the purpose + penalty (piecewise-linear, breakpoints at the two starts).""" + + def val(x): + return abs(x - sb) - abs(x - sa) + (pb - pa) + + if val(t0) < 0: + return t0 + lefts = [t0] + sorted(x for x in (sa, sb) if x > t0) + for i, left in enumerate(lefts): + right = lefts[i + 1] if i + 1 < len(lefts) else None + v = val(left) + if v < 0: + return left + probe = (left + right) / 2 if right is not None else left + 1.0 + slope = ((abs(probe - sb) - abs(probe - sa)) - (abs(left - sb) - abs(left - sa))) / (probe - left) + if slope < 0: + tc = left + v / (-slope) + if right is None or tc <= right: + return tc + return None + + +def find_background_file_for_time( + time, + *, + fido_client, + rid_lut=None, + window_past=None, + window_future=None, + min_duration=None, + require_same_elut=None, + purpose_penalty=None, + keywords=DEFAULT_BKG_KEYWORDS, + exclude_keywords=None, + exclude_flare_comment=None, +): + """Find the best quiet-time background CPD file applicable at ``time``. + + Strategy (see :func:`stixcore.io.RidLutManager.search_background_candidates`): + rank background requests from the RID LUT **nearest-in-time first** (closest + request start, past or future, within the separate ``window_past`` / + ``window_future`` bounds). For each candidate resolve the real ``sci_xray_cpd`` + L1 file(s) via FIDO, keep those whose filename carries the candidate rid, and + accept the first one that passes the "good background" checks (requested long + enough, attenuator out over the whole interval). + + Three optional stricter filters exist but default **off** (they move common + cases, not just edge cases): a ``purpose == "Background"`` preference + (``purpose_penalty`` > 0), dropping "elevated" backgrounds (``exclude_keywords``), + and dropping requests whose comment references a specific flare id + (``exclude_flare_comment``). See the ``[Processing]`` config keys. + + Alongside the winning file a validity interval ``[valid_from, valid_to]`` is + returned so a time-ordered caller (monthly flare processing) can reuse the + result for every later time inside the interval. ``valid_to`` is the earliest + time at which another candidate's effective distance would overtake the chosen + one, or the chosen start plus ``window_past`` when there is none — capped at + ``time + window_future`` so newly-reachable candidates are re-scanned. A + "no file found" result is cached until the next candidate could appear. + + Parameters + ---------- + time : `~astropy.time.Time` or str + The query time (e.g. a flare peak). + fido_client : `~stixpy.net.client.STIXClient` + Client used for the ``sci_xray_cpd`` L1 search. + rid_lut : `~astropy.table.Table`, optional + The RID LUT; defaults to ``RidLutManager.instance.rid_lut``. + window_past, window_future : `~astropy.units.Quantity`, optional + Search windows; default to the ``[Processing]`` config keys + ``flarelist_bkg_window_past_days`` / ``flarelist_bkg_window_future_days``. + min_duration : `~astropy.units.Quantity`, optional + Minimum requested integration time; defaults to + ``flarelist_bkg_min_duration_s``. + require_same_elut : bool, optional + If True, a candidate is only accepted when the ELUT active at its request + start (per `~stixcore.calibration.elut_manager.ELUTManager`) matches the + one active at ``time`` — i.e. the background was taken under the same + on-board ELUT configuration as the flare. Defaults to the ``[Processing]`` + config key ``flarelist_bkg_require_same_elut``. + purpose_penalty : `~astropy.units.Quantity`, optional + Distance penalty for non-``Background``-purpose candidates; defaults to the + ``[Processing]`` config key ``flarelist_bkg_purpose_penalty_days``. + keywords, exclude_keywords : tuple of str, optional + Positive / negative background keywords for the candidate search. + exclude_flare_comment : bool, optional + Drop candidates whose comment references a specific flare id. + + Returns + ------- + BackgroundSelection + """ + if window_past is None: + window_past = CONFIG.getfloat("Processing", "flarelist_bkg_window_past_days", fallback=30.0) * u.day + if window_future is None: + window_future = CONFIG.getfloat("Processing", "flarelist_bkg_window_future_days", fallback=7.0) * u.day + if min_duration is None: + min_duration = CONFIG.getfloat("Processing", "flarelist_bkg_min_duration_s", fallback=1200.0) * u.s + if require_same_elut is None: + require_same_elut = CONFIG.getboolean("Processing", "flarelist_bkg_require_same_elut", fallback=True) + if purpose_penalty is None: + purpose_penalty = CONFIG.getfloat("Processing", "flarelist_bkg_purpose_penalty_days", fallback=0.0) * u.day + if exclude_keywords is None: + exclude_keywords = ( + DEFAULT_BKG_EXCLUDE_KEYWORDS + if CONFIG.getboolean("Processing", "flarelist_bkg_exclude_elevated", fallback=False) + else () + ) + if exclude_flare_comment is None: + exclude_flare_comment = CONFIG.getboolean("Processing", "flarelist_bkg_exclude_flare_comment", fallback=False) + + t = time if isinstance(time, Time) else Time(time) + if rid_lut is None: + rid_lut = RidLutManager.instance.rid_lut + + # ELUT active at the flare time; only enforced when it can be resolved + flare_elut = _elut_id(t) if require_same_elut else None + + candidates = search_background_candidates( + rid_lut, + t, + window_past=window_past, + window_future=window_future, + keywords=keywords, + exclude_keywords=exclude_keywords, + exclude_flare_comment=exclude_flare_comment, + purpose_penalty=purpose_penalty, + ) + + t_jd = t.to_value("jd") + penalty_days = purpose_penalty.to_value(u.day) + + def _valid_to(chosen): + # the selection holds until another candidate's *effective* distance overtakes + # the chosen one's (accounts for the purpose penalty), or the chosen leaves the + # past window if none does. Capped at time + window_future so newly-reachable + # candidates get re-scanned. + pa = 0.0 if chosen.is_background else penalty_days + sa = chosen.start.to_value("jd") + edge = sa + window_past.to_value(u.day) + for j in candidates: + if j.rid == chosen.rid: + continue + pb = 0.0 if j.is_background else penalty_days + tc = _effective_crossover(t_jd, sa, pa, j.start.to_value("jd"), pb) + if tc is not None: + edge = min(edge, tc) + return min(Time(edge, format="jd", scale="utc"), t + window_future) + + for cand in candidates: + if (cand.end - cand.start) < min_duration: + logger.debug(f"bkg candidate rid {cand.rid}: requested duration below {min_duration}; skipping") + continue + if flare_elut is not None and _elut_id(cand.start) != flare_elut: + logger.debug(f"bkg candidate rid {cand.rid}: different ELUT configuration; skipping") + continue + try: + res = fido_client.search( + a.Time(cand.start, cand.end), + a.Instrument.stix, + a.stix.DataProduct.sci_xray_cpd, + a.Level("L1"), + ) + except Exception as e: + logger.warning(f"bkg candidate rid {cand.rid}: CPD search failed: {e}") + continue + if len(res) == 0: + continue + res.filter_for_latest_version() + url_to_path(res) + if "path" not in res.columns: + continue + for path in res["path"]: + if path is None: + continue + if _rid_from_filename(path) != cand.rid: + continue + try: + p = STIXPYProduct(path) + except Exception as e: + logger.warning(f"bkg candidate rid {cand.rid}: could not load {path}: {e}") + continue + rcr = p.data["rcr"] if "rcr" in p.data.colnames else None + if rcr is None or np.any(np.asarray(rcr) != 0): + logger.debug(f"bkg candidate rid {cand.rid}: attenuator in (rcr != 0); skipping {path}") + continue + valid_to = _valid_to(cand) + logger.info(f"selected background file {path} (rid {cand.rid}, {cand.side}) for {t.isot}") + return BackgroundSelection(path=Path(str(path)), rid=cand.rid, valid_from=t, valid_to=valid_to) + + later = [c.start for c in candidates if c.start > t] + valid_to = min(later) if later else t + window_future + logger.info(f"no usable background file found for {t.isot}") + return BackgroundSelection(path=None, rid=-1, valid_from=t, valid_to=valid_to) + + class FlareListManager: @property def flarelist(self): @@ -40,6 +369,229 @@ def flarelistname(self): def productCls(self): return self._product_cls + def _build_ql_month_timeline(self, *, start, end, fido_client, data_product, n_detectors, track_energy): + """Search + load all L1 QL files of ``data_product`` for ``[start, end)`` and + stack them into one time-sorted timeline with a per-bin count-rate column. + + Each daily file is opened exactly once (used for both the energy-table + lookup and the counts), so the whole month costs one open per day. + + Returns + ------- + (timeline, energy, date_to_eidx) + ``timeline`` : QTable with columns ``time``, ``counts`` (raw ``ct``, + ``(N, 5)``), ``counts_rate`` (``(N, 5)``) and, for the lightcurve, + ``rcr``. ``energy`` and ``date_to_eidx`` are only populated when + ``track_energy`` is True. + """ + energy = QTable() + energy_look_up = {} + date_to_eidx = {} + daily_tables = [] + + try: + res = fido_client.search(a.Time(start, end), a.Instrument.stix, data_product, a.Level("L1")) + if len(res) > 0: + res.filter_for_latest_version() + url_to_path(res) + except Exception as e: + logger.error(f"error searching L1 QL {data_product} files for month {start}: {e}") + return QTable(), energy, date_to_eidx + + if len(res) == 0 or "path" not in res.columns: + return QTable(), energy, date_to_eidx + + for path in res["path"]: + if path is None: + continue + try: + p = STIXPYProduct(path) + except Exception as e: + logger.warning(f"could not load QL product {path}: {e}") + continue + + energies = getattr(p, "_energies", None) + if energies is None: + logger.warning(f"QL product {path} has no energies table; skipping") + continue + + counts = p.data["counts"] + if counts.shape[1] != 5: + logger.warning(f"QL product {path} has {counts.shape[1]} channels (expected 5); skipping") + continue + + energy_delta = energies["e_high"] - energies["e_low"] + rate = compute_ql_count_rate( + counts, p.data["timedel"], p.data["triggers"], energy_delta, n_detectors=n_detectors + ) + + daily = QTable() + daily["time"] = p.data["time"] + daily["counts"] = counts + daily["counts_rate"] = rate + if "rcr" in p.data.colnames: + daily["rcr"] = np.asarray(p.data["rcr"]).astype(np.int16) + daily_tables.append(daily) + + if track_energy: + e_sub = QTable() + e_sub["channel"] = energies["channel"] + e_sub["e_low"] = energies["e_low"] + e_sub["e_high"] = energies["e_high"] + e_hash = frozenset(pd.core.util.hashing.hash_array(e_sub.as_array())) + if e_hash not in energy_look_up: + e_idx = len(energy_look_up) + energy_look_up[e_hash] = e_idx + e_sub["index"] = Column(e_idx, description="energy edge table index", dtype=np.int8) + energy = vstack([energy, e_sub]) + if e_idx > 0: + logger.warning(f"multiple energy ql-lc tables found for month {start}") + eidx = energy_look_up[e_hash] + for d in {t.to_datetime().date() for t in p.data["time"]}: + date_to_eidx.setdefault(d, eidx) + logger.info(f"loaded QL product {path} with {len(p.data)} bins and {len(energies)} energy channels") + + return build_month_timeline(daily_tables), energy, date_to_eidx + + def add_lc_bkg_columns(self, data, *, start, end, fido_client): + """Populate LC/BKG peak counts + rates, RCR and ``att_in`` on ``data`` from the + real L1 QL lightcurve + background products, and return the energy QTable. + + ``data`` must already carry ``flare_id`` and the astropy ``Time`` columns + ``start_UTC`` / ``end_UTC`` / ``peak_UTC``. Columns are added in place: + ``lc_peak``, ``lc_peak_rate``, ``lc_bgk_peak``, ``lc_bgk_peak_rate``, + ``rcr_at_peak``, ``rcr_max``, ``att_in``, ``energy_index``. + """ + n = len(data) + tol = CONFIG.getfloat("Processing", "flarelist_peak_max_dist_s", fallback=60.0) * u.s + + lc_timeline, energy, date_to_eidx = self._build_ql_month_timeline( + start=start, + end=end, + fido_client=fido_client, + data_product=a.stix.DataProduct.ql_lightcurve, + n_detectors=16, + track_energy=True, + ) + bkg_timeline, _, _ = self._build_ql_month_timeline( + start=start, + end=end, + fido_client=fido_client, + data_product=a.stix.DataProduct.ql_background, + n_detectors=1, + track_energy=False, + ) + + lc_peak = np.zeros((n, 5), dtype=np.int64) + lc_peak_rate = np.zeros((n, 5), dtype=np.float64) + lc_bgk_peak = np.zeros((n, 5), dtype=np.int64) + lc_bgk_peak_rate = np.zeros((n, 5), dtype=np.float64) + rcr_at_peak = np.full(n, -1, dtype=np.int8) + rcr_max = np.full(n, -1, dtype=np.int8) + energy_index = np.zeros(n, dtype=np.int8) + + rate_unit = u.ct / (u.s * u.keV) + lc_has = len(lc_timeline) > 0 + bkg_has = len(bkg_timeline) > 0 + if not lc_has: + logger.warning(f"No L1 QL lightcurve data found for month {start}") + if not bkg_has: + logger.warning(f"No L1 QL background data found for month {start}") + + for i, row in enumerate(data): + peak = row["peak_UTC"] + fid = row["flare_id"] + if lc_has: + j = nearest_bin_index(lc_timeline["time"], peak, tol) + if j is not None: + lc_peak[i] = lc_timeline["counts"][j].to_value(u.ct) + lc_peak_rate[i] = lc_timeline["counts_rate"][j].to_value(rate_unit) + rcr_at_peak[i] = int(lc_timeline["rcr"][j]) + rcr_max[i] = max_rcr_in_window( + lc_timeline["time"], + lc_timeline["rcr"], + row["start_UTC"], + row["end_UTC"], + fallback=int(rcr_at_peak[i]), + ) + energy_index[i] = date_to_eidx.get(peak.to_datetime().date(), 0) + else: + logger.warning(f"flare {fid}: no LC bin within {tol} of peak {peak.isot}") + if bkg_has: + k = nearest_bin_index(bkg_timeline["time"], peak, tol) + if k is not None: + lc_bgk_peak[i] = bkg_timeline["counts"][k].to_value(u.ct) + lc_bgk_peak_rate[i] = bkg_timeline["counts_rate"][k].to_value(rate_unit) + else: + logger.warning(f"flare {fid}: no BKG bin within {tol} of peak {peak.isot}") + + data["lc_peak"] = Column( + lc_peak * u.ct, + description="raw counts at the L1 QL lightcurve bin nearest the flare peak (5 energy channels)", + dtype=np.int64, + ) + data["lc_peak_rate"] = Column( + lc_peak_rate * rate_unit, + description="livetime-corrected count rate at the peak lightcurve bin (5 energy channels)", + ) + data["lc_bgk_peak"] = Column( + lc_bgk_peak * u.ct, + description="raw background counts at the L1 QL background bin nearest the flare peak (5 energy channels)", + dtype=np.int64, + ) + data["lc_bgk_peak_rate"] = Column( + lc_bgk_peak_rate * rate_unit, + description="livetime-corrected background count rate at the peak bin (5 energy channels)", + ) + data["rcr_at_peak"] = Column( + rcr_at_peak, description="rate control regime at the peak bin (>0 attenuator in)", dtype=np.int8 + ) + data["rcr_max"] = Column( + rcr_max, description="max rate control regime over the flare start..end window", dtype=np.int8 + ) + data["att_in"] = Column(rcr_max > 0, description="was attenuator in during flare (rcr_max > 0)") + data["energy_index"] = Column(energy_index, description="energy band index", dtype=np.int8) + + return energy + + def add_background_file_column(self, data, *, fido_client): + """Add ``bkg_file`` / ``bkg_rid`` columns: the best quiet-time background + CPD file for each flare peak. + + ``data`` rows are assumed to be peak-time ascending (as produced by the + source flare list), so each per-time background search + (:func:`find_background_file_for_time`) is cached with its validity + interval and only re-run once a flare peak crosses out of that period. + """ + n = len(data) + bkg_files = [""] * n + bkg_rids = np.full(n, -1, dtype=np.int64) + + primer = "" + baseurl = getattr(fido_client, "baseurl", None) + datapath = getattr(fido_client, "datapath", None) + if baseurl is not None and datapath is not None: + primer = baseurl.replace(datapath, "") + primer = primer[7:] if primer.startswith("file://") else primer + + selection = None + searches = 0 + for i, row in enumerate(data): + peak = row["peak_UTC"] + if selection is None or not (selection.valid_from <= peak <= selection.valid_to): + selection = find_background_file_for_time(peak, fido_client=fido_client) + searches += 1 + if selection.path is not None: + bkg_files[i] = str(selection.path).replace(primer, "") + bkg_rids[i] = selection.rid + + data["bkg_file"] = Column(bkg_files, description="path to the quiet-time background CPD file") + data["bkg_rid"] = Column( + bkg_rids, description="BSD request id of the selected background file (-1 if none)", dtype=np.int64 + ) + logger.info(f"background file search ran {searches}x for {n} flares") + return data + class SCFlareListManager(FlareListManager, metaclass=Singleton): """Manages a local copy of the flarelist provided by STIXCore or runs the flare detection @@ -188,15 +740,19 @@ 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") data["GOES_class"] = Column( mt["GOES_class"].astype(str), description="GOES class of the GOES XRS data at time of flare" @@ -220,15 +776,15 @@ def get_data(self, *, start, end, fido_client): "flare isn't visible to Earth", ) data["goes_min_flux_est"] = Column( - mt["goes_estimated_min_flux"].astype(float) * u.W / u.m**2, + (10 ** mt["goes_estimated_min_flux"].astype(float)) * u.W / u.m**2, description="min GOES flux estimate derived from STIX data", ) data["goes_max_flux_est"] = Column( - mt["goes_estimated_max_flux"].astype(float) * u.W / u.m**2, + (10 ** mt["goes_estimated_max_flux"].astype(float)) * u.W / u.m**2, description="max GOES flux estimate derived from STIX data", ) data["goes_mean_flux_est"] = Column( - mt["goes_estimated_mean_flux"].astype(float) * u.W / u.m**2, + (10 ** mt["goes_estimated_mean_flux"].astype(float)) * u.W / u.m**2, description="mean GOES flux estimate derived from STIX data", ) @@ -277,7 +833,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,15 +995,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"]] - 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") + + 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["GOES_class"] = Column( mt["GOES_class"].astype(str), description="GOES class of the GOES XRS data at time of flare" @@ -471,15 +1040,15 @@ def get_data(self, *, start, end, fido_client): "flare isn't visible to Earth", ) data["goes_min_flux_est"] = Column( - mt["goes_estimated_min_flux"].astype(float) * u.W / u.m**2, + (10 ** mt["goes_estimated_min_flux"].astype(float)) * u.W / u.m**2, description="min GOES flux estimate derived from STIX data", ) data["goes_max_flux_est"] = Column( - mt["goes_estimated_max_flux"].astype(float) * u.W / u.m**2, + (10 ** mt["goes_estimated_max_flux"].astype(float)) * u.W / u.m**2, description="max GOES flux estimate derived from STIX data", ) data["goes_mean_flux_est"] = Column( - mt["goes_estimated_mean_flux"].astype(float) * u.W / u.m**2, + (10 ** mt["goes_estimated_mean_flux"].astype(float)) * u.W / u.m**2, description="mean GOES flux estimate derived from STIX data", ) @@ -490,24 +1059,13 @@ def get_data(self, *, start, end, fido_client): # description="coarse flare location in y direction provided by" # "onboard algorithm. (0,0) represents disk center") - data["lc_peak"] = Column( - ( - np.vstack( - ( - mt["LC0_PEAK_COUNTS_4S"].value, - mt["LC1_PEAK_COUNTS_4S"].value, - mt["LC2_PEAK_COUNTS_4S"].value, - mt["LC3_PEAK_COUNTS_4S"].value, - mt["LC4_PEAK_COUNTS_4S"].value, - ) - ).T - * u.ct - ).astype(int), - description="counts in 4s peak window from quicklook lightcurve", - dtype=np.int64, + # background estimates carried over from the source flare list (distinct from the + # QL-derived at-peak background added below) + data["bkg_baseline"] = Column( + mt["LC0_BKG"].astype(float) * u.ct, + description="median value of the fitted baseline", ) - - data["lc_bgk_peak"] = Column( + data["bkg_quiet_period"] = Column( ( np.vstack( ( @@ -519,55 +1077,20 @@ def get_data(self, *, start, end, fido_client): ) ).T * u.ct - ).astype(int), - description="background counts in 4s peak windowfrom quicklook lightcurve", + ).astype(np.int64), + description="background counts per QL energy channel, median value for the most recent quiet period", dtype=np.int64, ) - data["energy_index"] = Column(0, description="energy band index", dtype=np.int8) + # LC/BKG peak counts+rates, RCR, att_in and energy_index are derived from the + # real L1 QL lightcurve + background products (the source CSV values are + # unreliable) using one monthly timeline per product built once. + energy = self.add_lc_bkg_columns(data, start=start, end=end, fido_client=fido_client) - data.add_index("flare_id") - - # add energy axis for the lightcurve peek 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. - # So we do the lookup also just grouped by peak day in order to save file lookups + # select the best quiet-time background data file for each flare peak + self.add_background_file_column(data, fido_client=fido_client) - energy_look_up = {} - data["peak_day"] = [d.datetime.day for d in data["peak_UTC"]] - data_by_day = data.group_by("peak_day") - - for day, flares in zip(data_by_day.groups.keys, data_by_day.groups): - time = flares["peak_UTC"][0] - lc_data = fido_client.search(a.Time(time, time), a.Instrument.stix, a.stix.DataProduct.ql_lightcurve) - lc_data.filter_for_latest_version() - url_to_path(lc_data) - - if len(lc_data) == 0: - logger.warning(f"No lightcurve data found for flare at time {time}") - continue - lc = Product(lc_data["path"][0]) - - energy_table_hash = frozenset(pd.core.util.hashing.hash_array(lc.energies.as_array())) - - # add the energy table to the energy table list if not already - # present and define a new index number - if energy_table_hash not in energy_look_up: - e_idx = len(energy_look_up.keys()) - energy_look_up[energy_table_hash] = e_idx - lc.energies["index"] = Column(e_idx, description="energy edge table index", dtype=np.int8) - energy = vstack([energy, lc.energies]) - if e_idx > 0: - logger.warning(f"multiple energy ql-lc tables found for month {start}") - - # add the energy index to the flare data to all flares of the same day - # https://docs.astropy.org/en/latest/table/modify_table.html#caveats - replace = data.loc[flares["flare_id"]] - replace["energy_index"] = energy_look_up[energy_table_hash] - data.loc[flares["flare_id"]] = replace - - del data["peak_day"] + data.add_index("flare_id") return data, control, energy diff --git a/stixcore/io/RidLutManager.py b/stixcore/io/RidLutManager.py index bd1f9d26..154956b8 100644 --- a/stixcore/io/RidLutManager.py +++ b/stixcore/io/RidLutManager.py @@ -1,23 +1,163 @@ +import re import sys import time import tempfile import urllib.request from datetime import date, datetime, timedelta +from collections import namedtuple import numpy as np +import astropy.units as u from astropy.io import ascii from astropy.table import Table from astropy.table.operations import unique, vstack +from astropy.time import Time from stixcore.config.config import CONFIG from stixcore.util.logging import get_logger from stixcore.util.singleton import Singleton -__all__ = ["RidLutManager"] +__all__ = ["RidLutManager", "BackgroundCandidate", "search_background_candidates"] logger = get_logger(__name__) +#: Keywords (case-insensitive) that mark a BSD request as a background/quiet +#: observation in the descriptive columns of the RID LUT. +DEFAULT_BKG_KEYWORDS = ("bkg", "quiet", "background", "non-flaring") + +#: Negative keywords: a candidate whose descriptive text contains any of these is +#: rejected (e.g. "elevated" background is not a clean quiet baseline). +DEFAULT_BKG_EXCLUDE_KEYWORDS = ("elevated",) + +#: Matches a specific flare id referenced in a comment, e.g. "for Flare 2309081508". +#: Such rows are flare data requests, not dedicated backgrounds. +_FLARE_REF_RE = re.compile(r"flare\s*\d{3,}") + +#: A single background-request candidate from the RID LUT. ``start``/``end``/``mid`` +#: are `~astropy.time.Time`, ``side`` is ``"past"``/``"future"`` relative to the query +#: time, and ``is_background`` is True when ``purpose == "Background"`` (a clean +#: background request, preferred over subject-only keyword matches). +BackgroundCandidate = namedtuple("BackgroundCandidate", ["rid", "start", "end", "mid", "side", "is_background"]) + + +def _col_as_lower_str(tbl, name): + """Return column ``name`` of ``tbl`` as a lower-cased ``str`` numpy array.""" + col = tbl[name] + try: + col = col.filled("") + except (AttributeError, TypeError): + pass + return np.char.lower(np.asarray(col, dtype=str)) + + +def search_background_candidates( + rid_lut, + time, + *, + window_past, + window_future, + keywords=DEFAULT_BKG_KEYWORDS, + exclude_keywords=DEFAULT_BKG_EXCLUDE_KEYWORDS, + exclude_flare_comment=True, + purpose_penalty=1.0 * u.day, +): + """Find background-request candidates in a RID LUT near ``time``. + + Rows are recognised as background requests by a case-insensitive keyword + match (``keywords``) over the ``subject``/``purpose``/``comment`` columns, + then filtered and ranked: + + * **exclude keywords** — a row whose text contains any of ``exclude_keywords`` + (default ``"elevated"``) is dropped. + * **flare-id comment** — if ``exclude_flare_comment`` a row whose comment + references a specific flare id (e.g. "for Flare 2309081508") is dropped. + * **ranking** — **nearest-in-time first** by ``|time - start|`` (past preferred + on a tie), but a request with ``purpose == "Background"`` is preferred over a + subject-only keyword match unless the latter is more than ``purpose_penalty`` + closer. This is implemented as an *effective* distance + ``|time - start| + purpose_penalty`` for non-Background rows. + + ``window_past`` / ``window_future`` remain separate bounds on how far a + candidate's start may lie in each direction. + + Parameters + ---------- + rid_lut : `~astropy.table.Table` + The RID LUT (as produced by `RidLutManager.read_rid_lut`). + time : `~astropy.time.Time` or str + The query time (e.g. a flare peak). + window_past, window_future : `~astropy.units.Quantity` + How far back / forward from ``time`` a candidate's start may lie. + keywords, exclude_keywords : tuple of str, optional + Positive / negative case-insensitive keywords. + exclude_flare_comment : bool, optional + Drop rows whose comment references a specific flare id. + purpose_penalty : `~astropy.units.Quantity`, optional + Distance penalty added to non-``Background``-purpose candidates so a clean + Background request wins unless a subject-only match is clearly closer. + + Returns + ------- + list of BackgroundCandidate + Ranked by ascending effective distance (past preferred on a tie). Empty + if the LUT has no matching rows in range. + """ + if rid_lut is None or len(rid_lut) == 0: + return [] + t = time if isinstance(time, Time) else Time(time) + + subj = _col_as_lower_str(rid_lut, "subject") + purp = _col_as_lower_str(rid_lut, "purpose") + comm = _col_as_lower_str(rid_lut, "comment") + haystack = np.char.add(np.char.add(subj, " "), np.char.add(purp, np.char.add(" ", comm))) + + mask = np.zeros(len(rid_lut), dtype=bool) + for kw in keywords: + mask |= np.char.find(haystack, kw.lower()) >= 0 + for kw in exclude_keywords: # negative keywords drop the row + mask &= np.char.find(haystack, kw.lower()) < 0 + if not np.any(mask): + return [] + + sub = rid_lut[mask] + purp_sub = purp[mask] + comm_sub = comm[mask] + starts = Time(np.asarray(sub["start_utc"], dtype=str), format="isot", scale="utc") + durations = np.asarray(sub["duration"], dtype=float) * u.s + ends = starts + durations + mids = starts + durations / 2 + rids = np.asarray(sub["unique_id"]).astype(np.int64) + + wp = window_past.to_value(u.day) + wf = window_future.to_value(u.day) + penalty = purpose_penalty.to_value(u.day) + kept = [] + for i in range(len(sub)): + if exclude_flare_comment and _FLARE_REF_RE.search(comm_sub[i]): + continue + d = (t - starts[i]).to_value(u.day) # > 0 starts before ``time`` (past), < 0 after (future) + if d >= 0: + if d > wp: + continue + side = "past" + else: + if -d > wf: + continue + side = "future" + is_bg = purp_sub[i].strip() == "background" + effective = abs(d) + (0.0 if is_bg else penalty) # prefer Background at equal-ish distance + kept.append( + ( + effective, + 0 if side == "past" else 1, + BackgroundCandidate(int(rids[i]), starts[i], ends[i], mids[i], side, is_bg), + ) + ) + + kept.sort(key=lambda x: (x[0], x[1])) + return [c for _, _, c in kept] + class RidLutManager(metaclass=Singleton): """Manages metadata for BSD requests @@ -78,6 +218,13 @@ def get_reason(self, rid): logger.warning("can't get request purpose: no request founds for rid: {rid}") return "" + def find_background_candidates(self, time, *, window_past, window_future, keywords=DEFAULT_BKG_KEYWORDS): + """Find background-request candidates near ``time`` (see + :func:`search_background_candidates`).""" + return search_background_candidates( + self.rid_lut, time, window_past=window_past, window_future=window_future, keywords=keywords + ) + def get_scaling_factor(self, rid): """Gets the trigger descaling factor connected to the BSD request. diff --git a/stixcore/io/product_processors/fits/processors.py b/stixcore/io/product_processors/fits/processors.py index 90630eea..b843f369 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 @@ -805,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"), @@ -1049,7 +1068,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) primary_hdu = fits.PrimaryHDU() @@ -1126,7 +1146,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) @@ -1139,7 +1160,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/io/tests/test_flarelistmanager.py b/stixcore/io/tests/test_flarelistmanager.py new file mode 100644 index 00000000..95b9dd93 --- /dev/null +++ b/stixcore/io/tests/test_flarelistmanager.py @@ -0,0 +1,609 @@ +from pathlib import Path + +import numpy as np +import pytest + +import astropy.units as u +from astropy.table import QTable, Table +from astropy.time import Time + +import stixcore.io.FlareListManager as flm_mod +from stixcore.io.FlareListManager import ( + BackgroundSelection, + FlareListManager, + build_month_timeline, + compute_ql_count_rate, + find_background_file_for_time, + max_rcr_in_window, + nearest_bin_index, +) +from stixcore.io.RidLutManager import RidLutManager, search_background_candidates + +RATE_UNIT = u.ct / (u.s * u.keV) + + +# --- helpers to build synthetic QL data --------------------------------------- + + +def _energies(): + e = QTable() + e["channel"] = np.arange(5, dtype=np.uint8) + e["e_low"] = [4, 10, 15, 25, 50] * u.keV + e["e_high"] = [10, 15, 25, 50, 84] * u.keV + return e + + +def _ql_data(t0, n, *, rcr=None, counts_value=100, with_rcr=True): + """Build a synthetic QL product ``.data`` QTable on a 4 s grid.""" + times = Time(t0) + np.arange(n) * 4 * u.s + data = QTable() + data["time"] = times + data["timedel"] = np.full(n, 4.0) * u.s + data["triggers"] = np.zeros(n) # zero triggers -> live_frac == 1 (deterministic) + data["counts"] = (np.full((n, 5), counts_value)).astype(int) * u.ct + if with_rcr: + data["rcr"] = np.zeros(n, dtype=np.ubyte) if rcr is None else np.asarray(rcr, dtype=np.ubyte) + return data + + +class FakeProduct: + def __init__(self, data, energies): + self.data = data + self._energies = energies + + +class FakeResponse: + """Minimal stand-in for a StixQueryResponse carrying a ``path`` column.""" + + def __init__(self, paths): + self._paths = list(paths) + + def __len__(self): + return len(self._paths) + + @property + def columns(self): + return ["path"] if self._paths else [] + + def filter_for_latest_version(self): + pass + + def __getitem__(self, key): + assert key == "path" + return self._paths + + +class FakeFido: + def __init__(self, lc_paths, bkg_paths): + self.lc_paths = lc_paths + self.bkg_paths = bkg_paths + + def search(self, time, instrument, data_product, level): + name = getattr(data_product, "value", str(data_product)) + if "background" in name: + return FakeResponse(self.bkg_paths) + return FakeResponse(self.lc_paths) + + +# --- compute_ql_count_rate ---------------------------------------------------- + + +def test_count_rate_units_and_shape(): + data = _ql_data("2024-06-15T12:00:00", 3) + ed = _energies()["e_high"] - _energies()["e_low"] + rate = compute_ql_count_rate(data["counts"], data["timedel"], data["triggers"], ed, n_detectors=16) + assert rate.shape == (3, 5) + assert rate.unit.is_equivalent(RATE_UNIT) + + +def test_count_rate_zero_triggers_is_deterministic(): + # zero triggers => live_frac == 1 => rate == counts / (timedel * energy_delta), + # independent of n_detectors, so LC and BKG agree exactly. + data = _ql_data("2024-06-15T12:00:00", 2, counts_value=200) + ed = _energies()["e_high"] - _energies()["e_low"] + lc = compute_ql_count_rate(data["counts"], data["timedel"], data["triggers"], ed, n_detectors=16) + bkg = compute_ql_count_rate(data["counts"], data["timedel"], data["triggers"], ed, n_detectors=1) + expected = data["counts"] / ((data["timedel"]).reshape(-1, 1) * ed) + assert u.allclose(lc, expected.to(RATE_UNIT)) + assert u.allclose(lc, bkg) + + +def test_count_rate_more_detectors_gives_smaller_rate_when_triggers_nonzero(): + # nonzero triggers: larger n_detectors -> lower trigger rate -> higher live_frac + # -> larger denominator -> smaller count rate. + data = _ql_data("2024-06-15T12:00:00", 1) + data["triggers"] = np.array([5000.0]) + ed = _energies()["e_high"] - _energies()["e_low"] + lc = compute_ql_count_rate(data["counts"], data["timedel"], data["triggers"], ed, n_detectors=16) + bkg = compute_ql_count_rate(data["counts"], data["timedel"], data["triggers"], ed, n_detectors=1) + assert np.all(lc.to_value(RATE_UNIT) < bkg.to_value(RATE_UNIT)) + + +# --- build_month_timeline ----------------------------------------------------- + + +def test_build_timeline_sorts_and_dedups_overlap(): + d1 = _ql_data("2024-06-15T12:00:00", 3) # 12:00:00, :04, :08 + d2 = _ql_data("2024-06-15T12:00:08", 3) # :08 (dup), :12, :16 + timeline = build_month_timeline([d2, d1]) # pass out of order + t = timeline["time"] + assert len(timeline) == 5 # 6 bins minus 1 duplicate at :08 + assert np.all(np.diff(t.jd) > 0) # strictly increasing + + +def test_build_timeline_empty(): + assert len(build_month_timeline([])) == 0 + assert len(build_month_timeline([None, QTable()])) == 0 + + +# --- nearest_bin_index -------------------------------------------------------- + + +def test_nearest_bin_exact(): + times = Time("2024-06-15T12:00:00") + np.arange(5) * 4 * u.s + assert nearest_bin_index(times, times[2], 2 * u.s) == 2 + + +def test_nearest_bin_within_tol(): + times = Time("2024-06-15T12:00:00") + np.arange(5) * 4 * u.s + target = times[3] + 1 * u.s + assert nearest_bin_index(times, target, 2 * u.s) == 3 + + +def test_nearest_bin_beyond_tol_returns_none(): + times = Time("2024-06-15T12:00:00") + np.arange(5) * 4 * u.s + target = times[-1] + 1 * u.h + assert nearest_bin_index(times, target, 60 * u.s) is None + + +def test_nearest_bin_empty_returns_none(): + assert nearest_bin_index(Time([], format="isot"), Time("2024-06-15T12:00:00"), 60 * u.s) is None + + +# --- max_rcr_in_window -------------------------------------------------------- + + +def test_max_rcr_over_window(): + times = Time("2024-06-15T12:00:00") + np.arange(6) * 4 * u.s + rcr = np.array([0, 0, 1, 2, 1, 0]) + assert max_rcr_in_window(times, rcr, times[1], times[4], fallback=-1) == 2 + + +def test_max_rcr_peak_only_window(): + times = Time("2024-06-15T12:00:00") + np.arange(6) * 4 * u.s + rcr = np.array([0, 0, 1, 2, 1, 0]) + assert max_rcr_in_window(times, rcr, times[0], times[0], fallback=-1) == 0 + + +def test_max_rcr_empty_window_returns_fallback(): + times = Time("2024-06-15T12:00:00") + np.arange(6) * 4 * u.s + rcr = np.array([0, 0, 1, 2, 1, 0]) + before = Time("2024-06-15T11:00:00") + assert max_rcr_in_window(times, rcr, before, before, fallback=7) == 7 + + +# --- add_lc_bkg_columns (integration, monkeypatched) -------------------------- + + +@pytest.fixture +def flare_data(): + # bin grid starts 2024-06-15T12:00:00, 4 s cadence; bins 5-7 attenuated (rcr=1) + peaks = Time( + [ + "2024-06-15T12:00:20", # bin 5 (attenuated), window covers attenuated bins + "2024-06-15T12:00:00", # bin 0, no attenuation + "2024-06-15T13:00:00", # far away -> beyond tolerance + ] + ) + data = QTable() + data["flare_id"] = [1, 2, 3] + data["start_UTC"] = peaks - 8 * u.s + data["end_UTC"] = peaks + 8 * u.s + data["peak_UTC"] = peaks + return data + + +def _patch_fido(monkeypatch): + n = 20 + rcr = np.zeros(n, dtype=np.ubyte) + rcr[5:8] = 1 + lc = FakeProduct(_ql_data("2024-06-15T12:00:00", n, rcr=rcr), _energies()) + bkg = FakeProduct(_ql_data("2024-06-15T12:00:00", n, with_rcr=False), _energies()) + + products = {"lc": lc, "bkg": bkg} + monkeypatch.setattr("stixcore.io.FlareListManager.STIXPYProduct", lambda path: products[path]) + return FakeFido(lc_paths=["lc"], bkg_paths=["bkg"]) + + +def test_add_lc_bkg_columns(flare_data, monkeypatch): + fido = _patch_fido(monkeypatch) + from datetime import date + + energy = FlareListManager().add_lc_bkg_columns( + flare_data, start=date(2024, 6, 1), end=date(2024, 7, 1), fido_client=fido + ) + + # shapes / units (QTable stores unit-bearing columns as float64 Quantity, + # same as the pre-existing lc_peak column; values stay integral counts) + assert flare_data["lc_peak"].shape == (3, 5) + assert flare_data["lc_peak"].unit == u.ct + assert np.all(flare_data["lc_peak"].value == np.round(flare_data["lc_peak"].value)) + assert flare_data["lc_peak_rate"].unit.is_equivalent(RATE_UNIT) + assert flare_data["lc_bgk_peak"].shape == (3, 5) + assert flare_data["lc_bgk_peak"].unit == u.ct + assert flare_data["att_in"].dtype == bool + assert flare_data["energy_index"].dtype == np.int8 + + # rcr semantics + assert np.all(flare_data["rcr_max"] >= flare_data["rcr_at_peak"]) + assert list(flare_data["att_in"]) == [True, False, False] + assert flare_data["rcr_max"][0] == 1 # window covers attenuated bins + assert flare_data["rcr_at_peak"][0] == 1 + assert flare_data["rcr_at_peak"][1] == 0 + + # far flare -> beyond tolerance -> zero filled, rcr -1 + assert np.all(flare_data["lc_peak"][2].to_value(u.ct) == 0) + assert flare_data["rcr_at_peak"][2] == -1 + assert flare_data["rcr_max"][2] == -1 + + # counts pulled from the real timeline for in-range flares + assert np.all(flare_data["lc_peak"][0].to_value(u.ct) == 100) + assert np.all(flare_data["lc_bgk_peak"][0].to_value(u.ct) == 100) + + # returned energy table schema + assert set(energy.colnames) == {"channel", "e_low", "e_high", "index"} + assert len(energy) == 5 + + +def test_add_lc_bkg_columns_no_files(flare_data, monkeypatch): + from datetime import date + + fido = FakeFido(lc_paths=[], bkg_paths=[]) + energy = FlareListManager().add_lc_bkg_columns( + flare_data, start=date(2024, 6, 1), end=date(2024, 7, 1), fido_client=fido + ) + + assert np.all(flare_data["lc_peak"].to_value(u.ct) == 0) + assert np.all(flare_data["lc_bgk_peak"].to_value(u.ct) == 0) + assert np.all(flare_data["rcr_at_peak"] == -1) + assert np.all(flare_data["rcr_max"] == -1) + assert not np.any(flare_data["att_in"]) + assert len(energy) == 0 + + +# --- background candidate search (RID LUT) ------------------------------------ + + +def _bkg_lut(rows): + """Build a minimal RID LUT ``Table`` from ``(rid, start_iso, duration_s)`` rows.""" + t = Table() + t["unique_id"] = [r[0] for r in rows] + t["start_utc"] = [r[1] for r in rows] + t["duration"] = [r[2] for r in rows] + t["subject"] = ["BKG quiet"] * len(rows) + t["purpose"] = ["Background"] * len(rows) + t["comment"] = [""] * len(rows) + return t + + +def test_candidates_nearest_in_time_regardless_of_side(): + lut = _bkg_lut( + [ + (3001, "2023-06-13T00:00:00", 3600), # past 2 d + (3003, "2023-06-18T00:00:00", 3600), # future 3 d + (3002, "2023-06-05T00:00:00", 3600), # past 10 d + (3004, "2023-04-15T00:00:00", 3600), # past 61 d -> outside 30 d window + ] + ) + t = Time("2023-06-15T00:00:00") + cands = search_background_candidates(lut, t, window_past=30 * u.day, window_future=7 * u.day) + # nearest start first regardless of side: 2 d past, 3 d future, 10 d past + assert [c.rid for c in cands] == [3001, 3003, 3002] + assert [c.side for c in cands] == ["past", "future", "past"] + + # widening the past window pulls in the legacy request, ranked by its (large) distance + cands90 = search_background_candidates(lut, t, window_past=90 * u.day, window_future=7 * u.day) + assert [c.rid for c in cands90] == [3001, 3003, 3002, 3004] + + +def test_candidates_tie_prefers_past(): + # equal distance past vs future -> past wins the tie + lut = _bkg_lut([(3001, "2023-06-10T00:00:00", 3600), (3003, "2023-06-20T00:00:00", 3600)]) + t = Time("2023-06-15T00:00:00") + cands = search_background_candidates(lut, t, window_past=30 * u.day, window_future=30 * u.day) + assert [c.rid for c in cands] == [3001, 3003] + assert cands[0].side == "past" + + +def _bkg_lut_full(rows): + """RID LUT from ``(rid, start_iso, duration_s, subject, purpose, comment)`` rows.""" + t = Table() + t["unique_id"] = [r[0] for r in rows] + t["start_utc"] = [r[1] for r in rows] + t["duration"] = [r[2] for r in rows] + t["subject"] = [r[3] for r in rows] + t["purpose"] = [r[4] for r in rows] + t["comment"] = [r[5] for r in rows] + return t + + +def test_candidates_exclude_elevated(): + lut = _bkg_lut_full( + [ + (1, "2023-06-14T00:00:00", 3600, "BKG elevated", "Background", ""), # closer but excluded + (2, "2023-06-10T00:00:00", 3600, "BKG quiet", "Background", ""), + ] + ) + cands = search_background_candidates( + lut, Time("2023-06-15T00:00:00"), window_past=30 * u.day, window_future=7 * u.day + ) + assert [c.rid for c in cands] == [2] + + +def test_candidates_exclude_flare_comment(): + lut = _bkg_lut_full( + [ + ( + 1, + "2023-06-14T00:00:00", + 3600, + "non-flaring AR?", + "Solar Flare", + "CL1 data request for Flare 2309081508", + ), # closer but flare-referenced + (2, "2023-06-10T00:00:00", 3600, "BKG quiet", "Background", ""), + ] + ) + cands = search_background_candidates( + lut, Time("2023-06-15T00:00:00"), window_past=30 * u.day, window_future=7 * u.day + ) + assert [c.rid for c in cands] == [2] + # keeping flare-referenced rows brings the closer one back to the front + keep = search_background_candidates( + lut, Time("2023-06-15T00:00:00"), window_past=30 * u.day, window_future=7 * u.day, exclude_flare_comment=False + ) + assert [c.rid for c in keep] == [1, 2] + + +def test_candidates_prefer_background_purpose(): + P = 1.0 * u.day + t = Time("2023-06-15T00:00:00") + # subject-only match 0.5 d closer than the Background request -> Background still wins (within penalty) + lut1 = _bkg_lut_full( + [ + (1, "2023-06-13T00:00:00", 3600, "BKG quiet", "Background", ""), # 2.0 d -> eff 2.0 + (2, "2023-06-13T12:00:00", 3600, "quiet region", "obs", ""), # 1.5 d -> eff 2.5 + ] + ) + c1 = search_background_candidates(lut1, t, window_past=30 * u.day, window_future=7 * u.day, purpose_penalty=P) + assert [c.rid for c in c1] == [1, 2] + assert c1[0].is_background + + # subject-only match clearly closer (>penalty) -> it wins + lut2 = _bkg_lut_full( + [ + (1, "2023-06-13T00:00:00", 3600, "BKG quiet", "Background", ""), # 2.0 d -> eff 2.0 + (2, "2023-06-14T12:00:00", 3600, "quiet region", "obs", ""), # 0.5 d -> eff 1.5 + ] + ) + c2 = search_background_candidates(lut2, t, window_past=30 * u.day, window_future=7 * u.day, purpose_penalty=P) + assert [c.rid for c in c2] == [2, 1] + assert not c2[0].is_background + + +def test_candidates_future_window_excludes_far_future(): + lut = _bkg_lut([(4001, "2023-06-25T00:00:00", 3600)]) # 10 d in the future + t = Time("2023-06-15T00:00:00") + assert search_background_candidates(lut, t, window_past=30 * u.day, window_future=7 * u.day) == [] + got = search_background_candidates(lut, t, window_past=30 * u.day, window_future=14 * u.day) + assert [c.rid for c in got] == [4001] + + +def test_candidates_keyword_match_only(): + lut = Table() + lut["unique_id"] = [1, 2] + lut["start_utc"] = ["2023-06-10T00:00:00", "2023-06-11T00:00:00"] + lut["duration"] = [3600, 3600] + lut["subject"] = ["Solar Flare", "some quiet interval"] # only row 2 matches + lut["purpose"] = ["flare", "obs"] + lut["comment"] = ["", ""] + cands = search_background_candidates( + lut, Time("2023-06-15T00:00:00"), window_past=30 * u.day, window_future=7 * u.day + ) + assert [c.rid for c in cands] == [2] + + +def test_candidates_empty_lut(): + assert ( + search_background_candidates(Table(), Time("2023-06-15"), window_past=30 * u.day, window_future=7 * u.day) == [] + ) + + +def test_ridlutmanager_singleton_find_background_candidates(): + # exercises the method against the shipped test LUT (rows 3001-3004) + cands = RidLutManager.instance.find_background_candidates( + Time("2023-06-15T00:00:00"), window_past=30 * u.day, window_future=7 * u.day + ) + # 3001 (5 d past) and 3003 (5 d future) are equidistant -> past wins tie, then 3002 (14 d past) + assert [c.rid for c in cands] == [3001, 3003, 3002] + + +# --- find_background_file_for_time -------------------------------------------- + + +def _cpd_name(rid): + return f"solo_L1_stix-sci-xray-cpd_20230101T000000-20230101T010000_V01_{rid:010d}-00001.fits" + + +class FakeCpdFido: + """Returns the same CPD file list for every search; the rid-in-filename filter + inside ``find_background_file_for_time`` selects the per-candidate file.""" + + def __init__(self, paths): + self._paths = list(paths) + self.searches = 0 + + def search(self, time, instrument, data_product, level): + self.searches += 1 + return FakeResponse(self._paths) + + +def _patch_cpd_products(monkeypatch, spec): + """``spec``: {rid: rcr_array_or_None}. Builds CPD filenames + FakeProducts and + monkeypatches STIXPYProduct to resolve them.""" + products = {} + paths = [] + for rid, rcr in spec.items(): + name = _cpd_name(rid) + paths.append(name) + with_rcr = rcr is not None + products[name] = FakeProduct(_ql_data("2023-06-10T00:00:00", 5, rcr=rcr, with_rcr=with_rcr), _energies()) + monkeypatch.setattr("stixcore.io.FlareListManager.STIXPYProduct", lambda path: products[path]) + return FakeCpdFido(paths) + + +def test_find_bkg_picks_closest_usable(monkeypatch): + # 3001 starts 5 d before, 3002 starts 7 d after -> 3001 is nearest in time + lut = _bkg_lut([(3001, "2023-06-10T00:00:00", 3600), (3002, "2023-06-22T00:00:00", 3600)]) + fido = _patch_cpd_products(monkeypatch, {3001: np.zeros(5), 3002: np.zeros(5)}) + t = Time("2023-06-15T00:00:00") + sel = find_background_file_for_time( + t, fido_client=fido, rid_lut=lut, min_duration=600 * u.s, require_same_elut=False + ) + assert sel.rid == 3001 # nearest-in-time wins + assert sel.path == Path(_cpd_name(3001)) + assert sel.valid_from == t + # valid until the midpoint to the next later start: (2023-06-10 + 2023-06-22) / 2 = 2023-06-16 + assert sel.valid_to == Time("2023-06-16T00:00:00") + + +def test_find_bkg_skips_attenuator_in(monkeypatch): + # 3001 (closer past) has attenuator in (rcr>0) -> fall through to 3002 + lut = _bkg_lut([(3001, "2023-06-10T00:00:00", 3600), (3002, "2023-06-05T00:00:00", 3600)]) + fido = _patch_cpd_products(monkeypatch, {3001: np.ones(5), 3002: np.zeros(5)}) + sel = find_background_file_for_time( + Time("2023-06-15T00:00:00"), + fido_client=fido, + rid_lut=lut, + min_duration=600 * u.s, + require_same_elut=False, + ) + assert sel.rid == 3002 + + +def test_find_bkg_skips_too_short(monkeypatch): + # 3001 requested only 300 s (< min_duration) -> skipped without even a search + lut = _bkg_lut([(3001, "2023-06-10T00:00:00", 300), (3002, "2023-06-05T00:00:00", 3600)]) + fido = _patch_cpd_products(monkeypatch, {3001: np.zeros(5), 3002: np.zeros(5)}) + sel = find_background_file_for_time( + Time("2023-06-15T00:00:00"), + fido_client=fido, + rid_lut=lut, + min_duration=600 * u.s, + require_same_elut=False, + ) + assert sel.rid == 3002 + + +def test_find_bkg_future_fallback(monkeypatch): + # no past candidate -> the future one is used + lut = _bkg_lut([(3003, "2023-06-18T00:00:00", 3600)]) + fido = _patch_cpd_products(monkeypatch, {3003: np.zeros(5)}) + sel = find_background_file_for_time( + Time("2023-06-15T00:00:00"), + fido_client=fido, + rid_lut=lut, + min_duration=600 * u.s, + require_same_elut=False, + ) + assert sel.rid == 3003 + + +def test_find_bkg_none_when_nothing_qualifies(monkeypatch): + # only candidate has attenuator in -> no file, but a valid period is still returned + lut = _bkg_lut([(3001, "2023-06-10T00:00:00", 3600)]) + fido = _patch_cpd_products(monkeypatch, {3001: np.ones(5)}) + t = Time("2023-06-15T00:00:00") + sel = find_background_file_for_time( + t, + fido_client=fido, + rid_lut=lut, + window_future=7 * u.day, + min_duration=600 * u.s, + require_same_elut=False, + ) + assert sel.path is None + assert sel.rid == -1 + assert sel.valid_from == t + assert sel.valid_to == t + 7 * u.day # no later candidate -> t + window_future + + +# --- same-ELUT criterion ------------------------------------------------------ + + +class _FakeELUTInstance: + def __init__(self, fn): + self._fn = fn + + def _find_elut_file(self, dt): + return self._fn(dt) + + +def _patch_elut(monkeypatch, fn): + import types + + monkeypatch.setattr(flm_mod, "ELUTManager", types.SimpleNamespace(instance=_FakeELUTInstance(fn))) + + +def test_find_bkg_requires_same_elut(monkeypatch): + # 3001 (closer past) is under a DIFFERENT ELUT than the flare; 3002 matches it + lut = _bkg_lut([(3001, "2023-06-10T00:00:00", 3600), (3002, "2023-06-05T00:00:00", 3600)]) + fido = _patch_cpd_products(monkeypatch, {3001: np.zeros(5), 3002: np.zeros(5)}) + # flare(15)->'A', 3001(10)->'B' (different), 3002(5)->'A' (same) + _patch_elut(monkeypatch, lambda dt: "B" if dt.day == 10 else "A") + t = Time("2023-06-15T00:00:00") + + sel = find_background_file_for_time( + t, fido_client=fido, rid_lut=lut, min_duration=600 * u.s, require_same_elut=True + ) + assert sel.rid == 3002 # 3001 skipped: different ELUT configuration + + sel_off = find_background_file_for_time( + t, fido_client=fido, rid_lut=lut, min_duration=600 * u.s, require_same_elut=False + ) + assert sel_off.rid == 3001 # criterion off -> closest past wins + + +def test_find_bkg_same_elut_skipped_when_flare_elut_unknown(monkeypatch): + # ELUT can't be resolved for the flare -> criterion cannot be enforced -> not restrictive + lut = _bkg_lut([(3001, "2023-06-10T00:00:00", 3600)]) + fido = _patch_cpd_products(monkeypatch, {3001: np.zeros(5)}) + _patch_elut(monkeypatch, lambda dt: None) + sel = find_background_file_for_time( + Time("2023-06-15T00:00:00"), fido_client=fido, rid_lut=lut, min_duration=600 * u.s, require_same_elut=True + ) + assert sel.rid == 3001 + + +# --- add_background_file_column (valid-period cache) -------------------------- + + +def test_add_background_file_column_uses_valid_period_cache(monkeypatch): + calls = [] + + def fake_find(time, *, fido_client, **kwargs): + calls.append(time) + # each selection stays valid for 2 days from the queried time + return BackgroundSelection(path=Path("bkg.fits"), rid=42, valid_from=time, valid_to=time + 2 * u.day) + + monkeypatch.setattr(flm_mod, "find_background_file_for_time", fake_find) + + data = QTable() + data["peak_UTC"] = Time( + ["2023-06-15T00:00:00", "2023-06-16T00:00:00", "2023-06-18T00:00:00"] # 3rd is beyond the 1st period + ) + FlareListManager().add_background_file_column(data, fido_client=object()) + + assert len(calls) == 2 # 1st + 3rd flare trigger a search; 2nd reuses the cache + assert list(data["bkg_rid"]) == [42, 42, 42] + assert list(data["bkg_file"]) == ["bkg.fits"] * 3 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..76f5c070 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(2023, 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..642ed876 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,16 +263,17 @@ 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 = flarelist_sdc.find_processing_months(phs) # fl_sdc_months = [] # TODO reactivate once flarelist processing is finalized @@ -287,9 +281,12 @@ def run_daily_pipeline(args): # 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 = fl_to_fl_files[2:-1] # fl_to_fl_files = [] + fl_to_fl_files = [f for f in fl_to_fl_files if "/2024/" in str(f[2])] + # all processing files should be terminated before the next step as the different # processing steeps might create new candidates # let each processing "task" run in its own process @@ -306,16 +303,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 +327,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..e38f5cf2 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 @@ -18,13 +19,15 @@ import astropy.units as u from astropy.coordinates import SkyCoord +from astropy.coordinates.representation import CartesianRepresentation from astropy.io import fits from astropy.table import Column, QTable from astropy.time import Time 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 stx_estimate_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 +42,7 @@ "FlareSOOPMixin", "FlareList", "FlarelistSDCLocImg", - "FlarePeekPreviewMixin", + "FlarePeakPreviewMixin", "FlarelistSC", "FlarelistSCLoc", "FlarelistSCLocImg", @@ -75,7 +78,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 @@ -85,18 +102,19 @@ 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", + location_time_colname="location_time_UTC", keep_all_flares=True, month=None, ): - data["flare_position"] = [SkyCoord(0, 0, frame="icrs", unit="deg") for i in range(0, len(data))] + anc_ephemeris_paths = [] + cpd_paths = [] + position_statuses = [] + position_messages = [] + solo_cartesian_list = [] - 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 @@ -108,12 +126,16 @@ def add_flare_position( day_asp_ephemeris_cache = dict() for i, row in enumerate(data): - if filter_function(row): + _anc_path = "" + _cpd_path = "" + _status = False + _message = "" + 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 - peak_time = row[peek_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: @@ -129,10 +151,28 @@ def add_flare_position( if len(anc_res) < 1: logger.warning(f"No ephemeris data found for flare at time {start_time} : {end_time}") - data[i]["_position_message"] = "no ephemeris data found" + _message = "no ephemeris data found" no_ephemeris += 1 + solo_cartesian_list.append( + ( + np.nan * u.km, + np.nan * u.km, + np.nan * u.km, + peak_time, + 0 * u.s, + np.nan * u.km, + np.nan * u.km, + np.nan * u.km, + 0, + np.nan, + ) + ) + anc_ephemeris_paths.append(_anc_path) + cpd_paths.append(_cpd_path) + position_statuses.append(_status) + position_messages.append(_message) continue - data[i]["anc_ephemeris_path"] = anc_res["path"][0] + _anc_path = str(anc_res["path"][0]) if start_time.datetime.hour < 2: start_time = start_time - 2 * u.hour @@ -145,8 +185,26 @@ def add_flare_position( if len(cpd_res) < 1: logger.warning(f"No CPD data found for flare at time {start_time} : {end_time}") - data[i]["_position_message"] = "no CPD data found" + _message = "no CPD data found" no_cpd += 1 + solo_cartesian_list.append( + ( + np.nan * u.km, + np.nan * u.km, + np.nan * u.km, + peak_time, + 0 * u.s, + np.nan * u.km, + np.nan * u.km, + np.nan * u.km, + 0, + np.nan, + ) + ) + anc_ephemeris_paths.append(_anc_path) + cpd_paths.append(_cpd_path) + position_statuses.append(_status) + position_messages.append(_message) continue if len(cpd_res) > 1: logger.debug(f"Many CPD data found for flare at time {start_time} : {end_time}") @@ -179,22 +237,206 @@ 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: one_cpd += 1 best_cpd_idx = 0 - data[i]["cpd_path"] = cpd_res["path"][best_cpd_idx] + _cpd_path = str(cpd_res["path"][best_cpd_idx]) + + try: + stixpy_cpd = STIXPYProduct(Path(_cpd_path)) + 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 + + _times = stixpy_cpd.data["time"] + _half_bin = stixpy_cpd.data["timedel"] / 2 + mask = (_times + _half_bin >= time_range.start) & (_times - _half_bin <= time_range.end) + data_at_peak = stixpy_cpd.data[mask] + energy_range = [4, 16] * u.keV + + 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 = (_times + _half_bin >= time_range.start) & (_times - _half_bin <= 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}" + ) - # do the calculations with stixpy + rcr_at_peak = data_at_peak["rcr"].max() + if rcr_at_peak > 0: + energy_range = [4, 25] * u.keV - data[i]["flare_position"] = SkyCoord(1, 1, frame="icrs", unit="deg") - data[i]["_position_status"] = True - data[i]["_position_message"] = "OK" + _, flare_loc, sidelobe, solo, img_time_range = stx_estimate_flare_location( + stixpy_cpd, time_range, energy_range + ) + + with SphericalScreen(solo, only_off_disk=True): + center_hgs = flare_loc.transform_to( + HeliographicStonyhurst(obstime=img_time_range.center) + ).cartesian + solo_cartesian_list.append( + ( + center_hgs.x, + center_hgs.y, + center_hgs.z, + img_time_range.center, + img_time_range.seconds, + solo.x, + solo.y, + solo.z, + rcr_at_peak, + sidelobe, + ) + ) + + _status = True + _message = "OK" + except Exception as e: + _status = False + _message = f"Error: {type(e)}" + logger.warning(f"Error calculating flare position for flare at time {start_time} : {end_time}: {e}") + solo_cartesian_list.append( + ( + np.nan * u.km, + np.nan * u.km, + np.nan * u.km, + peak_time, + 0 * u.s, + np.nan * u.km, + np.nan * u.km, + np.nan * u.km, + 0, + np.nan, + ) + ) + anc_ephemeris_paths.append(_anc_path) + cpd_paths.append(_cpd_path) + position_statuses.append(_status) + position_messages.append(_message) else: to_remove.append(i) + solo_cartesian_list.append( + ( + np.nan * u.km, + np.nan * u.km, + np.nan * u.km, + peak_time, + 0 * u.s, + np.nan * u.km, + np.nan * u.km, + np.nan * u.km, + 0, + np.nan, + ) + ) + anc_ephemeris_paths.append(_anc_path) + cpd_paths.append(_cpd_path) + position_statuses.append(False) + position_messages.append("flare did not pass the filter function") + + primer = fido_client.baseurl.replace(fido_client.datapath, "") + primer = primer[7:] if primer.startswith("file://") else primer + + data["anc_ephemeris_path"] = [v.replace(primer, "") for v in anc_ephemeris_paths] + data["anc_ephemeris_path"].info.description = "Path to the daily ancillary ephemeris file" + + data["cpd_path"] = [v.replace(primer, "") for v in cpd_paths] + data["cpd_path"].info.description = "Path to the CPD file used for flare location estimation" + + data["_position_status"] = position_statuses + data["_position_status"].info.description = "Status of the flare position calculation" + + data["_position_message"] = position_messages + data["_position_message"].info.description = "Message describing the status of the flare position calculation" + + flare_x, flare_y, flare_z, solo_times, duration, solo_x, solo_y, solo_z, rcr_at_peak, sidelobe = zip( + *solo_cartesian_list + ) + solo_times = Time(solo_times) + + hgs_coords = SkyCoord( + u.Quantity(flare_x), + u.Quantity(flare_y), + u.Quantity(flare_z), + frame=HeliographicStonyhurst(obstime=solo_times), + representation_type="cartesian", + ) + + solo_coords = SkyCoord( + u.Quantity(solo_x), + u.Quantity(solo_y), + u.Quantity(solo_z), + frame=HeliographicStonyhurst(obstime=solo_times), + representation_type="cartesian", + ) + + # hgc_coords = hgs_coords.transform_to(HeliographicCarrington(obstime=solo_times, observer="Earth")) + hp_coords = hgs_coords.transform_to(Helioprojective(obstime=solo_times, observer="Earth")) + + data["location_hgs"] = hgs_coords + data["location_hgs"].info.description = "Flare location in Heliographic Stonyhurst coordinates" + + data["solo_location_hgs"] = solo_coords + data["solo_location_hgs"].info.description = "SOLO location in Heliographic Stonyhurst coordinates" + + data["sidelobes_ratio"] = sidelobe + data["sidelobes_ratio"].info.description = "Ratio of sidelobes in the STIX image used to assess imaging quality" + + data["rcr_at_peak"] = rcr_at_peak + data[ + "rcr_at_peak" + ].info.description = "max rcr level at flare location estimation time range, > 0 attenuator in place" + + data["visible_from_earth"] = FlarePositionMixin.is_visible(hp_coords) + data[ + "visible_from_earth" + ].info.description = "Whether the flare location is visible from Earth (not occulted by the Sun)" + + data[location_time_colname] = solo_times + data[location_time_colname].info.description = "time center used for flare location estimation in UTC" + + data["location_duration"] = duration + data["location_duration"].info.description = "duration of the flare location estimation time range" + + ( + time_shift, + disc_size, + ) = zip( + *[ + (Spice.instance.get_earth_solo_time_shift(date=scet), Spice.instance.get_sun_disc_size(date=scet)) + for t in solo_times + for scet in (Spice.instance.datetime_to_scet(t),) + ] + ) + + data["time_shift"] = time_shift + data["time_shift"].info.description = "Time(Sun to Earth) - Time(Sun to S/C)" + + data["sun_disc_size"] = disc_size + data["sun_disc_size"].info.description = "Apparent photospheric solar radius" if not keep_all_flares: data.remove_rows(to_remove) @@ -203,23 +445,81 @@ 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) - len(to_remove)} flare locations found" ) + def on_serialize(self, data): + logger.warning( + "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) + if "solo_location_hgs" in data.colnames: + icrs = data["solo_location_hgs"].icrs + icrs_coord = SkyCoord(icrs.ra, icrs.dec, icrs.distance, frame="icrs") + col_idx = data.colnames.index("solo_location_hgs") + data.remove_column("solo_location_hgs") + data.add_column(icrs_coord, name="solo_location_icrs", index=col_idx) + super().on_serialize(data) + + def on_deserialize(self, data, *, location_time_colname=None, **kwargs): + logger.warning( + "FlarePositionMixin on_deserialize called, transforming location columns back to heliographic coordinates" + ) + time_col = location_time_colname or self.location_time_colname + if time_col not in data.colnames: + logger.warning(f"on_deserialize: column '{time_col}' not found, skipping location transform") + else: + obstime = Time(data[time_col]) + if "location_icrs" in data.colnames: + data["location_hgs"] = data["location_icrs"].transform_to(HeliographicStonyhurst(obstime=obstime)) + if "solo_location_icrs" in data.colnames: + data["solo_location_hgs"] = data["solo_location_icrs"].transform_to( + HeliographicStonyhurst(obstime=obstime) + ) + + super().on_deserialize(data, **kwargs) + + @classmethod + def is_visible(cls, coord): + """ + Returns whether the coordinate is on the visible side of the Sun. + This function is a modified version of PR#7118 + """ + + coord = coord.make_3d() + data = coord.cartesian + data_to_sun = coord.observer.radius * CartesianRepresentation(1, 0, 0) - data -class FlareSOOPMixin: + is_behind = data.x < 0 + # print(data.x.to(u.AU)) + is_beyond_limb = np.sqrt(1 - (data.x / data.norm()) ** 2) > coord.rsun / coord.observer.radius + # is_above_surface = data_to_sun.norm() >= coord.rsun + + is_on_near_side = data.dot(data_to_sun) >= 0 + + return is_behind | is_beyond_limb | (is_on_near_side) + + +class FlareSOOPMixin(_SerializeMixin): """_summary_""" @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) @@ -234,14 +534,22 @@ 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 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 +557,7 @@ class FlarePeekPreviewMixin: """ @classmethod - def add_peek_preview( + def add_peak_preview( cls, data, energies, @@ -257,7 +565,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 +574,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 +594,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 +603,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 +696,7 @@ def add_peek_preview( maps.append((map_with_erange, header)) - ppi = PeekPreviewImage( + ppi = PeakPreviewImage( control=QTable(), data=preview_data, month=month, @@ -407,18 +715,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 @@ -480,7 +788,7 @@ class FlarelistSDC(FlareList, FlareSOOPMixin): In L3 product format. """ - PRODUCT_PROCESSING_VERSION = 2 + PRODUCT_PROCESSING_VERSION = 3 NAME = "sdc" def __init__(self, *, service_type=0, service_subtype=0, ssid=2, data, month, **kwargs): @@ -536,7 +844,7 @@ class FlarelistSDCLoc(FlarelistSDC, FlarePositionMixin): In ANC product format. """ - PRODUCT_PROCESSING_VERSION = 2 + PRODUCT_PROCESSING_VERSION = 3 NAME = "sdcloc" def __init__(self, *, service_type=0, service_subtype=0, ssid=3, data, month, **kwargs): @@ -544,6 +852,7 @@ def __init__(self, *, service_type=0, service_subtype=0, ssid=3, data, month, ** self.name = FlarelistSDCLoc.NAME self.ssid = 3 + self.location_time_colname = "location_time_UTC" def enhance_from_product(self, in_prod: GenericProduct): pass @@ -558,10 +867,10 @@ 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, + keep_all_flares=True, month=month, ) @@ -570,7 +879,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 +898,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", @@ -681,6 +990,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 @@ -695,7 +1005,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 +1017,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 +1036,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", @@ -746,3 +1056,44 @@ def add_peek_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/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..73106521 --- /dev/null +++ b/stixcore/products/level3/processing.py @@ -0,0 +1,419 @@ +######################################################### +### 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.sun.constants as sun_const +import sunpy.time +import xrayvision.imaging +from stixpy.calibration.visibility import ( + calibrate_visibility, + create_meta_pixels, + create_visibility, +) +from stixpy.coordinates.frames import STIXImaging +from stixpy.coordinates.transforms import get_hpc_info +from sunpy.coordinates import HeliographicStonyhurst, SphericalScreen, frames +from sunpy.time import TimeRange +from xrayvision.imaging import vis_to_image + +from astropy import units as u +from astropy.coordinates import SkyCoord +from astropy.coordinates.representation import CartesianRepresentation +from astropy.time import Time + + +def get_rsun_obs(observer): + """ + Get the observed radius of the Sun from an observer location. + """ + + rsun_obs = ((sun_const.radius / (observer.spherical.distance - sun_const.radius)).decompose() * u.radian).to( + u.arcsec + ) + return rsun_obs + + +def get_distance_off_limb(coord): + theta_x = coord.Tx + theta_y = coord.Ty + rsun_obs = get_rsun_obs(coord.observer) + distance_off_limb = np.sqrt(theta_x**2 + theta_y**2) - rsun_obs + distance_r_sun = np.sqrt(theta_x**2 + theta_y**2) / rsun_obs + + return distance_off_limb, distance_r_sun + + +def generate_blank_map(date_obs, observer): + """ + Given a date and an observer create a blank map + + """ + data = np.full((12, 12), np.nan) + + # Define a reference coordinate and create a header using sunpy.map.make_fitswcs_header + skycoord = SkyCoord(0 * u.arcsec, 0 * u.arcsec, frame=frames.Helioprojective(observer=observer, obstime=date_obs)) + + # Scale set to the following for solar limb to be in the field of view + header = sunpy.map.make_fitswcs_header(data, skycoord, scale=[600, 600] * u.arcsec / u.pixel) + + # Use sunpy.map.Map to create the blank map + blank_map = sunpy.map.Map(data, header) + return blank_map + + +def is_visible(coord): + """ + Returns whether the coordinate is on the visible side of the Sun. + This function is a modified version of PR#7118 + """ + + coord = coord.make_3d() + data = coord.cartesian + data_to_sun = coord.observer.radius * CartesianRepresentation(1, 0, 0) - data + + is_behind = data.x < 0 + # print(data.x.to(u.AU)) + is_beyond_limb = np.sqrt(1 - (data.x / data.norm()) ** 2) > coord.rsun / coord.observer.radius + # is_above_surface = data_to_sun.norm() >= coord.rsun + + is_on_near_side = data.dot(data_to_sun) >= 0 + + return is_behind | is_beyond_limb | (is_on_near_side) + + +def stx_estimate_flare_location(cpd_sci, time_range, energy_range): + """ + Estimate the flare location using STIX imaging data. + + This function processes the imaging data from the STIX instrument on Solar Orbiter to estimate the location of a solar flare. + It is based on the IDL software `stx_estimate_flare_location`. + + It creates back-projected images in both STIX imaging coordinates and Helioprojective coordinates, and finds the maximum location of the pixel. + + Optionally, it plots the results showing the maximum pixel locations in both coordinate systems. + + Parameters + ---------- + cpd_sci : CPDProduct + the STIX pixel data product. + time_range : `sunpy.time.TimeRange` + The time range over which to estimate the flare location. + energy_range : `astropy.units.Quantity` + The energy range (e.g., in keV) for the analysis. + + Returns + ------- + max_stix : `astropy.coordinates.SkyCoord` + The estimated flare location in STIX imaging coordinates. + max_hpc : `astropy.coordinates.SkyCoord` + The estimated flare location in Helioprojective Cartesian coordinates. + + Notes + ----- + The function involves the following steps: + - Reading STIX pixel data and generating meta pixels for a given time and energy range. + - Creating visibility data from the meta pixels. + - Obtaining solar observer coordinates and converting them to the Heliographic Stonyhurst frame. + - Creating a back-projected image from the visibility data. + - Transforming the coordinates of the maximum pixel in the image to Helioprojective coordinates. + + """ + + meta_pixels_sci = create_meta_pixels( + cpd_sci, + time_range=[time_range.start, time_range.end], + energy_range=energy_range, + flare_location=[0, 0] * u.arcsec, + no_shadowing=True, + ) + + # create visibilities + vis = create_visibility(meta_pixels_sci) + vis_tr = TimeRange(vis.meta["time_range"]) + + roll, solo_xyz, pointing = get_hpc_info(vis_tr.start, vis_tr.end) + solo = frames.HeliographicStonyhurst(*solo_xyz, obstime=vis_tr.center, representation_type="cartesian") + + center_map = SkyCoord(0 * u.arcsec, 0 * u.arcsec, frame=frames.Helioprojective(observer=solo, obstime=solo.obstime)) + center_coord = center_map.transform_to(STIXImaging(obstime=vis_tr.start, obstime_end=vis_tr.end, observer=solo)) + + # get calibrated visibilities - use center of Sun as phase center + cal_vis = calibrate_visibility(vis, flare_location=center_coord) + + # order by sub-collimator e.g. 10a, 10b, 10c, 9a, 9b, 9c .... + isc_10_7 = [3, 20, 22, 16, 14, 32, 21, 26, 4, 24, 8, 28] + idx = np.argwhere(np.isin(cal_vis.meta["isc"], isc_10_7)).ravel() + + # only use subcolimators 7 - 10 + vis10_7 = cal_vis[idx] + + # set up image size + imsize = [512, 512] * u.pixel + + # to make sure the full Sun is within FOV - the 2.6 is taken to be the same as the IDL software + pixel = get_rsun_obs(solo) * 2.6 / imsize + + # get back projection image + bp_image = vis_to_image(vis10_7, imsize, pixel_size=pixel) + + # Make a sunpy map from the bp_image, in STIX imaging frame + header = sunpy.map.make_fitswcs_header( + bp_image, center_coord, telescope="STIX", observatory="Solar Orbiter", scale=pixel + ) + fd_bp_map = sunpy.map.Map((bp_image, header)) + + sidelobes_ratio = calculate_sidelobes_ratio(fd_bp_map) + + # Make a sunpy map from the bp_image, in HPC from STIX observer + hpc_ref = center_coord.transform_to(frames.Helioprojective(observer=solo, obstime=vis_tr.center)) + header_hp = sunpy.map.make_fitswcs_header(bp_image, hpc_ref, scale=pixel, rotation_angle=90 * u.deg + roll) + hp_map = sunpy.map.Map((bp_image, header_hp)) + + # get the position of the max pixel + max_pixel = np.argwhere(fd_bp_map.data == fd_bp_map.data.max()).ravel() * u.pixel + # get the world coord of the max pixel - (note WCS axes and array are reversed) + max_stix = fd_bp_map.pixel_to_world(max_pixel[1], max_pixel[0]) + + # get the coordinate of the max pixel in HPC - if coordinate is off limb, assume spherical screen for transform + with SphericalScreen(hp_map.observer_coordinate, only_off_disk=True): + max_hpc = max_stix.transform_to(hp_map.coordinate_frame) + + vis_time_range = TimeRange(vis.meta["time_range"][0], vis.meta["time_range"][1]) + return max_stix, max_hpc, sidelobes_ratio, solo, vis_time_range + + +def calculate_sidelobes_ratio(bp_nat_map, threshold=200 * u.arcsec): + """ + + Calculate the sidelobes ratio for a back-projected image map. + + The sidelobes ratio is a measure of the relative strength of the sidelobes compared to the main peak of the image. + This ratio helps determine the reliability of the flare location. A sidelobes ratio close to or above 0.9 suggests + that the flare location may not be reliable due to significant sidelobe interference. + + Parameters + ---------- + bp_nat_map : `sunpy.map.Map` + The back-projected image map (in natural units) to analyze. This map is typically generated from visibility data + and contains the image of the flare. + threshold : `astropy.units.Quantity`, optional + The angular separation threshold (in arcseconds) around the peak within which sidelobes are excluded from the calculation. + Default is 200 arcseconds. + + Returns + ------- + sidelobes_ratio : float + The ratio of the maximum sidelobe intensity to the peak intensity in the back-projected image. + A value close to 1 indicates significant sidelobes, potentially making the flare location unreliable. + + Notes + ----- + - This is based upon the methodology in the STIX-GSW IDL software. + """ + max_bp = np.max(bp_nat_map.data) + ind_max = np.unravel_index(np.argmax(bp_nat_map.data, axis=None), bp_nat_map.data.shape) + max_bp_coord = bp_nat_map.pixel_to_world(ind_max[1] * u.pix, ind_max[0] * u.pix) + + yy, xx = np.indices(bp_nat_map.data.shape) + world_coords = bp_nat_map.pixel_to_world(xx * u.pix, yy * u.pix) + + distance_wrt_peak = world_coords.separation(max_bp_coord) + + bp_image_masked = np.copy(bp_nat_map.data) + mask = distance_wrt_peak <= threshold + bp_image_masked[mask] = 0 + + sidelobes_ratio = np.max(bp_image_masked) / max_bp + + return sidelobes_ratio + + +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, solo + + +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..4d9dbf7d 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 + 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,7 +94,7 @@ 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: @@ -336,6 +339,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 +615,23 @@ 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(). + + 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().""" + deserialize = getattr(super(), "on_deserialize", None) + if deserialize is not None: + deserialize(data, **kwargs) + 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..07522646 --- /dev/null +++ b/stixcore/products/tests/test_flarelist.py @@ -0,0 +1,191 @@ +from datetime import date + +import numpy as np +import pytest +from sunpy.coordinates import HeliographicStonyhurst +from sunpy.time import TimeRange + +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 +from stixcore.products.level3.flarelist import ( + FlarelistSDCLoc, + calculate_overlap, + longest_constant_sequence, +) +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 + + 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=lon * u.deg, + lat=lat * u.deg, + radius=np.ones(N) * 1.0 * u.AU, + frame=HeliographicStonyhurst(obstime=peak_times), + ) + + data = QTable() + data["peak_UTC"] = peak_times + data["location_time_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 + + return data + + +@pytest.fixture +def written_fits(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 + 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() + + # read back via Product factory — calls on_deserialize internally + recovered = Product(fits_path) + + assert isinstance(recovered, FlarelistSDCLoc) + 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["location_time_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) + + +# --- longest_constant_sequence --- + + +def test_lcs_empty(): + length, start, state = longest_constant_sequence([]) + assert length == 0 + assert start is None + assert state is None + + +def test_lcs_single_element(): + length, start, state = longest_constant_sequence([5]) + assert length == 1 + assert start == 0 + assert state == 5 + + +def test_lcs_all_same(): + length, start, state = longest_constant_sequence([3, 3, 3, 3]) + assert length == 4 + assert start == 0 + assert state == 3 + + +def test_lcs_clear_winner(): + length, start, state = longest_constant_sequence([1, 2, 2, 2, 3, 3]) + assert length == 3 + assert start == 1 + assert state == 2 + + +def test_lcs_tie_prefers_lower_state(): + # two runs of length 2: state=1 at index 0, state=3 at index 2 + length, start, state = longest_constant_sequence([1, 1, 3, 3]) + assert length == 2 + assert start == 0 + assert state == 1 + + +def test_lcs_numpy_array(): + arr = np.array([0, 0, 1, 1, 1, 0]) + length, start, state = longest_constant_sequence(arr) + assert length == 3 + assert start == 2 + assert state == 1 + + +# --- calculate_overlap --- + + +def test_overlap_no_intersection(): + r1 = TimeRange("2024-01-01T00:00:00", "2024-01-01T01:00:00") + r2 = TimeRange("2024-01-01T02:00:00", "2024-01-01T03:00:00") + assert calculate_overlap(r1, r2) is None + + +def test_overlap_partial(): + r1 = TimeRange("2024-01-01T00:00:00", "2024-01-01T02:00:00") + r2 = TimeRange("2024-01-01T01:00:00", "2024-01-01T03:00:00") + result = calculate_overlap(r1, r2) + assert result is not None + assert result.start == Time("2024-01-01T01:00:00") + assert result.end == Time("2024-01-01T02:00:00") + + +def test_overlap_contained(): + r1 = TimeRange("2024-01-01T00:00:00", "2024-01-01T04:00:00") + r2 = TimeRange("2024-01-01T01:00:00", "2024-01-01T03:00:00") + result = calculate_overlap(r1, r2) + assert result is not None + assert result.start == Time("2024-01-01T01:00:00") + assert result.end == Time("2024-01-01T03:00:00") + + +def test_overlap_identical(): + r1 = TimeRange("2024-01-01T00:00:00", "2024-01-01T01:00:00") + result = calculate_overlap(r1, r1) + assert result is not None + assert result.start == r1.start + assert result.end == r1.end