8.3.3. Peak finding in a stellar image#
8.3.3.1. Data used in this notebook#
The image in this notebook is of the field of the star TIC 125489084 on a night when the moon was nearly full. The moonlight caused a smooth gradient across the background of the image.
As discussed in the section about background subtraction with this image, the best way to remove the background in this case is to use photutils to construct a 2D model of the background.
As we will see in a moment, DAOStarFinder does a good job identifying most of the sources in the image, since all of the sources are stars.
import warnings
from astropy.nddata import CCDData
from astropy.stats import sigma_clipped_stats, gaussian_sigma_to_fwhm, SigmaClip
from astropy.table import QTable
from astropy.visualization import simple_norm, SqrtStretch
from astropy.visualization.mpl_normalize import AsymmetricPercentileInterval, ImageNormalize
import matplotlib.pyplot as plt
import numpy as np
from photutils.aperture import CircularAperture, EllipticalAperture
from photutils.background import Background2D, MeanBackground
from photutils.centroids import centroid_2dg
from photutils.datasets import make_100gaussians_image, make_gaussian_sources_image
from photutils.detection import find_peaks, DAOStarFinder
plt.style.use('../photutils_notebook_style.mplstyle')
To begin, let’s create the image and subtract the background. Recall that this image contains only stars.
ccd_image = CCDData.read('TIC_125489084.01-S001-R055-C001-ip.fit.bz2')
sigma_clip = SigmaClip(sigma=3., maxiters=5)
bkg_estimator = MeanBackground()
bkg = Background2D(ccd_image, box_size=200, filter_size=(9, 9), mask=ccd_image.mask,
sigma_clip=sigma_clip, bkg_estimator=bkg_estimator)
# Calculate the 2D background subtraction, maintaining metadata, unit, and mask
ccd_2d_bkgdsub = ccd_image.subtract(bkg.background)
data_subtracted = ccd_2d_bkgdsub.data
WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / Equatorial coordinate system
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
mean, med, std = sigma_clipped_stats(data_subtracted, sigma=3.0, maxiters=5)
8.3.3.2. Source Detection with peak finding#
As in the previous two peak finding examples, we include a centroid function so that we can compare the sources detected by peak finding with the sources identified by the star-finding algorithms.
The centroid-finding, which we do here only so that we can compare positions found by this method to positions found by the star-finding methods, can generate many warnings when an attempt is made to fit a 2D Gaussian to the sources. The warnings are suppressed below to avoid cluttering the output.
with warnings.catch_warnings(action="ignore"):
sources_findpeaks = find_peaks(data_subtracted,
threshold=5. * std, box_size=29,
centroid_func=centroid_2dg)
print(sources_findpeaks)
id x_peak y_peak peak_value x_centroid y_centroid
--- ------ ------ ------------------ ------------------ ------------------
1 594 9 75.48432203308093 589.0952756642187 8.118290769552006
2 0 11 1198758.7748060143 8.77598977676826 12.176195213387953
3 351 11 246.65855271165825 353.70687383993516 7.443245887904401
4 2437 19 145.94909280988685 2436.400192873056 19.06115030572174
5 114 25 75.18232329013745 112.02970821066799 28.48959187723852
6 820 58 287.01603957515306 820.2761073463196 57.617455777410086
7 2950 83 136.11548539236028 2938.301322736696 82.35932871516538
8 3582 91 13919.254779455989 3581.7150437134924 91.02785627934037
9 0 97 407317.19684704806 7.252263374091754 101.53312061481788
10 1066 103 398.52729216961643 1068.65832037505 124.94316408227596
... ... ... ... ... ...
543 674 4023 240.61123785548745 672.8391203918304 4024.0354373236714
544 0 4029 176856.29746557667 8.696875311378543 4029.007264941508
545 905 4038 77.64652344995051 912.1633119304082 3971.38481059316
546 2093 4049 170.61674009204864 2091.7194224694895 4050.0541157416196
547 1640 4056 198.92298484050536 1640.011594256438 4055.5841322456877
548 0 4060 1242276.4335966578 7.61684134636393 4069.146720573554
549 2214 4071 78.92570788169812 2210.933634915861 4074.336715526592
550 623 4073 155.93380587902922 622.7529605095823 4074.5111686743962
551 0 4079 401068.7122964189 4.013596319533449 4084.810795390543
552 3994 4080 74.31587924706832 3990.967572055417 4082.849626388253
Length = 552 rows
The roughly 550 sources identified by find_peaks is much larger than the number detected by the star-finding approach.
Let’s take a look at the detections on the image to see what is going on.
# Set up the figure with subplots
fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))
# Plot the 2D-subtracted data
norm = ImageNormalize(ccd_2d_bkgdsub, interval=AsymmetricPercentileInterval(30, 99.5))
fitsplot = ax1.imshow(data_subtracted, cmap='viridis', origin='lower', norm=norm,
interpolation='nearest')
ax1.scatter(sources_findpeaks['x_peak'], sources_findpeaks['y_peak'], s=30, marker='o',
lw=1, alpha=0.7, facecolor='None', edgecolor='r')
ax1.set_title("find_peaks, no masking")
Text(0.5, 1.0, 'find_peaks, no masking')
8.3.3.2.1. Mask bad regions of the CCCD#
The column on the left side of this CCD is not actually exposed to light, so find_peaks finds a bunch of spurious sources. Recall that we can provide a mask to the source finder. As we did in the star-finding approach, we exclude the region within 50 pixels of the edge.
# Create a mask array
mask = np.zeros_like(ccd_2d_bkgdsub, dtype=bool)
# Set the mask for a border of 50 pixels around the image
mask[:50, :] = mask[-50:, :] = mask[:, :50] = mask[:, -50:] = 1.0
with warnings.catch_warnings(action="ignore"):
sources_findpeaks = find_peaks(data_subtracted, mask=mask,
threshold=5. * std, box_size=29,
centroid_func=centroid_2dg)
# Set up the figure with subplots
fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))
# Plot the 2D-subtracted data
norm = ImageNormalize(ccd_2d_bkgdsub, interval=AsymmetricPercentileInterval(30, 99.5))
fitsplot = ax1.imshow(data_subtracted, cmap='viridis', origin='lower', norm=norm,
interpolation='nearest')
ax1.scatter(sources_findpeaks['x_peak'], sources_findpeaks['y_peak'], s=30, marker='o',
lw=1, alpha=0.7, facecolor='None', edgecolor='r')
ax1.set_title("find_peaks, image edges masked");
This looks quite a bit more reasonable.
8.3.3.3. Comparing Detection Methods#
Let’s compare how peak finding compares to using DAOStarFinder on this image.
# Create the star finder
fwhm = 6.2
daofinder = DAOStarFinder(threshold=5 * std, fwhm=2 * fwhm)
sources_dao = daofinder(ccd_2d_bkgdsub.data, mask=mask)
print(f'''DAOStarFinder: {len(sources_dao)} sources
find_peaks: {len(sources_findpeaks)} sources''')
DAOStarFinder: 442 sources
find_peaks: 389 sources
In this example, many more sources are detected by DAOStarFinder than by find_peaks. A plot showing both detection methods illustrates that the star-finding approach is doing a better job here.
fig, ax1 = plt.subplots(1, 1, figsize=(8, 8))
# Plot the data
fitsplot = ax1.imshow(data_subtracted, norm=norm, cmap='Greys', origin='lower')
marker_size = 60
ax1.scatter(sources_findpeaks['x_peak'], sources_findpeaks['y_peak'], s=marker_size, marker='s',
lw=1, alpha=1, facecolor='None', edgecolor='r', label='Found by find_peaks')
ax1.scatter(sources_dao['xcentroid'], sources_dao['ycentroid'], s=2 * marker_size, marker='D',
lw=1, alpha=1, facecolor='None', edgecolor='#0077BB', label='Found by DAOfind')
# Add legend
ax1.legend(ncol=2, loc='lower center', bbox_to_anchor=(0.5, -0.20))
ax1.set_xlabel('X (pixels)')
ax1.set_ylabel('Y (pixels)')
ax1.set_title('Sources Found by Different Methods');
8.3.3.4. Summary#
In this case, DAOStarFinder performs better than find_peaks. This makes sense, since DAOStarFinder is tuned for finding stars, and all of the sources in this image are stars.