Reduction_CCDproc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Reduction package for 1.3m DFOT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ccdproc is an Astropy affiliated package for basic data reductions of CCD images. It provides the essential tools for the cleaning of CCD images."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Steps to follow:\n",
    "1. Slicing the images taken in kinetic mode to single frames.\n",
    "2. Updating the Julian Dates (JD) for each sliced image to the header for time information.\n",
    "3. Alignment of the images using ASTROALIGN\n",
    "4. Cleaning the images\n",
    "          \n",
    "     a) Bias correction.\n",
    "     b) Flat correction.\n",
    "     c) Cosmic Ray removal\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following packages are necessary to import:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib as mp\n",
    "import ccdproc,os,sys,time,random\n",
    "from astropy import units as u\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import fits\n",
    "from astropy.time import Time\n",
    "from glob import glob\n",
    "#from astroquery.astrometry_net import AstrometryNet\n",
    "from astropy.wcs import WCS\n",
    "import astroalign as aa\n",
    "from astropy.nddata import CCDData\n",
    "from astropy.stats import sigma_clipped_stats, SigmaClip\n",
    "from astropy.visualization import ImageNormalize, LogStretch\n",
    "from matplotlib.ticker import LogLocator\n",
    "from astropy.stats import SigmaClip, mad_std"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1) 2) Slicing the images from datacube & Updating the J.D. for each sliced image using header information "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we are slicing the images from the datacube (i.e. that are taken in kinetic mode) and \n",
    "updating their time in JD (Julian Date) by taking information from the header.\n",
    "\n",
    "    INPUT:\n",
    "\n",
    "    path: The directory where the images are kept (string)\n",
    "    filename: The first few characters and the extension of the images (string). Example:\n",
    "    j0946*fits, HD1104*fits etc.\n",
    "\n",
    "\n",
    "    OUTPUT:\n",
    "\n",
    "    1. The images sliced in the case of kinetic mode images.\n",
    "    2. The Julian date updated to the header of the fits files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0. 0. 0. 0. 0. 0. 0. 0. 0.]\n",
      "kinetic mode image with no. of files: 4.0\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_14_17_sliced_0.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_14_17_sliced_1.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_14_17_sliced_2.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_14_17_sliced_3.fits has been written\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: A 'JD' keyword already exists in this header.  Inserting duplicate keyword. [astropy.io.fits.header]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kinetic mode image with no. of files: 4.0\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_19_22_sliced_0.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_19_22_sliced_1.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_19_22_sliced_2.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_19_22_sliced_3.fits has been written\n",
      "kinetic mode image with no. of files: 3.0\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_23_25_sliced_0.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_23_25_sliced_1.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_23_25_sliced_2.fits has been written\n",
      "kinetic mode image with no. of files: 5.0\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_2_6_sliced_0.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_2_6_sliced_1.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_2_6_sliced_2.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_2_6_sliced_3.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_2_6_sliced_4.fits has been written\n",
      "kinetic mode image with no. of files: 5.0\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_8_12_sliced_0.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_8_12_sliced_1.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_8_12_sliced_2.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_8_12_sliced_3.fits has been written\n",
      "/home/vivek/phot_dfot/dimple/photometry/data/j0947_i_8_12_sliced_4.fits has been written\n"
     ]
    }
   ],
   "source": [
    "path='/home/vivek/phot_dfot/dimple/photometry/data/'\n",
    "filename='j0947_i*fits'\n",
    "\n",
    "files=sorted(glob(os.path.join(path,filename)))\n",
    "nof=np.zeros(len(files))\n",
    "for i in range(0,len(files)):\n",
    "    data=fits.open(files[i])\n",
    "    header=data[0].header\n",
    "    image=data[0].data\n",
    "    k=np.shape(image)\n",
    "    nof[i]=k[0]\n",
    "\n",
    "    check_header=header['ACQMODE']\n",
    "\n",
    "    if (check_header=='Single Scan'):\n",
    "        jd_up=image\n",
    "        time=header['DATE']\n",
    "        t=Time(time,format='isot',scale='utc')\n",
    "        time_jd=t.jd\n",
    "        header.insert(15,('JD',time_jd))\n",
    "        files[i]\n",
    "        mod_file_1=files[i].replace('.fits','')\n",
    "        fits.writeto(mod_file_1+'_sliced'+'.fits',jd_up,header,overwrite=True)\n",
    "\n",
    "        #print(files[i],t.jd,t.mjd,'single scan image')\n",
    "\n",
    "\n",
    "\n",
    "    elif (check_header=='Kinetics'):\n",
    "\n",
    "        print('kinetic mode image with no. of files:',nof[i])\n",
    "\n",
    "        name_of_file=files[i]\n",
    "        mod_file=name_of_file.replace('.fits','')\n",
    "        time=header['DATE']\n",
    "        t=Time(time,format='isot',scale='utc')\n",
    "        t_jd=t.jd\n",
    "        temp=int(nof[i])\n",
    "        mod_jd=np.zeros(temp)\n",
    "        exp_time=header['ACT']/86400  # for the 'day' from seconds calculation.\n",
    "        mod_jd[0]=t_jd   # first frame time from the header and the next ones from the 'ACT'\n",
    "        for j in range(1,temp):\n",
    "            mod_jd[j]=mod_jd[j-1]+exp_time\n",
    "\n",
    "\n",
    "        for k in range(0,len(mod_jd)):\n",
    "            mod_header=header\n",
    "            sliced_image=image[k]\n",
    "            time_jd=mod_jd[k]\n",
    "            mod_header.insert(15,('JD',time_jd))\n",
    "            fits.writeto(mod_file+'_sliced_%g'%k+'.fits',sliced_image,mod_header,overwrite=True)\n",
    "            print(mod_file+'_sliced_%g'%k+'.fits has been written')\n",
    "            continue\n",
    "\n",
    "\n",
    "    #os.system(\"mv *sliced* cleaned/\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Aligning the images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Astro-Align package is used to alighn the images\n",
    "\n",
    "\n",
    "INPUT:\n",
    "\n",
    "    path: The directory where the images are kept (string)\n",
    "    filename: The first few characters and the extension of the images (string). Example:\n",
    "    j0946*fits, HD1104*fits etc.\n",
    "    ref_image: Reference image which will be used to align the images.\n",
    "\n",
    "OUTPUT:\n",
    "\n",
    "    aligned images.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No. 0 done\n",
      "No. 1 done\n",
      "No. 2 done\n",
      "No. 3 done\n",
      "No. 4 done\n",
      "No. 5 done\n",
      "No. 6 done\n",
      "No. 7 done\n",
      "No. 8 done\n",
      "No. 9 done\n",
      "No. 10 done\n",
      "No. 11 done\n",
      "No. 12 done\n",
      "No. 13 done\n",
      "No. 14 done\n",
      "No. 15 done\n",
      "No. 16 done\n",
      "No. 17 done\n",
      "No. 18 done\n",
      "No. 19 done\n",
      "No. 20 done\n",
      "No. 21 done\n",
      "No. 22 done\n",
      "No. 23 done\n",
      "No. 24 done\n"
     ]
    }
   ],
   "source": [
    "filename=\"*sliced*fits\"\n",
    "ref_image=\"j0947_i_1_sliced.fits\"\n",
    "nfiles=sorted(glob(path+filename))\n",
    "image_data=fits.open(path+ref_image)\n",
    "reference_image=image_data[0].data\n",
    "for i in range(len(nfiles)):\n",
    "    image_data=fits.open(nfiles[i])\n",
    "    source_image=image_data[0].data\n",
    "    header=image_data[0].header\n",
    "    image_aligned,footprint=aa.register(source_image,reference_image)\n",
    "\n",
    "    aligned_file=nfiles[i].replace('.fits','')\n",
    "    fits.writeto(aligned_file+'_aligned'+'.fits',image_aligned,header,overwrite=True)\n",
    "\n",
    "    print('No. %i alignment done'%i)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cleaning the images"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This module is meant for cleaning the images. The tasks to be included are: bias correction, flat correction, trimming, overscan as well as the cosmic ray removal from the science cases. (For the time we are skipping the overscan and trimming part.)\n",
    "\n",
    "INPUT:\n",
    "\n",
    "    path: The directory where the images are kept (string)\n",
    "    filename: The first few characters and the extension of the images (string). Example:\n",
    "    j0946*fits, HD1104*fits etc.\n",
    "\n",
    "OUTPUT:\n",
    "\n",
    "\n",
    "    cleaned images in the new directory: path/cleaned"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# a) BIAS correction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]\n",
      "Master bias generated\n",
      " Mean and median of the masterbias:  515.4987641731899 515.5\n"
     ]
    }
   ],
   "source": [
    "gain = 2 * u.electron / u.adu  # gain and readout noise are properties of the CCD and will change for different CCDs.\n",
    "readnoise = 7.5 * u.electron\n",
    "\n",
    "#ra=input('Enter the RA of the source:   ')\n",
    "#dec=input('Enter the DEC of the source: ')\n",
    "\n",
    "\n",
    "bias_files = sorted(glob(os.path.join(path,'bias*.fits')))\n",
    "biaslist = []\n",
    "for i in range (0,len(bias_files)):\n",
    "    data= ccdproc.CCDData.read(bias_files[i],unit='adu')\n",
    "        #data = ccdproc.create_deviation(data, gain=gain, readnoise=readnoise)\n",
    "        #data= data-(data.uncertainty.array)\n",
    "    biaslist.append(data)\n",
    "    \n",
    "masterbias = ccdproc.combine(biaslist,method='average',sigma_clip=True, sigma_clip_low_thresh=5, sigma_clip_high_thresh=5,\n",
    "                             sigma_clip_func=np.ma.median, sigma_clip_dev_func=mad_std)\n",
    "masterbias.write('masterbias.fits', overwrite=True)\n",
    "mbias=ccdproc.CCDData.read('masterbias.fits',unit='adu')\n",
    "print('Master bias generated')\n",
    "print(\" Mean and median of the masterbias: \",np.mean(masterbias), np.median(masterbias))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# b) Flat correction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: using the unit adu passed to the FITS reader instead of the unit adu in the FITS file. [astropy.nddata.ccddata]\n",
      "Master flat generated\n",
      " Mean and median of the masterflat:  0.9992749321566833 0.9999320144129445\n"
     ]
    }
   ],
   "source": [
    "flat_files=sorted(glob(os.path.join(dir,'flat*.fits')))\n",
    "flatlist = []\n",
    "for j in range(0,len(flat_files)):\n",
    "    flat=ccdproc.CCDData.read(flat_files[j],unit='adu')\n",
    "    flat_bias_removed=ccdproc.subtract_bias(flat,masterbias)\n",
    "    flatlist.append(flat_bias_removed)\n",
    "\n",
    "    def inv_median(a):\n",
    "        return 1 / np.median(a)\n",
    "\n",
    "masterflat = ccdproc.combine(flatlist,method='median', scale=inv_median,\n",
    "                                 sigma_clip=True, sigma_clip_low_thresh=5, sigma_clip_high_thresh=5,\n",
    "                                 sigma_clip_func=np.ma.median, sigma_clip_dev_func=mad_std)\n",
    "masterflat.write('masterflat.fits', overwrite=True)\n",
    "mflat=ccdproc.CCDData.read('masterflat.fits',unit='adu')\n",
    "print('Master flat generated')\n",
    "print(\" Mean and median of the masterflat: \",np.mean(masterflat), np.median(masterflat))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# c) Cosmic Ray removal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename=\"*aligned*fits\"\n",
    "file_names = sorted(glob(os.path.join(dir,filename)))\n",
    "for i in range(0,len(file_names)):\n",
    "    image=ccdproc.CCDData.read(file_names[i],unit='adu')\n",
    "    header=fits.getheader(file_names[i],0)\n",
    "    bias_subtracted = ccdproc.subtract_bias(image, masterbias)\n",
    "    flat_corrected = ccdproc.flat_correct(bias_subtracted, masterflat)\n",
    "    cr_cleaned = ccdproc.cosmicray_lacosmic(flat_corrected,readnoise=7.5, sigclip=5,satlevel=65535,niter=20,cleantype='meanmask',gain_apply=True)\n",
    "        #print('Cosmic rays removed')\n",
    "    clean_file=file_names[i].replace('.fits','')\n",
    "    fits.writeto(clean_file+'_cleaned.fits',cr_cleaned,header,overwrite=True)\n",
    "    print('Image no-%i has been cleaned'%i)"
   ]
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}