diff options
| -rw-r--r-- | wcsutil/altwcs.py | 134 | ||||
| -rw-r--r-- | wcsutil/headerlet.py | 117 | ||||
| -rw-r--r-- | wcsutil/hstwcs.py | 126 | ||||
| -rw-r--r-- | wcsutil/wcscorr.py | 173 |
4 files changed, 305 insertions, 245 deletions
diff --git a/wcsutil/altwcs.py b/wcsutil/altwcs.py index b483e46..cbeeaef 100644 --- a/wcsutil/altwcs.py +++ b/wcsutil/altwcs.py @@ -1,37 +1,40 @@ from __future__ import division # confidence high +import os +import string -import os.path, string -import pywcs import numpy as np +import pywcs import pyfits -altwcskw = ['WCSAXES', 'CRVAL', 'CRPIX', 'PC', 'CDELT', 'CD', 'CTYPE', 'CUNIT', 'PV', 'PS'] +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' - + 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 + 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: @@ -40,11 +43,11 @@ def archiveWCS(fname, ext, wcskey=" ", wcsname=" ", clobber=False): if not parpasscheck(f, ext, wcskey, wcsname): closefobj(fname, f) return - + if isinstance(ext, int) or isinstance(ext, tuple): ext = [ext] - - if wcskey == " ": + + if wcskey == " ": # try getting the key from WCSNAME if not wcsname.strip(): wkey = next_wcskey(f[1].header) @@ -53,7 +56,7 @@ def archiveWCS(fname, ext, wcskey=" ", wcsname=" ", clobber=False): else: if wcskey not in available_wcskeys(f[1].header): # reuse wcsname - if not wcsname.strip(): + if not wcsname.strip(): wcsname = f[1].header["WCSNAME"+wcskey] wname = wcsname wkey = wcskey @@ -61,9 +64,9 @@ def archiveWCS(fname, ext, wcskey=" ", wcsname=" ", clobber=False): wkey = wcskey wname = wcsname else: - wkey = wcskey + wkey = wcskey wname = wcsname - + for e in ext: w = pywcs.WCS(f[e].header, fobj=f) hwcs = w.to_header() @@ -76,27 +79,27 @@ def archiveWCS(fname, ext, wcskey=" ", wcsname=" ", clobber=False): 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) + #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'_. + 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 + 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)) + 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. @@ -105,27 +108,27 @@ def restoreWCS(f, ext, wcskey=" ", wcsname=" ", clobber=False): 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: + 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: @@ -143,7 +146,7 @@ def restoreWCS(f, ext, wcskey=" ", wcsname=" ", clobber=False): print "Could not find alternate WCS with key %s in this file" % wcskey closefobj(f, fobj) return - wkey = wcskey + wkey = wcskey for e in ext: try: @@ -161,7 +164,7 @@ def restoreWCS(f, ext, wcskey=" ", wcsname=" ", clobber=False): closefobj(f, fobj) return #raise hwcs = nwcs.to_header() - + if nwcs.wcs.has_cd(): pc2cd(hwcs, key=wkey) for k in hwcs.keys(): @@ -187,7 +190,7 @@ def restoreWCS(f, ext, wcskey=" ", wcsname=" ", clobber=False): fobj[e].header.update(key='ORIENTAT', value=norient) else: continue - + if not clobber: fobj.writeto(name) closefobj(f, fobj) @@ -196,7 +199,7 @@ 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 @@ -210,21 +213,21 @@ def deleteWCS(fname, ext, wcskey=" ", wcsname=" "): 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." + print "Wcskey 'O' is reserved for the original WCS and should not be deleted." closefobj(fname, fobj) return - - if wcskey == " ": + + if wcskey == " ": # try getting the key from WCSNAME if wcsname == " ": print "Could not get a valid key from header" @@ -241,8 +244,8 @@ def deleteWCS(fname, ext, wcskey=" ", wcsname=" "): print "Could not find alternate WCS with key %s in this file" % wcskey closefobj(fname, fobj) return - wkey = wcskey - + wkey = wcskey + prexts = [] for i in ext: hdr = fobj[i].header @@ -266,9 +269,9 @@ def deleteWCS(fname, ext, wcskey=" ", wcsname=" "): #header operations def wcskeys(header): """ - Returns a list of characters used in the header for alternate + Returns a list of characters used in the header for alternate WCS description with WCSNAME keyword - + Parameters ---------- hdr: pyfits.Header @@ -280,7 +283,7 @@ def wcskeys(header): def wcsnames(header): """ Returns a dictionary of wcskey: WCSNAME pairs - + Parameters ---------- header: pyfits.Header @@ -290,13 +293,13 @@ def wcsnames(header): 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 @@ -304,9 +307,9 @@ def available_wcskeys(header): assert isinstance(header, pyfits.Header), "Requires a pyfits.Header object as input" all_keys = list(string.ascii_uppercase) used_keys = wcskeys(header) - try: + try: used_keys.remove("") - except ValueError: + except ValueError: pass [all_keys.remove(key) for key in used_keys] return all_keys @@ -314,10 +317,10 @@ def available_wcskeys(header): def next_wcskey(header): """ Returns next available character to be used for an alternate WCS - + Parameters ---------- - header: pyfits.Header + header: pyfits.Header """ assert isinstance(header, pyfits.Header), "Requires a pyfits.Header object as input" allkeys = available_wcskeys(header) @@ -328,10 +331,10 @@ def next_wcskey(header): def getKeyFromName(header, wcsname): """ - If WCSNAME is found in header, return its key, else return + 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 @@ -349,17 +352,17 @@ def getKeyFromName(header, wcsname): 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 + 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: @@ -370,7 +373,7 @@ def pc2cd(hdr, key=' '): val = 1. else: val = 0. - hdr.update(key='CD'+c+'%s' %key, value=val) + hdr.update(key='CD'+c+'%s' %key, value=val) return hdr def parpasscheck(fobj, ext, wcskey, wcsname, clobber=True): @@ -389,12 +392,12 @@ def parpasscheck(fobj, ext, wcskey, wcsname, clobber=True): 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 - - if wcskey == " ": + + if wcskey == " ": # try getting the key from WCSNAME """ if wcsname == " " or wcsname == "": @@ -421,7 +424,7 @@ def parpasscheck(fobj, ext, wcskey, wcsname, clobber=True): return True -def closefobj(fname,f): +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 @@ -429,6 +432,3 @@ def closefobj(fname,f): """ if isinstance(fname, str): f.close() - - -
\ No newline at end of file diff --git a/wcsutil/headerlet.py b/wcsutil/headerlet.py index 96dd71a..c580d1f 100644 --- a/wcsutil/headerlet.py +++ b/wcsutil/headerlet.py @@ -1,21 +1,24 @@ from __future__ import division +import logging import os import tarfile import tempfile import time -import numpy as np import warnings -import pyfits - from cStringIO import StringIO -from hstwcs import HSTWCS +import numpy as np +import pyfits + import altwcs +import wcscorr +from hstwcs import HSTWCS from mappings import basic_wcs +from pytools.fileutil import countExtn -import logging logger = logging.getLogger(__name__) + def isWCSIdentical(scifile, file2, verbose=False): """ Compares the WCS solution of 2 files. @@ -50,8 +53,8 @@ def isWCSIdentical(scifile, file2, verbose=False): logger.info("Starting isWCSIdentical: %s" % time.asctime()) result = True - numsci1 = max(countExt(scifile, 'SCI'), countExt(scifile, 'SIPWCS')) - numsci2 = max(countExt(file2, 'SCI'), countExt(file2, 'SIPWCS')) + numsci1 = max(countExtn(scifile), countExtn(scifile, 'SIPWCS')) + numsci2 = max(countExtn(file2), countExtn(file2, 'SIPWCS')) if numsci1 == 0 or numsci2 == 0 or numsci1 != numsci2: logger.info("Number of SCI and SIPWCS extensions do not match.") @@ -149,7 +152,12 @@ def createHeaderlet(fname, hdrname, destim=None, output=None, verbose=False): phdukw = {'IDCTAB': True, 'NPOLFILE': True, 'D2IMFILE': True} - fobj = pyfits.open(fname) + 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'] @@ -172,7 +180,7 @@ def createHeaderlet(fname, hdrname, destim=None, output=None, verbose=False): "Version of STWCS used to update the WCS") if 'O' in altkeys: altkeys.remove('O') - numsci = countExt(fname, 'SCI') + numsci = countExtn(fname) logger.debug("Number of 'SCI' extensions in file %s is %s" % (fname, numsci)) hdul = pyfits.HDUList() @@ -242,8 +250,8 @@ def createHeaderlet(fname, hdrname, destim=None, output=None, verbose=False): # temporary fix for pyfits ticket # 48 hdu._extver = e hdul.append(hdu) - numwdvarr = countExt(fname, 'WCSDVARR') - numd2im = countExt(fname, 'D2IMARR') + 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 @@ -259,7 +267,8 @@ def createHeaderlet(fname, hdrname, destim=None, output=None, verbose=False): if not output.endswith('_hdr.fits'): output = output + '_hdr.fits' hdul.writeto(output, clobber=True) - fobj.close() + if close_file: + fobj.close() return Headerlet(hdul) def applyHeaderlet(hdrfile, destfile, destim=None, hdrname=None, @@ -341,10 +350,12 @@ def mapFitsExt2HDUListInd(fname, extname): Map FITS extensions with 'EXTNAME' to HDUList indexes. """ - if isinstance(fname, basestring): + 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 @@ -352,26 +363,10 @@ def mapFitsExt2HDUListInd(fname, extname): if hdu.header.has_key('EXTNAME') and hdu.header['EXTNAME'] == extname: extver = hdu.header['EXTVER'] d[(extname, extver)] = f.index_of((extname, extver)) - f.close() + if close_file: + f.close() return d -def countExt(fname, extname): - """ - Return the number of extensions with EXTNAME in the file. - """ - - if isinstance(fname, basestring): - f = pyfits.open(fname) - else: - f = fname - numext = 0 - for hdu in f: - if hdu.header.has_key('EXTNAME') and hdu.header['EXTNAME'] == extname: - numext+=1 - logger.debug("file %s has %s extensions with extname=%s" - % (f.filename(), numext, extname)) - return numext - class Headerlet(pyfits.HDUList): """ @@ -412,22 +407,43 @@ class Headerlet(pyfits.HDUList): self.d2imerr = 0 self.axiscorr = 1 - def apply(self, dest, destim=None, hdrname=None, createheaderlet=True): + def apply(self, dest, createheaderlet=True): """ Apply this headerlet to a file. """ 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 + 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 = fobj['WCSCORR'] + except KeyError: + # The WCSCORR table needs to be created + wcscorr.init_wcscorr(fobj) + if createheaderlet: - # TODO: Currently this does nothing...the created headerlet - # isn't used--it should be given an appropriate name and - # appended to the end of the file. - orig_headerlet = createHeaderlet(dest, hdrname, destim) - fobj = pyfits.open(dest, mode='update') + # Create a headerlet for the original WCS data in the file, + # create an HDU from the original headerlet, and append it to + # the file + hdrname = fobj[0].header['ROOTNAME'] + '_orig' + orig_hlt = createHeaderlet(fobj, hdrname) + orig_hlt_hdu = HeaderletHDU.fromheaderlet(orig_hlt) + fobj.append(orig_hlt_hdu) + self._delDestWCS(fobj) updateRefFiles(self[0].header.ascard, fobj[0].header.ascard) - numsip = countExt(self, 'SIPWCS') + numsip = countExtn(self, 'SIPWCS') for idx in range(1, numsip + 1): fhdr = fobj[('SCI', idx)].header.ascard siphdr = self[('SIPWCS', idx)].header.ascard @@ -454,13 +470,14 @@ class Headerlet(pyfits.HDUList): #! Always attach these extensions last. Otherwise their headers may # get updated with the other WCS kw. - numwdvar = countExt(self, 'WCSDVARR') - numd2im = countExt(self, 'D2IMARR') + 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()) - fobj.close() + if close_dest: + fobj.close() else: logger.critical("Observation %s cannot be updated with headerlet " "%s" % (dest, self.hdrname)) @@ -486,7 +503,10 @@ class Headerlet(pyfits.HDUList): """ try: - droot = pyfits.getval(dest, 'ROOTNAME') + if not isinstance(dest, pyfits.HDUList): + droot = pyfits.getval(dest, 'ROOTNAME') + else: + droot = dest[0].header['ROOTNAME'] except KeyError: logger.debug("Keyword 'ROOTNAME' not found in destination file") droot = dest.split('.fits')[0] @@ -507,8 +527,7 @@ class Headerlet(pyfits.HDUList): Delete the WCS of a science file """ - logger.info("Deleting all WCSs of file %s" - % dest.fileinfo()['filename']) + logger.info("Deleting all WCSs of file %s" % dest.filename()) numext = len(dest) for idx in range(numext): @@ -523,8 +542,8 @@ class Headerlet(pyfits.HDUList): pass self._removeAltWCS(dest, ext=range(numext)) - numwdvarr = countExt(dest, 'WCSDVARR') - numd2im = countExt(dest, 'D2IMARR') + 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): @@ -691,13 +710,13 @@ class HeaderletHDU(pyfits.core._NonstandardExtHDU): elif len(members) > 1: warnings.warn('More than one file is contained in this ' 'only the headerlet file should be present.') - hlt_name = self.name + '_hdr.fits' + 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.') + 'headerlet file.' % hlt_name) hlt_info = members[0] hlt_file = t.extractfile(hlt_info) # hlt_file is a file-like object @@ -755,7 +774,7 @@ class HeaderletHDU(pyfits.core._NonstandardExtHDU): 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', phdu.header['HDRNAME'], + pyfits.Card('EXTNAME', cls._extension, 'name of the headerlet extension'), phdu.header.ascard['HDRNAME'], phdu.header.ascard['DATE'], diff --git a/wcsutil/hstwcs.py b/wcsutil/hstwcs.py index b4e9123..d16e174 100644 --- a/wcsutil/hstwcs.py +++ b/wcsutil/hstwcs.py @@ -17,22 +17,22 @@ 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 + 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 @@ -42,14 +42,14 @@ class HSTWCS(WCS): 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 + 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 != None: filename, hdr0, ehdr, phdu = getinput.parseSingleInput(f=fobj, ext=ext) self.filename = filename @@ -62,14 +62,14 @@ class HSTWCS(WCS): # 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() + 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 + # create a default HSTWCS object self.instrument = 'DEFAULT' WCS.__init__(self, minerr=self.minerr, key=self.wcskey) self.pc2cd() @@ -89,22 +89,22 @@ class HSTWCS(WCS): """ Populate the instrument specific attributes: - These can be in different headers but each instrument class has knowledge + 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)) @@ -115,7 +115,7 @@ class HSTWCS(WCS): raise else: pass - + else: raise KeyError, "Unsupported instrument - %s" %self.instrument @@ -132,7 +132,7 @@ class HSTWCS(WCS): 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 @@ -146,30 +146,30 @@ class HSTWCS(WCS): 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, + CX10, CX11, CY10, CY11, IDCSCALE, IDCTHETA, IDCXREF, IDCYREF, IDCV2REF, IDCV3REF """ if self.idctab in [None, '', ' ','N/A']: @@ -187,7 +187,7 @@ class HSTWCS(WCS): 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' @@ -197,7 +197,7 @@ class HSTWCS(WCS): model.cy = cy model.name = "sip" model.norder = header['A_ORDER'] - + refpix = {} refpix['XREF'] = header['IDCXREF'] refpix['YREF'] = header['IDCYREF'] @@ -206,23 +206,23 @@ class HSTWCS(WCS): 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, + CX10, CX11, CY10, CY11, IDCSCALE, IDCTHETA, IDCXREF, IDCYREF, IDCV2REF, IDCV3REF - + """ if self.date_obs == None: print 'date_obs not available\n' @@ -247,10 +247,10 @@ class HSTWCS(WCS): 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 @@ -259,7 +259,7 @@ class HSTWCS(WCS): 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) @@ -267,7 +267,7 @@ class HSTWCS(WCS): 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) @@ -282,7 +282,7 @@ class HSTWCS(WCS): 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) @@ -290,26 +290,26 @@ class HSTWCS(WCS): 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. + 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]), + 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'] @@ -321,51 +321,51 @@ class HSTWCS(WCS): 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 + """ + 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 + 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 + # 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) @@ -374,7 +374,7 @@ class HSTWCS(WCS): ev_old = None for i in xrange(20): x[1] = x[0] + 1.0 - y[1] = y[0] + y[1] = y[0] x[2] = x[0] y[2] = y[0] + 1.0 # Perform full transformation on pixel position @@ -382,11 +382,11 @@ class HSTWCS(WCS): # 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 + # 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] @@ -402,14 +402,14 @@ class HSTWCS(WCS): # 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 + + # 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)): + if ev < err or (np.abs(ev) > np.abs(ev_old)): break # remember error measurement from previous iteration ev_old = ev @@ -418,7 +418,7 @@ class HSTWCS(WCS): 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 @@ -432,7 +432,7 @@ class HSTWCS(WCS): 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. diff --git a/wcsutil/wcscorr.py b/wcsutil/wcscorr.py index 1ff10c5..9b4a39e 100644 --- a/wcsutil/wcscorr.py +++ b/wcsutil/wcscorr.py @@ -16,14 +16,15 @@ 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 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 will NOT overwrite any rows already present. - This function works on all SCI extensions at one time. + This function works on all SCI extensions at one time. """ # TODO: Create some sort of decorator or (for Python2.5) context for @@ -38,7 +39,8 @@ def init_wcscorr(input, force=False): # 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 extn.header.has_key('extname') and \ + extn.header['extname'] == 'WCSCORR': if not force: return else: @@ -49,10 +51,12 @@ def init_wcscorr(input, force=False): # 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) + 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('TROWS', numsci * 2, + comment='Number of updated rows in table') + wcsext.header.update('EXTNAME', 'WCSCORR', + comment='Table with WCS Update history') used_wcskeys = None wcs1 = stwcs.wcsutil.HSTWCS(fimg,ext=('SCI',1)) @@ -61,13 +65,15 @@ def init_wcscorr(input, force=False): 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'}) + 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 + 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) @@ -77,8 +83,9 @@ def init_wcscorr(input, force=False): if 'O' in used_wcskeys: wkey = 'O' else: - wkey = " " - wcshdr = stwcs.wcsutil.HSTWCS(fimg,ext=('SCI',extver),wcskey=wkey).wcs2header(idc2hdr=idc2header) + wkey = ' ' + wcs = stwcs.wcsutil.HSTWCS(fimg, ext=('SCI', extver), wcskey=wkey) + wcshdr = wcs.wcs2header(idc2hdr=idc2header) for key in DEFAULT_PRI_KEYS: prihdr_keys = [] @@ -99,26 +106,28 @@ def init_wcscorr(input, force=False): # 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") + if 'O' in used_wcskeys: + used_wcskeys.remove('O') # Now copy remaining alternate WCSs into table - for extver in xrange(1,numsci+1): - hdr = fimg['SCI',extver].header + for uwkey in used_wcskeys: if wkey == ' ': break - for uwkey in used_wcskeys: - wcshdr = stwcs.wcsutil.HSTWCS(fimg,ext=('SCI',extver),wcskey=uwkey).wcs2header() - if 'WCSNAME'+uwkey not in wcshdr: - idctab =fileutil.osfn(fimg[0].header['idctab']) + 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+'_'+fileutil.getDate() + wcsid = 'IDC_' + idcname + '_' + fileutil.getDate() else: - wcsid = wcshdr['WCSNAME'+uwkey] + wcsid = wcshdr['WCSNAME' + uwkey] # identify next empty row - rowind = find_wcscorr_row(wcsext.data,selections={'wcs_id':' '}) + rowind = find_wcscorr_row(wcsext.data, selections={'wcs_id': ''}) rows = np.where(rowind) if len(rows[0]) > 0: rownum = np.where(rowind)[0][0] @@ -134,7 +143,7 @@ def init_wcscorr(input, force=False): # 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] + 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] @@ -168,8 +177,9 @@ def find_wcscorr_row(wcstab, selections): 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. + """ + Update WCSCORR table with rows for each SCI extension to record the + newly updated WCS keyword values. """ if not isinstance(image, pyfits.HDUList): @@ -184,29 +194,35 @@ def archive_wcs_file(image, wcs_id=None): if close_image: fimg.close() -def update_wcscorr(fimg,hdr,selections=None): - """ Update WCSCORR table with a new row for this extension header. It +def update_wcscorr(fimg, hdr, selections=None): + """ + Update WCSCORR table with a new row 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). """ + # Now update the table... if selections is None: # define the WCS ID for this update wcs_key = altwcs.wcskeys(hdr)[-1] - selections = {'WCS_ID':'TWEAK_'+fileutil.getDate(),'EXTVER':hdr['extver'],'WCS_key':wcs_key} + selections = {'WCS_ID': 'TWEAK_' + fileutil.getDate(), + 'EXTVER': hdr['extver'], 'WCS_key': wcs_key} if selections == 'IDC': idcname = os.path.split(fileutil.osfn(fimg[0].header['idctab']))[1] - selections = {'WCS_ID':'IDC_'+idcname+'_'+fileutil.getDate(),'EXTVER':hdr['extver'],'WCS_key':wcs_key} + selections = {'WCS_ID': 'IDC_%s_%s' % (idcname, fileutil.getDate()), + 'EXTVER': hdr['extver'], 'WCS_key': wcs_key} # create new table for hdr and populate it with the newly updated values new_table = create_wcscorr() if hdr.has_key('extname'): extname = hdr['extname'] else: extname = 'PRIMARY' - if hdr.has_key('extver'): extver = hdr['extver'] - else: extver = 1 - extn = (extname,extver) - wcshdr = stwcs.wcsutil.HSTWCS(fimg,ext=extn).wcs2header() + if hdr.has_key('extver'): + extver = hdr['extver'] + else: + extver = 1 + extn = (extname, extver) + wcshdr = stwcs.wcsutil.HSTWCS(fimg, ext=extn).wcs2header() # ===> NOTE: need to add checks to insure that these IDs are unique # Update selection column values for key in selections: @@ -222,19 +238,21 @@ def update_wcscorr(fimg,hdr,selections=None): # Now, we need to merge this into the existing table old_table = fimg['WCSCORR'] - rowind = find_wcscorr_row(old_table.data,{'wcs_id':''}) + rowind = find_wcscorr_row(old_table.data, {'wcs_id': ''}) old_nrows = np.where(rowind)[0][0] # check to see if there is room for the new row if (old_nrows + 1) > old_table.data.shape[0]: # if not, create a new table with 'pad_rows' new empty rows - upd_table = pyfits.new_table(old_table.columns,nrows=old_table.data.shape[0]+pad_rows) + upd_table = pyfits.new_table(old_table.columns, + nrows=old_table.data.shape[0] + pad_rows) else: upd_table = old_table # Now, add for name in old_table.columns.names: - upd_table.data.field(name)[old_nrows:old_nrows+1] = new_table.data.field(name) - upd_table.header.update('TROWS',old_nrows+1) + upd_table.data.field(name)[old_nrows:old_nrows+1] = \ + new_table.data.field(name) + upd_table.header.update('TROWS', old_nrows + 1) # replace old extension with newly updated table extension fimg['WCSCORR'] = upd_table @@ -285,7 +303,8 @@ def restore_file_from_wcscorr(image, id='OPUS', wcskey=''): fimg.close() def create_wcscorr(descrip=False, numrows=1, padding=0): - """ Return the basic definitions for a WCSCORR table. + """ + 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). @@ -293,44 +312,66 @@ def create_wcscorr(descrip=False, numrows=1, padding=0): 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 - def_float64_zeros = np.array([0.0]*trows,dtype=np.float64) + # 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_range_1nrows = np.array(range(1,numrows+1),dtype=np.int16) - 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)} + 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)] + 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")) + 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 + 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'])) + 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"))) + 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) + 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 |