Aperture and simple PSF-fitting photometrymem

Contents:

Python script for various photometry tasks. See especially aperphot(), for basic aperture photometry. If you need something fancier, try PyRAF, DAOPHOT, etc.

phot.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.

INPUTS:

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.

retfull:

Also return arrays of target mask, sky mask, and frame used. This option is a memory hog!

OUTPUTS:

phot object.

EXAMPLE:

import astropy.io.fits
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 = astropy.io.fits.open(img)
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:   ", observation.bg

REQUIREMENTS:

scipy.interpolate, pyfits, numpy...

phot.centroid(im, mask=None, w=None, x=None, y=None)

Compute the centroid of an image with a specified binary mask projected upon it.

INPUT:

im – image array mask – binary mask, 0 in ignored regions and 1 in desired regions w is typically 1.0/u**2, where u is the uncertainty on im x,y are those generated by meshgrid.

OUTPUT:

(x0,y0) tuple of centroid location

phot.dms(d, delim=':')

Convert degrees, minutes, seconds to decimal degrees, and back.

EXAMPLES:

dms(‘150:15:32.8’) dms([7, 49]) dms(18.235097)

Also works for negative values.

SEE ALSO: hms

phot.egauss2d(param, x, y, z, ez=None)

Return the chi-squared error on a 2D gaussian fit. See gauss2d for more details on the input parameters.

z is the array of data to be fit ez is an optional array of one-sigma errors to the data in z.

phot.estbg(im, mask=None, bins=None, plotalot=False, rout=(3, 200), badval=nan, verbose=False)

Estimate the background value of a masked image via histogram fitting.

INPUTS:

im – numpy array. Input image.

OPTIONAL INPUTS:

mask – numpy array. logical mask, False/0 in regions to ignore bins – sequence. edges of bins to pass to HIST plotalot – bool. Plot the histogram and fit. rout – 2-tuple of (nsigma, niter) for analysis.removeoutliers.

Set to (Inf, 0) to not cut any outliers.

badval – value returned when things go wrong.

OUTPUT:

b, s_b – tuple of (background, error on background) from gaussian fit.

Note that the error is analagous to the standard deviation on the mean

COMMENTS:

The fit parameters appear to be robust across a fairly wide range of bin sizes.

phot.gauss2d(param, x, y)

Compute a gaussian distribution at the points x, y.

p is a three- or four-component array, list, or tuple:

z = [p4 +] p0/(p1*sqrt(4pi)) * exp(-r**2 / (2*p1**2))

where r = sqrt((x-p2)**2 + (y-p3)**2)

p[0] – Volume under the gaussian p[1] – one-sigma dispersion p[2] – X- offset of center p[3] – Y- offset of center p[4] – optional constant, vertical offset

x & y must be equal-sized arrays from numpy.meshgrid

NOTE: FWHM = 2*sqrt(2*ln(2)) * p1 ~ 2.3548*p1

SEE ALSO: egauss2d, numpy.meshgrid

phot.hess(color, mag, binsize, **kw)

Compute a hess diagram (surface-density CMD) on photometry data.

INPUT:

color mag binsize – width of bins, in magnitudes

OPTIONAL INPUT:

cbin= – set the centers of the color bins mbin= – set the centers of the magnitude bins

OUTPUT:

A 3-tuple consisting of:

Cbin – the centers of the color bins Mbin – the centers of the magnitude bins Hess – The Hess diagram array

phot.hms(d, delim=':')

Convert hours, minutes, seconds to decimal degrees, and back.

EXAMPLES:

hms(‘15:15:32.8’) hms([7, 49]) hms(18.235097)

Also works for negative values.

SEE ALSO: dms

phot.lmcextinct(ra, dec, **kw)

Use the Zaritsky & Harris (ZH) map to get A_V extinction in the LMC.

INPUT:

ra – decimal degrees of right ascension dec – decimal degrees of declination

OPTIONAL INPUT:

method=’griddata’ /’nearest’ – interpolation method map=’both’ /’hot’/’cool’ – Which ZH map to use hot=’lmc_hotav.fits’ – filename of hot map cool=’lmc_coolav.fits’ – filename of cool map null=0.0 – “no data” value in map verbose=False / True – print out some spam and eggs

EXAMPLE:

ra = hms(‘05:51:56.6’) dec = dms(‘-66:25:08.5’) lmcextinct(ra, dec)

If their map returns null, interpolate from the surroundings.

Note that these calculations are definitely _not_ optimized.

SEE ALSO:

hms, dms

phot.makemask(x, y, params, shape='circle')

Generate a binary (logical) mask with given shape and location.

INPUTS:

x = x-coodinate system (made with meshgrid) y = y-coodinate system (made with meshgrid)

params:

shape=’circle’:

params(1) = x-offset params(2) = y-offset params(3) = x-diameter params(4) = OPTIONAL y-diameter

shape=’quad’:

params: list of quadrants to include in the mask. The

midpoint for determining quadrants will be mid = (xmax+xmin)/2. Quadrants are: 0: x<midX and y<midY 1: x<midX and y>=midY 2: x>=midX and y<midY 3: x>=midX and y>=midY

