{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Image stacking, summing, median.\n",
    "\n",
    "Images of M13 taken with a 10 inch telescope.\n",
    "\n",
    "This dataset is downloaded from the astropy.stsci.edu server.\n",
    "\n",
    "Import numpy, matplotlib, pyplot, download_file, fits\n",
    "\n",
    "Download 5 files using 'download_file',  print the names of the downloaded files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.utils.data import download_file\n",
    "from astropy.io import fits\n",
    "image_list = [ download_file('http://data.astropy.org/tutorials/FITS-images/M13_blue_000'+n+'.fits', cache=True ) \\\n",
    "              for n in ['1','2','3','4','5'] ]\n",
    "for im_name in image_list:\n",
    "    print(im_name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a list containing the 5 data sets.\n",
    "\n",
    "Also create a 'final_image' containing zeros that is the same shape as our downloaded files.\n",
    "\n",
    "Sum the 5 images into the 'final_image' and divide by 5 to get an average.\n",
    "\n",
    "Print the mean value of the final image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_concat = []\n",
    "for image in image_list:\n",
    "    image_concat.append(fits.getdata(image))\n",
    "    \n",
    "    final_image = np.zeros(shape=image_concat[0].shape)\n",
    "\n",
    "for image in image_concat:\n",
    "    final_image += image\n",
    "    \n",
    "final_image = final_image/5.0\n",
    "print(np.mean(final_image))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot each of the 5 input images to we what we are working with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for image in image_list:\n",
    "    plt.imshow(fits.getdata(image), cmap='gray', origin=\"lower\")\n",
    "    plt.figure()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot a histogram of this average image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "image_hist = plt.hist(final_image.flatten(), 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the final, average image using a larger size figure and with min and max set from the histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15,15))\n",
    "plt.imshow(final_image, cmap='gray', origin='lower',vmin=450, vmax=600)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we use numpy to construct a 3D data cube out of the 5 images.\n",
    "\n",
    "Then we make an image consisting of the median value of each pixel by running median on axis=2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stack = np.dstack(image_concat)\n",
    "print(stack.shape)\n",
    "med_data = np.median(stack, axis=2)\n",
    "med_hist = plt.hist(med_data.flatten(), 100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot this median image the same way we did the average image for comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15,15))\n",
    "plt.imshow(med_data, cmap='gray', origin = 'lower', vmin=450, vmax=600)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show how we can write the average image out to disk."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outfile = 'stacked_M13_blue.fits'\n",
    "\n",
    "hdu = fits.PrimaryHDU(final_image)\n",
    "hdu.writeto(outfile, overwrite=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}