# Copyright 2014 - 2015 Anna-Lea Lesage for MASCARA
#
# This file is part of the mascara package.
# mascara is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# mascara is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with mascara. If not, see <http://www.gnu.org/licenses/>.
'''
.. module:: phot
.. moduleauthor:: Anna-Lea Lesage
:mod:`mascara.phot` provides a set of photometric routines:
- find_stars, an export from the DAOPHOT package, to find stars on an image
- GetPhot, the equivalent of the DAOPHOT APER routine.
- OnaPhot, an extensive version of GetPhot which performs quick image quality
analysis, and returns more informations
'''
__version__ = '15.10.12'
### *** Latest Changes Record *** ###
### 09.12.2014 Changed in Find the output of the function is No stars are found.
### Instead of returning just None, it is now returning 4 times None
### for (x, y, sharpness, roundness). This allows an easier condition
### on the success of star detection.
### 09.12.2014 Added a condition in find that if the std of the image is null, it
### returns only None. Goal is to skip empty images
### 07.01.2015 Found an error in the calculation of the mask. Changed around line 451.
### Should result in lower scatter in the light curve for small apertures.
### 09.01.2015 Added a getSkyValue function to get more robust estimates of the
### sky. It is based on the IDL mmm.pro function. In order to keep fast
### and robust, only the first pass is used. (Though the rest is coded)
### 13.01.2015 Added keywords precise_sky and onepass to OnaPhot to give the choise of
### the method used for calculating the background
### 13.01.2015 Added function makeMask, using meshgrid
### 09.03.2015 Rewrote the centroiding rountine in find to get accurate results!
### 26.05.2015 Added in the output of OnaPhot the max pixel value inside the aperture
### 01.09.2015 Updated GetPhot for the find_pointing routine: changed empty arrays into
### zero arrays.
### ******************************* ###
[docs]def find_stars(image, hmin, fwhm, nsigma=1.5, roundlim=[-1.,1.], sharplim=[0.2,1.]):
""" find_stars, adapted from the ASTROLIB-routine to find point sources on an image. To do so,
the image is convolved with a Gaussian kernel of fwhm as given. Any remaining sources
above the background are then potential stars...
Inputs:
- image, the 2D image to look at
- hmin, the noise threshold above which the object is qualified as a star. Typically hmin > 3*std(image)
- fhwm, the estimated full width half maximum of the stars on the image
Keywords:
- nsigma, radius of the convolution kernel. Default 1.5
- roundlim ([-1.,1.]): Threshold for the roundness criterion.
- sharplim ([0.2,1.]): Threshold for the sharpness criterion.
Outputs:
- xobs, yobs, the X and Y coordinates of the stars found on the CCD
- sharp, the sharpness of the stars
- round, the roundness of the stars
"""
from scipy.ndimage import convolve
import numpy as np
from math import sqrt
import bottleneck as bt
### Getting image information
ny, nx = image.shape
### Quick analyse to check is the image is not empty:
if image.std() == 0:
return None, None, None, None
################################
# Setting the convolution kernel
################################
sigmatofwhm = 2*sqrt(2*np.log(2))
radius = nsigma * fwhm / sigmatofwhm # Radius is 1.5 sigma
if radius < 1.0:
radius = 1.0
fwhm = sigmatofwhm/nsigma
print( "WARNING!!! Radius of convolution box smaller than one." )
print( "Setting the 'fwhm' to minimum value, %f, given the provided nsigma." %fwhm )
sigsq = (fwhm/sigmatofwhm)**2 # sigma squared
nhalf = int(radius) # Center of the kernel
nbox = 2*nhalf+1 # Number of pixels inside of convolution box
middle = nhalf # Index of central pixel
### Define x,y coordinates of the kernel
kerny, kernx = np.meshgrid(range(nbox),range(nbox))
### Compute the square of the distance to the center
g = np.float64((kernx-nhalf)**2+(kerny-nhalf)**2)
mask = g <= radius**2 # We make a mask to select the inner circle of radius "radius"
nmask = mask.sum() # The number of pixels in the mask within the inner circle.
g = np.exp(-0.5*g/sigsq) # We make the 2D gaussian profile
### Convolving the image with a kernel representing a gaussian
c = g*mask
### Normalisation of the Gaussian kernel
c[mask] = (c[mask] - c[mask].mean())/(c[mask].var() * nmask)
c1 = g[nhalf] # c1 will be used to the test the roundness
c1 = (c1-c1.mean())/((c1**2).sum() - c1.mean())
h = convolve(image,c,mode='constant',cval=0.0)
h[:nhalf,:] = 0 # Set the sides to zero in order to avoid border effects
h[-nhalf:,:] = 0
h[:,:nhalf] = 0
h[:,-nhalf:] = 0
###########################
### Find source candidates:
###########################
mask[middle,middle] = False # From now on we exclude the central pixel
nmask = mask.sum() # so the number of valid pixels is reduced by 1
goody,goodx = mask.nonzero() # "good" identifies position of valid pixels
### Identifying the point source candidates that stand above the background
indy,indx = (h >= hmin).nonzero()
nfound = indx.size # nfound is the number of candidates
### If we didn't find any stars, return x, y, sharpness, roundness = None
if nfound <= 0:
return None, None, None, None
# a (nfound, nmask) array of good positions in the mask, mask coordinate
offsetsx = np.resize(goodx-middle,(nfound,nmask))
offsetsx = indx + offsetsx.T
# a (nfound, nmask) array of good positions in the mask, mask coordinate
offsetsy = np.resize(goody-middle,(nfound,nmask))
offsetsy = indy + offsetsy.T
### Get values in the image with offsets
offsets_vals = h[offsetsy,offsetsx]
### Get pixel values in the images wihtout offsets
vals = h[indy,indx]
### Identifying the candidates that are local maxima (hence larger than offset)
ind_goodcandidates = ((vals - offsets_vals) > 0).all(axis=0)
nfound = ind_goodcandidates.sum() # update the number of candidates
if nfound <= 0:
return None, None, None, None
### Keep the local maximas
indx = indx[ind_goodcandidates]
indy = indy[ind_goodcandidates]
##############################################################
# Identifying the candidates that meet the sharpness criterion
##############################################################
d = h[indy,indx] # a (nfound) array of the intensity of good candidates
d_image = image[indy,indx]
offsetsx = offsetsx[:,ind_goodcandidates]
offsetsy = offsetsy[:,ind_goodcandidates]
offsets_vals = image[offsetsy,offsetsx]
sharpness = (d_image - offsets_vals.sum(0)/nmask) / d
ind_goodcandidates = (sharpness > sharplim[0]) * (sharpness < sharplim[1])
# a (nfound) array of candidates that meet the sharpness criterion
nfound = ind_goodcandidates.sum() # update the number of candidates
if nfound <= 0:
return None, None, None, None
## Update the good source candidates
indx = indx[ind_goodcandidates]
indy = indy[ind_goodcandidates]
# update sharpness with the good candidates
sharpness = sharpness[ind_goodcandidates]
##############################################################
# Identifying the candidates that meet the roundness criterion
##############################################################
temp = np.arange(nbox)-middle # make 1D indices of the kernel box
temp = np.resize(temp, (nbox,nbox)) # make 2D indices of the kernel box (for x or y)
offsetsx = np.resize(temp, (nfound,nbox,nbox))
# make 2D indices of the kernel box for x, repeated nfound times
offsetsy = np.resize(temp.T, (nfound,nbox,nbox))
# make 2D indices of the kernel box for y, repeated nfound times
offsetsx = (indx + offsetsx.swapaxes(0,-1)).swapaxes(0,-1) # make it relative to image coordinate
offsetsy = (indy + offsetsy.swapaxes(0,-1)).swapaxes(0,-1) # make it relative to image coordinate
offsets_vals = image[offsetsy,offsetsx]
# a (nfound, nbox, nbox) array of values (i.e. the kernel box values for each nfound candidate)
dx = (offsets_vals.sum(2)*c1).sum(1)
dy = (offsets_vals.sum(1)*c1).sum(1)
roundness = 2*(dx-dy)/(dx+dy)
ind_goodcandidates = (roundness > roundlim[0]) * (roundness < roundlim[1]) * (dx >= 0.) * (dy >= 0.)
# a (nfound) array of candidates that meet the roundness criterion
nfound = ind_goodcandidates.sum() # update the number of candidates
if nfound <= 0:
return None, None, None, None
## Update the good source candidates
indx = indx[ind_goodcandidates]
indy = indy[ind_goodcandidates]
sharpness = sharpness[ind_goodcandidates]
roundness = roundness[ind_goodcandidates]
#return indx, indy, sharpness, roundness
###Preprate centroiding
#g[0,:] = g[-1,:] = 1.
width = nhalf - np.abs(np.arange(nbox)-nhalf)+1
xwidth, ywidth = np.meshgrid(width, width)
sumgaussx = np.sum(g*xwidth,axis=1)
sumgaussy = np.sum(g*ywidth,axis=0)
P = np.sum(width)
sumgx = np.sum(width*sumgaussy)
sumgy = np.sum(width*sumgaussx)
vec = nhalf - np.arange(nbox)
dgdx = sumgaussy*vec
dgdy = sumgaussx*vec
x = np.array(indx, dtype=float)
y = np.array(indy, dtype=float)
for i in xrange(nfound):
temp = image[indy[i]-nhalf:indy[i]+nhalf+1, indx[i]-nhalf:indx[i]+nhalf+1]
### Find X centroid
sd = np.sum(temp*ywidth, 0)
sumofwsgsd = np.sum(width*sumgaussy*sd)
sumofwsd = np.sum(width*sd)
sumofwsddx = np.sum(width*sd*dgdx)
HeightX = (sumofwsgsd-sumgx*sumofwsd/P)/(np.sum(width*sumgaussy**2)-sumgx**2/P)
skylevel = (sumofwsd - HeightX*sumgx)/P
### Calculating the dx for centroiding
dx = (np.sum(width*dgdx*sumgaussy)-(sumofwsddx-np.sum(width*dgdx)*(HeightX*sumgx + \
skylevel*P)))/(HeightX*np.sum(width*dgdx**2)/sigsq)
if np.abs(dx) >= nhalf:
dx = 0.
x[i] += dx
### Find Y centroid
sd = np.sum(temp*xwidth,1)
sumofwsgsd = np.sum(width*sumgaussx*sd)
sumofwsd = np.sum(width*sd)
sumofwsddy = np.sum(width*sd*dgdy)
HeightY = (sumofwsgsd-sumgy*sumofwsd/P)/(np.sum(width*sumgaussx**2)-sumgy**2/P)
skylevel = (sumofwsd - HeightY*sumgy)/P
### Calculating the dy for centroiding:
dy = (np.sum(width*dgdy*sumgaussx)-(sumofwsddy-np.sum(width*dgdy)*(HeightY*sumgy + \
skylevel*P)))/(HeightY*np.sum(width*dgdy**2)/sigsq)
if np.abs(dy) >= nhalf:
dy = 0.
y[i] += dy
if nfound == 1:
return np.array(x), np.array(y), sharpness, roundness
### Remove possible duplicates of the same star
xf, yf = [],[]
x = x[np.isfinite(x)]
y = y[np.isfinite(y)]
index = list(np.arange(x.size))
for indi in index:
tmp = ((np.abs(x-x[indi]) < 2.*fwhm)*(np.abs(y-y[indi]) < 2.*fwhm)).nonzero()
xc, yc = x[tmp], y[tmp]
xf.append(xc.mean()), yf.append(yc.mean())
if len(tmp)>1:
for w in tmp[1:]:
if w > indi:
index.remove(w)
return np.array(xf), np.array(yf), sharpness, roundness
[docs]def GetPhot(im, xc, yc, apr=None, skyrad=None, phpadu=1, Flux=True):
""" Aperture Photometry based on DAOPHOT Aper procedure
Computes concentric aperture photometry for several user specified aperture radii
(apr). The sky is evaluated separately using the outer and inner sky radii (skyrad)
Inputs :
- im, input image array
- x, y the coordinates of the stars on the array
Keywords:
- apr, a single value (can be float) for the aperture radius
- skyrad, a tupple of 2 elements defining the radii of the sky annulus
- phpadu, the photon per ADU (optional)
Outputs:
- flux, the measured flux
- eflux, the errors on the flux
- sky, the measured background
- esky, the errors on the background
"""
from math import floor
import numpy as np
badval = [0, 63000]
dim = im.shape
if dim[0] == 1:
return np.zeros(1)
if len(dim) == 3:
ny, nx, ndim = dim
else:
ny, nx = dim
### Verify how many coordinates are given as input
Nstars = np.size(xc)
if Nstars == 1:
xc = np.array(xc)
yc = np.array(yc)
elif Nstars == 0:
return
### Verify that the aperture radius is given
napr = np.size(apr)
if apr == None:
apr = input('Give aperture radii: ')
napr = np.size(apr)
apr = np.array([apr])
# Some failsafes : ask for aperture radius if not given,
# ask for sky radii if skyrad is an integer, ask for phpadu if not given
if skyrad == None:
skyrad = input('Give inner and outer sky radii: ')
if phpadu == 0:
phpadu = input('Give the photon to ADU convertion: ')
if napr > 1:
flux = np.zeros([Nstars, napr])
erflux= np.zeros([Nstars, napr])
else:
flux = np.zeros(Nstars)
erflux= np.zeros(Nstars)
sky = np.zeros(Nstars)
ersky = np.zeros(Nstars)
# Set the limit to the edges of the CCD
lx = np.floor(xc - skyrad[1])
upx = np.floor(xc + skyrad[1])
ly = np.floor(yc - skyrad[1])
upy = np.floor(yc + skyrad[1])
# Attention to the edges of the CCD (remove the outer sky radius)
lx = lx*(lx > 0)
ly = ly*(ly > 0)
upx = upx*(upx < nx-0)
upy = upy*(upy < ny-0)
# Keep only the good values where the star is not on the edges or outside the CCD
dnx = upx-lx
dny = upy-ly
# Define the stars that lie on the edges or out of the CCD
edges = (lx == 0)+(upx == 0)+(ly == 0)+(upy == 0)
# Get the subpixel information in the subarray
dx = xc - lx
dy = yc - ly
# Defining the aperture areas
area = np.pi * apr**2
# Start aperture photometry
for i in xrange(Nstars):
# If the star is not on the edge? or if the star is on the edge...
try:
tmp = edges[i]
except IndexError:
tmp = edges
if tmp == False:
skymod = 0
apmag = np.zeros([napr])
magerr = np.zeros([napr])
error1 = np.zeros([napr])
error2 = np.zeros([napr])
error3 = np.zeros([napr])
try:
box = im[ly[i]:upy[i], lx[i]:upx[i]] # Extract the subimage centered on the star
except IndexError:
box = im[ly:upy, lx:upx]
# Define the distances of each pixel to the central pixel
try:
dxsq = (np.arange(dnx[i])-dx[i])** 2
rsq = np.zeros([dny[i], dnx[i]])
except IndexError:
dxsq = (np.arange(dnx)-dx)**2
rsq = np.zeros([dny,dnx])
try:
for j in xrange(0, int(dny[i])):
rsq[j,:] = dxsq + (j-dy[i])**2
except IndexError:
for j in xrange(0, int(dny)):
rsq[j,:] = dxsq + (j-dy)**2
r = np.sqrt(rsq) - 0.5
# Select the pixels inside the sky annulus
skypix = (rsq >= skyrad[0]**2)*(rsq <= skyrad[1]**2)
skyvalues = box[skypix]
# Skip the meanclip calculation of the sky for the moment
# Rough estimate of the sky mode
skymed = np.median(skyvalues) # Sky median
skymean = skyvalues.mean()
skymod = 3*skymed - 2*skymean # Evaluation of the mode... under the assumption
# of slight sky contamination
skyvar = skyvalues.var() # Sky variance
sigsq = skyvar/len(skyvalues) # Square of the standard error by mean sky brightness
# Go through the aperture for the star
for k, APR in enumerate(apr):
indexap = (r < APR)
inap = box[indexap]
rap = r[indexap]
fractn = (APR-rap) # fraction of pixel to count
fullpix = (fractn >= 1.).nonzero()[0]
fractn[fullpix] = 1.0
fcounts = fullpix.size
fracpix = (fractn < 1.).nonzero()[0]
factor = (area[k] - fcounts)/ sum(fractn[fracpix])
fractn[fracpix] = fractn[fracpix]*factor
apmag[k] = sum(inap*fractn)
if Flux:
# Keep the finite values:
gd = np.isfinite(apmag)
else:
tmp = np.abs(apmag - badval)*(np.abs(apmag-badval) > 0.01)
gd = tmp.nonzero()
apmag[gd] = apmag[gd] - skymod*area[gd]
# Add the errors
error1[gd] = area[gd]*skyvar # Scatter in the sky
error2[gd] = (apmag[gd])/phpadu # Random photon noise
error3[gd] = sigsq*area[gd]**2 # Uncertainties in mean sky brigtness
magerr[gd] = np.sqrt(error1[gd] + error2[gd]+error3[gd])
if not Flux:
good = apmag.nonzero()
if len(apmag[good]) > 0:
magerr[good] = 1.0857*magerr[good]/apmag[good]
# supposely log(10)/2.5=1.0857
apmag[good] = 25. - 2.5*np.log10(apmag[good])
# Store the results in ouput variables
sky[i] = skymod
ersky[i] = np.sqrt(skyvar)
if napr > 1:
flux[i,:] = apmag
erflux[i,:] = magerr
else:
flux[i] = apmag
erflux[i] = magerr
return flux, erflux, sky, ersky
[docs]def OnaPhot(*args, **kwargs):
''' OnaPhot is an aperture photometry routine based on the APER or GetPhot
routine, but developped for the following cases:
- when one star is treated at the time
- when the shape of the star is highly dependant of its position on the
CCD
OnaPhot adapt the shape of the aperture to the shape of the star.
In addition, it runs a diagnostic of the background.
OnaPhot accepts one to 10 different aperture radii.
Input:
- im, the 2d complete image
- xc, yc the coordinates of ONE star.
Keywords:
- apr, a vector, list or tuple of aperture radii
- skyrad, a tuple [inner sky rad, outer sky rad]
- mask, an 2D array. If not provided, it is calculated.
It may be interesting to provide mask, (via masphot.makeMask)
- precise_sky, to do a precise calculation of the sky background
- onepass, if only one pass of the getSkyValue should be done.
Outputs:
- flag, a diagnostic for the star noted in binary:
- 0 -> all went well
- 1 -> saturated pixel in the aperture
- 2 -> saturated pixel in the background
- 4 -> reserved (for astrometry)
- 8 -> reserved (for clouds)
- 16 -> Background <= 0
- 32 -> Flux or sky <= 0 or not finite
- 64 -> --
- 128 -> Star potentially outside the aperture
- 256 -> Sky annulus outside CCD
- flux, a vector of flux
- erflux, the corresponding errors
- sky, the sky background
- ersky, the corresponding errors
- xm, ym, the coordinates of the star
The output is expressed under the form:
[flag, flux0, erflux0, flux1, erflux1, ..., sky, ersky, xm, ym]
'''
from mascara.funcs.utils import percentile, Gfitting
from math import sqrt, floor
import numpy as np
import bottleneck as bt
### Define the 2 cases for centering the star:
if len(args) ==3:
im, xc, yc = args
centering = kwargs.pop('centering', False)
elif len(args) ==2:
im, xc = args
yc = xc
centering = kwargs.pop('centering', True)
### Reading in the other eywords
apr = kwargs.pop('apr', None)
skyrad = kwargs.pop('skyrad', None)
phpadu = kwargs.pop('phpadu', 0.62)
precise = kwargs.pop('precise_sky', False)
onepass = kwargs.pop('onepass', True)
ny, nx = im.shape
apr = np.array(apr)
napr = apr.size
if ny == 1:
return np.zeros(2*napr + 3 + 3)
badpixval = 63000.
flag = kwargs.pop('flag', 0)
flux = np.zeros([napr, 2])
# 1. Check if the outer skyrad is still in the CCD:
outskyradin = ((xc>=skyrad[1])*(xc<=nx-skyrad[1])*(yc>=skyrad[1])*(yc<=ny-skyrad[1]))
if outskyradin == False:
flag = 256.
flux = np.ravel(flux)
sky = np.zeros([2])
return np.hstack((flag, flux, sky, 0., xc, yc))
# 2. Check if the star is inside the aperture radius:
# If several apertures are given, use the largest one to center the star!
lx = floor(xc - skyrad[1])
upx = floor(xc + skyrad[1])
ly = floor(yc - skyrad[1])
upy = floor(yc + skyrad[1])
mask = makeMask(skyrad[1], xc, yc)
# 3. Recenter the star if keyword is set
xm, ym = np.copy(xc), np.copy(yc)
if centering:
from scipy.ndimage import center_of_mass
box = im[ly:upy, lx:upx]
thresh = percentile(box, 0.75)
cy, cx = center_of_mass((box > thresh)*(box<badpixval))
# 4. Update coordinates
xm, ym = xc + cx -skyrad[1], yc+cy-skyrad[1]
lx = np.floor(xm - skyrad[1])
upx = np.floor(xm + skyrad[1])
ly = np.floor(ym - skyrad[1])
upy = np.floor(ym + skyrad[1])
mask = makeMask(skyrad[1], xm, ym)
##### ####
box = im[ly:upy, lx:upx]
# 5. Check the sky background for irregularities.
outermask = ((mask+0.5)**2<=skyrad[1]**2)*((mask+0.5)**2>=skyrad[0]**2)
skydonut = box[outermask]
if precise:
skymod, skystd = getSkyValue(skydonut, onepass=onepass)
else:
skymed = bt.median(skydonut) # Sky median
skymean = bt.nanmean(skydonut)
skymod = 3*skymed - 2*skymean # Evaluation of the mode... under the assumption
# of slight sky contamination
skystd = bt.nanstd(skydonut)
### Flag negative sky measurements
if skymod < 0 :
flag = flag+16
skyvar = skystd**2 # Sky variance
sigsq = skyvar/len(skydonut)
outer = percentile(skydonut, 0.75)
indexap = (mask < bt.nanmax(apr))
inap = box[indexap]
amax = inap.max()
inner = percentile(inap, 0.9)
### Check if the aperture is brighter than the sky
if inner > outer:
### Flag if there is a bad pixel in the skydonut
if bt.nanmax(skydonut) >= badpixval:
flag = flag + 2
else:
sky = np.array([skymod, skystd])
flux = np.ravel(flux)
flag = 128 + flag
return np.hstack((flag, flux, sky, amax, xm, ym))
# Now do the photometry!
for k, APR in np.ndenumerate(apr):
area = np.pi * APR**2
indexap = mask < APR
inap = box[indexap]
## Flag if there is a bad pixel in the aperture
if bt.nanmax(inap) >= badpixval:
flag = flag + 1
rap = mask[indexap]
fractn = (APR-rap) # fraction of pixel to count
fullpix = (fractn >= 1.).nonzero()[0]
fractn[fullpix] = 1.0
fcounts = fullpix.size
fracpix = (fractn < 1.).nonzero()[0]
factor = (area - fcounts)/ bt.nansum(fractn[fracpix])
fractn[fracpix] = fractn[fracpix]*factor
apmag = bt.nansum(inap*fractn)
apmag = apmag - skymod*area
## Flag the negative stellar fluxes
if apmag < 0:
flag = flag + 32
error1 = area*skystd**2
error2 = apmag / phpadu
error3 = sigsq*area**2
magerr = np.sqrt(error1 + error2 + error3)
flux[k,:] = apmag, magerr
sky = np.array([skymod, skystd])
flux = np.ravel(flux)
return np.hstack((flag, flux, sky, amax, xm, ym))
def getSkyValue(sky, onepass= True):
''' getSkyValue calculates the background value even in the case of background
contamination.
Adapted from the mmm module of the idl daophot library
Inputs:
- sky
Keyords:
- onepass, if true, the function return the sky's mode and standard
deviation after only one pass.
Outputs:
- Skymod and sky standard deviation
'''
import numpy as np
import bottleneck as bt
minsky = 20.
nlast = sky.size -1
sky = bt.partsort(sky,1)
skymid = bt.nanmedian(sky)
cut1 = bt.nanmin([skymid-sky[0], sky[-1]-skymid])
cut2 = skymid+cut1
cut1 = skymid-cut1
## select the pixels between cut1 and cut2
good = ((sky <= cut2)*(sky>=cut1)).nonzero()[0]
Ng = good.size
## If no good pixel then...
if Ng == 0:
return skymid, 0.0
delta = sky[good] - skymid
Totdelta = bt.nansum(delta)
Totdeltasq = bt.nansum(delta**2)
Hsky = bt.nanmax(good)
Lsky = bt.nanmin(good)-1
## Compute mean and sigma from the first pass
skymed = 0.5*sky[(Lsky+Hsky+1)/2] + 0.5*sky[(Lsky+Hsky)/2+1]
skymn = Totdelta/(Hsky-Lsky)
sigma = sqrt(Totdeltasq/(Hsky-Lsky) - skymn**2)
skymn = skymn+skymid
if skymed < skymn:
skymod = 3.*skymed - 2.*skymn
else:
skymod = skymn
if onepass:
return skymod, sigma
## Rejection and recomputation loop
niter = 0
old = 0.
clamp = 0.
redo = True
while (niter < 20) and redo:
niter += 1
redo = False
print 'Iteration number %g' %niter
## Compute Chauvenet rejection criterion
r = np.log10(Hsky-Lsky)
r = bt.nanmax([2., (-0.1042*r+1.1695)*r + 0.8895])
## Compute the rejection limits
cut = r*sigma + 0.5*np.abs(skymn-skymod)
cut1 = skymod-cut
cut2 = skymod+cut
## recompute mean and sigma by adding/subtracting sky values at
## both end of the interval of acceptable values
newmin = Lsky
tst_min = sky[newmin+1] >= cut1
done = (newmin == -1)*tst_min
if not done:
done = (sky[newmin>0] > cut1) and tst_min
istep = 1 - 2*tst_min
while not done:
newmin += istep
done = (newmin == -1) or (newmin == nlast)
if not done:
done = (sky[newmin] <= cut1) and (sky[newmin+1] >= cut1)
if tst_min:
delta = sky[newmin+1: Lsky] - skymid
else:
delta = sky[Lsky+1:newmin] - skymid
Totdelta = Totdelta - istep * bt.nansum(delta)
Totdeltasq = Totdeltasq - istep * bt.nansum(delta**2)
Lsky = newmin
redo = True
newmax = Hsky
tst_max = sky[Hsky] <= cut2
done = (Hsky == nlast ) * tst_max
if not done:
done = (tst_max) and (sky[Hsky+1 < nlast] > cut2)
istep = -1 + 2*tst_max
while not done:
newmax += istep
done = (newmax == nlast) or (newmax == -1)
if not done:
done = (sky[newmax] <= cut2) * (sky[newmax+1] >= cut2)
if tst_max:
delta = sky[Hsky+1:newmax] = skymid
else:
delta = sky[newmax+1:Hsky] - skymid
Totdelta = Totdelta + istep * bt.nansum(delta)
Totdeltasq = Totdeltasq + istep * bt.nansum(delta**2)
Hsky = newmax
redo = True
nsky = Hsky-Lsky
if (nsky < minsky):
print 'Error - Outlier rejection left too few sky elements'
return skymod, sigma
skymn = Totdelta/nsky
sigma = sqrt(Totdeltasq / nsky - skymn**2)
skymn = skymn+skymid
### Calculating a robust median by averaging the central 20% pixels
center = (Hsky+1 + Lsky)/2.
side = np.round(0.2*(Hsky-Lsky))/2. + 0.25
J = np.round(center-side)
K = np.round(center+side)
skymed = bt.nansum(sky[J:K])/(K-J+1)
if skymed < skymn:
dmod = 3.*skymed - 2*skymn-skymod
else:
dmod = skymn - skymod
if dmod*old < 0:
clamp = 0.5*clamp
skymod = skymod + clamp*dmod
old = dmod
return skymod, sigma
def makeMask(skyrad, xc, yc):
from math import floor
import numpy as np
lx = floor(xc - skyrad)
ly = floor(yc - skyrad)
dx = xc - lx
dy = yc - ly
x = np.arange(0, 2*skyrad)
xsq, ysq = np.meshgrid(x, x)
mask = (xsq-dx)**2 + (ysq-dy)**2
return np.sqrt(mask)-0.5
[docs]def Photcore(im, xc, yc, apr=4, skyrad=[5,7], phpadu=1.):
''' Photcore contains just the core of the aperture photometric routine.
It is design to work only with one star on one image
Input:
- im, the 2d image
- xc, yc, 2 scalars giving the coordinates of the star. No verification
are done to check whether the star is inside the image or not.
It is assumed that the star has to be centered or well inside the
image.
Keywords:
- apr, the aperture radius, default = 4
- skyrad, 2 values defining the radius of the sky annulus.
Default = [5,7]
- phpadu, the photon to analog constant. Default=1
Output:
- a vector containing, the flux, the flux errors, the sky
and the sky errors.
'''
import numpy as np
xm, ym = np.copy(xc), np.copy(yc)
box = im[ym-skyrad[1]:ym+skyrad[1], xm-skyrad[1]:xm+skyrad[1]]
area = apr**2 *np.pi
xsq = (np.arange(-(skyrad[1]), skyrad[1],1))**2
mask = np.sqrt(np.array([xsq + (j)**2 for j in np.arange(-(skyrad[1]),skyrad[1])]))
innermask = np.copy(mask)
indexap = (innermask < apr)
inap = box[indexap]
rap = innermask[indexap]
fractn = (apr-rap) # fraction of pixel to count
fullpix = (fractn >= 1.).nonzero()[0]
fractn[fullpix] = 1.0
fcounts = fullpix.size
fracpix = (fractn < 1.).nonzero()[0]
factor = (area - fcounts)/ sum(fractn[fracpix])
fractn[fracpix] = fractn[fracpix]*factor
apmag = sum(inap*fractn)
skypix = (mask >= skyrad[0])*(mask <= skyrad[1])
skyvalues = box[skypix]
# Rough estimate of the sky mode
skymed = np.median(skyvalues) # Sky median
skymean = skyvalues.mean()
skymod = 3*skymed - 2*skymean # Evaluation of the mode... under the assumption
# of slight sky contamination
skyvar = skyvalues.var() # Sky variance
sigsq = skyvar/len(skyvalues) # Square of
apmag = apmag - skymod*area
error1 = area*sigsq**2
error2 = apmag / phpadu
error3 = skyvar*area**2
magerr = np.sqrt(error1 + error2+error3)
return np.hstack((apmag, magerr, skymed, np.sqrt(skyvar)))
