Plotting RHESSI Hard X-Ray emission light curve

Hi everyone. I tried plotting RHESSI Hard X-Ray emission light curves using the codes below but getting an error. I am beginner level Python programmer with no experience in RHESSI or sunpy. I’m just learning. Please guide to solve this. Thank you in advance.

Input:
import glob
import matplotlib
import matplotlib.pyplot as plt
import numpy as np

from sunpy import timeseries as ts
from sunpy.net import Fido
from sunpy.net import attrs as a
from sunpy.net.dataretriever.sources.rhessi import RHESSIClient as RC

Data File Directory

outdir=‘/Users/A/Desktop/Filename’

Search for the data with Fido

trange=a.Time(“2010-11-05 12:00:26”,“2010-11-05 17:00:00”)

result = Fido.search(trange, a.Instrument.rhessi, a.Physobs.summary_lightcurve)

Get the file from Fido - can obviously skip this step if already downloaded (Not Working)

fg15 = Fido.fetch(result,path=outdir)

If already downloaded can skip above step and just load back in (Data files are detected correctly and displayed as output since I downloaded the data files manually)

fg15=glob.glob(outdir+‘hsi20101105*.fits’)

Good idea to always cancatenate incase more than one file

rhessi = ts.TimeSeries(fg15,concatenate=True)

What did timeseries load in

rhessi.meta

Quick peek at what timeseries loaded in

rhessi.peek()

Output:

IndexError Traceback (most recent call last)
Cell In[6], line 6
1 # My Random Comments
2 # My Random Comments
3 # My Random Comments
4
5 # Good idea to always cancatenate in case more than one file
----> 6 rhessi = ts.TimeSeries(fg15,concatenate=True)
8 # What did timeseries load in
9 rhessi.meta

File ~\anaconda3\Lib\site-packages\sunpy\timeseries\timeseries_factory.py:430, in TimeSeriesFactory.call(self, silence_errors, *args, **kwargs)
410 “”"
411 Method for running the factory. Takes arbitrary arguments and keyword
412 arguments and passes them to a sequence of pre-registered types to
(…)
427 Extra keyword arguments are passed through to sunpy.io.read_file such as memmap for FITS files.
428 “”"
429 self.silence_errors = silence_errors
→ 430 new_timeseries = self._parse_args(*args, **kwargs)
432 # Concatenate the timeseries into one if specified.
433 concatenate = kwargs.get(‘concatenate’, False)

File ~\anaconda3\Lib\site-packages\sunpy\timeseries\timeseries_factory.py:328, in TimeSeriesFactory._parse_args(self, *args, **kwargs)
326 for arg in args:
327 try:
→ 328 all_ts += self._parse_arg(arg, **kwargs)
329 except (NoMatchError, MultipleMatchError, ValidationFunctionError):
330 if self.silence_errors:

File ~\anaconda3\Lib\site-packages\sunpy\util\functools.py:18, in seconddispatch..wrapper(*args, **kwargs)
17 def wrapper(*args, **kwargs):
—> 18 return dispatcher.dispatch(args[1].class)(*args, **kwargs)

File ~\anaconda3\Lib\site-packages\sunpy\timeseries\timeseries_factory.py:394, in TimeSeriesFactory._parse_path(self, path, **kwargs)
390 raise MultipleMatchError(“Multiple HDUs return multiple matching classes.”)
392 cls = types[0]
→ 394 data_header_unit_tuple = cls._parse_hdus(pairs)
395 all_ts += self._parse_arg(data_header_unit_tuple)
397 return all_ts

File ~\anaconda3\Lib\site-packages\sunpy\timeseries\sources\rhessi.py:233, in RHESSISummaryTimeSeries._parse_hdus(cls, hdulist)
223 @classmethod
224 def _parse_hdus(cls, hdulist):
225 “”"
226 Parses a RHESSI astropy.io.fits.HDUList from a FITS file.
227
(…)
231 A HDU list.
232 “”"
→ 233 header, d = parse_observing_summary_hdulist(hdulist)
234 # The time of dict d is astropy.time, but dataframe can only take datetime
235 d[‘time’] = d[‘time’].datetime

File ~\anaconda3\Lib\site-packages\sunpy\timeseries\sources\rhessi.py:73, in parse_observing_summary_hdulist(hdulist)
58 “”"
59 Parse a RHESSI observation summary file.
60
(…)
69 Returns a dictionary.
70 “”"
71 header = hdulist[0].header
—> 73 reference_time_ut = parse_time(hdulist[5].data.field(‘UT_REF’)[0],
74 format=‘utime’)
75 time_interval_sec = hdulist[5].data.field(‘TIME_INTV’)[0]
76 # label_unit = fits[5].data.field(‘DIM1_UNIT’)[0]
77 # labels = fits[5].data.field(‘DIM1_IDS’)

