diff options
| -rw-r--r-- | lib/__init__.py | 6 | ||||
| -rw-r--r-- | lib/mappings.py | 38 | ||||
| -rw-r--r-- | lib/utils.py | 2 | ||||
| -rw-r--r-- | updatewcs/apply_corrections.py | 24 | ||||
| -rw-r--r-- | updatewcs/corrections.py | 8 | ||||
| -rw-r--r-- | updatewcs/dgeo.py | 157 | ||||
| -rw-r--r-- | updatewcs/makewcs.py | 2 | ||||
| -rw-r--r-- | wcsutil/__init__.py | 44 | ||||
| -rw-r--r-- | wcsutil/instruments.py | 11 |
9 files changed, 145 insertions, 147 deletions
diff --git a/lib/__init__.py b/lib/__init__.py index c9dca75..4f383fe 100644 --- a/lib/__init__.py +++ b/lib/__init__.py @@ -1,8 +1,10 @@ #import all needed modules here to avoid relative imports -import mappings +#import mappings import utils import distortion import pywcs +from pytools import fileutil +DEGTORAD = fileutil.DEGTORAD +RADTODEG = fileutil.RADTODEG -
\ No newline at end of file diff --git a/lib/mappings.py b/lib/mappings.py deleted file mode 100644 index 4cc4327..0000000 --- a/lib/mappings.py +++ /dev/null @@ -1,38 +0,0 @@ -from pytools import fileutil - -# This dictionary maps an instrument into an instrument class -# The instrument class handles instrument specific keywords - -inst_mappings={'WFPC2': 'WFPC2WCS', - 'ACS': 'ACSWCS' - } - -DEGTORAD = fileutil.DEGTORAD -RADTODEG = fileutil.RADTODEG - -# A dictionary which lists the allowed corrections for each instrument. -# These are the default corrections applied also in the pipeline. -#Dgeo correction is applied separately. -allowed_corrections={'WFPC2': ['MakeWCS','CompSIP', 'VACorr', 'DGEOCorr'], - 'ACS': ['TDDCorr', 'MakeWCS','CompSIP', 'VACorr', 'DGEOCorr'] - } - -# A list of instrument specific keywords -# Every instrument class must have methods which define each of these -# as class attributes. -ins_spec_kw = ['ltv1', 'ltv2', 'parity', 'binned','vafactor', 'chip', - 'naxis1', 'naxis2', 'filter1', 'filter2'] - -# A list of keywords defined in the primary header. -# The HSTWCS class sets this as attributes -prim_hdr_kw = ['detector', 'offtab', 'idctab', 'date-obs', - 'pa_v3', 'ra_targ', 'dec_targ'] - -# These are the keywords which are archived before MakeWCS is run -basic_wcs = ['CD1_1', 'CD1_2', 'CD2_1', 'CD2_2', 'CRVAL1','CRVAL2','CTYPE1', 'CTYPE2', - 'CRPIX1', 'CRPIX2', 'CTYPE1', 'CTYPE2', 'ORIENTAT', 'NAXIS1', 'NAXIS2'] - -dgeo_vals = {'ACS': - {'naxis1':65, 'naxis2':33, 'extver':1, 'crpix1':33.5, 'crpix2':16.5, - 'cdelt1':1., 'cdelt2':1., 'crval1':2048., 'crval2':1024.} - }
\ No newline at end of file diff --git a/lib/utils.py b/lib/utils.py index b29c67c..efe8594 100644 --- a/lib/utils.py +++ b/lib/utils.py @@ -1,6 +1,6 @@ from pytools import parseinput, fileutil import pyfits -from mappings import basic_wcs +from wcsutil.mappings import basic_wcs def restoreWCS(fnames): """ diff --git a/updatewcs/apply_corrections.py b/updatewcs/apply_corrections.py index 6dc7cdf..4758c5a 100644 --- a/updatewcs/apply_corrections.py +++ b/updatewcs/apply_corrections.py @@ -1,13 +1,20 @@ import os import pyfits -from hstwcs.mappings import allowed_corrections +#import allowed_corrections import time from pytools import fileutil import os.path #Note: The order of corrections is important __docformat__ = 'restructuredtext' - + +# A dictionary which lists the allowed corrections for each instrument. +# These are the default corrections applied also in the pipeline. + +allowed_corrections={'WFPC2': ['MakeWCS','CompSIP', 'VACorr', 'DGEOCorr'], + 'ACS': ['TDDCorr', 'MakeWCS', 'CompSIP','VACorr', 'DGEOCorr'] + } + def setCorrections(fname, vacorr=True, tddcorr=True, dgeocorr=True): """ Purpose @@ -19,11 +26,11 @@ def setCorrections(fname, vacorr=True, tddcorr=True, dgeocorr=True): 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') - + # make a copy of this list ! + acorr = allowed_corrections[instrument][:] + if 'VACorr' in acorr and vacorr==False: acorr.remove('VACorr') + if 'TDDCorr' in acorr and tddcorr==False: acorr.remove('TDDCorr') + if 'DGEOCorr' in acorr and dgeocorr==False: acorr.remove('DGEOCorr') return acorr def applyTDDCorr(fname, utddcorr): @@ -103,8 +110,7 @@ def applyDgeoCorr(fname, udgeocorr): applyDGEOCorr = False if isOldStyleDGEO(fname, fdgeo0): - applyDGEOCorr = False - print 'udgeocorr', udgeocorr, applyDGEOCorr + applyDGEOCorr = False return (applyDGEOCorr and udgeocorr) def isOldStyleDGEO(fname, dgname): diff --git a/updatewcs/corrections.py b/updatewcs/corrections.py index 61f73e7..15fd106 100644 --- a/updatewcs/corrections.py +++ b/updatewcs/corrections.py @@ -19,9 +19,9 @@ class TDDCorr(object): Solution of the WFC :Parameters: - `owcs`: HSTWCS object + `ext_wcs`: HSTWCS object An extension HSTWCS object to be modified - `refwcs`: HSTWCS object + `ref_wcs`: HSTWCS object A reference HSTWCS object """ @@ -141,8 +141,8 @@ class CompSIP(object): Bkey="B_%d_%d" % (m,n-m) kw2update[Akey] = sipval[0,0] kw2update[Bkey] = sipval[1,0] - #kw2update['CTYPE1'] = 'RA---TAN-SIP' - #kw2update['CTYPE2'] = 'DEC--TAN-SIP' + kw2update['CTYPE1'] = 'RA---TAN-SIP' + kw2update['CTYPE2'] = 'DEC--TAN-SIP' return kw2update updateWCS = classmethod(updateWCS) diff --git a/updatewcs/dgeo.py b/updatewcs/dgeo.py index d15a485..70e6a98 100644 --- a/updatewcs/dgeo.py +++ b/updatewcs/dgeo.py @@ -1,6 +1,7 @@ import pyfits from pytools import fileutil -from hstwcs.mappings import dgeo_vals +#from hstwcs.mappings import dgeo_vals +import numpy class DGEOCorr(object): """ @@ -13,7 +14,7 @@ class DGEOCorr(object): ========= - Using extensions in the reference file create a WCSDVARR extension and add it to the file object. - - Add record-valued keywords which describe the looikup tables to the + - Add record-valued keywords which describe the lookup tables to the science extension header - Add a keyword 'DGEOFILE' to the science extension header, whose value is the reference file used to create the WCSVARR extension @@ -21,35 +22,38 @@ class DGEOCorr(object): If WCSDVARR extensions exist, subsequent updates will overwrite them. If not, they will be added to the file object. + It is assumed that the DGEO reference files were created to work with IDC tables + but will be applied with SIP coefficients. A transformation is applied to correct + for the fact that the lookup tables will be applied before the first order coefficients + which are in the CD matrix when the SIP convention is used. """ + def updateWCS(cls, fobj): """ :Parameters: `fobj`: pyfits object - Science file, for which a distortion correc tion in a DGEOFILE is available + Science file, for which a distortion correction in a DGEOFILE is available """ assert isinstance(fobj, pyfits.NP_pyfits.HDUList) - #self.fobj = fobj cls.applyDgeoCorr(fobj) dgfile = fobj[0].header['DGEOFILE'] - new_kw = {'DGEOFIEL': dgfile} + new_kw = {'DGEOFILE': dgfile} return new_kw updateWCS = classmethod(updateWCS) - def applyDgeoCorr(cls,fobj): + def applyDgeoCorr(cls, fobj): """ For each science extension in a pyfits file object: - create a WCSDVARR extension - update science header - add/update DGEOFILE keyword """ - - dgeover = 0 dgfile = fileutil.osfn(fobj[0].header['DGEOFILE']) instrument = fobj[0].header.get('INSTRUME', None) + # Map WCSDVARR EXTVER numbers to extension numbers wcsdvarr_ind = cls.getWCSIndex(fobj) for ext in fobj: try: @@ -58,22 +62,31 @@ class DGEOCorr(object): continue if extname == 'sci': extversion = ext.header['EXTVER'] - for ename in ['DX', 'DY']: - dgeover +=1 - - cls.addSciExtKw(ext.header, extver=dgeover, ename=ename) - hdu = cls.createDgeoHDU(instrument, dgeofile=dgfile, dgeover=dgeover,ename=ename, extver=extversion) + header = ext.header + # get the data arrays from the reference file and transform them for use with SIP + dx,dy = cls.getData(dgfile, extversion) + ndx, ndy = cls.transformData(header, dx,dy) + # determine EXTVER for the WCSDVARR extension from the DGEO file (EXTNAME, EXTVER) kw + wcsdvarr_x_version = 2 * extversion -1 + wcsdvarr_y_version = 2 * extversion + + for ename in zip(['DX', 'DY'], [wcsdvarr_x_version,wcsdvarr_y_version],[ndx, ndy]): + cls.addSciExtKw(header, wdvarr_ver=ename[1], dgeo_name=ename[0]) + hdu = cls.createDgeoHDU(header, dgeofile=dgfile, wdvarr_ver=ename[1],dgeo_name=ename[0], data=ename[2],extver=extversion) if wcsdvarr_ind: - fobj[wcsdvarr_ind[dgeover]] = hdu + fobj[wcsdvarr_ind[ename[1]]] = hdu else: fobj.append(hdu) - + + applyDgeoCorr = classmethod(applyDgeoCorr) - + def getWCSIndex(cls, fobj): """ - Returns a mapping of WCSDVARR 'EXTVER' and pyfits.HDUList indeces. - + If fobj has WCSDVARR extensions: + returns a mapping of their EXTVER kw are mapped to extension numbers + if fobj does not have WCSDVARR extensions: + an empty dictionary is returned. """ wcsd = {} for e in range(len(fobj)): @@ -88,17 +101,12 @@ class DGEOCorr(object): getWCSIndex = classmethod(getWCSIndex) - def addSciExtKw(cls, hdr, extver=None, ename=None): + def addSciExtKw(cls, hdr, wdvarr_ver=None, dgeo_name=None): """ - Adds kw to sci extension to define dgeo correction extensions - kw to be added to sci ext: - CPERRORj - CPDISj - DPj.EXTVER - DPj.NAXES - DPj.AXIS.i (i=DP1.NAXES) + Adds kw to sci extension to define WCSDVARR lookup table extensions + """ - if ename =='DX': + if dgeo_name =='DX': j=1 else: j=2 @@ -110,10 +118,10 @@ class DGEOCorr(object): dpaxis1 = 'DP%s.' %j+'AXIS.1' dpaxis2 = 'DP%s.' %j+'AXIS.2' keys = [cperror, cpdis, dpext, dpnaxes, dpaxis1, dpaxis2] - values = {cperror: 0.0, cpdis: 'Lookup', dpext: extver, dpnaxes: 2, + values = {cperror: 0.0, cpdis: 'Lookup', dpext: wdvarr_ver, dpnaxes: 2, dpaxis1: 1, dpaxis2: 2} - comments = {cperror: 'Maximum error of dgeo correction for axis %s' % (extver/2), + comments = {cperror: 'Maximum error of dgeo correction for axis %s' % (wdvarr_ver/2), cpdis: 'Prior distortion funcion type', dpext: 'Version number of WCSDVARR extension containing lookup distortion table', dpnaxes: 'Number of independent variables in distortion function', @@ -126,46 +134,75 @@ class DGEOCorr(object): addSciExtKw = classmethod(addSciExtKw) - def getDgeoData(cls, dgfile=None, ename=None, extver=1): + def getData(cls,dgfile, extver): """ - 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(cls, instrument, dgeofile=None, dgeover=1, ename=None,extver=1): + Get the data arrays from the reference DGEO files + """ + xdata = pyfits.getdata(dgfile, ext=('DX',extver)) + ydata = pyfits.getdata(dgfile, ext=('DY',extver)) + return xdata, ydata + getData = classmethod(getData) + + def transformData(cls, header, dx, dy): + """ + Transform the DGEO data arrays for use with SIP + """ + coeffs = cls.getCoeffs(header) + idcscale = header['IDCSCALE'] + sclcoeffs = numpy.linalg.inv(coeffs)/idcscale + ndx, ndy = numpy.dot(sclcoeffs, [dx.ravel(), dy.ravel()]) + ndx.shape = dx.shape + ndy.shape=dy.shape + return ndx, ndy + transformData = classmethod(transformData) + + def getCoeffs(cls, header): + """ + Return a matrix of the first order IDC coefficients. + """ + try: + ocx10 = header['OCX10'] + ocx11 = header['OCX11'] + ocy10 = header['OCY10'] + ocy11 = header['OCY11'] + except AttributeError: + print 'First order IDCTAB coefficients are not available.\n' + print 'Cannot convert SIP to IDC coefficients.\n' + return None + return numpy.array([[ocx11, ocx10], [ocy11,ocy10]], dtype=numpy.float32) + + getCoeffs = classmethod(getCoeffs) + + def createDgeoHDU(cls, sciheader, dgeofile=None, wdvarr_ver=1, dgeo_name=None,data = 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 = cls.createDgeoHdr(instrument, **dgeokw) - data = cls.getDgeoData(dgfile=dgeofile, ename=ename, extver=extver) + hdr = cls.createDgeoHdr(sciheader, dgeofile=dgeofile, wdvarr_ver=wdvarr_ver, dgeoname=dgeo_name, extver=extver) hdu=pyfits.ImageHDU(header=hdr, data=data) return hdu + createDgeoHDU = classmethod(createDgeoHDU) - def createDgeoHdr(cls, instr, **kw): + def createDgeoHdr(cls, sciheader, dgeofile, wdvarr_ver, dgeoname, extver): """ - Creates a header for the dgeo extension based on dgeo file + Creates a header for the WCSDVARR extension based on the DGEO reference 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_vals = dgeo_vals[instr] - naxis1 = kw['naxis1'] or instr_vals['naxis1'] - naxis2 = kw['naxis2'] or instr_vals['naxis2'] - extver = kw['extver'] or instr_vals['extver'] - crpix1 = kw['crpix1'] or instr_vals['crpix1'] - crpix2 = kw['crpix2'] or instr_vals['crpix2'] - cdelt1 = kw['cdelt1'] or instr_vals['cdelt1'] - cdelt2 = kw['cdelt2'] or instr_vals['cdelt2'] - crval1 = kw['crval1'] or instr_vals['crval1'] - crval2 = kw['crval2'] or instr_vals['crval2'] - + """ + dgeo_header = pyfits.getheader(dgeofile, ext=(dgeoname,extver)) + sci_naxis1 = sciheader['NAXIS1'] + sci_naxis2 = sciheader['NAXIS2'] + sci_crpix1 = sciheader['CRPIX1'] + sci_crpix2 = sciheader['CRPIX2'] + + naxis1 = dgeo_header['naxis1'] + naxis2 = dgeo_header['naxis2'] + extver = dgeo_header['extver'] + crpix1 = naxis1/2. + crpix2 = naxis2/2. + cdelt1 = sci_naxis1/naxis1 + cdelt2 = sci_naxis2/naxis2 + crval1 = sci_crpix1 + crval2 = sci_crpix2 keys = ['XTENSION','BITPIX','NAXIS','NAXIS1','NAXIS2', 'EXTNAME','EXTVER','PCOUNT','GCOUNT','CRPIX1', 'CDELT1','CRVAL1','CRPIX2','CDELT2','CRVAL2'] @@ -192,7 +229,7 @@ class DGEOCorr(object): 'NAXIS1': naxis1, 'NAXIS2': naxis2, 'EXTNAME': 'WCSDVARR', - 'EXTVER': extver, + 'EXTVER': wdvarr_ver, 'PCOUNT': 0, 'GCOUNT': 1, 'CRPIX1': crpix1, diff --git a/updatewcs/makewcs.py b/updatewcs/makewcs.py index 4b06d95..449cf2a 100644 --- a/updatewcs/makewcs.py +++ b/updatewcs/makewcs.py @@ -1,5 +1,5 @@ #from .. mappings import DEGTORAD, RADTODEG -from hstwcs.mappings import DEGTORAD, RADTODEG +from hstwcs import DEGTORAD, RADTODEG import numpy from math import sin, sqrt, pow, cos, asin, atan2,pi from hstwcs import utils diff --git a/wcsutil/__init__.py b/wcsutil/__init__.py index 8b06d7a..40a3649 100644 --- a/wcsutil/__init__.py +++ b/wcsutil/__init__.py @@ -6,10 +6,12 @@ import instruments from hstwcs.distortion import models import numpy as N from pytools import fileutil +from pytools.fileutil import DEGTORAD, RADTODEG #from .. mappings import inst_mappings, ins_spec_kw, DEGTORAD, RADTODEG, basic_wcs -from hstwcs.mappings import inst_mappings, ins_spec_kw, DEGTORAD, RADTODEG -from hstwcs.mappings import basic_wcs, prim_hdr_kw +import mappings +from mappings import inst_mappings, ins_spec_kw +from mappings import basic_wcs, prim_hdr_kw __docformat__ = 'restructuredtext' @@ -75,8 +77,6 @@ class HSTWCS(WCS): WCS.__init__(self, ehdr, fobj=fobj) self.setHDR0kw(hdr0, ehdr) - #self.detector = self.setDetector(hdr0) - self.setInstrSpecKw(hdr0, ehdr) self.pscale = self.setPscale() self.orientat = self.setOrient() @@ -90,33 +90,17 @@ class HSTWCS(WCS): def setHDR0kw(self, primhdr, ehdr): - if primhdr == None: - # we are given only an extension header - header = ehdr - elif ehdr == None: - header = primhdr - else: - hcards = primhdr.ascardlist() - hcards.extend(ehdr.ascardlist()) - header = pyfits.Header(cards = hcards) # Set attributes from kw defined in the primary header. - self.instrument = header.get('INSTRUME', None) - self.offtab = header.get('OFFTAB', None) - self.idctab = header.get('IDCTAB', None) - self.date_obs = header.get('DATE-OBS', None) - self.pav3 = header.get('PA_V3', None) - self.ra_targ = header.get('RA_TARG', None) - self.dec_targ = header.get('DEC_TARG', None) - self.detector = header.get('DETECTOR', None) - self.filename = header.get('FILENAME', "") - """ - def setDetector(self, header): - # Set detector attribute for instuments which have more than one detector - if self.instrument in ['ACS', 'WFC3']: - return header.get('DETECTOR', None) - else: - return None - """ + self.instrument = primhdr.get('INSTRUME', None) + self.offtab = primhdr.get('OFFTAB', None) + self.idctab = primhdr.get('IDCTAB', None) + self.date_obs = primhdr.get('DATE-OBS', None) + self.pav3 = primhdr.get('PA_V3', None) + self.ra_targ = primhdr.get('RA_TARG', None) + self.dec_targ = primhdr.get('DEC_TARG', None) + self.filename = primhdr.get('FILENAME', "") + self.detector = primhdr.get('DETECTOR', None) + def readIDCCoeffs(self, header): """ Reads in first order IDCTAB coefficients if present in the header diff --git a/wcsutil/instruments.py b/wcsutil/instruments.py index e2c2b0b..092597c 100644 --- a/wcsutil/instruments.py +++ b/wcsutil/instruments.py @@ -1,7 +1,7 @@ import pyfits import numpy as N #from .. mappings import ins_spec_kw -from hstwcs.mappings import ins_spec_kw, prim_hdr_kw +from mappings import ins_spec_kw, prim_hdr_kw class InstrWCS(object): """ @@ -29,6 +29,7 @@ class InstrWCS(object): self.set_binned() self.set_chip() self.set_parity() + self.set_extver() def set_filter1(self): try: @@ -77,7 +78,13 @@ class InstrWCS(object): self.chip = self.exthdr.get('CCDCHIP', 1) except: self.chip = 1 - + + def set_extver(self): + try: + self.extver = self.exthdr.get('EXTVER', 1) + except: + self.extver = 1 + def set_parity(self): self.parity = [[1.0,0.0],[0.0,-1.0]] |