-
Notifications
You must be signed in to change notification settings - Fork 1
/
stara.py
79 lines (67 loc) · 2.82 KB
/
stara.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
import astropy.units as u
import numpy as np
from skimage.measure import regionprops, label, regionprops_table
from skimage.filters import median
from skimage.morphology import disk, square, white_tophat
from skimage.util import invert
import sunpy.map
from astropy.table import QTable
from astropy.time import Time
@u.quantity_input
def stara(
smap,
circle_radius: u.deg = 100 * u.arcsec,
median_box: u.deg = 10 * u.arcsec,
threshold=6000,
limb_filter: u.percent = None,
):
"""
A method for automatically detecting sunspots in white-light data using morphological operations
Parameters
----------
smap : `sunpy.map.GenericMap`
The map to apply the algorithm to.
circle_radius : `astropy.units.Quantity`, optional
The angular size of the structuring element used in the
`skimage.morphology.white_tophat`. This is the maximum radius of
detected features.
median_box : `astropy.units.Quantity`, optional
The size of the structuring element for the median filter, features
smaller than this will be averaged out.
threshold : `int`, optional
The threshold used for detection, this will be subject to detector
degradation. The default is a reasonable value for HMI continuum images.
limb_filter : `astropy.units.Quantity`, optional
If set, ignore features close to the limb within a percentage of the
radius of the disk. A value of 10% generally filters out false
detections around the limb with HMI continuum images.
"""
data = invert(smap.data)
# Filter things that are close to limb to reduce false detections
if limb_filter is not None:
hpc_coords = sunpy.map.all_coordinates_from_map(smap)
r = np.sqrt(hpc_coords.Tx ** 2 + hpc_coords.Ty ** 2) / (
smap.rsun_obs - smap.rsun_obs * limb_filter
)
data[r > 1] = np.nan
# Median filter to remove detections based on hot pixels
m_pix = int((median_box / smap.scale[0]).to_value(u.pix))
med = median(data, square(m_pix), behavior="ndimage")
# Construct the pixel structuring element
c_pix = int((circle_radius / smap.scale[0]).to_value(u.pix))
circle = disk(c_pix / 2)
finite = white_tophat(med, circle)
finite[np.isnan(finite)] = 0 # Filter out nans
return finite > threshold
def get_regions(segmentation, smap, properties=("label", "centroid", "area", "min_intensity")):
labelled = label(segmentation)
if labelled.max() == 0:
return QTable()
regions = regionprops_table(
labelled, smap.data, properties=properties
)
regions["obstime"] = Time([smap.date] * regions["label"].size)
regions["center_coord"] = smap.pixel_to_world(
regions["centroid-0"] * u.pix, regions["centroid-1"] * u.pix
).heliographic_stonyhurst
return QTable(regions)