Redoing the r13221-13223 merge in the actual trunk now. This updates trunk to the setup_refactoring branch (however, coords, pysynphot, and pywcs are still being pulled from the astrolib setup_refactoring branch. Will have to do that separately then update the svn:externals)

git-svn-id: http://svn.stsci.edu/svn/ssb/stsci_python/stsci_python/trunk/stwcs@13225 fe389314-cf27-0410-b35b-8c050e845b92

10 files changed, 3006 insertions, 0 deletions

diff --git a/lib/stwcs/wcsutil/__init__.py b/lib/stwcs/wcsutil/__init__.py
new file mode 100644
index 0000000..9b7ed8c
--- /dev/null
+++ b/lib/stwcs/wcsutil/__init__.py
@@ -0,0 +1,35 @@
+from __future__ import division # confidence high
+
+from altwcs import *
+from hstwcs import HSTWCS
+
+__docformat__ = 'restructuredtext'
+
+
+def help():
+    print 'How to create an HSTWCS object:\n\n'
+    print """ \
+    1. Using a pyfits HDUList object and an extension number \n
+    Example:\n
+
+    fobj = pyfits.open('some_file.fits')\n
+    w = wcsutil.HSTWCS(fobj, 3)\n\n
+
+    2. Create an HSTWCS object using a qualified file name. \n
+    Example:\n
+    w = wcsutil.HSTWCS('j9irw4b1q_flt.fits[sci,1]')\n\n
+
+    3. Create an HSTWCS object using a file name and an extension number. \n
+    Example:\n
+    w = wcsutil.HSTWCS('j9irw4b1q_flt.fits', ext=2)\n\n
+
+    4. Create an HSTWCS object from WCS with key 'O'.\n
+    Example:\n
+
+    w = wcsutil.HSTWCS('j9irw4b1q_flt.fits', ext=2, wcskey='O')\n\n
+
+    5. Create a template HSTWCS object for a DEFAULT object.\n
+    Example:\n
+    w = wcsutil.HSTWCS(instrument='DEFAULT')\n\n
+    """
+
diff --git a/lib/stwcs/wcsutil/altwcs.py b/lib/stwcs/wcsutil/altwcs.py
new file mode 100644
index 0000000..d250b52
--- /dev/null
+++ b/lib/stwcs/wcsutil/altwcs.py
@@ -0,0 +1,451 @@
+from __future__ import division # confidence high
+import os
+import string
+
+import numpy as np
+import pywcs
+import pyfits
+
+altwcskw = ['WCSAXES', 'CRVAL', 'CRPIX', 'PC', 'CDELT', 'CD', 'CTYPE', 'CUNIT',
+            'PV', 'PS']
+
+# file operations
+def archiveWCS(fname, ext, wcskey=" ", wcsname=" ", clobber=False):
+    """
+    Copy the primary WCS to the hader as an alternate WCS
+    with wcskey and name WCSNAME. It loops over all extensions in 'ext'
+
+    Parameters
+    ----------
+    fname:  string or pyfits.HDUList
+            a file name or a file object
+    ext:    an int, a tuple, a python list of integers or a python list of tuples (e.g.('sci',1))
+            fits extensions to work with
+    wcskey: string "A"-"Z" or " "
+            if " ": get next available key if wcsname is also " " or try
+            to get a key from WCSNAME value
+    wcsname: string
+             Name of alternate WCS description
+    clobber: boolean
+             if Ture - overwrites a WCS with the same key
+
+    See Also
+    --------
+    wcsutils.restoreWCS: Copy an alternate WCS to the primary WCS
+
+    """
+
+    if isinstance(fname, str):
+        f = pyfits.open(fname, mode='update')
+    else:
+        f = fname
+
+    if not parpasscheck(f, ext, wcskey, wcsname):
+        closefobj(fname, f)
+        return
+
+    if isinstance(ext, int) or isinstance(ext, tuple):
+        ext = [ext]
+    
+    wcsext = 0
+    eindx = 0
+    for e in f:
+        if 'extname' in e.header and 'sci' in e.header['extname'].lower():
+            wcsext = eindx
+            break
+        eindx += 1
+
+    if wcskey == " ":
+        # try getting the key from WCSNAME
+        if not wcsname.strip():
+            wkey = next_wcskey(f[wcsext].header)
+        else:
+            wkey = getKeyFromName(f[wcsext].header, wcsname)
+    else:
+        if wcskey not in available_wcskeys(f[wcsext].header):
+            # reuse wcsname
+            if not wcsname.strip():
+                wcsname = f[wcsext].header["WCSNAME"+wcskey]
+                wname = wcsname
+                wkey = wcskey
+            else:
+                wkey = wcskey
+                wname = wcsname
+        else:
+            wkey = wcskey
+            wname = wcsname
+
+    for e in ext:
+        w = pywcs.WCS(f[e].header, fobj=f)
+        hwcs = w.to_header()
+        wcsnamekey = 'WCSNAME' + wkey
+        f[e].header.update(key=wcsnamekey, value=wcsname)
+        if w.wcs.has_cd():
+            pc2cd(hwcs)
+        for k in hwcs.keys():
+            key = k[:7] + wkey
+            f[e].header.update(key=key, value=hwcs[k])
+        #norient = np.rad2deg(np.arctan2(hwcs['CD1_2'],hwcs['CD2_2']))
+        #okey = 'ORIENT%s' % wkey
+        #f[e].header.update(key=okey, value=norient)
+    closefobj(fname, f)
+
+def restoreWCS(f, ext, wcskey=" ", wcsname=" ", clobber=False):
+    """
+    Copy a WCS with key "WCSKEY" to a primary WCS
+
+    Reads in a WCS defined with wcskey and saves it as the primary WCS.
+    If clobber is False, writes out new files whose names are the original
+    names with an attached 3 character string  _'WCSKEY'_.
+    Otherwise overwrites the files. Goes sequentially through the list of extensions
+    The WCS is restored from the 'SCI' extension but the primary WCS of all
+    extensions with the same EXTVER are updated.
+
+
+    Parameters
+    ----------
+    f:       string or pyfits.HDUList object
+             a file name or a file object
+    ext:     an int, a tuple, a python list of integers or a python list
+             of tuples (e.g.('sci',1))
+             fits extensions to work with
+    wcskey:  a charater
+             "A"-"Z" - Used for one of 26 alternate WCS definitions.
+             or " " - find a key from WCSNAMe value
+    wcsname: string (optional)
+             if given and wcskey is " ", will try to restore by WCSNAME value
+    clobber: boolean
+             A flag to define if the original files should be overwritten
+
+    See Also
+    --------
+    wcsutil.archiveWCS - copy the primary WCS as an alternate WCS
+
+    """
+    if isinstance(f, str):
+        if clobber:
+            fobj = pyfits.open(f, mode='update')
+        else:
+            fobj = pyfits.open(f)
+    else:
+        fobj = f
+
+    if not parpasscheck(fobj, ext, wcskey, wcsname):
+        closefobj(f, fobj)
+        return
+
+    if isinstance(ext, int) or isinstance(ext, tuple):
+        ext = [ext]
+
+    if not clobber:
+        name = (fobj.filename().split('.fits')[0] + '_%s_' + '.fits') %wcskey
+    else:
+        name = fobj.filename()
+
+    if wcskey == " ":
+        if wcsname.strip():
+            wkey = getKeyFromName(fobj[1].header, wcsname)
+            if not wkey:
+                print 'Could not get a key from wcsname %s .' % wcsname
+                closefobj(f, fobj)
+                return
+    else:
+        if wcskey not in wcskeys(fobj[1].header):
+            print "Could not find alternate WCS with key %s in this file" % wcskey
+            closefobj(f, fobj)
+            return
+        wkey = wcskey
+
+    for e in ext:
+        try:
+            extname = fobj[e].header['EXTNAME'].lower()
+        except KeyError:
+            continue
+        #Restore always from a 'SCI' extension but write it out to 'ERR' and 'DQ'
+        if extname == 'sci':
+            sciver = fobj[e].header['extver']
+            try:
+                nwcs = pywcs.WCS(fobj[e].header, fobj=fobj, key=wkey)
+            except KeyError:
+                print 'restoreWCS: Could not read WCS with key %s in file %s,  \
+                extension %d' % (wkey, fobj.filename(), e)
+                closefobj(f, fobj)
+                return #raise
+            hwcs = nwcs.to_header()
+
+            if nwcs.wcs.has_cd():
+                pc2cd(hwcs, key=wkey)
+            for k in hwcs.keys():
+                key = k[:-1]
+                if key in fobj[e].header.keys():
+                    fobj[e].header.update(key=key, value = hwcs[k])
+                else:
+                    continue
+            if wcskey == 'O' and fobj[e].header.has_key('TDDALPHA'):
+                fobj[e].header['TDDALPHA'] = 0.0
+                fobj[e].header['TDDBETA'] = 0.0
+            if fobj[e].header.has_key('ORIENTAT'):
+                norient = np.rad2deg(np.arctan2(hwcs['CD1_2'+'%s' %wkey],hwcs['CD2_2'+'%s' %wkey]))
+                fobj[e].header.update(key='ORIENTAT', value=norient)
+        elif extname in ['err', 'dq', 'sdq', 'time', 'samp']:
+            cextver = fobj[e].header['extver']
+            if cextver == sciver:
+                for k in hwcs.keys():
+                    key = k[:-1]
+                    fobj[e].header.update(key=key, value = hwcs[k])
+                if fobj[e].header.has_key('ORIENTAT'):
+                    norient = np.rad2deg(np.arctan2(hwcs['CD1_2'+'%s' %wkey],hwcs['CD2_2'+'%s' %wkey]))
+                    fobj[e].header.update(key='ORIENTAT', value=norient)
+        else:
+            continue
+
+    if not clobber:
+        fobj.writeto(name)
+    closefobj(f, fobj)
+
+def deleteWCS(fname, ext, wcskey=" ", wcsname=" "):
+    """
+    Delete an alternate WCS defined with wcskey.
+    If wcskey is " " try to get a key from WCSNAME.
+
+    Parameters
+    ----------
+    fname:   sting or a pyfits.HDUList object
+    ext:     an int, a tuple, a python list of integers or a python list of tuples (e.g.('sci',1))
+             fits extensions to work with
+    wcskey:  one of 'A'-'Z' or " "
+    wcsname: string
+             Name of alternate WCS description
+    """
+    if isinstance(fname, str):
+        fobj = pyfits.open(fname, mode='update')
+    else:
+        fobj = fname
+
+    if not parpasscheck(fobj, ext, wcskey, wcsname):
+        closefobj(fname, fobj)
+        return
+
+    if isinstance(ext, int) or isinstance(ext, tuple):
+        ext = [ext]
+
+    # Do not allow deleting the original WCS.
+    if wcskey == 'O':
+        print "Wcskey 'O' is reserved for the original WCS and should not be deleted."
+        closefobj(fname, fobj)
+        return
+
+    if wcskey == " ":
+        # try getting the key from WCSNAME
+        if wcsname == " ":
+            print "Could not get a valid key from header"
+            closefobj(fname, fobj)
+            return
+        else:
+            wkey = getKeyFromName(fobj[1].header, wcsname)
+            if not wkey:
+                print 'Could not get a key: wcsname "%s" not found in header.' % wcsname
+                closefobj(fname, fobj)
+                return
+    else:
+        if wcskey not in wcskeys(fobj[1].header):
+            print "Could not find alternate WCS with key %s in this file" % wcskey
+            closefobj(fname, fobj)
+            return
+        wkey = wcskey
+
+    prexts = []
+    for i in ext:
+        hdr = fobj[i].header
+        try:
+            w = pywcs.WCS(hdr, fobj, key=wkey)
+        except KeyError:
+            continue
+        hwcs = w.to_header()
+        if w.wcs.has_cd():
+            pc2cd(hwcs, key=wkey)
+        for k in hwcs.keys():
+            del hdr[k]
+            #del hdr['ORIENT'+wkey]
+        prexts.append(i)
+    if prexts != []:
+        print 'Deleted all instances of WCS with key %s in extensions' % wkey, prexts
+    else:
+        print "Did not find WCS with key %s in any of the extensions" % wkey
+    closefobj(fname, fobj)
+
+#header operations
+def wcskeys(header):
+    """
+    Returns a list of characters used in the header for alternate
+    WCS description with WCSNAME keyword
+
+    Parameters
+    ----------
+    hdr: pyfits.Header
+    """
+    assert isinstance(header, pyfits.Header), "Requires a pyfits.Header object as input"
+    names = header["WCSNAME*"]
+    return [key.split('WCSNAME')[1].upper() for key in names.keys()]
+
+def wcsnames(header):
+    """
+    Returns a dictionary of wcskey: WCSNAME pairs
+
+    Parameters
+    ----------
+    header: pyfits.Header
+    """
+    names = header["WCSNAME*"]
+    d = {}
+    for c in names:
+        d[c.key[-1]] = c.value
+    return d
+
+def available_wcskeys(header):
+    """
+    Returns a list of characters which are not used in the header
+    with WCSNAME keyword. Any of them can be used to save a new
+    WCS.
+
+    Parameters
+    ----------
+    header: pyfits.Header
+    """
+    assert isinstance(header, pyfits.Header), "Requires a pyfits.Header object as input"
+    all_keys = list(string.ascii_uppercase)
+    used_keys = wcskeys(header)
+    try:
+        used_keys.remove("")
+    except ValueError:
+        pass
+    [all_keys.remove(key) for key in used_keys]
+    return all_keys
+
+def next_wcskey(header):
+    """
+    Returns next available character to be used for an alternate WCS
+
+    Parameters
+    ----------
+    header: pyfits.Header
+    """
+    assert isinstance(header, pyfits.Header), "Requires a pyfits.Header object as input"
+    allkeys = available_wcskeys(header)
+    if allkeys != []:
+        return allkeys[0]
+    else:
+        return None
+
+def getKeyFromName(header, wcsname):
+    """
+    If WCSNAME is found in header, return its key, else return
+    None. This is used to update an alternate WCS
+    repeatedly and not generate new keys every time.
+
+    Parameters
+    ----------
+    header:  pyfits.Header
+    wcsname: str
+             Value of WCSNAME
+    """
+    wkey = None
+    names = wcsnames(header)
+    for item in names.items():
+        if item[1].lower() == wcsname.lower():
+            wkey = item[0]
+            break
+    return wkey
+
+def pc2cd(hdr, key=' '):
+    """
+    Convert a CD PC matrix to a CD matrix.
+
+    WCSLIB (and PyWCS) recognizes CD keywords as input
+    but converts them and works internally with the PC matrix.
+    to_header() returns the PC matrix even if the i nput was a
+    CD matrix. To keep input and output consistent we check
+    for has_cd and convert the PC back to CD.
+
+    Parameters
+    ----------
+    hdr: pyfits.Header
+
+    """
+    for c in ['1_1', '1_2', '2_1', '2_2']:
+        try:
+            val = hdr['PC'+c+'%s' % key]
+            del hdr['PC'+c+ '%s' % key]
+        except KeyError:
+            if c=='1_1' or c == '2_2':
+                val = 1.
+            else:
+                val = 0.
+        hdr.update(key='CD'+c+'%s' %key, value=val)
+    return hdr
+
+def parpasscheck(fobj, ext, wcskey, wcsname, clobber=True):
+
+    if not isinstance(fobj,pyfits.HDUList):
+        print "First parameter must be a fits file object or a file name."
+        return False
+    try:
+        assert (fobj.fileinfo(0)['filemode'] == 'update')
+    except AssertionError:
+        print "First parameter must be a file name or a file object opened in 'update' mode."
+        return False
+
+    if not isinstance(ext, int) and not isinstance(ext, tuple) \
+       and not isinstance(ext, list):
+        print "Ext must be integer, tuple, a list of int extension numbers, \
+        or a list of tuples representing a fits extension, for example ('sci', 1)."
+        return False
+
+    if len(wcskey) != 1:
+        print 'Parameter wcskey must be a character - one of "A"-"Z" or " "'
+        return False
+
+    wcsext = 0
+    eindx = 0
+    for e in fobj:
+        if 'extname' in e.header and 'sci' in e.header['extname'].lower():
+            wcsext = eindx
+            break
+        eindx += 1
+        
+
+    if wcskey == " ":
+        # try getting the key from WCSNAME
+        """
+        if wcsname == " " or wcsname == "":
+            #wkey = next_wcskey(f[1].header)
+            #if not wkey:
+            #    print "Could not get a valid key from header"
+            return False
+        """
+        if wcsname.strip():
+            wkey = getKeyFromName(fobj[wcsext].header, wcsname)
+            if wkey and not clobber:
+                print 'Wcsname %s is already used.' % wcsname
+                print 'Use "wcsutil.next_wcskey" to obtain a valid wcskey'
+                print 'or use "clobber=True" to overwrite the values.'
+                return False
+    else:
+        if wcskey not in available_wcskeys(fobj[wcsext].header):
+            if clobber==False:
+                print 'Wcskey %s is already used.' % wcskey
+                print 'Use "wcsutil.next_wcskey" to obtain a valid wcskey'
+                print 'or use "clobber=True" to overwrite the values.'
+                return False
+
+
+    return True
+
+def closefobj(fname, f):
+    """
+    Functions in this module accept as input a file name or a file object.
+    If the input was a file name (string) we close the object. If the user
+    passed a file object we leave it to the user to close it.
+    """
+    if isinstance(fname, str):
+        f.close()
diff --git a/lib/stwcs/wcsutil/convertwcs.py b/lib/stwcs/wcsutil/convertwcs.py
new file mode 100644
index 0000000..80276c3
--- /dev/null
+++ b/lib/stwcs/wcsutil/convertwcs.py
@@ -0,0 +1,118 @@
+import pyfits
+try:
+    import stwcs
+    from stwcs import wcsutil
+except:
+    stwcs = None
+
+from stsci.tools import fileutil
+
+OPUS_WCSKEYS = ['OCRVAL1','OCRVAL2','OCRPIX1','OCRPIX2',
+                'OCD1_1','OCD1_2','OCD2_1','OCD2_2',
+                'OCTYPE1','OCTYPE2']
+
+
+def archive_prefix_OPUS_WCS(fobj,extname='SCI'):
+    """ Identifies WCS keywords which were generated by OPUS and archived
+        using a prefix of 'O' for all 'SCI' extensions in the file
+
+        Parameters
+        ----------
+        fobj: string or pyfits.HDUList
+            Filename or pyfits object of a file
+
+    """
+    if stwcs is None:
+        print '====================='
+        print 'The STWCS package is needed to convert an old-style OPUS WCS to an alternate WCS'
+        print '====================='
+        raise ImportError
+
+
+    closefits = False
+    if isinstance(fobj,str):
+        # A filename was provided as input
+        fobj = pyfits.open(fobj,mode='update')
+        closefits=True
+
+    # Define the header
+    ext = ('sci',1)
+    hdr = fobj[ext].header
+
+    numextn = fileutil.countExtn(fobj)
+    extlist = []
+    for e in xrange(1,numextn+1):
+        extlist.append(('sci',e))
+
+    # Insure that the 'O' alternate WCS is present
+    if 'O' not in wcsutil.wcskeys(hdr):
+        # if not, archive the Primary WCS as the default OPUS WCS
+        wcsutil.archiveWCS(fobj,extlist, wcskey='O', wcsname='OPUS')
+
+    # find out how many SCI extensions are in the image
+    numextn = fileutil.countExtn(fobj,extname=extname)
+    if numextn == 0:
+        extname = 'PRIMARY'
+
+    # create HSTWCS object from PRIMARY WCS
+    wcsobj = wcsutil.HSTWCS(fobj,ext=ext,wcskey='O')
+    # get list of WCS keywords
+    wcskeys = wcsobj.wcs2header().keys()
+
+    # For each SCI extension...
+    for e in xrange(1,numextn+1):
+        # Now, look for any WCS keywords with a prefix of 'O'
+        for key in wcskeys:
+            okey = 'O'+key[:7]
+            hdr = fobj[(extname,e)].header
+            if hdr.has_key(okey):
+                # Update alternate WCS keyword with prefix-O OPUS keyword value
+                hdr[key] = hdr[okey]
+
+    if closefits:
+        fobj.close()
+
+def create_prefix_OPUS_WCS(fobj,extname='SCI'):
+    """ Creates alternate WCS with a prefix of 'O' for OPUS generated WCS values
+        to work with old MultiDrizzle.
+
+        Parameters
+        ----------
+        fobj: string or pyfits.HDUList
+            Filename or pyfits object of a file
+
+        Raises
+        ------
+        IOError:
+            if input PyFITS object was not opened in 'update' mode
+
+    """
+    # List of O-prefix keywords to create
+    owcskeys = OPUS_WCSKEYS
+
+    closefits = False
+    if isinstance(fobj,str):
+        # A filename was provided as input
+        fobj = pyfits.open(fobj,mode='update')
+        closefits=True
+    else:
+        # check to make sure this FITS obj has been opened in update mode
+        if fobj.fileinfo(0)['filemode'] != 'update':
+            print 'File not opened with "mode=update". Quitting...'
+            raise IOError
+
+    # check for existance of O-prefix WCS
+    if not fobj['sci',1].header.has_key(owcskeys[0]):
+
+        # find out how many SCI extensions are in the image
+        numextn = fileutil.countExtn(fobj,extname=extname)
+        if numextn == 0:
+            extname = ''
+        for extn in xrange(1,numextn+1):
+            hdr = fobj[(extname,extn)].header
+            for okey in owcskeys:
+                hdr.update(okey,hdr[okey[1:]+'O'])
+
+    # Close FITS image if we had to open it...
+    if closefits:
+        fobj.close()
diff --git a/lib/stwcs/wcsutil/getinput.py b/lib/stwcs/wcsutil/getinput.py
new file mode 100644
index 0000000..2f64f46
--- /dev/null
+++ b/lib/stwcs/wcsutil/getinput.py
@@ -0,0 +1,62 @@
+import pyfits
+from stsci.tools import irafglob, fileutil, parseinput
+
+def parseSingleInput(f=None, ext=None):
+    if isinstance(f, str):
+        # create an HSTWCS object from a filename
+        if ext != None:
+            filename = f
+            if isinstance(ext,tuple):
+                if ext[0] == '':
+                    extnum = ext[1] # handle ext=('',1)
+                else:
+                    extnum = ext
+            else:
+                extnum = int(ext)
+        elif ext == None:
+            filename, ext = fileutil.parseFilename(f)
+            ext = fileutil.parseExtn(ext)
+            if ext[0] == '':
+                extnum = int(ext[1]) #handle ext=('',extnum)
+            else:
+                extnum = ext
+        phdu = pyfits.open(filename)
+        hdr0 = pyfits.getheader(filename)
+        try:
+            ehdr = pyfits.getheader(filename, ext=extnum)
+        except (IndexError,KeyError):
+            print 'Unable to get extension.', extnum
+            raise
+
+    elif isinstance(f, pyfits.HDUList):
+        phdu = f
+        if ext == None:
+            extnum = 0
+        else:
+            extnum = ext
+        ehdr = f[extnum].header
+        hdr0 = f[0].header
+        filename = hdr0.get('FILENAME', "")
+        
+    else:
+        raise ValueError('Input must be a file name string or a pyfits file object')
+
+    return filename, hdr0, ehdr, phdu
+
+
+def parseMultipleInput(input):
+    if isinstance(input, str):
+        if input[0] == '@':
+            # input is an @ file
+            filelist = irafglob.irafglob(input)
+        else:
+            try:
+                filelist, output = parseinput.parseinput(input)
+            except IOError: raise
+    elif isinstance(input, list):
+        if isinstance(input[0], wcsutil.HSTWCS):
+            # a list of HSTWCS objects
+            return input
+        else:
+            filelist = input[:]
+    return filelist
\ No newline at end of file
diff --git a/lib/stwcs/wcsutil/headerlet.py b/lib/stwcs/wcsutil/headerlet.py
new file mode 100644
index 0000000..0318bf2
--- /dev/null
+++ b/lib/stwcs/wcsutil/headerlet.py
@@ -0,0 +1,898 @@
+from __future__ import division
+import logging
+import os
+import tarfile
+import tempfile
+import time
+import warnings
+from cStringIO import StringIO
+
+import numpy as np
+import pyfits
+
+import altwcs
+import wcscorr
+from hstwcs import HSTWCS
+from mappings import basic_wcs
+from stsci.tools.fileutil import countExtn
+
+module_logger = logging.getLogger('headerlet')
+
+import atexit
+atexit.register(logging.shutdown)
+
+def setLogger(logger, level, mode='w'):
+    formatter = logging.Formatter("%(asctime)s - %(name)s - %(levelname)s - %(message)s")
+    log_filename = 'headerlet.log'
+    fh = logging.FileHandler(log_filename, mode=mode)
+    fh.setLevel(logging.DEBUG)
+    fh.setFormatter(formatter)
+    logger.addHandler(fh)
+    logger.setLevel(level)
+    
+def isWCSIdentical(scifile, file2, verbose=False):
+    """
+    Compares the WCS solution of 2 files.
+
+    Parameters
+    ----------
+    scifile: file1
+    file2: file2
+    verbose: False or a python logging level
+             (one of 'INFO', 'DEBUG' logging levels)
+             (an integer representing a logging level)
+
+    Notes
+    -----
+    These can be 2 science observations or 2 headerlets
+    or a science observation and a headerlet. The two files
+    have the same WCS solution if the following are the same:
+
+    - rootname/destim
+    - primary WCS
+    - SIP coefficients
+    - NPOL distortion
+    - D2IM correction
+    - Velocity aberation
+
+    """
+    if verbose:
+        setLogger(module_logger, verbose)
+    else:
+        module_logger.setLevel(100)
+
+    module_logger.info("Starting isWCSIdentical: %s" % time.asctime())
+
+    result = True
+    numsci1 = max(countExtn(scifile), countExtn(scifile, 'SIPWCS'))
+    numsci2 = max(countExtn(file2), countExtn(file2, 'SIPWCS'))
+
+    if numsci1 == 0 or numsci2 == 0 or numsci1 != numsci2:
+        module_logger.info("Number of SCI and SIPWCS extensions do not match.")
+        result = False
+
+    if getRootname(scifile) != getRootname(file2):
+        module_logger.info('Rootnames do not match.')
+        result = False
+    try:
+        extname1 = pyfits.getval(scifile, 'EXTNAME', ext=('SCI', 1))
+    except KeyError:
+        extname1 = 'SIPWCS'
+    try:
+        extname2 = pyfits.getval(file2, 'EXTNAME', ext=('SCI', 1))
+    except KeyError:
+        extname2 = 'SIPWCS'
+
+    for i in range(1, numsci1 + 1):
+        w1 = HSTWCS(scifile, ext=(extname1, i))
+        w2 = HSTWCS(file2, ext=(extname2, i))
+        if not np.allclose(w1.wcs.crval, w2.wcs.crval, rtol=1e-7) or \
+        not np.allclose(w1.wcs.crpix, w2.wcs.crpix, rtol=1e-7)  or \
+        not np.allclose(w1.wcs.cd, w2.wcs.cd, rtol=1e-7) or \
+        not (np.array(w1.wcs.ctype) == np.array(w2.wcs.ctype)).all():
+            module_logger.info('Primary WCSs do not match')
+            result = False
+        if w1.sip or w2.sip:
+            if (w2.sip and not w1.sip) or (w1.sip and not w2.sip) or \
+               not np.allclose(w1.sip.a, w2.sip.a, rtol=1e-7) or \
+               not np.allclose(w1.sip.b, w2.sip.b, rtol=1e-7):
+                module_logger.info('SIP coefficients do not match')
+                result = False
+        if w1.cpdis1 or w2.cpdis1:
+            if w1.cpdis1 and not w2.cpdis1 or \
+                w2.cpdis1 and not w1.cpdis1 or \
+                not np.allclose(w1.cpdis1.data, w2.cpdis1.data):
+                module_logger.info('NPOL distortions do not match')
+                result = False
+        if w1.cpdis2 or w2.cpdis2:
+            if w1.cpdis2 and not w2.cpdis2 or \
+                w2.cpdis2 and not w1.cpdis2 or \
+                not np.allclose(w1.cpdis2.data, w2.cpdis2.data):
+                module_logger.info('NPOL distortions do not match')
+                result = False
+        if w1.det2im1 or w2.det2im1:
+            if w1.det2im1 and not w2.det2im1 or \
+                w2.det2im1 and not w1.det2im1 or\
+                not np.allclose(w1.det2im1.data, w2.det2im1.data):
+                module_logger.info('Det2Im corrections do not match')
+                result =  False
+        if w1.det2im2 or w2.det2im2:
+            if w1.det2im2 and not w2.det2im2 or \
+                w2.det2im2 and not w1.det2im2 or\
+                not np.allclose(w1.det2im2.data, w2.det2im2.data):
+                module_logger.info('Det2Im corrections do not match')
+                result = False
+        if w1.vafactor != w2.vafactor:
+            module_logger.info('VA factors do not match')
+            result = False
+
+    return result
+
+
+# TODO: It would be logical for this to be part of the Headerlet class, perhaps
+# as a classmethod
+def createHeaderlet(fname, hdrname, destim=None, output=None, verbose=False, logmode='w'):
+    """
+    Create a headerlet from a science observation
+
+    Parameters
+    ----------
+    fname: string
+           Name of file with science observation
+    hdrname: string
+           Name for the headerlet, stored in the primary header of the headerlet
+    destim: string
+           Destination image, stored in the primary header of the headerlet.
+           If None ROOTNAME is used of the science observation is used.
+           ROOTNAME has precedence, destim is used for observations without
+           ROOTNAME in the primary header
+    output: string
+           Save the headerlet to the given filename.
+    verbose: False or a python logging level
+             (one of 'INFO', 'DEBUG' logging levels)
+             (an integer representing a logging level)
+    logmode: 'w', 'a'
+            used internally for controling access to the log file
+    """
+
+    if verbose:
+        setLogger(module_logger, verbose, mode=logmode)
+    else:
+        module_logger.setLevel(100)
+
+    module_logger.info("Starting createHeaderlet: %s" % time.asctime())
+    phdukw = {'IDCTAB': True,
+            'NPOLFILE': True,
+            'D2IMFILE': True}
+    if not isinstance(fname, pyfits.HDUList):
+        fobj = pyfits.open(fname)
+        close_file = True
+    else:
+        fobj = fname
+        close_file = False
+    if destim is None:
+        try:
+            destim = fobj[0].header['ROOTNAME']
+        except KeyError:
+            module_logger.exception('Required keyword "DESTIM" not found')
+            print 'Please provide a value for the DESTIM keyword'
+            raise
+    if hdrname is None:
+        module_logger.critical("Required keyword 'HDRNAME' not given")
+        raise ValueError("Please provide a name for the headerlet, HDRNAME is "
+                         "a required parameter.")
+        
+    # get the version of STWCS used to create the WCS of the science file.
+    try:
+        upwcsver = fobj[0].header.ascard['STWCSVER']
+    except KeyError:
+        upwcsver = pyfits.Card("STWCSVER", " ",
+                               "Version of STWCS used to update the WCS")
+    
+    # get all keys for alternate WCSs
+    altkeys = altwcs.wcskeys(fobj[('SCI', 1)].header)
+    
+    if 'O' in altkeys:
+        altkeys.remove('O')
+    
+    numsci = countExtn(fname, 'SCI')
+    module_logger.debug("Number of 'SCI' extensions in file %s is %s"
+                 % (fname, numsci))
+    hdul = pyfits.HDUList()
+    phdu = pyfits.PrimaryHDU()
+    phdu.header.update('DESTIM', destim,
+                       comment='Destination observation root name')
+    phdu.header.update('HDRNAME', hdrname, comment='Headerlet name')
+    fmt="%Y-%m-%dT%H:%M:%S"
+    phdu.header.update('DATE', time.strftime(fmt),
+                       comment='Date FITS file was generated')
+    phdu.header.ascard.append(upwcsver)
+
+    refs = updateRefFiles(fobj[0].header.ascard, phdu.header.ascard, verbose=verbose)
+    phdukw.update(refs)
+    phdu.header.update(key='VAFACTOR',
+                       value=fobj[('SCI',1)].header.get('VAFACTOR', 1.))
+    hdul.append(phdu)
+
+    for e in range(1, numsci + 1):
+        hwcs = HSTWCS(fname, ext=('SCI', e))
+        h = hwcs.wcs2header(sip2hdr=True).ascard
+        for ak in altkeys:
+            awcs = HSTWCS(fname,ext=('SCI', e), wcskey=ak)
+            h.extend(awcs.wcs2header(idc2hdr=False).ascard)
+        h.append(pyfits.Card(key='VAFACTOR', value=hwcs.vafactor, 
+                             comment='Velocity aberration plate scale factor'))
+        h.insert(0, pyfits.Card(key='EXTNAME', value='SIPWCS',
+                                comment='Extension name'))
+        h.insert(1, pyfits.Card(key='EXTVER', value=e,
+                                comment='Extension version'))
+
+        fhdr = fobj[('SCI', e)].header.ascard
+        if phdukw['NPOLFILE']:
+            cpdis = fhdr['CPDIS*...']
+            for c in range(1, len(cpdis) + 1):
+                h.append(cpdis[c - 1])
+                dp = fhdr['DP%s.*...' % c]
+                h.extend(dp)
+
+                try:
+                    h.append(fhdr['CPERROR%s' % c])
+                except KeyError:
+                    pass
+
+            try:
+                h.append(fhdr['NPOLEXT'])
+            except KeyError:
+                pass
+
+        if phdukw['D2IMFILE']:
+            try:
+                h.append(fhdr['D2IMEXT'])
+            except KeyError:
+                pass
+
+            try:
+                h.append(fhdr['AXISCORR'])
+            except KeyError:
+                module_logger.exception("'D2IMFILE' kw exists but keyword 'AXISCORR' was not found in "
+                                 "%s['SCI',%d]" % (fname, e))
+                raise
+
+            try:
+                h.append(fhdr['D2IMERR'])
+            except KeyError:
+                h.append(pyfits.Card(key='DPERROR', value=0,
+                                     comment='Maximum error of D2IMARR'))
+
+        hdu = pyfits.ImageHDU(header=pyfits.Header(h))
+        # temporary fix for pyfits ticket # 48
+        hdu._extver = e
+        hdul.append(hdu)
+    numwdvarr = countExtn(fname, 'WCSDVARR')
+    numd2im = countExtn(fname, 'D2IMARR')
+    for w in range(1, numwdvarr + 1):
+        hdu = fobj[('WCSDVARR', w)].copy()
+        # temporary fix for pyfits ticket # 48
+        hdu._extver = w
+        hdul.append(hdu)
+    for d in range(1, numd2im + 1):
+        hdu = fobj[('D2IMARR', d)].copy()
+        # temporary fix for pyfits ticket # 48
+        hdu._extver = d
+        hdul.append(hdu)
+    if output is not None:
+        # write the headerlet to a file
+        if not output.endswith('_hdr.fits'):
+            output = output + '_hdr.fits'
+        hdul.writeto(output, clobber=True)
+    if close_file:
+        fobj.close()
+    return Headerlet(hdul,verbose=verbose, logmode='a')
+
+def applyHeaderlet(hdrfile, destfile, createheaderlet=True, hdrname=None,
+                   verbose=False):
+    """
+    Apply headerlet 'hdrfile' to a science observation 'destfile'
+
+    Parameters
+    ----------
+    hdrfile: string
+             Headerlet file
+    destfile: string
+             File name of science observation whose WCS solution will be updated
+    createheaderlet: boolean
+            True (default): before updating, create a headerlet with the
+            WCS old solution.
+    hdrname: string or None (default)
+            will be the value of the HDRNAME keyword in the headerlet generated
+            for the old WCS solution.  If not specified, a sensible default
+            will be used.  Not required if createheaderlet is False
+    verbose: False or a python logging level
+             (one of 'INFO', 'DEBUG' logging levels)
+             (an integer representing a logging level)
+    """
+    if verbose:
+        setLogger(module_logger, verbose)
+    else:
+        module_logger.setLevel(100)
+    module_logger.info("Starting applyHeaderlet: %s" % time.asctime())
+    hlet = Headerlet(hdrfile, verbose=verbose, logmode='a')
+    hlet.apply(destfile, createheaderlet=createheaderlet, hdrname=hdrname)
+
+def updateRefFiles(source, dest, verbose=False):
+    """
+    Update the reference files name in the primary header of 'dest'
+    using values from 'source'
+
+    Parameters
+    ----------
+    source: pyfits.Header.ascardlist
+    dest:   pyfits.Header.ascardlist
+    """
+    module_logger.info("Updating reference files")
+    phdukw = {'IDCTAB': True,
+            'NPOLFILE': True,
+            'D2IMFILE': True}
+
+    try:
+        wind = dest.index_of('HISTORY')
+    except KeyError:
+        wind = len(dest)
+    for key in phdukw.keys():
+        try:
+            value = source[key]
+            dest.insert(wind, value)
+        except KeyError:
+            # TODO: I don't understand what the point of this is.  Is it meant
+            # for logging purposes?  Right now it isn't used.
+            phdukw[key] = False
+    return phdukw
+
+def getRootname(fname):
+    """
+    returns the value of ROOTNAME or DESTIM
+    """
+
+    try:
+        rootname = pyfits.getval(fname, 'ROOTNAME')
+    except KeyError:
+        rootname = pyfits.getval(fname, 'DESTIM')
+    return rootname
+
+def mapFitsExt2HDUListInd(fname, extname):
+    """
+    Map FITS extensions with 'EXTNAME' to HDUList indexes.
+    """
+
+    if not isinstance(fname, pyfits.HDUList):
+        f = pyfits.open(fname)
+        close_file = True
+    else:
+        f = fname
+        close_file = False
+    d = {}
+    for hdu in f:
+        # TODO: Replace calls to header.has_key() with `in header` once
+        # pyfits refactoring branch is in production use
+        if hdu.header.has_key('EXTNAME') and hdu.header['EXTNAME'] == extname:
+            extver = hdu.header['EXTVER']
+            d[(extname, extver)] = f.index_of((extname, extver))
+    if close_file:
+        f.close()
+    return d
+
+
+class Headerlet(pyfits.HDUList):
+    """
+    A Headerlet class
+    Ref: http://stsdas.stsci.edu/stsci_python_sphinxdocs/stwcs/headerlet_def.html
+    """
+
+    def __init__(self, fobj, wcskeys=[], mode='copyonwrite', verbose=False, logmode='w'):
+        """
+        Parameters
+        ----------
+        fobj:  string
+                Name of headerlet file, file-like object, a list of HDU
+                instances, or an HDUList instance
+        wcskeys: python list
+                a list of wcskeys to be included in the headerlet
+                created from the old WCS solution before the
+                science file is updated. If empty: all alternate (if any)
+                WCSs are copied to the headerlet.
+        mode: string, optional
+                Mode with which to open the given file object
+        """
+        self.verbose = verbose
+        self.hdr_logger = logging.getLogger('headerlet.Headerlet')
+        if self.verbose:
+                setLogger(self.hdr_logger, self.verbose, mode=logmode)
+        else:
+            self.hdr_logger.setLevel(100)
+
+        self.hdr_logger.info("Creating a Headerlet object from wcskeys %s" % wcskeys)
+        self.wcskeys = wcskeys
+        if not isinstance(fobj, list):
+            fobj = pyfits.open(fobj, mode=mode)
+
+        super(Headerlet, self).__init__(fobj)
+        self.fname = self.filename()
+        self.hdrname = self[0].header["HDRNAME"]
+        self.stwcsver = self[0].header.get("STWCSVER", "")
+        self.destim = self[0].header["DESTIM"]
+        self.idctab = self[0].header.get("IDCTAB", "")
+        self.npolfile = self[0].header.get("NPOLFILE", "")
+        self.d2imfile = self[0].header.get("D2IMFILE", "")
+        self.vafactor = self[1].header.get("VAFACTOR", 1) #None instead of 1?
+        self.d2imerr = 0
+        self.axiscorr = 1
+        
+    def apply(self, dest, createheaderlet=True, hdrname=None, attach=True, createsummary=True):
+        """
+        Apply this headerlet to a file.
+        
+        Parameters
+        ----------
+        dest:    string
+                 Name of file to which to apply the WCS in the headerlet
+        hdrname: string
+                 A unique name for the headerlet
+        createheaderlet: boolean
+                 A flag which indicates if a headerlet should be created 
+                 from the old WCS and attached to the science file (default: True)
+        attach:  boolean, default: True
+                 By default the headerlet being applied will be attached 
+                 as an extension to the science file. Set 'attach' to False
+                 to disable this.
+        createsummary: boolean, default: True
+                 Set this to False to disable creating and updating of wcscorr table.
+                 This is used primarily for testing.
+        """
+        self.hverify()
+        if self.verify_dest(dest):
+            if not isinstance(dest, pyfits.HDUList):
+                fobj = pyfits.open(dest, mode='update')
+                close_dest = True
+            else:
+                fobj = dest
+                close_dest = False
+
+            # Create the WCSCORR HDU/table from the existing WCS keywords if
+            # necessary
+            if createsummary:
+                try:
+                    # TODO: in the pyfits refactoring branch if will be easier to
+                    # test whether an HDUList contains a certain extension HDU
+                    # without relying on try/except
+                    wcscorr_table = fobj['WCSCORR']
+                except KeyError:
+                    # The WCSCORR table needs to be created
+                    wcscorr.init_wcscorr(fobj)
+
+            orig_hlt_hdu = None
+            numhlt = countExtn(fobj, 'HDRLET')
+            if createheaderlet:
+                # Create a headerlet for the original WCS data in the file,
+                # create an HDU from the original headerlet, and append it to
+                # the file
+                if not hdrname:
+                    hdrname = fobj[0].header['ROOTNAME'] + '_orig'
+                orig_hlt = createHeaderlet(fobj, hdrname, verbose=self.verbose, logmode='a')
+                orig_hlt_hdu = HeaderletHDU.fromheaderlet(orig_hlt)
+                orig_hlt_hdu.update_ext_version(numhlt + 1)
+                numhlt += 1
+
+            self._delDestWCS(fobj)
+            refs = updateRefFiles(self[0].header.ascard, fobj[0].header.ascard, verbose=self.verbose)
+            numsip = countExtn(self, 'SIPWCS')
+            for idx in range(1, numsip + 1):
+                fhdr = fobj[('SCI', idx)].header.ascard
+                siphdr = self[('SIPWCS', idx)].header.ascard
+                # a minimal attempt to get the position of the WCS keywords group
+                # in the header by looking for the PA_APER kw.
+                # at least make sure the WCS kw are written befir the HISTORY kw
+                # if everything fails, append the kw to the header
+                try:
+                    wind = fhdr.index_of('PA_APER')
+                except KeyError:
+                    try:
+                        wind = fhdr.index_of('HISTORY')
+                    except KeyError:
+                        wind = len(fhdr)
+                self.hdr_logger.debug("Inserting WCS keywords at index %s" % wind)
+                # TODO: Drop .keys() when refactored pyfits comes into use
+                for k in siphdr.keys():
+                    if k not in ['XTENSION', 'BITPIX', 'NAXIS', 'PCOUNT',
+                                 'GCOUNT','EXTNAME', 'EXTVER', 'ORIGIN',
+                                 'INHERIT', 'DATE', 'IRAF-TLM']:
+                        fhdr.insert(wind, siphdr[k])
+                    else:
+                        pass
+
+            #! Always attach these extensions last. Otherwise their headers may
+            # get updated with the other WCS kw.
+            numwdvar = countExtn(self, 'WCSDVARR')
+            numd2im = countExtn(self, 'D2IMARR')
+            for idx in range(1, numwdvar + 1):
+                fobj.append(self[('WCSDVARR', idx)].copy())
+            for idx in range(1, numd2im + 1):
+                fobj.append(self[('D2IMARR', idx)].copy())
+
+            # Update the WCSCORR table with new rows from the headerlet's WCSs
+            if createsummary:
+                wcscorr.update_wcscorr(fobj, self, 'SIPWCS')
+
+            # Append the original headerlet
+            if createheaderlet and orig_hlt_hdu:
+                fobj.append(orig_hlt_hdu)
+            
+            if attach:
+                # Finally, append an HDU for this headerlet
+                new_hlt = HeaderletHDU.fromheaderlet(self)
+                new_hlt.update_ext_version(numhlt + 1)
+                fobj.append(new_hlt)
+
+            if close_dest:
+                fobj.close()
+        else:
+            self.hdr_logger.critical("Observation %s cannot be updated with headerlet "
+                            "%s" % (fobj.filename(), self.hdrname))
+            print "Observation %s cannot be updated with headerlet %s" \
+                  % (fobj.filename(), self.hdrname)
+
+
+    def hverify(self):
+        self.verify()
+        assert(self[0].header.has_key('DESTIM') and
+               self[0].header['DESTIM'].strip())
+        assert(self[0].header.has_key('HDRNAME') and
+               self[0].header['HDRNAME'].strip())
+        assert(self[0].header.has_key('STWCSVER'))
+
+
+    def verify_dest(self, dest):
+        """
+        verifies that the headerlet can be applied to the observation
+
+        DESTIM in the primary header of the headerlet must match ROOTNAME
+        of the science file (or the name of the destination file)
+        """
+
+        try:
+            if not isinstance(dest, pyfits.HDUList):
+                droot = pyfits.getval(dest, 'ROOTNAME')
+            else:
+                droot = dest[0].header['ROOTNAME']
+        except KeyError:
+            self.hdr_logger.debug("Keyword 'ROOTNAME' not found in destination file")
+            droot = dest.split('.fits')[0]
+        if droot == self.destim:
+            self.hdr_logger.debug("verify_destim() returned True")
+            return True
+        else:
+            self.hdr_logger.debug("verify_destim() returned False")
+            return False
+
+    def tofile(self, fname, destim=None, hdrname=None, clobber=False):
+        if not destim or not hdrname:
+            self.hverify()
+        self.writeto(fname, clobber=clobber)
+
+    def _delDestWCS(self, dest):
+        """
+        Delete the WCS of a science file
+        """
+
+        self.hdr_logger.info("Deleting all WCSs of file %s" % dest.filename())
+        numext = len(dest)
+
+        for idx in range(numext):
+            # Only delete WCS from extensions which may have WCS keywords
+            if dest[idx].__dict__.has_key('_xtn') and "IMAGE" in dest[idx]._xtn:
+                self._removeD2IM(dest[idx])
+                self._removeSIP(dest[idx])
+                self._removeLUT(dest[idx])
+                self._removePrimaryWCS(dest[idx])
+                self._removeIDCCoeffs(dest[idx])
+                try:
+                    del dest[idx].header.ascard['VAFACTOR']
+                except KeyError:
+                    pass
+                
+        self._removeRefFiles(dest[0])
+        self._removeAltWCS(dest, ext=range(numext))
+        numwdvarr = countExtn(dest, 'WCSDVARR')
+        numd2im = countExtn(dest, 'D2IMARR')
+        for idx in range(1, numwdvarr + 1):
+            del dest[('WCSDVARR', idx)]
+        for idx in range(1, numd2im + 1):
+            del dest[('D2IMARR', idx)]
+    
+    def _removeRefFiles(self, phdu):
+        """
+        phdu: Primary HDU
+        """
+        refkw = ['IDCTAB', 'NPOLFILE', 'D2IMFILE']
+        for kw in refkw:
+            try:
+                del phdu.header.ascard[kw]
+            except KeyError:
+                pass
+        
+    def _removeSIP(self, ext):
+        """
+        Remove the SIP distortion of a FITS extension
+        """
+
+        self.hdr_logger.debug("Removing SIP distortion from (%s, %s)"
+                     % (ext.name, ext._extver))
+        for prefix in ['A', 'B', 'AP', 'BP']:
+            try:
+                order = ext.header[prefix + '_ORDER']
+                del ext.header[prefix + '_ORDER']
+            except KeyError:
+                continue
+            for i in range(order + 1):
+                for j in range(order + 1):
+                    key = prefix + '_%d_%d' % (i, j)
+                    try:
+                        del ext.header[key]
+                    except KeyError:
+                        pass
+        try:
+            del ext.header['IDCTAB']
+        except KeyError:
+            pass
+
+    def _removeLUT(self, ext):
+        """
+        Remove the Lookup Table distortion of a FITS extension
+        """
+
+        self.hdr_logger.debug("Removing LUT distortion from (%s, %s)"
+                     % (ext.name, ext._extver))
+        try:
+            cpdis = ext.header['CPDIS*']
+        except KeyError:
+            return
+        try:
+            for c in range(1, len(cpdis) + 1):
+                del ext.header['DP%s.*...' % c]
+                del ext.header[cpdis[c - 1].key]
+            del ext.header['CPERR*']
+            del ext.header['NPOLFILE']
+            del ext.header['NPOLEXT']
+        except KeyError:
+            pass
+
+    def _removeD2IM(self, ext):
+        """
+        Remove the Detector to Image correction of a FITS extension
+        """
+
+        self.hdr_logger.debug("Removing D2IM correction from (%s, %s)"
+                     % (ext.name, ext._extver))
+        d2imkeys = ['D2IMFILE', 'AXISCORR', 'D2IMEXT', 'D2IMERR']
+        for k in d2imkeys:
+            try:
+                del ext.header[k]
+            except KeyError:
+                pass
+
+    def _removeAltWCS(self, dest, ext):
+        """
+        Remove Alternate WCSs of a FITS extension.
+        A WCS with wcskey 'O' is never deleted.
+        """
+        dkeys = altwcs.wcskeys(dest[('SCI', 1)].header)
+        self.hdr_logger.debug("Removing alternate WCSs with keys %s from %s"
+                     % (dkeys, dest.filename()))
+        for k in dkeys:
+            altwcs.deleteWCS(dest, ext=ext, wcskey=k)
+
+    def _removePrimaryWCS(self, ext):
+        """
+        Remove the primary WCS of a FITS extension
+        """
+
+        self.hdr_logger.debug("Removing Primary WCS from (%s, %s)"
+                     % (ext.name, ext._extver))
+        naxis = ext.header.ascard['NAXIS'].value
+        for key in basic_wcs:
+            for i in range(1, naxis + 1):
+                try:
+                    del ext.header.ascard[key + str(i)]
+                except KeyError:
+                    pass
+        try:
+            del ext.header.ascard['WCSAXES']
+        except KeyError:
+            pass
+
+    def _removeIDCCoeffs(self, ext):
+        """
+        Remove IDC coefficients of a FITS extension
+        """
+
+        self.hdr_logger.debug("Removing IDC coefficient from (%s, %s)"
+                     % (ext.name, ext._extver))
+        coeffs = ['OCX10', 'OCX11', 'OCY10', 'OCY11', 'IDCSCALE']
+        for k in coeffs:
+            try:
+                del ext.header.ascard[k]
+            except KeyError:
+                pass
+
+
+class HeaderletHDU(pyfits.core._NonstandardExtHDU):
+    """
+    A non-standard extension HDU for encapsulating Headerlets in a file.  These
+    HDUs have an extension type of HDRLET and their EXTNAME is derived from the
+    Headerlet's HDRNAME.
+
+    The data itself is a tar file containing a single file, which is the
+    Headerlet file itself.  The file name is derived from the HDRNAME keyword,
+    and should be in the form `<HDRNAME>_hdr.fits`.  If the COMPRESS keyword
+    evaluates to `True`, the tar file is compressed with gzip compression.
+
+    The Headerlet contained in the HDU's data can be accessed by the
+    `headerlet` attribute.
+    """
+
+    _xtn = _extension = 'HDRLET'
+
+    def __init__(self, data=None, header=None):
+        super(HeaderletHDU, self).__init__(data=data, header=header)
+        # TODO: This can be removed after the next pyfits release, but for now
+        # the _ExtensionHDU base class sets self._xtn = '' in its __init__().
+        self._xtn = self._extension
+        # For some reason _NonstandardExtHDU.__init__ sets self.name = None,
+        # even if it's already been set by the EXTNAME keyword in
+        # _ExtensionHDU.__init__() -_-;
+        if header and header.has_key('EXTNAME') and not self.name:
+            self.name = header['EXTNAME']
+        # self._extver, if set, is still preserved.  From
+        # _ExtensionHDU.__init__()...totally inconsistent.
+
+    def __getattr__(self, attr):
+        if attr == 'data':
+            size = self.size()
+            self._file.seek(self._datLoc)
+            self.__dict__[attr] = self._file.read(size)
+        elif attr == 'headerlet':
+            self._file.seek(self._datLoc)
+            s = StringIO()
+            # Read the data into a StringIO--reading directly from the file
+            # won't work (at least for gzipped files) due to problems deep
+            # within the gzip module that make it difficult to read gzip files
+            # embedded in another file
+            s.write(self._file.read(self.size()))
+            s.seek(0)
+            if self._header['COMPRESS']:
+                mode = 'r:gz'
+            else:
+                mode = 'r'
+            t = tarfile.open(mode=mode, fileobj=s)
+            members = t.getmembers()
+            if not len(members):
+                raise ValueError('The Headerlet contents are missing.')
+            elif len(members) > 1:
+                warnings.warn('More than one file is contained in this '
+                              'only the headerlet file should be present.')
+            hlt_name = self._header['HDRNAME'] + '_hdr.fits'
+            try:
+                hlt_info = t.getmember(hlt_name)
+            except KeyError:
+                warnings.warn('The file %s was missing from the HDU data.  '
+                              'Assuming that the first file in the data is '
+                              'headerlet file.' % hlt_name)
+                hlt_info = members[0]
+            hlt_file = t.extractfile(hlt_info)
+            # hlt_file is a file-like object
+            hlt = Headerlet(hlt_file, mode='readonly')
+            self.__dict__[attr] = hlt
+        else:
+            return pyfits.core._ValidHDU.__getattr__(self, attr)
+        try:
+            return self.__dict__[attr]
+        except KeyError:
+            raise AttributeError(attr)
+
+    @classmethod
+    def fromheaderlet(cls, headerlet, compress=False):
+        """
+        Creates a new HeaderletHDU from a given Headerlet object.
+
+        Parameters
+        ----------
+        headerlet : Headerlet
+            A valid Headerlet object.
+        compress : bool, optional
+            Gzip compress the headerlet data.
+        """
+
+        phdu = headerlet[0]
+        phduhdr = phdu.header
+        hlt_filename = phdu.header['HDRNAME'] + '_hdr.fits'
+
+        # TODO: As it stands there's no good way to write out an HDUList in
+        # memory, since it automatically closes the given file-like object when
+        # it's done writing.  I'd argue that if passed an open file handler it
+        # should not close it, but for now we'll have to write to a temp file.
+        fd, name = tempfile.mkstemp()
+        try:
+            f = os.fdopen(fd, 'rb+')
+            headerlet.writeto(f)
+            # The tar file itself we'll write in memory, as it should be
+            # relatively small
+            if compress:
+                mode = 'w:gz'
+            else:
+                mode = 'w'
+            s = StringIO()
+            t = tarfile.open(mode=mode, fileobj=s)
+            t.add(name, arcname=hlt_filename)
+            t.close()
+        finally:
+            os.remove(name)
+
+        cards = [
+            pyfits.Card('XTENSION', cls._extension, 'Headerlet extension'),
+            pyfits.Card('BITPIX',    8, 'array data type'),
+            pyfits.Card('NAXIS',     1, 'number of array dimensions'),
+            pyfits.Card('NAXIS1',    len(s.getvalue()), 'Axis length'),
+            pyfits.Card('PCOUNT',    0, 'number of parameters'),
+            pyfits.Card('GCOUNT',    1, 'number of groups'),
+            pyfits.Card('EXTNAME', cls._extension,
+                        'name of the headerlet extension'),
+            phdu.header.ascard['HDRNAME'],
+            phdu.header.ascard['DATE'],
+            pyfits.Card('SIPNAME', headerlet['SIPWCS', 1].header['WCSNAMEA'],
+                        'SIP distortion model name'),
+            phdu.header.ascard['NPOLFILE'],
+            phdu.header.ascard['D2IMFILE'],
+            pyfits.Card('COMPRESS', compress, 'Uses gzip compression')
+        ]
+
+        header = pyfits.Header(pyfits.CardList(cards))
+
+        hdu = cls(data=pyfits.DELAYED, header=header)
+        hdu._file = s
+        hdu._datLoc = 0
+        return hdu
+
+    @classmethod
+    def match_header(cls, header):
+        """
+        This is a class method used in the pyfits refactoring branch to
+        recognize whether or not this class should be used for instantiating
+        an HDU object based on values in the header.
+
+        It is included here for forward-compatibility.
+        """
+
+        card = header.ascard[0]
+        if card.key != 'XTENSION':
+            return False
+        xtension = card.value.rstrip()
+        return xtension == cls._extension
+
+    # TODO: Add header verification
+
+    def _summary(self):
+        # TODO: Perhaps make this more descriptive...
+        return '%-10s  %-12s  %4d' % (self.name, self.__class__.__name__,
+                                      len(self._header.ascard))
+
+# Monkey-patch pyfits to add minimal support for HeaderletHDUs
+# TODO: Get rid of this ASAP!!! (it won't be necessary with the pyfits
+# refactoring branch)
+if not hasattr(pyfits.Header._updateHDUtype, '_patched'):
+    __old_updateHDUtype = pyfits.Header._updateHDUtype
+    def __updateHDUtype(self):
+        if HeaderletHDU.match_header(self):
+            self._hdutype = HeaderletHDU
+        else:
+            __old_updateHDUtype(self)
+    __updateHDUtype._patched = True
+    pyfits.Header._updateHDUtype = __updateHDUtype
diff --git a/lib/stwcs/wcsutil/hstwcs.py b/lib/stwcs/wcsutil/hstwcs.py
new file mode 100644
index 0000000..11c1f42
--- /dev/null
+++ b/lib/stwcs/wcsutil/hstwcs.py
@@ -0,0 +1,451 @@
+from __future__ import division # confidence high
+
+import os.path
+from pywcs import WCS
+import pyfits
+import instruments
+from stwcs.distortion import models, coeff_converter
+import altwcs
+import numpy as np
+from stsci.tools import fileutil
+from stsci.tools.fileutil import DEGTORAD, RADTODEG
+
+import getinput
+import mappings
+from mappings import inst_mappings, ins_spec_kw
+from mappings import basic_wcs
+
+
+__docformat__ = 'restructuredtext'
+
+class HSTWCS(WCS):
+
+    def __init__(self, fobj=None, ext=None, minerr=0.0, wcskey=" "):
+        """
+        Create a WCS object based on the instrument.
+
+        In addition to basic WCS keywords this class provides
+        instrument specific information needed in distortion computation.
+
+        Parameters
+        ----------
+        fobj: string or PyFITS HDUList object or None
+                a file name, e.g j9irw4b1q_flt.fits
+                a fully qualified filename[EXTNAME,EXTNUM], e.g. j9irw4b1q_flt.fits[sci,1]
+                a pyfits file object, e.g pyfits.open('j9irw4b1q_flt.fits'), in which case the
+                user is responsible for closing the file object.
+        ext:  int or None
+                extension number
+                if ext==None, it is assumed the data is in the primary hdu
+        minerr: float
+                minimum value a distortion correction must have in order to be applied.
+                If CPERRja, CQERRja are smaller than minerr, the corersponding
+                distortion is not applied.
+        wcskey: str
+                A one character A-Z or " " used to retrieve and define an
+                alternate WCS description.
+        """
+
+        self.inst_kw = ins_spec_kw
+        self.minerr = minerr
+        self.wcskey = wcskey
+
+        if fobj is not None:
+            filename, hdr0, ehdr, phdu = getinput.parseSingleInput(f=fobj,
+                                                                   ext=ext)
+            self.filename = filename
+            instrument_name = hdr0.get('INSTRUME', 'DEFAULT')
+            if instrument_name in ['IRAF/ARTDATA','',' ','N/A']:
+                self.instrument = 'DEFAULT'
+            else:
+                self.instrument = instrument_name
+            WCS.__init__(self, ehdr, fobj=phdu, minerr=self.minerr,
+                         key=self.wcskey)
+            # If input was a pyfits HDUList object, it's the user's
+            # responsibility to close it, otherwise, it's closed here.
+            if not isinstance(fobj, pyfits.HDUList):
+                phdu.close()
+            self.setInstrSpecKw(hdr0, ehdr)
+            self.readIDCCoeffs(ehdr)
+            extname = ehdr.get('EXTNAME', '')
+            extnum = ehdr.get('EXTVER', None)
+            self.extname = (extname, extnum)
+        else:
+            # create a default HSTWCS object
+            self.instrument = 'DEFAULT'
+            WCS.__init__(self, minerr=self.minerr, key=self.wcskey)
+            self.pc2cd()
+            self.setInstrSpecKw()
+        self.setPscale()
+        self.setOrient()
+
+    def readIDCCoeffs(self, header):
+        """
+        Reads in first order IDCTAB coefficients if present in the header
+        """
+        coeffs = ['ocx10', 'ocx11', 'ocy10', 'ocy11', 'idcscale']
+        for c in coeffs:
+            self.__setattr__(c, header.get(c, None))
+
+    def setInstrSpecKw(self, prim_hdr=None, ext_hdr=None):
+        """
+        Populate the instrument specific attributes:
+
+        These can be in different headers but each instrument class has knowledge
+        of where to look for them.
+
+        Parameters
+        ----------
+        prim_hdr: pyfits.Header
+                  primary header
+        ext_hdr:  pyfits.Header
+                  extension header
+
+        """
+        if self.instrument in inst_mappings.keys():
+            inst_kl = inst_mappings[self.instrument]
+            inst_kl = instruments.__dict__[inst_kl]
+            insobj = inst_kl(prim_hdr, ext_hdr)
+
+            for key in self.inst_kw:
+                try:
+                    self.__setattr__(key, insobj.__getattribute__(key))
+                except AttributeError:
+                    # Some of the instrument's attributes are recorded in the primary header and
+                    # were already set, (e.g. 'DETECTOR'), the code below is a check for that case.
+                    if not self.__getattribute__(key):
+                        raise
+                    else:
+                        pass
+
+        else:
+            raise KeyError, "Unsupported instrument - %s" %self.instrument
+
+    def setPscale(self):
+        """
+        Calculates the plate scale from the CD matrix
+        """
+        try:
+            cd11 = self.wcs.cd[0][0]
+            cd21 = self.wcs.cd[1][0]
+            self.pscale = np.sqrt(np.power(cd11,2)+np.power(cd21,2)) * 3600.
+        except AttributeError:
+            print "This file has a PC matrix. You may want to convert it \n \
+            to a CD matrix, if reasonable, by running pc2.cd() method.\n \
+            The plate scale can be set then by calling setPscale() method.\n"
+            self.pscale = None
+
+    def setOrient(self):
+        """
+        Computes ORIENTAT from the CD matrix
+        """
+        try:
+            cd12 = self.wcs.cd[0][1]
+            cd22 = self.wcs.cd[1][1]
+            self.orientat = RADTODEG(np.arctan2(cd12,cd22))
+        except AttributeError:
+            print "This file has a PC matrix. You may want to convert it \n \
+            to a CD matrix, if reasonable, by running pc2.cd() method.\n \
+            The orientation can be set then by calling setOrient() method.\n"
+            self.pscale = None
+
+    def updatePscale(self, scale):
+        """
+        Updates the CD matrix with a new plate scale
+        """
+        self.wcs.cd = self.wcs.cd/self.pscale*scale
+        self.setPscale()
+
+    def readModel(self, update=False, header=None):
+        """
+        Reads distortion model from IDCTAB.
+
+        If IDCTAB is not found ('N/A', "", or not found on disk), then
+        if SIP coefficients and first order IDCTAB coefficients are present
+        in the header, restore the idcmodel from the header.
+        If not - assign None to self.idcmodel.
+
+        Parameters
+        ----------
+        header: pyfits.Header
+                fits extension header
+        update: boolean (False)
+                if True - record the following IDCTAB quantities as header keywords:
+                CX10, CX11, CY10, CY11, IDCSCALE, IDCTHETA, IDCXREF, IDCYREF,
+                IDCV2REF, IDCV3REF
+        """
+        if self.idctab in [None, '', ' ','N/A']:
+            #Keyword idctab is not present in header - check for sip coefficients
+            if header is not None and header.has_key('IDCSCALE'):
+                self._readModelFromHeader(header)
+            else:
+                print "Distortion model is not available: IDCTAB=None\n"
+                self.idcmodel = None
+        elif not os.path.exists(fileutil.osfn(self.idctab)):
+            if header is not None and header.has_key('IDCSCALE'):
+                self._readModelFromHeader(header)
+            else:
+                print 'Distortion model is not available: IDCTAB file %s not found\n' % self.idctab
+                self.idcmodel = None
+        else:
+            self.readModelFromIDCTAB(header=header, update=update)
+
+    def _readModelFromHeader(self, header):
+        # Recreate idc model from SIP coefficients and header kw
+        print 'Restoring IDC model from SIP coefficients\n'
+        model = models.GeometryModel()
+        cx, cy = coeff_converter.sip2idc(self)
+        model.cx = cx
+        model.cy = cy
+        model.name = "sip"
+        model.norder = header['A_ORDER']
+
+        refpix = {}
+        refpix['XREF'] = header['IDCXREF']
+        refpix['YREF'] = header['IDCYREF']
+        refpix['PSCALE'] = header['IDCSCALE']
+        refpix['V2REF'] = header['IDCV2REF']
+        refpix['V3REF'] = header['IDCV3REF']
+        refpix['THETA'] = header['IDCTHETA']
+        model.refpix = refpix
+
+        self.idcmodel = model
+
+
+    def readModelFromIDCTAB(self, header=None, update=False):
+        """
+        Read distortion model from idc table.
+
+        Parameters
+        ----------
+        header: pyfits.Header
+                fits extension header
+        update: boolean (False)
+                if True - save teh following as header keywords:
+                CX10, CX11, CY10, CY11, IDCSCALE, IDCTHETA, IDCXREF, IDCYREF,
+                IDCV2REF, IDCV3REF
+
+        """
+        if self.date_obs == None:
+            print 'date_obs not available\n'
+            self.idcmodel = None
+            return
+        if self.filter1 == None and self.filter2 == None:
+            'No filter information available\n'
+            self.idcmodel = None
+            return
+
+        self.idcmodel = models.IDCModel(self.idctab,
+                    chip=self.chip, direction='forward', date=self.date_obs,
+                    filter1=self.filter1, filter2=self.filter2,
+                    offtab=self.offtab, binned=self.binned)
+
+        if update:
+            if header==None:
+                print 'Update header with IDC model kw requested but header was not provided\n.'
+            else:
+                self._updatehdr(header)
+
+    def wcs2header(self, sip2hdr=False, idc2hdr=True):
+        """
+        Create a pyfits.Header object from WCS keywords.
+
+        If the original header had a CD matrix, return a CD matrix,
+        otherwise return a PC matrix.
+
+        Parameters
+        ----------
+        sip2hdr: boolean
+                 If True - include SIP coefficients
+        """
+        h = self.to_header()
+        if self.wcs.has_cd():
+            h = altwcs.pc2cd(h, key=self.wcskey)
+
+        if idc2hdr:
+            for card in self._idc2hdr():
+                h.update(card.key,value=card.value,comment=card.comment)
+        try:
+            del h.ascard['RESTFRQ']
+            del h.ascard['RESTWAV']
+        except KeyError: pass
+
+        if sip2hdr and self.sip:
+            for card in self._sip2hdr('a'):
+                h.update(card.key,value=card.value,comment=card.comment)
+            for card in self._sip2hdr('b'):
+                h.update(card.key,value=card.value,comment=card.comment)
+
+            try:
+                ap = self.sip.ap
+            except AssertionError:
+                ap = None
+            try:
+                bp = self.sip.bp
+            except AssertionError:
+                bp = None
+
+            if ap:
+                for card in self._sip2hdr('ap'):
+                    h.update(card.key,value=card.value,comment=card.comment)
+            if bp:
+                for card in self._sip2hdr('bp'):
+                    h.update(card.key,value=card.value,comment=card.comment)
+        return h
+
+
+    def _sip2hdr(self, k):
+        """
+        Get a set of SIP coefficients in the form of an array
+        and turn them into a pyfits.Cardlist.
+        k - one of 'a', 'b', 'ap', 'bp'
+        """
+
+        cards = pyfits.CardList()
+        korder = self.sip.__getattribute__(k+'_order')
+        cards.append(pyfits.Card(key=k.upper()+'_ORDER', value=korder))
+        coeffs = self.sip.__getattribute__(k)
+        ind = coeffs.nonzero()
+        for i in range(len(ind[0])):
+            card = pyfits.Card(key=k.upper()+'_'+str(ind[0][i])+'_'+str(ind[1][i]),
+                               value=coeffs[ind[0][i], ind[1][i]])
+            cards.append(card)
+        return cards
+
+    def _idc2hdr(self):
+        # save some of the idc coefficients
+        coeffs = ['ocx10', 'ocx11', 'ocy10', 'ocy11', 'idcscale']
+        cards = pyfits.CardList()
+        for c in coeffs:
+            try:
+                val = self.__getattribute__(c)
+            except AttributeError:
+                continue
+            cards.append(pyfits.Card(key=c, value=val))
+        return cards
+
+    def pc2cd(self):
+        self.wcs.cd = self.wcs.pc.copy()
+
+    def all_sky2pix(self,ra,dec,origin):
+        """
+        Performs full inverse transformation using iterative solution
+        on full forward transformation with complete distortion model.
+
+        NOTES
+        -----
+        We now need to find the position we want by iterative
+        improvement of an initial guess - the centre of the chip
+
+        The method is to derive an "effective CD matrix" and use that
+        to apply a correction until we are close enough (as defined by
+        the ERR variable)
+
+        Code from the drizzle task TRANBACK (dither$drizzle/tranback.f)
+        defined the algorithm for this implementation
+
+        """
+        from stwcs.distortion import utils
+
+        # Define some output arrays
+        xout = np.zeros(len(ra),dtype=np.float64)
+        yout = np.zeros(len(ra),dtype=np.float64)
+        # ... and internal arrays
+        x = np.zeros(3,dtype=np.float64)
+        y = np.zeros(3,dtype=np.float64)
+
+        # define delta for comparison
+        err = 0.0001
+
+        # Use linear WCS as frame in which to perform fit
+        # rather than on the sky
+        undistort = True
+        if self.sip is None:
+            # Only apply distortion if distortion coeffs are present.
+            undistort = False
+        wcslin = utils.output_wcs([self],undistort=undistort)
+
+        # We can only transform 1 position at a time
+        for r,d,n in zip(ra,dec,xrange(len(ra))):
+            # First guess for output
+            x[0],y[0] = self.wcs_sky2pix(r,d,origin)
+            # also convert RA,Dec into undistorted linear pixel positions
+            lx,ly = wcslin.wcs_sky2pix(r,d,origin)
+
+            # Loop around until we are close enough (max 20 iterations)
+            ev_old = None
+            for i in xrange(20):
+                x[1] = x[0] + 1.0
+                y[1] = y[0]
+                x[2] = x[0]
+                y[2] = y[0] + 1.0
+                # Perform full transformation on pixel position
+                rao,deco = self.all_pix2sky(x,y,origin)
+                # convert RA,Dec into linear pixel positions for fitting
+                xo,yo = wcslin.wcs_sky2pix(rao,deco,origin)
+
+                # Compute deltas between output and initial guess as
+                # an affine transform then invert that transformation
+                dxymat = np.array([[xo[1]-xo[0],yo[1]-yo[0]],
+                          [xo[2]-xo[0],yo[2]-yo[0]]],dtype=np.float64)
+
+                invmat = np.linalg.inv(dxymat)
+                # compute error in output position
+                dx = lx - xo[0]
+                dy = ly - yo[0]
+
+                # record the old position
+                xold = x[0]
+                yold = y[0]
+
+                # Update the initial guess position using the transform
+                x[0] = xold + dx*dxymat[0][0] + dy*dxymat[1][0]
+                y[0] = yold + dx*dxymat[0][1] + dy*dxymat[1][1]
+
+                # Work out the error vector length
+                ev = np.sqrt((x[0] - xold)**2 + (y[0] - yold)**2)
+
+                # initialize record of previous error measurement during 1st iteration
+                if ev_old is None:
+                    ev_old = ev
+
+                # Check to see whether we have reached the limit or
+                # the new error is greater than error from previous iteration
+                if ev < err or (np.abs(ev) > np.abs(ev_old)):
+                    break
+                # remember error measurement from previous iteration
+                ev_old = ev
+
+            xout[n] = x[0]
+            yout[n] = y[0]
+
+        return xout,yout
+
+    def _updatehdr(self, ext_hdr):
+        #kw2add : OCX10, OCX11, OCY10, OCY11
+        # record the model in the header for use by pydrizzle
+        ext_hdr.update('OCX10', self.idcmodel.cx[1,0])
+        ext_hdr.update('OCX11', self.idcmodel.cx[1,1])
+        ext_hdr.update('OCY10', self.idcmodel.cy[1,0])
+        ext_hdr.update('OCY11', self.idcmodel.cy[1,1])
+        ext_hdr.update('IDCSCALE', self.idcmodel.refpix['PSCALE'])
+        ext_hdr.update('IDCTHETA', self.idcmodel.refpix['THETA'])
+        ext_hdr.update('IDCXREF', self.idcmodel.refpix['XREF'])
+        ext_hdr.update('IDCYREF', self.idcmodel.refpix['YREF'])
+        ext_hdr.update('IDCV2REF', self.idcmodel.refpix['V2REF'])
+        ext_hdr.update('IDCV3REF', self.idcmodel.refpix['V3REF'])
+
+    def printwcs(self):
+        """
+        Print the basic WCS keywords.
+        """
+        print 'WCS Keywords\n'
+        print 'CD_11  CD_12: %r %r' % (self.wcs.cd[0,0],  self.wcs.cd[0,1])
+        print 'CD_21  CD_22: %r %r' % (self.wcs.cd[1,0],  self.wcs.cd[1,1])
+        print 'CRVAL    : %r %r' % (self.wcs.crval[0], self.wcs.crval[1])
+        print 'CRPIX    : %r %r' % (self.wcs.crpix[0], self.wcs.crpix[1])
+        print 'NAXIS    : %d %d' % (self.naxis1, self.naxis2)
+        print 'Plate Scale : %r' % self.pscale
+        print 'ORIENTAT : %r' % self.orientat
+
+
diff --git a/lib/stwcs/wcsutil/instruments.py b/lib/stwcs/wcsutil/instruments.py
new file mode 100644
index 0000000..997bdc8
--- /dev/null
+++ b/lib/stwcs/wcsutil/instruments.py
@@ -0,0 +1,321 @@
+from __future__ import division # confidence high
+
+import pyfits
+from mappings import ins_spec_kw
+    
+class InstrWCS(object):
+    """
+    A base class for instrument specific keyword definition.
+    It prvides a default implementation (modeled by ACS) for
+    all set_kw methods.
+    """
+    def __init__(self, hdr0=None, hdr=None):
+        self.exthdr = hdr   
+        self.primhdr = hdr0
+        self.set_ins_spec_kw()
+        
+    def set_ins_spec_kw(self):
+        """
+        This method MUST call all set_kw methods.
+        There should be a set_kw method for all kw listed in 
+        mappings.ins_spec_kw. TypeError handles the case when 
+        fobj='DEFAULT'.
+        """
+        self.set_idctab()
+        self.set_offtab()
+        self.set_date_obs()
+        self.set_ra_targ()
+        self.set_dec_targ()
+        self.set_pav3()
+        self.set_detector()
+        self.set_filter1()
+        self.set_filter2()
+        self.set_vafactor()
+        self.set_naxis1()
+        self.set_naxis2()
+        self.set_ltv1()
+        self.set_ltv2()
+        self.set_binned()
+        self.set_chip()
+        self.set_parity()
+        
+    def set_idctab(self):
+        try:
+            self.idctab = self.primhdr['IDCTAB']
+        except (KeyError, TypeError):
+            self.idctab = None
+            
+    def set_offtab(self):
+        try:
+            self.offtab = self.primhdr['OFFTAB']
+        except (KeyError, TypeError):
+            self.offtab = None
+
+    def set_date_obs(self):
+        try:
+            self.date_obs = self.primhdr['DATE-OBS']
+        except (KeyError, TypeError):
+            self.date_obs = None
+            
+    def set_ra_targ(self):
+        try:
+            self.ra_targ = self.primhdr['RA-TARG']
+        except (KeyError, TypeError):
+            self.ra_targ = None
+    
+    def set_dec_targ(self):
+        try:
+            self.dec_targ = self.primhdr['DEC-TARG']
+        except (KeyError, TypeError):
+            self.dec_targ = None
+
+    def set_pav3(self):
+        try:
+            self.pav3 = self.primhdr['PA_V3']
+        except (KeyError, TypeError):
+            self.pav3 = None
+            
+    def set_filter1(self):
+        try:
+            self.filter1 = self.primhdr['FILTER1']
+        except (KeyError, TypeError):
+            self.filter1 = None
+            
+    def set_filter2(self):
+        try:
+            self.filter2 = self.primhdr['FILTER2']
+        except (KeyError, TypeError):
+            self.filter2 = None
+            
+    def set_vafactor(self): 
+        try:
+            self.vafactor = self.exthdr['VAFACTOR']
+        except (KeyError, TypeError):
+            self.vafactor = 1
+            
+    def set_naxis1(self):
+        try:
+            self.naxis1 = self.exthdr['naxis1']
+        except (KeyError, TypeError):
+            try:
+                self.naxis1 = self.exthdr['npix1']
+            except (KeyError, TypeError):
+                self.naxis1 = None
+                
+    def set_naxis2(self):
+        try:
+            self.naxis2 = self.exthdr['naxis2']
+        except (KeyError, TypeError):
+            try:
+                self.naxis2 = self.exthdr['npix2']
+            except (KeyError, TypeError):
+                self.naxis2 = None
+                
+    def set_ltv1(self):
+        try:
+            self.ltv1 = self.exthdr['LTV1']
+        except (KeyError, TypeError):
+            self.ltv1 = 0.0
+        
+    def set_ltv2(self):
+        try:
+            self.ltv2 = self.exthdr['LTV2']
+        except (KeyError, TypeError):
+            self.ltv2 = 0.0
+            
+    def set_binned(self):
+        try:
+            self.binned = self.exthdr['BINAXIS1']
+        except (KeyError, TypeError):
+            self.binned = 1
+        
+    def set_chip(self):
+        try:
+            self.chip = self.exthdr['CCDCHIP']
+        except (KeyError, TypeError):
+            self.chip = 1
+                  
+    def set_parity(self):
+        self.parity = [[1.0,0.0],[0.0,-1.0]]
+    
+    def set_detector(self):
+        # each instrument has a different kw for detector and it can be 
+        # in a different header, so this is to be handled by the instrument classes 
+        self.detector = 'DEFAULT'
+        
+class ACSWCS(InstrWCS):
+    """
+    get instrument specific kw    
+    """
+    
+    def __init__(self, hdr0, hdr):
+        self.primhdr = hdr0
+        self.exthdr = hdr
+        InstrWCS.__init__(self,hdr0, hdr)
+        self.set_ins_spec_kw()
+        
+    def set_detector(self):
+        try:
+            self.detector = self.primhdr['DETECTOR']  
+        except KeyError:
+            print 'ERROR: Detector kw not found.\n'
+            raise
+        
+    def set_parity(self):
+        parity = {'WFC':[[1.0,0.0],[0.0,-1.0]],
+                'HRC':[[-1.0,0.0],[0.0,1.0]],
+                'SBC':[[-1.0,0.0],[0.0,1.0]]}
+        
+        if self.detector not in parity.keys():
+            parity = InstrWCS.set_parity(self)
+        else:
+            self.parity = parity[self.detector]
+    
+        
+class WFPC2WCS(InstrWCS):   
+
+
+    def __init__(self, hdr0, hdr):
+        self.primhdr = hdr0
+        self.exthdr = hdr
+        InstrWCS.__init__(self,hdr0, hdr)
+        self.set_ins_spec_kw()
+    
+    def set_filter1(self):
+        self.filter1 = self.primhdr.get('FILTNAM1', None)
+        if self.filter1 == " " or self.filter1 == None:
+            self.filter1 = 'CLEAR1'
+
+    def set_filter2(self):
+        self.filter2 = self.primhdr.get('FILTNAM2', None)
+        if self.filter2 == " " or self.filter2 == None:
+            self.filter2 = 'CLEAR2'
+            
+    
+    def set_binned(self):
+        mode = self.primhdr.get('MODE', 1)
+        if mode == 'FULL':
+            self.binned = 1
+        elif mode == 'AREA':
+            self.binned = 2
+
+    def set_chip(self):
+        self.chip = self.exthdr.get('DETECTOR', 1)
+    
+    def set_parity(self):
+        self.parity = [[-1.0,0.],[0.,1.0]]
+        
+    def set_detector(self):
+        try:
+            self.detector = self.exthdr['DETECTOR']
+        except KeyError:
+            print 'ERROR: Detector kw not found.\n'
+            raise
+    
+
+class WFC3WCS(InstrWCS):
+    """
+    Create a WFC3 detector specific class
+    """
+    
+    def __init__(self, hdr0, hdr):
+        self.primhdr = hdr0
+        self.exthdr = hdr
+        InstrWCS.__init__(self,hdr0, hdr)
+        self.set_ins_spec_kw()
+    
+    def set_detector(self):
+        try:
+            self.detector = self.primhdr['DETECTOR']  
+        except KeyError:
+            print 'ERROR: Detector kw not found.\n'
+            raise
+        
+    def set_filter1(self):
+        self.filter1 = self.primhdr.get('FILTER', None)
+        if self.filter1 == " " or self.filter1 == None:
+            self.filter1 = 'CLEAR'
+    
+    def set_filter2(self):
+        #Nicmos idc tables do not allow 2 filters.
+        self.filter2 = 'CLEAR'
+        
+    def set_parity(self):
+        parity = {'UVIS':[[-1.0,0.0],[0.0,1.0]], 
+          'IR':[[-1.0,0.0],[0.0,1.0]]}
+        
+        if self.detector not in parity.keys():
+            parity = InstrWCS.set_parity(self)
+        else:
+            self.parity = parity[self.detector]
+            
+class NICMOSWCS(InstrWCS):
+    """
+    Create a NICMOS specific class
+    """
+    
+    def __init__(self, hdr0, hdr):
+        self.primhdr = hdr0
+        self.exthdr = hdr
+        InstrWCS.__init__(self,hdr0, hdr)
+        self.set_ins_spec_kw()
+    
+    def set_parity(self):
+        self.parity = [[-1.0,0.],[0.,1.0]]
+        
+    def set_filter1(self):
+        self.filter1 = self.primhdr.get('FILTER', None)
+        if self.filter1 == " " or self.filter1 == None:
+            self.filter1 = 'CLEAR'
+
+    def set_filter2(self):
+        #Nicmos idc tables do not allow 2 filters.
+        self.filter2 = 'CLEAR'
+        
+    def set_chip(self):
+        self.chip = self.detector
+            
+    def set_detector(self):
+        try:
+            self.detector = self.primhdr['CAMERA']  
+        except KeyError:
+            print 'ERROR: Detector kw not found.\n'
+            raise
+    
+class STISWCS(InstrWCS):
+    """
+    A STIS specific class
+    """
+    
+    def __init__(self, hdr0, hdr):
+        self.primhdr = hdr0
+        self.exthdr = hdr
+        InstrWCS.__init__(self,hdr0, hdr)
+        self.set_ins_spec_kw()
+        
+    def set_parity(self):
+        self.parity = [[-1.0,0.],[0.,1.0]]
+        
+    def set_filter1(self):
+        self.filter1 = self.exthdr.get('OPT_ELEM', None)
+        if self.filter1 == " " or self.filter1 == None:
+            self.filter1 = 'CLEAR1'
+
+    def set_filter2(self):
+        self.filter2 = self.exthdr.get('FILTER', None)
+        if self.filter2 == " " or self.filter2 == None:
+            self.filter2 = 'CLEAR2'
+            
+    def set_detector(self):
+        try:
+            self.detector = self.primhdr['DETECTOR']  
+        except KeyError:
+            print 'ERROR: Detector kw not found.\n'
+            raise
+        
+    def set_date_obs(self):
+        try:
+            self.date_obs = self.exthdr['DATE-OBS']
+        except (KeyError, TypeError):
+            self.date_obs = None
+    
\ No newline at end of file
diff --git a/lib/stwcs/wcsutil/mappings.py b/lib/stwcs/wcsutil/mappings.py
new file mode 100644
index 0000000..24038bf
--- /dev/null
+++ b/lib/stwcs/wcsutil/mappings.py
@@ -0,0 +1,29 @@
+from __future__ import division # confidence high
+
+# This dictionary maps an instrument into an instrument class
+# The instrument class handles instrument specific keywords
+
+inst_mappings={'WFPC2': 'WFPC2WCS',
+                'ACS': 'ACSWCS',
+                'NICMOS': 'NICMOSWCS',
+                'STIS': 'STISWCS',
+                'WFC3': 'WFC3WCS',
+                'DEFAULT': 'InstrWCS'
+                }
+
+
+# A list of instrument specific keywords
+# Every instrument class must have methods which define each of these
+# as class attributes.
+ins_spec_kw = [ 'idctab', 'offtab', 'date_obs', 'ra_targ', 'dec_targ', 'pav3', \
+                'detector', 'ltv1', 'ltv2', 'parity', 'binned','vafactor', \
+                'chip', 'naxis1', 'naxis2', 'filter1', 'filter2']
+
+# A list of keywords defined in the primary header.
+# The HSTWCS class sets this as attributes 
+prim_hdr_kw = ['detector', 'offtab', 'idctab', 'date-obs', 
+              'pa_v3', 'ra_targ', 'dec_targ']
+
+# These are the keywords which are archived before MakeWCS is run
+basic_wcs = ['CD1_', 'CD2_', 'CRVAL', 'CTYPE', 'CRPIX', 'CTYPE', 'CDELT', 'CUNIT']
+            
diff --git a/lib/stwcs/wcsutil/mosaic.py b/lib/stwcs/wcsutil/mosaic.py
new file mode 100644
index 0000000..0c02265
--- /dev/null
+++ b/lib/stwcs/wcsutil/mosaic.py
@@ -0,0 +1,183 @@
+from __future__ import division
+import numpy as np
+from matplotlib import pyplot as plt
+import pyfits
+import string
+
+from stsci.tools import parseinput, irafglob
+from stwcs.distortion import utils
+from stwcs import updatewcs, wcsutil
+from stwcs.wcsutil import altwcs
+
+def vmosaic(fnames, outwcs=None, ref_wcs=None, ext=None, extname=None, undistort=True, wkey='V', wname='VirtualMosaic', plot=False, clobber=False):
+    """
+    Create a virtual mosaic using the WCS of the input images.
+
+    Parameters
+    ----------
+    fnames: a string or a list
+              a string or a list of filenames, or a list of wcsutil.HSTWCS objects
+    outwcs: an HSTWCS object
+             if given, represents the output tangent plane
+             if None, the output WCS is calculated from the input observations.
+    ref_wcs: an HSTWCS object
+             if output wcs is not given, this will be used as a reference for the
+             calculation of the output WCS. If ref_wcs is None and outwcs is None,
+             then the first observation in th einput list is used as reference.
+    ext:    an int, a tuple or a list
+              an int - represents a FITS extension, e.g. 0 is the primary HDU
+              a tuple - uses the notation (extname, extver), e.g. ('sci',1)
+              Can be a list of integers or tuples representing FITS extensions
+    extname: string
+              the value of the EXTNAME keyword for the extensions to be used in the mosaic
+    undistort: boolean (default: True)
+               undistort (or not) the output WCS
+    wkey:   string
+              default: 'V'
+              one character A-Z to be used to record the virtual mosaic WCS as
+              an alternate WCS in the headers of the input files.
+    wname:  string
+              default: 'VirtualMosaic
+              a string to be used as a WCSNAME value for the alternate WCS representign
+              the virtual mosaic
+    plot:   boolean
+              if True and matplotlib is installed will make a plot of the tangent plane
+              and the location of the input observations.
+    clobber: boolean
+              This covers the case when an alternate WCS with the requested key
+              already exists in the header of the input files.
+              if clobber is True, it will be overwritten
+              if False, it will compute the new one but will not write it to the headers.
+
+    Notes
+    -----
+    The algorithm is:
+    1. If output WCS is not given it is calculated from the input WCSes.
+       The first image is used as a reference, if no reference is given.
+       This represents the virtual mosaic WCS.
+    2. For each input observation/chip, an HSTWCS object is created
+       and its footprint on the sky is calculated (using only the four corners).
+    3. For each input observation the footprint is projected on the output
+       tangent plane and the virtual WCS is recorded in the header.
+    """
+    wcsobjects = readWCS(fnames, ext, extname)
+    if outwcs != None:
+        outwcs = outwcs.deepcopy()
+    else:
+        if ref_wcs != None:
+            outwcs = utils.output_wcs(wcsobjects, ref_wcs=ref_wcs, undistort=undistort)
+        else:
+            outwcs = utils.output_wcs(wcsobjects, undistort=undistort)
+    if plot:
+        outc=np.array([[0.,0], [outwcs.naxis1,0],
+                             [outwcs.naxis1, outwcs.naxis2],
+                             [0,outwcs.naxis2], [0,0]])
+        plt.plot(outc[:,0], outc[:,1])
+    for wobj in wcsobjects:
+        outcorners = outwcs.wcs_sky2pix(wobj.calcFootprint(),1)
+        if plot:
+            plt.plot(outcorners[:,0], outcorners[:,1])
+        objwcs = outwcs.deepcopy()
+        objwcs.wcs.crpix = objwcs.wcs.crpix - (outcorners[0])
+        updatehdr(wobj.filename, objwcs,wkey=wkey, wcsname=wname, ext=wobj.extname, clobber=clobber)
+    return outwcs
+
+def updatehdr(fname, wcsobj, wkey, wcsname, ext=1, clobber=False):
+    hdr = pyfits.getheader(fname, ext=ext)
+    all_keys = list(string.ascii_uppercase)
+    if wkey.upper() not in all_keys:
+        raise KeyError, "wkey must be one character: A-Z"
+    if wkey not in altwcs.available_wcskeys(hdr):
+        if not clobber:
+            raise ValueError, "wkey %s is already in use. Use clobber=True to overwrite it or specify a different key." %wkey
+        else:
+            altwcs.deleteWCS(fname, ext=ext, wcskey='V')
+    f = pyfits.open(fname, mode='update')
+
+    hwcs = wcs2header(wcsobj)
+    wcsnamekey = 'WCSNAME' + wkey
+    f[ext].header.update(key=wcsnamekey, value=wcsname)
+    for k in hwcs.keys():
+        f[ext].header.update(key=k[:7]+wkey, value=hwcs[k])
+
+    f.close()
+
+def wcs2header(wcsobj):
+
+    h = wcsobj.to_header()
+
+    if wcsobj.wcs.has_cd():
+        altwcs.pc2cd(h)
+    h.update('CTYPE1', 'RA---TAN')
+    h.update('CTYPE2', 'DEC--TAN')
+    norient = np.rad2deg(np.arctan2(h['CD1_2'],h['CD2_2']))
+    #okey = 'ORIENT%s' % wkey
+    okey = 'ORIENT'
+    h.update(key=okey, value=norient)
+    return h
+
+def readWCS(input, exts=None, extname=None):
+    if isinstance(input, str):
+        if input[0] == '@':
+            # input is an @ file
+            filelist = irafglob.irafglob(input)
+        else:
+            try:
+                filelist, output = parseinput.parseinput(input)
+            except IOError: raise
+    elif isinstance(input, list):
+        if isinstance(input[0], wcsutil.HSTWCS):
+            # a list of HSTWCS objects
+            return input
+        else:
+            filelist = input[:]
+    wcso = []
+    fomited = []
+    # figure out which FITS extension(s) to use
+    if exts == None and extname == None:
+        #Assume it's simple FITS and the data is in the primary HDU
+        for f in filelist:
+            try:
+                wcso = wcsutil.HSTWCS(f)
+            except AttributeError:
+                fomited.append(f)
+                continue
+    elif exts != None and validateExt(exts):
+        exts = [exts]
+        for f in filelist:
+            try:
+                wcso.extend([wcsutil.HSTWCS(f, ext=e) for e in exts])
+            except KeyError:
+                fomited.append(f)
+                continue
+    elif extname != None:
+        for f in filelist:
+            fobj = pyfits.open(f)
+            for i in range(len(fobj)):
+                try:
+                    ename = fobj[i].header['EXTNAME']
+                except KeyError:
+                    continue
+                if ename.lower() == extname.lower():
+                    wcso.append(wcsutil.HSTWCS(f,ext=i))
+                else:
+                    continue
+            fobj.close()
+    if fomited != []:
+        print "These files were skipped:"
+        for f in fomited:
+            print f
+    return wcso
+
+
+def validateExt(ext):
+    if not isinstance(ext, int) and not isinstance(ext, tuple) \
+       and not isinstance(ext, list):
+        print "Ext must be integer, tuple, a list of int extension numbers, \
+        or a list of tuples representing a fits extension, for example ('sci', 1)."
+        return False
+    else:
+        return True
+
+
+
diff --git a/lib/stwcs/wcsutil/wcscorr.py b/lib/stwcs/wcsutil/wcscorr.py
new file mode 100644
index 0000000..a6b1f94
--- /dev/null
+++ b/lib/stwcs/wcsutil/wcscorr.py
@@ -0,0 +1,458 @@
+import os,copy
+import pyfits
+import numpy as np
+
+from stsci.tools import fileutil
+import stwcs
+from stwcs.wcsutil import altwcs
+import convertwcs
+
+
+DEFAULT_WCS_KEYS = ['CRVAL1','CRVAL2','CRPIX1','CRPIX2',
+                    'CD1_1','CD1_2','CD2_1','CD2_2',
+                    'CTYPE1','CTYPE2']
+DEFAULT_PRI_KEYS = ['PA_V3']
+###
+### WCSEXT table related keyword archive functions
+###
+def init_wcscorr(input, force=False):
+    """
+    This function will initialize the WCSCORR table if it is not already present,
+    and look for WCS keywords with a prefix of 'O' as the original OPUS
+    generated WCS as the initial row for the table or use the current WCS
+    keywords as initial row if no 'O' prefix keywords are found.
+
+    This function will NOT overwrite any rows already present.
+
+    This function works on all SCI extensions at one time.
+    """
+
+    # TODO: Create some sort of decorator or (for Python2.5) context for
+    # opening a FITS file and closing it when done, if necessary
+    if not isinstance(input, pyfits.HDUList):
+        # input must be a filename, so open as PyFITS object
+        fimg = pyfits.open(input, mode='update')
+        need_to_close = True
+    else:
+        fimg = input
+        need_to_close = False
+
+    # Do not try to generate a WCSCORR table for a simple FITS file
+    if len(fimg) == 1:
+        return 
+    
+    # Verify that a WCSCORR extension does not already exist...
+    for extn in fimg:
+        if extn.header.has_key('extname') and \
+           extn.header['extname'] == 'WCSCORR':
+            if not force:
+                return
+            else:
+                del fimg['WCSCORR']
+    # define the primary columns of the WCSEXT table with initial rows for each
+    # SCI extension for the original OPUS solution
+    numsci = fileutil.countExtn(fimg)
+
+    # create new table with more rows than needed initially to make it easier to
+    # add new rows later
+    wcsext = create_wcscorr(numrows=numsci, padding=numsci * 4)
+    # Assign the correct EXTNAME value to this table extension
+    wcsext.header.update('TROWS', numsci * 2,
+                         comment='Number of updated rows in table')
+    wcsext.header.update('EXTNAME', 'WCSCORR',
+                         comment='Table with WCS Update history')
+    wcsext.header.update('EXTVER', 1)
+
+    used_wcskeys = None
+    wcs1 = stwcs.wcsutil.HSTWCS(fimg,ext=('SCI',1))
+    idc2header = True
+    if wcs1.idcscale is None:
+        idc2header = False
+    wcs_keywords = wcs1.wcs2header(idc2hdr=idc2header).keys()
+
+    # Now copy original OPUS values into table
+    for extver in xrange(1, numsci + 1):
+        rowind = find_wcscorr_row(wcsext.data,
+                                  {'WCS_ID': 'OPUS', 'EXTVER': extver,
+                                   'WCS_key':'O'})
+        # There should only EVER be a single row for each extension with OPUS values
+        rownum = np.where(rowind)[0][0]
+        #print 'Archiving OPUS WCS in row number ',rownum,' in WCSCORR table for SCI,',extver
+
+        hdr = fimg['SCI', extver].header
+        # define set of WCS keywords which need to be managed and copied to the table
+        if used_wcskeys is None:
+            used_wcskeys = altwcs.wcskeys(hdr)
+        # Check to see whether or not there is an OPUS alternate WCS present,
+        # if so, get its values directly, otherwise, archive the PRIMARY WCS
+        # as the OPUS values in the WCSCORR table
+        if 'O' not in used_wcskeys:
+            altwcs.archiveWCS(fimg,('SCI',extver),wcskey='O', wcsname='OPUS')
+        wkey = 'O'
+
+        wcs = stwcs.wcsutil.HSTWCS(fimg, ext=('SCI', extver), wcskey=wkey)
+        wcshdr = wcs.wcs2header(idc2hdr=idc2header)
+
+        for key in DEFAULT_PRI_KEYS:
+            prihdr_keys = []
+            if not hdr.has_key(key):
+                prihdr_keys.append(key)
+
+        if wcsext.data.field('CRVAL1')[rownum] != 0:
+            # If we find values for these keywords already in the table, do not
+            # overwrite them again
+            print 'WCS keywords already updated...'
+            break
+        for key in wcs_keywords:
+            if key in wcsext.data.names:
+                wcsext.data.field(key)[rownum] = wcshdr[(key+wkey)[:8]]
+        # Now get any keywords from PRIMARY header needed for WCS updates
+        for key in prihdr_keys:
+            wcsext.data.field(key)[rownum] = fimg[0].header[key]
+
+    # Now that we have archived the OPUS alternate WCS, remove it from the list
+    # of used_wcskeys
+    if 'O' in used_wcskeys:
+        used_wcskeys.remove('O')
+
+    # Now copy remaining alternate WCSs into table
+    # TODO: Much of this appears to be redundant with update_wcscorr; consider
+    # merging them...
+    for uwkey in used_wcskeys:
+        if wkey == ' ':
+            break
+        for extver in xrange(1, numsci + 1):
+            hdr = fimg['SCI', extver].header
+            wcs = stwcs.wcsutil.HSTWCS(fimg, ext=('SCI', extver),
+                                       wcskey=uwkey)
+            wcshdr = wcs.wcs2header()
+            if 'WCSNAME' + uwkey not in wcshdr:
+                idctab = fileutil.osfn(fimg[0].header['idctab'])
+                idcname = os.path.split(idctab)[-1]
+                idcname = idcname[:idcname.find('_')]
+                wcsid = 'IDC_' + idcname
+            else:
+                wcsid = wcshdr['WCSNAME' + uwkey]
+
+            # identify next empty row
+            rowind = find_wcscorr_row(wcsext.data, selections={'wcs_id': ''})
+            rows = np.where(rowind)
+            if len(rows[0]) > 0:
+                rownum = np.where(rowind)[0][0]
+            else:
+                print 'No available rows found for updating. '
+            #print 'Archiving current WCS row number ',rownum,' in WCSCORR table for SCI,',extver
+
+            # Update selection columns for this row with relevant values
+            wcsext.data.field('WCS_ID')[rownum] = wcsid
+            wcsext.data.field('EXTVER')[rownum] = extver
+            wcsext.data.field('WCS_key')[rownum] = uwkey
+
+            # Look for standard WCS keyword values
+            for key in wcs_keywords:
+                if key in wcsext.data.names:
+                    wcsext.data.field(key)[rownum] = wcshdr[key + uwkey]
+            # Now get any keywords from PRIMARY header needed for WCS updates
+            for key in prihdr_keys:
+                wcsext.data.field(key)[rownum] = fimg[0].header[key]
+
+    # Append this table to the image FITS file
+    fimg.append(wcsext)
+    # force an update now
+    # set the verify flag to 'warn' so that it will always succeed, but still
+    # tell the user if PyFITS detects any problems with the file as a whole
+    fimg.flush('warn')
+
+    if need_to_close:
+        fimg.close()
+
+
+def find_wcscorr_row(wcstab, selections):
+    """
+    Return an array of indices from the table (NOT HDU) 'wcstab' that matches the
+    selections specified by the user.
+
+    The row selection criteria must be specified as a dictionary with
+    column name as key and value(s) representing the valid desired row values.
+    For example, {'wcs_id':'OPUS','extver':2}.
+    """
+
+    mask = None
+    for i in selections:
+        icol = wcstab.field(i)
+        if isinstance(icol,np.chararray): icol = icol.rstrip()
+        bmask = (icol == selections[i])
+        if mask is None:
+            mask = bmask.copy()
+        else:
+            mask = np.logical_and(mask,bmask)
+        del bmask
+    return mask
+
+
+def archive_wcs_file(image, wcs_id=None):
+    """
+    Update WCSCORR table with rows for each SCI extension to record the
+    newly updated WCS keyword values.
+    """
+
+    if not isinstance(image, pyfits.HDUList):
+        fimg = pyfits.open(image, mode='update')
+        close_image = True
+    else:
+        fimg = image
+        close_image = False
+
+    update_wcscorr(fimg, wcs_id=wcs_id)
+
+    if close_image:
+        fimg.close()
+
+
+def update_wcscorr(dest, source=None, extname='SCI', wcs_id=None):
+    """
+    Update WCSCORR table with a new row or rows for this extension header. It
+    copies the current set of WCS keywords as a new row of the table based on
+    keyed WCSs as per Paper I Multiple WCS standard).
+
+    Parameters
+    ----------
+    dest : HDUList
+        The HDU list whose WCSCORR table should be appended to (the WCSCORR HDU
+        must already exist)
+    source : HDUList, optional
+        The HDU list containing the extension from which to extract the WCS
+        keywords to add to the WCSCORR table.  If None, the dest is also used
+        as the source.
+    extname : str, optional
+        The extension name from which to take new WCS keywords.  If there are
+        multiple extensions with that name, rows are added for each extension
+        version.
+    wcs_id : str, optional
+        The name of the WCS to add, as in the WCSNAMEa keyword.  If
+        unspecified, all the WCSs in the specified extensions are added.
+    """
+
+    if source is None:
+        source = dest
+
+    numext = fileutil.countExtn(source, extname)
+    if numext == 0:
+        raise ValueError('No %s extensions found in the source HDU list.'
+                         % extname)
+
+    # Current implementation assumes the same WCS keywords are in each
+    # extension version; if this should not be assumed then this can be
+    # modified...
+    wcs_keys = altwcs.wcskeys(source[(extname, 1)].header)
+    wcs_keys = filter(None, wcs_keys)
+    wcshdr = stwcs.wcsutil.HSTWCS(source, ext=(extname, 1)).wcs2header()
+    wcs_keywords = wcshdr.keys()
+
+    if 'O' in wcs_keys:
+        wcs_keys.remove('O') # 'O' is reserved for original OPUS WCS
+
+    # If we're looking for a particular wcs_id, test ahead of time that it's
+    # actually present in the specified extension headers
+    if wcs_id:
+        wcs_key = ''
+        for wcs_key in wcs_keys:
+            wcsname = source[(extname, 1)].header['WCSNAME' + wcs_key]
+            if wcs_id == wcsname:
+                break
+        else:
+            raise ValueError('A WCS with name %s was not found in the %s '
+                             'extension headers in the source HDU list.'
+                             % (wcs_id, extname))
+        wcs_keys = [wcs_key] # We're only interested in this one
+
+    # create new table for hdr and populate it with the newly updated values
+    new_table = create_wcscorr(numrows=0, padding=len(wcs_keys)*numext)
+    old_table = dest['WCSCORR']
+
+    idx = -1
+    for wcs_key in wcs_keys:
+        for extver in range(1, numext + 1):
+            extn = (extname, extver)
+            hdr = source[extn].header
+            wcsname = hdr['WCSNAME' + wcs_key]
+            selection = {'WCS_ID': wcsname, 'EXTVER': extver,
+                         'WCS_key': wcs_key}
+
+            # Ensure that an entry for this WCS is not already in the dest
+            # table; if so just skip it
+            rowind = find_wcscorr_row(old_table.data, selection)
+            if np.any(rowind):
+                continue
+
+            idx += 1
+
+            wcs = stwcs.wcsutil.HSTWCS(source, ext=extn, wcskey=wcs_key)
+            wcshdr = wcs.wcs2header()
+
+            # Update selection column values
+            for key, val in selection.iteritems():
+                new_table.data.field(key)[idx] = val
+
+            for key in wcs_keywords:
+                if key in new_table.data.names:
+                    new_table.data.field(key)[idx] = wcshdr[key + wcs_key]
+
+            prihdr = source[0].header
+            for key in DEFAULT_PRI_KEYS:
+                if key in new_table.data.names and prihdr.has_key(key):
+                    new_table.data.field(key)[idx] = prihdr[key]
+
+    # If idx was never incremented, no rows were added, so there's nothing else
+    # to do...
+    if idx < 0:
+        return
+
+    # Now, we need to merge this into the existing table
+    rowind = find_wcscorr_row(old_table.data, {'wcs_id':''})
+    old_nrows = np.where(rowind)[0][0]
+    new_nrows = new_table.data.shape[0]
+
+    # check to see if there is room for the new row
+    if (old_nrows + new_nrows) > old_table.data.shape[0]-1:
+        pad_rows = 2 * new_nrows
+        # if not, create a new table with 'pad_rows' new empty rows
+        upd_table = pyfits.new_table(old_table.columns,header=old_table.header,
+                                     nrows=old_table.data.shape[0]+pad_rows)
+    else:
+        upd_table = old_table
+        pad_rows = 0
+    # Now, add
+    for name in old_table.columns.names:
+        # reset the default values to ones specific to the row definitions
+        for i in range(pad_rows):
+            upd_table.data.field(name)[old_nrows+i] = old_table.data.field(name)[-1]
+        # Now populate with values from new table
+        upd_table.data.field(name)[old_nrows:old_nrows + new_nrows] = \
+                new_table.data.field(name)
+    upd_table.header.update('TROWS', old_nrows + new_nrows)
+    
+    # replace old extension with newly updated table extension
+    dest['WCSCORR'] = upd_table
+
+
+def restore_file_from_wcscorr(image, id='OPUS', wcskey=''):
+    """ Copies the values of the WCS from the WCSCORR based on ID specified by user.
+    The default will be to restore the original OPUS-derived values to the Primary WCS.
+    If wcskey is specified, the WCS with that key will be updated instead.
+    """
+
+    if not isinstance(image, pyfits.HDUList):
+        fimg = pyfits.open(image, mode='update')
+        close_image = True
+    else:
+        fimg = image
+        close_image = False
+    numsci = fileutil.countExtn(fimg)
+    wcs_table = fimg['WCSCORR']
+    orig_rows = (wcs_table.data.field('WCS_ID') == 'OPUS')
+    # create an HSTWCS object to figure out what WCS keywords need to be updated
+    wcsobj = stwcs.wcsutil.HSTWCS(fimg,ext=('sci',1))
+    wcshdr = wcsobj.wcs2header()
+    for extn in range(1,numsci+1):
+        # find corresponding row from table
+        ext_rows = (wcs_table.data.field('EXTVER') == extn)
+        erow = np.where(np.logical_and(ext_rows,orig_rows))[0][0]
+        for key in wcshdr:
+            if key in wcs_table.data.names: # insure that keyword is column in table
+                tkey = key
+
+                if 'orient' in key.lower():
+                    key = 'ORIENT'
+                if wcskey == '':
+                    skey = key
+                else:
+                    skey = key[:7]+wcskey
+                fimg['sci',extn].header.update(skey,wcs_table.data.field(tkey)[erow])
+        for key in DEFAULT_PRI_KEYS:
+            if key in wcs_table.data.names:
+                if wcskey == '':
+                    pkey = key
+                else:
+                    pkey = key[:7]+wcskey
+                fimg[0].header.update(pkey,wcs_table.data.field(key)[erow])
+
+    # close the image now that the update has been completed.
+    if close_image:
+        fimg.close()
+
+
+def create_wcscorr(descrip=False, numrows=1, padding=0):
+    """
+    Return the basic definitions for a WCSCORR table.
+    The dtype definitions for the string columns are set to the maximum allowed so
+    that all new elements will have the same max size which will be automatically
+    truncated to this limit upon updating (if needed).
+
+    The table is initialized with rows corresponding to the OPUS solution
+    for all the 'SCI' extensions.
+    """
+
+    trows = numrows + padding
+    # define initialized arrays as placeholders for column data
+    # TODO: I'm certain there's an easier way to do this... for example, simply
+    # define the column names and formats, then create an empty array using
+    # them as a dtype, then create the new table from that array.
+    def_float64_zeros = np.array([0.0] * trows, dtype=np.float64)
+    def_float64_ones = def_float64_zeros + 1.0
+    def_float_col = {'format': 'D', 'array': def_float64_zeros.copy()}
+    def_float1_col = {'format': 'D', 'array':def_float64_ones.copy()}
+    def_str40_col = {'format': '40A',
+                     'array': np.array([''] * trows, dtype='S40')}
+    def_str24_col = {'format': '24A',
+                     'array': np.array([''] * trows, dtype='S24')}
+    def_int32_col = {'format': 'J',
+                     'array': np.array([0]*trows,dtype=np.int32)}
+
+    # If more columns are needed, simply add their definitions to this list
+    col_names = [('CRVAL1',   def_float_col), ('CRVAL2',    def_float_col),
+                 ('CRPIX1',   def_float_col), ('CRPIX2',    def_float_col),
+                 ('CD1_1',    def_float_col), ('CD1_2',     def_float_col),
+                 ('CD2_1',    def_float_col), ('CD2_2',     def_float_col),
+                 ('CTYPE1',   def_str24_col), ('CTYPE2',    def_str24_col),
+                 ('ORIENTAT', def_float_col), ('PA_V3',     def_float_col),
+                 ('Delta_RA', def_float_col), ('Delta_Dec', def_float_col),
+                 ('RMS_RA',   def_float_col), ('RMS_Dec',   def_float_col),
+                 ('Delta_Orientat', def_float_col),
+                 ('Delta_Scale', def_float1_col),
+                 ('NMatch',   def_int32_col), ('Catalog',   def_str40_col)]
+
+    # Define selector columns
+    id_col = pyfits.Column(name='WCS_ID', format='40A',
+                           array=np.array(['OPUS'] * numrows + [''] * padding,
+                                          dtype='S24'))
+    extver_col = pyfits.Column(name='EXTVER', format='I',
+                               array=np.array(range(1, numrows + 1),
+                                              dtype=np.int16))
+    wcskey_col = pyfits.Column(name='WCS_key', format='A',
+                               array=np.array(['O'] * numrows + [''] * padding,
+                                              dtype='S'))
+    # create list of remaining columns to be added to table
+    col_list = [id_col, extver_col, wcskey_col] # start with selector columns
+    
+    for c in col_names:
+        cdef = copy.deepcopy(c[1])
+        col_list.append(pyfits.Column(name=c[0], format=cdef['format'],
+                        array=cdef['array']))
+
+    if descrip:
+        col_list.append(
+            pyfits.Column(name='Descrip', format='128A',
+                          array=np.array(
+                              ['Original WCS computed by OPUS'] * numrows,
+                              dtype='S128')))
+
+    # Now create the new table from the column definitions
+    newtab = pyfits.new_table(pyfits.ColDefs(col_list), nrows=trows)
+    # The fact that setting .name is necessary should be considered a bug in
+    # pyfits.
+    # TODO: Make sure this is fixed in pyfits, then remove this
+    newtab.name = 'WCSCORR'
+
+    return newtab
+