diff options
| -rw-r--r-- | lib/__init__.py | 47 | ||||
| -rw-r--r-- | lib/utils.py | 83 | ||||
| -rw-r--r-- | updatewcs/__init__.py | 45 | ||||
| -rw-r--r-- | updatewcs/corrections.py | 167 | ||||
| -rw-r--r-- | updatewcs/makewcs.py | 189 | ||||
| -rw-r--r-- | wcsutil/__init__.py | 115 | ||||
| -rw-r--r-- | wcsutil/instruments.py | 2 |
7 files changed, 299 insertions, 349 deletions
diff --git a/lib/__init__.py b/lib/__init__.py index 08e421b..c9dca75 100644 --- a/lib/__init__.py +++ b/lib/__init__.py @@ -1,47 +1,8 @@ -import wcsutil -from wcsutil import HSTWCS -from pytools import fileutil, parseinput -import pyfits - -#to avoid relative imports +#import all needed modules here to avoid relative imports import mappings -from mappings import basic_wcs +import utils import distortion import pywcs -def restoreWCS(fnames): - """ - Given a list of fits file names restore the original basic WCS kw - and write out the files. This overwrites the original files. - """ - 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 on multiextension fits files." - return - else: - fobj = pyfits.open(f, mode='update') - hdr0 = fobj[0].header - for ext in fobj: - try: - extname = ext.header['EXTNAME'].lower() - except KeyError: - continue - if extname in ['sci', 'err', 'sdq']: - hdr = ext.header - owcs = HSTWCS(hdr0, hdr) - #Changed restore to update the attributes ??? - # this may need to be changed - backup = owcs.restore() - if not backup: - print '\No archived keywords found.\n' - continue - else: - for kw in basic_wcs: - nkw = ('O'+kw)[:7] - if backup.has_key(kw): - hdr.update(kw, hdr[nkw]) - fobj.close() - - + +
\ No newline at end of file diff --git a/lib/utils.py b/lib/utils.py new file mode 100644 index 0000000..bec12f7 --- /dev/null +++ b/lib/utils.py @@ -0,0 +1,83 @@ +from pytools import parseinput, fileutil +import pyfits +from mappings import basic_wcs + +def restoreWCS(fnames): + """ + Given a list of fits file names restore the original basic WCS kw + and write out the files. This overwrites the original files. + + Affected keywords: + 'CD1_1', 'CD1_2', 'CD2_1', 'CD2_2', 'CRVAL1','CRVAL2','CTYPE1', 'CTYPE2', + 'CRPIX1', 'CRPIX2', 'CTYPE1', 'CTYPE2', 'ORIENTAT' + """ + 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: + fobj = pyfits.open(f, mode='update') + for ext in fobj: + try: + extname = ext.header['EXTNAME'].lower() + except KeyError: + continue + if extname in ['sci', 'err', 'sdq']: + hdr = ext.header + backup = get_archive(hdr) + if not backup: + #print 'No archived keywords found.\n' + continue + else: + for kw in basic_wcs: + nkw = ('O'+kw)[:7] + if backup.has_key(kw): + hdr.update(kw, hdr[nkw]) + fobj.close() + +def get_archive(header): + """ + Returns a dictionary with the archived basic WCS keywords. + """ + + backup = {} + for k in basic_wcs: + try: + nkw = ('O'+k)[:7] + backup[k] = header[nkw] + except KeyError: + pass + return backup + +def write_archive(header): + """ + Archives original WCS kw before recalculating them. + """ + backup_kw = get_archive(header) + if backup_kw != {}: + print "Archive already exists\n." + return + else: + cmt = 'archived value' + for kw in basic_wcs: + nkw = 'O'+kw + header.update(nkw[:7], header[kw], comment=cmt) + + +def diff_angles(a,b): + """ + Perform angle subtraction a-b taking into account + small-angle differences across 360degree line. + """ + + diff = a - b + + if diff > 180.0: + diff -= 360.0 + + if diff < -180.0: + diff += 360.0 + + return diff
\ No newline at end of file diff --git a/updatewcs/__init__.py b/updatewcs/__init__.py index db89d13..107a9f5 100644 --- a/updatewcs/__init__.py +++ b/updatewcs/__init__.py @@ -4,12 +4,13 @@ import pyfits 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 -#NB! the order of corrections matters +#Note: The order of corrections is important __docformat__ = 'restructuredtext' @@ -52,19 +53,30 @@ def updatewcs(input, vacorr=True, tddcorr=True, checkfiles=True): print 'No valid input, quitting ...\n' return for f in files: - instr = pyfits.getval(f, 'INSTRUME') + 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: - acorr = setCorrections(instr,vacorr=vacorr, tddcorr=tddcorr) + #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 - + #restore the original WCS keywords + utils.restoreWCS(f) makecorr(f, acorr) return files -def makecorr(fname, acorr): +def makecorr(fname, allowed_corr): """ Purpose ======= @@ -78,23 +90,26 @@ def makecorr(fname, acorr): """ f = pyfits.open(fname, mode='update') + #Determine the reference chip and make a copy of its restored header. nrefchip, nrefext = getNrefchip(f) primhdr = f[0].header - + ref_hdr = f[nrefext].header.copy() + utils.write_archive(ref_hdr) for extn in f: - refwcs = HSTWCS(primhdr, f[nrefext].header) - refwcs.archive_kw() - refwcs.readModel() + # Perhaps all ext headers should be corrected (to be consistent) if extn.header.has_key('extname') and extn.header['extname'].lower() == 'sci': + refwcs = HSTWCS(primhdr, ref_hdr, fobj=f) + refwcs.readModel(ref_hdr) hdr = extn.header - owcs = HSTWCS(primhdr, hdr) - owcs.archive_kw() - owcs.readModel() - for c in acorr: - owcs.__setattr__('DO'+c, 'PERFORM') + owcs = HSTWCS(primhdr, hdr, fobj=f) + utils.write_archive(hdr) + owcs.readModel(hdr) + for c in allowed_corr: corr_klass = corrections.__getattribute__(c) - corr_klass(owcs, refwcs) + kw2update = corr_klass.updateWCS(owcs, refwcs) + for kw in kw2update: + hdr.update(kw, kw2update[kw]) diff --git a/updatewcs/corrections.py b/updatewcs/corrections.py index e94be64..ac49833 100644 --- a/updatewcs/corrections.py +++ b/updatewcs/corrections.py @@ -2,51 +2,44 @@ import datetime import numpy from numpy import linalg from pytools import fileutil +from hstwcs.utils import diff_angles +import makewcs, dgeo -from makewcs import MakeWCS -from dgeo import DGEO +MakeWCS = makewcs.MakeWCS +DGEO = dgeo.DGEO class TDDCorr(object): """ Purpose ======= Apply time dependent distortion correction to SIP coefficients and basic - WCS keywords. Atthi stime it is applicable only to ACS/WFC data. + WCS keywords. It is applicable only to ACS/WFC data. Ref: Jay Anderson, ACS ISR 2007-08, Variation of the Distortion Solution of the WFC + :Parameters: + `owcs`: HSTWCS object + An extension HSTWCS object to be modified + `refwcs`: HSTWCS object + A reference HSTWCS object """ - def __init__(self, owcs=None, refwcs=None): - """ - :Parameters: - `owcs`: HSTWCS object - An extension HSTWCS object to be modified - `refwcs`: HSTWCS object - A reference HSTWCS object - """ - self.wcsobj = owcs - if self.wcsobj.DOTDDCorr == 'PERFORM': - self.updateWCS() - self.wcsobj.DOTDDCorr = 'COMPLETE' - else: - pass - - def updateWCS(self): + + def updateWCS(cls, wcs, refwcs): """ - Calculates alpha and beta for ACS/WFC data. - Writes 2 new kw to the extension header: TDDALPHA and TDDBETA + - sets TDDCORR to COMPLETE """ - if self.wcsobj.instrument == 'ACS' and self.wcsobj.detector == 'WFC': - newkw = ['TDDALPHA', 'TDDBETA'] - self.set_alpha_beta() - self.applyTDD() - self.updatehdr(newkeywords=newkw) - - else: - pass + + alpha, beta = cls.set_alpha_beta(wcs) + cls.apply_tdd2idc(wcs, alpha, beta) + newkw = {'TDDALPHA': alpha, 'TDDBETA':beta, 'TDDCORR': 'COMPLETE'} + + return newkw + updateWCS = classmethod(updateWCS) - def applyTDD(self): + def apply_tdd2idc(cls,wcs, alpha, beta): """ Applies TDD to the idctab coefficients of a ACS/WFC observation. This should be always the first correction. @@ -56,36 +49,36 @@ class TDDCorr(object): idcrad = fileutil.DEGTORAD(idctheta) mrotp = fileutil.buildRotMatrix(idctheta) mrotn = fileutil.buildRotMatrix(-idctheta) - abmat = numpy.array([[self.TDDBETA,self.TDDALPHA],[self.TDDALPHA,self.TDDBETA]]) - tdd_mat = numpy.array([[1+(self.TDDBETA/2048.), self.TDDALPHA/2048.],[self.TDDALPHA/2048.,1-(self.TDDBETA/2048.)]],numpy.float64) + abmat = numpy.array([[beta,alpha],[alpha,beta]]) + 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) - xshape, yshape = self.wcsobj.idcmodel.cx.shape, self.wcsobj.idcmodel.cy.shape - icxy = numpy.dot(abmat2,[self.wcsobj.idcmodel.cx.ravel(),self.wcsobj.idcmodel.cy.ravel()]) - self.wcsobj.idcmodel.cx = icxy[0] - self.wcsobj.idcmodel.cy = icxy[1] - self.wcsobj.idcmodel.cx.shape = xshape - self.wcsobj.idcmodel.cy.shape = yshape + xshape, yshape = wcs.idcmodel.cx.shape, wcs.idcmodel.cy.shape + icxy = numpy.dot(abmat2,[wcs.idcmodel.cx.ravel(),wcs.idcmodel.cy.ravel()]) + wcs.idcmodel.cx = icxy[0] + wcs.idcmodel.cy = icxy[1] + wcs.idcmodel.cx.shape = xshape + wcs.idcmodel.cy.shape = yshape + wcs.set() + + apply_tdd2idc = classmethod(apply_tdd2idc) - def set_alpha_beta(self): + def set_alpha_beta(cls, wcs): """ Compute the time dependent distortion skew terms default date of 2004.5 = 2004-7-1 """ dday = 2004.5 - year,month,day = self.wcsobj.date_obs.split('-') + year,month,day = wcs.date_obs.split('-') rdate = datetime.datetime(int(year),int(month),int(day)) rday = float(rdate.strftime("%j"))/365.25 + rdate.year alpha = 0.095 + 0.090*(rday-dday)/2.5 beta = -0.029 - 0.030*(rday-dday)/2.5 - self.TDDALPHA = alpha - self.TDDBETA = beta - - def updatehdr(self, newkeywords=None): - - for kw in newkeywords: - self.wcsobj.hdr.update(kw, self.__getattribute__(kw)) + return alpha, beta + + set_alpha_beta = classmethod(set_alpha_beta) + class VACorr(object): """ @@ -93,60 +86,44 @@ class VACorr(object): ======= Apply velocity aberation correction to WCS keywords. Modifies the CD matrix and CRVAL1/2 + """ - def __init__(self, owcs=None, refwcs=None): - self.wcsobj = owcs - self.refwcs = refwcs - if self.wcsobj.DOVACorr == 'PERFORM': - self.updateWCS() - self.wcsobj.DOVACorr == 'COMPLETE' - else: - pass - - def updateWCS(self): - if self.wcsobj.vafactor != 1: - self.wcsobj.cd = self.wcsobj.cd * self.wcsobj.vafactor + + def updateWCS(cls, wcs, refwcs): + if wcs.vafactor != 1: + wcs.cd = wcs.cd * wcs.vafactor + crval1 = wcs.crval[0] + crval2 = wcs.crval[1] + wcs.crval[0] = wcs.crval[0] + wcs.vafactor*diff_angles(wcs.crval[0], crval1) + wcs.crval[1] = wcs.crval[1] + wcs.vafactor*diff_angles(wcs.crval[1], crval2) - #self.wcsobj.crval[0] = self.wcsobj.ra_targ + self.wcsobj.vafactor*(self.wcsobj.crval[0] - self.wcsobj.ra_targ) - #self.wcsobj.crval[1] = self.wcsobj.dec_targ + self.wcsobj.vafactor*(self.wcsobj.crval[1] - self.wcsobj.dec_targ) - ref_backup = self.refwcs.restore() - crval1 = ref_backup['CRVAL1'] - crval2 = ref_backup['CRVAL2'] - self.wcsobj.crval[0] = crval1 + self.wcsobj.vafactor*diff_angles(self.wcsobj.crval[0], crval1) - self.wcsobj.crval[1] = crval2 + self.wcsobj.vafactor*diff_angles(self.wcsobj.crval[1], crval2) - - kw2update={'CD1_1': self.wcsobj.cd[0,0], 'CD1_2':self.wcsobj.cd[0,1], - 'CD2_1':self.wcsobj.cd[1,0], 'CD2_2':self.wcsobj.cd[1,1], - 'CRVAL1':self.wcsobj.crval[0], 'CRVAL2':self.wcsobj.crval[1]} - self.updatehdr(newkeywords=kw2update) + kw2update={'CD1_1': wcs.cd[0,0], 'CD1_2':wcs.cd[0,1], + 'CD2_1':wcs.cd[1,0], 'CD2_2':wcs.cd[1,1], + 'CRVAL1':wcs.crval[0], 'CRVAL2':wcs.crval[1]} + wcs.set() else: pass - - def updatehdr(self, newkeywords=None): - for kw in newkeywords: - self.wcsobj.hdr.update(kw, newkeywords[kw]) + return kw2update + + updateWCS = classmethod(updateWCS) + class CompSIP(object): """ Purpose ======= Compute SIP coefficients from idc table coefficients. - """ - def __init__(self, owcs, refwcs): - self.wcsobj = owcs - self.updateWCS() - self.DOCOMPSIP = 'COMPLETE' - - def updateWCS(self): + + def updateWCS(cls, wcs, refwcs): kw2update = {} - order = self.wcsobj.idcmodel.norder + order = wcs.idcmodel.norder kw2update['A_ORDER'] = order kw2update['B_ORDER'] = order - pscale = self.wcsobj.idcmodel.refpix['PSCALE'] + pscale = wcs.idcmodel.refpix['PSCALE'] - cx = self.wcsobj.idcmodel.cx - cy = self.wcsobj.idcmodel.cy + cx = wcs.idcmodel.cx + cy = wcs.idcmodel.cy matr = numpy.array([[cx[1,1],cx[1,0]], [cy[1,1],cy[1,0]]], dtype=numpy.float) imatr = linalg.inv(matr) @@ -163,26 +140,8 @@ class CompSIP(object): Bkey="B_%d_%d" % (m,n-m) kw2update[Akey] = sipval[0,0] kw2update[Bkey] = sipval[1,0] - self.updatehdr(newkw=kw2update) - - def updatehdr(self, newkw=None): - if not newkw: return - for kw in newkw.keys(): - self.wcsobj.hdr.update(kw, newkw[kw]) - - -def diff_angles(a,b): - """ - Perform angle subtraction a-b taking into account - small-angle differences across 360degree line. - """ + return kw2update - diff = a - b + updateWCS = classmethod(updateWCS) - if diff > 180.0: - diff -= 360.0 - if diff < -180.0: - diff += 360.0 - - return diff diff --git a/updatewcs/makewcs.py b/updatewcs/makewcs.py index 451da13..ddb4d06 100644 --- a/updatewcs/makewcs.py +++ b/updatewcs/makewcs.py @@ -3,6 +3,7 @@ from hstwcs.mappings import DEGTORAD, RADTODEG import numpy from numpy import math from math import sin, sqrt, pow, cos, asin, atan2,pi +from hstwcs import utils class MakeWCS(object): """ @@ -17,74 +18,73 @@ class MakeWCS(object): - update extension header """ - - def __init__(self, owcs, refwcs): - self.wcsobj = owcs - self.refwcs = refwcs - self.updateWCS() - self.DOMakeWCS = 'COMPLETE' - - def updateWCS(self): + + def updateWCS(cls, wcs, refwcs): """ recomputes the basic WCS kw """ - backup_wcs = self.wcsobj.restore() - ltvoff, offshift = self.getOffsets(backup_wcs) + ltvoff, offshift = cls.getOffsets(wcs) - self.uprefwcs(ltvoff, offshift) - self.upextwcs(ltvoff, offshift) - self.updatehdr() - - def upextwcs(self, loff, offsh): + cls.uprefwcs(wcs, refwcs, ltvoff, offshift) + cls.upextwcs(wcs, refwcs, ltvoff, offshift) + + kw2update = {'CD1_1': wcs.cd[0,0], + 'CD1_2': wcs.cd[0,1], + 'CD2_1': wcs.cd[1,0], + 'CD2_2': wcs.cd[1,1], + 'CRVAL1': wcs.crval[0], + 'CRVAL2': wcs.crval[1], + 'CRPIX1': wcs.crpix[0], + 'CRPIX2': wcs.crpix[1], + } + return kw2update + + updateWCS = classmethod(updateWCS) + + def upextwcs(self, wcs, refwcs, loff, offsh): """ updates an extension wcs """ ltvoffx, ltvoffy = loff[0], loff[1] offshiftx, offshifty = offsh[0], offsh[1] - backup_wcs = self.wcsobj.restore() - ltv1 = self.wcsobj.ltv1 - ltv2 = self.wcsobj.ltv2 - - + ltv1 = wcs.ltv1 + ltv2 = wcs.ltv2 if ltv1 != 0. or ltv2 != 0.: - offsetx = backup_wcs['CRPIX1'] - ltv1 - self.wcsobj.idcmodel.refpix['XREF'] - offsety = backup_wcs['CRPIX2'] - ltv2 - self.wcsobj.idcmodel.refpix['YREF'] - fx,fy = self.wcsobj.idcmodel.shift(self.wcsobj.idcmodel.cx,self.wcsobj.idcmodel.cy,offsetx,offsety) + offsetx = backup_wcs['CRPIX1'] - ltv1 - wcs.idcmodel.refpix['XREF'] + offsety = backup_wcs['CRPIX2'] - ltv2 - wcs.idcmodel.refpix['YREF'] + fx,fy = wcs.idcmodel.shift(wcs.idcmodel.cx,wcs.idcmodel.cy,offsetx,offsety) else: - fx, fy = self.wcsobj.idcmodel.cx, self.wcsobj.idcmodel.cy + fx, fy = wcs.idcmodel.cx, wcs.idcmodel.cy - ridcmodel = self.refwcs.idcmodel - v2 = self.wcsobj.idcmodel.refpix['V2REF'] - v3 = self.wcsobj.idcmodel.refpix['V3REF'] - v2ref = self.refwcs.idcmodel.refpix['V2REF'] + ridcmodel = refwcs.idcmodel + v2 = wcs.idcmodel.refpix['V2REF'] + v3 = wcs.idcmodel.refpix['V3REF'] + v2ref = refwcs.idcmodel.refpix['V2REF'] - v3ref = self.refwcs.idcmodel.refpix['V3REF'] - R_scale = self.refwcs.idcmodel.refpix['PSCALE']/3600.0 + v3ref = refwcs.idcmodel.refpix['V3REF'] + R_scale = refwcs.idcmodel.refpix['PSCALE']/3600.0 off = sqrt((v2-v2ref)**2 + (v3-v3ref)**2)/(R_scale*3600.0) - - # Here we must include the PARITY if v3 == v3ref: theta=0.0 else: - theta = atan2(self.wcsobj.parity[0][0]*(v2-v2ref),self.wcsobj.parity[1][1]*(v3-v3ref)) + theta = atan2(wcs.parity[0][0]*(v2-v2ref), wcs.parity[1][1]*(v3-v3ref)) - if self.refwcs.idcmodel.refpix['THETA']: theta += self.refwcs.idcmodel.refpix['THETA']*pi/180.0 + if refwcs.idcmodel.refpix['THETA']: theta += refwcs.idcmodel.refpix['THETA']*pi/180.0 dX=(off*sin(theta)) + offshiftx dY=(off*cos(theta)) + offshifty px = numpy.array([[dX,dY]]) - newcrval = self.refwcs.pixel2world_fits(px)[0] - newcrpix = numpy.array([self.wcsobj.idcmodel.refpix['XREF'] + ltvoffx, - self.wcsobj.idcmodel.refpix['YREF'] + ltvoffy]) - - self.wcsobj.crval = newcrval - self.wcsobj.crpix = newcrpix + newcrval = refwcs.pixel2world_fits(px)[0] + newcrpix = numpy.array([wcs.idcmodel.refpix['XREF'] + ltvoffx, + wcs.idcmodel.refpix['YREF'] + ltvoffy]) + wcs.crval = newcrval + wcs.crpix = newcrpix # Account for subarray offset # Angle of chip relative to chip - if self.wcsobj.idcmodel.refpix['THETA']: - dtheta = self.wcsobj.idcmodel.refpix['THETA'] - self.refwcs.idcmodel.refpix['THETA'] + if wcs.idcmodel.refpix['THETA']: + dtheta = wcs.idcmodel.refpix['THETA'] - refwcs.idcmodel.refpix['THETA'] else: dtheta = 0.0 @@ -96,8 +96,7 @@ class MakeWCS(object): delYX=fy[1,1]/R_scale/3600. delXY=fx[1,0]/R_scale/3600. delYY=fy[1,0]/R_scale/3600. - - + # Convert to radians rr=dtheta*pi/180.0 @@ -110,93 +109,74 @@ class MakeWCS(object): px = numpy.array([[dX+dXX,dY+dYX]]) # Transform to sky coordinates - - wc = self.refwcs.pixel2world_fits(px) + wc = refwcs.pixel2world_fits(px) a = wc[0,0] b = wc[0,1] px = numpy.array([[dX+dXY,dY+dYY]]) - wc = self.refwcs.pixel2world_fits(px) + wc = refwcs.pixel2world_fits(px) c = wc[0,0] d = wc[0,1] - + # Calculate the new CDs and convert to degrees - cd11 = diff_angles(a,newcrval[0])*cos(newcrval[1]*pi/180.0) - cd12 = diff_angles(c,newcrval[0])*cos(newcrval[1]*pi/180.0) - cd21 = diff_angles(b,newcrval[1]) - cd22 = diff_angles(d,newcrval[1]) - + cd11 = utils.diff_angles(a,newcrval[0])*cos(newcrval[1]*pi/180.0) + cd12 = utils.diff_angles(c,newcrval[0])*cos(newcrval[1]*pi/180.0) + cd21 = utils.diff_angles(b,newcrval[1]) + cd22 = utils.diff_angles(d,newcrval[1]) cd = numpy.array([[cd11, cd12], [cd21, cd22]]) - self.wcsobj.cd = cd + wcs.cd = cd + + upextwcs = classmethod(upextwcs) - def uprefwcs(self, loff, offsh): + def uprefwcs(self, wcs, refwcs, loff, offsh): """ Updates the reference chip """ ltvoffx, ltvoffy = loff[0], loff[1] offshift = offsh - backup_refwcs = self.refwcs.restore() - dec = backup_refwcs['CRVAL2'] - + dec = refwcs.crval[1] # Get an approximate reference position on the sky - rref = numpy.array([[self.refwcs.idcmodel.refpix['XREF']+ltvoffx, - self.refwcs.idcmodel.refpix['YREF']+ltvoffy]]) - - - crval = self.refwcs.pixel2world_fits(rref) + rref = numpy.array([[refwcs.idcmodel.refpix['XREF']+ltvoffx, + refwcs.idcmodel.refpix['YREF']+ltvoffy]]) + crval = refwcs.pixel2world_fits(rref) # Convert the PA_V3 orientation to the orientation at the aperture # This is for the reference chip only - we use this for the # reference tangent plane definition # It has the same orientation as the reference chip - pv = troll(self.wcsobj.pav3,dec,self.refwcs.idcmodel.refpix['V2REF'], - self.refwcs.idcmodel.refpix['V3REF']) + pv = troll(wcs.pav3,dec,refwcs.idcmodel.refpix['V2REF'], + refwcs.idcmodel.refpix['V3REF']) # Add the chip rotation angle - if self.refwcs.idcmodel.refpix['THETA']: - pv += self.refwcs.idcmodel.refpix['THETA'] + if refwcs.idcmodel.refpix['THETA']: + pv += refwcs.idcmodel.refpix['THETA'] # Set values for the rest of the reference WCS - self.refwcs.crval = crval[0] - self.refwcs.crpix = numpy.array([0.0,0.0])+offsh - parity = self.refwcs.parity - R_scale = self.refwcs.idcmodel.refpix['PSCALE']/3600.0 + refwcs.crval = crval[0] + refwcs.crpix = numpy.array([0.0,0.0])+offsh + parity = refwcs.parity + R_scale = refwcs.idcmodel.refpix['PSCALE']/3600.0 cd11 = parity[0][0] * cos(pv*pi/180.0)*R_scale cd12 = parity[0][0] * -sin(pv*pi/180.0)*R_scale cd21 = parity[1][1] * sin(pv*pi/180.0)*R_scale cd22 = parity[1][1] * cos(pv*pi/180.0)*R_scale rcd = numpy.array([[cd11, cd12], [cd21, cd22]]) - self.refwcs.cd = rcd - self.refwcs.set() - - + refwcs.cd = rcd + refwcs.set() - def updatehdr(self): - """ - Keywords to update: - CD1_1, CD1_2, CD2_1, CD2_2, CRPIX1, CRPIX2, CRVAL1, CRVAL2 - """ - - self.wcsobj.hdr.update('CD1_1', self.wcsobj.cd[0,0]) - self.wcsobj.hdr.update('CD1_2', self.wcsobj.cd[0,1]) - self.wcsobj.hdr.update('CD2_1', self.wcsobj.cd[1,0]) - self.wcsobj.hdr.update('CD2_2', self.wcsobj.cd[1,1]) - self.wcsobj.hdr.update('CRVAL1', self.wcsobj.crval[0]) - self.wcsobj.hdr.update('CRVAL2', self.wcsobj.crval[1]) - self.wcsobj.hdr.update('CRPIX1', self.wcsobj.crpix[0]) - self.wcsobj.hdr.update('CRPIX2', self.wcsobj.crpix[1]) - self.wcsobj.hdr.update('ORIENTAT', self.wcsobj.orientat) + uprefwcs = classmethod(uprefwcs) - def getOffsets(self, backup_wcs): - ltv1 = self.wcsobj.ltv1 - ltv2 = self.wcsobj.ltv2 + + def getOffsets(cls, wcs): + ltv1 = wcs.ltv1 + ltv2 = wcs.ltv2 - offsetx = backup_wcs['CRPIX1'] - ltv1 - self.wcsobj.idcmodel.refpix['XREF'] - offsety = backup_wcs['CRPIX2'] - ltv2 - self.wcsobj.idcmodel.refpix['YREF'] + offsetx = wcs.crpix[0] - ltv1 - wcs.idcmodel.refpix['XREF'] + offsety = wcs.crpix[1] - ltv2 - wcs.idcmodel.refpix['YREF'] - shiftx = self.wcsobj.idcmodel.refpix['XREF'] + ltv1 - shifty = self.wcsobj.idcmodel.refpix['YREF'] + ltv2 + shiftx = wcs.idcmodel.refpix['XREF'] + ltv1 + shifty = wcs.idcmodel.refpix['YREF'] + ltv2 if ltv1 != 0. or ltv2 != 0.: ltvoffx = ltv1 + offsetx ltvoffy = ltv2 + offsety @@ -212,21 +192,8 @@ class MakeWCS(object): offshift = numpy.array([offshiftx, offshifty]) return ltvoff, offshift -def diff_angles(a,b): - """ - Perform angle subtraction a-b taking into account - small-angle differences across 360degree line. - """ - - diff = a - b - - if diff > 180.0: - diff -= 360.0 - - if diff < -180.0: - diff += 360.0 + getOffsets = classmethod(getOffsets) - return diff def troll(roll, dec, v2, v3): """ Computes the roll angle at the target position based on: diff --git a/wcsutil/__init__.py b/wcsutil/__init__.py index 46ad42e..9482a48 100644 --- a/wcsutil/__init__.py +++ b/wcsutil/__init__.py @@ -16,34 +16,34 @@ class HSTWCS(SIP): """ Purpose ======= - Create a wcs object based on the instrument. + Create a WCS object based on the instrument. It has all basic WCS kw as attribbutes (set by pywcs). It also uses the primary and extension header to define instrument specific attributes needed by the correction classes. """ - def __init__(self, hdr0, hdr): + def __init__(self, hdr0, ehdr, fobj=None): """ :Parameters: `hdr0`: Pyfits Header primary header - `hdr`: Pyfits Header + `ehdr`: Pyfits Header extension header + `fobj`: PyFITS HDUList object or None + pyfits file object """ - self.hdr = hdr - self.hdr0 = hdr0 - SIP.__init__(self, self.hdr) + SIP.__init__(self, ehdr, fobj=fobj) self.setHDR0kw(hdr0) - self.detector = self.setDetector() + self.detector = self.setDetector(hdr0) self.inst_kw = ins_spec_kw - self.setInstrSpecKw() + self.setInstrSpecKw(hdr0, ehdr) self.pscale = self.setPscale() self.orientat = self.setOrient() def setHDR0kw(self, primhdr): - #sets attributes from kw defined in the primary header + # Set attributes from kw defined in the primary header. self.instrument = primhdr.get('INSTRUME', None) self.offtab = primhdr.get('OFFTAB', None) self.idctab = primhdr.get('IDCTAB', None) @@ -53,34 +53,35 @@ class HSTWCS(SIP): self.dec_targ = primhdr.get('DEC_TARG', None) - def setDetector(self): - #sets detector attribute for instuments which have more than one detector + def setDetector(self, header): + # Set detector attribute for instuments which have more than one detector if self.instrument in ['ACS', 'WFC3']: - return self.hdr0.get('DETECTOR', None) - - - def setInstrSpecKw(self): - #Based on the instrument kw creates an instalnce of an instrument WCS class - #and sets attributes from instrument specific kw + return header.get('DETECTOR', None) + else: + return None + + def setInstrSpecKw(self, prim_hdr, ext_hdr): + # Based on the instrument kw creates an instalnce of an instrument WCS class + # and sets attributes from instrument specific kw if self.instrument in inst_mappings.keys(): inst_kl = inst_mappings[self.instrument] inst_kl = instruments.__dict__[inst_kl] - insobj = inst_kl(self.hdr0, self.hdr) + insobj = inst_kl(prim_hdr, ext_hdr) for key in self.inst_kw: self.__setattr__(key, insobj.__getattribute__(key)) else: raise KeyError, "Unsupported instrument - %s" %self.instrument def setPscale(self): - #Calculates the plate scale from the cd matrix - #this may not be needed now and shoufl probably be done after all + # Calculates the plate scale from the cd matrix + # this may not be needed now and shoufl probably be done after all # corrections cd11 = self.cd[0][0] cd21 = self.cd[1][0] return N.sqrt(N.power(cd11,2)+N.power(cd21,2)) * 3600. def setOrient(self): - #same as setPscale + # Recompute ORIENTAT cd12 = self.cd[0][1] cd22 = self.cd[1][1] return RADTODEG(N.arctan2(cd12,cd22)) @@ -91,7 +92,7 @@ class HSTWCS(SIP): def verifyKw(self): """verify that all required kw have meaningful values""" - def readModel(self): + def readModel(self, header): """ Purpose ======= @@ -106,61 +107,27 @@ class HSTWCS(SIP): offtab=self.offtab, binned=self.binned) - self.updatehdr() + self.updatehdr(header) - def updatehdr(self, newkeywords=None): + def updatehdr(self, ext_hdr, newkeywords=None): #kw2add : OCX10, OCX11, OCY10, OCY11 # record the model in the header for use by pydrizzle - self.hdr.update('OCX10', self.idcmodel.cx[1,0]) - self.hdr.update('OCX11', self.idcmodel.cx[1,1]) - self.hdr.update('OCY10', self.idcmodel.cy[1,0]) - self.hdr.update('OCY11', self.idcmodel.cy[1,1]) - self.hdr.update('IDCSCALE', self.idcmodel.refpix['PSCALE']) - self.hdr.update('IDCTHETA', self.idcmodel.refpix['THETA']) - self.hdr.update('IDCXREF', self.idcmodel.refpix['XREF']) - self.hdr.update('IDCYREF', self.idcmodel.refpix['YREF']) - self.hdr.update('IDCV2REF', self.idcmodel.refpix['V2REF']) - self.hdr.update('IDCV3REF', self.idcmodel.refpix['V3REF']) - self.hdr.update('IDCXSIZE', self.idcmodel.refpix['XSIZE']) - self.hdr.update('IDCYSIZE', self.idcmodel.refpix['YSIZE']) - self.hdr.update('IDCXDELT', self.idcmodel.refpix['XDELTA']) - self.hdr.update('IDCYDELT', self.idcmodel.refpix['YDELTA']) - self.hdr.update('CENTERED', self.idcmodel.refpix['centered']) + 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']) + ext_hdr.update('IDCXSIZE', self.idcmodel.refpix['XSIZE']) + ext_hdr.update('IDCYSIZE', self.idcmodel.refpix['YSIZE']) + ext_hdr.update('IDCXDELT', self.idcmodel.refpix['XDELTA']) + ext_hdr.update('IDCYDELT', self.idcmodel.refpix['YDELTA']) + ext_hdr.update('CENTERED', self.idcmodel.refpix['centered']) - def restore(self): - """ - restores basic wcs keywords from archive - this should be modified to always return a populated - dictionary, although the values may be None. - """ - - backup = {} - for k in basic_wcs: - try: - nkw = ('O'+k)[:7] - #backup[k] = self.hdr.__getitem__('O'+k) - backup[k] = self.hdr[nkw] - #self.__setattr__(k, self.hdr.__getitem__('O'+k)) - self.__setattr__(k, self.hdr[nkw]) - except KeyError: - print 'Keyword %s not found' % k - - return backup - - def archive_kw(self): - """ - archive original WCS kw before recalculating them. - """ - backup_kw = self.restore() - if backup_kw != {}: - return - else: - # keep the archived kw 8 characters long - cmt = 'archived value' - for kw in basic_wcs: - nkw = 'O'+kw - self.hdr.update(nkw[:7], self.hdr[kw], comment=cmt) - - +
\ No newline at end of file diff --git a/wcsutil/instruments.py b/wcsutil/instruments.py index 2790d8a..a80ec31 100644 --- a/wcsutil/instruments.py +++ b/wcsutil/instruments.py @@ -71,13 +71,11 @@ class ACSWCS(InstrWCS): InstrWCS.__init__(self,hdr0, hdr) self.set_ins_spec_kw() - 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]]} detector = self.primhdr.get('detector', None) - print 'detector', detector self.parity = parity[detector] |