Aperture Photometry (photutils.aperture)

Introduction

In Photutils, the aperture_photometry() function is the main tool to perform aperture photometry on an astronomical image for a given set of apertures.

Photutils provides several apertures defined in pixel or sky coordinates. The aperture classes that are defined in pixel coordinates are:

Each of these classes has a corresponding variant defined in celestial coordinates:

To perform aperture photometry with sky-based apertures, one will need to specify a WCS transformation.

Users can also create their own custom apertures (see Defining Your Own Custom Apertures).

Creating Aperture Objects

The first step in performing aperture photometry is to create an aperture object. An aperture object is defined by a position (or a list of positions) and parameters that define its size and possibly, orientation (e.g., an elliptical aperture).

We start with an example of creating a circular aperture in pixel coordinates using the CircularAperture class:

>>> from photutils import CircularAperture
>>> positions = [(30., 30.), (40., 40.)]
>>> apertures = CircularAperture(positions, r=3.)

The positions should be either a single tuple of (x, y), a list of (x, y) tuples, or an array with shape Nx2, where N is the number of positions. The above example defines two circular apertures located at pixel coordinates (30, 30) and (40, 40) with a radius of 3 pixels.

Creating an aperture object in celestial coordinates is similar. One first uses the SkyCoord class to define celestial coordinates and then the SkyCircularAperture class to define the aperture object:

>>> from astropy import units as u
>>> from astropy.coordinates import SkyCoord
>>> from photutils import SkyCircularAperture
>>> positions = SkyCoord(l=[1.2, 2.3] * u.deg, b=[0.1, 0.2] * u.deg,
...                      frame='galactic')
>>> apertures = SkyCircularAperture(positions, r=4. * u.arcsec)

Note

Sky apertures are not defined completely in celestial coordinates. They simply use celestial coordinates to define the central position, and the remaining parameters are converted to pixels using the pixel scale of the image at the central position. Projection distortions are not taken into account. If the apertures were defined completely in celestial coordinates, their shapes would not be preserved when converting to pixel coordinates.

Performing Aperture Photometry

After the aperture object is created, we can then perform the photometry using the aperture_photometry() function. We start by defining the apertures as described above:

>>> positions = [(30., 30.), (40., 40.)]
>>> apertures = CircularAperture(positions, r=3.)

and then we call the aperture_photometry() function with the data and the apertures:

>>> import numpy as np
>>> from photutils import aperture_photometry
>>> data = np.ones((100, 100))
>>> phot_table = aperture_photometry(data, apertures)
>>> print(phot_table)    
 id xcenter ycenter  aperture_sum
      pix     pix
--- ------- ------- -------------
  1    30.0    30.0 28.2743338823
  2    40.0    40.0 28.2743338823

This function returns the results of the photometry in an Astropy QTable. In this example, the table has four columns, named 'id', 'xcenter', 'ycenter', and 'aperture_sum'.

Since all the data values are 1.0, the aperture sums are equal to the area of a circle with a radius of 3:

>>> print(np.pi * 3. ** 2)    
28.2743338823

Aperture and Pixel Overlap

The overlap of the apertures with the data pixels can be handled in different ways. For the default method (method='exact'), the exact intersection of the aperture with each pixel is calculated. The other options, 'center' and 'subpixel', are faster, but with the expense of less precision. For 'center', a pixel is considered to be entirely in or out of the aperture depending on whether its center is in or out of the aperture. For 'subpixel', pixels are divided into a number of subpixels, which are in or out of the aperture based on their centers. For this method, the number of subpixels needs to be set with the subpixels keyword.

This example uses the 'subpixel' method where pixels are resampled by a factor of 5 (subpixels=5) in each dimension:

>>> phot_table = aperture_photometry(data, apertures, method='subpixel',
...                                  subpixels=5)
>>> print(phot_table)    
 id xcenter ycenter aperture_sum
      pix     pix
--- ------- ------- ------------
  1    30.0    30.0        27.96
  2    40.0    40.0        27.96

Note that the results differ from the true value of 28.274333 (see above).

For the 'subpixel' method, the default value is subpixels=5, meaning that each pixel is equally divided into 25 smaller pixels (this is the method and subsampling factor used in SourceExtractor.). The precision can be increased by increasing subpixels, but note that computation time will be increased.

Multiple Apertures at Each Position

While the Aperture objects support multiple positions, they must have a fixed shape, e.g. radius, size, and orientation.

To perform photometry in multiple apertures at each position, one may input a list of aperture objects to the aperture_photometry() function.

Suppose that we wish to use three circular apertures, with radii of 3, 4, and 5 pixels, on each source:

>>> radii = [3., 4., 5.]
>>> apertures = [CircularAperture(positions, r=r) for r in radii]
>>> phot_table = aperture_photometry(data, apertures)
>>> print(phot_table)    
 id xcenter ycenter aperture_sum_0 aperture_sum_1 aperture_sum_2
      pix     pix
--- ------- ------- -------------- -------------- --------------
  1    30.0    30.0  28.2743338823  50.2654824574  78.5398163397
  2    40.0    40.0  28.2743338823  50.2654824574  78.5398163397

For multiple apertures, the output table column names are appended with the positions index.

Other apertures have multiple parameters specifying the aperture size and orientation. For example, for elliptical apertures, one must specify a, b, and theta:

>>> from photutils import EllipticalAperture
>>> a = 5.
>>> b = 3.
>>> theta = np.pi / 4.
>>> apertures = EllipticalAperture(positions, a, b, theta)
>>> phot_table = aperture_photometry(data, apertures)
>>> print(phot_table)    
 id xcenter ycenter  aperture_sum
      pix     pix
--- ------- ------- -------------
  1    30.0    30.0 47.1238898038
  2    40.0    40.0 47.1238898038

Again, for multiple apertures one should input a list of aperture objects, each with identical positions:

>>> a = [5., 6., 7.]
>>> b = [3., 4., 5.]
>>> theta = np.pi / 4.
>>> apertures = [EllipticalAperture(positions, a=ai, b=bi, theta=theta)
...              for (ai, bi) in zip(a, b)]
>>> phot_table = aperture_photometry(data, apertures)
>>> print(phot_table)    
 id xcenter ycenter aperture_sum_0 aperture_sum_1 aperture_sum_2
      pix     pix
--- ------- ------- -------------- -------------- --------------
  1    30.0    30.0  47.1238898038  75.3982236862  109.955742876
  2    40.0    40.0  47.1238898038  75.3982236862  109.955742876

Background Subtraction

Global Background Subtraction

aperture_photometry() assumes that the data have been background-subtracted. If bkg is an array representing the background of the data (determined by Background2D or an external function), simply do:

>>> phot_table = aperture_photometry(data - bkg, apertures)  

Local Background Subtraction

Suppose we want to perform the photometry in a circular aperture with a radius of 3 pixels and estimate the local background level around each source with a circular annulus of inner radius 6 pixels and outer radius 8 pixels. We start by defining the apertures:

>>> from photutils import CircularAnnulus
>>> apertures = CircularAperture(positions, r=3)
>>> annulus_apertures = CircularAnnulus(positions, r_in=6., r_out=8.)

We then perform the photometry in both apertures:

>>> apers = [apertures, annulus_apertures]
>>> phot_table = aperture_photometry(data, apers)
>>> print(phot_table)    
 id xcenter ycenter aperture_sum_0 aperture_sum_1
      pix     pix
--- ------- ------- -------------- --------------
  1    30.0    30.0  28.2743338823  87.9645943005
  2    40.0    40.0  28.2743338823  87.9645943005

Note that we cannot simply subtract the aperture sums because the apertures have different areas.

To calculate the mean local background within the circular annulus aperture, we need to divide its sum by its area, which can be calculated using the area() method:

>>> bkg_mean = phot_table['aperture_sum_1'] / annulus_apertures.area()

The background sum within the circular aperture is then the mean local background times the circular aperture area:

>>> bkg_sum = bkg_mean * apertures.area()
>>> final_sum = phot_table['aperture_sum_0'] - bkg_sum
>>> phot_table['residual_aperture_sum'] = final_sum
>>> print(phot_table['residual_aperture_sum'])    
residual_aperture_sum
---------------------
    -7.1054273576e-15
    -7.1054273576e-15

The result here should be zero because all of the data values are 1.0 (the tiny difference from 0.0 is due to numerical precision).

Error Estimation

If and only if the error keyword is input to aperture_photometry(), the returned table will include a 'aperture_sum_err' column in addition to 'aperture_sum'. 'aperture_sum_err' provides the propagated uncertainty associated with 'aperture_sum'.

For example, suppose we have previously calculated the error on each pixel’s value and saved it in the array error:

>>> error = 0.1 * data
>>> phot_table = aperture_photometry(data, apertures, error=error)
>>> print(phot_table)    
 id xcenter ycenter  aperture_sum aperture_sum_err
      pix     pix
--- ------- ------- ------------- ----------------
  1    30.0    30.0 28.2743338823   0.531736155272
  2    40.0    40.0 28.2743338823   0.531736155272

'aperture_sum_err' values are given by:

\[\Delta F = \sqrt{\sum_{i \in A} \sigma_{\mathrm{tot}, i}^2}\]

where \(\Delta F\) is source_sum_err, \(A\) are the non-masked pixels in the aperture, and \(\sigma_{\mathrm{tot}, i}\) is the input error array.

