Audience: Advanced beginner / intermediate (requires basic knowledge of scipp.DataArray)

Objectives: Constructing binned data from an event list, modifying binning in order to extract different quantities, and basic masking and filtering of event and binned data.

RHESSI Solar Flares#

This tutorial covers the basics of binned data in Scipp by analyzing the list of solar flares recorded by NASA’s RHESSI small explorer [Lin et al].

The input data has been constructed from the official flare list. It is available as a HDF5 file in Scipp’s own format and can be downloaded and accessed directly via scipp.data as shown below.

Attention

The tutorial data has been filtered and modified. It should not be used for any actual scientific analyses!

See docs/tutorials/prepare_data_rhessi.py in the scipp source for details.

2 Spatial Distribution#

2.1 Bin by x and y#

Plot the spatial distribution of flares.

Plotting the event list flares would yield a scatter plot which is not particularly useful. A better approach is computing and plotting the density as a function of x and y. This is commonly done by histogramming the events. But Scipp offers an alternative: binning.

Scipp’s ‘binned data’ is similar to a histogram, except that the individual events are preserved. They are simply collected into bins as defined by bin-edge coordinates.

Define bin-edges for x and y, use scipp.bin to create binned data from flares, and plot the result.

Tip:

Use scipp.arange or scipp.linspace to construct the edges. Make sure to use the correct unit!
Use sc.bin(flares, edges=[<edge_x>, <edge_y>]).
Turn your binned data into a histogram before plotting using <binned>.bins.sum().
See Plotting 2D data for the relevant parts of the plotting API.

[9]:

# YOUR CODE

Solution

Bin spatially.

spatial = flares.bin(y=sc.linspace('y', -1200, 1200, 100, unit='asec'),
                     x=sc.linspace('x', -1200, 1200, 100, unit='asec'))

Histogram and plot. The arguments to plot are optional but improve the result here.

spatial.hist().plot(aspect='equal',
                    norm='log',
                    labels={
                        'x': 'x [asec]',
                        'y': 'y [asec]'
                    },
                    cmap='inferno')

2.2 Remove Outliers#

The plot shows a lot of outliers. They are not the non-solar events from before because those are not visible in the plot as they were masked out. Instead, the outliers are caused by RHESSI’s electronics or analysis software glitching out and assigning bad positions to the flares. This includes the circular shape even though it looks deceptively Sun-like.

The flare list does not include enough information to exclude all such bad positions. But the instrument can only detect x-rays for \(x \in [-1000~\text{asec}, 1000~\text{asec}]\) and \(y \in [-600~\text{asec}, 600~\text{asec}]\). So to first order, all events outside that range should be removed.

Instead of using a mask as before, use binning this time. Create new bin-edges for x and y with the proper limits and bin the data with them, this will remove all events outside the valid range. Plot the result. (You can either bin the original flares data array or apply new bins to the previously binned array; both via scipp.bin).

[10]:

# YOUR CODE

Solution

Bin into a more narrow range. Bin sizes are chosen such that bins are square.

spatial = flares.bin(y=sc.linspace('y', -600, 600, 90, unit='asec'),
                     x=sc.linspace('x', -1000, 1000, 150, unit='asec'))

And plot like before.

spatial.hist().plot(aspect='equal',
                    norm='log',
                    labels={
                        'x': 'x [asec]',
                        'y': 'y [asec]'
                    },
                    cmap='inferno')

2.3 Correct for Detector Efficiency#

In this tutorial, we assume that the instrument consists of a 3x3 grid of detectors which each record x-rays from distinct directions. (The reality is more complicated, of course. See the wiki for more information.) Furthermore, the tutorial data has been manipulated to simulate different efficiencies of the individual detectors.

The efficiency is available as an attribute:

[11]:

efficiency = flares.attrs['detector_efficiency'].value
efficiency

[11]:

scipp.DataArray (1.38 KB)

Dimensions:
- x: 3
- y: 3
Coordinates: (2)
- x
  (x [bin-edge])
  float64
  4.84813681109535984e-06rad
  -1000.000, -300.0, 300.0, 1000.000
```
Values:
array([-1000.,  -300.,   300.,  1000.])
```
- y
  (y [bin-edge])
  float64
  4.84813681109535984e-06rad
  -600.0, -200.0, 200.0, 600.0
```
Values:
array([-600., -200.,  200.,  600.])
```

Data: (1)

(x, y)

float64

𝟙

0.9, 0.95, ..., 0.94, 0.91

Values:
array([[0.9 , 0.95, 0.77],
       [0.98, 0.89, 0.82],
       [0.93, 0.94, 0.91]])

Normalize the data by dividing by efficiency. You first need to bin into the correct bins as defined by the coordinates of efficiency.

[12]:

# YOUR CODE

Solution

coarse_spatial = spatial.bin(x=efficiency.coords['x'],
                             y=efficiency.coords['y'])
corrected = coarse_spatial / efficiency

You can plot the corrected data as before. But it makes sense to return to smaller bins in order to resolve the distribution properly. This can be done using scipp.bin with finer edges because the data was only binned and not histogrammed.

[13]:

# YOUR CODE

Solution

Bin in the same way as before but using corrected as input

