Commit fcf59a4b authored by Maxime Charpentier's avatar Maxime Charpentier

Add new file

def aperphot(fn, timekey=None, pos=[0,0], dap=[2,4,6], mask=None, verbose=False, nanval=999, resamp=None, retfull=False, ignorenan=True):
"""Do aperture photometry on a specified file.
pos : 2-sequence
center of apertures (as if indexing: fn[pos[0], pos[1]])
dap : 3-sequence
Photometry aperture DIAMETERS:
-- target aperture (within which to sum flux)
-- inner sky aperture
-- outer sky aperture
resamp : int
Factor by which to interpolate frame before measuring photometry
(in essence, does partial-pixel aperture photometry)
Aperture masking:
If no mask is passed in, use the star's input position and
aperture diameters to create binary pixel masks for stellar and
background photometry. If a mask is passed in, stellar and
background aperture masks are generated from where all input
mask elements equal 1 and 2, respectively.
Also return arrays of target mask, sky mask, and frame used.
This option is a memory hog!
:class:`phot` object.
from astropy import wcs
import numpy as np
from phot import aperphot
img='62_z_CDFs_goods_stamp_img.fits' #path to the image
RA = 52.9898239
DEC = -27.7143114
hdulist =
w = wcs.WCS(hdulist['PRIMARY'].header)
world = np.array([[RA, DEC]])
pix = w.wcs_world2pix(world,1) # Pixel coordinates of (RA, DEC)
print "Pixel Coordinates: ", pix[0,0], pix[0,1]
#call aperture function
observation=aperphot(img, timekey=None, pos=[pix[0,0], pix[0,1]], dap=[4,8,12], resamp=2, retfull=False)
# Print outputs
print "Aperture flux:", observation.phot
print "Background: ",
scipy.interpolate, pyfits, numpy...
# 2009-09-14 10:49 IJC: Created
# 2010-01-15 14:20 IJC: Added numpy "_string" check
# 2011-12-29 12:01 IJMC: Added peak pixel values to photometry report.
# 2012-01-25 11:26 IJMC: Adding "resamp" option -- thanks to
# K. Stevenson and J. Harrington of UCF for
# the suggestion.
# 2012-02-26 11:53 IJMC: Now return 'ntarg' and 'nsky' -- number of pixels used.
# 2012-06-07 08:27 IJMC: 'peak' values are now corrected for the
# resampling factor.
# 2012-07-03 10:35 IJMC: Fixed a key bug: frames were not
# correctly background-subtracted when
# applying partial-pixel resampling.
# 2012-10-19 13:41 IJMC: Documented 'retfull' option; changed default.
# 2013-03-20 09:21 IJMC: More error-checking for saving header
# keywords. Thanks to A. Weigel @
# ETH-Zurich for catching this!
# 2014-05-07 16:05: added "ephot" by J. Xavier Prochaska, UCSC
# 2014-08-28 09:32 IJMC: Added better checking for user-input
# masks and resamp>1. Corrected x/y
# indexing for non-square inputs arrays.
# 2014-09-06 14:41 IJMC: Added better checking for dap & pos.
# 2014-09-08 10:10 IJMC: Fix masking & nan-handling.
from numpy import meshgrid, median,isfinite,sort,ndarray,string_
import numpy as np
from import fits as pyfits
import pyfits
from analysis import fixval
from os import path
from scipy import interpolate
thisobs = phot()
if pos is None:
pos = 0, 0
if dap is None:
dap = 1, 2, 3
x0, y0 = pos
dap_targ, dap_skyinner, dap_skyouter = dap
if resamp is None or resamp<1:
resamp = 1
resamp = float(resamp)
# Determine size:
if isinstance(fn,str):
nx = pyfits.getval(fn, 'NAXIS2')
ny = pyfits.getval(fn, 'NAXIS1')
elif isinstance(fn,ndarray):
nx,ny = fn.shape
nx0, ny0 = nx, ny
nx = ((nx - 1)*resamp + 1.) # Avoid resampling at pixel locations
ny = ((ny - 1)*resamp + 1.) # outside the original boundaries.
# Generate or load masks:
if mask==None:
xx,yy = meshgrid(np.arange(ny)/resamp, np.arange(nx)/resamp)
mask_targ = makemask(xx, yy, (y0, x0, dap_targ))
mask_s1 = makemask(xx, yy, (y0,x0, dap_skyinner))
mask_s2 = makemask(xx, yy, (y0,x0, dap_skyouter))
mask_sky = mask_s2 - mask_s1
userInputMask = False
mask_targ = mask==1
mask_sky = mask==2
userInputMask = True
# Load data frame:
thisobs = phot()
if pos is None:
pos = 0, 0
if dap is None:
dap = 1, 2, 3
if isinstance(fn,ndarray):
frame = fn
elif isinstance(fn, str) or isinstance(fn,string_):
if not path.isfile(fn):
print "file %s not found! exiting..." % fn
return thisobs
frame = pyfits.getdata(fn)
fixval(frame, nanval)
# Resample data frame
badval = -12345.6789
if resamp>1:
frame0 = frame.copy()
mask_good0 = np.isfinite(frame0)
frame0[True-mask_good0] = badval
xx0 = range(nx0)
yy0 = range(ny0)
x1,y1 = np.arange(nx)/resamp, np.arange(ny)/resamp
rectspline = interpolate.fitpack2.RectBivariateSpline(xx0, yy0, frame0, kx=1, ky=1, s=0)
frame = rectspline(x1, y1)
if userInputMask and \
((frame.shape<>mask_targ.shape) or (frame.shape<>mask_sky.shape)):
print "User-entered mask did not have the correct size for the resampling"
print " input used (resamp=%1.3f). Beware!" % resamp
if ignorenan:
mask_good = np.isfinite(frame)
mask_targ = mask_targ * mask_good
mask_sky = mask_sky * mask_good
#from pylab import *
# Measure background and aperture photometry
thisbg, thisebg = estbg(frame, mask=mask_sky, plotalot=verbose, rout=[3,99], verbose=verbose)
if not np.isfinite(thisbg):
thisbg = 0.
thisphot = ((frame - thisbg)[mask_targ]).sum() /resamp/resamp
peak = frame.max()
peak_targ = frame[mask_targ].max()
peak_annulus = frame[mask_sky].max()
thisobs.bgstr='phot.estbg: SDOM on bg histogram mean & dispersion after outlier rejection'
thisobs.photstr='by-hand background-subtracted aperture photometry'
thisobs.ntarg = mask_targ.sum()/resamp/resamp
thisobs.nsky = mask_sky.sum()/resamp/resamp
thisobs.peak = peak
thisobs.peak_targ = peak_targ
thisobs.peak_annulus = peak_annulus
thisobs.peakstr = 'peak pixel value in frame'
thisobs.peak_targstr = 'peak pixel value in target aperture'
thisobs.peak_annulusstr = 'peak pixel value in sky annulus'
thisobs.position = pos
thisobs.positionstr = 'user-specified, zero-indexed pixel coordinates.'
if isinstance(fn, str):
header = pyfits.getheader(fn)
if not timekey==None:
if timekey in header:
thisobs.timestr='heliocentric modified julian date'
if 'object' in header: thisobs.object = header['object']
if 'exptime' in header: thisobs.exptime = header['exptime']
thisobs.aper = dap
thisobs.aperstr = 'target, inner, outer aperture diameters, in pixels.'
thisobs.resamp = resamp
# Simple stats :: JXP 2014 May 6
var = thisphot + np.sqrt(thisobs.nsky)*
thisobs.ephot = np.sqrt(var)
if retfull:
thisobs.mask_targ = mask_targ
thisobs.mask_sky = mask_sky
thisobs.frame = frame
if verbose>0:
from pylab import figure, colorbar
from nsdata import imshow
figure(); imshow(frame*mask_targ); colorbar()
figure(); imshow(frame*mask_sky); colorbar()
return thisobs
\ No newline at end of file
Markdown is supported
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment