diff options
Diffstat (limited to 'lib/stwcs/wcsutil')
| -rw-r--r-- | lib/stwcs/wcsutil/__init__.py | 35 | ||||
| -rw-r--r-- | lib/stwcs/wcsutil/altwcs.py | 451 | ||||
| -rw-r--r-- | lib/stwcs/wcsutil/convertwcs.py | 118 | ||||
| -rw-r--r-- | lib/stwcs/wcsutil/getinput.py | 62 | ||||
| -rw-r--r-- | lib/stwcs/wcsutil/headerlet.py | 898 | ||||
| -rw-r--r-- | lib/stwcs/wcsutil/hstwcs.py | 451 | ||||
| -rw-r--r-- | lib/stwcs/wcsutil/instruments.py | 321 | ||||
| -rw-r--r-- | lib/stwcs/wcsutil/mappings.py | 29 | ||||
| -rw-r--r-- | lib/stwcs/wcsutil/mosaic.py | 183 | ||||
| -rw-r--r-- | lib/stwcs/wcsutil/wcscorr.py | 458 |
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 + |