notebook_GMOS_a03_12amp_process_science_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Processing GMOS-N and GMOS-S science data\n",
    "\n",
    "Performs overscan correction, trimming, bias subtraction, flat-field correction, and cosmic ray removal on raw Gemini science data downloaded from the Gemini archive.  Uses raw bias data downloaded from the Gemini archive and processed flatfield files created by <tt>notebook_GMOS_a01_12amp_process_raw_twilight_flats.ipynb</tt> and <tt>notebook_GMOS_a02_12amp_combine_twilight_flatfield_images.ipynb</tt>.\n",
    "\n",
    "## Before running notebook:\n",
    "\n",
    "- Install <tt>fpack</tt> and <tt>funpack</tt> if not already installed (https://heasarc.gsfc.nasa.gov/fitsio/fpack/); change paths to local installation location in functions <tt>compress_file_fpack</tt> and <tt>decompress_file_fpack</tt> in this notebook<br><br>\n",
    "\n",
    "- Requires <tt>cosmics_py3.py</tt> to be present in same directory as notebook; uses the L.A.Cosmic cosmic ray removal algorithm developed by Pieter van Dokkum (https://arxiv.org/abs/astro-ph/0108003)<br><br>\n",
    "\n",
    "- If not already done, run <tt>notebook_GMOS_a01_12amp_process_raw_twilight_flats.ipynb</tt> and <tt>notebook_GMOS_a02_12amp_combine_twilight_flatfield_images</tt> to prepare flatfield images closest to the given night; place these files in the main directory where filenames should follow the format of <tt>nYYYYMMDD.instr.twiskyflat.?.chipN.fits.fz</tt>, where \"<tt>YYYYMMDD</tt>\" is the date associated with the median-combined flatfield data, \"<tt>instr</tt>\" indicates the instrument used (\"<tt>gmosn</tt>\" for GMOS-N and \"<tt>gmoss</tt>\" for GMOS-S), \"<tt>?</tt>\" is the filter associated with that flatfield file, and \"<tt>N</tt>\" is the chip number of that flatfield file (i.e., 1, 2, or 3).<br><br>\n",
    "\n",
    "- Download science data for a particular night from https://archive.gemini.edu/searchform; unpack archive, separate individual observations by filter, and place data files into subdirectories named <tt>rawfits_?</tt> in the main data directory where \"<tt>?</tt>\" is the filter associated with that subdirectory<br><br>\n",
    "\n",
    "- Go to https://archive.gemini.edu/searchform; specify UTC Date of corresponding flatfield data, \"GMOS-N\" or \"GMOS-S\" for Instrument, \"BIAS\" as Obs. Type, \"2x2\" for Binning, and \"Full Frame\" for ROI; search for bias frames corresponding to each night of flatfield data; download data and place in a subdirectory named <tt>rawfits_bias</tt> in the main data directory<br><br>\n",
    "\n",
    "- Run <tt>process_gemini_science_directory(data_path,observatory,date_id)</tt> where <tt>data_path</tt> is the main data directory, observatory is either \"<tt>GeminiN</tt>\" or \"<tt>GeminiS</tt>, and <tt>date_id</tt> is the date of observations in YYYYMMDD format\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import required packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "from astropy.io import fits\n",
    "from astropy.io.fits import getheader\n",
    "from astropy.modeling import models\n",
    "from astropy import nddata\n",
    "from astropy import units as u\n",
    "import ccdproc as cp\n",
    "from ccdproc import CCDData\n",
    "import numpy as np\n",
    "import math, glob, subprocess\n",
    "import cosmics_py3\n",
    "import datetime\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define file-handling functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_directory(path):\n",
    "    if not os.path.exists(path):\n",
    "        os.mkdir(path)\n",
    "    else:\n",
    "        print('Directory {:s} already exists.'.format(path))\n",
    "    return None\n",
    "\n",
    "def decompress_file_bzip2(filename):\n",
    "    cmd = ['bzip2','-d',filename]\n",
    "    process = subprocess.call(cmd)\n",
    "    return None\n",
    "\n",
    "def compress_file_fpack(filename):\n",
    "    process = subprocess.call(['/Users/hhsieh/Astro/tools/cfitsio/fpack',filename])\n",
    "    os.remove(filename)\n",
    "    return None\n",
    "\n",
    "def decompress_file_fpack(filename):\n",
    "    process = subprocess.call(['/Users/hhsieh/Astro/tools/cfitsio/funpack',filename])\n",
    "    os.remove(filename)\n",
    "    return None\n",
    "\n",
    "def decompress_directory_bzip2(file_path):\n",
    "    print('>>> Starting decompression of bz2 files in {:s}...'.format(file_path))\n",
    "    os.chdir(file_path)\n",
    "    for bz2_file in sorted(glob.glob('*.bz2')):\n",
    "        decompress_file_bzip2(bz2_file)\n",
    "    print('>>> Decompression of bz2 files complete.')\n",
    "    return None\n",
    "\n",
    "def compress_directory_fpack(file_path):\n",
    "    print('>>> Starting fpack compression of FITS files in {:s}...'.format(file_path))\n",
    "    os.chdir(file_path)\n",
    "    for fits_file in glob.glob('*.fits'):\n",
    "        compress_file_fpack(fits_file)\n",
    "    print('>>> fpack compression of FITS files complete.')\n",
    "    return None\n",
    "    \n",
    "def decompress_directory_fpack(file_path):\n",
    "    print('>>> Starting decompression of fz files in {:s}...'.format(file_path))\n",
    "    os.chdir(file_path)\n",
    "    for fz_file in glob.glob('*.fz'):\n",
    "        decompress_file_fpack(fz_file)\n",
    "    print('>>> Decompression of fz files complete.')\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define function to retrieve date of FITS files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_date_id():\n",
    "    # obtain YYYYMMDD date from filename of GMOS FITS files being processed (regardless of current compression state)\n",
    "    fits_date = ''\n",
    "    if len(glob.glob('*.fits')) > 0:\n",
    "        for raw_mef_file in sorted(glob.glob('*.fits')):\n",
    "            fits_date = raw_mef_file[1:9]\n",
    "    elif len(glob.glob('*.fits.fz')) > 0:\n",
    "        for raw_mef_file in sorted(glob.glob('*.fits.fz')):\n",
    "            fits_date = raw_mef_file[1:9]\n",
    "    elif len(glob.glob('*.fits.bz2')) > 0:\n",
    "        for raw_mef_file in sorted(glob.glob('*.fits.bz2')):\n",
    "            fits_date = raw_mef_file[1:9]\n",
    "    return fits_date\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define functions to construct image statistics files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_stats_files(fits_date,imgtype):\n",
    "    for extid in range(2,10):\n",
    "        ext_stats = 'n{:s}.{:s}.{:02d}.stats'.format(fits_date,imgtype,extid+1)\n",
    "        with open(ext_stats,'w') as of:\n",
    "            of.write('# Extension {:02d}              NPIX        MEAN     STDDEV         MIN         MAX\\n'.format(extid+1))\n",
    "    return None\n",
    "        \n",
    "def compile_stats_files(fits_date,imgtype):\n",
    "    print('Starting image statistics output...')\n",
    "    output_stats_filename = 'n{:s}.{:s}.stats'.format(fits_date,imgtype)\n",
    "    with open(output_stats_filename,'w') as output_stats_file:\n",
    "        for stats_filename in glob.glob('*.??.stats'):\n",
    "            with open(stats_filename,'r') as stats_file:\n",
    "                for line in stats_file:\n",
    "                    output_stats_file.write(line)\n",
    "            os.remove(stats_filename)\n",
    "    print('>>> Image statistics output complete.')\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define header editing functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_fits_header_line(fits_file,param_name,param_value,param_comment):\n",
    "    # Add line to FITS header\n",
    "    with fits.open(fits_file) as hdulist:\n",
    "        data = hdulist[0].data\n",
    "        hdr = hdulist[0].header\n",
    "    hdr[param_name] = (param_value,param_comment)\n",
    "    fits.writeto(fits_file,data,hdr,overwrite=True)\n",
    "    return None\n",
    "\n",
    "def fix_gain_rdnoise_gmosn_hamamatsu():\n",
    "    # Fix gain and read noise values in FITS headers\n",
    "    for fits_file in sorted(glob.glob('*.???.??.fits')):\n",
    "        with fits.open(fits_file) as hdulist:\n",
    "            data = hdulist[0].data\n",
    "            hdr = hdulist[0].header\n",
    "        hdr['GAIN'] = 1.63\n",
    "        hdr['RDNOISE'] = 4.14\n",
    "        fits.writeto(fits_file,data,hdr,overwrite=True)\n",
    "    print('>>> Gain and read noise updated in image headers.')\n",
    "    return None\n",
    "\n",
    "def fix_gain_rdnoise_gmoss_hamamatsu():\n",
    "    # Fix gain and read noise values in FITS headers\n",
    "    for fits_file in sorted(glob.glob('*.???.??.fits')):\n",
    "        with fits.open(fits_file) as hdulist:\n",
    "            data = hdulist[0].data\n",
    "            hdr = hdulist[0].header\n",
    "        hdr['GAIN'] = 1.83\n",
    "        hdr['RDNOISE'] = 3.98\n",
    "        fits.writeto(fits_file,data,hdr,overwrite=True)\n",
    "    print('>>> Gain and read noise updated in image headers.')\n",
    "    return None\n",
    "\n",
    "def fix_gmos_header_info(instr_id,filter_id):\n",
    "    # Fix header information for PDS submission\n",
    "    for fits_file in sorted(glob.glob('*.otzfc.fits')):\n",
    "        with fits.open(fits_file) as hdulist:\n",
    "            hdr,data = hdulist[0].header,hdulist[0].data\n",
    "\n",
    "        # Change keywords for CDi_j parameters (conflict with other WCS keywords)\n",
    "        cd1_1,cd1_1_comment = hdr['CD1_1'],hdr.comments['CD1_1']\n",
    "        cd1_2,cd1_2_comment = hdr['CD1_2'],hdr.comments['CD1_2']\n",
    "        cd2_1,cd2_1_comment = hdr['CD2_1'],hdr.comments['CD2_1']\n",
    "        cd2_2,cd2_2_comment = hdr['CD2_2'],hdr.comments['CD2_2']\n",
    "        hdr['x_CD1_1'] = (cd1_1,cd1_1_comment)\n",
    "        hdr['x_CD1_2'] = (cd1_2,cd1_2_comment)\n",
    "        hdr['x_CD2_1'] = (cd2_1,cd2_1_comment)\n",
    "        hdr['x_CD2_2'] = (cd2_2,cd2_2_comment)\n",
    "        del hdr['CD1_1']\n",
    "        del hdr['CD1_2']\n",
    "        del hdr['CD2_1']\n",
    "        del hdr['CD2_2']\n",
    "        \n",
    "        # Change keyword for EPOCH parameter\n",
    "        target_epoch,target_epoch_comment = hdr['EPOCH'],hdr.comments['EPOCH']\n",
    "        hdr['TRGEPOCH'] = (target_epoch,target_epoch_comment)\n",
    "        del hdr['EPOCH']\n",
    "        \n",
    "        # Remove RADVEL parameter (not correctly populated)\n",
    "        del hdr['RADVEL']\n",
    "        \n",
    "        # Change EXPTIME value to integer format\n",
    "        exptime = int(round(hdr['EXPTIME']))\n",
    "        hdr['EXPTIME'] = exptime\n",
    "        \n",
    "        # Add central filter wavelength keyword\n",
    "        if filter_id == 'g':   hdr['CENTWAVE'] = 475\n",
    "        elif filter_id == 'r': hdr['CENTWAVE'] = 630\n",
    "        elif filter_id == 'i': hdr['CENTWAVE'] = 780\n",
    "        elif filter_id == 'z': hdr['CENTWAVE'] = 876\n",
    "        \n",
    "        # Remove checksums\n",
    "        if 'CHECKSUM' in hdr: del hdr['CHECKSUM']\n",
    "        if 'DATASUM'  in hdr: del hdr['DATASUM']\n",
    "        \n",
    "        fits.writeto(fits_file,data,hdr,overwrite=True,checksum=True)\n",
    "        \n",
    "    print('>>> Header information fixed for PDS submission.')\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define function for splitting multi-extension FITS files into individual extension files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def split_extensions(imgtype):\n",
    "    # split 12-element GMOS multi-extension FITS files into individual elements\n",
    "    # return: dimensions of extensions; also writes individual extension files to working directory\n",
    "    for raw_mef_file in glob.glob('*.fits'):\n",
    "        with fits.open(raw_mef_file) as hdulist:\n",
    "            print('Splitting {:s}...'.format(raw_mef_file))\n",
    "            fits_date = raw_mef_file[1:9]\n",
    "            fits_id   = raw_mef_file[11:14]\n",
    "            hdr1 = getheader(raw_mef_file,0)\n",
    "            for extid in range(2,10):\n",
    "                extension = hdulist[extid+1].data\n",
    "                hdr2 = hdulist[extid+1].header\n",
    "                radecsys = hdr2['RADECSYS']\n",
    "                hdr2['RADESYSA'] = radecsys\n",
    "                del hdr2['RADECSYS']\n",
    "                imstat_npix = extension.size\n",
    "                imstat_min,imstat_max    = np.min(extension),np.max(extension)\n",
    "                imstat_mean,imstat_stdev = np.mean(extension),np.std(extension)\n",
    "                ext_filename = 'n' + fits_date + '.' + fits_id + '.{:02d}.fits'.format(extid+1)\n",
    "                outputfilename = 'n{:s}.{:s}.{:02d}.stats'.format(fits_date,imgtype,extid+1)\n",
    "                with open(outputfilename,'a') as of:\n",
    "                    of.write('{:s}   {:>8d}    {:>8.2f}    {:>7.2f}    {:>8.2f}    {:>8.2f}\\n'.format(ext_filename,imstat_npix,imstat_mean,imstat_stdev,imstat_min,imstat_max))\n",
    "                dimensions = extension.shape\n",
    "                dimension1 = dimensions[0]\n",
    "                dimension2 = dimensions[1]\n",
    "                hdr = hdr1 + hdr2\n",
    "                fits.writeto(ext_filename,extension,hdr)\n",
    "    print('>>> Multi-extension fits file splitting complete.')\n",
    "    return dimension1, dimension2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define function for performing overscan correction and trimming"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def overscan_and_trim():\n",
    "    # Do overscan correction and trim images in current directory\n",
    "    print('Starting overscan correction and trimming...')\n",
    "    for fits_file in sorted(glob.glob('*.???.??.fits')):\n",
    "        ot_fits_file = fits_file[0:16] + '.ot.fits'\n",
    "        fits_data = CCDData.read(fits_file,unit=u.adu)\n",
    "        file_ext = fits_file[14:16]\n",
    "        # Overscan correction\n",
    "        if ((file_ext == '03') or (file_ext=='05') or (file_ext=='07') or (file_ext=='09')):\n",
    "            o_fits_data  = cp.subtract_overscan(fits_data, fits_section='[257:288,1:2112]', overscan_axis=1, add_keyword={'oscansub': True, 'calstat': 'O'}, model=models.Polynomial1D(1))\n",
    "        else:\n",
    "            o_fits_data  = cp.subtract_overscan(fits_data, fits_section='[1:32,1:2112]', overscan_axis=1, add_keyword={'oscansub': True, 'calstat': 'O'}, model=models.Polynomial1D(1))\n",
    "        # Image trimming\n",
    "        if ((file_ext == '03') or (file_ext=='07')):\n",
    "            ot_fits_data = cp.trim_image(o_fits_data,fits_section='[1:256,65:2110]',add_keyword={'trimmed': True, 'calstat': 'OT'})\n",
    "        elif ((file_ext == '04') or (file_ext == '08')):\n",
    "            ot_fits_data = cp.trim_image(o_fits_data,fits_section='[33:282,65:2110]',add_keyword={'trimmed': True, 'calstat': 'OT'})\n",
    "        elif ((file_ext == '05') or (file_ext =='09')):\n",
    "            ot_fits_data = cp.trim_image(o_fits_data,fits_section='[7:256,65:2110]',add_keyword={'trimmed': True, 'calstat': 'OT'})\n",
    "        else: # file_ext == '02' or file_ext == '10'\n",
    "            ot_fits_data = cp.trim_image(o_fits_data,fits_section='[33:288,65:2110]',add_keyword={'trimmed': True, 'calstat': 'OT'})\n",
    "        ot_fits_data.write(ot_fits_file)\n",
    "        os.remove(fits_file)\n",
    "    print('>>> Overscan correction and trimming complete.')\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define bias correction functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_bias_data(path_rawbias,instr_id,date_id,reduction_log):\n",
    "    print('\\n{:s} - Processing bias-frame data in {:s}...'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),path_rawbias))\n",
    "    reduction_log.write('\\n{:s} - Processing bias-frame science data in {:s}...\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),path_rawbias))\n",
    "    \n",
    "    # decompress raw bias fits files\n",
    "    print('\\n{:s} - Decompressing, filtering, and splitting bias files...'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S')))\n",
    "    reduction_log.write('{:s} - Decompressing, filtering, and splitting bias files...\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S')))\n",
    "    os.chdir(path_rawbias)\n",
    "    decompress_directory_bzip2(path_rawbias)       # decompress downloaded GMOS files\n",
    "    decompress_directory_fpack(path_rawbias)       # decompress previously fpacked GMOS files\n",
    "    create_stats_files(date_id,'bias')       # create files to record image statistics\n",
    "    split_extensions('bias')                 # split bias MEFs into individual FITS files\n",
    "    if instr_id == 'gmosn':\n",
    "        fix_gain_rdnoise_gmosn_hamamatsu()\n",
    "    elif instr_id == 'gmoss':\n",
    "        fix_gain_rdnoise_gmoss_hamamatsu()\n",
    "    compile_stats_files(date_id,'bias')      # collect image stats of bias files together and delete individual files\n",
    "    \n",
    "    # perform overscan corrections and trim data for bias frames\n",
    "    os.chdir(path_rawbias)\n",
    "    overscan_and_trim()\n",
    "    bias_median_combine(date_id)\n",
    "    \n",
    "    print('\\n{:s} - Processing bias-frame data in {:s} complete.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),path_rawbias))\n",
    "    reduction_log.write('\\n{:s} - Processing bias-frame science data in {:s} complete.\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),path_rawbias))\n",
    "    return reduction_log\n",
    "\n",
    "def bias_median_combine(dateid):\n",
    "    print('Starting median combination of bias frames...')\n",
    "    for extid in range(2,10):\n",
    "        bias_list = []\n",
    "        ext_file_format = '*.{:02d}.ot.fits'.format(extid+1)\n",
    "        output_filename = 'n' + dateid + '.bias.{:02d}.fits'.format(extid+1)\n",
    "        for fits_file in sorted(glob.glob(ext_file_format)):\n",
    "            fits_data = CCDData.read(fits_file)\n",
    "            bias_list.append(fits_data)\n",
    "            os.remove(fits_file)\n",
    "        master_bias = cp.combine(bias_list,method='median')\n",
    "        master_bias.write(output_filename)\n",
    "    print('>>> Median combination of bias frames complete.')\n",
    "    return None\n",
    "\n",
    "def bias_correct(dateid,cwd_raw_bias):\n",
    "    print('Starting bias subtraction...')\n",
    "    for extid in range(2,10):\n",
    "        bias_filename = cwd_raw_bias + 'n' + dateid + '.bias.{:02d}.fits'.format(extid+1)\n",
    "        ext_file_format = '*.{:02d}.ot.fits'.format(extid+1)\n",
    "        for fits_file in sorted(glob.glob(ext_file_format)):\n",
    "            fits_data = CCDData.read(fits_file)\n",
    "            bias_data = CCDData.read(bias_filename)\n",
    "            fits_date_imageid = fits_file[0:13]\n",
    "            output_filename = fits_date_imageid + '.{:02d}.otz.fits'.format(extid+1)\n",
    "            bias_corrected_data = cp.subtract_bias(fits_data,bias_data,add_keyword={'zerocorr':True,'calstat':'OTZ'})\n",
    "            bias_corrected_data.write(output_filename)\n",
    "            os.remove(fits_file)\n",
    "    print('>>> Bias subtraction complete.')\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define function to concatenate amp sections for each GMOS chip"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def concatenate_gmos_amps():\n",
    "    print('>>> Starting adjacent amp area concatenation...')\n",
    "    for fits_file in sorted(glob.glob('*.03.otz.fits')):\n",
    "        file_prefix = fits_file[0:13]\n",
    "        ext1_filename = file_prefix + '.03.otz.fits'\n",
    "        ext2_filename = file_prefix + '.04.otz.fits'\n",
    "        hdr = getheader(ext1_filename,0)\n",
    "        with fits.open(ext1_filename) as hdul1, fits.open(ext2_filename) as hdul2:\n",
    "            ext1_data = hdul1[0].data\n",
    "            ext2_data = hdul2[0].data\n",
    "        chip1_data = np.hstack((ext1_data,ext2_data))\n",
    "        output_filename = file_prefix + '.chip1.otz.fits'\n",
    "        fits.writeto(output_filename,chip1_data,hdr)\n",
    "        os.remove(ext1_filename)\n",
    "        os.remove(ext2_filename)\n",
    "\n",
    "        ext1_filename = file_prefix + '.05.otz.fits'\n",
    "        ext2_filename = file_prefix + '.06.otz.fits'\n",
    "        ext3_filename = file_prefix + '.07.otz.fits'\n",
    "        ext4_filename = file_prefix + '.08.otz.fits'\n",
    "        hdr = getheader(ext1_filename,0)\n",
    "        with fits.open(ext1_filename) as hdul1, fits.open(ext2_filename) as hdul2, fits.open(ext3_filename) as hdul3, fits.open(ext4_filename) as hdul4:        \n",
    "            ext1_data = hdul1[0].data\n",
    "            ext2_data = hdul2[0].data\n",
    "            ext3_data = hdul3[0].data\n",
    "            ext4_data = hdul4[0].data\n",
    "        chip1_data = np.hstack((ext1_data,ext2_data,ext3_data,ext4_data))\n",
    "        output_filename = file_prefix + '.chip2.otz.fits'\n",
    "        fits.writeto(output_filename,chip1_data,hdr)\n",
    "        os.remove(ext1_filename)\n",
    "        os.remove(ext2_filename)\n",
    "        os.remove(ext3_filename)\n",
    "        os.remove(ext4_filename)\n",
    "\n",
    "        ext1_filename = file_prefix + '.09.otz.fits'\n",
    "        ext2_filename = file_prefix + '.10.otz.fits'\n",
    "        hdr = getheader(ext1_filename,0)\n",
    "        with fits.open(ext1_filename) as hdul1, fits.open(ext2_filename) as hdul2:\n",
    "            ext1_data = hdul1[0].data\n",
    "            ext2_data = hdul2[0].data\n",
    "        chip1_data = np.hstack((ext1_data,ext2_data))\n",
    "        output_filename = file_prefix + '.chip3.otz.fits'\n",
    "        fits.writeto(output_filename,chip1_data,hdr)\n",
    "        os.remove(ext1_filename)\n",
    "        os.remove(ext2_filename)\n",
    "    print('>>> Adjacent amp area concatenation complete.')\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define functions for flatfield correction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_flatfield_files(data_path,instr_id,filter_id):\n",
    "    ff_chip1_filename,ff_chip2_filename,ff_chip3_filename = 'none','none','none'\n",
    "    os.chdir(data_path)\n",
    "    decompress_directory_fpack(data_path)\n",
    "    flatfield_filename_pattern = 'n*.' + instr_id + '.twiskyflat.' + filter_id + '*.fits'\n",
    "    if len(glob.glob(flatfield_filename_pattern)) != 3:\n",
    "        print('Incorrect number of flatfield files found.')\n",
    "    else:\n",
    "        for flatfield_file in sorted(glob.glob(flatfield_filename_pattern)):\n",
    "            if flatfield_file[-10:-5] == 'chip1': ff_chip1_filename = flatfield_file\n",
    "            if flatfield_file[-10:-5] == 'chip2': ff_chip2_filename = flatfield_file\n",
    "            if flatfield_file[-10:-5] == 'chip3': ff_chip3_filename = flatfield_file\n",
    "    return ff_chip1_filename,ff_chip2_filename,ff_chip3_filename\n",
    "\n",
    "def flatfield_correction(data_path,flatfield_chip1_file,flatfield_chip2_file,flatfield_chip3_file):\n",
    "    print('Starting flatfield correction...')\n",
    "    for fits_file in sorted(glob.glob('*.chip1.otz.fits')):\n",
    "        fits_data = CCDData.read(fits_file)\n",
    "        flat_data = CCDData.read(data_path+flatfield_chip1_file,unit=u.adu)\n",
    "        fits_date_imageid = fits_file[0:13]\n",
    "        output_filename = fits_date_imageid + '.chip1.otzf.fits'\n",
    "        flat_corrected_data = cp.flat_correct(fits_data,flat_data,add_keyword={'flatcorr': True, 'calstat': 'OTZF'})\n",
    "        flat_corrected_data.write(output_filename)\n",
    "        add_fits_header_line(output_filename,'BIAS_1',fits_date_imageid+'.bias.03.fits','1st bias frame used')\n",
    "        add_fits_header_line(output_filename,'BIAS_2',fits_date_imageid+'.bias.04.fits','2nd bias frame used')\n",
    "        add_fits_header_line(output_filename,'BIAS_3','n/a','3rd bias frame used')\n",
    "        add_fits_header_line(output_filename,'BIAS_4','n/a','4th bias frame used')\n",
    "        add_fits_header_line(output_filename,'FLATUSED',flatfield_chip1_file,'Flatfield frame used')\n",
    "        os.remove(fits_file)\n",
    "\n",
    "    for fits_file in sorted(glob.glob('*.chip2.otz.fits')):\n",
    "        fits_data = CCDData.read(fits_file)\n",
    "        flat_data = CCDData.read(data_path+flatfield_chip2_file,unit=u.adu)\n",
    "        fits_date_imageid = fits_file[0:13]\n",
    "        output_filename = fits_date_imageid + '.chip2.otzf.fits'\n",
    "        flat_corrected_data = cp.flat_correct(fits_data,flat_data,add_keyword={'flatcorr': True, 'calstat': 'OTZF'})\n",
    "        flat_corrected_data.write(output_filename)\n",
    "        add_fits_header_line(output_filename,'BIAS_1',fits_date_imageid+'.bias.05.fits','1st bias frame used')\n",
    "        add_fits_header_line(output_filename,'BIAS_2',fits_date_imageid+'.bias.06.fits','2nd bias frame used')\n",
    "        add_fits_header_line(output_filename,'BIAS_3',fits_date_imageid+'.bias.07.fits','3rd bias frame used')\n",
    "        add_fits_header_line(output_filename,'BIAS_4',fits_date_imageid+'.bias.08.fits','4th bias frame used')\n",
    "        add_fits_header_line(output_filename,'FLATUSED',flatfield_chip2_file,'Flatfield frame used')\n",
    "        os.remove(fits_file)\n",
    "\n",
    "    for fits_file in sorted(glob.glob('*.chip3.otz.fits')):\n",
    "        fits_data = CCDData.read(fits_file)\n",
    "        flat_data = CCDData.read(data_path+flatfield_chip3_file,unit=u.adu)\n",
    "        fits_date_imageid = fits_file[0:13]\n",
    "        output_filename = fits_date_imageid + '.chip3.otzf.fits'\n",
    "        flat_corrected_data = cp.flat_correct(fits_data,flat_data,add_keyword={'flatcorr': True, 'calstat': 'OTZF'})\n",
    "        flat_corrected_data.write(output_filename)\n",
    "        add_fits_header_line(output_filename,'BIAS_1',fits_date_imageid+'.bias.09.fits','1st bias frame used')\n",
    "        add_fits_header_line(output_filename,'BIAS_2',fits_date_imageid+'.bias.10.fits','2nd bias frame used')\n",
    "        add_fits_header_line(output_filename,'BIAS_3','n/a','3rd bias frame used')\n",
    "        add_fits_header_line(output_filename,'BIAS_4','n/a','4th bias frame used')\n",
    "        add_fits_header_line(output_filename,'FLATUSED',flatfield_chip3_file,'Flatfield frame used')\n",
    "        os.remove(fits_file)\n",
    "\n",
    "    print('>>> Flatfield correction complete.')\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define function to perform cosmic ray detection and removal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cosmicray_cleaning():\n",
    "    print('Starting cosmic ray correction...')\n",
    "    for fits_file in sorted(glob.glob('*.chip1.otzf.fits')):\n",
    "        print('Cosmic ray cleaning for {:s}'.format(fits_file))\n",
    "        file_prefix = fits_file[0:19]\n",
    "        hdr = getheader(fits_file,0)\n",
    "        fits_gain = hdr['gain']\n",
    "        fits_readnoise = hdr['rdnoise']\n",
    "        array,header = cosmics_py3.fromfits(fits_file,verbose=False)\n",
    "        c = cosmics_py3.cosmicsimage(array,gain=fits_gain,readnoise=fits_readnoise,sigclip=5.0,sigfrac=0.3,objlim=5.0,verbose=False)\n",
    "        c.run(maxiter=4)\n",
    "        output_file = file_prefix + '.otzfc.fits'\n",
    "        cosmics_py3.tofits(output_file,c.cleanarray,header,verbose=False)\n",
    "        add_fits_header_line(output_file,'CRCORR',True,'Cosmic ray removal performed? (T/F)')\n",
    "        add_fits_header_line(output_file,'CALSTAT','OTZFC','Image calibration status')\n",
    "        add_fits_header_line(output_file,'CRR_GAIN',fits_gain,'CR removal parameter: gain')\n",
    "        add_fits_header_line(output_file,'CRR_NOIS',fits_readnoise,'CR removal parameter: read noise')\n",
    "        add_fits_header_line(output_file,'CRR_SIGC',5.0,'CR removal parameter: sigclip')\n",
    "        add_fits_header_line(output_file,'CRR_SIGF',0.3,'CR removal parameter: sigfrac')\n",
    "        add_fits_header_line(output_file,'CRR_OBJL',5.0,'CR removal parameter: objlim')\n",
    "        os.remove(fits_file)\n",
    "\n",
    "    for fits_file in sorted(glob.glob('*.chip2.otzf.fits')):\n",
    "        print('Cosmic ray cleaning for {:s}'.format(fits_file))\n",
    "        file_prefix = fits_file[0:19]\n",
    "        hdr = getheader(fits_file,0)\n",
    "        fits_gain = hdr['gain']\n",
    "        fits_readnoise = hdr['rdnoise']\n",
    "        array,header = cosmics_py3.fromfits(fits_file,verbose=False)\n",
    "        c = cosmics_py3.cosmicsimage(array,gain=fits_gain,readnoise=fits_readnoise,sigclip=5.0,sigfrac=0.3,objlim=5.0,verbose=False)\n",
    "        c.run(maxiter=4)\n",
    "        output_file = file_prefix + '.otzfc.fits'\n",
    "        cosmics_py3.tofits(output_file,c.cleanarray,header,verbose=False)\n",
    "        add_fits_header_line(output_file,'CRCORR',True,'Cosmic ray removal performed? (T/F)')\n",
    "        add_fits_header_line(output_file,'CALSTAT','OTZFC','Image calibration status')\n",
    "        add_fits_header_line(output_file,'CRR_GAIN',fits_gain,'CR removal parameter: gain')\n",
    "        add_fits_header_line(output_file,'CRR_NOIS',fits_readnoise,'CR removal parameter: read noise')\n",
    "        add_fits_header_line(output_file,'CRR_SIGC',5.0,'CR removal parameter: sigclip')\n",
    "        add_fits_header_line(output_file,'CRR_SIGF',0.3,'CR removal parameter: sigfrac')\n",
    "        add_fits_header_line(output_file,'CRR_OBJL',5.0,'CR removal parameter: objlim')\n",
    "        os.remove(fits_file)\n",
    "\n",
    "    for fits_file in sorted(glob.glob('*.chip3.otzf.fits')):\n",
    "        print('Cosmic ray cleaning for {:s}'.format(fits_file))\n",
    "        file_prefix = fits_file[0:19]\n",
    "        hdr = getheader(fits_file,0)\n",
    "        fits_gain = hdr['gain']\n",
    "        fits_readnoise = hdr['rdnoise']\n",
    "        array,header = cosmics_py3.fromfits(fits_file,verbose=False)\n",
    "        c = cosmics_py3.cosmicsimage(array,gain=fits_gain,readnoise=fits_readnoise,sigclip=5.0,sigfrac=0.3,objlim=5.0,verbose=False)\n",
    "        c.run(maxiter=4)\n",
    "        output_file = file_prefix + '.otzfc.fits'\n",
    "        cosmics_py3.tofits(output_file,c.cleanarray,header,verbose=False)\n",
    "        add_fits_header_line(output_file,'CRCORR',True,'Cosmic ray removal performed? (T/F)')\n",
    "        add_fits_header_line(output_file,'CALSTAT','OTZFC','Image calibration status')\n",
    "        add_fits_header_line(output_file,'CRR_GAIN',fits_gain,'CR removal parameter: gain')\n",
    "        add_fits_header_line(output_file,'CRR_NOIS',fits_readnoise,'CR removal parameter: read noise')\n",
    "        add_fits_header_line(output_file,'CRR_SIGC',5.0,'CR removal parameter: sigclip')\n",
    "        add_fits_header_line(output_file,'CRR_SIGF',0.3,'CR removal parameter: sigfrac')\n",
    "        add_fits_header_line(output_file,'CRR_OBJL',5.0,'CR removal parameter: objlim')\n",
    "        os.remove(fits_file)\n",
    "        \n",
    "    print('>>> Cosmic ray correction complete.')\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define functions for sorting final output files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clean_up_science_data(path_rawfits,path_procfits,path_stats):\n",
    "    # Clean up data\n",
    "    move_reduced_science_data(path_rawfits,path_procfits,path_stats)\n",
    "    compress_directory_fpack(path_rawfits)\n",
    "    return None\n",
    "\n",
    "def clean_up_common_data(data_path,path_rawbias,path_procfits,path_stats):\n",
    "    # Clean up data\n",
    "    move_reduced_bias_data(path_rawbias,path_procfits,path_stats)\n",
    "    compress_directory_fpack(data_path)\n",
    "    compress_directory_fpack(path_rawbias)\n",
    "    if os.path.exists(path_procfits):\n",
    "        compress_directory_fpack(path_procfits)\n",
    "    return None\n",
    "\n",
    "def move_reduced_science_data(rawfitsdir,flatdir,path_stats):\n",
    "    os.chdir(rawfitsdir)\n",
    "    for fits_file in sorted(glob.glob('*.otzfc.fits')):\n",
    "        os.rename(rawfitsdir+fits_file,flatdir+fits_file)\n",
    "    for stats_file in sorted(glob.glob('*.stats')):\n",
    "        os.rename(rawfitsdir+stats_file,path_stats+stats_file)\n",
    "    print('>>> Processed data moved to output directory.')\n",
    "    return None\n",
    "\n",
    "def move_reduced_bias_data(rawbiasdir,flatdir,path_stats):\n",
    "    os.chdir(rawbiasdir)\n",
    "    if os.path.exists(flatdir):\n",
    "        for fits_file in sorted(glob.glob('*.bias.??.fits')):\n",
    "            os.rename(rawbiasdir+fits_file,flatdir+fits_file)\n",
    "    for stats_file in sorted(glob.glob('*.stats')):\n",
    "        os.rename(rawbiasdir+stats_file,path_stats+stats_file)\n",
    "    print('>>> Processed bias data moved to output directory.')\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define functions for science data processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_science_data(data_path,path_rawfits,path_rawbias,path_procfits,path_stats,instr_id,filter_id,reduction_log):\n",
    "    print('\\n{:s} - Processing {:s}-band science data in {:s}...'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),filter_id,data_path))\n",
    "    reduction_log.write('\\n{:s} - Processing {:s}-band science data in {:s}...\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),filter_id,data_path))\n",
    "    \n",
    "    flatfield_chip1_filename,flatfield_chip2_filename,flatfield_chip3_filename = find_flatfield_files(data_path,instr_id,filter_id)\n",
    "    \n",
    "    if flatfield_chip1_filename != 'none' and flatfield_chip2_filename != 'none' and flatfield_chip3_filename != 'none':\n",
    "        if not os.path.isdir(path_procfits): os.mkdir(path_procfits)\n",
    "        if not os.path.isdir(path_stats):    os.mkdir(path_stats)\n",
    "        # Process raw science fits files\n",
    "        print('\\n{:s} - Decompressing and splitting data files...'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S')))\n",
    "        reduction_log.write('{:s} - Decompressing and splitting data files...\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S')))\n",
    "        os.chdir(path_rawfits)\n",
    "        decompress_directory_bzip2(path_rawfits)       # decompress downloaded GMOS files\n",
    "        decompress_directory_fpack(path_rawfits)       # decompress previously fpacked GMOS files\n",
    "        date_id = get_date_id()\n",
    "        create_stats_files(date_id,'scidata.{:s}'.format(filter_id))    # create files to record image statistics\n",
    "        dim1, dim2 = split_extensions('scidata') # split science MEFs into individual FITS files\n",
    "        if instr_id == 'gmosn':\n",
    "            fix_gain_rdnoise_gmosn_hamamatsu()\n",
    "        elif instr_id == 'gmoss':\n",
    "            fix_gain_rdnoise_gmoss_hamamatsu()\n",
    "        compile_stats_files(date_id,'scidata.{:s}'.format(filter_id))   # collect image stats of science files together and delete individual files\n",
    "            \n",
    "        # Process science frames\n",
    "        os.chdir(path_rawfits)\n",
    "        overscan_and_trim()                                                    # Perform overscan correction and trim data\n",
    "        bias_correct(date_id,path_rawbias)                                      # Perform bias correction\n",
    "        concatenate_gmos_amps()                                                # Join science frames\n",
    "        flatfield_correction(data_path,flatfield_chip1_filename,flatfield_chip2_filename,flatfield_chip3_filename)  # Do flatfield corrections\n",
    "        cosmicray_cleaning()                                                   # Do cosmic ray cleaning\n",
    "    \n",
    "        print('{:s} - Processing science data in {:s} complete.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),data_path))\n",
    "        reduction_log.write('{:s} - Processing science data in {:s} complete.\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),data_path))\n",
    "            \n",
    "    else:\n",
    "        print('{:s} - Appropriate flatfield images for {:s} not found.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),path_rawfits))        \n",
    "        reduction_log.write('{:s} - Appropriate flatfield images for {:s} not found.\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),path_rawfits))        \n",
    "            \n",
    "    return reduction_log\n",
    "\n",
    "def process_gemini_science_directory(data_path,observatory,date_id):\n",
    "    if observatory == 'GeminiN':   instr_id = 'gmosn'\n",
    "    elif observatory == 'GeminiS': instr_id = 'gmoss'\n",
    "    else: instr_id = 'not recognized'\n",
    "    \n",
    "    if instr_id == 'gmosn' or instr_id == 'gmoss':\n",
    "        os.chdir(data_path)\n",
    "        reduction_logfile = data_path + 'log_reduction_{:s}_{:s}_{:s}.txt'.format(date_id,instr_id,datetime.datetime.today().strftime('%Y%m%d_%H%M%S'))\n",
    "        path_rawbias  = data_path + 'rawfits_bias/'\n",
    "        path_procfits = data_path + 'flatfits_pyt/'\n",
    "        path_stats    = data_path + 'stats/'\n",
    "        with open(reduction_logfile,'w') as reduction_log:\n",
    "            if os.path.exists(path_procfits):\n",
    "                print('{:s} - Science data in {:s} already processed.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),data_path))\n",
    "            elif not os.path.exists(path_rawbias):\n",
    "                print('{:s} - Bias directory in {:s} not found.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),data_path))\n",
    "                reduction_log.write('{:s} - Bias directory in {:s} not found.\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),data_path))\n",
    "            elif len(os.listdir(path_rawbias)) == 0:\n",
    "                print('{:s} - Bias directory in {:s} is empty.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),data_path))\n",
    "                reduction_log.write('{:s} - Bias directory in {:s} is empty.\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),data_path))\n",
    "            elif len(glob.glob('rawfits_?/')) == 0:\n",
    "                print('{:s} - Raw science data folder(s) in {:s} not found.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),data_path))\n",
    "                reduction_log.write('{:s} - Raw science data folder(s) in {:s} not found.\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),data_path))\n",
    "            else:\n",
    "                print('{:s} - Starting processing of data in {:s}...'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),data_path))\n",
    "                reduction_log.write('{:s} - Starting processing of data in {:s}...\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),data_path))\n",
    "                for dir_rawfits in sorted(glob.glob('rawfits_?/')):\n",
    "                    os.chdir(data_path + dir_rawfits)\n",
    "                    date_id = get_date_id()\n",
    "                reduction_log = process_bias_data(path_rawbias,instr_id,date_id,reduction_log)\n",
    "                os.chdir(data_path)\n",
    "                for dir_rawfits in sorted(glob.glob('rawfits_?/')):\n",
    "                    filter_id = dir_rawfits[-2:-1]\n",
    "                    path_rawfits  = data_path + dir_rawfits\n",
    "                    reduction_log = process_science_data(data_path,path_rawfits,path_rawbias,path_procfits,path_stats,instr_id,filter_id,reduction_log)\n",
    "                    clean_up_science_data(path_rawfits,path_procfits,path_stats)\n",
    "                clean_up_common_data(data_path,path_rawbias,path_procfits,path_stats)\n",
    "                print('{:s} - Processing of data in {:s} complete.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),data_path))\n",
    "                reduction_log.write('{:s} - Processing of data in {:s} complete.\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),data_path))\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform science data processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# date_id: date of observations in YYYYMMDD format\n",
    "# observatory: 'GeminiN' or 'GeminiS'\n",
    "\n",
    "process_gemini_science_directory(data_path,observatory,date_id)\n"
   ]
  }
 ],
 "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.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}