Source code for sunpy.timeseries.sources.rhessi

"""
This module provides a RHESSI `~sunpy.timeseries.TimeSeries` source.
"""
import itertools
from collections import OrderedDict

import numpy as np
from pandas import DataFrame

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.metadata import MetaDict

__all__ = ['RHESSISummaryTimeSeries']


def uncompress_countrate(compressed_countrate):
    """
    Convert the compressed count rate inside of observing summary file from a
    compressed byte to a true count rate.

    Parameters
    ----------
    compressed_countrate : `byte` array
        A compressed count rate returned from an observing summary file.

    References
    ----------
    `Hsi_obs_summ_decompress.pro <https://hesperia.gsfc.nasa.gov/ssw/hessi/idl/qlook_archive/hsi_obs_summ_decompress.pro>`_
    """

    # Ensure uncompressed counts are between 0 and 255
    if (compressed_countrate.min() < 0) or (compressed_countrate.max() > 255):
        raise ValueError(
            f'Exepected uncompressed counts {compressed_countrate} to in range 0-255')

    # TODO Must be a better way than creating entire lookup table on each call
    ll = np.arange(0, 16, 1)
    lkup = np.zeros(256, dtype='int')
    _sum = 0
    for i in range(0, 16):
        lkup[16 * i:16 * (i + 1)] = ll * 2 ** i + _sum
        if i < 15:
            _sum = lkup[16 * (i + 1) - 1] + 2 ** i

    return lkup[compressed_countrate]


def parse_observing_summary_hdulist(hdulist):
    """
    Parse a RHESSI observation summary file.

    Parameters
    ----------
    hdulist : `list`
        The HDU list from the fits file.

    Returns
    -------
    out : `dict`
        Returns a dictionary.
    """
    header = hdulist[0].header

    reference_time_ut = parse_time(hdulist[5].data.field('UT_REF')[0],
                                   format='utime')
    time_interval_sec = hdulist[5].data.field('TIME_INTV')[0]
    # label_unit = fits[5].data.field('DIM1_UNIT')[0]
    # labels = fits[5].data.field('DIM1_IDS')
    labels = ['3 - 6 keV', '6 - 12 keV', '12 - 25 keV', '25 - 50 keV',
              '50 - 100 keV', '100 - 300 keV', '300 - 800 keV',
              '800 - 7000 keV', '7000 - 20000 keV']

    # The data stored in the fits file are "compressed" countrates stored as
    # one byte
    compressed_countrate = np.array(hdulist[6].data.field('countrate'))

    countrate = uncompress_countrate(compressed_countrate)
    dim = np.array(countrate[:, 0]).size

    time_array = parse_time(reference_time_ut) + \
        TimeDelta(time_interval_sec * np.arange(dim) * u.second)

    #  TODO generate the labels for the dict automatically from labels
    data = {'time': time_array, 'data': countrate, 'labels': labels}

    return header, data


[docs]class RHESSISummaryTimeSeries(GenericTimeSeries): """ RHESSI X-ray Summary lightcurve TimeSeries. The RHESSI mission consists of a single spin-stabilized spacecraft in a low-altitude orbit inclined 38 degrees to the Earth's equator. The only instrument on board is a set of 9 Germanium spectrometers with the ability to obtain high fidelity solar spectra from X rays (down to 3 keV) to gamma rays (1 MeV). Each spectrometer is coupled to a set of grids with different pitches which enable fourier-style imaging as the spacecraft spins. RHESSI provides summary lightcurves in the following passbands: * 3 - 6 keV * 6 - 12 keV * 12 - 25 keV * 25 - 50 keV * 50 - 100 keV * 100 - 300 keV * 300 - 800 keV * 800 - 7000 keV * 7000 - 20000 keV RHESSI was launched on 5th February 2002. Examples -------- >>> import sunpy.data.sample # doctest: +REMOTE_DATA >>> import sunpy.timeseries >>> rhessi = sunpy.timeseries.TimeSeries(sunpy.data.sample.RHESSI_TIMESERIES) # doctest: +REMOTE_DATA >>> rhessi.peek() # doctest: +SKIP References ---------- * `RHESSI Homepage. <https://hesperia.gsfc.nasa.gov/rhessi3/index.html>`_ * `Mission Paper. <https://doi.org/10.1023/A:1022428818870>`_ """ # Class attributes used to specify the source class of the TimeSeries # and a URL to the mission website. _source = 'rhessi' _url = "https://hesperia.gsfc.nasa.gov/rhessi3/index.html" _peek_title = "RHESSI Observing Summary Count Rate"
[docs] def plot(self, axes=None, columns=None, **kwargs): """ Plots RHESSI count rate light curve. 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) # These are a matplotlib version of the default RHESSI color cycle default_colors = ('black', 'tab:pink', 'tab:green', 'tab:cyan', 'tab:olive', 'tab:red', 'tab:blue', 'tab:orange', 'tab:brown') colors = kwargs.pop('colors', default_colors) for color, (item, frame) in zip(itertools.cycle(colors), self._data[columns].items()): axes.plot(self.to_dataframe()[columns].index, frame.values, color=color, label=item, **kwargs) axes.set_yscale("log") axes.set_ylabel('Count Rate s$^{-1}$ detector$^{-1}$') axes.legend() self._setup_x_axis(axes) return axes
@classmethod def _parse_file(cls, filepath): """ Parses rhessi FITS data files to create TimeSeries. 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 RHESSI `astropy.io.fits.HDUList` from a FITS file. Parameters ---------- hdulist : `astropy.io.fits.HDUList` A HDU list. """ header, d = parse_observing_summary_hdulist(hdulist) # The time of dict `d` is astropy.time, but dataframe can only take datetime d['time'] = d['time'].datetime header = MetaDict(OrderedDict(header)) data = DataFrame(d['data'], columns=d['labels'], index=d['time']) # Add the units data units = OrderedDict([('3 - 6 keV', u.ct / u.s / u.Unit('detector')), ('6 - 12 keV', u.ct / u.s / u.Unit('detector')), ('12 - 25 keV', u.ct / u.s / u.Unit('detector')), ('25 - 50 keV', u.ct / u.s / u.Unit('detector')), ('50 - 100 keV', u.ct / u.s / u.Unit('detector')), ('100 - 300 keV', u.ct / u.s / u.Unit('detector')), ('300 - 800 keV', u.ct / u.s / u.Unit('detector')), ('800 - 7000 keV', u.ct / u.s / u.Unit('detector')), ('7000 - 20000 keV', u.ct / u.s / u.Unit('detector'))]) # Todo: check units used. https://hesperia.gsfc.nasa.gov/ssw/hessi/doc/guides/hessi_data_access.htm return data, header, units
[docs] @classmethod def is_datasource_for(cls, **kwargs): """ Determines if the file corresponds to a RHESSI X-ray Summary `~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('telescop', '').startswith('HESSI')