Aim: Demonstrate photometry on a series of bias and flat field corrected images of a Near Earth Asteroid.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import astropy.io.fits as fits
## make matplotlib appear in the notebook rather than in a new window
%matplotlib inline
In [2]:
datadir = ''
objname = '2016HO3'
In [3]:
def plotfits(imno):
img = fits.open(datadir+objname+'_{0:02d}.fits'.format(imno))[0].data
f = plt.figure(figsize=(10,12))
#im = plt.imshow(img, cmap='hot')
im = plt.imshow(img[480:580, 460:600], cmap='hot')
plt.clim(1800, 2800)
plt.colorbar(im, fraction=0.034, pad=0.04)
plt.savefig("figure{0}.png".format(imno))
plt.show()
In [4]:
numb = 1
plotfits(numb)
In [5]:
numb = 2
plotfits(numb)
Select part of the image for ease of display
In [6]:
partimg = fits.open(datadir+objname+'_01.fits')[0].data[480:580, 460:600]
Define starting values
In [7]:
targcen = np.array([22,42]) ## target center
compcen = np.array([75,125]) ## comparison center
Aperture photometry set up
In [8]:
searchr = 6 ## search box size
ap_r = 2 ## aperture radius
sky_inner = 3
sky_outer = 5
Calculate Center of Mass (CoM) defined as: $\bar{x} = \frac{\sum A_i x_i}{\sum A_i }$, $\bar{y} = \frac{\sum A_i y_i}{\sum A_i }$.
In [9]:
def cent_weight(n):
"""
Assigns centroid weights
"""
wghts=np.zeros((n),np.float)
for i in range(n):
wghts[i]=float(i-n/2)+0.5
return wghts
def calc_CoM(psf, weights):
"""
Finds Center of Mass of image
"""
cent=np.zeros((2),np.float)
temp=sum(sum(psf) - min(sum(psf) ))
print(temp)
cent[1]=sum(( sum(psf) - min(sum(psf)) ) * weights)/temp
cent[0]=sum(( sum(psf.T) - min(sum(psf.T)) ) *weights)/temp
return cent
Use centroiding algorithm to find the actual centers of the targe and comparison.
In [10]:
## Cut a box between search limits, centered around targcen
targbox = partimg[targcen[0]-searchr : targcen[0]+searchr, targcen[1]-searchr : targcen[1]+searchr]
weights = cent_weight(len(targbox))
tcenoffset = calc_CoM(targbox, weights)
print(tcenoffset)
tcenter = targcen + tcenoffset
In [11]:
plt.plot(sum(targbox))
plt.show()
In [12]:
compbox = partimg[compcen[0]-searchr : compcen[0]+searchr, compcen[1]-searchr : compcen[1]+searchr]
compw = cent_weight(len(compbox))
ccenoffset = calc_CoM(compbox,compw)
ccenter = compcen + ccenoffset
In [13]:
print(tcenter)
In [14]:
def circle(npix, r1):
"""
Builds a circle
"""
pup=np.zeros((npix,npix),np.int)
for i in range(npix):
for j in range(npix):
r=np.sqrt((float(i-npix/2)+0.5)**2+(float(j-npix/2)+0.5)**2)
if r<=r1:
pup[i,j]=1
return pup
In [15]:
def annulus(npix, r_inner,r_outer=-1.):
"""
Builds an annulus
"""
pup=np.zeros((npix,npix),np.int)
for i in range(npix):
for j in range(npix):
r=np.sqrt((float(i-npix/2)+0.5)**2+(float(j-npix/2)+0.5)**2)
if ((r<=r_outer)&(r>=r_inner)):
pup[i,j]=1
return pup
Create mask
In [16]:
circmask = circle(searchr*2, ap_r)
annmask = annulus(searchr*2, sky_inner, sky_outer)
Define new regions where the target and comparison are centered.
In [17]:
newtarg = partimg[int(round(tcenter[0]))-searchr : int(round(tcenter[0]))+searchr, int(round(tcenter[1]))-searchr : int(round(tcenter[1]))+searchr]
newcomp = partimg[int(round(ccenter[0]))-searchr : int(round(ccenter[0]))+searchr, int(round(ccenter[1]))-searchr : int(round(ccenter[1]))+searchr]
Place mask on region
In [18]:
targaper = newtarg * circmask
compaper = newcomp * circmask
Place mask on sky annulus slice.
In [19]:
targann = newtarg * annmask
compann = newcomp * annmask
a. Display image with target and comparison centers before and after centroiding
In [20]:
im = plt.imshow(partimg, cmap='hot')
plt.clim(1800, 2800)
plt.scatter(targcen[1], targcen[0], c='g', marker='x')
plt.scatter(compcen[1], compcen[0], c='g', marker='x')
plt.scatter(tcenter[1], tcenter[0], c='b', marker='x')
plt.scatter(ccenter[1], ccenter[0], c='b', marker='x')
plt.show()
b. Disply image with aperture mask and sky annulus
In [21]:
im = plt.imshow(targaper, cmap='hot')
plt.clim(1800, 2800)
plt.show()
In [22]:
im = plt.imshow(targann, cmap='hot')
plt.clim(1800, 2800)
plt.show()
Calculate Signal-to-Noise Ratio. CCD noise = sqrt(signal + background + dark current + read noise)
In [23]:
def calcsnr(target, bg):
signal = target - bg
noise = np.sqrt(signal + bg)
snr = signal / noise
return snr, noise
Sum all flux inside target and comparison apertures and divide by number of pixels to get average count per pixel.
In [24]:
targc = np.sum(targaper) / np.sum(circmask)
targbg= np.sum(targann) / np.sum(annmask)
compc = np.sum(compaper) / np.sum(circmask)
compbg= np.sum(compann) / np.sum(annmask)
In [25]:
snr, noise = calcsnr(targc, targbg)
print(snr)
In [26]:
snr, noise = calcsnr(compc, compbg)
print(snr)
In [27]:
#Try a range of aperture sizes
apertures = np.arange(1, 10, 1)
snrold = 0
for aper in apertures:
apermask = circle(searchr*2, aper)
targc = np.sum(apermask*newtarg) / np.sum(apermask)
snr, noise = calcsnr(targc, targbg)
if snr > snrold:
bestaper = aper
snrold = snr
snr, noise = calcsnr(targc, targbg)
In [28]:
targc = circle(searchr*2, ap_r)*newtarg
targskyc = annulus(searchr*2, sky_inner, sky_outer)*newtarg
compc = circle(searchr*2, ap_r)*newcomp
compskyc = annulus(searchr*2, sky_inner, sky_outer)*newcomp
ratio = np.sum(compc)/np.sum(targc)
sigmaratio = ratio*np.sqrt((np.sum(targc)/np.sum(targskyc))**2 + (np.sum(compc)/np.sum(compc))**2)
deltamag = -2.5*np.log10(ratio)
sigmamag = 2.5*sigmaratio/(ratio*np.log(10))
refmag = 19.4
mag = refmag - deltamag
print("Measured Magnitude = {:0.3f} ± {:0.3f}".format(mag, sigmamag))
In [ ]: