#
# Copyright 2026 by University College Dublin. All rights reserved.
#
# Developed by: Derek O'Callaghan
# University College Dublin
# https://www.ucd.ie/
#
# Builds on:
# Gamma-ray Data Tools - Core Components (https://github.com/USRA-STI/gdt-core)
# Gamma-ray Data Tools - Fermi mission components (https://github.com/USRA-STI/gdt-fermi/)
#
# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except
# in compliance with the License. You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software distributed under the License
# is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or
# implied. See the License for the specific language governing permissions and limitations under the
# License.
#
import warnings
import astropy.io.fits as fits
import numpy as np
from scipy.stats import chi2
import scipy.interpolate
import healpy as hp
from astropy.coordinates import get_sun, SkyCoord, angular_separation
from astropy.coordinates.representation import CartesianRepresentation
from astropy.units import Quantity
from gdt.core.coords import Quaternion
from gdt.missions.fermi.gbm.localization import GbmHealPix
from gdt.missions.gifts.detectors import GiftsDetectors
from gdt.missions.gifts.frame import *
from gdt.missions.gifts.headers import HealpixHeaders
from gdt.missions.gifts.time import Time
__all__ = ['GiftsHealPix']
[docs]
class GiftsHealPix(GbmHealPix):
"""Class for GIFTS HEALPix localization files.
"""
[docs]
@classmethod
def from_data(cls, prob_arr, sig_arr, trigtime=None, headers=None, filename=None,
quaternion=None, scpos=None):
"""Create a GiftsHealPix object from a healpix array.
Args:
prob_arr (np.array):
The HEALPix array containing the probability/pixel
sig_arr (np.array):
The HEALPix array containing the significance/pixel
trigtime (float, optional):
The time corresponding to the localization
headers (:class:`~gdt.core.headers.FileHeaders`, optional):
The file headers
filename (str, optional): The filename
quaternion (:class:`~gdt.core.coords.Quaternion`, optional):
The associated spacecraft quaternion used to determine the
detector pointings in equatorial coordinates
scpos (astropy.coordinates.representation.CartesianRepresentation, optional):
The associated spacecraft position in Earth inertial coordinates
used to determine the geocenter location in equatorial
coordinates
Returns:
(:class:`GiftsHealPix`)
"""
# TODO: add sig_arr
obj = super(GbmHealPix, cls).from_data(prob_arr=prob_arr,
trigtime=trigtime,
filename=filename)
obj._sig = obj._assert_sig(sig_arr)
# Due to hardcoded HealpixHeaders, FermiFrame and GbmDetectors in GbmHealPix.from_data()
# need to replicate it here
if headers is not None:
if not isinstance(headers, HealpixHeaders):
raise TypeError('headers must be a HealpixHeaders object')
else:
headers = cls._none_default_headers()
obj._headers = headers
if quaternion is not None:
if not isinstance(quaternion, Quaternion):
raise TypeError('quaternion must be a Quaternion object')
obj._quat = quaternion
if scpos is not None:
if not isinstance(scpos, CartesianRepresentation):
raise TypeError('scpos must be a CartesianRepresentation object')
obj._scpos = scpos
# Two trigtime representations are used, float (input)
# and Time instance (derived)
# if we have a trigtime, calculate sun position
if trigtime is not None:
trigtime = Time(trigtime, format='gifts')
obj._sun_loc = get_sun(trigtime)
elif obj._headers[1]['SUN_RA'] is not None:
obj._sun_loc = SkyCoord(obj._headers[1]['SUN_RA'],
obj._headers[1]['SUN_DEC'], unit='deg',
frame='gcrs')
else:
obj._sun_loc = None
# create the spacecraft frame with the info that we have at hand
obj._frame = GiftsFrame(obstime=trigtime,
quaternion=obj._quat, obsgeoloc=obj._scpos,
detectors=GiftsDetectors)
# if have scpos, create geocenter location and earth radius
# if not, then try to pull from header
if obj._scpos is not None:
obj._geo_loc = obj._frame.geocenter
obj._geo_rad = obj._frame.earth_angular_radius
elif obj._headers[1]['GEO_RA'] is not None:
obj._geo_loc = SkyCoord(obj._headers[1]['GEO_RA'],
obj._headers[1]['GEO_DEC'], unit='deg',
frame='gcrs')
obj._geo_rad = Quantity(obj._headers[1]['GEO_RAD'], unit='deg')
# if have trigtime and quaternion, create detector pointings
# if not, then try to pull from header
# TODO: not sure I actually need the full block here
# TODO: at a minimum, remove duplication so that setattr, SkyCoord only called once
# Seems to be rotating the detector coordinates according to the quaternion
# I'm inclined to just have the first block and remove the second,
# i.e. a quaternion will always be provided, even if it's no rotation
# In fact, the latter can be the default
if (trigtime is not None) and (quaternion is not None):
for det in obj._frame.detectors:
pointing = (det.azimuth, det.elevation)
det_coord = SkyCoord(*pointing, frame=obj._frame).gcrs[0]
setattr(obj, det.name + '_pointing', det_coord)
elif obj._headers[1]['N0_RA'] is not None:
for det in GiftsDetectors:
ra_key = det.name + '_RA'
dec_key = det.name + '_DEC'
det_coord = SkyCoord(obj._headers[1][ra_key],
obj._headers[1][dec_key], unit='deg',
frame='gcrs')
setattr(obj, det.name + '_pointing', det_coord)
# build headers
obj._headers = obj._build_headers(obj.trigtime, obj.nside)
return obj
[docs]
@classmethod
def open(cls, file_path, **kwargs):
"""Open a GIFTS HEALPix FITS file and return the GiftsHealPix object
Args:
file_path (str): The file path of the FITS file
Returns:
(:class:`GiftsHealPix`)
"""
# As GbmHealPix.open() hardcodes the header class
# call this first, then replace with the corresponding
# GIFTS headers
obj = GbmHealPix.open(file_path, **kwargs)
headers = HealpixHeaders.from_file_headers(obj.headers)
# Need to retain sig, and TODO: add it to GiftsHealPix.from_data(),
# especially as GbmHealPix.from_data() is largely replicated
sig = obj._sig
obj.close()
obj = cls.from_data(prob_arr=obj.prob, sig_arr=sig, trigtime=obj._trigtime,
quaternion=obj._quat, scpos=obj._scpos, filename=obj.filename,
headers=headers)
# Unsure why sig isn't set in from_data(), so have to explicitly set it
# here as in GbmHealPix.open()
#obj._sig = sig
return obj
def _build_headers(self, trigtime, nside):
# Remove quaternion to avoid GbmDetector error in parent
quaternion = self._quat
self._quat = None
headers = super()._build_headers(trigtime, nside)
# Restore quaternion
self._quat = quaternion
if self.quaternion is not None:
headers['HEALPIX']['COMMENT'][1] = 'QUAT: ' + np.array2string(self.quaternion.scalar_last)
for det in GiftsDetectors:
pointing = getattr(self, det.name + '_pointing')
headers['HEALPIX'][det.name + '_RA'] = pointing.ra.value
headers['HEALPIX'][det.name + '_DEC'] = pointing.dec.value
return headers
@staticmethod
def _none_default_headers():
hdr = HealpixHeaders()
hdr[0]['TRIGTIME'] = None
for hk in [f"{k}_{c}" for k in ["SUN", "GEO"] + [d.name for d in GiftsDetectors] for c in ["RA", "DEC"]] + ["GEO_RAD"]:
hdr[1][hk] = None
return hdr