Local coordinate aperture photometry for input to DSLR data worksheet



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)
%matplotlib inline

from PythonPhot import aper


Photometry of a list of FITS files, creating a table of times and instrumental magnitudes

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)
        dataframe[fits_file[0:fits_file.rindex(suffix)]] = [time] + instr_mags

Single image+coordinate photometry, returning a time and instrumental magnitude. Invoked by multi_file_photometry()

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], 

    return time, mag[0]

Display an image with target and reference stars annotated, to sanity check local coordinates

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)

    if annotate:
        for designation in coord_map:
            xy = coord_map[designation]
            annotate_image(fig.axes[0], designation, xy, aperture_size)

Annotate plot axis with coordinate positions and designations. Invoked by show_image()

In [ ]:
def annotate_image(axis, designation, xy, aperture_size):
    axis.plot(xy[0], xy[1], 'o', markersize=aperture_size, 
              markeredgecolor='r', markerfacecolor='none', 
    offsetbox = TextArea(designation, minimumdescent=False)

    ab = AnnotationBbox(offsetbox, xy,
                        xybox=(-20, 40+random()*10-10),
                        boxcoords="offset points",


Change these to suit your environment

File settings

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"

Map of object designations to local coordinates

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),
 "52 795":(1657.205,1312.103),
 "52 708":(1777.295,1584.834),
 # ** END **

Aperture radii and gain

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.


Find B, G, R files in the FITS file directory

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])

Aperture location sanity check by visual inspection

Arbitrarily choose the first G FITS 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)

Aperture photometry

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)

# Display photometry table

In [ ]: