diff options
| -rw-r--r-- | updatewcs/__init__.py | 129 | ||||
| -rw-r--r-- | updatewcs/apply_corrections.py | 124 | ||||
| -rw-r--r-- | updatewcs/corrections.py | 27 | ||||
| -rw-r--r-- | updatewcs/dgeo.py | 85 | ||||
| -rw-r--r-- | updatewcs/makewcs.py | 1 |
5 files changed, 199 insertions, 167 deletions
diff --git a/updatewcs/__init__.py b/updatewcs/__init__.py index 2f77949..323b0e6 100644 --- a/updatewcs/__init__.py +++ b/updatewcs/__init__.py @@ -2,20 +2,21 @@ import os import pyfits #from .. wcsutil import HSTWCS from hstwcs.wcsutil import HSTWCS -from hstwcs.mappings import allowed_corrections + #from .. mappings import allowed_corrections from hstwcs import utils import corrections, makewcs import dgeo import time from pytools import parseinput, fileutil +import apply_corrections #Note: The order of corrections is important __docformat__ = 'restructuredtext' -def updatewcs(input, vacorr=True, tddcorr=True, checkfiles=True): +def updatewcs(input, vacorr=True, tddcorr=True, dgeocorr=True, checkfiles=True): """ Purpose ======= @@ -39,7 +40,7 @@ def updatewcs(input, vacorr=True, tddcorr=True, checkfiles=True): `vacorr`: boolean If True, vecocity aberration correction will be applied `tddcorr`: boolean - If True, time dependen correction will be applied to the distortion model + If True, time dependent distortion correction will be applied `checkfiles`: boolean If True, the format of the input files will be checked, geis and waiver fits files will be converted to MEF format. @@ -53,24 +54,7 @@ def updatewcs(input, vacorr=True, tddcorr=True, checkfiles=True): print 'No valid input, quitting ...\n' return for f in files: - instrument = pyfits.getval(f, 'INSTRUME') - if instrument == 'ACS' and pyfits.getval(f, 'DETECTOR')== 'WFC': - try: - tdc = pyfits.getval(f, 'TDDCORR') - if tdc == 'PERFORM': - tddcorr = True - else: - tddcorr = False - except KeyError: - tddcorr = True - try: - #define which corrections will be performed - acorr = setCorrections(instrument,vacorr=vacorr, tddcorr=tddcorr) - except KeyError: - print 'Unsupported instrument %s ' %instr - print 'Removing %s from list of processed files\n' % f - files.remove(f) - continue + acorr = apply_corrections.setCorrections(f, vacorr=vacorr, tddcorr=tddcorr,dgeocorr=dgeocorr) #restore the original WCS keywords utils.restoreWCS(f) makecorr(f, acorr) @@ -80,7 +64,7 @@ def makecorr(fname, allowed_corr): """ Purpose ======= - Applies corrections to a single file + Applies corrections to the WCS of a single file :Parameters: `fname`: string @@ -106,99 +90,18 @@ def makecorr(fname, allowed_corr): utils.write_archive(hdr) ext_wcs.readModel(hdr) for c in allowed_corr: - corr_klass = corrections.__getattribute__(c) - kw2update = corr_klass.updateWCS(ext_wcs, ref_wcs) - for kw in kw2update: - hdr.update(kw, kw2update[kw]) - + if c != 'DGEOCorr': + corr_klass = corrections.__getattribute__(c) + kw2update = corr_klass.updateWCS(ext_wcs, ref_wcs) + for kw in kw2update: + hdr.update(kw, kw2update[kw]) + + if 'DGEOCorr' in allowed_corr: + kw2update = dgeo.DGEOCorr.updateWCS(f) + for kw in kw2update: + f[1].header.update(kw, kw2update[kw]) - - #always do dgeo correction - if applyDgeoCorr(fname): - dgeo.DGEO(f) f.close() - - -def setCorrections(instrument, vacorr=True, tddcorr=True): - """ - Purpose - ======= - Creates a list of corrections to be applied to a file. - based on user input paramters and allowed corrections - for the instrument, which are defined in mappings.py. - """ - acorr = allowed_corrections[instrument] - if 'VACorr' in acorr and not vacorr: acorr.remove('VACorr') - if 'TDDCorr' in acorr and not tddcorr: acorr.remove('TDDCorr') - if 'DGEOCorr' in acorr and not dgeocorr: acorr.remove('DGEOCorr') - - return acorr - - - -def applyDgeoCorr(fname): - """ - Purpose - ======= - Adds dgeo extensions to files based on the DGEOFILE keyword in the primary - header. This is a default correction and will always run in the pipeline. - The file used to generate the extensions is - recorded in the DGEOFILE keyword in each science extension. - If 'DGEOFILE' in the primary header is different from 'DGEOFILE' in the - extension header and the file exists on disk and is a 'new type' dgeofile, - then the dgeo extensions will be updated. - """ - applyDGEOCorr = True - try: - # get DGEOFILE kw from primary header - fdgeo0 = pyfits.getval(fname, 'DGEOFILE') - fdgeo0 = fileutil.osfn(fdgeo0) - if not fileutil.findFile(fdgeo0): - print 'Kw DGEOFILE exists in primary header but file %s not found\n' % fdgeo0 - print 'DGEO correction will not be applied\n' - applyDGEOCorr = False - return applyDGEOCorr - try: - # get DGEOFILE kw from first extension header - fdgeo1 = pyfits.getval(fname, 'DGEOFILE', ext=1) - fdgeo1 = fileutil.osfn(fdgeo1) - if fdgeo1 and fileutil.findFile(fdgeo1): - if fdgeo0 != fdgeo1: - applyDGEOCorr = True - else: - applyDGEOCorr = False - else: - # dgeo file defined in first extension may not be found - # but if a valid kw exists in the primary header, dgeo should be applied. - applyDGEOCorr = True - except KeyError: - # the case of DGEOFILE kw present in primary header but missing - # in first extension header - applyDGEOCorr = True - except KeyError: - - print 'DGEOFILE keyword not found in primary header' - applyDGEOCorr = False - - if isOldStyleDGEO(fname, fdgeo0): - applyDGEOCorr = False - - return applyDGEOCorr - -def isOldStyleDGEO(fname, dgname): - # checks if the file defined in a DGEOFILE kw is a full size - # (old style) image - - sci_naxis1 = pyfits.getval(fname, 'NAXIS1', ext=1) - sci_naxis2 = pyfits.getval(fname, 'NAXIS2', ext=1) - dg_naxis1 = pyfits.getval(dgname, 'NAXIS1', ext=1) - dg_naxis2 = pyfits.getval(dgname, 'NAXIS2', ext=1) - if sci_naxis1 <= dg_naxis1 or sci_naxis2 <= dg_naxis2: - print 'Only full size (old style) XY file was found.' - print 'DGEO correction will not be applied.\n' - return True - else: - return False def getNrefchip(fobj): """ diff --git a/updatewcs/apply_corrections.py b/updatewcs/apply_corrections.py new file mode 100644 index 0000000..6dc7cdf --- /dev/null +++ b/updatewcs/apply_corrections.py @@ -0,0 +1,124 @@ +import os +import pyfits +from hstwcs.mappings import allowed_corrections +import time +from pytools import fileutil +import os.path +#Note: The order of corrections is important + +__docformat__ = 'restructuredtext' + +def setCorrections(fname, vacorr=True, tddcorr=True, dgeocorr=True): + """ + Purpose + ======= + Creates a list of corrections to be applied to a file. + based on user input paramters and allowed corrections + for the instrument, which are defined in mappings.py. + """ + instrument = pyfits.getval(fname, 'INSTRUME') + tddcorr = applyTDDCorr(fname, tddcorr) + dgeocorr = applyDgeoCorr(fname, dgeocorr) + acorr = allowed_corrections[instrument] + if 'VACorr' in acorr and not vacorr: acorr.remove('VACorr') + if 'TDDCorr' in acorr and not tddcorr: acorr.remove('TDDCorr') + if 'DGEOCorr' in acorr and not dgeocorr: acorr.remove('DGEOCorr') + + return acorr + +def applyTDDCorr(fname, utddcorr): + """ + The default value of tddcorr for all ACS images is True. + This correction will be performed if all conditions below are True: + - the user did not turn it off on the command line + - the detector is WFC + - the idc table specified in the primary header is available. + """ + + instrument = pyfits.getval(fname, 'INSTRUME') + try: + detector = pyfits.getval(fname, 'DETECTOR') + except KeyError: + detector = None + + if instrument == 'ACS' and detector == 'WFC' and utddcorr == True: + tddcorr = True + try: + idctab = pyfits.getval(fname, 'IDCTAB') + except KeyError: + tddcorr = False + #print "***IDCTAB keyword not found - not applying TDD correction***\n" + if os.path.exists(fileutil.osfn(idctab)): + tddcorr = True + else: + tddcorr = False + #print "***IDCTAB file not found - not applying TDD correction***\n" + else: + tddcorr = False + + return tddcorr + +def applyDgeoCorr(fname, udgeocorr): + """ + Purpose + ======= + Adds dgeo extensions to files based on the DGEOFILE keyword in the primary + header. This is a default correction and will always run in the pipeline. + The file used to generate the extensions is + recorded in the DGEOFILE keyword in each science extension. + If 'DGEOFILE' in the primary header is different from 'DGEOFILE' in the + extension header and the file exists on disk and is a 'new type' dgeofile, + then the dgeo extensions will be updated. + """ + applyDGEOCorr = True + try: + # get DGEOFILE kw from primary header + fdgeo0 = pyfits.getval(fname, 'DGEOFILE') + fdgeo0 = fileutil.osfn(fdgeo0) + if not fileutil.findFile(fdgeo0): + print 'Kw DGEOFILE exists in primary header but file %s not found\n' % fdgeo0 + print 'DGEO correction will not be applied\n' + applyDGEOCorr = False + return applyDGEOCorr + try: + # get DGEOFILE kw from first extension header + fdgeo1 = pyfits.getval(fname, 'DGEOFILE', ext=1) + fdgeo1 = fileutil.osfn(fdgeo1) + if fdgeo1 and fileutil.findFile(fdgeo1): + if fdgeo0 != fdgeo1: + applyDGEOCorr = True + else: + applyDGEOCorr = False + else: + # dgeo file defined in first extension may not be found + # but if a valid kw exists in the primary header, dgeo should be applied. + applyDGEOCorr = True + except KeyError: + # the case of DGEOFILE kw present in primary header but missing + # in first extension header + applyDGEOCorr = True + except KeyError: + + print 'DGEOFILE keyword not found in primary header' + applyDGEOCorr = False + + if isOldStyleDGEO(fname, fdgeo0): + applyDGEOCorr = False + print 'udgeocorr', udgeocorr, applyDGEOCorr + return (applyDGEOCorr and udgeocorr) + +def isOldStyleDGEO(fname, dgname): + # checks if the file defined in a DGEOFILE kw is a full size + # (old style) image + + sci_naxis1 = pyfits.getval(fname, 'NAXIS1', ext=1) + sci_naxis2 = pyfits.getval(fname, 'NAXIS2', ext=1) + dg_naxis1 = pyfits.getval(dgname, 'NAXIS1', ext=1) + dg_naxis2 = pyfits.getval(dgname, 'NAXIS2', ext=1) + if sci_naxis1 <= dg_naxis1 or sci_naxis2 <= dg_naxis2: + print 'Only full size (old style) XY file was found.' + print 'DGEO correction will not be applied.\n' + return True + else: + return False + diff --git a/updatewcs/corrections.py b/updatewcs/corrections.py index dc4eeff..4699b1b 100644 --- a/updatewcs/corrections.py +++ b/updatewcs/corrections.py @@ -6,7 +6,7 @@ from hstwcs.utils import diff_angles import makewcs, dgeo MakeWCS = makewcs.MakeWCS -DGEO = dgeo.DGEO +DGEOCorr = dgeo.DGEOCorr class TDDCorr(object): """ @@ -32,9 +32,11 @@ class TDDCorr(object): - sets TDDCORR to COMPLETE """ - alpha, beta = cls.set_alpha_beta(ext_wcs) + alpha, beta = cls.compute_alpha_beta(ext_wcs) cls.apply_tdd2idc(ext_wcs, alpha, beta) - newkw = {'TDDALPHA': alpha, 'TDDBETA':beta, 'TDDCORR': 'COMPLETE'} + ext_wcs.idcmodel.refpix['TDDALPHA'] = alpha + ext_wcs.idcmodel.refpix['TDDBETA'] = beta + newkw = {'TDDALPHA': alpha, 'TDDBETA':beta} return newkw updateWCS = classmethod(updateWCS) @@ -45,11 +47,8 @@ class TDDCorr(object): This should be always the first correction. """ theta_v2v3 = 2.234529 - idctheta = theta_v2v3 - idcrad = fileutil.DEGTORAD(idctheta) - mrotp = fileutil.buildRotMatrix(idctheta) - mrotn = fileutil.buildRotMatrix(-idctheta) - abmat = numpy.array([[beta,alpha],[alpha,beta]]) + mrotp = fileutil.buildRotMatrix(theta_v2v3) + mrotn = fileutil.buildRotMatrix(-theta_v2v3) tdd_mat = numpy.array([[1+(beta/2048.), alpha/2048.],[alpha/2048.,1-(beta/2048.)]],numpy.float64) abmat1 = numpy.dot(tdd_mat, mrotn) abmat2 = numpy.dot(mrotp,abmat1) @@ -59,11 +58,11 @@ class TDDCorr(object): ext_wcs.idcmodel.cy = icxy[1] ext_wcs.idcmodel.cx.shape = xshape ext_wcs.idcmodel.cy.shape = yshape - ext_wcs.wcs.set() + #ext_wcs.wcs.set() apply_tdd2idc = classmethod(apply_tdd2idc) - def set_alpha_beta(cls, ext_wcs): + def compute_alpha_beta(cls, ext_wcs): """ Compute the time dependent distortion skew terms default date of 2004.5 = 2004-7-1 @@ -77,7 +76,7 @@ class TDDCorr(object): return alpha, beta - set_alpha_beta = classmethod(set_alpha_beta) + compute_alpha_beta = classmethod(compute_alpha_beta) class VACorr(object): @@ -92,10 +91,8 @@ class VACorr(object): def updateWCS(cls, ext_wcs, ref_wcs): if ext_wcs.vafactor != 1: ext_wcs.wcs.cd = ext_wcs.wcs.cd * ext_wcs.vafactor - crval1 = ext_wcs.wcs.crval[0] - crval2 = ext_wcs.wcs.crval[1] - ext_wcs.wcs.crval[0] = ext_wcs.wcs.crval[0] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[0], crval1) - ext_wcs.wcs.crval[1] = ext_wcs.wcs.crval[1] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[1], crval2) + ext_wcs.wcs.crval[0] = ref_wcs.wcs.crval[0] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[0], ref_wcs.wcs.crval[0]) + ext_wcs.wcs.crval[1] = ref_wcs.wcs.crval[1] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[1], ref_wcs.wcs.crval[1]) kw2update={'CD1_1': ext_wcs.wcs.cd[0,0], 'CD1_2':ext_wcs.wcs.cd[0,1], 'CD2_1':ext_wcs.wcs.cd[1,0], 'CD2_2':ext_wcs.wcs.cd[1,1], diff --git a/updatewcs/dgeo.py b/updatewcs/dgeo.py index 32e1a38..d15a485 100644 --- a/updatewcs/dgeo.py +++ b/updatewcs/dgeo.py @@ -2,7 +2,7 @@ import pyfits from pytools import fileutil from hstwcs.mappings import dgeo_vals -class DGEO(object): +class DGEOCorr(object): """ Purpose ======= @@ -18,8 +18,11 @@ class DGEO(object): - Add a keyword 'DGEOFILE' to the science extension header, whose value is the reference file used to create the WCSVARR extension + If WCSDVARR extensions exist, subsequent updates will overwrite them. + If not, they will be added to the file object. + """ - def __init__(self, fobj): + def updateWCS(cls, fobj): """ :Parameters: `fobj`: pyfits object @@ -27,11 +30,16 @@ class DGEO(object): """ assert isinstance(fobj, pyfits.NP_pyfits.HDUList) - self.fobj = fobj - self.applyDgeoCorr() - + #self.fobj = fobj + cls.applyDgeoCorr(fobj) + dgfile = fobj[0].header['DGEOFILE'] - def applyDgeoCorr(self): + new_kw = {'DGEOFIEL': dgfile} + return new_kw + + updateWCS = classmethod(updateWCS) + + def applyDgeoCorr(cls,fobj): """ For each science extension in a pyfits file object: - create a WCSDVARR extension @@ -40,9 +48,10 @@ class DGEO(object): """ dgeover = 0 - dgfile = fileutil.osfn(self.fobj[0].header['DGEOFILE']) - wcsdvarr_ind = self.getwcsindex() - for ext in self.fobj: + dgfile = fileutil.osfn(fobj[0].header['DGEOFILE']) + instrument = fobj[0].header.get('INSTRUME', None) + wcsdvarr_ind = cls.getWCSIndex(fobj) + for ext in fobj: try: extname = ext.header['EXTNAME'].lower() except KeyError: @@ -52,33 +61,34 @@ class DGEO(object): for ename in ['DX', 'DY']: dgeover +=1 - self.addSciExtKw(ext.header, extver=dgeover, ename=ename) - hdu = self.createDgeoHDU(dgeofile=dgfile, dgeover=dgeover,ename=ename, extver=extversion) + cls.addSciExtKw(ext.header, extver=dgeover, ename=ename) + hdu = cls.createDgeoHDU(instrument, dgeofile=dgfile, dgeover=dgeover,ename=ename, extver=extversion) if wcsdvarr_ind: - self.fobj[wcsdvarr_ind[dgeover]] = hdu + fobj[wcsdvarr_ind[dgeover]] = hdu else: - self.fobj.append(hdu) - self.updateDGEOkw(ext.header) + fobj.append(hdu) - def getwcsindex(self): + applyDgeoCorr = classmethod(applyDgeoCorr) + + def getWCSIndex(cls, fobj): """ - Returns the index of a WCSDVARR extension in a pyfits file object if it exists. - If it exists subsequent updates will overwrite it. If not, it will be - added to the file object. + Returns a mapping of WCSDVARR 'EXTVER' and pyfits.HDUList indeces. + """ wcsd = {} - for e in range(len(self.fobj)): + for e in range(len(fobj)): try: - ename = self.fobj[e].header['EXTNAME'] + ename = fobj[e].header['EXTNAME'] except KeyError: continue if ename == 'WCSDVARR': - wcsd[self.fobj[e].header['EXTVER']] = e + wcsd[fobj[e].header['EXTVER']] = e + return wcsd - return self.fobj.index_of(('WCSDVARR', dgeover)) + getWCSIndex = classmethod(getWCSIndex) - def addSciExtKw(self, hdr, extver=None, ename=None): + def addSciExtKw(cls, hdr, extver=None, ename=None): """ Adds kw to sci extension to define dgeo correction extensions kw to be added to sci ext: @@ -114,36 +124,37 @@ class DGEO(object): for key in keys: hdr.update(key=key, value=values[key], comment=comments[key], before='HISTORY') - dgfile = self.fobj[0].header['DGEOFILE'] - hdr.update(key='DGEOFILE', value=dgfile, comment='DGEOFILE used to create WCSDVARR extensions', before='HISTORY') - - - def getDgeoData(self, dgfile=None, ename=None, extver=1): + addSciExtKw = classmethod(addSciExtKw) + + def getDgeoData(cls, dgfile=None, ename=None, extver=1): """ Given a dgeo file name, creates an array to be used as a data array in the dgeo extension. """ return pyfits.getdata(dgfile, ext=(ename,extver)) + + getDgeoData = classmethod(getDgeoData) - def createDgeoHDU(self, dgeofile=None, dgeover=1, ename=None,extver=1): + def createDgeoHDU(cls, instrument, dgeofile=None, dgeover=1, ename=None,extver=1): """ Creates an HDU to be added to the file object. """ dgeokw = {'naxis1':None, 'naxis2':None, 'extver':dgeover, 'crpix1':None, 'crpix2':None, 'cdelt1':None, 'cdelt2':None, 'crval1':None, 'crval2':None} - hdr = self.createDgeoHdr(**dgeokw) - data = self.getDgeoData(dgfile=dgeofile, ename=ename, extver=extver) + hdr = cls.createDgeoHdr(instrument, **dgeokw) + data = cls.getDgeoData(dgfile=dgeofile, ename=ename, extver=extver) hdu=pyfits.ImageHDU(header=hdr, data=data) return hdu + createDgeoHDU = classmethod(createDgeoHDU) - def createDgeoHdr(self, **kw): + def createDgeoHdr(cls, instr, **kw): """ Creates a header for the dgeo extension based on dgeo file and sci extension header. **kw = {'naxis1':None, 'naxis2':None, 'extver':None, 'crpix1':None, 'crpix2':None, 'cdelt1':None, 'cdelt2':None, 'crval1':None, 'crval2':None} """ - instr = self.fobj[0].header['INSTRUME'] + #instr = self.fobj[0].header['INSTRUME'] instr_vals = dgeo_vals[instr] naxis1 = kw['naxis1'] or instr_vals['naxis1'] naxis2 = kw['naxis2'] or instr_vals['naxis2'] @@ -199,8 +210,6 @@ class DGEO(object): hdr = pyfits.Header(cards=cdl) return hdr - - def updateDGEOkw(self, hdr): - dgfile = self.fobj[0].header['DGEOFILE'] - hdr.update(key='DGEOFILE', value=dgfile, comment='DGEOFILE used to create WCSDVARR extensions', before='HISTORY') -
\ No newline at end of file + + createDgeoHdr = classmethod(createDgeoHdr) + diff --git a/updatewcs/makewcs.py b/updatewcs/makewcs.py index 3df2fd4..1f3e640 100644 --- a/updatewcs/makewcs.py +++ b/updatewcs/makewcs.py @@ -163,7 +163,6 @@ class MakeWCS(object): cd22 = parity[1][1] * cos(pv*pi/180.0)*R_scale rcd = numpy.array([[cd11, cd12], [cd21, cd22]]) - print 'cd shape', rcd.shape, ref_wcs.wcs.cd.shape ref_wcs.wcs.cd = rcd ref_wcs.wcs.set() |