2. Centering#
This part will teach you:
Basic ways to find the centroid (center of mass) position using few methods.
Important
Before proceed, read this photutils/centroids tutorial and this photutils/detection tutorial. Follow the tutorials before starting this lecture note.
# Ignore this cell if you encounter errors
%load_ext version_information
import time
now = time.strftime("%Y-%m-%d %H:%M:%S (%Z = GMT%z)")
print(f"This notebook was generated at {now} ")
vv = %version_information astropy, numpy, scipy, matplotlib, photutils, version_information
for i, pkg in enumerate(vv.packages):
print(f"{i} {pkg[0]:10s} {pkg[1]:s}")
This notebook was generated at 2023-05-09 15:36:23 (KST = GMT+0900)
0 Python 3.10.10 64bit [Clang 14.0.6 ]
1 IPython 8.13.2
2 OS macOS 13.1 arm64 arm 64bit
3 astropy 5.2.2
4 numpy 1.24.3
5 scipy 1.10.1
6 matplotlib 3.7.1
7 photutils 1.6.0
8 version_information 1.0.4
# %matplotlib notebook
from IPython.core.interactiveshell import InteractiveShell
from IPython import get_ipython
%config InlineBackend.figure_format = 'retina'
InteractiveShell.ast_node_interactivity = 'last_expr'
ipython = get_ipython()
from pathlib import Path
import numpy as np
from astropy.io import fits
from astropy.time import Time
from astropy import units as u
from astropy.nddata import CCDData, Cutout2D
from astropy.stats import sigma_clipped_stats
from matplotlib import pyplot as plt
from matplotlib import rcParams
plt.style.use('default')
rcParams.update({'font.size':12})
from photutils.centroids import centroid_com
import ysfitsutilpy as yfu
import ysphotutilpy as ypu
import warnings
warnings.filterwarnings('ignore', category=UserWarning, append=True)
import _tool_visualization as vis
DATAPATH = Path('../../Tutorial_Data')
TMPDIR = Path('tmp')
TMPDIR.mkdir(exist_ok=True)
2.1. Prepare Data#
Cut the image with
astropy.nddata.Cutout2D.I cut the image centered at
position=(270, 320)and size ofsize=(100, 100)pixels.
allfits = list(DATAPATH.glob("*p4179*.fits"))
allfits.sort()
ccd = CCDData.read(allfits[0])
cut = Cutout2D(ccd, position=(270,320), size=(95,100))
fig, axs = plt.subplots(1, 1, figsize=(3, 4), sharex=False, sharey=False, gridspec_kw=None)
vis.norm_imshow(axs, cut.data, zscale=True)
plt.tight_layout()
WARNING: FITSFixedWarning: RADECSYS= 'ICRS ' / Telescope Coordinate System
the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs]
2.2. Finding Centroid#
There are multiple practical ways to find the center of an object. Here, I will introduce four simplest ways.
The center-of-mass of the image.
Using Gaussian function.
Find the center by DAOPHOT-like algorithm
Find the center by SExtractor-like algorithm
The center of mass is easy to calculate, but in python, scipy provides a simple function center_of_mass. In astronomical images, people will just use photutils.centroids.centroid_com.
To fit Gaussian function to the image, there are several ways to do this. First, decide whether you want a circular Gaussian (4 parameters = x-center, y-center, amplitude, sigma), or an elliptical Gaussian (6 parameters = same as circular case except sigma_x, sigma_y, theta). Then you can directly fit 2-D Gaussian to the image. When background level is non-zero, you need one more constant parameter.
Due to the mathematical characteristic of Gaussian function, you can also prove this: Say \(I(x, y) = G_x(x | \mu_x, \sigma_x) \times G_y(y | \mu_y, \sigma_y)\). Here \(G_x\) and \(G_y\) are the Gaussian fitting results for the image is summed along the y- and x-direction, respectively. \(\mu\) and \(\sigma\) are the center and standard deviation of fitted Gaussians, respectively. For example, \(G_x(x| \mu_x, \sigma_y)\) is the fitted Gaussian to data \(D(x) = \sum_{y=1}{N_y} D(x, y)\). Then \(I(x, y)\) is a Gaussian with center \((\mu_x, \mu_y)\). This holds true only if \(\mathrm{min} \{I\} = 0\) is guaranteed, which is usually not the case in real images.
For several reasons, Gaussian fitting is not always the best choice. Although it’s a personal taste, but I prefer using DAOFIND or Extractor-like algorithms over Gaussian fitting.
DAOPHOT’s DAOFIND algorithm is that described in Stetson 1987 and you have already used it in photutils/detection tutorial. This is basically used when you are dealing with circular stars. If you are studying extended sources (nebula, galaxy, non-sidereal objects, etc), it may not work properly or you may undergo hard time fine-tuning the parameters.
SExtractor-like algorithm is basically used to detect general objects and also has power to calculate proper/useful center information. You will study this later in extended sources lecture note.
2.2.1. Simple Center-of-Mass#
Consider you want to find the center of mass of the star in the cropped image above. This is acutally not an easy task, because real data is contaminated by cosmic-rays, defects (dead pixels, hot pixels, …), background objects, uneven sky, etc.
To improve the robustness of center-of-mass, I here used the most elementary level algorithm I could think of.
Do sigma-clipping to the cutout data. Get statistics of the clipped data.
Set the threshold as median plus 3-sigma.
Get centroid by center of mass algorithm using only the pixels above that threshold level.
avg, med, std = sigma_clipped_stats(cut.data) # by default, 3-sigma 5-iteration.
thresh_3sig = med + 3 * std
mask_3sig = (cut.data < thresh_3sig)
center = centroid_com(data=cut.data, mask=mask_3sig)
print(center)
print(cut.to_original_position(center))
fig, axs = plt.subplots(1, 1, figsize=(3, 4), sharex=False, sharey=False, gridspec_kw=None)
vis.norm_imshow(axs, mask_3sig.astype(int))
vis.norm_imshow(axs, cut.data, alpha=0.4, zscale=True)
axs.plot(*center, 'rx')
plt.tight_layout()
[53.09601666 41.04309455]
(273.09601666127566, 314.04309454723943)
Thus, (273.10, 314.04) is the found centroid.
In the code, I overplotted the image onto the mask_3sig. The faint background-like feature is the original image, and the highlight means the pixels used for the centroid calculation. Red cross is the thus found centroid.
Practice
In realistic centroiding, this method should not be used. This is becasue virtually all of the photometric information is contained in only very small area around the true center, so we need to do centroiding using, e.g., cbox (shorthand for centroiding box) size of \(\sim\) 1-3 FWHM to minimize error sources. What you have to do is to cutout only a small region near the center of the object and do centroiding. This is not an easy task to realize, so it is better to use pre-made functions.
Consider you have to know the center of mass of this object (e.g., you cannot use Gaussian fitting, DAOPHOT, etc). How would you calculate it robustly? Give at least three possible algorithms/methods/improvements you can think of. Possible hints: You can use mask, bad-pixel rejection/interpolation, etc. You may have an idea which of them will work best only after you encounter an enormous amount of real data…
2.2.2. DAOFIND#
As you already have learnt it from the photutils tutorial, I will quickly go through the detection part:
from photutils.detection import DAOStarFinder
from photutils.aperture import CircularAperture
avg, med, std = sigma_clipped_stats(ccd.data) # by default, 3-sigma 5-iteration.
finder = DAOStarFinder(threshold=5.*std, fwhm=4, exclude_border=True)
sources = finder(ccd.data - med)
for col in sources.colnames:
sources[col].info.format = "%d" if col in ('id', 'npix') else '%.2f'
sources.pprint(max_width=76) # astropy Table
fig, axs = plt.subplots(1, 1, figsize=(8, 5), sharex=False, sharey=False, gridspec_kw=None)
vis.norm_imshow(axs, ccd.data, zscale=True)
pos = np.transpose((sources['xcentroid'], sources['ycentroid']))
aps = CircularAperture(pos, r=10.)
aps.plot(color='blue', lw=1.5, alpha=0.5)
plt.tight_layout()
plt.show();
id xcentroid ycentroid sharpness roundness1 ... sky peak flux mag
--- --------- --------- --------- ---------- ... ---- ------- ----- -----
1 175.60 15.85 0.36 0.18 ... 0.00 82.02 1.43 -0.39
2 572.53 45.07 0.39 -0.14 ... 0.00 358.09 6.22 -1.98
3 788.41 46.58 0.38 -0.32 ... 0.00 570.82 9.35 -2.43
4 150.12 62.00 0.41 -0.21 ... 0.00 424.43 7.18 -2.14
5 68.21 76.65 0.34 -0.56 ... 0.00 68.67 1.19 -0.19
6 430.40 86.52 0.44 -0.33 ... 0.00 298.27 4.92 -1.73
7 514.85 131.52 0.45 -0.10 ... 0.00 169.61 2.83 -1.13
8 406.48 152.47 0.36 -0.25 ... 0.00 255.01 4.51 -1.63
9 455.38 182.60 0.49 -0.60 ... 0.00 128.59 2.32 -0.91
10 733.55 185.89 0.40 -0.26 ... 0.00 438.38 7.39 -2.17
11 8.49 234.83 0.43 -0.35 ... 0.00 189.61 3.38 -1.32
12 852.94 245.47 0.32 -0.30 ... 0.00 259.96 4.89 -1.72
13 989.26 245.81 0.24 -0.30 ... 0.00 98.44 1.93 -0.72
14 366.73 259.55 0.50 -0.30 ... 0.00 297.76 4.56 -1.65
15 99.22 311.19 0.41 -0.23 ... 0.00 1100.25 18.58 -3.17
16 273.13 314.14 0.42 -0.21 ... 0.00 1972.81 32.66 -3.79
17 835.50 351.49 0.22 -0.46 ... 0.00 71.58 1.54 -0.47
18 683.60 352.10 0.42 -0.32 ... 0.00 2583.85 42.66 -4.08
19 979.67 417.63 0.42 -0.20 ... 0.00 4548.72 73.96 -4.67
20 398.62 428.02 0.42 -0.25 ... 0.00 233.51 3.57 -1.38
21 354.83 450.93 0.37 -0.60 ... 0.00 81.11 1.42 -0.38
22 499.84 500.13 0.36 -0.91 ... 0.00 31.97 1.00 -0.00
23 823.93 526.23 0.35 -0.17 ... 0.00 71.80 1.40 -0.37
24 988.46 539.78 0.45 -0.33 ... 0.00 376.02 6.17 -1.98
25 892.30 619.67 0.27 -0.06 ... 0.00 127.36 2.52 -1.00
26 885.56 625.93 0.38 -0.13 ... 0.00 301.92 5.34 -1.82
27 204.33 649.61 0.35 -0.27 ... 0.00 110.41 2.00 -0.75
28 913.12 652.41 0.38 -0.18 ... 0.00 358.40 6.17 -1.98
29 118.52 661.34 0.26 -0.63 ... 0.00 56.86 1.27 -0.26
30 958.73 664.12 0.42 -0.35 ... 0.00 1001.99 16.08 -3.02
31 961.49 689.83 0.59 -0.59 ... 0.00 112.69 1.68 -0.56
32 834.10 751.03 0.45 0.14 ... 0.00 186.52 3.45 -1.35
33 378.47 753.85 0.36 -0.35 ... 0.00 62.34 1.24 -0.23
34 581.31 772.72 0.36 -0.23 ... 0.00 193.50 3.23 -1.27
35 197.80 775.26 0.50 -0.05 ... 0.00 121.26 1.84 -0.66
36 700.40 775.54 0.39 -0.02 ... 0.00 521.16 8.82 -2.36
37 667.09 819.09 0.56 -0.24 ... 0.00 87.42 1.44 -0.40
38 277.67 893.89 0.42 -0.34 ... 0.00 96.85 1.75 -0.61
39 729.21 915.12 0.42 -0.34 ... 0.00 384.59 6.78 -2.08
40 813.20 930.60 0.37 -0.31 ... 0.00 1043.05 17.96 -3.14
41 561.69 949.95 0.50 -0.15 ... 0.00 101.47 1.53 -0.46
42 373.76 965.13 0.40 -0.28 ... 0.00 238.18 4.18 -1.55
Practice
Find the astropy table documentation by yourself and try other print options.
Note that DAOStarFinder has multiple parameters that is required a priori. Among them, fwhm is the FWHM of the Gaussian kernel, which must be given. Also, the bounds for roundness and sharpness related parameters have to be set properly especially when field stars are “trailed” (e.g., your telescope could not follow the stars). sigma_radius=1.5 is the factor that determines the size of the kernel. The kernel size is calculated by _StarFinderKernel.
Now if you are interested in all of these, you are good to go. To find the object that is nearest to our expected location, you may do distance calculation:
sources["dist"] = np.sqrt(np.sum((pos - np.array([270, 320]))**2, axis=1))
sources.sort("dist")
sources[0]
| id | xcentroid | ycentroid | sharpness | roundness1 | roundness2 | npix | sky | peak | flux | mag | dist |
|---|---|---|---|---|---|---|---|---|---|---|---|
| int64 | float64 | float64 | float64 | float32 | float64 | int64 | float64 | float64 | float64 | float64 | float64 |
| 16 | 273.13 | 314.14 | 0.42 | -0.21 | -0.06 | 25 | 0.00 | 1972.81 | 32.66 | -3.79 | 6.643679381248131 |
Thus, (273.13, 314.14) is the found centroid.
2.3. Accuracy of the Centroid#
In this simple case, the two methods gave very similar center location (differ by < 0.1 pix).
You have learned how to calculate the uncertainty in the centroid location from mathematics in theory class. However, in reality, the uncertainty due to the sky subtraction, instrumental defects (bad pixels or problematic flat correction, etc), cosmic-rays, etc can be larger than the centroid uncertainty. How accurate this centroid should be?
This question can be rephrased: “How accurate photometry (or astrometry) do you need?” It is difficult to give a general answer to this question because it really depends on what you are looking for. In case you only seek for 0.1 mag accuracy and aperture radius is about 3 FWHM = 6 pixels, change in centroid of 0.1 pixel order will not affect the final photometry, because (1) vast majority of flux will already have been included in 3 FWHM radius and (2) 0.1 pixel shift is much smaller than 6 pixel radius (or 12 pixel diameter). In this specific case, the partial-pixel summation algorithm would give more error than the centroid error of < 0.1 pixel.