Reproduce IRAF's phot task 'merr' value with aperture_photometry

[Gabriel-p]

The phot task performs aperture photometry, returning (among many other things) the magnitude and magnitude uncertainty for each star. This uncertainty is called merr and is obtained as:

merr = 1.0857 * error / flux

where

error = sqrt(flux / epadu + area * stdev**2 +  area**2 * stdev**2 / nsky)

with

stdev: standard deviation (of the sky value)
epadu: Gain (electrons  per adu)
nsky: number of sky pixels

How can I reproduce this value using the aperture_sum_err returned by aperture_photometry() (shown here)?

This is what I have so far, but I'm not sure what to use for bkg_error.

# Assuming 'hdu_data' and 'gain' are known, what do I use for 'bkg_error'?
error = calc_total_error(hdu_data, bkg_error, gain)

# Define apertures based on given 'positions'  and 'aper_rad' value
apertures = CircularAperture(positions, r=aper_rad)
annulus_apertures = CircularAnnulus(positions, r_in=aper_rad + 5, r_out=aper_rad + 10.)
apers = [apertures, annulus_apertures]

# Aperture photometry
phot_table = aperture_photometry(hdu_data, apers, error=error)
bkg_mean = phot_table['aperture_sum_1'] / annulus_apertures.area()
bkg_sum = bkg_mean * apertures.area()
phot_table['flux_fit'] = phot_table['aperture_sum_0'] - bkg_sum

# Obtain 'merr'
phot_table['merr'] = 1.0857 *phot_table['aperture_sum_err_0'] / phot_table['flux_fit']

[Vb2341]

To get IRAF's merr nearly exactly, you actually need to do local background computation, with local background error estimation as well. This is tricky using PhotUtils alone, but I've actually implemented it rather recently (though it seems @larrybradley also did it very recently here. Anyway, I've simply just taken it a step further and computed the DAOPHOT style errors by implementing the same formula in Python.

The code lives over at https://github.com/spacetelescope/wfc3_photometry/blob/master/photometry_tools/photometry_tools.py. An example is in the docstring, but for your posted code I think it should be something like:

from photometry_tools import iraf_style_photometry

photometry_tbl = iraf_style_photometry(
        apertures, annulus_apertures, hdu_data,
        error_array=None, bg_method='mean') # median and mode are also valid options

This should get you errors more in line with what IRAF would produce. There are a few small caveats, but generally it seems the magnitudes agree to ~0.01 mags on average (1%) and the errors agree to ~.001 mag (the fainter stars are worse, though I think this is more an IRAF limitation).

[Gabriel-p]

Thank you very much for your answer @Vb2341. One small note: your compute_phot_error function seems to assume an epadu (Gain) of 1. with no possibility to pass a different value from iraf_style_photometry. Is this for a reason or you just forgot to add the epadu parameter to iraf_style_photometry?

[Vb2341]

@Gabriel-p No problem!

The reason I didn't really think about it was mainly because I was developing these tools for HST photometry, where our images are already typically in electrons rather than DN. I can easily add it as a parameter, but I'm not sure if I'll have to apply the gain factor anywhere else. If the image is in DN, then I think it should just be the error, but I could be wrong. I suppose if the results don't match DAOPHOT afterwards, then the answer will come out of there. I'll add the epadu to the iraf_style_photometry definition within the hour.

Edit: It's added: spacetelescope/wfc3_photometry@4cdd941

[Gabriel-p]

Thanks again @Vb2341! I'll check this out, compare with IRAF, and report back. If the results are close enough, I'll close the issue.