Getting Started with Photutils¶

The following example uses Photutils to find sources in an astronomical image and then perform circular aperture photometry on them.

We start by loading an image from the bundled datasets and selecting a subset of the image:

>>> import numpy as np
>>> from photutils.datasets import load_star_image
>>> hdu = load_star_image()  
>>> image = hdu.data[500:700, 500:700].astype(float)  

We then subtract a rough estimate of the background, calculated using the image median:

>>> image -= np.median(image)  

In the remainder of this example, we assume that the data is background-subtracted.

Photutils supports several source detection algorithms. For this example, we use DAOStarFinder to detect the stars in the image. We set the detection threshold at the 3-sigma noise level, estimated using the median absolute deviation (mad_std) of the image. The parameters of the detected sources are returned as an Astropy Table:

>>> from photutils.detection import DAOStarFinder
>>> from astropy.stats import mad_std
>>> bkg_sigma = mad_std(image)  
>>> daofind = DAOStarFinder(fwhm=4.0, threshold=3.0 * bkg_sigma)  
>>> sources = daofind(image)  
>>> for col in sources.colnames:  
...     sources[col].info.format = '%.8g'  # for consistent table output
>>> print(sources)  
 id xcentroid ycentroid  sharpness  ... sky peak    flux       mag
--- --------- ---------- ---------- ... --- ---- --------- -----------
  1 182.83866 0.16767019 0.85099873 ...   0 3824 2.8028346  -1.1189937
  2 189.20431 0.26081353  0.7400477 ...   0 4913 3.8729185  -1.4700959
  3 5.7946491  2.6125424 0.39589731 ...   0 7752 4.1029107  -1.5327302
  4 36.847063  1.3220228 0.29594528 ...   0 8739 7.4315818  -2.1777032
  5 3.2565602   5.418952 0.35985495 ...   0 6935 3.8126298  -1.4530616
...       ...        ...        ... ... ...  ...       ...         ...
147 197.24864  186.16647 0.31211532 ...   0 8302 7.5814629  -2.1993825
148 124.31327  188.30523  0.5362742 ...   0 6702 6.6358543  -2.0547421
149 24.257207  194.71494 0.44169546 ...   0 8342 3.2671037  -1.2854073
150    116.45  195.05923 0.67080547 ...   0 3299 2.8775221  -1.1475467
151 18.958086  196.34207 0.56502139 ...   0 3854 2.3835296 -0.94305138
152 111.52575  195.73192 0.45827852 ...   0 8109 7.9278607    -2.24789
Length = 152 rows

Using the source locations (i.e., the xcentroid and ycentroid columns), we now define circular apertures centered at these positions with a radius of 4 pixels and compute the sum of the pixel values within the apertures. The aperture_photometry() function returns an Astropy QTable with the results of the photometry:

>>> from photutils.aperture import aperture_photometry, CircularAperture
>>> positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
>>> apertures = CircularAperture(positions, r=4.0)  
>>> phot_table = aperture_photometry(image, apertures)  
>>> for col in phot_table.colnames:  
...     phot_table[col].info.format = '%.8g'  # for consistent table output
>>> print(phot_table)  
 id  xcenter   ycenter   aperture_sum
       pix       pix
--- --------- ---------- ------------
  1 182.83866 0.16767019    18121.759
  2 189.20431 0.26081353    29836.515
  3 5.7946491  2.6125424    331979.82
  4 36.847063  1.3220228    183705.09
  5 3.2565602   5.418952    349468.98
...       ...        ...          ...
148 124.31327  188.30523    45084.874
149 24.257207  194.71494    355778.01
150    116.45  195.05923    31232.912
151 18.958086  196.34207    162076.26
152 111.52575  195.73192    82795.715
Length = 152 rows

The sum of the pixel values within the apertures are given in the aperture_sum column.

Finally, we plot the image and the defined apertures:

>>> import matplotlib.pyplot as plt
>>> plt.imshow(image, cmap='gray_r', origin='lower')
>>> apertures.plot(color='blue', lw=1.5, alpha=0.5)

(Source code, png, hires.png, pdf, svg)

Navigation

Getting Started with Photutils¶