diff options
| -rw-r--r-- | updatewcs/__init__.py | 10 | ||||
| -rw-r--r-- | updatewcs/corrections.py | 58 | ||||
| -rw-r--r-- | updatewcs/makewcs.py | 122 | ||||
| -rw-r--r-- | wcsutil/__init__.py | 14 |
4 files changed, 103 insertions, 101 deletions
diff --git a/updatewcs/__init__.py b/updatewcs/__init__.py index 107a9f5..2f77949 100644 --- a/updatewcs/__init__.py +++ b/updatewcs/__init__.py @@ -99,15 +99,15 @@ def makecorr(fname, allowed_corr): for extn in f: # 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) + ref_wcs = HSTWCS(primhdr, ref_hdr, fobj=f) + ref_wcs.readModel(ref_hdr) hdr = extn.header - owcs = HSTWCS(primhdr, hdr, fobj=f) + ext_wcs = HSTWCS(primhdr, hdr, fobj=f) utils.write_archive(hdr) - owcs.readModel(hdr) + ext_wcs.readModel(hdr) for c in allowed_corr: corr_klass = corrections.__getattribute__(c) - kw2update = corr_klass.updateWCS(owcs, refwcs) + kw2update = corr_klass.updateWCS(ext_wcs, ref_wcs) for kw in kw2update: hdr.update(kw, kw2update[kw]) diff --git a/updatewcs/corrections.py b/updatewcs/corrections.py index ac49833..dc4eeff 100644 --- a/updatewcs/corrections.py +++ b/updatewcs/corrections.py @@ -25,21 +25,21 @@ class TDDCorr(object): A reference HSTWCS object """ - def updateWCS(cls, wcs, refwcs): + def updateWCS(cls, ext_wcs, ref_wcs): """ - Calculates alpha and beta for ACS/WFC data. - Writes 2 new kw to the extension header: TDDALPHA and TDDBETA - sets TDDCORR to COMPLETE """ - alpha, beta = cls.set_alpha_beta(wcs) - cls.apply_tdd2idc(wcs, alpha, beta) + alpha, beta = cls.set_alpha_beta(ext_wcs) + cls.apply_tdd2idc(ext_wcs, alpha, beta) newkw = {'TDDALPHA': alpha, 'TDDBETA':beta, 'TDDCORR': 'COMPLETE'} return newkw updateWCS = classmethod(updateWCS) - def apply_tdd2idc(cls,wcs, alpha, beta): + def apply_tdd2idc(cls,ext_wcs, alpha, beta): """ Applies TDD to the idctab coefficients of a ACS/WFC observation. This should be always the first correction. @@ -53,23 +53,23 @@ class TDDCorr(object): 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 = 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() + xshape, yshape = ext_wcs.idcmodel.cx.shape, ext_wcs.idcmodel.cy.shape + icxy = numpy.dot(abmat2,[ext_wcs.idcmodel.cx.ravel(),ext_wcs.idcmodel.cy.ravel()]) + ext_wcs.idcmodel.cx = icxy[0] + ext_wcs.idcmodel.cy = icxy[1] + ext_wcs.idcmodel.cx.shape = xshape + ext_wcs.idcmodel.cy.shape = yshape + ext_wcs.wcs.set() apply_tdd2idc = classmethod(apply_tdd2idc) - def set_alpha_beta(cls, wcs): + def set_alpha_beta(cls, ext_wcs): """ Compute the time dependent distortion skew terms default date of 2004.5 = 2004-7-1 """ dday = 2004.5 - year,month,day = wcs.date_obs.split('-') + year,month,day = ext_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 @@ -89,18 +89,18 @@ class VACorr(object): """ - 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) + 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) - 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() + 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], + 'CRVAL1':ext_wcs.wcs.crval[0], 'CRVAL2':ext_wcs.wcs.crval[1]} + ext_wcs.wcs.set() else: pass return kw2update @@ -115,15 +115,15 @@ class CompSIP(object): Compute SIP coefficients from idc table coefficients. """ - def updateWCS(cls, wcs, refwcs): + def updateWCS(cls, ext_wcs, ref_wcs): kw2update = {} - order = wcs.idcmodel.norder + order = ext_wcs.idcmodel.norder kw2update['A_ORDER'] = order kw2update['B_ORDER'] = order - pscale = wcs.idcmodel.refpix['PSCALE'] + pscale = ext_wcs.idcmodel.refpix['PSCALE'] - cx = wcs.idcmodel.cx - cy = wcs.idcmodel.cy + cx = ext_wcs.idcmodel.cx + cy = ext_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) diff --git a/updatewcs/makewcs.py b/updatewcs/makewcs.py index ddb4d06..3df2fd4 100644 --- a/updatewcs/makewcs.py +++ b/updatewcs/makewcs.py @@ -19,72 +19,72 @@ class MakeWCS(object): """ - def updateWCS(cls, wcs, refwcs): + def updateWCS(cls, ext_wcs, ref_wcs): """ recomputes the basic WCS kw """ - ltvoff, offshift = cls.getOffsets(wcs) + ltvoff, offshift = cls.getOffsets(ext_wcs) - 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], + cls.uprefwcs(ext_wcs, ref_wcs, ltvoff, offshift) + cls.upextwcs(ext_wcs, ref_wcs, ltvoff, offshift) + + 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], + 'CRVAL1': ext_wcs.wcs.crval[0], + 'CRVAL2': ext_wcs.wcs.crval[1], + 'CRPIX1': ext_wcs.wcs.crpix[0], + 'CRPIX2': ext_wcs.wcs.crpix[1], } return kw2update updateWCS = classmethod(updateWCS) - def upextwcs(self, wcs, refwcs, loff, offsh): + def upextwcs(self, ext_wcs, ref_wcs, loff, offsh): """ updates an extension wcs """ ltvoffx, ltvoffy = loff[0], loff[1] offshiftx, offshifty = offsh[0], offsh[1] - ltv1 = wcs.ltv1 - ltv2 = wcs.ltv2 + ltv1 = ext_wcs.ltv1 + ltv2 = ext_wcs.ltv2 if ltv1 != 0. or ltv2 != 0.: - 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) + offsetx = backup_wcs['CRPIX1'] - ltv1 - ext_wcs.idcmodel.refpix['XREF'] + offsety = backup_wcs['CRPIX2'] - ltv2 - ext_wcs.idcmodel.refpix['YREF'] + fx,fy = ext_wcs.idcmodel.shift(ext_wcs.idcmodel.cx,ext_wcs.idcmodel.cy,offsetx,offsety) else: - fx, fy = wcs.idcmodel.cx, wcs.idcmodel.cy + fx, fy = ext_wcs.idcmodel.cx, ext_wcs.idcmodel.cy - ridcmodel = refwcs.idcmodel - v2 = wcs.idcmodel.refpix['V2REF'] - v3 = wcs.idcmodel.refpix['V3REF'] - v2ref = refwcs.idcmodel.refpix['V2REF'] + ridcmodel = ref_wcs.idcmodel + v2 = ext_wcs.idcmodel.refpix['V2REF'] + v3 = ext_wcs.idcmodel.refpix['V3REF'] + v2ref = ref_wcs.idcmodel.refpix['V2REF'] - v3ref = refwcs.idcmodel.refpix['V3REF'] - R_scale = refwcs.idcmodel.refpix['PSCALE']/3600.0 + v3ref = ref_wcs.idcmodel.refpix['V3REF'] + R_scale = ref_wcs.idcmodel.refpix['PSCALE']/3600.0 off = sqrt((v2-v2ref)**2 + (v3-v3ref)**2)/(R_scale*3600.0) if v3 == v3ref: theta=0.0 else: - theta = atan2(wcs.parity[0][0]*(v2-v2ref), wcs.parity[1][1]*(v3-v3ref)) + theta = atan2(ext_wcs.parity[0][0]*(v2-v2ref), ext_wcs.parity[1][1]*(v3-v3ref)) - if refwcs.idcmodel.refpix['THETA']: theta += refwcs.idcmodel.refpix['THETA']*pi/180.0 + if ref_wcs.idcmodel.refpix['THETA']: theta += ref_wcs.idcmodel.refpix['THETA']*pi/180.0 dX=(off*sin(theta)) + offshiftx dY=(off*cos(theta)) + offshifty px = numpy.array([[dX,dY]]) - 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 + newcrval = ref_wcs.all_pix2sky_fits(px)[0] + newcrpix = numpy.array([ext_wcs.idcmodel.refpix['XREF'] + ltvoffx, + ext_wcs.idcmodel.refpix['YREF'] + ltvoffy]) + ext_wcs.wcs.crval = newcrval + ext_wcs.wcs.crpix = newcrpix # Account for subarray offset # Angle of chip relative to chip - if wcs.idcmodel.refpix['THETA']: - dtheta = wcs.idcmodel.refpix['THETA'] - refwcs.idcmodel.refpix['THETA'] + if ext_wcs.idcmodel.refpix['THETA']: + dtheta = ext_wcs.idcmodel.refpix['THETA'] - ref_wcs.idcmodel.refpix['THETA'] else: dtheta = 0.0 @@ -109,13 +109,13 @@ class MakeWCS(object): px = numpy.array([[dX+dXX,dY+dYX]]) # Transform to sky coordinates - wc = refwcs.pixel2world_fits(px) + wc = ref_wcs.all_pix2sky_fits(px) a = wc[0,0] b = wc[0,1] px = numpy.array([[dX+dXY,dY+dYY]]) - wc = refwcs.pixel2world_fits(px) + wc = ref_wcs.all_pix2sky_fits(px) c = wc[0,0] d = wc[0,1] @@ -125,58 +125,60 @@ class MakeWCS(object): cd21 = utils.diff_angles(b,newcrval[1]) cd22 = utils.diff_angles(d,newcrval[1]) cd = numpy.array([[cd11, cd12], [cd21, cd22]]) - wcs.cd = cd + ext_wcs.wcs.cd = cd #??? upextwcs = classmethod(upextwcs) - def uprefwcs(self, wcs, refwcs, loff, offsh): + def uprefwcs(self, ext_wcs, ref_wcs, loff, offsh): """ Updates the reference chip """ ltvoffx, ltvoffy = loff[0], loff[1] offshift = offsh - dec = refwcs.crval[1] + dec = ref_wcs.wcs.crval[1] # Get an approximate reference position on the sky - rref = numpy.array([[refwcs.idcmodel.refpix['XREF']+ltvoffx, - refwcs.idcmodel.refpix['YREF']+ltvoffy]]) + rref = numpy.array([[ref_wcs.idcmodel.refpix['XREF']+ltvoffx, + ref_wcs.idcmodel.refpix['YREF']+ltvoffy]]) - crval = refwcs.pixel2world_fits(rref) + crval = ref_wcs.all_pix2sky_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(wcs.pav3,dec,refwcs.idcmodel.refpix['V2REF'], - refwcs.idcmodel.refpix['V3REF']) + pv = troll(ext_wcs.pav3,dec,ref_wcs.idcmodel.refpix['V2REF'], + ref_wcs.idcmodel.refpix['V3REF']) # Add the chip rotation angle - if refwcs.idcmodel.refpix['THETA']: - pv += refwcs.idcmodel.refpix['THETA'] + if ref_wcs.idcmodel.refpix['THETA']: + pv += ref_wcs.idcmodel.refpix['THETA'] # Set values for the rest of the reference WCS - refwcs.crval = crval[0] - refwcs.crpix = numpy.array([0.0,0.0])+offsh - parity = refwcs.parity - R_scale = refwcs.idcmodel.refpix['PSCALE']/3600.0 + ref_wcs.wcs.crval = crval[0] + ref_wcs.wcs.crpix = numpy.array([0.0,0.0])+offsh + parity = ref_wcs.parity + R_scale = ref_wcs.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]]) - refwcs.cd = rcd - refwcs.set() + print 'cd shape', rcd.shape, ref_wcs.wcs.cd.shape + ref_wcs.wcs.cd = rcd + ref_wcs.wcs.set() uprefwcs = classmethod(uprefwcs) - def getOffsets(cls, wcs): - ltv1 = wcs.ltv1 - ltv2 = wcs.ltv2 + def getOffsets(cls, ext_wcs): + ltv1 = ext_wcs.ltv1 + ltv2 = ext_wcs.ltv2 - offsetx = wcs.crpix[0] - ltv1 - wcs.idcmodel.refpix['XREF'] - offsety = wcs.crpix[1] - ltv2 - wcs.idcmodel.refpix['YREF'] + offsetx = ext_wcs.wcs.crpix[0] - ltv1 - ext_wcs.idcmodel.refpix['XREF'] + offsety = ext_wcs.wcs.crpix[1] - ltv2 - ext_wcs.idcmodel.refpix['YREF'] - shiftx = wcs.idcmodel.refpix['XREF'] + ltv1 - shifty = wcs.idcmodel.refpix['YREF'] + ltv2 + shiftx = ext_wcs.idcmodel.refpix['XREF'] + ltv1 + shifty = ext_wcs.idcmodel.refpix['YREF'] + ltv2 if ltv1 != 0. or ltv2 != 0.: ltvoffx = ltv1 + offsetx ltvoffy = ltv2 + offsety diff --git a/wcsutil/__init__.py b/wcsutil/__init__.py index 9482a48..f9e9c18 100644 --- a/wcsutil/__init__.py +++ b/wcsutil/__init__.py @@ -1,5 +1,5 @@ #from .. pywcs.sip import SIP -from pywcs.sip import SIP +from pywcs import WCS import pyfits import instruments #from .. distortion import models @@ -12,7 +12,7 @@ from hstwcs.mappings import basic_wcs, prim_hdr_kw __docformat__ = 'restructuredtext' -class HSTWCS(SIP): +class HSTWCS(WCS): """ Purpose ======= @@ -32,7 +32,7 @@ class HSTWCS(SIP): `fobj`: PyFITS HDUList object or None pyfits file object """ - SIP.__init__(self, ehdr, fobj=fobj) + WCS.__init__(self, ehdr, fobj=fobj) self.setHDR0kw(hdr0) self.detector = self.setDetector(hdr0) self.inst_kw = ins_spec_kw @@ -76,14 +76,14 @@ class HSTWCS(SIP): # 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] + cd11 = self.wcs.cd[0][0] + cd21 = self.wcs.cd[1][0] return N.sqrt(N.power(cd11,2)+N.power(cd21,2)) * 3600. def setOrient(self): # Recompute ORIENTAT - cd12 = self.cd[0][1] - cd22 = self.cd[1][1] + cd12 = self.wcs.cd[0][1] + cd22 = self.wcs.cd[1][1] return RADTODEG(N.arctan2(cd12,cd22)) def verifyPa_V3(self): |