Source code for sunpy.timeseries.sources.fermi_gbm

"""
This module FERMI GBM `~sunpy.timeseries.TimeSeries` source.
"""
from collections import OrderedDict

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import astropy.units as u
from astropy.time import TimeDelta

import sunpy.io
from sunpy.time import parse_time
from sunpy.timeseries.timeseriesbase import GenericTimeSeries
from sunpy.util.decorators import deprecate_positional_args_since
from sunpy.util.metadata import MetaDict
from sunpy.visualization import peek_show

__all__ = ['GBMSummaryTimeSeries']


[docs] class GBMSummaryTimeSeries(GenericTimeSeries): """ Fermi/GBM Summary lightcurve TimeSeries. The Gamma-ray Burst Monitor (GBM) is an instrument on board Fermi. It is meant to detect gamma-ray bursts but also detects solar flares. It consists of 12 Sodium Iodide (NaI) scintillation detectors and 2 Bismuth Germanate (BGO) scintillation detectors. The NaI detectors cover from a few keV to about 1 MeV and provide burst triggers and locations. The BGO detectors cover the energy range from about 150 keV to about 30 MeV. This summary lightcurve makes use of the CSPEC (daily version) data set which consists of the counts accumulated every 4.096 seconds in 128 energy channels for each of the 14 detectors. Note that the data is re-binned from the original 128 into the following 8 pre-determined energy channels. The rebinning method treats the counts in each of the original 128 channels as all having the energy of the average energy of that channel. For example, the counts in an 14.5--15.6 keV original channel would all be accumulated into the 15--25 keV rebinned channel. * 4-15 keV * 15-25 keV * 25-50 keV * 50-100 keV * 100-300 keV * 300-800 keV * 800-2000 keV Examples -------- >>> import sunpy.timeseries >>> import sunpy.data.sample # doctest: +REMOTE_DATA >>> gbm = sunpy.timeseries.TimeSeries(sunpy.data.sample.GBM_TIMESERIES, source='GBMSummary') # doctest: +REMOTE_DATA >>> gbm.peek() # doctest: +SKIP References ---------- * `Fermi Mission Homepage <https://fermi.gsfc.nasa.gov>`_ * `Fermi GBM Homepage <https://fermi.gsfc.nasa.gov/science/instruments/gbm.html>`_ * `Fermi Science Support Center <https://fermi.gsfc.nasa.gov/ssc/>`_ * `Fermi Data Product <https://fermi.gsfc.nasa.gov/ssc/data/access/>`_ * `GBM Instrument Papers <https://gammaray.nsstc.nasa.gov/gbm/publications/instrument_journal_gbm.html>`_ """ # Class attributes used to specify the source class of the TimeSeries # and a URL to the mission website. _source = 'gbmsummary' _url = "https://gammaray.nsstc.nasa.gov/gbm/#"
[docs] def plot(self, axes=None, columns=None, **kwargs): """ Plots the GBM timeseries. Parameters ---------- axes : `matplotlib.axes.Axes`, optional The axes on which to plot the TimeSeries. Defaults to current axes. columns : list[str], optional If provided, only plot the specified columns. **kwargs : `dict` Additional plot keyword arguments that are handed to `~matplotlib.axes.Axes.plot` functions. Returns ------- `~matplotlib.axes.Axes` The plot axes. """ axes, columns = self._setup_axes_columns(axes, columns) for d in columns: axes.plot(self._data.index, self._data[d], label=d, **kwargs) axes.set_yscale("log") axes.set_ylabel('Counts/s/keV') axes.legend() self._setup_x_axis(axes) return axes
[docs] @peek_show @deprecate_positional_args_since("4.1") def peek(self, *, title=None, columns=None, **kwargs): """ Displays the GBM timeseries by calling `~sunpy.timeseries.sources.fermi_gbm.GBMSummaryTimeSeries.plot`. .. plot:: import sunpy.timeseries import sunpy.data.sample gbm = sunpy.timeseries.TimeSeries(sunpy.data.sample.GBM_TIMESERIES, source='GBMSummary') gbm.peek() Parameters ---------- title : `str`, optional The title of the plot. columns : list[str], optional If provided, only plot the specified columns. **kwargs : `dict` Additional plot keyword arguments that are handed to `~matplotlib.axes.Axes.plot` functions. """ if title is None: title = 'Fermi GBM Summary data ' + str(self.meta.get('DETNAM').values()) fig, ax = plt.subplots() axes = self.plot(axes=ax, columns=columns, **kwargs) axes.set_title(title) return fig
@classmethod def _parse_file(cls, filepath): """ Parses a GBM CSPEC FITS file. Parameters ---------- filepath : `str` The path to the file you want to parse. """ hdus = sunpy.io.read_file(filepath) return cls._parse_hdus(hdus) @classmethod def _parse_hdus(cls, hdulist): """ Parses a GBM CSPEC `astropy.io.fits.HDUList`. Parameters ---------- hdulist : `str` The path to the file you want to parse. """ header = MetaDict(OrderedDict(hdulist[0].header)) # these GBM files have three FITS extensions. # extn1 - this gives the energy range for each of the 128 energy bins # extn2 - this contains the data, e.g. counts, exposure time, time of observation # extn3 - eclipse times? energy_bins = hdulist[1].data count_data = hdulist[2].data # rebin the 128 energy channels into some summary ranges # 4-15 keV, 15 - 25 keV, 25-50 keV, 50-100 keV, 100-300 keV, 300-800 keV, 800 - 2000 keV # put the data in the units of counts/s/keV summary_counts = _bin_data_for_summary(energy_bins, count_data) # get the time information in datetime format with the correct MET adjustment met_ref_time = parse_time('2001-01-01 00:00') # Mission elapsed time gbm_times = met_ref_time + TimeDelta(count_data['time'], format='sec') gbm_times.precision = 9 gbm_times = gbm_times.isot.astype('datetime64') column_labels = ['4-15 keV', '15-25 keV', '25-50 keV', '50-100 keV', '100-300 keV', '300-800 keV', '800-2000 keV'] # Add the units data units = OrderedDict([('4-15 keV', u.ct / u.s / u.keV), ('15-25 keV', u.ct / u.s / u.keV), ('25-50 keV', u.ct / u.s / u.keV), ('50-100 keV', u.ct / u.s / u.keV), ('100-300 keV', u.ct / u.s / u.keV), ('300-800 keV', u.ct / u.s / u.keV), ('800-2000 keV', u.ct / u.s / u.keV)]) return pd.DataFrame(summary_counts, columns=column_labels, index=gbm_times), header, units
[docs] @classmethod def is_datasource_for(cls, **kwargs): """ Determines if the file corresponds to a GBM summary lightcurve `~sunpy.timeseries.TimeSeries`. """ # Check if source is explicitly assigned if 'source' in kwargs.keys(): if kwargs.get('source', ''): return kwargs.get('source', '').lower().startswith(cls._source) # Check if HDU defines the source instrument if 'meta' in kwargs.keys(): return kwargs['meta'].get('INSTRUME', '').startswith('GBM')
def _bin_data_for_summary(energy_bins, count_data): """ Rebin the 128 energy channels into some summary ranges and put the data in the units of counts/s/keV. Bin ranges used: * 4-15 keV * 15-25 keV * 25-50 keV * 50-100 keV * 100-300 keV * 300-800 keV * 800-2000 keV Parameters ---------- energy_bins : `numpy.ndarray` The array of energy bins to rebin. count_data : `numpy.ndarray` The array of count data to rebin. """ # list of energy bands to sum between ebands = [4, 15, 25, 50, 100, 300, 800, 2000] e_center = (energy_bins['e_min'] + energy_bins['e_max']) / 2 indices = [np.searchsorted(e_center, e) for e in ebands] summary_counts = [] for ind_start, ind_end in zip(indices[:-1], indices[1:]): # sum the counts in the energy bands, and find counts/s/keV summed_counts = np.sum(count_data["counts"][:, ind_start:ind_end], axis=1) energy_width = (energy_bins["e_max"][ind_end - 1] - energy_bins["e_min"][ind_start]) summary_counts.append(summed_counts/energy_width/count_data["exposure"]) return np.array(summary_counts).T