{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start, as always, by loading the libraries we will use. These are mostly libraries within the astropy package, but also note the loading of the ccdproc package, which contain many of the core functions for processing CCD data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from astropy.io import fits\n",
    "from astropy.visualization import hist\n",
    "from astropy import units as u\n",
    "from astropy.stats import mad_std\n",
    "from pathlib import Path\n",
    "from astropy.nddata import CCDData\n",
    "import ccdproc as ccdp\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we need to set the directories. The first is the \"raw directory,\" where the raw images reside. If you're doing this on speedy, that is already pointing to the right place. If you want to do the reduction on your own computer, you will need to adjust that to point to wherever you put the raw data files. The second is the \"reduceddir,\" where the results of your work will be stored. You can leave that one as-is and it will just store it in a sub-directory of your current location. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rawdir=Path('/home/hugh/ngc55-raw/')\n",
    "reduceddir=Path('ngc55-reduced/')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This will make the reduced directory if it doesn't already exist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reduceddir.mkdir(exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we make an \"ImageFileCollection\" which, as the name suggests, is a data type that includes a list/collection of all images in a directory. In this case, let's make a list of all the raw images. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ifc_raw = ccdp.ImageFileCollection(rawdir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's list a summary of the properties of the images in that directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ifc_raw.summary['file', 'imagetyp', 'filter', 'exptime', 'naxis1', 'naxis2']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this list ist truncated (it has the \"...\" in the middle of the list). This likely means that we're missing some image types. Let's figure out what image types exist in this ImageFileCollection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_imagetypes=set(h['imagetyp'] for h in ifc_raw.headers())\n",
    "all_imagetypes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the results here, we can see that there are different image types for the different kinds of image. Using what you found here, edit this next block of code to define the names for the bias frames, flat frames, and science frames. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "science_imagetyp=''\n",
    "bias_imagetype=''\n",
    "flat_imagetyp=''"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will make a list of the bias frames (using the files_filtered method), combine those biases into a single averaged bias frame, update the header to indicate that it was a combined image, and write that image out to the reduced data directory. \n",
    "\n",
    "Before your run this code block, take a look at the options for ccdp.combine. Using the documentation that you read through, explain what each of the options for ccdp.combine do (i.e. \"method\", \"sigma_clip\", etc.). Include this description in your lab report. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_biases = ifc_raw.files_filtered(include_path=True, imagetyp=\"Bias Frame\")\n",
    "\n",
    "combined_bias = ccdp.combine(raw_biases,\n",
    "                             method='average',\n",
    "                             sigma_clip=True, \n",
    "                             sigma_clip_low_thresh=5, \n",
    "                             sigma_clip_high_thresh=5,\n",
    "                             sigma_clip_func=np.ma.median, \n",
    "                             sigma_clip_dev_func=mad_std,\n",
    "                             unit='adu',\n",
    "                             mem_limit=16e9\n",
    "                            )\n",
    "\n",
    "combined_bias.meta['combined'] = True\n",
    "combined_bias.write(reduceddir / 'combined_bias.fit', overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An important note: throughout this data reduction process, you may get warnings (often repeated warnings) about the sofware making adjustments to headers. For these data, the hardware cameras at iTelescope use some non-standard entries for the date. This is fine, but astropy freaks out about this. Don't be alarmed about warnings of the form:\n",
    "\n",
    ">INFO:astropy:using the unit adu passed to the FITS reader instead of the unit adu in the FITS file.\n",
    "\n",
    ">WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 60468.566390 from DATE-OBS'. [astropy.wcs.wcs]\n",
    "\n",
    ">WARNING:astropy:FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 60468.566390 from DATE-OBS'."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let's read in that combined bias image we just made and also make a list of images that are registered as flat fields. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combinedbias= CCDData.read(reduceddir / 'combined_bias.fit', unit=\"adu\")\n",
    "lightfiles=ifc_raw.files_filtered(imagetyp=flat_imagetyp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we will subtract the combined bias from each of the images that are labeled as a \"Flat\" and write that image to the reduced data directory"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ccd, file_name in ifc_raw.ccds(imagetyp='FLAT',\n",
    "                               return_fname=True\n",
    "                                ):\n",
    "    ccd=ccdp.subtract_bias(ccd,combinedbias)\n",
    "    ccd.write(reduceddir / file_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's set up an ImageFileCollection for the reduced data directory, since we now have files in there we want to work with (i.e. the bias-subtracted Flats)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ifc_reduced=ccdp.ImageFileCollection(reduceddir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's make a list of all the different filter names that exist for the flat field images and print that list out. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flat_filters=set(h['filter'] for h in ifc_reduced.headers(imagetyp=flat_imagetyp))\n",
    "flat_filters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're going to scale each flat field by the inverse median, so let's define a function that we can use when we need it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def inv_median(a):\n",
    "    return 1 / np.median(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK. Now, for each filter, we will combine the flat field images to make a combined flat. Note that we are careful to ONLY combine blue flat images with other blue flat images, green flat images with other green flat images, etc. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for filt in flat_filters:\n",
    "    to_combine = ifc_reduced.files_filtered(imagetyp=flat_imagetyp, filter=filt, include_path=True)\n",
    "    combined_flat = ccdp.combine(to_combine,\n",
    "                                 method='average', scale=inv_median,\n",
    "                                 sigma_clip=True, sigma_clip_low_thresh=5, sigma_clip_high_thresh=5,\n",
    "                                 sigma_clip_func=np.ma.median, signma_clip_dev_func=mad_std,\n",
    "                                 mem_limit=350e6\n",
    "                                )\n",
    "\n",
    "    combined_flat.meta['combined'] = True\n",
    "    flat_file_name = 'combined_flat_filter_{}.fit'.format(filt.replace(\"''\", \"p\"))\n",
    "    combined_flat.write(reduceddir / flat_file_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, you're going to have a long list of error messages above this all about 'datfix'. Don't worry about those!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK. Now, we're going to flat- and bias-correct all of our science images. First, let's see all of our science images:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ifc_reduced=ccdp.ImageFileCollection(reduceddir)\n",
    "sciences = ifc_raw.summary[ifc_raw.summary['imagetyp']==science_imagetyp]\n",
    "sciences['date-obs', 'file', 'imagetyp', 'object', 'filter', 'exptime']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we will set up a dictionary for the flats, since we want to make sure that each science image is matched up with the correct filter (i.e. we don't want to correct a science image taken through a Green filter with a Red flat). We'll also just find the combined bias and read that into a variable to use in the loop below when we bias subtract the science images. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "combined_flats = {ccd.header['filter']: ccd for ccd in ifc_reduced.ccds(imagetyp=flat_imagetyp, combined=True)}\n",
    "\n",
    "# There is only on bias, so no need to set up a dictionary.\n",
    "combined_bias = ifc_reduced.ccds(imagetyp='bias', combined=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK. Here's the main show. We're going to bias correct and flat correct each of the science images. Why is it important to bias correct your images before you flat field them? Include the answer to this question in your lab writeup."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_reds = []\n",
    "light_ccds = []\n",
    "for light, file_name in ifc_raw.ccds(imagetyp=science_imagetyp, return_fname=True):\n",
    "    light_ccds.append(light)\n",
    " \n",
    "    # Note that the first argument in the remainder of the ccdproc calls is\n",
    "    # the *reduced* image, so that the calibration steps are cumulative.\n",
    "    reduced = ccdp.subtract_bias(light, combinedbias)\n",
    "    \n",
    "    good_flat = combined_flats[reduced.header['filter']]\n",
    "    reduced = ccdp.flat_correct(reduced, good_flat)\n",
    "    all_reds.append(reduced)\n",
    "    reduced.write(reduceddir / file_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will re-create the ImageFileCollection, since now there are new files (i.e. the reduced flats) in the directory that we want to use in our final step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ifc_reduced=ccdp.ImageFileCollection(reduceddir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, for each filter, we will average together the science images. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "science_filters=set(h['filter'] for h in ifc_reduced.headers(imagetyp=science_imagetyp))\n",
    "for filt in science_filters:\n",
    "    to_combine = ifc_reduced.files_filtered(imagetyp=science_imagetyp, filter=filt, include_path=True)\n",
    "    combined_science = ccdp.combine(to_combine,\n",
    "                                 method='average',\n",
    "                                 sigma_clip=True, sigma_clip_low_thresh=5, sigma_clip_high_thresh=5,\n",
    "                                 sigma_clip_func=np.ma.median, signma_clip_dev_func=mad_std,\n",
    "                                 mem_limit=350e6\n",
    "                                )\n",
    "\n",
    "    combined_science.meta['combined'] = True\n",
    "    science_file_name = 'combined_science_filter_{}.fit'.format(filt.replace(\"''\", \"p\"))\n",
    "    combined_science.write(reduceddir / science_file_name)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}