{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calibrating bias images"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The purpose of calibrating bias images is threefold:\n",
    "\n",
    "+ Subtract overscan if you have decided your science will be better if you\n",
    "subtract overscan. See [this discussion of overscan](01-08-Overscan.ipynb) for some guidance.\n",
    "+ Trim the overscan region off of the image if it is present, regardless of\n",
    "whether or not you have chosen to subtract the overscan.\n",
    "+ Combine the bias images into a \"combined\" bias to be used in calibrating the\n",
    "rest of the images. The purpose of combining several images is to reduce as much\n",
    "as possible the read noise in the combined bias.\n",
    "\n",
    "The approach in this notebook will be to reduce a single image, look at the\n",
    "effects the reduction step had on that image and then demonstrate how to\n",
    "calibrate a folder containing several images of that type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pathlib import Path\n",
    "import os\n",
    "\n",
    "from astropy.nddata import CCDData\n",
    "from astropy.visualization import hist\n",
    "import ccdproc as ccdp\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from convenience_functions import show_image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use custom style for larger fonts and figures\n",
    "plt.style.use('guide.mplstyle')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data for these examples"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See the [Preface notebook](00-00-Preface.ipynb) for download links for all data."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 1: With overscan subtraction"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Decide where to put your Example 1 calibrated images\n",
    "Though it is possible to overwrite your raw data with calibrated images, this is not recommended. Here we create a folder called `example1-reduced` that will contain\n",
    "the calibrated data and create it if it doesn't exist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calibrated_data = Path('.', 'example1-reduced')\n",
    "calibrated_data.mkdir(exist_ok=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make an image file collection for the raw data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "example_cryo_path = Path('example-cryo-LFC')\n",
    "files = ccdp.ImageFileCollection(example_cryo_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files.summary['file', 'imagetyp', 'filter', 'exptime', 'naxis1', 'naxis2']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "darks_only = ccdp.ImageFileCollection(example_cryo_path / 'darks')\n",
    "darks_only.summary['file', 'imagetyp', 'exptime']"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determine overscan region for the LFC Chip 0"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please see the discussion of this camera in\n",
    "[the Overscan notebook](01-08-Overscan.ipynb#case-1-cryogenically-cooled-large-format-camera-lfc-at-palomar) for the appropriate overscan region\n",
    "to use for this camera. Note, in particular, that it differs from the value\n",
    "given in the `BIASSEC` keyword in the header of the images."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The astropy affiliated package [ccdproc](https://ccdproc.readthedocs.io) provides two\n",
    "useful functions here:\n",
    "\n",
    "+ `subtract_overscan` for subtracting the overscan from the image, and\n",
    "+ `trim_image` for trimming off the overscan.\n",
    "\n",
    "First, let's see the values of `BIASSEC`, which sometimes (but do not always)\n",
    "indicate that there is is overscan and which part of the chip is the overscan,\n",
    "as well as the values of `CCDSEC`, which is sometimes but not always present, and indicates which\n",
    "part of the chip light hit.\n",
    "\n",
    "Note that neither of these are standard; sometimes, for example, `trimsec` is\n",
    "used instead of `ccdsec`, and there are likely other variants. Some images may\n",
    "have neither keyword in the header. This does not necessarily indicate that\n",
    "ovserscan isn't present. The best advice is to carefully check the documentation\n",
    "for the camera you are using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files.summary['file', 'imagetyp', 'biassec', 'ccdsec', 'datasec'][0]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fits header claims the overscan extends from the 2049$^{th}$ column to the\n",
    "end of the image (this is one-based indexing) and that the part of the image\n",
    "exposed to light extends over all rows and from the first column to the\n",
    "2048$^{th}$ column (again, this is one-indexed)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### FITS vs. Python indexing\n",
    "\n",
    "There are two differences between FITS and Python in terms of indexing:\n",
    "\n",
    "+ Python indexes are zero-based (i.e., numbering starts at zero), while FITS indexes\n",
    "are one-based (i.e., numbering starts at one).\n",
    "+ The *order* of the indexes is swapped.\n",
    "\n",
    "For example, the **FITS** representation of the part of the chip exposed to\n",
    "light is `[1:2048,1:4128]`. To access that part of the data from a NumPy array\n",
    "in **Python**, switch the order so that the indexing looks like this: `[0:4128,\n",
    "0:2048]` (or, more compactly `[:, :2048]`). Note that the *ending* indexes given\n",
    "here for Python are correct because the second part of a range (after the colon)\n",
    "is *not included* in the array slice. For example, `0:2048` starts at 0 (the\n",
    "first pixel) and goes up to but does not include 2048, so the last pixel included\n",
    "is `2047` (the 2048$^{th}$ pixel)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As discussed in [the Overscan notebook](01-08-Overscan.ipynb#case-1-cryogenically-cooled-large-format-camera-lfc-at-palomar), the useful\n",
    "overscan region for this camera starts at the 2055$^{th}$ column, not column\n",
    "2049 as indicated by the `BIASSEC` keyword in the header. This situation is not\n",
    "unusual; column 2049 is the first of the columns masked from light by the manufacturer, \n",
    "but there is some leakage into this region from the rest of the CCD.\n",
    "\n",
    "If you are going to use overscan you need to carefully examine the overscan in a few\n",
    "representative images to understand which part of the overscan to use.\n",
    "\n",
    "In what follows, we will use for the overscan the region (Python/NumPy indexing)\n",
    "`[:, 2055:]`."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subtract and then trim the overscan (one sample image)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using `subtract_overscan` is reasonably concise, as shown in the cell\n",
    "below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_biases = files.files_filtered(include_path=True, imagetyp='BIAS')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "first_bias = CCDData.read(raw_biases[0], unit='adu')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bias_overscan_subtracted = ccdp.subtract_overscan(first_bias, overscan=first_bias[:, 2055:], median=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we trim off the full overscan region (not just the part we used for\n",
    "subtracting overscan)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trimmed_bias = ccdp.trim_image(bias_overscan_subtracted[:, :2048])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(20, 10))\n",
    "\n",
    "show_image(first_bias.data, cmap='gray', ax=ax1, fig=fig)\n",
    "ax1.set_title('Raw bias')\n",
    "show_image(trimmed_bias.data, cmap='gray', ax=ax2, fig=fig)\n",
    "ax2.set_title('Bias, overscan subtracted and trimmed')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Discussion\n",
    "\n",
    "Visually, the images look nearly identical before and after calibration. The\n",
    "only prominent difference is a shift in the pixel values, as one would expect\n",
    "from subtracting the same value from each pixel in an image. It simply shifts\n",
    "the zero point.\n",
    "\n",
    "There is one other important difference between the images: the input image\n",
    "uses 32MB of memory while the calibrated, overscan-subtracted image uses roughly\n",
    "128MB. The input image is stored as unsigned 16-bit integers; the calibrated\n",
    "image is stored as floating point numbers, which default in Python to 64-bit\n",
    "floats. The memory size is also the size the files will have when written to\n",
    "disk (ignoring any compression). You can reduce the memory and disk footprint by\n",
    "changing the `dtype` of the image: `trimmed_bias.data = trimmed_bias.data.astype('float32')`. It is best\n",
    "to do this just before writing the image out because arithmetic operations on\n",
    "the image may convert its `dtype` back to `float64`.\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Processing the folder of bias images for the LFC Chip 0"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Processing each of the bias images individually would be tedious, at best.\n",
    "Instead, we can use the [`ImageFileCollection`](https://ccdproc.readthedocs.io/en/latest/ccdproc/image_management.html) we created above to\n",
    "loop over only the bias images, saving each in the folder `calibrated_data`. In\n",
    "this example, the files are saved uncompressed because the Python library for\n",
    "compressing gzip files is extremely slow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ccd, file_name in files.ccds(imagetyp='BIAS',            # Just get the bias frames\n",
    "                                 ccd_kwargs={'unit': 'adu'}, # CCDData requires a unit for the image if \n",
    "                                                             # it is not in the header\n",
    "                                 return_fname=True           # Provide the file name too.\n",
    "                                ):\n",
    "        # Subtract the overscan\n",
    "    ccd = ccdp.subtract_overscan(ccd, overscan=ccd[:, 2055:], median=True)\n",
    "    \n",
    "    # Trim the overscan\n",
    "    ccd = ccdp.trim_image(ccd[:, :2048])\n",
    "    \n",
    "    # Save the result\n",
    "    ccd.write(calibrated_data / file_name)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's check that we really did get the images we expect by creating an\n",
    "[`ImageFileCollection`](https://ccdproc.readthedocs.io/en/latest/ccdproc/image_management.html) for the reduced folder and displaying the size\n",
    "of each image. We are expecting the images to be 2048 × 4128, and that there\n",
    "will be the same number of reduced bias images as input bias images (six)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reduced_images = ccdp.ImageFileCollection(calibrated_data)\n",
    "reduced_images.summary['file', 'imagetyp', 'naxis1', 'naxis2']"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2: No overscan subtraction, but trim the images"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are not subtracting overscan then the only manipulation you may need to\n",
    "do is trimming the overscan from the images. If there is no overscan region in\n",
    "your images then even that is unnecessary."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Decide where to put your calibrated Example 2 images"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Though it is possible to overwrite your raw data with calibrated images, this is not recommended. Here we create a folder called `example2-reduced` that will contain\n",
    "the calibrated data and create it if it doesn't exist."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calibrated_data = Path('.', 'example2-reduced')\n",
    "calibrated_data.mkdir(exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = ccdp.ImageFileCollection('example-thermo-electric')\n",
    "files.summary['file', 'imagetyp', 'filter', 'exptime', 'naxis1', 'naxis2']"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Determine overscan region for this camera"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please see the discussion of this camera in [the Overscan notebook](01-08-Overscan.ipynb#case-2-thermo-electrically-cooled-apogee-aspen-cg16m) for\n",
    "a discussion of the overscan region of this camera. The overscan for this camera\n",
    "is not useful but should be trimmed out at this stage.\n",
    "\n",
    "These headers have some information in the keywords `BIASSEC` and `TRIMSEC`\n",
    "indicating, in the FITS numbering convention, the overscan region and the\n",
    "science region of the chip."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files.summary['file', 'imagetyp', 'biassec', 'trimsec'][0]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on this, and the decision not to subtract overscan for this camera, we\n",
    "will only need to trim the overscan region off of the images. See the discussion\n",
    "at [FITS vs. Python indexing](#fits-vs-python-indexing), above, for some details about the\n",
    "difference between FITS and Python indexing. Essentially, to get Python indexes\n",
    "from FITS, reverse the order and subtract one."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Trim the overscan (one sample image)\n",
    "\n",
    "The function `trim_image` from [ccdproc](https://ccdproc.readthedocs.io) removes a\n",
    "portion of the image and updates the image metadata as needed.\n",
    "\n",
    "Below we get the first bias image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "raw_biases = files.files_filtered(include_path=True, imagetyp='BIAS')\n",
    "\n",
    "first_bias = CCDData.read(raw_biases[0], unit='adu')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two ways of specifying the region to trim. One is to slice the image in\n",
    "Python; the other is to use the `fits_section` argument to `trim_image`.\n",
    "\n",
    "The cell below uses a FITS-style section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trimmed_bias_fits = ccdp.trim_image(first_bias, fits_section='[1:4096, :]')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cell below does the same trimming as the one above, but with Python-style\n",
    "slicing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trimmed_bias_python = ccdp.trim_image(first_bias[:, :4096])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.testing.assert_allclose(trimmed_bias_python, trimmed_bias_fits)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Processing the folder of bias images for Example 2"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in [Example 1](#example-1-with-overscan-subtraction), above, we can use the\n",
    "[`ImageFileCollection`](https://ccdproc.readthedocs.io/en/latest/ccdproc/image_management.html) we created  to loop over only the bias images,\n",
    "saving each in the folder `calibrated_data`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ccd, file_name in files.ccds(imagetyp='BIAS',            # Just get the bias frames\n",
    "                                 return_fname=True           # Provide the file name too.\n",
    "                                ):    \n",
    "    # Trim the overscan\n",
    "    ccd = ccdp.trim_image(ccd[:, :4096])\n",
    "    \n",
    "    # Save the result\n",
    "    ccd.write(calibrated_data / file_name)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3: No overscan at all\n",
    "\n",
    "If there is no overscan then there is, in principle, nothing to be done with the\n",
    "bias frames. It may be convenient to copy them to the directory with the rest of\n",
    "your reduced images. The code below does that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calibrated_data = Path('.', 'example3-reduced')\n",
    "calibrated_data.mkdir(exist_ok=True)\n",
    "\n",
    "biases = files.files_filtered(imagetyp='BIAS', include_path=True)\n",
    "\n",
    "import shutil\n",
    "\n",
    "for bias in biases:\n",
    "    shutil.copy(bias, calibrated_data)"
   ]
  }
 ],
 "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}