diff options
| -rw-r--r-- | updatewcs/__init__.py | 2 | ||||
| -rw-r--r-- | updatewcs/apply_corrections.py | 38 | ||||
| -rw-r--r-- | updatewcs/dgeo.py | 86 | ||||
| -rw-r--r-- | wcsutil/__init__.py | 55 | ||||
| -rw-r--r-- | wcsutil/instruments.py | 3 |
5 files changed, 132 insertions, 52 deletions
diff --git a/updatewcs/__init__.py b/updatewcs/__init__.py index 00eda0f..ca1dd7a 100644 --- a/updatewcs/__init__.py +++ b/updatewcs/__init__.py @@ -1,9 +1,7 @@ import os import pyfits -#from .. wcsutil import HSTWCS from stwcs.wcsutil import HSTWCS -#from .. mappings import allowed_corrections from stwcs import utils import corrections, makewcs import dgeo, det2im diff --git a/updatewcs/apply_corrections.py b/updatewcs/apply_corrections.py index 5e55767..665d0f4 100644 --- a/updatewcs/apply_corrections.py +++ b/updatewcs/apply_corrections.py @@ -1,9 +1,9 @@ import os import pyfits -#import allowed_corrections import time from pytools import fileutil import os.path + #Note: The order of corrections is important __docformat__ = 'restructuredtext' @@ -24,15 +24,26 @@ def setCorrections(fname, vacorr=True, tddcorr=True, dgeocorr=True, d2imcorr=Tru ======= 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. + for the instrument. """ instrument = pyfits.getval(fname, 'INSTRUME') # make a copy of this list ! acorr = allowed_corrections[instrument][:] + + #First check if idctab was updated + if not foundIDCTAB(fname):# and not foundSIP(fname): + acorr.remove('TDDCorr') + acorr.remove('MakeWCS') + acorr.remove('CompSIP') + elif not foundIDCTAB(fname) and foundSIP(fname): + print 'IDCTAB not found, using SIP coefficients\n' + + #check if idctab is present on disk if 'VACorr' in acorr and vacorr==False: acorr.remove('VACorr') if 'TDDCorr' in acorr: tddcorr = applyTDDCorr(fname, tddcorr) if tddcorr == False: acorr.remove('TDDCorr') + if 'DGEOCorr' in acorr: dgeocorr = applyDgeoCorr(fname, dgeocorr) if dgeocorr == False: acorr.remove('DGEOCorr') @@ -41,6 +52,29 @@ def setCorrections(fname, vacorr=True, tddcorr=True, dgeocorr=True, d2imcorr=Tru if d2imcorr == False: acorr.remove('DET2IMCorr') return acorr +def foundIDCTAB(fname): + idctab_found = True + try: + idctab = fileutil.osfn(pyfits.getval(fname, 'IDCTAB')) + if idctab == 'N/A' or idctab == "": + idctab_found = False + if os.path.exists(idctab): + idctab_found = True + else: + idctab_found = False + except KeyError: + idctab_found = False + return idctab_found + +def foundSIP(fname): + sip_found = False + try: + idcscale = pyfits.getval(fname, 'IDCSCALE') + sip_found = True + except KeyError: + sip_found = False + return sip_found + def applyTDDCorr(fname, utddcorr): """ The default value of tddcorr for all ACS images is True. diff --git a/updatewcs/dgeo.py b/updatewcs/dgeo.py index 5dcb5ad..9242a11 100644 --- a/updatewcs/dgeo.py +++ b/updatewcs/dgeo.py @@ -61,21 +61,24 @@ class DGEOCorr(object): continue if extname == 'sci': extversion = ext.header['EXTVER'] - ccdchip = cls.get_ccdchip(fobj, extversion) + ccdchip = cls.get_ccdchip(fobj, extname='SCI', extver=extversion) binned = utils.getBinning(fobj, extversion) header = ext.header # get the data arrays from the reference file and transform them for use with SIP dx,dy = cls.getData(dgfile, ccdchip) - ndx, ndy = cls.transformData(header, dx,dy) + idccoeffs = cls.getIDCCoeffs(header) + if idccoeffs != None: + dx, dy = cls.transformData(dx,dy, idccoeffs) + # Determine EXTVER for the WCSDVARR extension from the DGEO file (EXTNAME, EXTVER) kw. # This is used to populate DPj.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]) + for ename in zip(['DX', 'DY'], [wcsdvarr_x_version,wcsdvarr_y_version],[dx, dy]): + cls.addSciExtKw(header, wdvarr_ver=ename[1], dgeo_extname=ename[0]) hdu = cls.createDgeoHDU(header, dgeofile=dgfile, \ - wdvarr_ver=ename[1], dgeo_name=ename[0], data=ename[2],ccdchip=ccdchip, binned=binned) + wdvarr_ver=ename[1], dgeo_extname=ename[0], data=ename[2],ccdchip=ccdchip, binned=binned) if wcsdvarr_ind: fobj[wcsdvarr_ind[ename[1]]] = hdu else: @@ -105,12 +108,12 @@ class DGEOCorr(object): getWCSIndex = classmethod(getWCSIndex) - def addSciExtKw(cls, hdr, wdvarr_ver=None, dgeo_name=None): + def addSciExtKw(cls, hdr, wdvarr_ver=None, dgeo_extname=None): """ Adds kw to sci extension to define WCSDVARR lookup table extensions """ - if dgeo_name =='DX': + if dgeo_extname =='DX': j=1 else: j=2 @@ -146,7 +149,7 @@ class DGEOCorr(object): dgf = pyfits.open(dgfile) for ext in dgf: dgextname = ext.header.get('EXTNAME',"") - dgccdchip = ext.header.get('CCDCHIP',0) + dgccdchip = ext.header.get('CCDCHIP',1) if dgextname == 'DX' and dgccdchip == ccdchip: xdata = ext.data.copy() continue @@ -159,47 +162,51 @@ class DGEOCorr(object): return xdata, ydata getData = classmethod(getData) - def transformData(cls, header, dx, dy): + def transformData(cls, dx, dy, coeffs): """ Transform the DGEO data arrays for use with SIP """ - coeffs = cls.getCoeffs(header) - idcscale = header['IDCSCALE'] - sclcoeffs = np.linalg.inv(coeffs/idcscale) - ndx, ndy = np.dot(sclcoeffs, [dx.ravel(), dy.ravel()]) + ndx, ndy = np.dot(coeffs, [dx.ravel(), dy.ravel()]) ndx.shape = dx.shape ndy.shape=dy.shape return ndx, ndy + transformData = classmethod(transformData) - def getCoeffs(cls, header): + def getIDCCoeffs(cls, header): """ - Return a matrix of the first order IDC coefficients. + Return a matrix of the scaled first order IDC coefficients. """ try: ocx10 = header['OCX10'] ocx11 = header['OCX11'] ocy10 = header['OCY10'] ocy11 = header['OCY11'] - except AttributeError: + coeffs = np.array([[ocx11, ocx10], [ocy11,ocy10]], dtype=np.float32) + except KeyError: print 'First order IDCTAB coefficients are not available.\n' print 'Cannot convert SIP to IDC coefficients.\n' return None - return np.array([[ocx11, ocx10], [ocy11,ocy10]], dtype=np.float32) + try: + idcscale = header['IDCSCALE'] + except KeyError: + idcscale = 1 + + return np.linalg.inv(coeffs/idcscale) - getCoeffs = classmethod(getCoeffs) + getIDCCoeffs = classmethod(getIDCCoeffs) - def createDgeoHDU(cls, sciheader, dgeofile=None, wdvarr_ver=1, dgeo_name=None,data = None, ccdchip=1, binned=1): + def createDgeoHDU(cls, sciheader, dgeofile=None, wdvarr_ver=1, dgeo_extname=None,data = None, ccdchip=1, binned=1): """ Creates an HDU to be added to the file object. """ - hdr = cls.createDgeoHdr(sciheader, dgeofile=dgeofile, wdvarr_ver=wdvarr_ver, dgeoname=dgeo_name, ccdchip=ccdchip, binned=binned) + hdr = cls.createDgeoHdr(sciheader, dgeofile=dgeofile, wdvarr_ver=wdvarr_ver, dg_extname=dgeo_extname, ccdchip=ccdchip, binned=binned) hdu=pyfits.ImageHDU(header=hdr, data=data) return hdu createDgeoHDU = classmethod(createDgeoHDU) - def createDgeoHdr(cls, sciheader, dgeofile, wdvarr_ver, dgeoname, ccdchip, binned): + def createDgeoHdr(cls, sciheader, dgeofile, wdvarr_ver, dg_extname, ccdchip, binned): """ Creates a header for the WCSDVARR extension based on the DGEO reference file and sci extension header. The goal is to always work in image coordinates @@ -209,9 +216,15 @@ class DGEOCorr(object): """ dgf = pyfits.open(dgeofile) for ext in dgf: - dgextname = ext.header.get('EXTNAME', "") - dgccdchip = ext.header.get('CCDCHIP', 0) - if dgextname == dgeoname and dgccdchip == ccdchip: + #for i in range(len(dgf)): + try: + dgextname = ext.header['EXTNAME'] + dgextver = ext.header['EXTVER'] + except KeyError: + continue + #dgccdchip = ext.header.get('CCDCHIP', 0) + dgccdchip = cls.get_ccdchip(dgf, extname=dgextname, extver=dgextver) + if dgextname == dg_extname and dgccdchip == ccdchip: dgeo_header = ext.header break else: @@ -219,7 +232,7 @@ class DGEOCorr(object): dgf.close() naxis = pyfits.getval(dgeofile, ext=1, key='NAXIS') - ccdchip = dgeo_header['CCDCHIP'] + ccdchip = dgextname #dgeo_header['CCDCHIP'] kw = { 'NAXIS': 'Size of the axis', 'CRPIX': 'Coordinate system reference pixel', @@ -274,23 +287,22 @@ class DGEOCorr(object): createDgeoHdr = classmethod(createDgeoHdr) - def get_ccdchip(cls, fobj, extver): + def get_ccdchip(cls, fobj, extname, extver): """ - Currently this is only needed for ACS, since this is the only instrument - with DGEO correction. The rest is here for completeness. + Given a science file or dgeo file determine CCDCHIP """ - dgfile = fileutil.osfn(fobj[0].header['DGEOFILE']) - if fobj[0].header['INSTRUME'] == 'ACS': - return fobj['SCI', extver].header['CCDCHIP'] - elif obj[0].header['INSTRUME'] == 'WFC3': - return fobj['SCI', extver].header['CCDCHIP'] + ccdchip = 1 + if fobj[0].header['INSTRUME'] == 'ACS' and fobj[0].header['DETECTOR'] == 'WFC': + ccdchip = fobj[extname, extver].header['CCDCHIP'] + elif fobj[0].header['INSTRUME'] == 'WFC3' and fobj[0].header['DETECTOR'] == 'UVIS': + ccdchip = fobj[extname, extver].header['CCDCHIP'] elif fobj[0].header['INSTRUME'] == 'WFPC2': - return fobj[0].header['DETECTOR'] + ccdchip = fobj[extname, extver].header['DETECTOR'] elif fobj[0].header['INSTRUME'] == 'STIS': - return fobj['SCI', extver].header['DETECTOR'] + ccdchip = fobj[extname, extver].header['DETECTOR'] elif fobj[0].header['INSTRUME'] == 'NICMOS': - return fobj['SCI', extver].header['CAMERA'] - + ccdchip = fobj[extname, extver].header['CAMERA'] + return ccdchip get_ccdchip = classmethod(get_ccdchip)
\ No newline at end of file diff --git a/wcsutil/__init__.py b/wcsutil/__init__.py index 12d302b..4284cf5 100644 --- a/wcsutil/__init__.py +++ b/wcsutil/__init__.py @@ -1,17 +1,15 @@ -#from .. pywcs.sip import SIP -from pywcs import WCS, DistortionLookupTable +import os.path +from pywcs import WCS import pyfits import instruments -#from .. distortion import models -from stwcs.distortion import models +from stwcs.distortion import models, coeff_converter import numpy as np from pytools import fileutil from pytools.fileutil import DEGTORAD, RADTODEG -#from .. mappings import inst_mappings, ins_spec_kw, DEGTORAD, RADTODEG, basic_wcs import mappings from mappings import inst_mappings, ins_spec_kw -from mappings import basic_wcs, prim_hdr_kw +from mappings import basic_wcs __docformat__ = 'restructuredtext' @@ -181,6 +179,45 @@ class HSTWCS(WCS): def readModel(self, update=False, header=None): + + if self.idctab == None: + #Keyword idctab is not present in header - check for sip coefficients + if header.has_key('IDCSCALE'): + self.readModelFromHeader(header) + else: + print 'Distortion model is not available\n' + return + elif not os.path.exists(fileutil.osfn(self.idctab)): + if header.has_key('IDCSCALE'): + self.readModelFromHeader(header) + else: + print 'Distortion model is not available\n' + return + else: + self.readModelFromIDCTAB(header=header, update=update) + + def readModelFromHeader(self, header): + # Recreate idc model from SIP coefficients and header kw + model = models.GeometryModel() + cx, cy = coeff_converter.sip2idc(self) + model.cx = cx + model.cy = cy + model.name = "sip" + model.norder = header['A_ORDER'] + + refpix = {} + refpix['XREF'] = header['IDCXREF'] + refpix['YREF'] = header['IDCYREF'] + refpix['PSCALE'] = header['IDCSCALE'] + refpix['V2REF'] = header['IDCV2REF'] + refpix['V3REF'] = header['IDCV3REF'] + refpix['THETA'] = header['IDCTHETA'] + model.refpix = refpix + + self.idcmodel = model + + + def readModelFromIDCTAB(self, header=hdr, update=False): """ Purpose ======= @@ -189,11 +226,11 @@ class HSTWCS(WCS): If header is provided and update is True, some IDC model kw will be recorded in the header. """ - if self.idctab == None or self.date_obs == None: - print 'idctab or date_obs not available\n' + if self.date_obs == None: + print 'date_obs not available\n' self.idcmodel = None return - if self.filter1 ==None and self.filter2 == None: + if self.filter1 == None and self.filter2 == None: 'No filter information available\n' self.idcmodel = None return diff --git a/wcsutil/instruments.py b/wcsutil/instruments.py index afdf1a9..df92d54 100644 --- a/wcsutil/instruments.py +++ b/wcsutil/instruments.py @@ -1,6 +1,5 @@ import pyfits -#from .. mappings import ins_spec_kw -from mappings import ins_spec_kw, prim_hdr_kw +from mappings import ins_spec_kw class InstrWCS(object): """ |