In [ ]:
import os
from random import random
from collections import OrderedDict
import numpy as np
import pandas as pd
from astropy.io import fits
from astropy.visualization import astropy_mpl_style
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from matplotlib.patches import Circle
from matplotlib.offsetbox import (TextArea, DrawingArea, OffsetImage, AnnotationBbox)
plt.style.use(astropy_mpl_style)
%matplotlib inline
from PythonPhot import aper
In [ ]:
def multi_file_photometry(fits_root, fits_files, data_index, coords, dataframe,
aperture_radius, inner_sky_radius, outer_sky_radius,
gain=1, zeropoint=0, suffix='.fit'):
for fits_file in fits_files:
fits_file_path = os.path.join(fits_root, fits_file)
hdus = fits.open(fits_file_path)
instr_mags = []
for x, y in coords:
time, mag = aperture_photometry(hdus[data_index], x, y,
aperture_radius, inner_sky_radius, outer_sky_radius,
gain=gain, zeropoint=zeropoint)
instr_mags.append(mag)
dataframe[fits_file[0:fits_file.rindex(suffix)]] = [time] + instr_mags
In [ ]:
def aperture_photometry(hdu, x, y,
aperture_radius, inner_sky_radius, outer_sky_radius,
gain=1, zeropoint=0):
image_data = hdu.data
time = hdu.header[time_name]
mag, magerr, flux, fluxerr, sky, skyerr, badflag, outstr = \
aper.aper(image_data, x, y, phpadu=gain,
apr=aperture_radius, zeropoint=zeropoint,
skyrad=[inner_sky_radius, outer_sky_radius],
exact=True)
return time, mag[0]
In [ ]:
def show_image(image_data, coord_map, aperture_size, annotate=True, vmin=10, vmax=200, figx=20, figy=10):
fig = plt.figure(figsize=(figx, figy))
plt.imshow(image_data, cmap='gray', vmin=vmin, vmax=vmax)
plt.gca().invert_yaxis()
plt.colorbar()
if annotate:
for designation in coord_map:
xy = coord_map[designation]
annotate_image(fig.axes[0], designation, xy, aperture_size)
plt.show()
In [ ]:
def annotate_image(axis, designation, xy, aperture_size):
axis.plot(xy[0], xy[1], 'o', markersize=aperture_size,
markeredgecolor='r', markerfacecolor='none',
markeredgewidth=2)
offsetbox = TextArea(designation, minimumdescent=False)
ab = AnnotationBbox(offsetbox, xy,
xybox=(-20, 40+random()*10-10),
xycoords='data',
boxcoords="offset points",
arrowprops=dict(arrowstyle="->"))
axis.add_artist(ab)
In [ ]:
# Instrumental magnitude output file path
instr_mag_file_root = "/Users/david/_photometry_working"
instr_mag_csv_file = "instr_mags.csv"
# FITS file directory
fits_root = "/Users/david/_photometry_working"
# B, G, and R FITS file prefixes to identify files,
# e.g. stk-median-g matches stk-median-g1.fit, stk-median-g2.fit, ...
fits_prefixes = ["stk-median-b", "stk-median-g", "stk-median-r"]
# FITS file data HDU index
data_index = 0
# Time column name
time_name = "JD"
In [ ]:
# Ordered dictionary of object names/IDs to local coordinates
position_map = OrderedDict({
# ** Your name-tuple pairs go here **
"eta Car":(2374.478,814.152),
"45 (CHK)":(2352.016,1629.581),
"46":(2585.264,600.724),
"47":(2512.030,1096.505),
"51":(2664.810,1040.385),
"52 795":(1657.205,1312.103),
"52 708":(1777.295,1584.834),
"55":(2482.270,1252.879)
# ** END **
})
In [ ]:
# Aperture radii
measurement_aperture = 13
inner_sky_annulus = 16
outer_sky_annulus = 21
# ph/ADU
# Note: PythonPhot's aperture photometry function takes a phadu parameter.
# Assumption: this is photons/ADU or e-/ADU, i.e. gain.
gain=1.67
In [ ]:
files = os.listdir(fits_root)
fits_files = []
for fits_prefix in fits_prefixes:
fits_files += sorted([file for file in files if fits_prefix in file])
In [ ]:
fits_file = fits_files[0]
hdus = fits.open(os.path.join(fits_root, fits_file))
image_data = hdus[data_index].data
median = np.median(image_data)
show_image(image_data, position_map, measurement_aperture, annotate=True, vmin=10, vmax=median*4)
In [ ]:
# Create empty table with time and object headers
pd.options.display.float_format = '{:,.6f}'.format
instr_mag_df = pd.DataFrame()
names = [name for name in position_map]
instr_mag_df['name'] = [time_name] + names
instr_mag_df.set_index('name', inplace=True)
In [ ]:
# Carry out photometry on B, G, R FITS files, yielding instrumental magnitudes
positions = position_map.values()
multi_file_photometry(fits_root, fits_files, data_index, positions, instr_mag_df,
measurement_aperture, inner_sky_annulus, outer_sky_annulus, gain=gain)
In [ ]:
# Save photometry table as CSV
instr_mag_csv_path = os.path.join(instr_mag_file_root, instr_mag_csv_file)
instr_mag_df.T.to_csv(instr_mag_csv_path)
# Display photometry table
instr_mag_df.T
In [ ]: