From b8858a93d9ab830f46a9ad12f24929516f235e33 Mon Sep 17 00:00:00 2001 From: Nicky Hochmuth Date: Fri, 27 Feb 2026 16:47:47 +0100 Subject: [PATCH] skip QLSpectraReshapeError in LB to L0 step --- processing_history.sqlite | Bin 16384 -> 0 bytes stixcore/processing/LBtoL0.py | 5 + stixcore/util/scripts/COMPONENT_STORAGE.md | 211 ---------------- .../scripts/SKYCOORD_FITS_SERIALIZATION.md | 225 ------------------ stixcore/util/scripts/SUMMARY.md | 94 -------- stixcore/util/scripts/skycoord.py | 178 -------------- stixcore/util/scripts/timecolumn.py | 68 ------ 7 files changed, 5 insertions(+), 776 deletions(-) delete mode 100644 processing_history.sqlite delete mode 100644 stixcore/util/scripts/COMPONENT_STORAGE.md delete mode 100644 stixcore/util/scripts/SKYCOORD_FITS_SERIALIZATION.md delete mode 100644 stixcore/util/scripts/SUMMARY.md delete mode 100644 stixcore/util/scripts/skycoord.py delete mode 100644 stixcore/util/scripts/timecolumn.py diff --git a/processing_history.sqlite b/processing_history.sqlite deleted file mode 100644 index f28f98f7fe3d15eabaf74fc25cfdd1e8838f9060..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 16384 zcmeI2PjBNy6u|AIZJM-67tm4`4!F0iwn=6@{_hK+1XReTOV-P7B_t~}E<{R`Z0)X= z0|yr&#Dyc$Ubyoa_y8QZaO8r7#0P;lPN|wmb~2?%Tj`DJwv)WEfAgDHPp!O%AGi8Z zpd5`(j=PbTgqt`9)EY(Jq~u1LGUaX>?qNfW9U2!PQw0Z zxT7R9>ks{LH~M1xL2hO9-aYvs>V7r|!lwgxLh{4l>1i+o57XWJ`AogWyRD`&^XBA= zFmn5w*X6>d>&oBlCg)6wdHC?JRsYZdLJ`Pduf)567ccYw@wkT3mH7U%{)6Jzw_=58~*T6XvQE$vH-+ zQId66J;~Vbb>RtqY6TWvylwt$zP9qGRQaX+{rYd~udaEeABw*h-(3B)@FTp84Go|H zG=K)sz>PEzyK5!swO8Jjx6X=^dC3cn2I&7icdKmSe*I_%Z zM?o0X$Ba_{aCAKOlWvREd3>GNb%=kz)#x-^dk?*OI2w%n7W4Vw>czRMa4Db16Pr+v zB(2K(luY`DYWbfKwZ{JiNI0kpC>}K&YCBA;nU-prR_vDM&|8A^mKeP?O;0RFJ(V#+ z?KFDJp?b}-n688bRlqXRu5Fvr+kokMR3(i>v*#tgl| z2zPJlEX}>)Fj_NUh8AQ_VGcc{zEIvn>Wk#PK1-kUDWR$KI;U5errCDvuFRp=1nD&~ zdToYYOWdmlGt%5^j%wCS-D0*8yUTOvnIJtAqldrG%l8~InGIDnHI1G*M6W@ZWkcRC z&7p@ol~CrWV)W_^y{cM}_c}|@9O`IVO*c%z26gu z(weQ=jF@LeZcVxayVY6oY=duCgVEHEC9+RlzxDs+>M+B}Mk%uFmPakvB^@1ZceZIz z5GP-oxa(hV5pfG&{X(2b#!bW#hp9EgGEH0M&;RnqUlM#_Lj!044WI!ufCkV28bAYR z01co4G=K(fu7S7Yyep?vjlhQ^@cDoFg;e>e@@*xk>{hnQf0zF#|6G2taa8&B=1v;P z(Eu7i184vZpaC?12G9T+Km#|%z*ZrtT9Z?2*48p3cwJa2GlJJ%6*D7vebZ`Y1g{|~ zWJd5hot4Z8UJJ9F8NusKmNFxFO-DYyvpOlA;LcvE!iik`#yMK?zo9K$955-VSfDcC EAE2aCFaQ7m diff --git a/stixcore/processing/LBtoL0.py b/stixcore/processing/LBtoL0.py index 251ef7b5..75a99d44 100644 --- a/stixcore/processing/LBtoL0.py +++ b/stixcore/processing/LBtoL0.py @@ -9,6 +9,7 @@ from stixcore.idb.manager import IDBManager from stixcore.io.product_processors.fits.processors import FitsL0Processor from stixcore.io.RidLutManager import RidLutManager +from stixcore.products.level0.quicklookL0 import QLSpectraReshapeError from stixcore.products.level0.scienceL0 import NotCombineException from stixcore.products.product import Product from stixcore.util.logging import get_logger @@ -103,6 +104,10 @@ def process_tm_type(files, tm_type, processor, spice_kernel_path, config, idbm): if level0: fits_files = processor.write_fits(level0) all_files.extend(fits_files) + except QLSpectraReshapeError as e: + logger.warning(str(e)) + if CONFIG.getboolean("Logging", "stop_on_error", fallback=False): + raise e except Exception as e: logger.error( f"Error processing file {file} for {levelb.service_type}, {levelb.service_subtype}, {levelb.ssid}", diff --git a/stixcore/util/scripts/COMPONENT_STORAGE.md b/stixcore/util/scripts/COMPONENT_STORAGE.md deleted file mode 100644 index 61764e5c..00000000 --- a/stixcore/util/scripts/COMPONENT_STORAGE.md +++ /dev/null @@ -1,211 +0,0 @@ -# Storing Coordinate Components Instead of SkyCoord Objects - -## Overview - -Instead of storing SkyCoord objects in FITS files, store the individual coordinate components (Tx, Ty for Helioprojective, or lon, lat for HeliographicStonyhurst) along with observer position and time. This approach is simpler, more transparent, and avoids many serialization issues. - -## Benefits - -1. **Simplicity**: No need for `fits.connect._encode_mixins()` or ICRS conversion -2. **Transparency**: FITS files contain plain columns with clear names and units -3. **No Distance Issues**: No need to handle distance/unit sphere complications -4. **Direct Reconstruction**: Create SkyCoord directly from stored components -5. **Frame Agnostic**: Works for any coordinate frame -6. **Better Inspection**: Easy to view coordinates with standard FITS tools - -## What to Store - -For **Helioprojective** coordinates: -- `hp_tx`: Tx coordinate (arcsec) -- `hp_ty`: Ty coordinate (arcsec) -- `observer_hgs_x`: Observer x position in HeliographicStonyhurst (km) -- `observer_hgs_y`: Observer y position in HeliographicStonyhurst (km) -- `observer_hgs_z`: Observer z position in HeliographicStonyhurst (km) -- `peak_UTC` (or relevant time): Observation time - -For **HeliographicStonyhurst** coordinates: -- `hgs_lon`: Longitude (deg) -- `hgs_lat`: Latitude (deg) -- `hgs_radius`: Radius (km) - optional, can default to solar radius -- `peak_UTC`: Observation time - -## Example: Writing to FITS - -```python -from astropy.table import QTable -from astropy.io import fits -from astropy.io.fits import table_to_hdu -import astropy.units as u -from sunpy.coordinates import Helioprojective, HeliographicStonyhurst -from stixpy.coordinates.transforms import get_hpc_info - -# Get coordinates and observer position -coord_hp = ... # Your Helioprojective coordinate -roll, solo_xyz, pointing = get_hpc_info(start_utc, end_utc) - -# Extract components -hp_tx = coord_hp.Tx -hp_ty = coord_hp.Ty -obs_x = solo_xyz[0] -obs_y = solo_xyz[1] -obs_z = solo_xyz[2] - -# Create table with component columns -data = QTable() -data["peak_UTC"] = [peak_utc] -data["hp_tx"] = [hp_tx] -data["hp_ty"] = [hp_ty] -data["observer_hgs_x"] = [obs_x] -data["observer_hgs_y"] = [obs_y] -data["observer_hgs_z"] = [obs_z] - -# Write to FITS - no special encoding needed! -data_enc = fits.connect._encode_mixins(data) # Still handles Time, Quantity -data_hdu = table_to_hdu(data_enc) -hdul = fits.HDUList([fits.PrimaryHDU(), data_hdu]) -hdul.writeto('coordinates.fits', overwrite=True) -``` - -## Example: Reading from FITS - -```python -from astropy.table import QTable -from astropy.coordinates import SkyCoord -from sunpy.coordinates import Helioprojective, HeliographicStonyhurst - -# Read back -data = QTable.read('coordinates.fits', astropy_native=True) - -# Extract components (they have units!) -hp_tx = data['hp_tx'][0] -hp_ty = data['hp_ty'][0] -obs_x = data['observer_hgs_x'][0] -obs_y = data['observer_hgs_y'][0] -obs_z = data['observer_hgs_z'][0] -obstime = data['peak_UTC'][0] - -# Reconstruct observer (NO SPICE!) -observer = HeliographicStonyhurst( - x=obs_x, y=obs_y, z=obs_z, - obstime=obstime, - representation_type='cartesian' -) - -# Reconstruct coordinate -coord_hp = SkyCoord( - hp_tx, hp_ty, - frame=Helioprojective(obstime=obstime, observer=observer) -) - -# Done! coord_hp is now a proper Helioprojective coordinate -``` - -## Comparison with SkyCoord Storage - -### Old Approach (storing SkyCoord objects): -1. Convert coordinate to ICRS -2. Drop distance to make frames compatible -3. Use `concatenate()` to create arrays -4. Encode with `fits.connect._encode_mixins()` -5. Write to FITS -6. Read back -7. Add distance back -8. Convert to desired frame - -**Problems:** -- Complex serialization -- Distance handling issues -- Frame conversion complications -- Not transparent in FITS files - -### New Approach (storing components): -1. Extract Tx, Ty (or lon, lat) -2. Extract observer x, y, z -3. Store in table columns -4. Write to FITS -5. Read back -6. Create SkyCoord from components - -**Advantages:** -- Simple and direct -- No frame conversions needed -- No distance issues -- Transparent FITS structure -- Fewer steps - -## Multiple Coordinate Frames - -You can store multiple coordinate representations: - -```python -# Store both Helioprojective and HeliographicStonyhurst -data["hp_tx"] = [coord_hp.Tx] -data["hp_ty"] = [coord_hp.Ty] -data["hgs_lon"] = [coord_hgs.lon] -data["hgs_lat"] = [coord_hgs.lat] -data["hgs_radius"] = [coord_hgs.radius] -data["observer_hgs_x"] = [obs_x] -data["observer_hgs_y"] = [obs_y] -data["observer_hgs_z"] = [obs_z] -``` - -Each can be reconstructed independently: - -```python -# Reconstruct Helioprojective -coord_hp = SkyCoord(hp_tx, hp_ty, frame=Helioprojective(obstime, observer)) - -# Reconstruct HeliographicStonyhurst -coord_hgs = SkyCoord(hgs_lon, hgs_lat, hgs_radius, - frame=HeliographicStonyhurst(obstime=obstime)) -``` - -## Column Naming Conventions - -Suggested naming: -- **Helioprojective**: `hp_tx`, `hp_ty` -- **HeliographicStonyhurst**: `hgs_lon`, `hgs_lat`, `hgs_radius` -- **HeliographicCarrington**: `hgc_lon`, `hgc_lat`, `hgc_radius` -- **Observer position**: `observer_hgs_x`, `observer_hgs_y`, `observer_hgs_z` -- **Time**: `peak_UTC`, `start_UTC`, `end_UTC` - -## Units - -All units are automatically preserved by `astropy.table.QTable`: -- Angles: `astropy.units.arcsec`, `astropy.units.deg` -- Distances: `astropy.units.km`, `astropy.units.AU` -- Times: `astropy.time.Time` objects - -When reading with `astropy_native=True`, units are automatically restored. - -## Integration with Existing STIXCore Code - -This approach can coexist with existing code: - -```python -# If you have existing SkyCoord columns, convert to components -if 'flare_location' in data.colnames: - # Extract components - data['hp_tx'] = [coord.Tx for coord in data['flare_location']] - data['hp_ty'] = [coord.Ty for coord in data['flare_location']] - # Remove old column - data.remove_column('flare_location') - -# Or keep both for backwards compatibility -data['flare_location'] = coords_array # Old format -data['hp_tx'] = [coord.Tx for coord in coords_array] # New format -data['hp_ty'] = [coord.Ty for coord in coords_array] -``` - -## Summary - -**Key Points:** -1. Store coordinate components (Tx, Ty, lon, lat) as separate columns -2. Store observer position as separate x, y, z columns -3. Store observation time -4. No SkyCoord serialization needed - just plain columns with units -5. Reconstruction is simple: `SkyCoord(Tx, Ty, frame=...)` -6. No SPICE needed when reading - all data is in the FITS file -7. FITS files are transparent and easy to inspect - -This approach provides the simplest path to storing and retrieving astronomical coordinates in FITS files while maintaining full scientific accuracy. diff --git a/stixcore/util/scripts/SKYCOORD_FITS_SERIALIZATION.md b/stixcore/util/scripts/SKYCOORD_FITS_SERIALIZATION.md deleted file mode 100644 index bdd9d414..00000000 --- a/stixcore/util/scripts/SKYCOORD_FITS_SERIALIZATION.md +++ /dev/null @@ -1,225 +0,0 @@ -# SkyCoord FITS Serialization Guide - -## Problem - -When trying to serialize Astropy `SkyCoord` objects to FITS files, you may encounter two types of errors: - -1. **Single coordinate with complex frame**: - ``` - ('cannot represent an object', ) - ``` - -2. **Multiple coordinates (mixed types)**: - ``` - Column 'flare_location_stix' contains unsupported object types or mixed types: {dtype('O')} - ``` - -## Root Cause - -The issue occurs because: -1. FITS format doesn't natively support complex Python objects like `SkyCoord` -2. Astropy's `fits.connect._encode_mixins()` can serialize **simple** coordinate frames (ICRS, FK5, Galactic) -3. **Complex frames with observer information** (HeliographicStonyhurst, Helioprojective, STIXImaging) contain additional metadata that cannot be automatically serialized - -## Solution 1: Convert to ICRS Frame (Recommended) - -Convert your coordinates to the ICRS frame before storing. This is the approach used in STIXCore's flarelist. - -```python -from astropy.io import fits -from astropy.io.fits import table_to_hdu -from astropy.table import QTable -from astropy.coordinates import SkyCoord - -# Your original coordinates in complex frames -coord_hp = SkyCoord(0 * u.deg, 0 * u.deg, - frame=Helioprojective(obstime=peak_utc, observer=solo)) -coord_stix = coord_hp.transform_to(STIXImaging(obstime=start_utc, - obstime_end=end_utc, - observer=solo)) - -# Convert to ICRS for serialization -coord_stix_icrs = coord_stix.transform_to('icrs') -coord_hp_icrs = coord_hp.transform_to('icrs') - -# Create table with ICRS coordinates -data = QTable() -data["flare_location_stix"] = [coord_stix_icrs] -data["flare_location_hp"] = [coord_hp_icrs] - -# Serialize using fits.connect._encode_mixins -data_enc = fits.connect._encode_mixins(data) -data_hdu = table_to_hdu(data_enc) -data_hdu.name = "DATA" - -# Write to FITS -primary_hdu = fits.PrimaryHDU() -hdul = fits.HDUList([primary_hdu, data_hdu]) -hdul.writeto('output.fits', overwrite=True) - -# Read back with astropy_native=True to reconstruct SkyCoord -tr = QTable.read('output.fits', astropy_native=True) -# tr['flare_location_stix'] is now a SkyCoord in ICRS frame -``` - -### Why This Works - -- ICRS is a simple, time-independent reference frame -- No observer position or observation time metadata to serialize -- `fits.connect._encode_mixins()` knows how to handle ICRS coordinates -- The coordinates can be fully reconstructed when reading back with `astropy_native=True` - -## Solution 2: Store Components Separately - -For maximum control, store RA/Dec as separate columns: - -```python -# Store coordinate components -data = QTable() -data["stix_ra"] = [coord_stix_icrs.ra] -data["stix_dec"] = [coord_stix_icrs.dec] - -# Write as before -data_enc = fits.connect._encode_mixins(data) -data_hdu = table_to_hdu(data_enc) -data_hdu.name = "DATA" - -primary_hdu = fits.PrimaryHDU() -hdul = fits.HDUList([primary_hdu, data_hdu]) -hdul.writeto('output.fits', overwrite=True) - -# Read back and reconstruct -tr = QTable.read('output.fits', astropy_native=True) -coord_reconstructed = SkyCoord(ra=tr['stix_ra'][0], - dec=tr['stix_dec'][0], - frame='icrs') -``` - -### Pros and Cons - -**Solution 1 (ICRS conversion)**: -- ✅ Preserves full SkyCoord structure -- ✅ Automatic reconstruction when reading -- ✅ Used in STIXCore flarelist code -- ⚠️ Loses original frame information (must reconvert if needed) - -**Solution 2 (Separate components)**: -- ✅ Maximum control over serialization -- ✅ Simpler FITS structure -- ✅ Easy to understand and debug -- ⚠️ Manual reconstruction required -- ⚠️ Need separate columns for each coordinate - -## Key Points - -1. **Always use `fits.connect._encode_mixins()`** before writing QTables with complex types to FITS -2. **Convert to simple frames** (ICRS, FK5, Galactic) before serialization -3. **Use `astropy_native=True`** when reading to reconstruct SkyCoord objects -4. **Store observer/frame metadata separately** if you need to reconstruct the original frame - -## Example in STIXCore - -The flarelist code always uses ICRS frame for stored coordinates: - -```python -# From stixcore/products/level3/flarelist.py:94 -data["flare_position"] = [SkyCoord(0, 0, frame="icrs", unit="deg") - for i in range(0, len(data))] -``` - -## Converting Back from ICRS (Without SPICE!) - -Yes, you can convert ICRS coordinates back to HeliographicStonyhurst or Helioprojective! The key is to **store the observer position in the FITS file** so you don't need SPICE for reconstruction. - -### Storing Data (Write Time) - -```python -from sunpy.coordinates import HeliographicStonyhurst, Helioprojective -from astropy.coordinates import SkyCoord - -# Get observer position from SPICE (only needed once at write time) -roll, solo_xyz, pointing = get_hpc_info(start_utc, end_utc) -# solo_xyz is a Quantity vector [x, y, z] in km - -# Convert coordinates to ICRS for storage -coord_icrs = coord_stix.transform_to('icrs') - -# For multiple coordinates, drop distance and create a SkyCoord array using concatenate (not a list) -# Distance must be dropped to make frames equivalent for concatenation -from astropy.coordinates import concatenate -coord_icrs_1_nodist = SkyCoord(ra=coord_icrs_1.ra, dec=coord_icrs_1.dec, frame='icrs') -coord_icrs_2_nodist = SkyCoord(ra=coord_icrs_2.ra, dec=coord_icrs_2.dec, frame='icrs') -coords_array = concatenate([coord_icrs_1_nodist, coord_icrs_2_nodist]) # if you have multiple - -# Store everything needed for reconstruction -data = QTable() -data["peak_UTC"] = [peak_utc] -data["flare_location_icrs"] = coords_array # or [coord_icrs] for single coord -# Store observer position as Cartesian vector (no SPICE needed later!) -data["observer_hgs_xyz"] = [solo_xyz] # [x, y, z] in km - -# Write to FITS -data_enc = fits.connect._encode_mixins(data) -data_hdu = table_to_hdu(data_enc) -# ... write as shown above -``` - -### Reading and Converting Back (Read Time - NO SPICE!) - -```python -from sunpy.coordinates import HeliographicStonyhurst, Helioprojective - -# Read back from FITS -tr = QTable.read('output.fits', astropy_native=True) - -# Extract coordinate and time -coord_icrs = tr['flare_location_icrs'][0] -obstime = tr['peak_UTC'][0] -xyz = tr['observer_hgs_xyz'][0] # [x, y, z] vector with units - -# Reconstruct observer from stored Cartesian vector -# NO SPICE REQUIRED! -observer = HeliographicStonyhurst( - x=xyz[0], y=xyz[1], z=xyz[2], - obstime=obstime, - representation_type='cartesian' -) - -# Convert back to Helioprojective -# Add distance to the coordinate (stored coords are on unit sphere without distance) -coord_icrs_with_dist = SkyCoord( - ra=coord_icrs.ra, dec=coord_icrs.dec, - distance=observer.radius, # Use Sun's distance from observer - frame='icrs' -) -coord_hp = coord_icrs_with_dist.transform_to(Helioprojective(obstime=obstime, observer=observer)) -``` - -### Important Notes - -- **Store observer position in FITS** - This eliminates the need for SPICE when reading -- **Store observation times** - ICRS has no time dependence, store `peak_UTC`, `start_UTC`, etc. -- **Distance information may vary** - ICRS preserves sky position; radial distance can change slightly -- **Self-contained FITS files** - Everything needed for coordinate transformations is in the file - -### Best Practice: Complete Storage - -Store these columns for full reconstruction without external dependencies: - -```python -data = QTable() -data["peak_UTC"] = [obstime] # For time-dependent transformations -data["start_UTC"] = [start_time] # Optional: integration window -data["end_UTC"] = [end_time] # Optional: integration window -data["flare_location_icrs"] = [coord_icrs] # Sky position -# Observer position as Cartesian vector (eliminates SPICE dependency) -data["observer_hgs_xyz"] = [solo_xyz] # [x, y, z] in km -# Optional: store distance if needed -data["flare_distance"] = [coord.distance] # Optional: 3D position -``` - -This approach makes your FITS files **self-contained** - no SPICE kernels needed for reading and transforming coordinates! - -## Working Example - -See [skycoord.py](stixcore/util/scripts/skycoord.py) for a complete working example that demonstrates both methods and round-trip conversion. diff --git a/stixcore/util/scripts/SUMMARY.md b/stixcore/util/scripts/SUMMARY.md deleted file mode 100644 index 110a71fd..00000000 --- a/stixcore/util/scripts/SUMMARY.md +++ /dev/null @@ -1,94 +0,0 @@ -# SkyCoord FITS Serialization - Quick Reference - -## The Problem -Complex coordinate frames (HeliographicStonyhurst, Helioprojective, STIXImaging) cannot be directly serialized to FITS because they contain observer metadata. - -## The Solution -Convert to ICRS frame and store observer position separately. - -## Quick Code Template - -### Writing FITS (with observer data) - -```python -from astropy.io import fits -from astropy.io.fits import table_to_hdu -from astropy.table import QTable -from astropy.coordinates import SkyCoord - -# 1. Get observer position (from SPICE, only needed at write time) -roll, solo_xyz, pointing = get_hpc_info(start_utc, end_utc) -# solo_xyz is a Quantity vector [x, y, z] in km - -# 2. Convert coordinates to ICRS -coord_icrs = coord_stix.transform_to('icrs') - -# For multiple coordinates, drop distance and create SkyCoord array using concatenate -# Distance must be dropped to make frames equivalent for concatenation -from astropy.coordinates import concatenate -coord_icrs_1_nodist = SkyCoord(ra=coord_icrs_1.ra, dec=coord_icrs_1.dec, frame='icrs') -coord_icrs_2_nodist = SkyCoord(ra=coord_icrs_2.ra, dec=coord_icrs_2.dec, frame='icrs') -coords_array = concatenate([coord_icrs_1_nodist, coord_icrs_2_nodist]) # or [coord_icrs] for single - -# 3. Create table with all necessary data -data = QTable() -data["peak_UTC"] = [peak_utc] -data["flare_location"] = coords_array # SkyCoord array, not list -# Store observer as Cartesian vector (no SPICE needed later!) -data["observer_hgs_xyz"] = [solo_xyz] # [x, y, z] in km - -# 4. Write to FITS -data_enc = fits.connect._encode_mixins(data) -data_hdu = table_to_hdu(data_enc) -data_hdu.name = "DATA" - -primary_hdu = fits.PrimaryHDU() -hdul = fits.HDUList([primary_hdu, data_hdu]) -hdul.writeto('output.fits', overwrite=True) -``` - -### Reading FITS and Converting Back (NO SPICE!) - -```python -from sunpy.coordinates import HeliographicStonyhurst, Helioprojective -from astropy.table import QTable - -# 1. Read from FITS -data = QTable.read('output.fits', astropy_native=True) - -# 2. Extract data -coord_icrs = data['flare_location'][0] -obstime = data['peak_UTC'][0] -xyz = data['observer_hgs_xyz'][0] # [x, y, z] vector with units - -# 3. Reconstruct observer from Cartesian vector (NO SPICE!) -observer = HeliographicStonyhurst( - x=xyz[0], y=xyz[1], z=xyz[2], - obstime=obstime, - representation_type='cartesian' -) - -# 4. Convert to desired frame -# Add distance to the coordinate (stored coords are on unit sphere without distance) -coord_icrs_with_dist = SkyCoord( - ra=coord_icrs.ra, dec=coord_icrs.dec, - distance=observer.radius, # Use Sun's distance from observer - frame='icrs' -) -coord_hp = coord_icrs_with_dist.transform_to( - Helioprojective(obstime=obstime, observer=observer) -) -``` - -## Key Benefits - -✅ **No SPICE dependency for reading** - All observer data stored in FITS -✅ **Self-contained files** - Everything needed is in the FITS file -✅ **Standard approach** - Used in STIXCore flarelist code -✅ **Preserves sky position** - RA/Dec accurately maintained -✅ **Enables coordinate transformations** - Can convert to any frame - -## Files - -- **[SKYCOORD_FITS_SERIALIZATION.md](SKYCOORD_FITS_SERIALIZATION.md)** - Complete documentation -- **[skycoord.py](skycoord.py)** - Working example with tests diff --git a/stixcore/util/scripts/skycoord.py b/stixcore/util/scripts/skycoord.py deleted file mode 100644 index e341558b..00000000 --- a/stixcore/util/scripts/skycoord.py +++ /dev/null @@ -1,178 +0,0 @@ -from pathlib import Path - -import numpy as np -from stixpy.coordinates.transforms import STIXImaging, get_hpc_info -from sunpy.coordinates import HeliographicStonyhurst, Helioprojective - -import astropy.units as u -from astropy.coordinates import SkyCoord -from astropy.io import fits -from astropy.io.fits import table_to_hdu -from astropy.table.table import QTable -from astropy.time import Time - -from stixcore.config.config import CONFIG -from stixcore.ephemeris.manager import Spice, SpiceKernelManager - -if __name__ == "__main__": - _spm = SpiceKernelManager(Path(CONFIG.get("Paths", "spice_kernels"))) - spicemeta = _spm.get_latest_mk_and_pred() - - Spice.instance = Spice(spicemeta) - - now = Time.now() - - start_utc = now - 100 * u.d - peak_utc = now - 98 * u.d - end_utc = now - 96 * u.d - - start_utc_2 = now - 10 * u.d - peak_utc_2 = now - 9 * u.d - end_utc_2 = now - 9 * u.d - - roll_1, solo_xyz_1, pointing_1 = get_hpc_info(start_utc, end_utc) - solo_1 = HeliographicStonyhurst(*solo_xyz_1, obstime=peak_utc, representation_type="cartesian") - center_hpc = SkyCoord(0 * u.deg, 0 * u.deg, frame=Helioprojective(obstime=peak_utc, observer=solo_1)) - coord_stix = center_hpc.transform_to(STIXImaging(obstime=start_utc, obstime_end=end_utc, observer=solo_1)) - - coord_hp = coord_stix.transform_to(Helioprojective(observer=solo_1)) - - roll_2, solo_xyz_2, pointing_2 = get_hpc_info(start_utc_2, end_utc_2) - solo_2 = HeliographicStonyhurst(*solo_xyz_2, obstime=peak_utc_2, representation_type="cartesian") - center_hpc = SkyCoord(0 * u.deg, 0 * u.deg, frame=Helioprojective(obstime=peak_utc_2, observer=solo_2)) - coord_stix_2 = center_hpc.transform_to(STIXImaging(obstime=start_utc_2, obstime_end=end_utc_2, observer=solo_2)) - - coord_hp_2 = coord_stix_2.transform_to(Helioprojective(observer=solo_2)) - - # ========================================================================== - # APPROACH: Store individual coordinate components instead of SkyCoord - # ========================================================================== - # This is simpler and more transparent than storing SkyCoord objects: - # - No need for ICRS conversion or distance handling - # - FITS files contain plain columns with units (hp_tx, hp_ty, observer_x, etc.) - # - Reconstruction is straightforward: SkyCoord(Tx, Ty, frame=Helioprojective(...)) - # - Works for any coordinate frame (Helioprojective, HeliographicStonyhurst, etc.) - # - No need for fits.connect._encode_mixins() complexity - - # Store Helioprojective coordinates as separate Tx, Ty columns - hp_tx = [coord_hp.Tx, coord_hp_2.Tx] - hp_ty = [coord_hp.Ty, coord_hp_2.Ty] - - # Store observer position as separate x, y, z columns (or as single vector) - observer_x = [solo_xyz_1[0], solo_xyz_2[0]] - observer_y = [solo_xyz_1[1], solo_xyz_2[1]] - observer_z = [solo_xyz_1[2], solo_xyz_2[2]] - - data = QTable() - data["start_UTC"] = [start_utc, start_utc_2] - data["peak_UTC"] = [peak_utc, peak_utc_2] - data["end_UTC"] = [end_utc, end_utc_2] - data["counts"] = 100 * u.s - data["duration"] = 10 * u.s - data["lc_peak"] = [u.Quantity(range(1, 10))] - - # Store Helioprojective coordinates as separate columns with units - data["hp_tx"] = hp_tx - data["hp_ty"] = hp_ty - - # Store observer position as separate columns with units - data["observer_hgs_x"] = observer_x - data["observer_hgs_y"] = observer_y - data["observer_hgs_z"] = observer_z - - # Using fits.connect._encode_mixins (recommended approach used in STIXCore) - # This properly serializes complex types like SkyCoord - data_enc = fits.connect._encode_mixins(data) - data_hdu = table_to_hdu(data_enc) - data_hdu.name = "DATA" - - primary_hdu = fits.PrimaryHDU() - hdul = fits.HDUList([primary_hdu, data_hdu]) - hdul.writeto("test_sky_coord.fits", overwrite=True) - - # Read back - everything needed is stored in the FITS file - # No SkyCoord objects involved - just plain columns with units! - tr_conv = QTable.read("test_sky_coord.fits", astropy_native=True) - - print("\n" + "=" * 60) - print("READING BACK FROM FITS - NO SKYCOORD OBJECTS") - print("=" * 60) - - # Get original coordinates for comparison - original_coords_hp = [coord_hp, coord_hp_2] - - # Process all rows - for i in range(len(tr_conv)): - # Read stored Helioprojective coordinates (just plain columns with units!) - hp_tx_stored = tr_conv["hp_tx"][i] - hp_ty_stored = tr_conv["hp_ty"][i] - obstime = tr_conv["peak_UTC"][i] - - # Read observer position from separate columns - obs_x = tr_conv["observer_hgs_x"][i] - obs_y = tr_conv["observer_hgs_y"][i] - obs_z = tr_conv["observer_hgs_z"][i] - - # Reconstruct observer from stored components (NO SPICE!) - observer = HeliographicStonyhurst(x=obs_x, y=obs_y, z=obs_z, obstime=obstime, representation_type="cartesian") - observer_spherical = observer.represent_as("spherical") - - # Reconstruct Helioprojective SkyCoord from stored Tx, Ty - # This is the only place we create a SkyCoord - from plain stored values! - coord_hp_reconstructed = SkyCoord( - hp_tx_stored, hp_ty_stored, frame=Helioprojective(obstime=obstime, observer=observer) - ) - - print(f"\nRow {i + 1}:") - print(f" Stored Tx, Ty: {hp_tx_stored:.2f}, {hp_ty_stored:.2f}") - print( - f" Reconstructed Helioprojective: Tx={coord_hp_reconstructed.Tx:.2f}, Ty={coord_hp_reconstructed.Ty:.2f}" - ) - print(f" Observer: x={obs_x:.2e}, y={obs_y:.2e}, z={obs_z:.2e}") - print( - f" Observer (Spherical): Lon={observer_spherical.lon:.2f}, Lat={observer_spherical.lat:.2f}, " - f"R={observer_spherical.distance.to(u.AU):.3f}" - ) - - # Verify round-trip accuracy - original_hp = original_coords_hp[i] - - # Tolerance for angular comparison (1 arcsecond is reasonable for stored Tx/Ty) - tolerance = 1 * u.arcsec - - # Compare Tx/Ty directly - assert np.abs(coord_hp_reconstructed.Tx - original_hp.Tx) < tolerance, ( - f"Row {i + 1} HP Tx mismatch: {coord_hp_reconstructed.Tx} vs {original_hp.Tx}" - ) - assert np.abs(coord_hp_reconstructed.Ty - original_hp.Ty) < tolerance, ( - f"Row {i + 1} HP Ty mismatch: {coord_hp_reconstructed.Ty} vs {original_hp.Ty}" - ) - - print(" ✓ Round-trip verification passed:") - print(f" - Original HP: Tx={original_hp.Tx:.2f}, Ty={original_hp.Ty:.2f}") - print(f" - Reconstructed HP: Tx={coord_hp_reconstructed.Tx:.2f}, Ty={coord_hp_reconstructed.Ty:.2f}") - print( - f" - Difference: ΔTx={np.abs(coord_hp_reconstructed.Tx - original_hp.Tx):.4f}, " - f"ΔTy={np.abs(coord_hp_reconstructed.Ty - original_hp.Ty):.4f}" - ) - - print("\n" + "=" * 60) - print("✓ All conversions completed without SPICE!") - print("=" * 60) - print("\nBenefits of storing components instead of SkyCoord:") - print(" • No need for SkyCoord serialization - just plain columns with units") - print(" • FITS files are more transparent and easier to inspect") - print(" • Simple reconstruction: just create SkyCoord from stored Tx, Ty") - print(" • Works with any coordinate frame (Helioprojective, HeliographicStonyhurst, etc.)") - print(" • Observer position stored as separate x, y, z columns") - print(" • All units are preserved automatically by astropy.table") - print("\nWhat's stored in FITS:") - print(" • hp_tx, hp_ty: Helioprojective coordinates (arcsec)") - print(" • observer_hgs_x, observer_hgs_y, observer_hgs_z: Observer position (km)") - print(" • peak_UTC: Observation time") - print("\nReconstruction is simple:") - print(" 1. Read Tx, Ty from columns (they have units!)") - print(" 2. Read observer x, y, z from columns") - print(" 3. Create observer: HeliographicStonyhurst(x=x, y=y, z=z, obstime=obstime)") - print(" 4. Create coordinate: SkyCoord(Tx, Ty, frame=Helioprojective(obstime, observer))") - print(" 5. Done! No SPICE, no ICRS conversion, no distance issues.") diff --git a/stixcore/util/scripts/timecolumn.py b/stixcore/util/scripts/timecolumn.py deleted file mode 100644 index 19d9fac5..00000000 --- a/stixcore/util/scripts/timecolumn.py +++ /dev/null @@ -1,68 +0,0 @@ -from pathlib import Path -from datetime import datetime - -from stixpy.coordinates.transforms import STIXImaging, get_hpc_info -from sunpy.coordinates import HeliographicStonyhurst, Helioprojective - -import astropy.units as u -from astropy.coordinates import SkyCoord -from astropy.table.table import QTable -from astropy.time import Time - -from stixcore.config.config import CONFIG -from stixcore.ephemeris.manager import Spice, SpiceKernelManager -from stixcore.io.product_processors.fits.processors import FitsL3Processor -from stixcore.products.level3.flarelist import FlarelistSDC -from stixcore.products.product import Product - -if __name__ == "__main__": - _spm = SpiceKernelManager(Path(CONFIG.get("Paths", "spice_kernels"))) - spicemeta = _spm.get_latest_mk_and_pred() - - Spice.instance = Spice(spicemeta) - - l3_fits_writer = FitsL3Processor(Path(".")) - - data = QTable() - start_utc = Time(datetime.now()) - 100 * u.d - peak_utc = Time(datetime.now()) - 98 * u.d - end_utc = Time(datetime.now()) - 96 * u.d - - start_utc_2 = Time(datetime.now()) - 10 * u.d - peak_utc_2 = Time(datetime.now()) - 9 * u.d - end_utc_2 = Time(datetime.now()) - 9 * u.d - - data["start_UTC"] = [start_utc, start_utc_2] - data["peak_UTC"] = [peak_utc, peak_utc_2] - data["end_UTC"] = [end_utc, end_utc_2] - data["counts"] = 100 * u.s - - data["duration"] = 10 * u.s - data["lc_peak"] = [u.Quantity(range(1, 10))] - - control = QTable() - month = datetime(2025, 1, 1, 0, 0, 0, 0) - p_new = FlarelistSDC(data=data, month=month, control=control, energy=None) - - roll, solo_xyz, pointing = get_hpc_info(start_utc, end_utc) - solo = HeliographicStonyhurst(*solo_xyz, obstime=peak_utc, representation_type="cartesian") - center_hpc = SkyCoord(0 * u.deg, 0 * u.deg, frame=Helioprojective(obstime=peak_utc, observer=solo)) - coord_stix = center_hpc.transform_to(STIXImaging(obstime=start_utc, obstime_end=end_utc, observer=solo)) - - coord_hp = coord_stix.transform_to(Helioprojective(observer=solo)) - - roll, solo_xyz, pointing = get_hpc_info(start_utc_2, end_utc_2) - solo = HeliographicStonyhurst(*solo_xyz, obstime=peak_utc_2, representation_type="cartesian") - center_hpc = SkyCoord(0 * u.deg, 0 * u.deg, frame=Helioprojective(obstime=peak_utc_2, observer=solo)) - coord_stix_2 = center_hpc.transform_to(STIXImaging(obstime=start_utc_2, obstime_end=end_utc_2, observer=solo)) - - coord_hp_2 = coord_stix_2.transform_to(Helioprojective(observer=solo)) - - data["flare_location_stix"] = [coord_stix, coord_stix] - data["flare_location_hp"] = [coord_hp, coord_hp_2] - - f = l3_fits_writer.write_fits(p_new) - - p_in = Product(f[0]) - - print("all done")