photutils Aperture Photometry¶

import matplotlib.pyplot as plt
import numpy as np
from astropy.io import fits
from astropy.stats import sigma_clip, sigma_clipped_stats
from photutils import DAOStarFinder, Background2D, SigmaClip, SExtractorBackground
from photutils import CircularAperture as CircAp
from photutils import CircularAnnulus as CircAn
from photutils import aperture_photometry as APPHOT

#%%
hdu = fits.open('HST_Tutorial/M13.fits')
img = hdu[0].data[900:1200, 900:1200]
# if image value < 10^(-6), replace the pixel as 10^(-6)
img[img < 1.e-6] = 1.e-6

# multiply 10**4 and integrize to mimic ADU
img *= 1.e4
img = img.astype(int)
ronoise = 5  # electrons
gain = 7     # e/ADU

FWHM   = 2.5
sky_th = 500   # sky_th * sky_sigma will be used for detection lower limit
sky_a, sky_m, sky_s  = sigma_clipped_stats(img) # 3-sigma, 5 iters
thresh = sky_th * sky_s

find   = DAOStarFinder(fwhm=FWHM, threshold=thresh,
                       sharplo=0.2, sharphi=1.0,  # default values 
                       roundlo=-1.0, roundhi=1.0, # default values
                       sigma_radius=1.5,          # default values
                       ratio=1.0,                 # 1.0: circular gaussian
                       exclude_border=True)       # To exclude sources near edges
found = find(img)

N_star = len(found)
# save XY coordinates and aperture:
coord = (found['xcentroid'], found['ycentroid']) 
apert = CircAp(positions=coord, r=max(3, FWHM))
annul = CircAn(positions=coord, r_in=4*FWHM, r_out=6*FWHM)

def mag_inst(flux, ferr):
    ''' Returns magnitude from flux.
    '''
    m_inst = -2.5 * np.log10(flux)
    merr   = 2.5/np.log(10) * ferr / flux
    return m_inst, merr

def sky_fit(all_sky, method='mode', sky_nsigma=3, sky_iter=5, \
            mode_option='sex', med_factor=2.5, mean_factor=1.5):
    '''
    Estimate sky from given sky values.

    Parameters
    ----------
    all_sky : ~numpy.ndarray
        The sky values as numpy ndarray format. It MUST be 1-d for proper use.
    method : {"mean", "median", "mode"}, optional
        The method to estimate sky value. You can give options to "mode"
        case; see mode_option.
        "mode" is analogous to Mode Estimator Background of photutils.
    sky_nsigma : float, optinal
        The input parameter for sky sigma clipping.
    sky_iter : float, optinal
        The input parameter for sky sigma clipping.
    mode_option : {"sex", "IRAF", "MMM"}, optional.
        sex  == (med_factor, mean_factor) = (2.5, 1.5)
        IRAF == (med_factor, mean_factor) = (3, 2)
        MMM  == (med_factor, mean_factor) = (3, 2)

    Returns
    -------
    sky : float
        The estimated sky value within the all_sky data, after sigma clipping.
    std : float
        The sample standard deviation of sky value within the all_sky data,
        after sigma clipping.
    nsky : int
        The number of pixels which were used for sky estimation after the
        sigma clipping.
    nrej : int
        The number of pixels which are rejected after sigma clipping.
    -------

    '''
    sky = all_sky.copy()
    if method == 'mean':
        return np.mean(sky), np.std(sky, ddof=1)

    elif method == 'median':
        return np.median(sky), np.std(sky, ddof=1)

    elif method == 'mode':
        sky_clip   = sigma_clip(sky, sigma=sky_nsigma, iters=sky_iter)
        sky_clipped= sky[np.invert(sky_clip.mask)]
        nsky       = np.count_nonzero(sky_clipped)
        mean       = np.mean(sky_clipped)
        med        = np.median(sky_clipped)
        std        = np.std(sky_clipped, ddof=1)
        nrej       = len(all_sky) - len(sky_clipped)

        if nrej < 0:
            raise ValueError('nrej < 0: check the code')

        if nrej > nsky: # rejected > survived
            raise Warning('More than half of the pixels rejected.')

        if mode_option == 'IRAF':
            if (mean < med):
                sky = mean
            else:
                sky = 3 * med - 2 * mean

        elif mode_option == 'MMM':
            sky = 3 * med - 2 * mean

        elif mode_option == 'sex':
            if (mean - med) / std > 0.3:
                sky = med
            else:
                sky = (2.5 * med) - (1.5 * mean)
        else:
            raise ValueError('mode_option not understood')

        return sky, std, nsky, nrej