

In [ ]:

    
from astropy import nddata
from astropy.modeling import models
import astropy.units as u
from astropy.io import fits
import matplotlib.pyplot as pl
%matplotlib inline
import numpy as np

import ccdproc



In [ ]:

    
nddata.conf.warn_unsupported_correlated = False



In [ ]:

    
path = "/Users/adrian/projects/triand-rrlyrae/data/mdm-fall-2015/n3/"
images = ccdproc.ImageFileCollection(path, keywords='*')

CCD Properties (Echelle at MDM)



In [ ]:

    
gain = 1.3 * u.electron / u.adu
readnoise = 7.9 * u.electron

Function to subtract overscan and trim a list of images



In [ ]:

    
def oscan_and_trim(image_list):
    """
    Remove overscan and trim a list of images. The original list is 
    replaced by a list of images with the changes applied.
    """
    for idx, img in enumerate(image_list):
        oscan = ccdproc.subtract_overscan(img, img[:, 306:3314], add_keyword={'oscan_sub': True, 'calstat': 'O'}, model=models.Polynomial1D(1))
        image_list[idx] = ccdproc.trim_image(oscan[:, :300], add_keyword={'trimmed': True, 'calstat': 'OT'})



In [ ]:

    
def med_over_images(masked_arr, axis=0):
    """
    Calculate median pixel value along specified axis
    
    Uses bottleneck.nanmedian for speed
    """
    
    dat = masked_arr.data.copy()
    dat[masked_arr.mask] = np.NaN
    return np.nanmedian(dat, axis=axis)

Make a master bias

First load all bias images



In [ ]:

    
bias_list = []
for hdu, fname in images.hdus(imagetyp='bias', return_fname=True):
    meta = hdu.header
    meta['filename'] = fname
    bias_list.append(ccdproc.CCDData(hdu.data*u.adu, meta=meta))



In [ ]:

    
print("{} bias frames".format(len(bias_list)))

Now subtract overscan and trim



In [ ]:

    
oscan_and_trim(bias_list)

Combine the bias frames by averaging



In [ ]:

    
biases = ccdproc.Combiner(bias_list)
master_bias = biases.median_combine()



In [ ]:

    
bias_mean = np.mean(master_bias.data)
bias_std = np.std(master_bias.data)

pl.figure(figsize=(15, 15))
pl.imshow(master_bias, vmax=bias_mean + 4*bias_std, vmin=bias_mean - 4*bias_std, cmap='Greys_r')

Make a master flat



In [ ]:

    
flats = []
for flat, fname in images.hdus(imagetyp='flat', return_fname=True):
    meta = flat.header
    meta['filename'] = fname
    flats.append(ccdproc.CCDData(flat.data, meta=meta, unit="adu"))

Now subtract overscan and trim



In [ ]:

    
oscan_and_trim(flats)

Correct for bias



In [ ]:

    
for i,flat in enumerate(flats):
    flats[i] = ccdproc.subtract_bias(flat, master_bias)

Combine flats with sigma clipping



In [ ]:

    
flat_combiner = ccdproc.Combiner(flats)



In [ ]:

    
# flat_combiner.sigma_clipping(func=med_over_images)
flat_combiner.sigma_clipping(func=lambda x,axis: np.nanmedian(x,axis))



In [ ]:

    
master_flat = flat_combiner.median_combine()



In [ ]:

    
flat_mean = np.mean(master_flat.data)
flat_std = np.std(master_flat.data)

pl.figure(figsize=(15, 15))
pl.imshow(master_flat, vmax=flat_mean + 5*flat_std, vmin=flat_mean - 4*flat_std, cmap='Greys_r')



In [ ]:

    
master_flat_e = ccdproc.gain_correct(master_flat, gain=gain)



In [ ]:

    
obj_img = ccdproc.CCDData.read("/Users/adrian/projects/triand-rrlyrae/data/mdm-fall-2015/n1/m110315.0036.fit", 
                               unit=u.adu)
obj_img = [obj_img]

oscan_and_trim(obj_img)
obj_img = obj_img[0]



In [ ]:

    
obj_img = ccdproc.subtract_bias(obj_img, master_bias)



In [ ]:

    
obj_img_corr = ccdproc.flat_correct(obj_img, master_flat)



In [ ]:

    
d = obj_img_corr

_mean = np.mean(d.data)
_std = np.std(d.data)

pl.figure(figsize=(15, 15))
pl.imshow(d, vmax=_mean + 5*_std, vmin=_mean - 5*_std, cmap='Greys_r')



In [ ]:

    
from scipy.signal import medfilt2d

d = medfilt2d(obj_img_corr.data, kernel_size=(3,3))

_mean = np.mean(d.data)
_std = np.std(d.data)

pl.figure(figsize=(15, 15))
pl.imshow(d, vmax=_mean + 5*_std, vmin=_mean - 5*_std, cmap='Greys_r')



In [ ]:

    
sub_d = d[:,50:100]



In [ ]:

    
x = np.arange(sub_d.shape[1])
pl.plot(sub_d[1024], marker=None, drawstyle='steps')
pl.axvline(x[sub_d[1024].argmax()])
pl.axvline(x[sub_d[1024].argmax()-7], alpha=0.4)
pl.axvline(x[sub_d[1024].argmax()+7], alpha=0.4)



In [ ]:

    
amps = np.zeros(sub_d.shape[0])
for i,row in enumerate(sub_d):
    sky = np.median(row)
    idx = row.argmax()
    amps[i] = (row[idx-7:idx+7] - sky).sum()



In [ ]:

    
pl.figure(figsize=(12,8))
pl.plot(amps, marker=None, drawstyle='steps')
pl.xlim(1700, 1000)



In [ ]: