#!/usr/bin/env python
# coding: utf-8

# # Combining flats
# 
# There is one step in combining flats that is different from most other image
# combination: the flats should be scaled to a common value before combining them.
# This is particularly important if the flats are twilight flats in which the
# average image value typically changes significantly as the images are being
# taken.
# 
# Flats are typically grouped by filter when combining them. That is, one combined
# flat is produced for each filter in which flats were taken.
# 
# Combination will be done for each of the two examples in the previous notebook.

# In[ ]:


from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt

import ccdproc as ccdp
from astropy.stats import mad_std
from astropy.visualization import hist

from convenience_functions import show_image


# In[ ]:


# Use custom style for larger fonts and figures
plt.style.use('guide.mplstyle')


# ## Example 1
# 
# We begin by setting up an image collection for the reduced data. This data is
# from chip 0 of the cryogenically-cooled Large Format Camera at Palomar
# Observatory.

# In[ ]:


calibrated_path = Path('example1-reduced')

flat_imagetyp = 'flatfield'

ifc = ccdp.ImageFileCollection(calibrated_path)


# We'll first check what filters are present.

# In[ ]:


flat_filters = set(h['filter'] for h in ifc.headers(imagetyp=flat_imagetyp))
flat_filters


# These flats are dome flats, essentially pictures of a screen in the dome
# illuminated by a light source, so you would not expect there to be much variable
# in the typical pixel value between different exposures. There is typically
# *some* variation, though, so we graph it below.

# In[ ]:


median_count = [np.median(data) for data in ifc.data(imagetyp=flat_imagetyp)]
mean_count = [np.mean(data) for data in ifc.data(imagetyp=flat_imagetyp)]
plt.plot(median_count, label='median')
plt.plot(mean_count, label='mean')
plt.xlabel('Image number')
plt.ylabel('Count (ADU)')
plt.title('Pixel value in calibrated flat frames')
plt.legend()
print(median_count)


# Although this is less frame-to-frame variation than we will see in Example 2, it
# is about 5%. If we were to combine these without scaling the flats to a common
# value then the images with higher counts would effectively get more weight than
# the images.
# 
# There is a substantial difference between the mean and median of this data.
# Typically it is better to use the median because extreme values do not affect
# the median as much as the mean.
# 
# To scale the frames so that they have the same median value, we need to define a
# function that can calculate the inverse of the median given the data.

# In[ ]:


def inv_median(a):
    return 1 / np.median(a)


# This function is passed into the `scale` argument of `combine` below. One
# combined flat is created for each filter in the data.

# In[ ]:


for filt in flat_filters:
    to_combine = ifc.files_filtered(imagetyp=flat_imagetyp, filter=filt, include_path=True)
    combined_flat = ccdp.combine(to_combine,
                                 method='average', scale=inv_median,
                                 sigma_clip=True, sigma_clip_low_thresh=5, sigma_clip_high_thresh=5,
                                 sigma_clip_func=np.ma.median, signma_clip_dev_func=mad_std,
                                 mem_limit=350e6
                                )

    combined_flat.meta['combined'] = True
    dark_file_name = 'combined_flat_filter_{}.fit'.format(filt.replace("''", "p"))
    combined_flat.write(calibrated_path / dark_file_name)


# ### Discussion of Example 1

# We will begin by checking that the right number of combined flats have been
# created. There were two filters, `g'` and `i'` in the raw data so there should
# be two combined flats. We need to refresh the `ImageFileCollection` for the
# reduced data because new files, our flats, have been added to them.

# In[ ]:


ifc.refresh()
ifc.files_filtered(imagetyp=flat_imagetyp, combined=True)


# That looks good. The two flats are displayed below.

# In[ ]:


fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 20), tight_layout=True)

for ccd, axis in zip(ifc.ccds(imagetyp=flat_imagetyp, combined=True), axes):
    show_image(ccd.data, cmap='gray', fig=fig, ax=axis)
    title = "Filter {}".format(ccd.header['filter'])
    axis.set_title(title)


# The first thing to notice is that the flats are different in these two filters.
# That is expected because one of the elements in the optical path, the filter, is
# different.
# 
# The pattern of electronics in the flat images is because this is a
# backside-illuminated CCD. The light-detecting pixels are on the under side of
# the chip and the light needs to pass through the chip to reach the sensor. The
# small dark spots are places where the chip wasn't thinned uniformly.
# 
# Compare this with [Example 2](#discussion-of-example-2) below, which shows a flat
# taken with a frontside-illuminated camera.

# ## Example 2
# 
# The data in this example is from a thermoelectrically-cooled Andor Aspen CG16M.
# These flats are twilight flats, taken just after sunset.

# In[ ]:


calibrated_path = Path('example2-reduced')

flat_imagetyp = 'flat'

ifc = ccdp.ImageFileCollection(calibrated_path)


# In[ ]:


flat_filters = set(h['filter'] for h in ifc.headers(imagetyp=flat_imagetyp))
flat_filters


# Twilight flats typically differ more frame-to-frame than dome flats, as
# shown in the figure below.

# In[ ]:


median_count = [np.median(hdu.data) for hdu in ifc.hdus(imagetyp=flat_imagetyp)]
mean_count = [np.mean(data) for data in ifc.data(imagetyp=flat_imagetyp)]
plt.plot(median_count, label='median')
plt.plot(mean_count, label='mean')
plt.xlabel('Image number')
plt.ylabel('Count (ADU)')
plt.title('Pixel value in calibrated flat frames')
plt.legend()
print(median_count)


# In[ ]:


def inv_median(a):
    return 1 / np.median(a)

for filt in flat_filters:
    to_combine = ifc.files_filtered(imagetyp=flat_imagetyp, filter=filt, include_path=True)
    combined_flat = ccdp.combine(to_combine,
                                 method='average', scale=inv_median,
                                 sigma_clip=True, sigma_clip_low_thresh=5, sigma_clip_high_thresh=5,
                                 sigma_clip_func=np.ma.median, signma_clip_dev_func=mad_std,
                                 mem_limit=350e6
                                )

    combined_flat.meta['combined'] = True
    dark_file_name = 'combined_flat_filter_{}.fit'.format(filt.replace("''", "p"))
    combined_flat.write(calibrated_path / dark_file_name)


# ### Discussion of Example 2

# We expect only one combined flat because there was only one filter. The
# `ImageFileCollection` is refreshed before we query it because the combined flats
# were added after the collection was created.

# In[ ]:


ifc.refresh()
ifc.files_filtered(imagetyp=flat_imagetyp, combined=True)


# The combined flat is displayed below.

# In[ ]:


show_image(combined_flat, cmap='gray', figsize=(10, 20))


# This flat looks very different than the one in [Example 1](#discussion-of-example-1)
# because this CCD is frontside-illuminated and the previous one is backside-illuminated. 
# That means the sensor is on the top of the chip and the light does
# not pass through the sensor chip itself to reach the sensors. Though only one
# filter is shown here, the flat field looks slightly different through other
# filters on this camera.