Photometry using astropy photutils and AAVSO comparison stars

1 Import dependencies and setup matplotlib

In [1]:
import requests, math, glob
import pandas as pd
import numpy as np
from photutils import DAOStarFinder
from astropy.stats import mad_std
from astropy.io import fits
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
import astropy.units as u
import matplotlib.pyplot as plt
from photutils import aperture_photometry, CircularAperture
from astroquery.simbad import Simbad
from datetime import datetime
import warnings
warnings.filterwarnings('ignore')

%matplotlib inline
plt.style.use('seaborn')
pd.options.mode.chained_assignment = None 

2 Define input file, star name and comparason star range¶

In [2]:
FITS_FILE = '/home/dokeeffe/Pictures/CalibratedLight/2017-10-02/KIC08462852-2017-10-02-20-56-22Light_005.fits'
STAR_NAME = 'KIC08462852'
BRIGHTEST_COMPARISON_STAR_MAG = 11.0
DIMMEST_COMPARISON_STAR_MAG = 13.0

3 Define a function to download comparison stars from AAVSO.

This will return a tuple containing an array and the chart ID.

In [3]:
def get_comp_stars(ra,dec,filter_band='V',field_of_view=18.5):
    result = []
    vsp_template = 'https://www.aavso.org/apps/vsp/api/chart/?format=json&fov={}&maglimit=18.5&ra={}&dec={}'
    print(vsp_template.format(field_of_view, ra, dec))
    r = requests.get(vsp_template.format(field_of_view, ra, dec))
    chart_id = r.json()['chartid']
    print('Downloaded Comparison Star Chart ID {}'.format(chart_id))
    for star in r.json()['photometry']:
        comparison = {}
        comparison['auid'] = star['auid']
        comparison['ra'] = star['ra']
        comparison['dec'] = star['dec']
        for band in star['bands']:
            if band['band'] == filter_band:
                comparison['vmag'] = band['mag']
                comparison['error'] = band['error']
        result.append(comparison)
    return result, chart_id

4 Download comparison stars and search simbad for our target.

Use astroquery to locate the RA/DEC of our target star.

Here we also define results which will be a collection of data progressivly enriched as we go

In [4]:
astroquery_results = Simbad.query_object(STAR_NAME)
TARGET_RA = str(astroquery_results[0]['RA'])
TARGET_DEC = str(astroquery_results[0]['DEC']).replace('+','').replace('-','')
results, chart_id = get_comp_stars(TARGET_RA, TARGET_DEC)
print('{} comp stars found'.format(len(results)))
results.append({'auid': 'target', 'ra': TARGET_RA, 'dec': TARGET_DEC})

results
https://www.aavso.org/apps/vsp/api/chart/?format=json&fov=18.5&maglimit=18.5&ra=20 06 15.4553&dec=44 27 24.793
Downloaded Comparison Star Chart ID X21328CWX
5 comp stars found
Out[4]:
[{'auid': '000-BLS-551',
  'dec': '44:22:48.0',
  'error': 0.054,
  'ra': '20:06:48.09',
  'vmag': 11.263},
 {'auid': '000-BML-045',
  'dec': '44:27:05.0',
  'error': 0.036,
  'ra': '20:06:36.68',
  'vmag': 12.415},
 {'auid': '000-BLS-555',
  'dec': '44:30:54.0',
  'error': 0.05,
  'ra': '20:06:21.25',
  'vmag': 12.789},
 {'auid': '000-BML-046',
  'dec': '44:35:41.8',
  'error': 0.048,
  'ra': '20:06:53.14',
  'vmag': 12.983},
 {'auid': '000-BML-047',
  'dec': '44:25:53.7',
  'error': 0.038,
  'ra': '20:06:00.39',
  'vmag': 13.327},
 {'auid': 'target', 'dec': '44 27 24.793', 'ra': '20 06 15.4553'}]

5 Extract all sources from image

In [5]:
# extract sources from image and add details to comp_stars
fwhm = 3.0
source_snr = 20
hdulist = fits.open(FITS_FILE)
data = hdulist[0].data.astype(float)
header = hdulist[0].header
wcs = WCS(header)
bkg_sigma = mad_std(data)    
daofind = DAOStarFinder(fwhm=fwhm, threshold=source_snr*bkg_sigma)    
sources = daofind(data)

sources
Out[5]:
<Table length=115>
idxcentroidycentroidsharpnessroundness1roundness2npixskypeakfluxmag
int64float64float64float64float64float64float64float64float64float64float64
1645.556580297.699083270450.454248555190.1746450322460.35553948023925.00.070365.293778128.269908551-3.62831100903
21309.1659823111.59958439960.446728725220.2537016144780.49747590310225.00.07287.614090362.85692255025-1.13974616759
31375.6579289911.99959568790.4965712464110.1859798474860.52549303743725.00.03368.496644031.05632128534-0.0594900777583
41122.9478215116.71589949090.460102848306-0.08966978266260.4727670861925.00.08200.645219463.19746747319-1.26201533809
51362.464514326.96197756560.4338037914440.1129425017510.42973115343525.00.021372.89932418.4265184611-2.31412044147
6611.03721443328.78760458110.4016858425060.1679216624560.41934408149925.00.05478.051407261.98099565096-0.742208805242
7358.18949385734.46178149110.503866132614-0.009481791565240.2661832946525.00.090850.628901434.2481732847-3.83659353055
8449.53959859282.10691777620.3793093254030.09805237015240.29422574310825.00.07390.863115372.73814589578-1.09364146198
9964.42383054983.98941471410.3853663893080.1206510008410.34791377159125.00.03227.591622441.11273974981-0.115984005994
101492.3411787988.0765522970.4033273158030.05610768190010.48617835936525.00.04414.124170451.67920723014-0.562760738787
.................................
106827.0799326341179.020255250.397659638526-0.2582973218530.7089859902225.00.05463.5716491.83315370736-0.657997203678
10737.70682835361181.085616580.491163497203-0.9000168656090.59673712403625.00.04461.481389461.29027419392-0.276705028177
10836.69505325051181.684362430.438245603628-0.4695660883090.74391564656725.00.04282.613768421.46004410038-0.410914934433
1091614.103243821208.141460270.4474985680290.2999194280330.42518684544325.00.03437.808923851.01760469367-0.0189477536934
1101198.99602061224.783851650.5049921450830.09047486346330.64985304922325.00.04715.865641241.49606300697-0.437374710765
111823.2822151181226.317686570.472764055771-0.2600892230410.66320868119725.00.011893.65540794.12646557823-1.53894556676
112754.380683891230.616205760.459792598755-0.4044745360980.48456900518125.00.04524.829135371.32212938817-0.303184896983
11330.93644894761238.42374320.514502694717-0.4074840644840.77266308855525.00.06368.251121691.93519483193-0.716811738856
114152.1075442221241.393311040.439041849586-0.8243671613860.62668987288825.00.04610.070897411.35224097488-0.32763522889
115630.4345297091245.659728820.382036296985-0.3447071196570.69681239235825.00.03931.927006941.17175061378-0.172087974203

6 Find the sources that correspond to our target and comparison stars

Any source within 4 arc seconds of a comparison star is considered a match. Here the results will be converted to a pandas dataframe and will now contain the x,y coords of the target and comparison stars.

In [6]:
for star in results:
    star_coord = SkyCoord(star['ra'],star['dec'], unit=(u.hourangle, u.deg))
    xy = SkyCoord.to_pixel(star_coord, wcs=wcs, origin=1)
    x = xy[0].item(0)
    y = xy[1].item(0)
    for source in sources:
        if(source['xcentroid']-4 < x < source['xcentroid']+4) and source['ycentroid']-4 < y < source['ycentroid']+4:
            star['x'] = x
            star['y'] = y
            star['peak'] = source['peak']
results = pd.DataFrame(results)
results
Out[6]:
auid dec error peak ra vmag x y
0 000-BLS-551 44:22:48.0 0.054 19038.106969 20:06:48.09 11.263 388.415326 817.537270
1 000-BML-045 44:27:05.0 0.036 8153.566262 20:06:36.68 12.415 491.879800 617.615307
2 000-BLS-555 44:30:54.0 0.050 5919.023955 20:06:21.25 12.789 628.526243 440.515839
3 000-BML-046 44:35:41.8 0.048 4777.952321 20:06:53.14 12.983 365.886101 203.792448
4 000-BML-047 44:25:53.7 0.038 4538.973603 20:06:00.39 13.327 797.954483 683.615847
5 target 44 27 24.793 NaN 13457.943342 20 06 15.4553 NaN 672.395700 607.642881

7 Plot the image annotated with the comparsion stars

This is just to check everything looks ok

In [7]:
# plot the image with overlay
results = results.query('x > 0 and y > 0') 

hdulist = fits.open(FITS_FILE)
data = hdulist[0].data.astype(float)
fig=plt.figure(figsize=(18, 16), dpi= 80, facecolor='w', edgecolor='k')
positions = (results['x'], results['y'])    
position_apps = CircularAperture(positions, r=20.)    
# target_app = CircularAperture(target_xy, r=20.)    
plt.imshow(data, cmap='gray_r', origin='lower', vmin=0, vmax=2500)
position_apps.plot(color='red', lw=1.5, alpha=0.5)
# target_app.plot(color='blue', lw=1.5, alpha=0.5)
to_plot = results.query('peak < 50000 and vmag < 20 and vmag > 0') 
for to_annotate in results.iterrows():
    plt.annotate('{}'.format(to_annotate[1]['auid']),
        xy=(to_annotate[1]['x'], to_annotate[1]['y']), xycoords='data',
        xytext=(-150, 130), textcoords='offset points', size='16',
        arrowprops=dict(arrowstyle="->"))