import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from astropy.io import fits

data = ('../data/Observing Project/raw_data/ngc7789_b_30_1.fit')
image_data = fits.getdata(data)
image_hdr  = fits.getheader(data)
image_data

fig, axis = plt.subplots(figsize=(10,10))
implot = axis.imshow(image_data,cmap='gray',origin='lower',vmin=1500,vmax=2500)
axis.grid(False)
plt.show()

print( image_hdr['EXPOSURE'] )

hdu  = fits.PrimaryHDU(data=image_data, header=image_hdr)
hdu.writeto('my_new_fits_file.fits',clobber=True)

import ccdproc # the ccdproc library
from reduction_tools import Combiner, Reduction # Stu's custom tools for combining and `reducing` CCD data.

data_dir = '../data/Observing Project/raw_data'
destination_dir = '../data/Observing Project/calib_data'

import os
#make destination directory if it does not exist
if not os.path.exists(destination_dir):
    os.mkdir(destination_dir)

images = ccdproc.ImageFileCollection(location=data_dir, keywords=['imagetyp', 'exposure', 'filter'])
print(images.summary_info)

#changes headers
from astropy.io import fits
hdu_list = fits.open('../data/Observing Project/raw_data/ngc7789_v_30_3.fit',mode='update')
data = hdu_list[0].data
header = hdu_list[0].header

print(header['IMAGETYP'])
header['IMAGETYP'] = "LIGHT"

hdu_list.flush()

# define the Reduction task
copy_file_task = Reduction(input_image_collection=images, 
                           destination=destination_dir, 
                           apply_to={'imagetyp': 'BIAS'},
                           copy_only=True)

# the line above simply creates the Reduction task object. To run it, call it's `action` function, like so
copy_file_task.action()

reduced_images = ccdproc.ImageFileCollection(location=destination_dir, keywords=['imagetyp', 'exposure', 'filter', 'master'])

make_master_bias = Combiner(image_source=reduced_images,
                            file_name_base='master_bias',
                            combine_method='average',
                            apply_to={'imagetyp':'BIAS'},
                            destination=destination_dir)
make_master_bias.action()

os.listdir(destination_dir)

# use os library to glue directories and filenames together easily
filepath = os.path.join(destination_dir, 'master_bias.fit')
print(filepath)

# open the file and store data in array called 'masterbias'
masterbias = fits.getdata(filepath)

# What's the mean level in the bias image?
print('The masterbias has mean level of {:.1f} counts'.format( masterbias.mean() ))

# display it
fig, axis = plt.subplots(figsize=(5,5))
# custom scaling - see companion notebook
vmin,vmax = np.percentile(masterbias, [1,99])
implot = axis.imshow(masterbias,cmap='gray',origin='lower',vmin=vmin,vmax=vmax)
axis.grid(False)
plt.show()

#It is best to combine the dark frames using the median, rather than the mean. Dark frames are often long exposures, which can be affected by cosmic rays. 
#Cosmic rays hitting the CCD will excite also excite electrons. Taking the median of the master dark will help remove cosmic rays from the master dark frame.

dark_preprocessing = Reduction(bias_subtract=True,
                               dark_subtract=False,
                               flat_correct=False,
                               master_source=reduced_images,
                               input_image_collection=images,
                               destination=destination_dir,
                               apply_to={'imagetyp': 'dark'})
dark_preprocessing.action()

#make a master dark

make_master_dark = Combiner(file_name_base='master_dark',
                            combine_method='median',
                            image_source=reduced_images,
                            apply_to={'imagetyp':'DARK'},
                            group_by='exposure',
                            destination=destination_dir)
make_master_dark.action()

os.listdir(destination_dir)

filepath = os.path.join(destination_dir,"master_dark_exposure_300.0.fit")

masterdark_im = fits.getdata(filepath)

# What's the mean level in the dark image?
print('The masterdark has mean level of {:.1f} counts'.format( masterdark_im.mean() ))

# display it
fig, axis = plt.subplots(figsize=(5,5))
# custom scaling - see companion notebook
vmin,vmax = np.percentile(masterdark_im, [1,99])
implot = axis.imshow(masterdark_im,cmap='gray',origin='lower',vmin=vmin,vmax=vmax)
axis.grid(False)
plt.show()

dark_fits_filename = os.path.join(destination_dir, 'dark_300_1.fit') # use os.path.join to glue together directory and filename
img_hdr = fits.getheader(dark_fits_filename) # read in FITS header using astropy's fits library
print(img_hdr['EXPOSURE'])

#process flats
flat_reduction = Reduction(bias_subtract=True,
                           dark_subtract=False,
                           flat_correct=False,
                           master_source=reduced_images,
                           input_image_collection=images,
                           destination=destination_dir,
                           apply_to={'imagetyp':'FLAT'})
flat_reduction.action()

#create master flat
make_master_flat = Combiner(file_name_base='master_flat',
                            group_by='filter',
                            image_source=reduced_images,
                            destination=destination_dir,
                            apply_to={'imagetyp':'FLAT'})
make_master_flat.action()

#plot master flat
filepath = os.path.join(destination_dir,"master_flat_filter_V.fit")

masterflat_im = fits.getdata(filepath)

# What's the mean level in the dark image?
print('The masterflat has mean level of {:.1f} counts'.format( masterflat_im.mean() ))

# display it
fig, axis = plt.subplots(figsize=(6,6))
# custom scaling - see companion notebook
vmin,vmax = np.percentile(masterflat_im, [1,99])
implot = axis.imshow(masterflat_im,cmap='gray',origin='lower',vmin=vmin,vmax=vmax)
axis.grid(False)
plt.show()

#subtract master bias
#subtract master dark
#divide by master flat
#save calib image
data_reduction = Reduction(bias_subtract=True,
                           dark_subtract=True,
                           flat_correct=True,
                           master_source=reduced_images,
                           input_image_collection=images,
                           destination=destination_dir,
                           apply_to={'imagetyp':'LIGHT'})
data_reduction.action()

image_filepath=os.path.join(destination_dir, "MED_b_30_aligned.fits")
m51_image=fits.getdata(image_filepath)

fig, axis = plt.subplots(figsize=(15,6))
implot = axis.imshow(m51_image,cmap='gray',origin='lower',vmin=200,vmax=1500)
axis.grid(False)
plt.show()