"""
Tools for aperture photometry with non native bg/error methods
This function serves to ease the computation of photometric magnitudes
and errors using PhotUtils by replicating DAOPHOT's photometry and
error methods. The formula for DAOPHOT's error is:
err = sqrt (Poisson_noise / epadu + area * stdev**2 + area**2 * stdev**2 / nsky)
Which gives a magnitude error:
mag_err = 1.0857 * err / flux
Where epadu is electrons per ADU (gain), area is the photometric
aperture area, stdev is the uncertainty in the sky measurement
and nsky is the sky annulus area. To get the uncertainty in the sky
we must use a custom background tool, which also enables computation of
the mean and median of the sky as well (more robust statistics).
All the stats are sigma clipped. These are calculated by the
functions in aperture_stats_tbl.
.. note::
Currently, the background computations will fully include a pixel that has ANY overlap with the
background aperture (the annulus). This is to simplify the computation of the median, as a weighted
median is nontrivial, and slower. Copied from https://grit.stsci.edu/HLA/software/blob/master/HLApipeline/HLApipe/scripts/photometry_tools.py
Authors
-------
- Varun Bajaj, January 2018
Use
---
::
from photometry_tools import iraf_style_photometry
phot_aps = CircularAperture((sources['xcentroid'], sources['ycentroid']),r=10.)
bg_aps = CircularAnnulus((sources['xcentroid'], sources['ycentroid']), r_in=13., r_out=15.)
Simplest call:
::
photometry_tbl = iraf_style_photometry(phot_aps, bg_aps, data)
Pass in pixelwise error and set background to mean
::
photometry_tbl = iraf_style_photometry(phot_aps, bg_aps, data, error_array=data_err, bg_method='mean')
Can also set the gain (if image units are DN)
::
photometry_tbl = iraf_style_photometry(phot_aps, bg_aps, data, epadu=2.5)
Classes and Functions
---------------------
"""
import numpy as np
import math
from astropy.table import Table
from drizzlepac.haputils import constants
from drizzlepac.haputils.background_median import aperture_stats_tbl
from photutils.aperture import aperture_photometry
[docs]
def iraf_style_photometry(phot_apertures, bg_apertures, data, photflam, photplam, error_array=None,
bg_method='mode', epadu=1.0):
"""
Computes photometry with PhotUtils apertures, with IRAF formulae
Parameters
----------
phot_apertures : photutils PixelAperture object (or subclass)
The PhotUtils apertures object to compute the photometry. i.e. the object returned via CirularAperture.
bg_apertures : photutils PixelAperture object (or subclass)
The phoutils aperture object to measure the background in. i.e. the object returned via CircularAnnulus.
data : array
The data for the image to be measured.
photflam : float
inverse sensitivity, in ergs/cm2/angstrom/electron
photplam : float
Pivot wavelength, in angstroms
error_array : array
(Optional) The array of pixelwise error of the data. If none, the Poisson noise term in the error computation
will just be the square root of the flux/epadu. If not none, the aperture_sum_err column output by
aperture_photometry (divided by epadu) will be used as the Poisson noise term.
bg_method: string
{'mean', 'median', 'mode'}, optional. The statistic used to calculate the background. All measurements are
sigma clipped. Default value is 'mode'. NOTE: From DAOPHOT, mode = 3 * median - 2 * mean.
epadu : float
(optional) Gain in electrons per adu (only use if image units aren't e-). Default value is 1.0
Returns
-------
An astropy Table with columns as follows:
X-Center Y-Center RA DEC ID MagAp1 MagErrAp1 MagAp2 MagErrAp2 MSkyAp2 StdevAp2 FluxAp2 CI Flags
"""
if bg_method not in ['mean', 'median', 'mode']:
raise ValueError('Invalid background method, choose either \
mean, median, or mode')
phot = aperture_photometry(data, phot_apertures, error=error_array)
bg_phot = aperture_stats_tbl(data, bg_apertures, sigma_clip=True)
names = ['X-Center', 'Y-Center', 'ID']
x, y = phot_apertures[0].positions.T
final_stacked = np.stack([x, y, phot["id"].data], axis=1)
# n_aper = 0
name_list = 'Flux', 'FluxErr', 'Mag', 'MagErr'
for aper_string in ['Ap1', 'Ap2']:
for col_name in name_list:
names.append("{}{}".format(col_name, aper_string))
# for item in list(phot.keys()):
# if item.startswith("aperture_sum_") and not item.startswith("aperture_sum_err_"):
# aper_size_arcsec = phot_apertures[n_aper].r * platescale
# for name in name_list:
# names.append("{}_{:.2f}".format(name, aper_size_arcsec))
# n_aper += 1
for aperCtr in range(0, 2):
ap_area = phot_apertures[aperCtr].area
bg_method_name = 'aperture_{}'.format(bg_method)
flux = phot['aperture_sum_{}'.format(aperCtr)] - bg_phot[bg_method_name] * ap_area
# Need to use variance of the sources
# for Poisson noise term in error computation.
#
# This means error needs to be squared.
# If no error_array error = flux ** .5
if error_array is not None:
flux_error = compute_phot_error(phot['aperture_sum_err_{}'.format(aperCtr)] ** 2.0, bg_phot, bg_method,
ap_area, epadu)
else:
flux_error = compute_phot_error(flux, bg_phot, bg_method, ap_area, epadu)
# When photflam is 0.0, mag is set to -9999.0.
mag = convert_flux_to_abmag(flux, photflam, photplam)
# NOTE: Magnitude error calculation comes from computing d(ABMAG)/d(flux).
# See https://iraf.net/forum/viewtopic.php?showtopic=83932 for details.
if math.isclose(photflam, 0.0, abs_tol=constants.TOLERANCE):
mag_err = mag
else:
mag_err = 1.0857 * flux_error / flux
# Set magnitude errors to infinity if flux is negative
mag_err[np.logical_not(flux>0)]=np.inf
# Build the final data table
stacked = np.stack([flux, flux_error, mag, mag_err], axis=1)
final_stacked = np.concatenate([final_stacked, stacked], axis=1)
# Build final output table
final_tbl = Table(data=final_stacked, names=names,
dtype=[np.float64, np.float64, np.int64, np.float64, np.float64, np.float64, np.float64,
np.float64, np.float64, np.float64, np.float64])
# add sky and std dev columns from background calculation subroutine
final_tbl.add_column(bg_phot[bg_method_name])
final_tbl.rename_column(bg_method_name, 'MSkyAp2')
final_tbl.add_column(bg_phot['aperture_std'])
final_tbl.rename_column('aperture_std', 'StdevAp2')
return final_tbl
[docs]
def compute_phot_error(flux_variance, bg_phot, bg_method, ap_area, epadu=1.0):
"""Computes the flux errors using the DAOPHOT style computation. Originally taken from
https://github.com/spacetelescope/wfc3_photometry/blob/527815bf580d0361754281b608008e539ed84368/photometry_tools/photometry_with_errors.py#L137
Parameters
----------
flux_variance : array
flux values
bg_phot : array
background brightness values.
bg_method : string
background method
ap_area : array
the area of the aperture in square pixels
epadu : float
(optional) Gain in electrons per adu (only use if image units aren't e-). Default value is 1.0
Returns
-------
flux_error : array
an array of flux errors
"""
bg_variance_terms = (ap_area * bg_phot['aperture_std'] ** 2.) * (1. + ap_area / bg_phot['aperture_area'])
variance = flux_variance / epadu + bg_variance_terms
flux_error = variance ** .5
return flux_error
[docs]
def convert_flux_to_abmag(in_flux, photflam, photplam):
"""converts flux (in units of electrons/sec) to ABMAG
Parameters
----------
in_flux : astropy.table.column.Column object
flux values to convert to ABMAG, in electrons/second
photflam : float
inverse sensitivity, in ergs/cm2/angstrom/electron
photplam : float
pivot wavelength, in angstroms
Returns
-------
abmag : astropy.table.column.Column object
input flux values converted to ABMAG
"""
# Avoid taking the log10 of zero
if math.isclose(photflam, 0.0, abs_tol=constants.TOLERANCE):
# Create a column of data from an existing column of data
# and set each value = -9999.0
abmag = constants.FLAG + in_flux * 0.0
else:
# convert flux from units of electrons/second to ergs/cm2/angstrom/second
f_lambda = in_flux * photflam
# Convert f_lambda to STMAG
stmag = -2.5 * np.log10(f_lambda) - 21.10
# Convert STMAG to ABMAG
abmag = stmag - 5.0 * np.log10(photplam) + 18.6921
return abmag