Local coordinate aperture photometry for input to DSLR data worksheet

Definitions

Imports


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

Functions

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)
            instr_mags.append(mag)
        
        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], 
                          exact=True)

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

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', 
              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)

Inputs

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),
 "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 **
})

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.
gain=1.67

Outputs

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)
instr_mag_df.T.to_csv(instr_mag_csv_path)

# Display photometry table
instr_mag_df.T

In [ ]: