Giter VIP home page Giter VIP logo

gdt-fermi's People

Contributors

adamgoldstein-usra avatar billcleveland-usra avatar joshuarwood avatar

Stargazers

avatar

Watchers

avatar avatar

Forkers

gifts-6u sumanbala2210-usra joshuarwood

gdt-fermi's Issues

Add datetime conversion with timezone info for Time.py

For example, when converting an MET to datetime, I use

from gdt.missions.fermi.time import Time
t1 = Time(trigtime, format='fermi', scale='utc').datetime

This object does not have timezone information and when comparing to another datetime object, e.g.,

t2 = datetime.datetime(y, m, d, dt, tzinfo=datetime.timezone.utc),

I receive this error: TypeError: can't compare offset-naive and offset-aware datetimes

However, in the GBM Data Tools, when using
from gbm.time import Met
t = Met(met)

This automatically returns a datetime object with the UTC timezone specified:
self.__time.utc.to_datetime(datetime.timezone.utc)

It would be nice to have this capability in the GDT without having to call
from gdt.missions.fermi.time import Time
t1 = Time(trigtime, format='fermi', scale='utc')
t1 = t1.to_datetime(datetime.timezone.utc)

Add the ability to handle relative time and override CREATOR keyword

Currently when the GbmTte class is being used to generate a new file, there is no method to specify that the time should be relative to trigger time. In addition, if 'CREATOR' is specified within the given header, it should override the one used in the PRIMARY header (and not show up in the extensions).

The current version of GbmTte fails to make the required TZEROn keywords for EVENTS and GTI.

`Attribute Error` when calling `SpectralFitterPgstat.fit()`

When performing a spectral fit (I've attached a portion of the code here), I get an attribute error:

from gdt.core.spectra.fitting import SpectralFitterPgstat
from gdt.core.spectra.functions import PowerLaw, Comptonized, Band
from gdt.missions.fermi.gbm.response import GbmRsp2
from gdt.missions.fermi.gbm.collection import GbmDetectorCollection

# Load a Response (.rsp) file for each NaI and BGO detector
nai_rsps = [GbmRsp2.open(str(gbm_path)+'/glg_cspec_'+nai+'_bn'+str(bn)+'_v01.rsp2') for nai in nais]
bgo_rsps = [GbmRsp2.open(str(gbm_path)+'/glg_cspec_'+bgo+'_bn'+str(bn)+'_v01.rsp2') for bgo in bgos]

# Create a collection from the list of our files
all_rsps = nai_rsps + bgo_rsps
rsps = GbmDetectorCollection.from_list(all_rsps)

# Initialize spectral fitter with the PHAs, backgrounds, and responses:
specfitter = SpectralFitterPgstat(phas, bkgds.to_list(), rsps.to_list(), method='TNC')

# Instantiate spectral functions
# a power law, cut-off power law, and a Band function
pl  = PowerLaw()
comp = Comptonized()
band  = Band()
models = [pl, comp, band]

for model in models:
    specfitter.fit(model, options={'maxiter': 1000})

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[68], line 5
      2 print ('Fitting with {} model ...'.format(model.name))
      4 # Perform fit 
----> 5 specfitter.fit(model, options={'maxiter': 1000})


File ~/anaconda3/envs/gdt_devel/lib/python3.12/site-packages/gdt/core/spectra/fitting.py:769, in SpectralFitter._fold_model(self, function, params)
    766 stat = np.zeros(self.num_sets)
    767 for i in range(self.num_sets):
    768     # fold model through response and convert to raw model counts
--> 769     model = self._rsp[i].drm.fold_spectrum(function, params, channel_mask=self._chan_masks[i])
    771     # perform likelihood calculation for one dataset
    772     stat[i] = self._eval_stat(i, model)

AttributeError: 'GbmRsp2' object has no attribute 'drm'

Feature request: add HTTPS support for HEASARC archive access

FTP access to HEASARC is currently unreliable due to frequent outages at the HEASARC FTP servers. HTTPS access currently appears more reliable. We should add a method to access archival HEASARC data through the HTTPS servers. I can start working on this in May 2024.

`AttributeError` when calling `GbmHealpix._repr_html()`

The following occurs in a notebook:

from gdt.missions.fermi.gbm.localization import GbmHealPix

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

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
File ~/anaconda3/envs/gifts/lib/python3.10/site-packages/IPython/core/formatters.py:343, in BaseFormatter.__call__(self, obj)
    341     method = get_real_method(obj, self.print_method)
    342     if method is not None:
--> 343         return method()
    344     return None
    345 else:

File ~/gitrepos/gifts-6u/gifts-software/gdt/core/file.py:217, in FitsFileContextManager._repr_html_(self)
    215 s += '<table>'
    216 s += '<tr><th>No.</th><th>Name</th><th>Ver</th><th>Type</th><th>Cards</th><th>Dimensions</th></tr>'
--> 217 for row in self._hdulist.info(False):
    218     s += f'<tr><td>{row[0]}</td><td>{row[1]}</td><td>{row[2]}</td><td>{row[3]}</td><td>{row[4]}</td>' \
    219          f'<td>{row[5]}</td></tr>'
    220 s += '</table>'

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

_hdulist is instantiated in FitsFileContextManager.open():

    def open(cls, file_path: Union[str, Path], mode: str = 'readonly', memmap: bool = None):
        """Open a FITS file.
            
        Args:
            file_path (str): The file path
            mode (str): The file access mode
            memmap (bool): If True, memory map when reading the file
        
        Returns:
            (:class:`FitsFileContextManager`)
        """
        path = Path(file_path)
        obj = cls()
        obj._hdulist = fits.open(path, mode=mode, memmap=memmap)

which is called by GbmHealpix.open():

obj = super().open(file_path, **kwargs)

However, the object is recreated in GbmHealpix.from_data(), without _hdulist:

obj = super().from_data(prob_arr, trigtime=trigtime, filename=filename)

This local FitsFileContextManager._repr_html() fix resolves the issue, if it's acceptable I'll create a PR:

diff --git a/src/gdt/core/file.py b/src/gdt/core/file.py
index c35202c..1f8f2cd 100644
--- a/src/gdt/core/file.py
+++ b/src/gdt/core/file.py
@@ -214,7 +214,7 @@ class FitsFileContextManager(AbstractContextManager):
         s = f'<p>&lt{self.__class__.__name__}(filename=<b>"{self.filename}"</b>) at {hex(id(self))}&gt</p>'
         s += '<table>'
         s += '<tr><th>No.</th><th>Name</th><th>Ver</th><th>Type</th><th>Cards</th><th>Dimensions</th></tr>'
-        for row in self._hdulist.info(False):
+        for row in self.hdulist.info(False):
             s += f'<tr><td>{row[0]}</td><td>{row[1]}</td><td>{row[2]}</td><td>{row[3]}</td><td>{row[4]}</td>' \
                  f'<td>{row[5]}</td></tr>'
         s += '</table>'

SuperFunction() does not work properly when more than two functions are used

SuperFunction() does not work properly when more than two functions are used. Here is example code to reproduce the issue:

import matplotlib.pyplot as plt

from gdt.missions.fermi.gbm.phaii import GbmPhaii
from gdt.missions.fermi.gbm.response import GbmRsp2
from gdt.missions.fermi.gbm.collection import GbmDetectorCollection
from gdt.core.background.fitter import BackgroundFitter
from gdt.core.background.binned import Polynomial
from gdt.core.spectra.fitting import SpectralFitterPstat
from gdt.core.spectra.functions import SmoothlyBrokenPowerLaw, PowerLaw, GaussLine
from gdt.core.plot.model import ModelFit

if __name__ == "__main__":

    b0_phaii = GbmPhaii.open('./data/glg_cspec_b0_bn221009553_v00.pha')
    b1_phaii = GbmPhaii.open('./data/glg_cspec_b1_bn221009553_v00.pha')
    n4_phaii = GbmPhaii.open('./data/glg_cspec_n4_bn221009553_v00.pha')
    n6_phaii = GbmPhaii.open('./data/glg_cspec_n6_bn221009553_v00.pha')
    n8_phaii = GbmPhaii.open('./data/glg_cspec_n8_bn221009553_v00.pha')
    cspecs = GbmDetectorCollection.from_list([b0_phaii, b1_phaii, n4_phaii, n6_phaii, n8_phaii])

    b0_rsp2 = GbmRsp2.open('./data/glg_cspec_b0_bn221009553_v13.rsp2')
    b1_rsp2 = GbmRsp2.open('./data/glg_cspec_b1_bn221009553_v13.rsp2')
    n4_rsp2 = GbmRsp2.open('./data/glg_cspec_n4_bn221009553_v13.rsp2')
    n6_rsp2 = GbmRsp2.open('./data/glg_cspec_n6_bn221009553_v13.rsp2')
    n8_rsp2 = GbmRsp2.open('./data/glg_cspec_n8_bn221009553_v13.rsp2')
    rsps = GbmDetectorCollection.from_list([b0_rsp2, b1_rsp2, n4_rsp2, n6_rsp2, n8_rsp2])


    bkgd_times = [(-150.0, -2.0), (1620.0, 2370.0)]
    poly_orders = [4, 4, 3, 4, 4]
    src_time = (310, 320)
    erange_nai = (45.0, 900.0)
    erange_bgo = (400.0, 39_000.0)

    backfitters = [BackgroundFitter.from_phaii(cspec, Polynomial, time_ranges=bkgd_times) for cspec in cspecs]
    backfitters = GbmDetectorCollection.from_list(backfitters, dets=cspecs.detector())

    [backfitter.fit(order=poly_order) for backfitter, poly_order in zip(backfitters, poly_orders)]

    bkgds = backfitters.interpolate_bins(cspecs.data()[0].tstart, cspecs.data()[0].tstop)
    bkgds = GbmDetectorCollection.from_list(bkgds, dets=cspecs.detector())

    bkgd_specs = bkgds.integrate_time(*src_time)

    phas = cspecs.to_pha(time_ranges=src_time, nai_kwargs={'energy_range':erange_nai}, bgo_kwargs={'energy_range':erange_bgo})
    rsps_interp = [rsp.interpolate(pha.tcent) for rsp, pha in zip(rsps, phas)]

    print('\n----- SBPL + PL Fit -----')

    SBPL_PL = SmoothlyBrokenPowerLaw()+PowerLaw()
    SBPL_PL.default_values = [0.06, 12.0, 9.0, 0.5, -4.4, 100.0, 0.06, -1.7, 87.0]

    specfitter = SpectralFitterPstat(phas, bkgds.to_list(), rsps_interp, method='L-BFGS-B')
    specfitter.fit(SBPL_PL, options={'maxiter': 10_000})

    modelplot = ModelFit(fitter=specfitter)
    plt.savefig('./SBPL_PL_cnts.png')

    modelplot = ModelFit(fitter=specfitter, view='nufnu')
    modelplot.nufnu_spectrum(num_samples=500)
    modelplot.ax.grid(which='both')
    plt.savefig('./SBPL_PL_nuFnu.png')

    print('\n----- SBPL + PL + Gline Fit -----')

    # SBPL_PL_Gline = SmoothlyBrokenPowerLaw()+PowerLaw()+GaussLine()
    # SBPL_PL_Gline = (SmoothlyBrokenPowerLaw()+PowerLaw())+GaussLine()
    SBPL_PL = SmoothlyBrokenPowerLaw()+PowerLaw()
    SBPL_PL_Gline = SBPL_PL + GaussLine()

    SBPL_PL_Gline.default_values = [0.06, 12.0, 9.0, 0.5, -4.4, 100.0, 0.06, -1.7, 87.0, 0.11, 8_975.0, 2_649.0]

    specfitter = SpectralFitterPstat(phas, bkgds.to_list(), rsps_interp, method='L-BFGS-B')
    specfitter.fit(SBPL_PL_Gline, options={'maxiter': 10_000})

    modelplot = ModelFit(fitter=specfitter)
    plt.savefig('./SBPL_PL_Gline_cnts.png')

    modelplot = ModelFit(fitter=specfitter, view='nufnu')
    modelplot.nufnu_spectrum(num_samples=500)
    modelplot.ax.grid(which='both')
    plt.savefig('./SBPL_PL_Gline_nuFnu.png')

Here is the failure trace:

Traceback (most recent call last):
  File "/Users/stephenlesage/Fermi-GBM/projects/gaussian_line/python_sims/test.py", line 87, in <module>
    modelplot = ModelFit(fitter=specfitter, view='nufnu')
  File "/Users/stephenlesage/Environments/GDTtools/lib/python3.9/site-packages/gdt/core/plot/model.py", line 78, in __init__
    self.set_fit(fitter, resid=resid)
  File "/Users/stephenlesage/Environments/GDTtools/lib/python3.9/site-packages/gdt/core/plot/model.py", line 215, in set_fit
    self.nufnu_spectrum()
  File "/Users/stephenlesage/Environments/GDTtools/lib/python3.9/site-packages/gdt/core/plot/model.py", line 173, in nufnu_spectrum
    self._plot_spectral_model(**kwargs)
  File "/Users/stephenlesage/Environments/GDTtools/lib/python3.9/site-packages/gdt/core/plot/model.py", line 288, in _plot_spectral_model
    label=comps[i], color=self.colors[i+1],
IndexError: list index out of range

I tried the following syntax suggestions and got the same error message:

  1. SBPL_PL_Gline = SmoothlyBrokenPowerLaw()+PowerLaw()+GaussLine()
  2. SBPL_PL_Gline = (SmoothlyBrokenPowerLaw()+PowerLaw())+GaussLine()
  3. SBPL_PL = SmoothlyBrokenPowerLaw()+PowerLaw()
    SBPL_PL_Gline = SBPL_PL + GaussLine()

It appears to be an issue in model.py, specifically that the number of components "num_comp" does not reflect the number of function components "comps" being fed to SuperFunction(). Implementing print statements reveals:

num_comp = 3
comps = ['SmoothlyBrokenPowerLaw+PowerLaw', 'GaussLine']

I can provide the data files to run the example code if necessary.

Attribute error when using GbmDetectorCollection with TTE data: cannot find detector information

When I convert TTE data to PHAII or PHA files and group several detectors using GbmDetectorCollection, there seems to be missing values for local arguments. This procedure works for other PHAII types like CTIME and CSPEC. I tried converting to PHAII and then to PHA and I get the same error during the conversion to PHA. May it is an issue specific to the opening/reading of TTE files and converting them to PHAII/PHA? Might be related to issue #13?

Screenshot 2024-01-26 at 8 54 28 AM

Screenshot 2024-01-26 at 8 54 54 AM

When converting directly to PHA from TTE:

Screenshot 2024-01-26 at 9 43 20 AM

Screenshot 2024-01-26 at 9 43 39 AM

GbmHealPix.multiply results in a map that cannot be plotted

Using the multiply() member of the GbmHealPix class results in a map that cannot be plotted. Here is example code to reproduce the issue:

import matplotlib.pyplot as plt

from gdt.core.plot.sky import EquatorialPlot
from gdt.missions.fermi.gbm.localization import GbmHealPix

# open a file
g = GbmHealPix.open(path)

# multiple
g2 = g.multiply(g, g)

# try to plot
skyplot = EquatorialPlot()
skyplot.add_localization(g2) # <-- fails at this step
plt.show()

where path is a path to any valid GBM localization. Here is the failure trace:

Traceback (most recent call last):
  File "multiply_issue.py", line 18, in <module>
    skyplot.add_localization(g2)
  File "pyvenv-3.8.9/lib/python3.8/site-packages/gdt/core/plot/sky.py", line 323, in add_localization
    detectors = [det.name for det in hpx.frame.detectors]
AttributeError: 'NoneType' object has no attribute 'detectors'

It appears to be a result of the frame object not being generated. I suspect the multiply function might not be preserving scpos and quaternion fields. I'll take a closer look tomorrow. At the moment I can fix the error by adding the following line after multiplication:

g2 = GbmHealPix.from_data(g.prob, trigtime=g.trigtime, scpos=g.scpos, quaternion=g.quaternion)

This forces the creation of the frame object, but seems counterintuitive to me.

Possible bug with GbmDetectorCollection

I tried to convert TTE detector data to a PHAII format so that I can group detectors for spectral fitting using GbmDetectorCollection. I follow the examples in the GBM Data Tools and GDT documentation https://fermi.gsfc.nasa.gov/ssc/data/analysis/rmfit/gbm_data_tools/gdt-docs/notebooks/SpectralFitting.html but get errors when trying to fit a background to a list of detectors. However, working through the error seems to show an empty list after using GbmDetectorCollection. The conversion beforehand from TTE to PHAII seems to work fine, but the collection element does not seem to populate a list after. I could be doing something wrong, but I am struggling to see the issue and would appreciate some feedback. I am using a MacOS (Monterrey 12.6.5). Attached is a screenshot of the executed code.

Screen Shot 2023-05-08 at 12 38 11 PM

Screen Shot 2023-05-08 at 12 41 02 PM

Unable to open FITS files when certain `GbmHeader` keywords aren't in the corresponding header

I've noticed that certain triggers/burst FITS files don't contain all GbmHeader subclass keywords, preventing the subclass from opening the file. This can be illustrated with Tcat.open() and a local copy of glg_tcat_all_bn180112842_v00.fit, but it also happens in other situations e.g. GbmTte.open():

In [1]: from astropy.io import fits
   ...: from gdt.missions.fermi.gbm.headers import TcatHeader
   ...: from gdt.missions.fermi.gbm.tcat import Tcat

In [2]: tcat_file = "./glg_tcat_all_bn180112842_v00.fit"

In [3]: Tcat.open(tcat_file)
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Input In [3], in <cell line: 1>()
----> 1 Tcat.open(tcat_file)

File ~/gitrepos/gifts-6u/gdt-fermi/src/gdt/missions/fermi/gbm/tcat.py:88, in Tcat.open(cls, file_path, **kwargs)
     86 obj = super().open(file_path, **kwargs)
     87 hdrs = [hdu.header for hdu in obj.hdulist]
---> 88 obj._headers = TcatHeaders.from_headers(hdrs)        
     89 return obj

File ~/gitrepos/gifts-6u/gdt-core/src/gdt/core/headers.py:192, in FileHeaders.from_headers(cls, headers)
    190             hidx += 1
    191         else:
--> 192             obj[i][key] = headers[i][key]
    194 return obj

File ~/anaconda3/envs/gifts/lib/python3.10/site-packages/astropy/io/fits/header.py:170, in Header.__getitem__(self, key)
    167 else:
    168     keyword = key
--> 170 card = self._cards[self._cardindex(key)]
    172 if card.field_specifier is not None and keyword == card.rawkeyword:
    173     # This is RVKC; if only the top-level keyword was specified return
    174     # the raw value, not the parsed out float value
    175     return card.rawvalue

File ~/anaconda3/envs/gifts/lib/python3.10/site-packages/astropy/io/fits/header.py:1771, in Header._cardindex(self, key)
   1768         indices = self._rvkc_indices.get(keyword, None)
   1770 if not indices:
-> 1771     raise KeyError(f"Keyword {keyword!r} not found.")
   1773 try:
   1774     return indices[n]

KeyError: "Keyword 'INFILE01' not found."

The FITS file can be opened directly, where it can be seen that INFILE01 isn't in the header:

In [5]: sorted(TcatHeader().keys())
Out[5]: 
['ADC_HI',
 'ADC_LO',
 'CHAN_HI',
 'CHAN_LO',
 'CLASS',
 'CREATOR',
 'DATE',
 'DATE-END',
 'DATE-OBS',
 'DEC_OBJ',
 'DEC_SCX',
 'DEC_SCZ',
 'DET_MASK',
 'EQUINOX',
 'ERR_RAD',
 'FILENAME',
 'FILETYPE',
 'GCN_FLAG',
 'GEO_LAT',
 'GEO_LONG',
 'HISTORY',
 'HISTORY',
 'HISTORY',
 'HISTORY',
 'INFILE01',
 'INSTRUME',
 'LOC_ENRG',
 'LOC_SRC',
 'LOC_VER',
 'MJDREFF',
 'MJDREFI',
 'OBJECT',
 'OBJ_CLAS',
 'OBSERVER',
 'ORIGIN',
 'PHI',
 'RADECSYS',
 'RA_OBJ',
 'RA_SCX',
 'RA_SCZ',
 'RELIABLT',
 'TELESCOP',
 'THETA',
 'TIMESYS',
 'TIMEUNIT',
 'TRIGSCAL',
 'TRIGTIME',
 'TRIG_ALG',
 'TRIG_SIG',
 'TSTART',
 'TSTOP']

In [7]: tcatf = fits.open("./glg_tcat_all_bn180112842_v00.fit")

In [8]: 'INFILE01' in tcatf["PRIMARY"].header.keys()
Out[8]: False

In [9]: tcatf["PRIMARY"].header
Out[9]: 
SIMPLE  =                    T / file does conform to FITS standard             
BITPIX  =                    8 / number of bits per data pixel                  
NAXIS   =                    0 / number of data axes                            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
CREATOR = 'DOL+GBM_TRIGDAT_Reader v1.28' / Software and version creating file   
FILETYPE= 'TRIGGER ENTRY'      / Name for this type of FITS file                
TELESCOP= 'GLAST   '           / Name of mission/satellite                      
INSTRUME= 'GBM     '           / Specific instrument used for observation       
OBSERVER= 'Meegan  '           / GLAST Burst Monitor P.I.                       
ORIGIN  = 'GIOC    '           / Name of organization making file               
DATE    = '2018-01-12T20:22:25' / file creation date (YYYY-MM-DDThh:mm:ss UT)   
DATE-OBS= '2018-01-12T20:10:14' / Date of start of observation                  
DATE-END= '2018-01-12T20:20:28' / Date of end of observation                    
TIMESYS = 'TT      '           / Time system used in time keywords              
TIMEUNIT= 's       '           / Time since MJDREF, used in TSTART and TSTOP    
MJDREFI =                51910 / MJD of GLAST reference epoch, integer part     
MJDREFF = 7.428703703703703D-4 / MJD of GLAST reference epoch, fractional part  
TSTART  =     537480614.201896 / Time of first TRIGDAT INFO packet              
TSTOP   =     537481228.613448 / Time of last TRIGDAT INFO packet               
FILENAME= 'glg_tcat_all_bn180112842_v00.fit' / Name of this file                
TRIGTIME=     537480753.340436 / Trigger time relative to MJDREF, double precisi
OBJECT  = 'GRB180112842'       / Burst name in standard format, yymmddfff       
RADECSYS= 'FK5     '           / Stellar reference frame                        
EQUINOX =               2000.0 / Equinox for RA and Dec                         
RA_OBJ  =             173.9200 / Calculated RA of burst                         
DEC_OBJ =              23.0600 / Calculated Dec of burst                        
ERR_RAD =               6.6700 / Calculated Location Error Radius               
TRIG_SIG=                  4.9 / Trigger significance (sigma)                   
LOC_SRC = 'Fermi, GBM'         / Mission/Instrument providing the localization. 
LOC_VER = '3       '           / Version string for localization software.      
LOC_ENRG= '(50.0,300.0)'       / [keV] Energy range used for localization.      
GEO_LONG=              45.6500 / [deg] Spacecraft geographical east longitude   
GEO_LAT =              25.4833 / [deg] Spacecraft geographical north latitude   
RA_SCX  =             298.1209 / [deg] Pointing of spacecraft x-axis: RA        
DEC_SCX =              13.4377 / [deg] Pointing of spacecraft x-axis: Dec       
RA_SCZ  =              94.5567 / [deg] Pointing of spacecraft z-axis: RA        
DEC_SCZ =              75.3901 / [deg] Pointing of spacecraft z-axis: Dec       
THETA   =              65.0000 / [deg] Angle from spacecraft z-axis (zenith)    
PHI     =             243.0000 / [deg] Angle from spacecraft +x-axis toward +y  
CLASS   = 'GRB     '           / Classification of trigger.                     
OBJ_CLAS= 'GRB     '           / GBM internal classification.                   
RELIABLT=               0.9529 / Reliability of classification.                 
TRIGSCAL=                  256 / [ms] Triggered timescale                       
TRIG_ALG=                    9 / Triggered algorithm number                     
CHAN_LO =                    3 / Trigger channel: low                           
CHAN_HI =                    4 / Trigger channel: high                          
ADC_LO  =                  259 / Trigger channel: low (ADC: 0 - 4095)           
ADC_HI  =                 1352 / Trigger channel: high (ADC: 0 - 4095)          
DET_MASK= '00000011000000'     / Triggered detectors: (0-13)                    
GCN_FLAG= 'No      '           / Was this file used for auto GCN Notice?        
CHECKSUM= 'ZEmKc9jHZCjHb9jH'   / HDU checksum updated 2018-01-12T20:22:25       
DATASUM = '         0'         / data unit checksum updated 2018-01-12T20:22:25

I'm currently using a local gdt-core workaround as follows, where any expected keywords missing from the header are ignored. Although I can create a PR with this modification, I'm unsure whether you always want to reject any FITS files with missing keywords, so I've held off for now and am logging this issue instead.

diff --git a/src/gdt/core/headers.py b/src/gdt/core/headers.py
index 65e39a4..37c61e3 100644
--- a/src/gdt/core/headers.py
+++ b/src/gdt/core/headers.py
@@ -180,7 +180,8 @@ class FileHeaders():
         for i in range(num_headers):
             cidx = 0
             hidx = 0
-            for key in obj[i].keys():
+            #for key in obj[i].keys():
+            for key in [k for k in obj[i].keys() if k in headers[i]]:
                 if (key == 'COMMENT'):
                     obj[i][key][cidx] = headers[i][key][cidx]
                     cidx += 1

GEOS version dependency error

I received the following error during install. We may need to add a note about GEOS version dependency in the install instructions.

(gdt_public) Tachyon> pip install astro-gdt-fermi
Collecting astro-gdt-fermi
  Downloading astro_gdt_fermi-2.0.0-py3-none-any.whl (73 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 73.1/73.1 kB 1.2 MB/s eta 0:00:00
Collecting astro-gdt>=2.0.0 (from astro-gdt-fermi)
  Downloading astro_gdt-2.0.1-py3-none-any.whl (491 kB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 491.4/491.4 kB 6.6 MB/s eta 0:00:00
Collecting pyproj>=1.9.6 (from astro-gdt-fermi)
  Downloading pyproj-3.5.0-cp310-cp310-macosx_10_9_x86_64.whl (8.5 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 8.5/8.5 MB 10.9 MB/s eta 0:00:00
Collecting numpy>=1.17.3 (from astro-gdt-fermi)
  Downloading numpy-1.24.3-cp310-cp310-macosx_10_9_x86_64.whl (19.8 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 19.8/19.8 MB 13.3 MB/s eta 0:00:00
Collecting scipy>=1.1.0 (from astro-gdt-fermi)
  Using cached scipy-1.10.1-cp310-cp310-macosx_10_9_x86_64.whl (35.1 MB)
Collecting matplotlib>=3.7.1 (from astro-gdt-fermi)
  Downloading matplotlib-3.7.1-cp310-cp310-macosx_10_12_x86_64.whl (7.4 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 7.4/7.4 MB 13.9 MB/s eta 0:00:00
Collecting astropy>=3.1 (from astro-gdt-fermi)
  Downloading astropy-5.2.2-cp310-cp310-macosx_10_9_x86_64.whl (7.0 MB)
     ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 7.0/7.0 MB 8.3 MB/s eta 0:00:00
Collecting healpy>=1.12.4 (from astro-gdt-fermi)
  Using cached healpy-1.16.2-cp310-cp310-macosx_10_9_x86_64.whl (7.0 MB)
Collecting cartopy>=0.21.1 (from astro-gdt-fermi)
  Using cached Cartopy-0.21.1.tar.gz (10.9 MB)
  Installing build dependencies ... done
  Getting requirements to build wheel ... error
  error: subprocess-exited-with-error
  
  × Getting requirements to build wheel did not run successfully.
  │ exit code: 1
  ╰─> [1 lines of output]
      GEOS version 3.6.1 is installed, but cartopy requires at least version 3.7.2.
      [end of output]
  
  note: This error originates from a subprocess, and is likely not a problem with pip.
error: subprocess-exited-with-error

× Getting requirements to build wheel did not run successfully.
│ exit code: 1
╰─> See above for output.

note: This error originates from a subprocess, and is likely not a problem with pip.

Time.gbm_bn does not work in Astropy 6.0

Example:

from gdt.missions.fermi.time import Time
Time.now().gbm_bn

Results in the error

AttributeError: 'GbmBurstNumber' object has no attribute 'mask_if_needed'

This appears to be due to a change in Astropy 6.0. I have verified this works in Astropy 5.2.2.

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

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:

gdt-fermi/src/gdt/missions/fermi/gbm/localization.py

Lines 213 to 219 in d4bb2b5

if (trigtime is not None) and (scpos is not None) and \
(quaternion is not None):
obj._scpos = scpos
obj._quat = Quaternion(quaternion)
obj._frame = SpacecraftFrame(obstime=Time(trigtime, format='fermi'),
quaternion=obj._quat,

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

gdt-fermi/src/gdt/missions/fermi/gbm/localization.py

Lines 348 to 361 in d4bb2b5

# quaternion and scpos are stored as comments
try:
headers[1]['COMMENT'][0] = obj.hdulist[1].header['COMMENT'][0]
headers[1]['COMMENT'][1] = obj.hdulist[1].header['COMMENT'][1]
except:
headers[1]['COMMENT'][0] = ''
headers[1]['COMMENT'][1] = ''
try:
scpos_comment = headers[1]['COMMENT'][0]
scpos = scpos_comment.split('[')[1].split(']')[0]
scpos = np.array([float(el) for el in scpos.split()])
except:
scpos = None

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:

gdt-fermi/src/gdt/missions/fermi/gbm/localization.py

Lines 236 to 242 in d4bb2b5

for det in GbmDetectors:
ra_key = det.name.upper() + '_RA'
dec_key = det.name.upper() + '_DEC'
det_coord = SkyCoord(obj._headers[1][ra_key],
obj._headers[1][dec_key], unit='deg',
frame='gcrs')
setattr(obj, det.name.lower() + '_pointing', det_coord)

UnboundLocalError with Trigdat._reconcile_timescales.mask

When loading a Trigdat file with Trigdat.open():

File ~/gitrepos/gifts-6u/gifts-ground-software/gdt/missions/fermi/gbm/trigdat.py:171, in Trigdat.open(cls, file_path, **kwargs)
    168 obj._rates = obj._data['RATE'].reshape(-1, 14, 8)
    170 # store the position history
--> 171 idx, dt = obj._time_indices(1024)
    172 eic = obj._data['EIC'][idx]
    173 quat = obj._data['SCATTITD'][idx]

File ~/gitrepos/gifts-6u/gifts-ground-software/gdt/missions/fermi/gbm/trigdat.py:664, in Trigdat._time_indices(self, time_res)
    662 cnt = len(idx)
    663 # reconcile 8 s and 1 s data
--> 664 idx = self._reconcile_timescales(back_idx, idx)
    666 # reconcile 8 s + 1 s and 256 ms data
    667 if time_res <= 256:

File ~/gitrepos/gifts-6u/gifts-ground-software/gdt/missions/fermi/gbm/trigdat.py:703, in Trigdat._reconcile_timescales(self, idx1, idx2)
    700 end_times2 = self._data['ENDTIME'][idx2]
    702 # find where bracketing timescale ends and inserted timescale begins
--> 703 if mask.sum():
    704     mask = end_times1 >= start_times2[0]
    705     start_idx =  (np.where(mask))[0][0]

UnboundLocalError: cannot access local variable 'mask' where it is not associated with a value

This appears to be related to the changes in #22, where mask isn't declared prior to usage:

gdt-fermi/src/gdt/missions/fermi/gbm/trigdat.py

Lines 700 to 712 in ad6c3c3

end_times2 = self._data['ENDTIME'][idx2]
# find where bracketing timescale ends and inserted timescale begins
if mask.sum():
mask = end_times1 >= start_times2[0]
start_idx = (np.where(mask))[0][0]
idx = np.concatenate((idx1[0:start_idx], idx2))
# find where inserted timescale ends and bracketing timescale begins again
mask = start_times1 >= end_times2[-1]
if mask.sum():
end_idx = (np.where(mask))[0][0]
idx = np.concatenate((idx, idx1[end_idx:]))

It looks like moving the mask initialization to before the first if mask.sum() statement fixes it, i.e.:

diff --git a/src/gdt/missions/fermi/gbm/trigdat.py b/src/gdt/missions/fermi/gbm/trigdat.py
index f23b9bc..348d08f 100644
--- a/src/gdt/missions/fermi/gbm/trigdat.py
+++ b/src/gdt/missions/fermi/gbm/trigdat.py
@@ -700,8 +700,8 @@ class Trigdat(FitsFileContextManager):
         end_times2 = self._data['ENDTIME'][idx2]
 
         # find where bracketing timescale ends and inserted timescale begins
+        mask = end_times1 >= start_times2[0]
         if mask.sum():
-            mask = end_times1 >= start_times2[0]
             start_idx =  (np.where(mask))[0][0]
             idx = np.concatenate((idx1[0:start_idx], idx2))

Size of McIlwain Plots

The current GBM BA tools produces larger McIlwain plots, which are easier to read.

In GDT, on line 67 of the script gdt/core/plot/earthplot.py, the latitudinal range is divided by 30.
In the BA tools, on line 76 of gbm/plot/earthplot.py, the latitudinal range is divided by 15.

The longitudinal and latitudinal ranges are not set as attributes of the EarthPlot class, so I suggest I change to the code. This change is not necessary, just a visibility issue.

standard_title() displays incorrect McIlwain L when plotting

Using standard_title() displays incorrect McIlwain L when plotting. Here is the code to reproduce the issue for burst number 231003905:

import matplotlib.pyplot as plt

from gdt.missions.fermi.time import Time
from gdt.missions.fermi.plot import FermiEarthPlot
from gdt.missions.fermi.gbm.finders import TriggerFtp
from gdt.missions.fermi.gbm.trigdat import Trigdat

# download trigdat file
bn = "231003905"
TriggerFtp(bn).get_trigdat(".")

trigdat = Trigdat.open("glg_trigdat_all_bn231003905_v01.fit")
trigtime = Time(trigdat.trigtime, format='fermi')

mcilwain_plot = FermiEarthPlot(mcilwain=True)
mcilwain_plot.add_spacecraft_frame(trigdat.poshist, trigtime=trigtime)
mcilwain_plot.standard_title()

plt.gcf().set_figheight(5)
plt.show()

In this case, the McIlwain L should be 1.56 but instead the title reports a value of 1.09. 1.09 is clearly inconsistent with the color code corresponding to the spacecraft position (see attached screenshot).
Screenshot 2023-11-20 at 2 17 02 PM

The source of the issue is the fact that standard_plot() does not account for East versus West when getting longitude from the coordinates of the Fermi spacecraft icon on line 344 of src/gdt/missions/fermi/plot.py:

if self.spacecraft is not None:
    coord = self.spacecraft.coordinates
    title = 'Latitude, East Longitude: ({0}, {1})\n'.format(*coord)
    mcl = calc_mcilwain_l(float(coord[0][:-1]), float(coord[1][:-1])) # this line is the issue. coord[1] corresponds to "99.70W"                                                                                                                      
                                                                      # and using coord[1][:-1] removes the "W" before float conversion

A potential fix is:

if self.spacecraft is not None:
    coord = self.spacecraft.coordinates
    title = 'Latitude, East Longitude: ({0}, {1})\n'.format(*coord)
    lat = float(coord[0][:-1])
    lon = float(coord[1][:-1]) * (-1 if "W" in coord[1] else 1)
    mcl = calc_mcilwain_l(lat, lon)

Recommend Projects

  • React photo React

    A declarative, efficient, and flexible JavaScript library for building user interfaces.

  • Vue.js photo Vue.js

    🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.

  • Typescript photo Typescript

    TypeScript is a superset of JavaScript that compiles to clean JavaScript output.

  • TensorFlow photo TensorFlow

    An Open Source Machine Learning Framework for Everyone

  • Django photo Django

    The Web framework for perfectionists with deadlines.

  • D3 photo D3

    Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

  • javascript

    JavaScript (JS) is a lightweight interpreted programming language with first-class functions.

  • web

    Some thing interesting about web. New door for the world.

  • server

    A server is a program made to process requests and deliver data to clients.

  • Machine learning

    Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.

  • Game

    Some thing interesting about game, make everyone happy.

Recommend Org

  • Facebook photo Facebook

    We are working to build community through open source technology. NB: members must have two-factor auth.

  • Microsoft photo Microsoft

    Open source projects and samples from Microsoft.

  • Google photo Google

    Google ❤️ Open Source for everyone.

  • D3 photo D3

    Data-Driven Documents codes.