search

USRA-STI / gdt-fermi

Gamma-ray Data Tools - Fermi mission components
Apache License 2.0
2 stars 3 forks source link

GBM localization HEALPix FITS without COMMENT quaternion, spacecraft position can't be plotted with `EquatorialPlot` #14

Open derekocallaghan opened 8 months ago

derekocallaghan commented 8 months ago

Certain GBM localizations can't be plotted with EquatorialPlot, similar to the documentation example. For example, using a local copy of glg_healpix_all_bn180112842_v00.fit:

from gdt.core.coords import Quaternion, SpacecraftFrame
from gdt.core.plot.sky import EquatorialPlot
from gdt.missions.fermi.gbm.localization import GbmHealPix
from gdt.missions.fermi.gbm.trigdat import Trigdat
from gdt.missions.fermi.time import Time

gbm_healpix = GbmHealPix.open("./glg_healpix_all_bn180112842_v00.fit")

eqplot = EquatorialPlot()
eqplot.add_localization(gbm_healpix)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [22], in <cell line: 2>()
      1 eqplot = EquatorialPlot()
----> 2 eqplot.add_localization(gbm_healpix)

File ~/gitrepos/gifts-6u/gifts-software/gdt/core/plot/sky.py:323, in SkyPlot.add_localization(self, hpx, gradient, clevels, sun, earth, detectors, galactic_plane)
    320     clevels = [0.997, 0.955, 0.687]
    322 if detectors == 'all':
--> 323     detectors = [det.name for det in hpx.frame.detectors]
    324 else:
    325     if isinstance(detectors, str):

AttributeError: 'NoneType' object has no attribute 'detectors'

The opened HEALPix currently has no frame:

gbm_healpix.frame is None
True

Looking at GbmHealpix.from_data(), a frame is only created when trigtime, scpos, and quaternion are specified:

https://github.com/USRA-STI/gdt-fermi/blob/d4bb2b563fdda47aa420a76cef4820621004424b/src/gdt/missions/fermi/gbm/localization.py#L213-L219

GbmHealpix.open() expects quaternion and scpos to be stored in the FITS COMMENT header:

https://github.com/USRA-STI/gdt-fermi/blob/d4bb2b563fdda47aa420a76cef4820621004424b/src/gdt/missions/fermi/gbm/localization.py#L348-L361

However, these are missing from this localization. Consequently, no frame is created, and the detectors can't be plotted due to the error above:

gbm_healpix.headers[1]

PIXTYPE = 'HEALPIX '           / HEALPIX pixelization                           
ORDERING= 'NESTED  '           / Pixel ordering scheme, either RING or NESTED   
COORDSYS= 'C       '           / Ecliptic, Galactic or Celestial (equatorial)   
EXTNAME = 'HEALPIX '           / name of this binary table extension            
NSIDE   =                  128 / Resolution parameter of HEALPIX                
FIRSTPIX=                    0 / First pixel # (0 based)                        
LASTPIX =               196607 / Last pixel # (0 based)                         
INDXSCHM= 'IMPLICIT'           / Indexing: IMPLICIT or EXPLICIT                 
OBJECT  = 'GRB180112842'       / Sky coverage, either FULLSKY or PARTIAL        
SUN_RA  =    294.1123920790658 / RA of Sun                                      
SUN_DEC =   -21.58754321200521 / Dec of Sun                                     
GEO_RA  =    281.2035873034539 / RA of Geocenter relative to Fermi              
GEO_DEC =   -25.49632766235672 / Dec of Geocenter relative to Fermi             
GEO_RAD =                 67.5 / Radius of the Earth                            
N0_RA   =    23.59186671612186 / RA pointing for detector n0                    
N0_DEC  =    69.77534070454293 / Dec pointing for detector n0                   
N1_RA   =    357.5663742293567 / RA pointing for detector n1                    
N1_DEC  =    48.44799927945594 / Dec pointing for detector n1                   
N2_RA   =    356.0821947483266 / RA pointing for detector n2                    
N2_DEC  =    2.031334870128302 / Dec pointing for detector n2                   
N3_RA   =    244.0724750337902 / RA pointing for detector n3                    
N3_DEC  =    57.86808851985533 / Dec pointing for detector n3                   
N4_RA   =    239.7218063400303 / RA pointing for detector n4                    
N4_DEC  =     11.7706142076168 / Dec pointing for detector n4                   
N5_RA   =    301.5350637607208 / RA pointing for detector n5                    
N5_DEC  =    13.30445107975222 / Dec pointing for detector n5                   
N6_RA   =    136.2278681786678 / RA pointing for detector n6                    
N6_DEC  =    60.84211079751814 / Dec pointing for detector n6                   
N7_RA   =    151.0283993825707 / RA pointing for detector n7                    
N7_DEC  =    36.85396162561901 / Dec pointing for detector n7                   
N8_RA   =    174.3529606601047 / RA pointing for detector n8                    
N8_DEC  =   -2.660678205035072 / Dec pointing for detector n8                   
N9_RA   =    76.50429475456934 / RA pointing for detector n9                    
N9_DEC  =    30.72804082095637 / Dec pointing for detector n9                   
NA_RA   =    60.20741564206021 / RA pointing for detector na                    
NA_DEC  =   -12.53686185800487 / Dec pointing for detector na                   
NB_RA   =    121.9461723321663 / RA pointing for detector nb                    
NB_DEC  =   -13.36910124803439 / Dec pointing for detector nb                   
B0_RA   =    298.0882301495907 / RA pointing for detector b0                    
B0_DEC  =     13.4533919264521 / Dec pointing for detector b0                   
B1_RA   =    118.0882301495907 / RA pointing for detector b1                    
B1_DEC  =   -13.45339192645209 / Dec pointing for detector b1                   
COMMENT                                                                         
COMMENT