In the example above, it is assumed that the error keyword specifies the total error – either it includes Poisson noise due to individual sources or such noise is irrelevant. However, it is often the case that one has calculated a smooth “background-only error” array, which by design doesn’t include increased noise on bright pixels. To include Poisson noise from the sources, we can use the calc_total_error() function.

Let’s assume we have a background-only image called bkg_error. If our data are in units of electrons/s, we would use the exposure time as the effective gain:

>>> from photutils.utils import calc_total_error
>>> effective_gain = 500   # seconds
>>> error = calc_total_error(data, bkg_error, effective_gain)    
>>> phot_table = aperture_photometry(data - bkg, apertures, error=error)    

Note

In cases where the error array is slowly varying across the image, it is not necessary to sum the error from every pixel in the aperture individually. Instead, we can approximate the error as being roughly constant across the aperture and simply take the value of \(\sigma\) at the center of the aperture. This can be done by setting the keyword pixelwise_errors=False. In this case the flux error is

\[\Delta F = \sigma \sqrt{A}\]

where \(\sigma\) is the error at the center of the aperture and \(A\) is the area of the aperture.

Pixel Masking

Pixels can be ignored/excluded (e.g., bad pixels) from the aperture photometry by providing an image mask via the mask keyword:

>>> data = np.ones((5, 5))
>>> aperture = CircularAperture((2, 2), 2.)
>>> mask = np.zeros_like(data, dtype=bool)
>>> data[2, 2] = 100.   # bad pixel
>>> mask[2, 2] = True
>>> t1 = aperture_photometry(data, aperture, mask=mask)
>>> print(t1['aperture_sum'])    
 aperture_sum
 -------------
 11.5663706144

The result is very different if a mask image is not provided:

>>> t2 = aperture_photometry(data, aperture)
>>> print(t2['aperture_sum'])    
aperture_sum
-------------
111.566370614

Aperture Photometry Using Sky Coordinates

As mentioned in Creating Aperture Objects, performing photometry using apertures defined in celestial coordinates simply requires defining a “sky” aperture at positions defined by a SkyCoord object. Here we show an example of photometry on real data using a SkyCircularAperture.

We start by loading a Spitzer 4.5 micron image of a region of the Galactic plane:

>>> from photutils import datasets
>>> hdu = datasets.load_spitzer_image()   
Downloading http://data.astropy.org/photometry/spitzer_example_image.fits [Done]
>>> catalog = datasets.load_spitzer_catalog()   
Downloading http://data.astropy.org/photometry/spitzer_example_catalog.xml [Done]

The catalog contains (among other things) the Galactic coordinates of the sources in the image as well as the PSF-fitted fluxes from the official Spitzer data reduction. We define the apertures positions based on the existing catalog positions:

>>> positions = SkyCoord(catalog['l'], catalog['b'], frame='galactic')   
>>> apertures = SkyCircularAperture(positions, r=4.8 * u.arcsec)   

Now perform the photometry in these apertures using the hdu. The hdu object is a FITS HDU that contains the data and a header describing the WCS transformation of the image. The WCS includes the coordinate frame of the image and the projection from celestial to pixel coordinates. The aperture_photometry function uses the WCS information to automatically convert the apertures defined in celestial coordinates into pixel coordinates:

>>> phot_table = aperture_photometry(hdu, apertures)    

The Spitzer catalog also contains the official fluxes for the sources, so we can compare to our fluxes. Because the Spitzer catalog fluxes are in units of mJy and the data are in units of MJy/sr, we need to convert units before comparing the results. The image data have a pixel scale of 1.2 arcsec/pixel.

>>> import astropy.units as u
>>> factor = (1.2 * u.arcsec) ** 2 / u.pixel
>>> fluxes_catalog = catalog['f4_5']   
>>> converted_aperture_sum = (phot_table['aperture_sum'] *
...                           factor).to(u.mJy / u.pixel)   

Finally, we can plot the comparison of the photometry:

>>> import matplotlib.pyplot as plt
>>> plt.scatter(fluxes_catalog, converted_aperture_sum.value)
>>> plt.xlabel('Spitzer catalog PSF-fit fluxes ')
>>> plt.ylabel('Aperture photometry fluxes')

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

../_images/aperture-1.png

Despite using different methods, the two catalogs are in good agreement. The aperture photometry fluxes are based on a circular aperture with a radius of 4.8 arcsec. The Spitzer catalog fluxes were computed using PSF photometry. Therefore, differences are expected between the two measurements.

Aperture Masks

All PixelAperture objects have a to_mask() method that returns a list of ApertureMask objects, one for each aperture position. The ApertureMask object contains a cutout of the aperture mask and a slices object that provides the locations where the mask is to be applied. It also provides a to_image() method to obtain an image of the mask in a 2D array of the given shape, a cutout() method to create a cutout from the input data over the mask bounding box, and an apply() method to apply the aperture mask to the input data to create a mask-weighted data cutout. All of these methods properly handle the cases of partial or no overlap of the aperture mask with the data.