spatial = corrected.bin(y=sc.linspace('y', -600, 600, 90, unit='asec'),
                        x=sc.linspace('x', -1000, 1000, 150, unit='asec'))

And plot like before.

spatial.hist().plot(aspect='equal',
                    norm='log',
                    labels={
                        'x': 'x [asec]',
                        'y': 'y [asec]'
                    },
                    cmap='inferno')

3 Time Series#

Now, we want to look at the distribution of flares over time. The temporal distribution can be obtained like the spatial one using binning.

Select one of the time coordinates (e.g. ‘peak_time’) and bin with an appropriate bin size. Plot the result.

Important: Use the pre-binned data from the previous task in order to include the detector normalization.

Tip:

The time coordinates are event-coordinates, that is, they are not defined at the top level of the binned data (i.e. per bin) but inside of the bins (i.e. per event). You can access them using <binned>.bins.coords['peak_time'].
You can add binning by time using <binned>.bin(peak_time=<n_bins>) and produce a three-dimensional array. But here, we are more interested in a one-dimensional distribution. So we need to erase the binning in x and y which can be achieved, e.g., using <binned>.bins.concat('x').

Alternatively, the lower level function scipp.binning.make_binned can be used to bin in time and erase spatial binning at the same time.

[14]:

# YOUR CODE

Solution

Bin in time and keep the existing spatial binning. transpose is used to put the time dimension on the slider in the plot.

temporal_and_spatial = spatial.bin(peak_time=200)
to_plot = temporal_and_spatial.hist()
to_plot.transpose(['peak_time', 'y', 'x']).plot(vmin=0.0 * sc.units.one,
                                                vmax=0.66 * sc.units.one,
                                                labels={
                                                    'x': 'x [asec]',
                                                    'y': 'y [asec]'
                                                },
                                                cmap='inferno')

Or remove the existing binning, yielding a 1D-variable.

temporal = spatial.bins.concat('x').bins.concat('y').bin(peak_time=200)
temporal.hist().plot(ylabel='number of flares')

Note that it is important to first remove the binning in x and y by using concat. Binning in time first and then removing other binning would use too much memory. It is possible to do everything in one step using a lower-level API:

from scipp.binning import make_binned

time = spatial.bins.coords['peak_time']
min_time = time.min().value
max_time = time.max().value
step = (max_time - min_time) / 200
time_edges = sc.arange('peak_time',
                       min_time,
                       max_time,
                       step,
                       unit=time.bins.unit)
temporal = make_binned(spatial, edges=[time_edges], erase=('x', 'y'))
temporal.hist().plot(ylabel='number of flares')

3.1 Flare Durations#

Another interesting quantity to look at is the duration of flares. The duration was computed earlier and should already be stored as an event coordinate.

Plot the duration as a function of time.

Tip:

Construct a new data array from the previous result: duration = <binned_by_time>.copy(). And assign new data using duration.bins.data = <duration_data>.
Previously, we used .hist() (which is an alias for .bins.sum()) to make histograms because the data was given as counts. Now the data is seconds which should be averaged instead of summed. So use .bins.mean() instead.

[15]:

# YOUR CODE

Solution

Create a new data array containing the duration as its data.

duration = temporal.copy()
duration.name = 'duration'
duration.bins.data = duration.bins.coords.pop('duration')

Plot duration vs time.

duration.bins.mean().plot()

4 Energy Bands#

The flares were recorded in several non-overlapping energy bands. Those are identified by the ‘min_energy’ and ‘max_energy’ attributes in flares. Since the bands do not overlap, it is sufficient to label them with ‘min_energy’ for simplicity.

In this section, we want to split the temporal distribution from above into the different energy bands. Group the temporal distribution by ‘min_energy’ to obtain two-dimensional data. Use scipp.group instead of scipp.bin this time because every event has exactly one of a set of possible energies instead of a value in a range.

[16]:

# YOUR CODE

Solution

Skipping the lowest energy band here because those events are not confirmed to be solar flares. This is optional.

grouped_by_energy = temporal.group('min_energy')['min_energy', 1:]

Plot the result. A 2D plot is not very useful here, so split the data by ‘min_energy’ and either plot each energy in a separate plot or combine them into a dictionary and plot that:

lines = {'<name>': <data_array>}
sc.plot(lines)

[17]:

# YOUR CODE

Solution

lines = {
    f"min_energy={grouped_by_energy['min_energy', i].attrs['min_energy'].value}":
    grouped_by_energy['min_energy', i].bins.sum()
    for i in range(grouped_by_energy.sizes['min_energy'])
}
sc.plot(lines, ylabel='number of flares')

References#

Lin, R., Dennis, B., Hurford, G. et al. The Reuven Ramaty High-Energy Solar Spectroscopic Imager (RHESSI). Sol Phys 210, 3–32 (2002). doi:10.1023/A:1022428818870

RHESSI Solar Flares

Contents

RHESSI Solar Flares#

1 Loading and preprocessing data#

1.1 Load Flare List#

1.2 Inspect the Loaded Data#

1.3 Compute Duration#

1.4 Create Masks#

2 Spatial Distribution#

2.1 Bin by x and y#

2.2 Remove Outliers#

2.3 Correct for Detector Efficiency#

3 Time Series#

3.1 Flare Durations#

4 Energy Bands#

References#