Note that the burst localization used in the documentation example does contain the required COMMENT header, resulting in frame creation and a successful plot:

docs_gbm_healpix = GbmHealPix.open("./glg_healpix_all_bn190915240_v00.fit")
docs_gbm_healpix.headers[1]

PIXTYPE = 'HEALPIX '           / HEALPIX pixelization                           
ORDERING= 'NESTED  '           / Pixel ordering scheme, either RING or NESTED   
COORDSYS= 'C       '           / Ecliptic, Galactic or Celestial (equatorial)   
EXTNAME = 'HEALPIX '           / name of this binary table extension            
NSIDE   =                  128 / Resolution parameter of HEALPIX                
FIRSTPIX=                    0 / First pixel # (0 based)                        
LASTPIX =               196608 / Last pixel # (0 based)                         
INDXSCHM= 'IMPLICIT'           / Indexing: IMPLICIT or EXPLICIT                 
OBJECT  = 'GRB190915240'       / Sky coverage, either FULLSKY or PARTIAL        
SUN_RA  =    172.5011935415178 / RA of Sun                                      
SUN_DEC =     3.23797213866954 / Dec of Sun                                     
GEO_RA  =    319.8312390218318 / RA of Geocenter relative to Fermi              
GEO_DEC =    17.40612934717674 / Dec of Geocenter relative to Fermi             
GEO_RAD =     67.2950460311874 / Radius of the Earth                            
N0_RA   =    146.5959532829778 / RA pointing for detector n0                    
N0_DEC  =    36.96759511828569 / Dec pointing for detector n0                   
N1_RA   =    177.8627043286993 / RA pointing for detector n1                    
N1_DEC  =    37.86505650718371 / Dec pointing for detector n1                   
N2_RA   =    235.2132411789626 / RA pointing for detector n2                    
N2_DEC  =    32.53574484506952 / Dec pointing for detector n2                   
N3_RA   =    141.1380239565408 / RA pointing for detector n3                    
N3_DEC  =   -11.86175670304888 / Dec pointing for detector n3                   
N4_RA   =    148.7150050771425 / RA pointing for detector n4                    
N4_DEC  =    -57.7151024222104 / Dec pointing for detector n4                   
N5_RA   =    204.6679387636932 / RA pointing for detector n5                    
N5_DEC  =   -14.17362482787578 / Dec pointing for detector n5                   
N6_RA   =    103.6451496378668 / RA pointing for detector n6                    
N6_DEC  =    19.90603710061072 / Dec pointing for detector n6                   
N7_RA   =    82.17216524135398 / RA pointing for detector n7                    
N7_DEC  =    4.861082974313736 / Dec pointing for detector n7                   
N8_RA   =    53.77139403781888 / RA pointing for detector n8                    
N8_DEC  =   -31.16408912405502 / Dec pointing for detector n8                   
N9_RA   =    76.25309367118034 / RA pointing for detector n9                    
N9_DEC  =    65.37497944156361 / Dec pointing for detector n9                   
NA_RA   =    329.2307441352281 / RA pointing for detector na                    
NA_DEC  =    56.85788605114446 / Dec pointing for detector na                   
NB_RA   =    24.77833487844497 / RA pointing for detector nb                    
NB_DEC  =    13.78283005529111 / Dec pointing for detector nb                   
B0_RA   =    203.0460127070913 / RA pointing for detector b0                    
B0_DEC  =    -17.1448614419065 / Dec pointing for detector b0                   
B1_RA   =    23.04601270709126 / RA pointing for detector b1                    
B1_DEC  =    17.14486144190651 / Dec pointing for detector b1                   
COMMENT SCPOS: [-5039500.  4254000. -2067500.]                                  
COMMENT QUAT: [-0.223915  0.447149  0.860062 -0.101055]
docs_gbm_healpix.frame

<SpacecraftFrame: 1 frames;
 obstime=[590219102.911008]
 obsgeoloc=[(-5039500., 4254000., -2067500.) m]
 obsgeovel=[(0., 0., 0.) m / s]
 quaternion=[(x, y, z, w) [-0.223915,  0.447149,  0.860062, -0.101055]]>

I've been using a local workaround where I retrieve a trigtime SpacecraftFrame from the corresponding Trigdat and specifiy it to the GbmHealpix, e.g.:

trigdat = Trigdat.open("./glg_trigdat_all_bn180112842_v01.fit")
gbm_healpix._frame = trigdat.poshist.at(Time(gbm_healpix.trigtime, format="fermi"))
gbm_healpix.frame

<SpacecraftFrame: 1 frames;
 obstime=[537480753.340436]
 obsgeoloc=[(-1210312.37296423, 6115468.77964168, 2972968.74576547) m]
 obsgeovel=[(0., 0., 0.) m / s]
 quaternion=[(x, y, z, w) [ 0.10279574,  0.0747444 ,  0.51519948, -0.84759413]]>

This will ensure that EquatorialPlot.add_localization() will work as expected.

However, GbmHealpix.from_data() also loads the detector coordinates from the header, similar to how this was previously done in GBM Data Tools. Perhaps these should be used as a backup in situations where no frame is created due to empty COMMENT headers:

https://github.com/USRA-STI/gdt-fermi/blob/d4bb2b563fdda47aa420a76cef4820621004424b/src/gdt/missions/fermi/gbm/localization.py#L236-L242

derekocallaghan commented 4 months ago

Hi @AdamGoldstein-USRA,

Just wanted to check whether the suggested GbmHealpix.from_data() workaround is okay with you? If so, I'll create a PR.

Thanks