notebook_GMOS_a01_12amp_process_raw_twilight_flats.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Initial processing of GMOS-N and GMOS-S flatfield files\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> below<br><br>\n",
    "\n",
    "- Go to https://archive.gemini.edu/searchform; specify UTC Date range of interest, \"GMOS-N\" or \"GMOS-S\" for Instrument, \"Twilight\" as Target Name, \"2x2\" for Binning, and \"Full Frame\" for ROI, and specify filter of interest; search for flatfield images; download data<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<br><br>\n",
    "\n",
    "- Unpack compressed archives into a new directory, putting raw flatfield data in a subdirectory called \"rawfits_twiskyflat\" and raw bias data in a subdirectory called \"rawfits_bias\"<br><br>\n",
    "\n",
    "- Run <tt>process_flatfield_directory(data_path)</tt> where data_path is location of directory containing \"<tt>rawfits_twiskyflat</tt>\" and \"<tt>rawfits_bias</tt>\" subdirectories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "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 glob, subprocess\n",
    "import datetime\n"
   ]
  },
  {
   "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",
    "    process = subprocess.call(['bzip2','-d',filename])\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 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",
    "\n",
    "def get_date_id():\n",
    "    # obtain YYYYMMDD date from filename of GMOS FITS files being processed\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(2,10):\n",
    "        ext_stats = 'n' + fits_date + '.{:02d}.stats'.format(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,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",
    "    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",
    "        os.remove(stats_filename)\n",
    "    output_stats_file.close()\n",
    "    print('>>> Image statistics output complete.')\n",
    "    return None\n",
    "\n",
    "def split_extensions():\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' + fits_date + '.{:02d}.stats'.format(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,dimension2 = dimensions[0],dimensions[1]\n",
    "                hdr = hdr1 + hdr2\n",
    "                fits.writeto(ext_filename,extension,hdr)\n",
    "    return dimension1,dimension2\n",
    "\n",
    "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 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",
    "\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 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 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",
    "    \n",
    "def concatenate_gmos_amps():\n",
    "    print('>>> Starting adjacent amp area concatenation...')\n",
    "    for fits_file in 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",
    "        \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",
    "            \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",
    "        \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",
    "        \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",
    "\n",
    "\n",
    "def process_flatfield_files(path_rawfits,path_rawbias,path_procfits):\n",
    "    \n",
    "    print('\\n')\n",
    "    create_directory(path_procfits)\n",
    "    \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",
    "    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)                   #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",
    "    # 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",
    "    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)                   #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",
    "    overscan_and_trim()\n",
    "    bias_median_combine(date_id)\n",
    "\n",
    "    # perform overscan corrections and trim data for science frames\n",
    "    os.chdir(path_rawfits)\n",
    "    overscan_and_trim()\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",
    "        os.rename(fits_file,path_procfits+fits_file)\n",
    "    os.chdir(path_rawbias)\n",
    "    for fits_file in glob.glob('*.bias.??.fits'):\n",
    "        os.rename(fits_file,path_procfits+fits_file)\n",
    "\n",
    "    compress_directory_fpack(path_rawfits)\n",
    "    compress_directory_fpack(path_rawbias)\n",
    "    compress_directory_fpack(path_procfits)\n",
    "        \n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_flatfield_directory(data_path):\n",
    "    os.chdir(data_path)\n",
    "    path_rawfits  = data_path + 'rawfits_twiskyflat/'\n",
    "    path_rawbias  = data_path + 'rawfits_bias/'\n",
    "    path_procfits = data_path + 'procfits_twiskyflat/'\n",
    "    \n",
    "    if os.path.exists(path_procfits):    \n",
    "        print('{:s} - Data in {:s} already processed.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "    elif not os.path.exists(path_rawfits):\n",
    "        print('{:s} - Raw flatfield directory ({:s}) not found.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "    elif not os.path.exists(path_rawbias):\n",
    "        print('{:s} - Bias directory ({:s}) not found.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\n",
    "    elif len(os.listdir(path_rawfits)) == 0:\n",
    "        print('{:s} - Raw flatfield directory in {:s} is empty.'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),datadir))\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'),datadir))\n",
    "    else:\n",
    "        process_flatfield_files(path_rawfits,path_rawbias,path_procfits)\n",
    "\n",
    "    return None\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "process_flatfield_directory(data_path)\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
}