notebook_GMOS_01_6amp_process_raw_twilight_flats.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.6.11 | packaged by conda-forge | (default, Aug  5 2020, 20:19:23) \n",
      "[GCC Clang 10.0.1 ]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import os, sys\n",
    "sys.path.append('/Users/hhsieh/anaconda3/envs/astroconda/lib/python3.6/site-packages')\n",
    "sys.path.append('/opt/local/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages')\n",
    "print(sys.version)\n",
    "print('')\n",
    "import matplotlib.pyplot as plt\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, ImageFileCollection\n",
    "import numpy as np\n",
    "import glob, os, bz2, subprocess\n",
    "import sqlite3\n",
    "from sqlite3 import Error\n",
    "import datetime\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Done.\n"
     ]
    }
   ],
   "source": [
    "def create_directory(path):\n",
    "    if not os.path.exists(path):\n",
    "        os.mkdir(path)\n",
    "        print('Directory {:s} successfully created.'.format(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",
    "    fpack_command  = '/Users/hhsieh/Astro/tools/cfitsio/fpack'\n",
    "    delete_command = '/bin/rm'\n",
    "    fpack_cmd = [fpack_command,filename]\n",
    "    process = subprocess.call(fpack_cmd)\n",
    "    delete_cmd = [delete_command,filename]\n",
    "    process = subprocess.call(delete_cmd)\n",
    "    return None\n",
    "\n",
    "def decompress_file_fpack(filename):\n",
    "    funpack_command = '/Users/hhsieh/Astro/tools/cfitsio/funpack'\n",
    "    delete_command  = '/bin/rm'\n",
    "    funpack_cmd = [funpack_command,filename]\n",
    "    process = subprocess.call(funpack_cmd)\n",
    "    delete_cmd = [delete_command,filename]\n",
    "    process = subprocess.call(delete_cmd)\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 glob.glob('*.bz2'):\n",
    "        decompress_file_bzip2(bz2_file)\n",
    "        #cmd = [bzip1,bzip2,bz2_file]\n",
    "        #process = subprocess.call(cmd)\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",
    "    #fpack_file = '/Users/hhsieh/Astro/tools/cfitsio/fpack'\n",
    "    #delete_file = '/bin/rm'\n",
    "    os.chdir(file_path)\n",
    "    for fits_file in glob.glob('*.fits'):\n",
    "        compress_file_fpack(fits_file)\n",
    "        #fpack_cmd = [fpack_file,fits_file]\n",
    "        #process = subprocess.call(fpack_cmd)\n",
    "        #delete_cmd = [delete_file,fits_file]\n",
    "        #process = subprocess.call(delete_cmd)\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",
    "    #funpack_file = '/Users/hhsieh/Astro/tools/cfitsio/funpack'\n",
    "    #delete_file = '/bin/rm'\n",
    "    os.chdir(file_path)\n",
    "    for fz_file in glob.glob('*.fz'):\n",
    "        decompress_file_fpack(fz_file)\n",
    "        #funpack_cmd = [funpack_file,fz_file]\n",
    "        #process = subprocess.call(funpack_cmd)\n",
    "        #delete_cmd = [delete_file,fz_file]\n",
    "        #process = subprocess.call(delete_cmd)\n",
    "    print('>>> Decompression of fz files complete.')\n",
    "    return None\n",
    "\n",
    "\n",
    "#def decompress_data_bzip2():\n",
    "#    print('>>> Starting decompression of bz2 files...')\n",
    "#    bzip1 = 'bzip2'\n",
    "#    bzip2 = '-d'\n",
    "#    for bz2_file in glob.glob('*.bz2'):\n",
    "#        cmd = [bzip1,bzip2,bz2_file]\n",
    "#        process = subprocess.call(cmd)\n",
    "#    print('>>> Decompression of bz2 files complete.')\n",
    "#\n",
    "#def decompress_data_fpack():\n",
    "#    print('>>> Starting decompression of fz files...')\n",
    "#    funpack_file = '/Users/hhsieh/Astro/tools/cfitsio/funpack'\n",
    "#    delete_file = '/bin/rm'\n",
    "#    for fz_file in glob.glob('*.fz'):\n",
    "#        funpack_cmd = [funpack_file,fz_file]\n",
    "#        process = subprocess.call(funpack_cmd)\n",
    "#        delete_cmd = [delete_file,fz_file]\n",
    "#        process = subprocess.call(delete_cmd)\n",
    "#    print('>>> Decompression of fz files complete.')\n",
    "#\n",
    "#def compress_data_fpack():\n",
    "#    print('>>> Starting fpack compression of fits files...')\n",
    "#    fpack_file = '/Users/hhsieh/Astro/tools/cfitsio/fpack'\n",
    "#    delete_file = '/bin/rm'\n",
    "#    for fits_file in glob.glob('*.fits'):\n",
    "#        fpack_cmd = [fpack_file,fits_file]\n",
    "#        process = subprocess.call(fpack_cmd)\n",
    "#        delete_cmd = [delete_file,fits_file]\n",
    "#        process = subprocess.call(delete_cmd)\n",
    "#    print('>>> Compression of fits files complete.')\n",
    "\n",
    "def get_date_id():\n",
    "    # obtain YYYYMMDD date from filename of GMOS FITS files being processed\n",
    "    # return: YYYYMMDD date id\n",
    "    for raw_mef_file in glob.glob('*.fits'):\n",
    "        fits_date = raw_mef_file[1:9]\n",
    "    return fits_date\n",
    "\n",
    "def create_stats_files(fits_date):\n",
    "    for extid in range(0,6):\n",
    "        ext_stats = 'n' + fits_date + '.{:02d}.stats'.format(extid+1)\n",
    "        outputfile = open(ext_stats,'w')\n",
    "        outputfile.write('# Extension {:02d}               NPIX        MEAN     STDDEV         MIN         MAX\\n'.format(extid+1))\n",
    "        outputfile.close()\n",
    "    return None\n",
    "        \n",
    "def compile_stats_files(fits_date,outputfile_path):\n",
    "    print('Starting image statistics output...')\n",
    "    output_stats_filename = outputfile_path + 'n' + fits_date + '.stats'\n",
    "    output_stats_file = open(output_stats_filename,'w')\n",
    "    cmd1 = '/bin/rm'\n",
    "    for stats_filename in glob.glob('*.??.stats'):\n",
    "        stats_file = open(stats_filename,'r')\n",
    "        for line in stats_file:\n",
    "            output_stats_file.write(line)\n",
    "        stats_file.close()\n",
    "        cmd = [cmd1,stats_filename]\n",
    "        process = subprocess.call(cmd)\n",
    "    output_stats_file.close()\n",
    "    print('>>> Image statistics output complete.')\n",
    "    return None\n",
    "\n",
    "def remove_test_images():\n",
    "    # remove test images from downloaded Gemini data tarball\n",
    "    delete_cmd = '/bin/rm'\n",
    "    for fits_file in glob.glob('*.fits'):\n",
    "        hdulist = fits.open(fits_file)\n",
    "        hdr = hdulist[1].header\n",
    "        naxis2 = hdr['naxis2']\n",
    "        if(naxis2 < 200):\n",
    "            delete_file = [delete_cmd,fits_file]\n",
    "            print('Removing {:s}...'.format(fits_file))\n",
    "            process = subprocess.call(delete_file)\n",
    "    return None\n",
    "\n",
    "def split_extensions():\n",
    "    # split 6-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",
    "        hdulist = fits.open(raw_mef_file)\n",
    "        print('Splitting {:s}...'.format(raw_mef_file))\n",
    "        fits_date = raw_mef_file[1:9]\n",
    "        fits_id   = raw_mef_file[10:14]\n",
    "        hdr1 = getheader(raw_mef_file,0)\n",
    "        for extid in range(0,6):\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 = np.min(extension)\n",
    "            imstat_max = np.max(extension)\n",
    "            imstat_mean = np.mean(extension)\n",
    "            imstat_stdev = np.std(extension)\n",
    "            ext_filename = 'n' + fits_date + '.' + fits_id + '.{:02d}.fits'.format(extid+1)\n",
    "            outputfilename = 'n' + fits_date + '.{:02d}.stats'.format(extid+1)\n",
    "            outputfile = open(outputfilename,'a')\n",
    "            outputfile.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",
    "            outputfile.close()\n",
    "            dimensions = extension.shape\n",
    "            dimension1 = dimensions[0]\n",
    "            dimension2 = dimensions[1]\n",
    "            hdr = hdr1 + hdr2\n",
    "            fits.writeto(ext_filename,extension,hdr)\n",
    "        hdulist.close()\n",
    "    return dimension1, dimension2\n",
    "    print('{:s} - >>> Multi-extension fits file splitting complete.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S')))\n",
    "    return None\n",
    "    \n",
    "def remove_wrong_dimensions(dimension1,dimension2):\n",
    "    # identify bias files with correct binning and remove others\n",
    "    # return: none; retains bias files with correct binning and removes others\n",
    "    for raw_mef_file in glob.glob('*.fits'):\n",
    "        hdulist = fits.open(raw_mef_file)\n",
    "        extension01 = hdulist[1].data\n",
    "        bias_dimensions = extension01.shape\n",
    "        bias_dimension1 = bias_dimensions[0]\n",
    "        bias_dimension2 = bias_dimensions[1]\n",
    "        hdulist.close()\n",
    "        delete_file = '/bin/rm'\n",
    "        if (bias_dimension1 != dimension1) or (bias_dimension2 != dimension2):\n",
    "            delete_cmd = [delete_file,raw_mef_file]\n",
    "            process = subprocess.call(delete_cmd)\n",
    "    print('{:s} - >>> Bias frame selection complete.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S')))\n",
    "    return None\n",
    "\n",
    "def overscan_and_trim_gmosn():\n",
    "    # Do overscan correction and trim images in current directory\n",
    "    print('Starting overscan correction and trimming...')\n",
    "    for fits_file in glob.glob('*.????.??.fits'):\n",
    "        ot_fits_file = fits_file[0:17] + '.ot.fits'\n",
    "        fits_data = CCDData.read(fits_file,unit=u.adu)\n",
    "        file_ext = fits_file[15:17]\n",
    "        # Overscan correction\n",
    "        if ((file_ext == '01') or (file_ext == '03')):\n",
    "            o_fits_data  = cp.subtract_overscan(fits_data, fits_section='[1:32,1:2304]', overscan_axis=1, add_keyword={'oscansub': True, 'calstat': 'O'}, model=models.Polynomial1D(1))\n",
    "        else: # if file_ext == 4 or 5\n",
    "            o_fits_data  = cp.subtract_overscan(fits_data, fits_section='[513:544,1:2304]', overscan_axis=1, add_keyword={'oscansub': True, 'calstat': 'O'}, model=models.Polynomial1D(1))\n",
    "        # Image trimming for Ext #s 1,3,4,5\n",
    "        if ((file_ext == '01') or (file_ext == '03')):\n",
    "            ot_fits_data = cp.trim_image(o_fits_data,fits_section='[33:543,3:2260]',add_keyword={'trimmed': True, 'calstat': 'OT'})\n",
    "            ot_fits_data.write(ot_fits_file)\n",
    "        elif ((file_ext == '04') or (file_ext == '05')):\n",
    "            ot_fits_data = cp.trim_image(o_fits_data,fits_section='[2:512,3:2260]',add_keyword={'trimmed': True, 'calstat': 'OT'})\n",
    "            ot_fits_data.write(ot_fits_file)\n",
    "        delete_file = '/bin/rm'\n",
    "        delete_cmd = [delete_file,fits_file]\n",
    "        process = subprocess.call(delete_cmd)\n",
    "    print('{:s} - >>> Overscan correction and trimming complete.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S')))\n",
    "    return None\n",
    "\n",
    "def bias_median_combine(dateid):\n",
    "    print('Starting median combination of bias frames...')\n",
    "    for extid in range(0,6):\n",
    "        if((extid == 0) or (extid == 2) or (extid == 3) or (extid == 4)):\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 glob.glob(ext_file_format):\n",
    "                fits_data = CCDData.read(fits_file)\n",
    "                bias_list.append(fits_data)\n",
    "                delete_file = '/bin/rm'\n",
    "                delete_cmd = [delete_file,fits_file]\n",
    "                process = subprocess.call(delete_cmd)\n",
    "            master_bias = cp.combine(bias_list,method='median')\n",
    "            master_bias.write(output_filename)\n",
    "    print('{:s} - >>> Median combination of bias frames complete.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S')))\n",
    "    return None\n",
    "\n",
    "def bias_correct(dateid,cwd_raw_bias):\n",
    "    print('Starting bias subtraction...')\n",
    "    delete_cmd = '/bin/rm'\n",
    "    for extid in range(0,6):\n",
    "        if((extid == 0) or (extid == 2) or (extid == 3) or (extid == 4)):\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 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:14]\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",
    "                delete_file = [delete_cmd,fits_file]\n",
    "                process = subprocess.call(delete_file)\n",
    "    print('{:s} - >>> Bias subtraction complete.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S')))\n",
    "    return None\n",
    "\n",
    "def concatenate_gmos_amps():\n",
    "    print('>>> Starting adjacent amp area concatenation...')\n",
    "    delete_cmd = '/bin/rm'\n",
    "    move_cmd = 'mv'\n",
    "    \n",
    "    for fits_file in glob.glob('*.01.otz.fits'):\n",
    "        file_prefix = fits_file[0:14]\n",
    "        output_filename = file_prefix + '.chip1.otz.fits'\n",
    "        move_file = [move_cmd,fits_file,output_filename]\n",
    "        process = subprocess.call(move_file)\n",
    "\n",
    "    for fits_file in glob.glob('*.03.otz.fits'):\n",
    "        file_prefix = fits_file[0:14]\n",
    "        ext1_filename = file_prefix + '.04.otz.fits'\n",
    "        ext2_filename = file_prefix + '.03.otz.fits'\n",
    "        hdr = getheader(ext1_filename,0)\n",
    "        hdulist1 = fits.open(ext1_filename)\n",
    "        ext1_data = hdulist1[0].data\n",
    "        hdulist2 = fits.open(ext2_filename)\n",
    "        ext2_data = hdulist2[0].data\n",
    "        hdulist1.close()\n",
    "        hdulist2.close()\n",
    "        chip1_data = np.hstack((ext1_data,ext2_data))\n",
    "        output_filename = file_prefix + '.chip2.otz.fits'\n",
    "        fits.writeto(output_filename,chip1_data,hdr)\n",
    "        delete_file = [delete_cmd,ext1_filename]\n",
    "        process = subprocess.call(delete_file)\n",
    "        delete_file = [delete_cmd,ext2_filename]\n",
    "        process = subprocess.call(delete_file)\n",
    "        \n",
    "    for fits_file in glob.glob('*.05.otz.fits'):\n",
    "        file_prefix = fits_file[0:14]\n",
    "        output_filename = file_prefix + '.chip3.otz.fits'\n",
    "        move_file = [move_cmd,fits_file,output_filename]\n",
    "        process = subprocess.call(move_file)\n",
    "        \n",
    "    print('{:s} - >>> Adjacent amp area concatenation complete.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S')))\n",
    "    return None\n",
    "\n",
    "\n",
    "def process_flatfield_files(path_rawfits,path_rawbias,path_procfits,instr_id):\n",
    "    \n",
    "    print('\\n')\n",
    "    create_directory(path_procfits)\n",
    "    \n",
    "    print('\\n{:s} - Decompressing and splitting data files...'.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",
    "    remove_test_images()\n",
    "    date_id = get_date_id()\n",
    "    create_stats_files(date_id)                   #create files to record image statistics of bias files\n",
    "    dim1,dim2 = split_extensions()                #split science MEFs into individual FITS files\n",
    "    compile_stats_files(date_id,path_procfits)    #collect image statistics of science files together and delete individual files\n",
    "\n",
    "    print('\\n{:s} - Decompressing, filtering, and splitting bias files...'.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",
    "    remove_test_images()\n",
    "    remove_wrong_dimensions(dim1,dim2)            #find and retain bias frames matching data file binning\n",
    "    create_stats_files(date_id)                   #create files to record image statistics of bias files\n",
    "    split_extensions()                            #split bias MEFs into individual FITS files\n",
    "    compile_stats_files(date_id+'.bias',path_procfits)          #collect image statistics 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",
    "    if instr_id == 'gmosn':\n",
    "        overscan_and_trim_gmosn()\n",
    "    if instr_id == 'gmoss':\n",
    "        overscan_and_trim_gmoss()\n",
    "    bias_median_combine(date_id)\n",
    "        \n",
    "    # perform overscan corrections and trim data for science frames\n",
    "    os.chdir(path_rawfits)\n",
    "    if instr_id == 'gmosn':\n",
    "        overscan_and_trim_gmosn()\n",
    "    if instr_id == 'gmoss':\n",
    "        overscan_and_trim_gmoss()\n",
    "    bias_correct(date_id,path_rawbias)\n",
    "\n",
    "    # join science frames\n",
    "    os.chdir(path_rawfits)\n",
    "    concatenate_gmos_amps()\n",
    "    \n",
    "    # move processed files to processed directory\n",
    "    os.chdir(path_rawfits)\n",
    "    for fits_file in glob.glob('*.chip?.otz.fits'):\n",
    "        cmd_mv = 'mv'\n",
    "        move_file_to_flat_dir = [cmd_mv,fits_file,path_procfits]\n",
    "        process = subprocess.call(move_file_to_flat_dir)\n",
    "    os.chdir(path_rawbias)\n",
    "    for fits_file in glob.glob('*.bias.??.fits'):\n",
    "        cmd_mv = 'mv'\n",
    "        move_file_to_flat_dir = [cmd_mv,fits_file,path_procfits]\n",
    "        process = subprocess.call(move_file_to_flat_dir)\n",
    "\n",
    "    compress_directory_fpack(path_rawfits)\n",
    "    compress_directory_fpack(path_rawbias)\n",
    "    compress_directory_fpack(path_procfits)\n",
    "        \n",
    "    return None\n",
    "        \n",
    "print('Done.')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_flatfield_directory(basewd):\n",
    "    os.chdir(basewd)\n",
    "    for datadir in glob.glob('ut*_gemini?_twiskyflat'):\n",
    "        path_rawfits  = basewd + datadir + '/rawfits_twiskyflat/'\n",
    "        path_rawbias  = basewd + datadir + '/rawfits_bias/'\n",
    "        path_procfits = basewd + datadir + '/procfits_twiskyflat/'\n",
    "    \n",
    "        instr_id = 'not recognized'\n",
    "        if datadir[-12:-11] == 'N': instr_id = 'gmosn'\n",
    "        if datadir[-12:-11] == 'S': instr_id = 'gmoss'\n",
    "    \n",
    "        if not os.path.exists(path_procfits):\n",
    "            if os.path.exists(path_rawfits):\n",
    "                if len(os.listdir(path_rawfits)) > 0:\n",
    "                    if os.path.exists(path_rawbias):\n",
    "                        if len(os.listdir(path_rawbias)) > 0:\n",
    "                            if instr_id != 'not recognized':\n",
    "                                if instr_id == 'gmosn':\n",
    "                                    process_flatfield_files(path_rawfits,path_rawbias,path_procfits,instr_id)\n",
    "                                else:\n",
    "                                    print('{:s} - GMOS-S reduction code not yet available.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S')))\n",
    "                            else:\n",
    "                                print('{:s} - Instrument for {:s} not recognized.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "                        else:\n",
    "                            print('{:s} - Bias directory in {:s} is empty.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "                    else:\n",
    "                        print('{:s} - Bias directory in {:s} not found.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "                else:\n",
    "                    print('{:s} - Flatfield directory in {:s} is empty.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "            else:\n",
    "                print('{:s} - Flatfield directory in {:s} not found.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "        else:\n",
    "            print('{:s} - Data in {:s} already processed.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "    return None\n",
    "\n",
    "# Set base working directory and subdirectory paths\n",
    "\n",
    "#basewd = '/volumes/Fantom12a/BackupData/gemini/flatfields_twilight/GMOS-N_TwilightFlats_20111015_20170202/GMOS-N_Twilight_Flats_g/'\n",
    "#basewd = '/volumes/Fantom12a/BackupData/gemini/flatfields_twilight/GMOS-N_TwilightFlats_20111015_20170202/GMOS-N_Twilight_Flats_r/'\n",
    "#basewd = '/volumes/Fantom12a/BackupData/gemini/flatfields_twilight/GMOS-N_TwilightFlats_20111015_20170202/GMOS-N_Twilight_Flats_i/'\n",
    "#basewd = '/volumes/Fantom12a/BackupData/gemini/flatfields_twilight/GMOS-S_TwilightFlats_20140601_20170301/GMOS-S_Twilight_Flats_g/'\n",
    "#basewd = '/volumes/Fantom12a/BackupData/gemini/flatfields_twilight/GMOS-S_TwilightFlats_20140601_20170301/GMOS-S_Twilight_Flats_r/'\n",
    "#basewd = '/volumes/Fantom12a/BackupData/gemini/flatfields_twilight/GMOS-S_TwilightFlats_20140601_20170301/GMOS-S_Twilight_Flats_i/'\n",
    "#rawfits_dir   = '/rawfits_twiskyflat/'\n",
    "#rawbias_dir   = '/rawfits_bias/'\n",
    "#procfits_dir  = '/procfits_twiskyflat/'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Directory /Users/hhsieh/Data/176P/ut20051229_gemini2_twiskyflat/ut20051229_geminiN_twiskyflat/procfits_twiskyflat/ successfully created.\n",
      "\n",
      "2022-07-15 22:49:44 - Decompressing and splitting data files...\n",
      ">>> Starting decompression of bz2 files in /Users/hhsieh/Data/176P/ut20051229_gemini2_twiskyflat/ut20051229_geminiN_twiskyflat/rawfits_twiskyflat/...\n",
      ">>> Decompression of bz2 files complete.\n",
      ">>> Starting decompression of fz files in /Users/hhsieh/Data/176P/ut20051229_gemini2_twiskyflat/ut20051229_geminiN_twiskyflat/rawfits_twiskyflat/...\n",
      ">>> Decompression of fz files complete.\n",
      "Removing N20051229S0329.fits...\n",
      "Splitting N20051229S0338.fits...\n"
     ]
    },
    {
     "ename": "IndexError",
     "evalue": "list index out of range",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-4-10e066bc523a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0;31m##process_flatfield_directory('/volumes/Fantom12a/BackupData/gemini/flatfields_twilight/GMOS-S_TwilightFlats_20140601_20150101/GMOS-S_Twilight_Flats_i/')\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 8\u001b[0;31m \u001b[0mprocess_flatfield_directory\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'/Users/hhsieh/Data/176P/ut20051229_gemini2_twiskyflat/'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     10\u001b[0m \u001b[0;31m#os.chdir(basewd)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-3-e01af29dd99f>\u001b[0m in \u001b[0;36mprocess_flatfield_directory\u001b[0;34m(basewd)\u001b[0m\n\u001b[1;32m     17\u001b[0m                             \u001b[0;32mif\u001b[0m \u001b[0minstr_id\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;34m'not recognized'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m                                 \u001b[0;32mif\u001b[0m \u001b[0minstr_id\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'gmosn'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 19\u001b[0;31m                                     \u001b[0mprocess_flatfield_files\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpath_rawfits\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpath_rawbias\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpath_procfits\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0minstr_id\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     20\u001b[0m                                 \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     21\u001b[0m                                     \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'{:s} - GMOS-S reduction code not yet available.'\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdatetime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdatetime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtoday\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstrftime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'%Y-%m-%d %H:%M:%S'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-2-ac7e93e76252>\u001b[0m in \u001b[0;36mprocess_flatfield_files\u001b[0;34m(path_rawfits, path_rawbias, path_procfits, instr_id)\u001b[0m\n\u001b[1;32m    310\u001b[0m     \u001b[0mdate_id\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_date_id\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    311\u001b[0m     \u001b[0mcreate_stats_files\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdate_id\u001b[0m\u001b[0;34m)\u001b[0m                   \u001b[0;31m#create files to record image statistics of bias files\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 312\u001b[0;31m     \u001b[0mdim1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mdim2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msplit_extensions\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m                \u001b[0;31m#split science MEFs into individual FITS files\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    313\u001b[0m     \u001b[0mcompile_stats_files\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdate_id\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpath_procfits\u001b[0m\u001b[0;34m)\u001b[0m    \u001b[0;31m#collect image statistics of science files together and delete individual files\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    314\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-2-ac7e93e76252>\u001b[0m in \u001b[0;36msplit_extensions\u001b[0;34m()\u001b[0m\n\u001b[1;32m    154\u001b[0m         \u001b[0mhdr1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgetheader\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mraw_mef_file\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    155\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mextid\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m6\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 156\u001b[0;31m             \u001b[0mextension\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mhdulist\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mextid\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    157\u001b[0m             \u001b[0mhdr2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mhdulist\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mextid\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mheader\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    158\u001b[0m             \u001b[0mradecsys\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mhdr2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'RADECSYS'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/astroconda/lib/python3.6/site-packages/astropy/io/fits/hdu/hdulist.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m    316\u001b[0m             \u001b[0;31m# Raise a more helpful IndexError if the file was not fully read.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    317\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_read_all\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 318\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    319\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    320\u001b[0m                 raise IndexError('HDU not found, possibly because the index '\n",
      "\u001b[0;32m~/anaconda3/envs/astroconda/lib/python3.6/site-packages/astropy/io/fits/hdu/hdulist.py\u001b[0m in \u001b[0;36m__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m    312\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    313\u001b[0m             return self._try_while_unread_hdus(super().__getitem__,\n\u001b[0;32m--> 314\u001b[0;31m                                                self._positive_index_of(key))\n\u001b[0m\u001b[1;32m    315\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mIndexError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    316\u001b[0m             \u001b[0;31m# Raise a more helpful IndexError if the file was not fully read.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/envs/astroconda/lib/python3.6/site-packages/astropy/io/fits/hdu/hdulist.py\u001b[0m in \u001b[0;36m_try_while_unread_hdus\u001b[0;34m(self, func, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1109\u001b[0m         \u001b[0;32mwhile\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1110\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1111\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1112\u001b[0m             \u001b[0;32mexcept\u001b[0m \u001b[0mException\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1113\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_read_next_hdu\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: list index out of range"
     ]
    }
   ],
   "source": [
    "#process_flatfield_directory('/volumes/Fantom12a/BackupData/gemini/flatfields_twilight/GMOS-N_TwilightFlats_20111015_20170202/GMOS-N_Twilight_Flats_g/')\n",
    "#process_flatfield_directory('/volumes/Fantom12a/BackupData/gemini/flatfields_twilight/GMOS-N_TwilightFlats_20111015_20170202/GMOS-N_Twilight_Flats_r/')\n",
    "#process_flatfield_directory('/volumes/Fantom12a/BackupData/gemini/flatfields_twilight/GMOS-N_TwilightFlats_20111015_20170202/GMOS-N_Twilight_Flats_i/')\n",
    "##process_flatfield_directory('/volumes/Fantom12a/BackupData/gemini/flatfields_twilight/GMOS-S_TwilightFlats_20140601_20150101/GMOS-S_Twilight_Flats_g/')\n",
    "##process_flatfield_directory('/volumes/Fantom12a/BackupData/gemini/flatfields_twilight/GMOS-S_TwilightFlats_20140601_20150101/GMOS-S_Twilight_Flats_r/')\n",
    "##process_flatfield_directory('/volumes/Fantom12a/BackupData/gemini/flatfields_twilight/GMOS-S_TwilightFlats_20140601_20150101/GMOS-S_Twilight_Flats_i/')\n",
    "\n",
    "process_flatfield_directory('/Users/hhsieh/Data/176P/ut20051229_gemini2_twiskyflat/')\n",
    "\n",
    "#os.chdir(basewd)\n",
    "#for datadir in glob.glob('ut*_gemini?_twiskyflat'):\n",
    "#    path_rawfits  = basewd + datadir + rawfits_dir\n",
    "#    path_rawbias  = basewd + datadir + rawbias_dir\n",
    "#    path_procfits = basewd + datadir + procfits_dir\n",
    "#    \n",
    "#    instr_id = 'not recognized'\n",
    "#    if datadir[-12:-11] == 'N': instr_id = 'gmosn'\n",
    "#    if datadir[-12:-11] == 'S': instr_id = 'gmoss'\n",
    "#    \n",
    "#    if not os.path.exists(path_procfits):\n",
    "#        if os.path.exists(path_rawfits):\n",
    "#            if len(os.listdir(path_rawfits)) > 0:\n",
    "#                if os.path.exists(path_rawbias):\n",
    "#                    if len(os.listdir(path_rawbias)) > 0:\n",
    "#                        if instr_id != 'not recognized':\n",
    "#                            if instr_id == 'gmosn':\n",
    "#                                process_flatfield_files(path_rawfits,path_rawbias,path_procfits,instr_id)\n",
    "#                            else:\n",
    "#                                print('{:s} - GMOS-S reduction code not yet available.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S')))\n",
    "#                        else:\n",
    "#                            print('{:s} - Instrument for {:s} not recognized.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "#                    else:\n",
    "#                        print('{:s} - Bias directory in {:s} is empty.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "#                else:\n",
    "#                    print('{:s} - Bias directory in {:s} not found.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "#            else:\n",
    "#                print('{:s} - Flatfield directory in {:s} is empty.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "#        else:\n",
    "#            print('{:s} - Flatfield directory in {:s} not found.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "#    else:\n",
    "#        print('{:s} - Data in {:s} already processed.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "\n",
    "print('Done.')\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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}