diff options
| -rw-r--r-- | lib/utils.py | 200 | ||||
| -rw-r--r-- | updatewcs/__init__.py | 98 | ||||
| -rw-r--r-- | updatewcs/apply_corrections.py | 41 |
3 files changed, 91 insertions, 248 deletions
diff --git a/lib/utils.py b/lib/utils.py index 6f63db0..29ba5f3 100644 --- a/lib/utils.py +++ b/lib/utils.py @@ -1,147 +1,5 @@ from __future__ import division # confidence high -from pytools import parseinput, fileutil -import pyfits -import pywcs -import numpy as np -import string - -def restoreWCS(fnames, wcskey, clobber=False): - """ - Purpose - ======= - 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 representing _'WCSKEY'_. - Otherwise overwrites the files. The WCS is restored from the 'SCI' - extension but the primary WCS of all extensions are updated. - - WCS keywords: - 'CD1_1', 'CD1_2', 'CD2_1', 'CD2_2', - 'CRVAL*', - 'CTYPE*', - 'CRPIX*', - 'CDELT*', - 'CUNIT*', - 'ORIENTAT' - ? - 'TDDALPHA', 'TDDBETA' - 'A_x_x', B_x_x' - SIP coefficients - 'CPERROR*', 'CPDIS*', 'DP*', - - `fnames`: a python list of file names, a string of comma separated file names, - an @file - `wcskey`: a charater - Used for one of 26 alternate WCS definitions. - `clobber`: boolean - A flag to define if the original files should be overwritten - """ - files = parseinput.parseinput(fnames)[0] - for f in files: - isfits, ftype = fileutil.isFits(f) - if not isfits or (isfits and ftype == 'waiver'): - print "RestoreWCS works only with true fits files." - return - else: - if clobber: - print 'Overwriting original files\n' - fobj = pyfits.open(f, mode='update') - name = f - else: - fobj = pyfits.open(f) - name = (f.split('.fits')[0] + '_%s_' + '.fits') %wcskey - for e in range(len(fobj)): - 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=wcskey) - except: - print 'utils.restoreWCS: Could not read WCS with key %s in file %s, \ - extension %d' % (wcskey, f, e) - return #raise - hwcs = nwcs.to_header() - if nwcs.wcs.has_pc(): - for c in ['1_1', '1_2', '2_1', '2_2']: - del hwcs['CD'+c+wcskey] - elif nwcs.wcs.has_cd(): - for c in ['1_1', '1_2', '2_1', '2_2']: - hwcs.update(key='CD'+c+wcskey, value=hwcs['PC'+c+wcskey]) - del hwcs['PC'+c] - 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'): - cd12 = 'CD1_2%s' % wcskey - cd22 = 'CD2_2%s' % wcskey - norient = np.rad2deg(np.arctan2(hwcs[cd12],hwcs[cd22])) - fobj[e].header.update(key='ORIENTAT', value=norient) - elif extname in ['err', 'dq', 'sdq']: - 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'): - cd12 = 'CD1_2%s' % wcskey - cd22 = 'CD2_2%s' % wcskey - norient = np.rad2deg(np.arctan2(hwcs[cd12],hwcs[cd22])) - fobj[e].header.update(key='ORIENTAT', value=norient) - else: - continue - - if not clobber: - fobj.writeto(name) - fobj.close() - -def archiveWCS(fname, ext, wcskey, wcsname=" ", clobber=False): - """ - Copy the primary WCS to an alternate WCS - with wcskey and name WCSNAME. - """ - assert len(wcskey) == 1 - if wcskey == " " and clobber==False: - print "Please provide a valid wcskey for this WCS." - print 'Use "utils.next_wcskey" to obtain a valid wcskey.' - print 'Or use utils.restoreWCS a specific WCS to the primary WCS.' - print 'WCS was NOT archived.' - return - if (wcskey not in available_wcskeys(pyfits.getheader(fname, ext=ext))) and clobber==False: - print 'Wcskey %s is already used in this header.' % wcskey - print 'Use "utils.next_wcskey" to obtain a valid wcskey' - print 'or use "clobber=True" to overwrite the values.' - return - - f = pyfits.open(fname, mode='update') - w = pywcs.WCS(f[ext].header, fobj=f) - hwcs = w.to_header() - wkey = 'WCSNAME' + wcskey - f[ext].header.update(key=wkey, value=wcsname) - if w.wcs.has_pc(): - for c in ['CD1_1', 'CD1_2', 'CD2_1', 'CD2_2']: - del hwcs[c] - elif w.wcs.has_cd(): - for c in ['1_1', '1_2', '2_1', '2_2']: - hwcs.update(key='CD'+c, value=hwcs['PC'+c]) - del hwcs['PC'+c] - for k in hwcs.keys(): - key = k+wcskey - f[ext].header.update(key=key, value = hwcs[k]) - norient = np.rad2deg(np.arctan2(hwcs['CD1_2'],hwcs['CD2_2'])) - okey = 'ORIENT%s' % wcskey - f[ext].header.update(key=okey, value=norient) - f.close() - - def diff_angles(a,b): """ Perform angle subtraction a-b taking into account @@ -168,61 +26,3 @@ def getBinning(fobj, extver=1): binned = fobj['SCI', extver].header.get('BINAXIS',1) return binned -def wcsnames(header): - """ - Purpose - ======= - Return a dictionary of wcskey: WCSNAME pairs - """ - names = header["WCSNAME*"] - d = {} - for c in names: - d[c.key[-1]] = c.value - return d - - -def wcskeys(header): - """ - Purpose - ======= - Returns a list of characters used in the header for alternate - WCS description via WCSNAME keyword - - `header`: pyfits.Header - """ - names = header["WCSNAME*"] - return [key.split('WCSNAME')[1] for key in names.keys()] - -def available_wcskeys(header): - """ - Purpose - ======= - 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. - - `header`: pyfits.Header - """ - 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): - """ - Purpose - ======= - Returns next available character to be used for an alternate WCS - - `header`: pyfits.Header - """ - allkeys = available_wcskeys(header) - if allkeys != []: - return allkeys[0] - else: - return None -
\ No newline at end of file diff --git a/updatewcs/__init__.py b/updatewcs/__init__.py index b69928f..cce9c8a 100644 --- a/updatewcs/__init__.py +++ b/updatewcs/__init__.py @@ -2,9 +2,11 @@ from __future__ import division # confidence high import os import pyfits +import numpy as np +from stwcs import wcsutil from stwcs.wcsutil import HSTWCS +import pywcs -from stwcs import utils import corrections, makewcs import dgeo, det2im from pytools import parseinput, fileutil @@ -14,7 +16,7 @@ import apply_corrections __docformat__ = 'restructuredtext' -__version__ = '0.5' +__version__ = '0.8' def updatewcs(input, vacorr=True, tddcorr=True, dgeocorr=True, d2imcorr=True, checkfiles=True, wcskey=" ", wcsname=" ", clobber=False): @@ -77,7 +79,7 @@ def updatewcs(input, vacorr=True, tddcorr=True, dgeocorr=True, d2imcorr=True, tddcorr=tddcorr,dgeocorr=dgeocorr, d2imcorr=d2imcorr) #restore the original WCS keywords - utils.restoreWCS(f, wcskey='O', clobber=True) + wcsutil.restoreWCS(f, wcskey='O', clobber=True) makecorr(f, acorr, wkey=wcskey, wname=wcsname, clobber=False) return files @@ -106,27 +108,28 @@ def makecorr(fname, allowed_corr, wkey=" ", wname=" ", clobber=False): f = pyfits.open(fname, mode='update') #Determine the reference chip and create the reference HSTWCS object nrefchip, nrefext = getNrefchip(f) - ref_wcs = HSTWCS(fobj=f, ext=nrefext) - ref_wcs.readModel(update=True,header=f[nrefext].header) - ref_wcs.copyWCS(header=f[nrefext].header, wcskey='O', wcsname='OPUS', clobber=True) + rwcs = HSTWCS(fobj=f, ext=nrefext) + rwcs.readModel(update=True,header=f[nrefext].header) + wcsutil.archiveWCS(f, 'O', wcsname='OPUS', ext=nrefext, clobber=True) if 'DET2IMCorr' in allowed_corr: det2im.DET2IMCorr.updateWCS(f) - + + # get a wcskey and wcsname from the first extension header + idcname = fileutil.osfn(rwcs.idctab) + key, name = getKeyName(f[1].header, wkey, wname, idcname) + for i in range(len(f))[1:]: - # Perhaps all ext headers should be corrected (to be consistent) extn = f[i] if extn.header.has_key('extname'): extname = extn.header['extname'].lower() if extname == 'sci': - sciextver = extn.header['extver'] - ref_wcs.restore(f[nrefext].header, wcskey="O") - + ref_wcs = rwcs.deepcopy() hdr = extn.header ext_wcs = HSTWCS(fobj=f, ext=i) - ext_wcs.copyWCS(header=hdr, wcskey='O', wcsname='OPUS', clobber=True) + wcsutil.archiveWCS(f, "O", wcsname="OPUS", ext=[i], clobber=True) ext_wcs.readModel(update=True,header=hdr) for c in allowed_corr: if c != 'DGEOCorr' and c != 'DET2IMCorr': @@ -134,23 +137,15 @@ def makecorr(fname, allowed_corr, wkey=" ", wname=" ", clobber=False): kw2update = corr_klass.updateWCS(ext_wcs, ref_wcs) for kw in kw2update: hdr.update(kw, kw2update[kw]) - - if wkey is not None: - # archive the updated primary WCS - if wkey == " " : - idcname = os.path.split(fileutil.osfn(ext_wcs.idctab))[1] - wname = ''.join(['IDC_',idcname.split('_idc.fits')[0]]) - wkey = getKey(hdr, wname) - #in this case clobber = true, to allow updatewcs to be run repeatedly - ext_wcs.copyWCS(header=hdr, wcskey=wkey, wcsname=wname, clobber=True) - else: - #clobber is set to False as a warning to users - ext_wcs.copyWCS(header=hdr, wcskey=wkey, wcsname=wname, clobber=False) - + #if wkey is None, do not archive the primary WCS + if key is not None: + wcsutil.archiveWCS(f, key, wcsname=name, ext=[i], clobber=True) elif extname in ['err', 'dq', 'sdq', 'samp', 'time']: cextver = extn.header['extver'] if cextver == sciextver: - ext_wcs.copyWCS(header=extn.header, wcskey=" ", wcsname=" ") + hdr = f[('SCI',sciextver)].header + w = pywcs.WCS(hdr, f) + copyWCS(w, extn.header, key, name) else: continue @@ -158,22 +153,49 @@ def makecorr(fname, allowed_corr, wkey=" ", wname=" ", clobber=False): kw2update = dgeo.DGEOCorr.updateWCS(f) for kw in kw2update: f[1].header.update(kw, kw2update[kw]) - + f.close() -def getKey(header, wcsname): +def getKeyName(hdr, wkey, wname, idcname): + if wkey is not None: # archive the primary WCS + if wkey == " ": + if wname == " " : + # get the next available key and use the IDCTABLE name as WCSNAME + idcname = os.path.split(idcname)[1] + name = ''.join(['IDC_',idcname.split('_idc.fits')[0]]) + key = wcsutil.getKeyFromName(hdr, name) + if not key: + key = wcsutil.next_wcskey(hdr) + else: + #try to get a key from WCSNAME + # if not - get the next availabble key + name = wname + key = wcsutil.getKeyFromName(hdr, wname) + if not wkey: + key = wcsutil.next_wcskey(hdr) + else: + key = wkey + name = wname + return key, name + +def copyWCS(w, hdr, wkey, wname): """ - If WCSNAME is found in header, return its key, else return - the next available key. This is used to update a specific WCS - repeatedly and not generate new keys every time. + This is a convenience function to copy a WCS object + to a header as a primary WCS. It is used only to copy the + WCS of the 'SCI' extension to the headers of 'ERR', 'DQ', 'SDQ', + 'TIME' or 'SAMP' extensions. """ - wkey = utils.next_wcskey(header) - names = utils.wcsnames(header) - for item in names.items(): - if item[1] == wcsname: - wkey = item[0] - return wkey - + hwcs = w.to_header() + + if w.wcs.has_cd(): + wcsutil.pc2cd(hwcs) + for k in hwcs.keys(): + key = k+wkey + hdr.update(key=key, value=hwcs[k]) + norient = np.rad2deg(np.arctan2(hwcs['CD1_2'],hwcs['CD2_2'])) + okey = 'ORIENT%s' % wkey + hdr.update(key=okey, value=norient) + def getNrefchip(fobj): """ This handles the fact that WFPC2 subarray observations diff --git a/updatewcs/apply_corrections.py b/updatewcs/apply_corrections.py index af1a081..4eee126 100644 --- a/updatewcs/apply_corrections.py +++ b/updatewcs/apply_corrections.py @@ -48,7 +48,17 @@ def setCorrections(fname, vacorr=True, tddcorr=True, dgeocorr=True, d2imcorr=Tru if 'MakeWCS' in acorr: acorr.remove('MakeWCS') if 'CompSIP' in acorr: acorr.remove('CompSIP') - + # A new IDCTAB means all previously computed WCS's are invalid + # We are deleting all of them except the original OPUS WCS. + if 'MakeWCS' in acorr and newIDCTAB(fname): + print "New IDCTAB file detected. This invalidates all WCS's." + print "Deleting all previous WCS's" + keys = utils.wcskeys(pyfits.getheader(fname, ext=1)) + if 'O' in keys: + keys.remove('O') + for key in keys: + utils.deleteWCS(fname, key) + if 'VACorr' in acorr and vacorr==False: acorr.remove('VACorr') if 'TDDCorr' in acorr: tddcorr = applyTDDCorr(fname, tddcorr) @@ -63,18 +73,16 @@ def setCorrections(fname, vacorr=True, tddcorr=True, dgeocorr=True, d2imcorr=Tru return acorr def foundIDCTAB(fname): - idctab_found = True try: idctab = fileutil.osfn(pyfits.getval(fname, 'IDCTAB')) - if idctab == 'N/A' or idctab == "": - idctab_found = False - if os.path.exists(idctab): - idctab_found = True - else: - idctab_found = False except KeyError: - idctab_found = False - return idctab_found + return False + if idctab == 'N/A' or idctab == "": + return False + if os.path.exists(idctab): + return True + else: + return False def applyTDDCorr(fname, utddcorr): """ @@ -214,4 +222,17 @@ def applyD2ImCorr(fname, d2imcorr): print 'D2IMFILE keyword not found in primary header' applyD2IMCorr = False return applyD2IMCorr + +def newIDCTAB(fname): + #When this is called we know there's a kw IDCTAB in the header + idctab = fileutil.osfn(pyfits.getval(fname, 'IDCTAB')) + try: + #check for the presence of IDCTAB in the first extension + oldidctab = fileutil.osfn(pyfits.getval(fname, 'IDCTAB', ext=1)) + except KeyError: + return False + if idctab == oldidctab: + return False + else: + return True
\ No newline at end of file |