GradPak_skysub.py

#!/usr/bin/python
#
# History:
#     v1 - A. Eigenbrot Nov. 2014
#
###################################################
"""
**************
GradPak_skysub
**************

This is a script used to do sky subtraction on multispec files produced by the WIYN Bench Spectrograph and the GradPak IFU. To make use of this script the spectra need to be reduced usig the gradpak_sizes.iraf aperture table. This table defines different beam numbers for each fiber size such that the beam number is fibersize(in microns)/100. The sky fibers are assigned beam numbers as fibersize/100\*11.

General Usage
-------------
First, make sure you have *skysub = no* in your **dohydra** parameter file so that **dohydra** doesn't do (incorrect) sky subtraction first. After that using this module is very straight forward. Simply run::

    >python GradPak_skysub.py INPUT_FILE.ms_lin.fits OUTPUT_FILE.ms_s_lin.fits

That's it! If you follow the above naming convention and have your extracted, linearly sampled spectra in \*.ms_lin.fits files then you can omit the output file argument and GradPak_skysub will assume it to be of the form above (.ms_s_lin.fits).

When you run GradPak_skysub you will be shown the interactive IRAF **skysub** graphs 5 times; one for each fiber size. In these windows you can delete specific fibers with the 'd' key (use 'r' to redraw the plot after you do this). Once you are satisified with the selection of sky fibers hit 'q' to be taken to the next fiber size. Once all 5 sizes are done the output file will be generated and you'll be thrown back to your command prompt.

Individual functions are documented below, but only for completeness; you should never have to call any of them directly from within python.

Functions
---------
"""
import sys
import os

#Load the IRAF packages we'll need
try:
    current_dir = os.getcwd()
    if os.getlogin() == 'Arthur':
            os.chdir('/Users/Arthur/Ureka/iraf/local')
    from pyraf import iraf
    os.chdir(current_dir)
    iraf.imred(_doprint=0)
    iraf.hydra(_doprint=0)
except Exception as e:
    print "Failure: could not find pyraf/iraf"
    sys.exit(1)

def skysub(imagename,fibersize):
    """Take a multispec file with GradPak beam numbers and do a sky subtraction
    of a single fiber size. 

    Actual sky subtraction is done with IRAF's hydra.skysub. The fiber size is
    an integer that is the actual size of the fibers (in microns) divided by
    100. So the 200 micron fibers have fibersize=2.

    Parameters
    ----------

    imagename : str
        Name of the fits file to do sky subtraction on. It should be a file that can be recognized by the hydra package.

    fibersize : {2,3,4,5,6} int
        Fiber size to perform sky subtraction on. It is given in units of arcseconds rounded to the nearest integer.

    Returns
    -------
    skyname : str
        Name of the temporary sky-subtracted file generated by skysub. This file contains *only* fibers of the size specified.

    Notes
    -----
    This is a helper function that should probably not be called directly unless you know what you're doing. The main result is a .ms file that contains sky subtracted spectra for only that particular fiber size. For this function to work you need to have set *apidtab = gradpak_sizes.iraf* in **dohydra**.

    """
    skyname = '{:}.ms_s{:n}_lin.fits'.\
        format(imagename.split('.ms_lin.fits')[0],fibersize)
    print "Generating {} micron sky sub image {}".format(fibersize*100,skyname)

    iraf.skysub(imagename,
                output=skyname,
                objbeam='{},{}'.format(fibersize,fibersize*11),
                skybeam='{}'.format(fibersize*11),
                skyedit='yes',
                combine='average',
                reject='avsigclip',
                scale='none',
                savesky='yes',
                logfile='spool.txt')

    return skyname

def recombine(imagelist,outputimage):
    """Recombine a group of multispec files by aperture
    
    This uses the scombine task with group set to apertures. Because the
    aperture information is preserved when doing sky subtraction all this does
    is concatenate the individual .ms files into a single file.
    
    Parameters
    ----------
    imagelist : list of str
        List of individual file names of the temporary, single-size .ms files created by skysub

    outputimage : str
        Name of the resulting image that has all apertures from the input list

    Returns
    -------
    None
        The result is a multispec file that is a concatination of all the files in the input list.

    Notes
    -----
    To work as intended, each file in the input list should contain a unique set of aperture identifications. Any repeating apertures will be averaged together, which probably isn't what you want.

    """
    print 'Combining appertures from:'
    for image in imagelist: print '\t'+image

    iraf.scombine(','.join(imagelist),outputimage,
                  apertures='',
                  group='apertures',
                  combine='average', #These are irrelivant b/c
                  reject='avsigclip',# we aren't actually combining anything
                  first='yes', #Very important
                  w1='INDEF',w2='INDEF',dw='INDEF',nw='INDEF',
                  weight='',
                  log='no',
                  gain=0.438,
                  rdnoise=3.9,
                  logfile='spool.txt')

    return

def cleanup(imagelist):
    """Delete all the intermediate, single-fiber-size, sky-subtracted images

    Parameters
    ----------
    imagelist : list of str
        List of files to be deleted

    Returns
    -------
    None
        All files in the input list will be deleted

    """
    print 'cleaning intermediates:'
    for image in imagelist:
        print '\t'+image
        os.remove(image)

    return
                  
def main():
    """Do everything

    Take the input file and do a sky subtraction for each of the GradPak fiber
    sizes. Then combine everything back into the output file name.
    """

    if len(sys.argv) < 2 or len(sys.argv) > 3:
        #We don't want to overwhelm the user with verbosity
        return "The request was made, but it was not good"

    imagename = sys.argv[1]
    
    try:
        outputname = sys.argv[2]
    except IndexError:
        outputname = '{}.ms_s_lin.fits'.format(imagename.split('.ms_lin')[0])
        print "Using default output "+outputname

    #We do this check now because scombine doesn't have a clobber option
    if os.path.exists(outputname):
        print "Warning, you are about to overwrite an existing {}".\
            format(outputname)
        clobber = raw_input('Continue? (Y/n): ')
        if clobber.lower() == 'n':
            sys.exit()
        else:
            os.remove(outputname)
    
    #Fuck this line is cool
    skylist = [skysub(imagename,f) for f in range(2,7)]
    recombine(skylist,outputname)
    cleanup(skylist)

    return 0

if __name__ == '__main__':
    sys.exit(main())