diff options
Diffstat (limited to 'updatewcs/dgeo.py')
| -rw-r--r-- | updatewcs/dgeo.py | 157 |
1 files changed, 97 insertions, 60 deletions
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, |