IndexError: list index out of range

Hello @TheNovusCoder,

I need a bit more information about your setup. What version of astropy and sunpy do you have installed?

Can you also see if your code works if you remove the fg15=glob.glob(outdir+‘hsi20101105*.fits’) line? I can get the download to work and plot but I bypassed the glob search. I suspect that might be pulling in a file that the Timeseries class can not handle. If you can provide the output from that glob search as well, so we can check.

For future reference you can use backticks ` to put all your code in a block that will be neatly formatted.

Hi @nabobalis . Thanks for your quick reply. I appreciate it very much.

I am using Jupyter Notebook. I am using SunPy 3.1.0 and Astopy 6.1.0. I am using the old versions since it is better for me to run some codes in this without error. I hope this is not an issue.

For downloading I used this:
fg15 = Fido.fetch(result,path=outdir)

And got this error:

--------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[8], line 2
      1 #Get the file from Fido - can obviously skip this step if already downloaded (Not Working)
----> 2 fg15 = Fido.fetch(result,path=outdir)

File ~\anaconda3\Lib\site-packages\sunpy\net\fido_factory.py:384, in UnifiedDownloaderFactory.fetch(self, path, max_conn, progress, overwrite, downloader, *query_results, **kwargs)
    381     raise ValueError("wait is not a valid keyword argument to Fido.fetch.")
    383 if downloader is None:
--> 384     downloader = Downloader(max_conn=max_conn, progress=progress, overwrite=overwrite)
    385 elif not isinstance(downloader, parfive.Downloader):
    386     raise TypeError("The downloader argument must be a parfive.Downloader object.")

File ~\anaconda3\Lib\site-packages\sunpy\util\parfive_helpers.py:30, in Downloader.__init__(self, *args, **kwargs)
     27 if headers:
     28     kwargs['headers'] = headers
---> 30 super().__init__(*args, **kwargs)

TypeError: Downloader.__init__() got an unexpected keyword argument 'headers'

I downloaded the fits file data from here (Index of /hessidata) for 5th November 2010 manually. Got the link from Fido.search function. Since, I can’t get my code to download the data, I downloaded them manually into my pc , then from there I just loaded them. I tried using the os and path modules separately to do so, now without glob.glob function as per your instruction. When I tried to plot I’m still getting the same error.

Output:

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
Cell In[39], line 2
      1 # Good idea to always cancatenate incase more than one file
----> 2 rhessi = ts.TimeSeries(fg15,concatenate=True)
      4 # What did timeseries load in
      5 rhessi.meta

File ~\anaconda3\Lib\site-packages\sunpy\timeseries\timeseries_factory.py:430, in TimeSeriesFactory.__call__(self, silence_errors, *args, **kwargs)
    410 """
    411 Method for running the factory. Takes arbitrary arguments and keyword
    412 arguments and passes them to a sequence of pre-registered types to
   (...)
    427 Extra keyword arguments are passed through to `sunpy.io.read_file` such as `memmap` for FITS files.
    428 """
    429 self.silence_errors = silence_errors
--> 430 new_timeseries = self._parse_args(*args, **kwargs)
    432 # Concatenate the timeseries into one if specified.
    433 concatenate = kwargs.get('concatenate', False)

File ~\anaconda3\Lib\site-packages\sunpy\timeseries\timeseries_factory.py:328, in TimeSeriesFactory._parse_args(self, *args, **kwargs)
    326 for arg in args:
    327     try:
--> 328         all_ts += self._parse_arg(arg, **kwargs)
    329     except (NoMatchError, MultipleMatchError, ValidationFunctionError):
    330         if self.silence_errors:

File ~\anaconda3\Lib\site-packages\sunpy\util\functools.py:18, in seconddispatch.<locals>.wrapper(*args, **kwargs)
     17 def wrapper(*args, **kwargs):
---> 18     return dispatcher.dispatch(args[1].__class__)(*args, **kwargs)

File ~\anaconda3\Lib\site-packages\sunpy\timeseries\timeseries_factory.py:394, in TimeSeriesFactory._parse_path(self, path, **kwargs)
    390             raise MultipleMatchError("Multiple HDUs return multiple matching classes.")
    392         cls = types[0]
--> 394         data_header_unit_tuple = cls._parse_hdus(pairs)
    395         all_ts += self._parse_arg(data_header_unit_tuple)
    397 return all_ts

File ~\anaconda3\Lib\site-packages\sunpy\timeseries\sources\rhessi.py:233, in RHESSISummaryTimeSeries._parse_hdus(cls, hdulist)
    223 @classmethod
    224 def _parse_hdus(cls, hdulist):
    225     """
    226     Parses a RHESSI `astropy.io.fits.HDUList` from a FITS file.
    227 
   (...)
    231         A HDU list.
    232     """
--> 233     header, d = parse_observing_summary_hdulist(hdulist)
    234     # The time of dict `d` is astropy.time, but dataframe can only take datetime
    235     d['time'] = d['time'].datetime

File ~\anaconda3\Lib\site-packages\sunpy\timeseries\sources\rhessi.py:73, in parse_observing_summary_hdulist(hdulist)
     58 """
     59 Parse a RHESSI observation summary file.
     60 
   (...)
     69     Returns a dictionary.
     70 """
     71 header = hdulist[0].header
---> 73 reference_time_ut = parse_time(hdulist[5].data.field('UT_REF')[0],
     74                                format='utime')
     75 time_interval_sec = hdulist[5].data.field('TIME_INTV')[0]
     76 # label_unit = fits[5].data.field('DIM1_UNIT')[0]
     77 # labels = fits[5].data.field('DIM1_IDS')

IndexError: list index out of range

I am also attaching the glob output for your reference.
The output from the glob search:
['/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_003400_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_005400_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_012720_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_023440_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_030320_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_040900_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_043900_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_055720_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_061440_001.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_073700_001.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_075040_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_092500_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_092620_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_110200_003.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_122240_003.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_123800_005.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_132240_005.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_135520_005.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_141340_004.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_154920_003.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_172520_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_190100_003.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_203700_001.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_221240_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_234700_002.fits', '/Users/A/Desktop/RHESSI_DATA/Goddard\\hsi_20101105_234820_002.fits']

The downloading function is not working separately as well just to plot. May I see your code? I am not sure how you made it to work.

I am using Jupyter Notebook. I am using SunPy 3.1.0 and Astopy 6.1.0. I am using the old versions since it is better for me to run some codes in this without error. I hope this is not an issue.

The problem with this is that these old versions of sunpy are not meant to be used where you have to update versions of astropy or other packages. They are just not compatible. You have installed the latest version of astropy for example here. That is not an old version.

For example:

TypeError: Downloader.init() got an unexpected keyword argument ‘headers’

Is due the older version of sunpy and a not compatible version of parfive, which handles the downloads underneath.

I think you cannot fix this without updating sunpy or by downgrading every other software package.

Regarding your file loading problem. From what I can tell, one of the files is not a valid file that is able to be read by sunpy. If you loop over each file like:

for file in glob.glob(outdir+"hsi20101105*.fits"):
   print(file)
    ts.TimeSeries(file)

and see if all of the error or none of them error.
It is also possible that the concatenate is also failing instead of loading the file.

I will say that if there is a bug here, we can’t patch an old version of sunpy, the best we can do is patch the release version which is 6.0. So you would have to update if we have discovered a bug and fixed it.

The downloading function is not working separately as well just to plot. May I see your code? I am not sure how you made it to work.

My code was the following:

from sunpy import timeseries as ts
from sunpy.net import Fido
from sunpy.net import attrs as a

trange=a.Time("2010-11-05 12:00:26","2010-11-05 17:00:00")

result = Fido.search(trange, a.Instrument.rhessi, a.Physobs.summary_lightcurve)
fg15 = Fido.fetch(result)
rhessi = ts.TimeSeries(fg15)
rhessi.peek()

Here’s where things went awry for you. You have “incorrectly” manually downloaded RHESSI Level-0 FITS files (which contain individual photon events), but the file that needs to be downloaded is the RHESSI observing-summary FITS file (which contain processed lightcurves). Here’s the file for the day you want:

https://hesperia.gsfc.nasa.gov/hessidata/metadata/catalog/hsi_obssumm_20101105_051.fits

2 Likes

Thank you @nabobalis, I managed to plot the graph nicely. Thank you @ayshih , I noticed it. I am using sunpy 6.0.1 and astropy 6.1.2 now after updating.

But I noticed the graph is weird, most probably with background noise. My SunPy plot looks different from the one in the RHESSI Quicklook browser. I just focused on 2 energy levels (12-25keV and 25-50keV) but the plots are smooth in the sample (12-25keV and 25-50keV) without the extremely protruding line on top. I could not find any resources in SunPy documentation to remove the background noise from the data. I need help in this part.

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

from sunpy import timeseries as ts
from sunpy.net import Fido
from sunpy.net import attrs as a

outdir='/Users/A/Desktop/RHESSI_CORRRECT_DATA/'

trange=a.Time("2010-11-05 12:00:26","2010-11-05 17:00:00")

result = Fido.search(trange, a.Instrument.rhessi, a.Physobs.summary_lightcurve)

fg15 = Fido.fetch(result,path=outdir)
rhessi = ts.TimeSeries(fg15)
rhessi.peek()

trhessi=rhessi.truncate(trange.start.iso,trange.end.iso)
print(trhessi)

#Plot using pandas and matplotlib to overplot some interesting time ranges
trhessi_df=trhessi.to_dataframe()
trhessi_tims=trhessi_df.index
trhessi_weak=trhessi_df['12 - 25 keV']
trhessi_strong=trhessi_df['25 - 50 keV']

# Can also overplot some interesting times
inttime=a.Time("2010-11-05 13:08:00","2010-11-05 13:09:30")

# And now doing a more manual of this truncated data
# Some extra lines to make it look nicer, particularly in the time labelling

plt.rcParams.update({'font.size': 18,'font.family':"sans-serif",\
                         'font.sans-serif':"Arial",'mathtext.default':"regular"})

fig,ax = plt.subplots(figsize=(10, 8))
plt.plot(trhessi_tims,trhessi_weak,drawstyle='steps-post',marker=None,color='red',lw=2,label='$12 - 25 keV$')
plt.plot(trhessi_tims,trhessi_strong,drawstyle='steps-post',marker=None,color='aqua',lw=2,label='$25 - 50 keV$')


ax.set_ylabel("RHESSI Count Rates [$s^{-1}$ detector $^{-1}$]")
ax.set_xlabel("Start Time "+trange.start.iso[:-4])
ax.set_yscale("log")
ax.set_ylim([1e0,1e5])
ax.set_xlim([trange.start.datetime,trange.end.datetime])
# precisely control the x time labels
myFmt = matplotlib.dates.DateFormatter('%H:%M')
majorx= matplotlib.dates.MinuteLocator(interval=60)
minorx= matplotlib.dates.MinuteLocator(interval=20)
ax.set_title('Plot of \nRHESSI Hard X-Ray Emission \nagainst \nStart Time 2010-11-05 12:00:26')
ax.xaxis.set_major_locator(majorx)
ax.xaxis.set_minor_locator(minorx)
ax.xaxis.set_major_formatter(myFmt)


ax.axvspan(inttime.start.datetime,inttime.end.datetime,color='black',alpha=0.2, label='interesting time range')

plt.legend()
plt.show()

I am adding the sample RHESSI QUicklook browser plot here. (Please focus on 12-25keV and 25-50keV)
sample rhessi quicklook browser

What you are seeing is the difference between the normal quicklook rates (your plot using sunpy) and the “corrected” quicklook rates (the image you are comparing against). If you choose “uncorrected” instead of “corrected” in the RHESSI Browser, you will see this plot:
hsi_20101105_123800_rate
Note that it matches what you see in your plot using sunpy.

This is not “background noise”, but rather that RHESSI has attenuators that are moved in and out to moderate the incident flux. These attenuators have the greatest effect at the lowest energies, and you can see the different attenuator states indicated as magenta bars at the top labeled “A0”, “A1”, and “A3” . The normal or “uncorrected” quicklook lightcurves are merely showing the attenuators doing their job.

It is not possible to accurately correct the rates for the changing attenuator states without a spectral model, which formally requires a careful analysis of the flare spectrum and the underlying background. Instead, for quicklook rates, the “corrected” quicklook rates are calculated using a heuristic algorithm that essentially scales each lightcurve segment so that the overall result looks continuous. However, this correction approach is intended only for qualitative use because it is not scientifically robust: the algorithm can get misled if there is a real transition/burst close to an attenuator-state change, resulting in short artifacts or even long-term scaling errors.

sunpy does not currently generate the “corrected” quicklook rates from the data in the FITS file, so there’s no easy way for you to reproduce the plot you want. You are welcome to submit a feature request to GitHub!

1 Like

Thank you @ayshih . Now, I understand it clearly. Sure, I will request for the feature. It would be helpful to have it in SunPy as well.