Note
Go to the end to download the full example code.
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.
import matplotlib.pyplot as plt
from scipy import ndimage
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 10% 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 bool
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
fig = plt.figure()
ax = fig.add_subplot(projection=aiamap_mask)
aiamap_mask.plot(axes=ax)
plt.colorbar()
plt.show()

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.
Now we will make a second sunpy map with this smoothed data.
aiamap2 = sunpy.map.Map(data2, aiamap.meta)
The function scipy.ndimage.label
counts the number of contiguous regions
in an image.
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 10% of the maximum level.
fig = plt.figure()
ax = fig.add_subplot(projection=aiamap)
aiamap.plot(axes=ax)
ax.contour(labels)
plt.figtext(0.3, 0.2, f'Number of regions = {n}', color='white')
plt.show()

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