OPTIONAL INPUTS:

shape=: desired mask shape. Currently only ‘circle’ is valid.

OUTPUTS:

mask = NxM grided representation of specified mask

where NxM are the size of the x,y input meshgrids

EXAMPLE:

x=arange(32); y=arange(64) xx,yy = meshgrid(x,y) m = makemask(xx,yy,(24, 53, 10)) m[53,24] # —-> True

phot.psffit(psf, frame, loc=None, w=None, scale=100, dframe=9, xoffs=None, yoffs=None, verbose=False)

Conduct simple PRF model-fitting at a specified grid of positions

INPUTS:

psf – model PSF (or PRF), supersampled by factor ‘scale’

frame – science frame to which psf will be fit. best if it’s odd-sized

OPTIONS:

loc – (x,y) integers, star location in data frame (e.g., data[x,y])

w – [array] weights of pixels in data frame, (typ. 1/sigma^2)

scale – supersampling level of PSF

dframe – [odd int] diameter of square box around target

location. Errors may occur if this size is larger than your binned-down PRF model.

xoffs – [int array] subpixel x offsets to test.

yoffs – [int array] subpixel y offsets to test.

RETURNS:

modelpsf, data, chisq, background, fluxscale, xoffset, yoffset, chisq[ii,jj],background[ii,jj], fluxscale[ii,jj]

EXAMPLE:

import k2
import phot
import pylab as py

fn = 'kplr060018142-2014044044430_lpd-targ.fits'
f = pyfits.open('kplr060018142-2014044044430_lpd-targ.fits')
image = f[1].data['flux'][0]
prf = k2.loadPRF(fn)
out = phot.psffit(prf, image, (33, 23), scale=50, dframe=7, verbose=True)

py.figure()
py.subplot(132)
py.imshow(out[0])
py.title('Best-fit Model PRF')
py.colorbar()
py.subplot(131)
py.imshow(out[1])
py.title('Observed Data')
py.clim([out[0].min(), out[0].max()])
py.colorbar()
py.subplot(133)
py.imshow(out[1] - out[0])
py.title('Data - Model')
py.colorbar()

phot.psffiterr(xyoffset, psf, frame, w=None, scale=100, dframe=9, loc=None, verbose=False)

USAGE:

fit = optimize.fmin(psffiterr, [0,0], args=(psf, frame,weights,100,13, loc),xtol=0.5,ftol=1, full_output=True)

modelout = psffit(psf, frame, loc, weights, scale=100, dframe=13, xoffs=[fit[0][0]], yoffs=[fit[0][1]])

phot.readogle(filename, **kw)

Read in an OGLE-II photometry map file and output the data.

Returns a 6-tuple with each element an array:

0 – RA. Strings of Right Ascension values. 1 – DEC. Strings of Declination values 2 – x_ref. OGLE x-coordinate, pixels 3 – y_ref. OGLE y-coordinate, pixels 4 – v_mag. OGLE V magnitude 5 – i_mag. OGLE I magnitude

If you like, use HMS and DMS to convert the RA and DEC values returned.

phot.sens(texp, mag, lin=380000.0, tdead=15.0, rn=8.0, csky=12.0, zp=22.8, z=1.0, h=1280.0, d=1.0)

Calculate sensitiviy for various exposure times.

texp – 1D array of times mag – target magnitude npix – number of pixels read out

lin – linearity limit of the detector (ADU) tdead – detector dead time (reset + readout) between integrations rn – read noise, in ADU csky – sky background per second, in ADU zp – photometric zero point (magnitude for 1 ADU/sec) z – airmass h – observatory altitude, in meters d – telescope mirror diameter, in meters

TBD: Flat-field noise!

For now, I impose an artificial floor to the number of pixels needed – it won’t spread it out over less than 8.

phot.snr(mag=20, itime=1.0, read=24.5, sky=8.43, npix=24.0, zero=26.44, dark=0.0)

Calculate SNR of a photometric CCD observation of given parameters.

The default keyword parameters are for the CTIO Blanco MOSAIC imager, operating in V band.

phot.subreg(fn, center=None, dim=None, verbose=False)

Load a subsection of a list of 2D FITS files.

INPUTS:

fn – (str) filename or (list) list of filenames to load center – (2-tuple) (x0,y0) center of region to load. dim – (2-tuple) (dx,dy) sizes of region to load.

OUTPUT:

region – (array) (N x dx x dy)-sized stack of subregions

NOTE: center can also be a list of tuples (1 per file), but dim cannot.

phot.subreg2(fn, center=None, dim=None, verbose=False)

Load a subsection of a list of 2D FITS files.

INPUTS:

fn – str, list of str, or 2- or 3D Numpy array.

1. filename to load, OR
2. list of filenames to load, OR
3. Numpy array of data, already loaded.

center – (2-tuple) (x0,y0) center of region to load.

dim – (2-tuple) (dx,dy) sizes of region to load.

OUTPUT:

region – (array)

(N x dx x dy)-sized stack of subregions. Note that this will always be 3D!

NOTE: center can also be a list of tuples (1 per file), but dim cannot.