Finding bright regions with ndimage

How you can to find the brightest regions in an AIA image and count the approximate number of regions of interest using ndimage.

# sphinx_gallery_thumbnail_number = 2

from scipy import ndimage
import matplotlib.pyplot as plt

import sunpy.map
from sunpy.data.sample import AIA_193_IMAGE

We start with the sample data

aiamap_mask = sunpy.map.Map(AIA_193_IMAGE)
aiamap = sunpy.map.Map(AIA_193_IMAGE)

First we make a mask, which tells us which regions are bright. We choose the criterion that the data should be at least 5% of the maximum value. Pixels with intensity values greater than this are included in the mask, while all other pixels are excluded.

mask = aiamap.data < aiamap.max() * 0.10

Mask is a boolean array. It can be used to modify the original map object without modifying the data. Once this mask attribute is set, we can plot the image again.

aiamap_mask.mask = mask
plt.figure()
aiamap.plot()
plt.colorbar()
plt.show()
../../../_images/sphx_glr_image_bright_regions_gallery_example_001.png

Only the brightest pixels remain in the image. However, these areas are artificially broken up into small regions. We can solve this by applying some smoothing to the image data. Here we apply a 2D Gaussian smoothing function to the data.

data2 = ndimage.gaussian_filter(aiamap.data * ~mask, 14)

The issue with the filtering is that it create pixels where the values are small (<100), so when we go on later to label this array, we get one large region which encompasses the entire array. If you want to see, just remove this line.

data2[data2 < 100] = 0

Now we will make a second SunPy map with this smoothed data.

aiamap2 = sunpy.map.Map(data2, aiamap.meta)

The function label from the scipy.ndimage module, counts the number of contiguous regions in an image.

labels, n = ndimage.label(aiamap2.data)

Finally, we plot the smoothed bright image data, along with the estimate of the number of distinct regions. We can see that approximately 6 distinct hot regions are present above the 5% of the maximum level.

plt.figure()
ax = plt.subplot(projection=aiamap)
aiamap.plot()
plt.contour(labels)
plt.figtext(0.3, 0.2, f'Number of regions = {n}', color='white')
plt.show()
../../../_images/sphx_glr_image_bright_regions_gallery_example_002.png

Total running time of the script: ( 0 minutes 1.211 seconds)

Gallery generated by Sphinx-Gallery