Let’s start by creating an aperture object:

>>> from photutils import CircularAperture
>>> positions = [(30., 30.), (40., 40.)]
>>> apertures = CircularAperture(positions, r=3.)

Now let’s create a list of ApertureMask objects using the to_mask() method:

>>> masks = aperture.to_mask(method='center')

We can now create an image with of the first aperture mask at its position:

>>> mask = masks[0]
>>> image = mask.to_image(shape=((200, 200)))

We can also create a cutout from a data image over the mask domain:

>>> data_cutout = mask.cutout(data)

We can also create a mask-weighted cutout from the data. Here the circular aperture mask has been applied to the data:

>>> data_cutout_aper = mask.apply(data)

Defining Your Own Custom Apertures

The aperture_photometry() function can perform aperture photometry in arbitrary apertures. This function accepts any Aperture-derived objects, such as CircularAperture. This makes it simple to extend functionality: a new type of aperture photometry simply requires the definition of a new Aperture subclass.

All PixelAperture subclasses must define a _slices property, to_mask() and plot() methods, and optionally an area() method. All SkyAperture subclasses must implement only a to_pixel() method.

  • _slices: A property defining the minimal bounding box slices for the aperture at each position.
  • to_mask(): A method to return a list of ApertureMask objects, one for each aperture position.
  • area(): A method to return the exact analytical area (in pixels**2) of the aperture.
  • plot(): A method to plot the aperture on a matplotlib.axes.Axes instance.

See Also

  1. IRAF’s APPHOT specification [PDF] (Sec. 3.3.5.8 - 3.3.5.9)
  2. SourceExtractor Manual [PDF] (Sec. 9.4 p. 36)

Reference/API

This subpackage contains modules and packages for identifying sources in an astronomical image.

Functions

aperture_photometry(data, apertures[, ...]) Perform aperture photometry on the input data by summing the flux within the given aperture(s).

Classes

Aperture Abstract base class for all apertures.
ApertureMask(mask, bbox_slice) Class for an aperture mask.
CircularAnnulus Circular annulus aperture(s), defined in pixel coordinates.
CircularAperture Circular aperture(s), defined in pixel coordinates.
CircularMaskMixin Mixin class to create masks for circular and circular-annulus aperture objects.
EllipticalAnnulus Elliptical annulus aperture(s), defined in pixel coordinates.
EllipticalAperture Elliptical aperture(s), defined in pixel coordinates.
EllipticalMaskMixin Mixin class to create masks for elliptical and elliptical-annulus aperture objects.
PixelAperture Abstract base class for 2D apertures defined in pixel coordinates.
RectangularAnnulus Rectangular annulus aperture(s), defined in pixel coordinates.
RectangularAperture Rectangular aperture(s), defined in pixel coordinates.
RectangularMaskMixin Mixin class to create masks for rectangular or rectangular-annulus aperture objects.
SkyAperture Abstract base class for 2D apertures defined in celestial coordinates.
SkyCircularAnnulus Circular annulus aperture(s), defined in sky coordinates.
SkyCircularAperture Circular aperture(s), defined in sky coordinates.
SkyEllipticalAnnulus Elliptical annulus aperture(s), defined in sky coordinates.
SkyEllipticalAperture Elliptical aperture(s), defined in sky coordinates.
SkyRectangularAnnulus Rectangular annulus aperture(s), defined in sky coordinates.
SkyRectangularAperture Rectangular aperture(s), defined in sky coordinates.

Class Inheritance Diagram

Inheritance diagram of photutils.aperture.core.Aperture, photutils.aperture.core.ApertureMask, photutils.aperture.circle.CircularAnnulus, photutils.aperture.circle.CircularAperture, photutils.aperture.circle.CircularMaskMixin, photutils.aperture.ellipse.EllipticalAnnulus, photutils.aperture.ellipse.EllipticalAperture, photutils.aperture.ellipse.EllipticalMaskMixin, photutils.aperture.core.PixelAperture, photutils.aperture.rectangle.RectangularAnnulus, photutils.aperture.rectangle.RectangularAperture, photutils.aperture.rectangle.RectangularMaskMixin, photutils.aperture.core.SkyAperture, photutils.aperture.circle.SkyCircularAnnulus, photutils.aperture.circle.SkyCircularAperture, photutils.aperture.ellipse.SkyEllipticalAnnulus, photutils.aperture.ellipse.SkyEllipticalAperture, photutils.aperture.rectangle.SkyRectangularAnnulus, photutils.aperture.rectangle.SkyRectangularAperture