diff options
Diffstat (limited to 'distortion')
| -rw-r--r-- | distortion/__init__.py | 0 | ||||
| -rw-r--r-- | distortion/coeff_converter.py | 141 | ||||
| -rw-r--r-- | distortion/models.py | 361 | ||||
| -rw-r--r-- | distortion/mutil.py | 663 | ||||
| -rw-r--r-- | distortion/utils.py | 204 |
5 files changed, 0 insertions, 1369 deletions
diff --git a/distortion/__init__.py b/distortion/__init__.py deleted file mode 100644 index e69de29..0000000 --- a/distortion/__init__.py +++ /dev/null diff --git a/distortion/coeff_converter.py b/distortion/coeff_converter.py deleted file mode 100644 index f2eb4ad..0000000 --- a/distortion/coeff_converter.py +++ /dev/null @@ -1,141 +0,0 @@ -from __future__ import division # confidence high - -import numpy as np -import pyfits -import pywcs - -def sip2idc(wcs): - """ - Converts SIP style coefficients to IDCTAB coefficients. - - :Parameters: - `wcs`: pyfits.Header or pywcs.WCS object - """ - if isinstance(wcs,pyfits.Header): - ocx10 = wcs.get('OCX10', None) - ocx11 = wcs.get('OCX11', None) - ocy10 = wcs.get('OCY10', None) - ocy11 = wcs.get('OCY11', None) - order = wcs.get('A_ORDER', None) - sipa, sipb = _read_sip_kw(wcs) - if None in [ocx10, ocx11, ocy10, ocy11, sipa, sipb]: - print 'Cannot convert SIP to IDC coefficients.\n' - return None, None - elif isinstance(wcs,pywcs.WCS): - try: - ocx10 = wcs.ocx10 - ocx11 = wcs.ocx11 - ocy10 = wcs.ocy10 - ocy11 = wcs.ocy11 - except AttributeError: - print 'First order IDCTAB coefficients are not available.\n' - print 'Cannot convert SIP to IDC coefficients.\n' - return None, None - try: - sipa = wcs.sip.a - sipb = wcs.sip.b - except AttributeError: - print 'SIP coefficients are not available.' - print 'Cannot convert SIP to IDC coefficients.\n' - return None, None - try: - order = wcs.sip.a_order - except AttributeError: - print 'SIP model order unknown, exiting ...\n' - return None, None - - else: - print 'Input to sip2idc must be a PyFITS header or a wcsutil.HSTWCS object\n' - return - - - if None in [ocx10, ocx11, ocy10, ocy11]: - print 'First order IDC coefficients not found, exiting ...\n' - return - idc_coeff = np.array([[ocx11, ocx10], [ocy11, ocy10]]) - cx = np.zeros((order+1,order+1), dtype=np.double) - cy = np.zeros((order+1,order+1), dtype=np.double) - for n in range(order+1): - for m in range(order+1): - if n >= m and n>=2: - sipval = np.array([[sipa[m,n-m]],[sipb[m,n-m]]]) - idcval = np.dot(idc_coeff, sipval) - cx[n,m] = idcval[0] - cy[n,m] = idcval[1] - - cx[1,0] = ocx10 - cx[1,1] = ocx11 - cy[1,0] = ocy10 - cy[1,1] = ocy11 - - return cx, cy - -def _read_sip_kw(header): - """ - Reads SIP header keywords and returns an array of coefficients. - - If no SIP header keywords are found, None is returned. - """ - if header.has_key("A_ORDER"): - if not header.has_key("B_ORDER"): - raise ValueError( - "A_ORDER provided without corresponding B_ORDER " - "keyword for SIP distortion") - - m = int(header["A_ORDER"]) - a = np.zeros((m+1, m+1), np.double) - for i in range(m+1): - for j in range(m-i+1): - a[i, j] = header.get("A_%d_%d" % (i, j), 0.0) - - m = int(header["B_ORDER"]) - b = np.zeros((m+1, m+1), np.double) - for i in range(m+1): - for j in range(m-i+1): - b[i, j] = header.get("B_%d_%d" % (i, j), 0.0) - elif header.has_key("B_ORDER"): - raise ValueError( - "B_ORDER provided without corresponding A_ORDER " - "keyword for SIP distortion") - else: - a = None - b = None - - return a , b - - -""" -def idc2sip(wcsobj, idctab = None): - if isinstance(wcs,pywcs.WCS): - try: - cx10 = wcsobj.ocx10 - cx11 = wcsobj.cx11 - cy10 = wcsobj.cy10 - cy11 = wcsobj.cy11 - except AttributeError: - print - try: - order = wcs.sip.a_order - except AttributeError: - print 'SIP model order unknown, exiting ...\n' - return - else: - print 'Input to sip2idc must be a PyFITS header or a wcsutil.HSTWCS object\n' - return - - if None in [ocx10, ocx11, ocy10, ocy11]: - print 'First order IDC coefficients not found, exiting ...\n' - return - idc_coeff = np.array([[wcsobj.cx11, wcsobj.cx10], [wcsobj.cy11, wcsobj.cy10]]) - cx = numpy.zeros((order+1,order+1), dtype=numpy.double) - cy = numpy.zeros((order+1,order+1), dtype=numpy.double) - for n in range(order+1): - for m in range(order+1): - if n >= m and n>=2: - sipval = numpy.array([[wcsobj.sip.a[n,m]],[wcsobj.sip.b[n,m]]]) - idcval = numpy.dot(idc_coeff, sipval) - cx[m,n-m] = idcval[0] - cy[m,n-m] = idcval[1] - - return cx, cy -"""
\ No newline at end of file diff --git a/distortion/models.py b/distortion/models.py deleted file mode 100644 index 96d5d72..0000000 --- a/distortion/models.py +++ /dev/null @@ -1,361 +0,0 @@ -from __future__ import division # confidence high - -import types -# Import PyDrizzle utility modules -import mutil -import numpy as np -import mutil -from mutil import combin - -yes = True -no = False - -################# -# -# -# Geometry/Distortion Classes -# -# -################# - -class GeometryModel: - """ - Base class for Distortion model. - There will be a separate class for each type of - model/filetype used with drizzle, i.e., IDCModel and - DrizzleModel. - - Each class will know how to apply the distortion to a - single point and how to convert coefficients to an input table - suitable for the drizzle task. - - Coefficients will be stored in CX,CY arrays. - """ - # - # - # - # - # - # - # - NORDER = 3 - - def __init__(self): - " This will open the given file and determine its type and norder." - - # Method to read in coefficients from given table and - # populate the n arrays 'cx' and 'cy'. - # This will be different for each type of input file, - # IDCTAB vs. drizzle table. - - # Set these up here for all sub-classes to use... - # But, calculate norder and cx,cy arrays in detector specific classes. - self.cx = None - self.cy = None - self.refpix = None - self.norder = self.NORDER - # Keep track of computed zero-point for distortion coeffs - self.x0 = None - self.y0 = None - - # default values for these attributes - self.direction = 'forward' - - self.pscale = 1.0 - - def shift(self,cx,cy,xs,ys): - """ - Shift reference position of coefficients to new center - where (xs,ys) = old-reference-position - subarray/image center. - This will support creating coeffs files for drizzle which will - be applied relative to the center of the image, rather than relative - to the reference position of the chip. - """ - - _cxs = np.zeros(shape=cx.shape,dtype=cx.dtype) - _cys = np.zeros(shape=cy.shape,dtype=cy.dtype) - _k = self.norder + 1 - # loop over each input coefficient - for m in xrange(_k): - for n in xrange(_k): - if m >= n: - # For this coefficient, shift by xs/ys. - _ilist = range(m, _k) - # sum from m to k - for i in _ilist: - _jlist = range(n, i - (m-n)+1) - # sum from n to i-(m-n) - for j in _jlist: - _cxs[m,n] += cx[i,j]*combin(j,n)*combin((i-j),(m-n))*pow(xs,(j-n))*pow(ys,((i-j)-(m-n))) - _cys[m,n] += cy[i,j]*combin(j,n)*combin((i-j),(m-n))*pow(xs,(j-n))*pow(ys,((i-j)-(m-n))) - return _cxs,_cys - - def convert(self, tmpname, xref=None,yref=None,delta=yes): - """ - Open up an ASCII file, output coefficients in drizzle - format after converting them as necessary. - First, normalize these coefficients to what drizzle expects - Normalize the coefficients by the MODEL/output plate scale. - - 16-May-2002: - Revised to work with higher order polynomials by John Blakeslee. - 27-June-2002: - Added ability to shift coefficients to new center for support - of subarrays. - """ - cx = self.cx/self.pscale - cy = self.cy/self.pscale - x0 = self.refpix['XDELTA'] + cx[0,0] - y0 = self.refpix['YDELTA'] + cy[0,0] - #xr = self.refpix['XREF'] - #yr = self.refpix['YREF'] - xr = self.refpix['CHIP_XREF'] - yr = self.refpix['CHIP_YREF'] - - - - ''' - if xref != None: - # Shift coefficients for use with drizzle - _xs = xref - self.refpix['XREF'] + 1.0 - _ys = yref - self.refpix['YREF'] + 1.0 - - - if _xs != 0 or _ys != 0: - cxs,cys= self.shift(cx, cy, _xs, _ys) - cx = cxs - cy = cys - - # We only want to apply this shift to coeffs - # for subarray images. - if delta == no: - cxs[0,0] = cxs[0,0] - _xs - cys[0,0] = cys[0,0] - _ys - - # Now, apply only the difference introduced by the distortion.. - # i.e., (undistorted - original) shift. - x0 += cxs[0,0] - y0 += cys[0,0] - ''' - self.x0 = x0 #+ 1.0 - self.y0 = y0 #+ 1.0 - - # Now, write out the coefficients into an ASCII - # file in 'drizzle' format. - lines = [] - - - lines.append('# Polynomial distortion coefficients\n') - lines.append('# Extracted from "%s" \n'%self.name) - lines.append('refpix %f %f \n'%(xr,yr)) - if self.norder==3: - lines.append('cubic\n') - elif self.norder==4: - lines.append('quartic\n') - elif self.norder==5: - lines.append('quintic\n') - else: - raise ValueError, "Drizzle cannot handle poly distortions of order %d"%self.norder - - str = '%16.8f %16.8g %16.8g %16.8g %16.8g \n'% (x0,cx[1,1],cx[1,0],cx[2,2],cx[2,1]) - lines.append(str) - str = '%16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cx[2,0],cx[3,3],cx[3,2],cx[3,1],cx[3,0]) - lines.append(str) - if self.norder>3: - str = '%16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cx[4,4],cx[4,3],cx[4,2],cx[4,1],cx[4,0]) - lines.append(str) - if self.norder>4: - str = '%16.8g %16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cx[5,5],cx[5,4],cx[5,3],cx[5,2],cx[5,1],cx[5,0]) - lines.append(str) - lines.append("\n") - - str = '%16.8f %16.8g %16.8g %16.8g %16.8g \n'% (y0,cy[1,1],cy[1,0],cy[2,2],cy[2,1]) - lines.append(str) - str = '%16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cy[2,0],cy[3,3],cy[3,2],cy[3,1],cy[3,0]) - lines.append(str) - if self.norder>3: - str = '%16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cy[4,4],cy[4,3],cy[4,2],cy[4,1],cy[4,0]) - lines.append(str) - if self.norder>4: - str = '%16.8g %16.8g %16.8g %16.8g %16.8g %16.8g \n'% (cy[5,5],cy[5,4],cy[5,3],cy[5,2],cy[5,1],cy[5,0]) - lines.append(str) - - output = open(tmpname,'w') - output.writelines(lines) - output.close() - - - def apply(self, pixpos,scale=1.0,order=None): - """ - Apply coefficients to a pixel position or a list of positions. - This should be the same for all coefficients tables. - Return the geometrically-adjusted position - in arcseconds from the reference position as a tuple (x,y). - - Compute delta from reference position - """ - - """ - scale actually is a ratio of pscale/self.model.pscale - what is pscale? - """ - if self.cx == None: - return pixpos[:,0],pixpos[:,1] - - if order is None: - order = self.norder - - # Apply in the same way that 'drizzle' would... - _cx = self.cx / (self.pscale * scale) - _cy = self.cy / (self.pscale * scale) - _convert = no - _p = pixpos - - # Do NOT include any zero-point terms in CX,CY here - # as they should not be scaled by plate-scale like rest - # of coeffs... This makes the computations consistent - # with 'drizzle'. WJH 17-Feb-2004 - _cx[0,0] = 0. - _cy[0,0] = 0. - - if isinstance(_p,types.ListType) or isinstance(_p,types.TupleType): - _p = np.array(_p,dtype=np.float64) - _convert = yes - - dxy = _p - (self.refpix['XREF'],self.refpix['YREF']) - # Apply coefficients from distortion model here... - c = _p * 0. - for i in range(order+1): - for j in range(i+1): - c[:,0] = c[:,0] + _cx[i][j] * pow(dxy[:,0],j) * pow(dxy[:,1],(i-j)) - c[:,1] = c[:,1] + _cy[i][j] * pow(dxy[:,0],j) * pow(dxy[:,1],(i-j)) - xc = c[:,0] - yc = c[:,1] - - # Convert results back to same form as original input - if _convert: - xc = xc.tolist() - yc = yc.tolist() - # If a single tuple was input, return just a single tuple - if len(xc) == 1: - xc = xc[0] - yc = yc[0] - - return xc,yc - - def setPScaleCoeffs(self,pscale): - self.cx[1,1] = pscale - self.cy[1,0] = pscale - - self.refpix['PSCALE'] = pscale - self.pscale = pscale - - -class IDCModel(GeometryModel): - """ - This class will open the IDCTAB, select proper row based on - chip/direction and populate cx,cy arrays. - We also need to read in SCALE, XCOM,YCOM, XREF,YREF as well. - """ - def __init__(self, idcfile, date=None, chip=1, direction='forward', - filter1='CLEAR1',filter2='CLEAR2',offtab=None, binned=1): - GeometryModel.__init__(self) - # - # Norder must be derived from the coeffs file itself, - # then the arrays can be setup. Thus, it needs to be - # done in the sub-class, not in the base class. - # Read in table. - # Populate cx,cy,scale, and other variables here. - # - self.name = idcfile - self.cx,self.cy,self.refpix,self.norder = mutil.readIDCtab(idcfile, - chip=chip,direction=direction,filter1=filter1,filter2=filter2, - date=date, offtab=offtab) - - if self.refpix.has_key('empty_model') and self.refpix['empty_model']: - pass - else: - self.refpix['PSCALE'] = self.refpix['PSCALE'] * binned - self.cx = self.cx * binned - self.cy = self.cy * binned - self.refpix['XREF'] = self.refpix['XREF'] / binned - self.refpix['YREF'] = self.refpix['YREF'] / binned - self.refpix['XSIZE'] = self.refpix['XSIZE'] / binned - self.refpix['YSIZE'] = self.refpix['YSIZE'] / binned - - self.pscale = self.refpix['PSCALE'] - - -class WCSModel(GeometryModel): - """ - This class sets up a distortion model based on coefficients - found in the image header. - """ - def __init__(self,header,rootname): - GeometryModel.__init__(self) - - - if header.has_key('rootname'): - self.name = header['rootname'] - else: - self.name = rootname - # Initialize all necessary distortion arrays with - # default model... - #self.cx,self.cy,self.refpix,self.order = mutil.defaultModel() - - # Read in values from header, and update distortion arrays. - self.cx,self.cy,self.refpix,self.norder = mutil.readWCSCoeffs(header) - - self.pscale = self.refpix['PSCALE'] - - - -class DrizzleModel(GeometryModel): - """ - This class will read in an ASCII Cubic - drizzle coeffs file and populate the cx,cy arrays. - """ - - def __init__(self, idcfile, scale = None): - GeometryModel.__init__(self) - # - # We now need to read in the file, populate cx,cy, and - # other variables as necessary. - # - self.name = idcfile - self.cx,self.cy,self.refpix,self.norder = mutil.readCubicTable(idcfile) - - # scale is the ratio wcs.pscale/model.pscale. - # model.pscale for WFPC2 is passed from REFDATA. - # This is needed for WFPC2 binned data. - - if scale != None: - self.pscale = scale - else: - self.pscale = self.refpix['PSCALE'] - - """ - The above definition looks wrong. - In one case it's a ratio in the other it's pscale. - - """ - -class TraugerModel(GeometryModel): - """ - This class will read in the ASCII Trauger coeffs - file, convert them to SIAF coefficients, then populate - the cx,cy arrays. - """ - NORDER = 3 - - def __init__(self, idcfile,lam): - GeometryModel.__init__(self) - self.name = idcfile - self.cx,self.cy,self.refpix,self.norder = mutil.readTraugerTable(idcfile,lam) - self.pscale = self.refpix['PSCALE'] - # - # Read in file here. - # Populate cx,cy, and other variables. - # - - diff --git a/distortion/mutil.py b/distortion/mutil.py deleted file mode 100644 index 2934d17..0000000 --- a/distortion/mutil.py +++ /dev/null @@ -1,663 +0,0 @@ -from __future__ import division # confidence high - -from pytools import fileutil -import numpy as np -import string -import calendar - -# Set up IRAF-compatible Boolean values -yes = True -no = False - -# This function read the IDC table and generates the two matrices with -# the geometric correction coefficients. -# -# INPUT: FITS object of open IDC table -# OUTPUT: coefficient matrices for Fx and Fy -# -#### If 'tabname' == None: This should return a default, undistorted -#### solution. -# - -def readIDCtab (tabname, chip=1, date=None, direction='forward', - filter1=None,filter2=None, offtab=None): - - """ - Read IDCTAB, and optional OFFTAB if sepcified, and generate - the two matrices with the geometric correction coefficients. - - If tabname == None, then return a default, undistorted solution. - If offtab is specified, dateobs also needs to be given. - - """ - - # Return a default geometry model if no IDCTAB filename - # is given. This model will not distort the data in any way. - if tabname == None: - print 'Warning: No IDCTAB specified! No distortion correction will be applied.' - return defaultModel() - - # Implement default values for filters here to avoid the default - # being overwritten by values of None passed by user. - if filter1 == None or filter1.find('CLEAR') == 0: - filter1 = 'CLEAR' - if filter2 == None or filter2.find('CLEAR') == 0: - filter2 = 'CLEAR' - - # Insure that tabname is full filename with fully expanded - # IRAF variables; i.e. 'jref$mc41442gj_idc.fits' should get - # expanded to '/data/cdbs7/jref/mc41442gj_idc.fits' before - # being used here. - # Open up IDC table now... - try: - ftab = fileutil.openImage(tabname) - except: - err_str = "------------------------------------------------------------------------ \n" - err_str += "WARNING: the IDCTAB geometric distortion file specified in the image \n" - err_str += "header was not found on disk. Please verify that your environment \n" - err_str += "variable ('jref'/'uref'/'oref'/'nref') has been correctly defined. If \n" - err_str += "you do not have the IDCTAB file, you may obtain the latest version \n" - err_str += "of it from the relevant instrument page on the STScI HST website: \n" - err_str += "http://www.stsci.edu/hst/ For WFPC2, STIS and NICMOS data, the \n" - err_str += "present run will continue using the old coefficients provided in \n" - err_str += "the Dither Package (ca. 1995-1998). \n" - err_str += "------------------------------------------------------------------------ \n" - raise IOError,err_str - - #First thing we need, is to read in the coefficients from the IDC - # table and populate the Fx and Fy matrices. - - if ftab['PRIMARY'].header.has_key('DETECTOR'): - detector = ftab['PRIMARY'].header['DETECTOR'] - else: - if ftab['PRIMARY'].header.has_key('CAMERA'): - detector = str(ftab['PRIMARY'].header['CAMERA']) - else: - detector = 1 - # First, read in TDD coeffs if present - phdr = ftab['PRIMARY'].header - instrument = phdr['INSTRUME'] - if instrument == 'ACS' and detector == 'WFC': - skew_coeffs = read_tdd_coeffs(phdr) - else: - skew_coeffs = None - - # Set default filters for SBC - if detector == 'SBC': - if filter1 == 'CLEAR': - filter1 = 'F115LP' - filter2 = 'N/A' - if filter2 == 'CLEAR': - filter2 = 'N/A' - - # Read FITS header to determine order of fit, i.e. k - norder = ftab['PRIMARY'].header['NORDER'] - if norder < 3: - order = 3 - else: - order = norder - - fx = np.zeros(shape=(order+1,order+1),dtype=np.float64) - fy = np.zeros(shape=(order+1,order+1),dtype=np.float64) - - #Determine row from which to get the coefficients. - # How many rows do we have in the table... - fshape = ftab[1].data.shape - colnames = ftab[1].data._names - row = -1 - - # Loop over all the rows looking for the one which corresponds - # to the value of CCDCHIP we are working on... - for i in xrange(fshape[0]): - - try: - # Match FILTER combo to appropriate row, - #if there is a filter column in the IDCTAB... - if 'FILTER1' in colnames and 'FILTER2' in colnames: - - filt1 = ftab[1].data.field('FILTER1')[i] - if filt1.find('CLEAR') > -1: filt1 = filt1[:5] - - filt2 = ftab[1].data.field('FILTER2')[i] - if filt2.find('CLEAR') > -1: filt2 = filt2[:5] - else: - if 'OPT_ELEM' in colnames: - filt1 = ftab[1].data.field('OPT_ELEM') - if filt1.find('CLEAR') > -1: filt1 = filt1[:5] - else: - filt1 = filter1 - - if 'FILTER' in colnames: - _filt = ftab[1].data.field('FILTER')[i] - if _filt.find('CLEAR') > -1: _filt = _filt[:5] - if 'OPT_ELEM' in colnames: - filt2 = _filt - else: - filt1 = _filt - filt2 = 'CLEAR' - else: - filt2 = filter2 - except: - # Otherwise assume all rows apply and compare to input filters... - filt1 = filter1 - filt2 = filter2 - - if 'DETCHIP' in colnames: - detchip = ftab[1].data.field('DETCHIP')[i] - if not str(detchip).isdigit(): - detchip = 1 - else: - detchip = 1 - - if 'DIRECTION' in colnames: - direct = string.strip(string.lower(ftab[1].data.field('DIRECTION')[i])) - else: - direct = 'forward' - - if filt1 == filter1.strip() and filt2 == filter2.strip(): - if direct == direction.strip(): - if int(detchip) == int(chip) or int(detchip) == -999: - row = i - break - - joinstr = ',' - if 'CLEAR' in filter1: - f1str = '' - joinstr = '' - else: - f1str = filter1.strip() - if 'CLEAR' in filter2: - f2str = '' - joinstr = '' - else: - f2str = filter2.strip() - filtstr = (joinstr.join([f1str,f2str])).strip() - if row < 0: - err_str = '\nProblem finding row in IDCTAB! Could not find row matching:\n' - err_str += ' CHIP: '+str(detchip)+'\n' - err_str += ' FILTERS: '+filstr+'\n' - ftab.close() - del ftab - raise LookupError,err_str - else: - print '- IDCTAB: Distortion model from row',str(row+1),'for chip',detchip,':',filtstr - - # Read in V2REF and V3REF: this can either come from current table, - # or from an OFFTAB if time-dependent (i.e., for WFPC2) - theta = None - if 'V2REF' in colnames: - v2ref = ftab[1].data.field('V2REF')[row] - v3ref = ftab[1].data.field('V3REF')[row] - else: - # Read V2REF/V3REF from offset table (OFFTAB) - if offtab: - v2ref,v3ref,theta = readOfftab(offtab, date, chip=detchip) - else: - v2ref = 0.0 - v3ref = 0.0 - - if theta == None: - if 'THETA' in colnames: - theta = ftab[1].data.field('THETA')[row] - else: - theta = 0.0 - - refpix = {} - refpix['XREF'] = ftab[1].data.field('XREF')[row] - refpix['YREF'] = ftab[1].data.field('YREF')[row] - refpix['XSIZE'] = ftab[1].data.field('XSIZE')[row] - refpix['YSIZE'] = ftab[1].data.field('YSIZE')[row] - refpix['PSCALE'] = round(ftab[1].data.field('SCALE')[row],8) - refpix['V2REF'] = v2ref - refpix['V3REF'] = v3ref - refpix['THETA'] = theta - refpix['XDELTA'] = 0.0 - refpix['YDELTA'] = 0.0 - refpix['DEFAULT_SCALE'] = yes - refpix['centered'] = no - refpix['skew_coeffs'] = skew_coeffs - - # Now that we know which row to look at, read coefficients into the - # numeric arrays we have set up... - # Setup which column name convention the IDCTAB follows - # either: A,B or CX,CY - if 'CX10' in ftab[1].data._names: - cxstr = 'CX' - cystr = 'CY' - else: - cxstr = 'A' - cystr = 'B' - - for i in xrange(norder+1): - if i > 0: - for j in xrange(i+1): - xcname = cxstr+str(i)+str(j) - ycname = cystr+str(i)+str(j) - fx[i,j] = ftab[1].data.field(xcname)[row] - fy[i,j] = ftab[1].data.field(ycname)[row] - - ftab.close() - del ftab - - # If CX11 is 1.0 and not equal to the PSCALE, then the - # coeffs need to be scaled - - if fx[1,1] == 1.0 and abs(fx[1,1]) != refpix['PSCALE']: - fx *= refpix['PSCALE'] - fy *= refpix['PSCALE'] - - # Return arrays and polynomial order read in from table. - # NOTE: XREF and YREF are stored in Fx,Fy arrays respectively. - return fx,fy,refpix,order -# -# -# Time-dependent skew correction coefficients (only ACS/WFC) -# -# -def read_tdd_coeffs(phdr): - ''' Read in the TDD related keywords from the PRIMARY header of the IDCTAB - ''' - skew_coeffs = {} - skew_coeffs['TDDORDER'] = 0 - skew_coeffs['TDD_DATE'] = "" - skew_coeffs['TDD_A'] = None - skew_coeffs['TDD_B'] = None - - if phdr.has_key("TDDORDER"): - n = int(phdr["TDDORDER"]) - else: - print 'TDDORDER kw not present, using default TDD correction' - return None - - a = np.zeros((n+1,), np.float64) - b = np.zeros((n+1,), np.float64) - for i in range(n+1): - a[i] = phdr.get(("TDD_A%d" % i), 0.0) - b[i] = phdr.get(("TDD_B%d" % i), 0.0) - if (a==0).all() and (b==0).all(): - print 'Warning: TDD_A and TDD_B coeffiecients have values of 0, \n \ - but TDDORDER is %d.' % TDDORDER - - skew_coeffs['TDDORDER'] = n - skew_coeffs['TDD_DATE'] = phdr['TDD_DATE'] - skew_coeffs['TDD_A'] = a - skew_coeffs['TDD_B'] = b - - return skew_coeffs - -def readOfftab(offtab, date, chip=None): - - -#Read V2REF,V3REF from a specified offset table (OFFTAB). -# Return a default geometry model if no IDCTAB filenam e -# is given. This model will not distort the data in any way. - - if offtab == None: - return 0.,0. - - # Provide a default value for chip - if chip: - detchip = chip - else: - detchip = 1 - - # Open up IDC table now... - try: - ftab = fileutil.openImage(offtab) - except: - raise IOError,"Offset table '%s' not valid as specified!" % offtab - - #Determine row from which to get the coefficients. - # How many rows do we have in the table... - fshape = ftab[1].data.shape - colnames = ftab[1].data._names - row = -1 - - row_start = None - row_end = None - - v2end = None - v3end = None - date_end = None - theta_end = None - - num_date = convertDate(date) - # Loop over all the rows looking for the one which corresponds - # to the value of CCDCHIP we are working on... - for ri in xrange(fshape[0]): - i = fshape[0] - ri - 1 - if 'DETCHIP' in colnames: - detchip = ftab[1].data.field('DETCHIP')[i] - else: - detchip = 1 - - obsdate = convertDate(ftab[1].data.field('OBSDATE')[i]) - - # If the row is appropriate for the chip... - # Interpolate between dates - if int(detchip) == int(chip) or int(detchip) == -999: - if num_date <= obsdate: - date_end = obsdate - v2end = ftab[1].data.field('V2REF')[i] - v3end = ftab[1].data.field('V3REF')[i] - theta_end = ftab[1].data.field('THETA')[i] - row_end = i - continue - - if row_end == None and (num_date > obsdate): - date_end = obsdate - v2end = ftab[1].data.field('V2REF')[i] - v3end = ftab[1].data.field('V3REF')[i] - theta_end = ftab[1].data.field('THETA')[i] - row_end = i - continue - - if num_date > obsdate: - date_start = obsdate - v2start = ftab[1].data.field('V2REF')[i] - v3start = ftab[1].data.field('V3REF')[i] - theta_start = ftab[1].data.field('THETA')[i] - row_start = i - break - - ftab.close() - del ftab - - if row_start == None and row_end == None: - print 'Row corresponding to DETCHIP of ',detchip,' was not found!' - raise LookupError - elif row_start == None: - print '- OFFTAB: Offset defined by row',str(row_end+1) - else: - print '- OFFTAB: Offset interpolated from rows',str(row_start+1),'and',str(row_end+1) - - # Now, do the interpolation for v2ref, v3ref, and theta - if row_start == None or row_end == row_start: - # We are processing an observation taken after the last calibration - date_start = date_end - v2start = v2end - v3start = v3end - _fraction = 0. - theta_start = theta_end - else: - _fraction = float((num_date - date_start)) / float((date_end - date_start)) - - v2ref = _fraction * (v2end - v2start) + v2start - v3ref = _fraction * (v3end - v3start) + v3start - theta = _fraction * (theta_end - theta_start) + theta_start - - return v2ref,v3ref,theta - -def readWCSCoeffs(header): - - #Read distortion coeffs from WCS header keywords and - #populate distortion coeffs arrays. - - # Read in order for polynomials - _xorder = header['a_order'] - _yorder = header['b_order'] - order = max(max(_xorder,_yorder),3) - - fx = np.zeros(shape=(order+1,order+1),dtype=np.float64) - fy = np.zeros(shape=(order+1,order+1),dtype=np.float64) - - # Read in CD matrix - _cd11 = header['cd1_1'] - _cd12 = header['cd1_2'] - _cd21 = header['cd2_1'] - _cd22 = header['cd2_2'] - _cdmat = np.array([[_cd11,_cd12],[_cd21,_cd22]]) - _theta = np.arctan2(-_cd12,_cd22) - _rotmat = np.array([[np.cos(_theta),np.sin(_theta)], - [-np.sin(_theta),np.cos(_theta)]]) - _rCD = np.dot(_rotmat,_cdmat) - _skew = np.arcsin(-_rCD[1][0] / _rCD[0][0]) - _scale = _rCD[0][0] * np.cos(_skew) * 3600. - _scale2 = _rCD[1][1] * 3600. - - # Set up refpix - refpix = {} - refpix['XREF'] = header['crpix1'] - refpix['YREF'] = header['crpix2'] - refpix['XSIZE'] = header['naxis1'] - refpix['YSIZE'] = header['naxis2'] - refpix['PSCALE'] = _scale - refpix['V2REF'] = 0. - refpix['V3REF'] = 0. - refpix['THETA'] = RADTODEG(_theta) - refpix['XDELTA'] = 0.0 - refpix['YDELTA'] = 0.0 - refpix['DEFAULT_SCALE'] = yes - refpix['centered'] = yes - - - # Set up template for coeffs keyword names - cxstr = 'A_' - cystr = 'B_' - # Read coeffs into their own matrix - for i in xrange(_xorder+1): - for j in xrange(i+1): - xcname = cxstr+str(j)+'_'+str(i-j) - if header.has_key(xcname): - fx[i,j] = header[xcname] - - # Extract Y coeffs separately as a different order may - # have been used to fit it. - for i in xrange(_yorder+1): - for j in xrange(i+1): - ycname = cystr+str(j)+'_'+str(i-j) - if header.has_key(ycname): - fy[i,j] = header[ycname] - - # Now set the linear terms - fx[0][0] = 1.0 - fy[0][0] = 1.0 - - return fx,fy,refpix,order - - -def readTraugerTable(idcfile,wavelength): - - # Return a default geometry model if no coefficients filename - # is given. This model will not distort the data in any way. - if idcfile == None: - return fileutil.defaultModel() - - # Trauger coefficients only result in a cubic file... - order = 3 - numco = 10 - a_coeffs = [0] * numco - b_coeffs = [0] * numco - indx = _MgF2(wavelength) - - ifile = open(idcfile,'r') - # Search for the first line of the coefficients - _line = fileutil.rAsciiLine(ifile) - while string.lower(_line[:7]) != 'trauger': - _line = fileutil.rAsciiLine(ifile) - # Read in each row of coefficients,split them into their values, - # and convert them into cubic coefficients based on - # index of refraction value for the given wavelength - # Build X coefficients from first 10 rows of Trauger coefficients - j = 0 - while j < 20: - _line = fileutil.rAsciiLine(ifile) - if _line == '': continue - _lc = string.split(_line) - if j < 10: - a_coeffs[j] = float(_lc[0])+float(_lc[1])*(indx-1.5)+float(_lc[2])*(indx-1.5)**2 - else: - b_coeffs[j-10] = float(_lc[0])+float(_lc[1])*(indx-1.5)+float(_lc[2])*(indx-1.5)**2 - j = j + 1 - - ifile.close() - del ifile - - # Now, convert the coefficients into a Numeric array - # with the right coefficients in the right place. - # Populate output values now... - fx = np.zeros(shape=(order+1,order+1),dtype=np.float64) - fy = np.zeros(shape=(order+1,order+1),dtype=np.float64) - # Assign the coefficients to their array positions - fx[0,0] = 0. - fx[1] = np.array([a_coeffs[2],a_coeffs[1],0.,0.],dtype=np.float64) - fx[2] = np.array([a_coeffs[5],a_coeffs[4],a_coeffs[3],0.],dtype=np.float64) - fx[3] = np.array([a_coeffs[9],a_coeffs[8],a_coeffs[7],a_coeffs[6]],dtype=np.float64) - fy[0,0] = 0. - fy[1] = np.array([b_coeffs[2],b_coeffs[1],0.,0.],dtype=np.float64) - fy[2] = np.array([b_coeffs[5],b_coeffs[4],b_coeffs[3],0.],dtype=np.float64) - fy[3] = np.array([b_coeffs[9],b_coeffs[8],b_coeffs[7],b_coeffs[6]],dtype=np.float64) - - # Used in Pattern.computeOffsets() - refpix = {} - refpix['XREF'] = None - refpix['YREF'] = None - refpix['V2REF'] = None - refpix['V3REF'] = None - refpix['XDELTA'] = 0. - refpix['YDELTA'] = 0. - refpix['PSCALE'] = None - refpix['DEFAULT_SCALE'] = no - refpix['centered'] = yes - - return fx,fy,refpix,order - - -def readCubicTable(idcfile): - # Assumption: this will only be used for cubic file... - order = 3 - # Also, this function does NOT perform any scaling on - # the coefficients, it simply passes along what is found - # in the file as is... - - # Return a default geometry model if no coefficients filename - # is given. This model will not distort the data in any way. - if idcfile == None: - return fileutil.defaultModel() - - ifile = open(idcfile,'r') - # Search for the first line of the coefficients - _line = fileutil.rAsciiLine(ifile) - - _found = no - while _found == no: - if _line[:7] in ['cubic','quartic','quintic'] or _line[:4] == 'poly': - found = yes - break - _line = fileutil.rAsciiLine(ifile) - - # Read in each row of coefficients, without line breaks or newlines - # split them into their values, and create a list for A coefficients - # and another list for the B coefficients - _line = fileutil.rAsciiLine(ifile) - a_coeffs = string.split(_line) - - x0 = float(a_coeffs[0]) - _line = fileutil.rAsciiLine(ifile) - a_coeffs[len(a_coeffs):] = string.split(_line) - # Scale coefficients for use within PyDrizzle - for i in range(len(a_coeffs)): - a_coeffs[i] = float(a_coeffs[i]) - - _line = fileutil.rAsciiLine(ifile) - b_coeffs = string.split(_line) - y0 = float(b_coeffs[0]) - _line = fileutil.rAsciiLine(ifile) - b_coeffs[len(b_coeffs):] = string.split(_line) - # Scale coefficients for use within PyDrizzle - for i in range(len(b_coeffs)): - b_coeffs[i] = float(b_coeffs[i]) - - ifile.close() - del ifile - # Now, convert the coefficients into a Numeric array - # with the right coefficients in the right place. - # Populate output values now... - fx = np.zeros(shape=(order+1,order+1),dtype=np.float64) - fy = np.zeros(shape=(order+1,order+1),dtype=np.float64) - # Assign the coefficients to their array positions - fx[0,0] = 0. - fx[1] = np.array([a_coeffs[2],a_coeffs[1],0.,0.],dtype=np.float64) - fx[2] = np.array([a_coeffs[5],a_coeffs[4],a_coeffs[3],0.],dtype=np.float64) - fx[3] = np.array([a_coeffs[9],a_coeffs[8],a_coeffs[7],a_coeffs[6]],dtype=np.float64) - fy[0,0] = 0. - fy[1] = np.array([b_coeffs[2],b_coeffs[1],0.,0.],dtype=np.float64) - fy[2] = np.array([b_coeffs[5],b_coeffs[4],b_coeffs[3],0.],dtype=np.float64) - fy[3] = np.array([b_coeffs[9],b_coeffs[8],b_coeffs[7],b_coeffs[6]],dtype=np.float64) - - # Used in Pattern.computeOffsets() - refpix = {} - refpix['XREF'] = None - refpix['YREF'] = None - refpix['V2REF'] = x0 - refpix['V3REF'] = y0 - refpix['XDELTA'] = 0. - refpix['YDELTA'] = 0. - refpix['PSCALE'] = None - refpix['DEFAULT_SCALE'] = no - refpix['centered'] = yes - - return fx,fy,refpix,order - -def factorial(n): - """ Compute a factorial for integer n. """ - m = 1 - for i in range(int(n)): - m = m * (i+1) - return m - -def combin(j,n): - """ Return the combinatorial factor for j in n.""" - return (factorial(j) / (factorial(n) * factorial( (j-n) ) ) ) - - -def defaultModel(): - """ This function returns a default, non-distorting model - that can be used with the data. - """ - order = 3 - - fx = np.zeros(shape=(order+1,order+1),dtype=np.float64) - fy = np.zeros(shape=(order+1,order+1),dtype=np.float64) - - fx[1,1] = 1. - fy[1,0] = 1. - - # Used in Pattern.computeOffsets() - refpix = {} - refpix['empty_model'] = yes - refpix['XREF'] = None - refpix['YREF'] = None - refpix['V2REF'] = 0. - refpix['XSIZE'] = 0. - refpix['YSIZE'] = 0. - refpix['V3REF'] = 0. - refpix['XDELTA'] = 0. - refpix['YDELTA'] = 0. - refpix['PSCALE'] = None - refpix['DEFAULT_SCALE'] = no - refpix['THETA'] = 0. - refpix['centered'] = yes - return fx,fy,refpix,order - -# Function to compute the index of refraction for MgF2 at -# the specified wavelength for use with Trauger coefficients -def _MgF2(lam): - _sig = pow((1.0e7/lam),2) - return np.sqrt(1.0 + 2.590355e10/(5.312993e10-_sig) + - 4.4543708e9/(11.17083e9-_sig) + 4.0838897e5/(1.766361e5-_sig)) - - -def convertDate(date): - """ Converts the DATE-OBS date string into an integer of the - number of seconds since 1970.0 using calendar.timegm(). - - INPUT: DATE-OBS in format of 'YYYY-MM-DD'. - OUTPUT: Date (integer) in seconds. - """ - - _dates = date.split('-') - _val = 0 - _date_tuple = (int(_dates[0]), int(_dates[1]), int(_dates[2]), 0, 0, 0, 0, 0, 0) - - return calendar.timegm(_date_tuple) diff --git a/distortion/utils.py b/distortion/utils.py deleted file mode 100644 index a99e227..0000000 --- a/distortion/utils.py +++ /dev/null @@ -1,204 +0,0 @@ -from __future__ import division # confidence high -import os -import numpy as np -import pywcs -import pyfits -from stwcs import wcsutil -from numpy import sqrt, arctan2 -from pytools import fileutil - -def output_wcs(list_of_wcsobj, ref_wcs=None, owcs=None, undistort=True): - """ - Create an output WCS. - - Parameters - ---------- - list_of_wcsobj: Python list - a list of HSTWCS objects - ref_wcs: an HSTWCS object - to be used as a reference WCS, in case outwcs is None. - if ref_wcs is None (default), the first member of the list - is used as a reference - outwcs: an HSTWCS object - the tangent plane defined by this object is used as a reference - undistort: boolean (default-True) - a flag whether to create an undistorted output WCS - """ - fra_dec = np.vstack([w.calcFootprint() for w in list_of_wcsobj]) - - # This new algorithm may not be strictly necessary, but it may be more - # robust in handling regions near the poles or at 0h RA. - crval1,crval2 = computeFootprintCenter(fra_dec) - - crval = np.array([crval1,crval2], dtype=np.float64) # this value is now zero-based - if owcs is None: - if ref_wcs == None: - ref_wcs = list_of_wcsobj[0].deepcopy() - if undistort: - outwcs = undistortWCS(ref_wcs) - else: - outwcs = ref_wcs.deepcopy() - outwcs.wcs.crval = crval - outwcs.wcs.set() - outwcs.pscale = sqrt(outwcs.wcs.cd[0,0]**2 + outwcs.wcs.cd[1,0]**2)*3600. - outwcs.orientat = arctan2(outwcs.wcs.cd[0,1],outwcs.wcs.cd[1,1]) * 180./np.pi - else: - outwcs = owcs.deepcopy() - outwcs.pscale = sqrt(outwcs.wcs.cd[0,0]**2 + outwcs.wcs.cd[1,0]**2)*3600. - outwcs.orientat = arctan2(outwcs.wcs.cd[0,1],outwcs.wcs.cd[1,1]) * 180./np.pi - - tanpix = outwcs.wcs.s2p(fra_dec, 0)['pixcrd'] - - outwcs.naxis1 = int(np.ceil(tanpix[:,0].max() - tanpix[:,0].min())) - outwcs.naxis2 = int(np.ceil(tanpix[:,1].max() - tanpix[:,1].min())) - crpix = np.array([outwcs.naxis1/2., outwcs.naxis2/2.], dtype=np.float64) - outwcs.wcs.crpix = crpix - outwcs.wcs.set() - tanpix = outwcs.wcs.s2p(fra_dec, 0)['pixcrd'] - - # shift crpix to take into account (floating-point value of) position of - # corner pixel relative to output frame size: no rounding necessary... - newcrpix = np.array([crpix[0]+tanpix[:,0].min(), crpix[1]+ - tanpix[:,1].min()]) - - newcrval = outwcs.wcs.p2s([newcrpix], 1)['world'][0] - outwcs.wcs.crval = newcrval - outwcs.wcs.set() - return outwcs - -def computeFootprintCenter(edges): - """ Geographic midpoint in spherical coords for points defined by footprints. - Algorithm derived from: http://www.geomidpoint.com/calculation.html - - This algorithm should be more robust against discontinuities at the poles. - """ - alpha = fileutil.DEGTORAD(edges[:,0]) - dec = fileutil.DEGTORAD(edges[:,1]) - - xmean = np.mean(np.cos(dec)*np.cos(alpha)) - ymean = np.mean(np.cos(dec)*np.sin(alpha)) - zmean = np.mean(np.sin(dec)) - - crval1 = fileutil.RADTODEG(np.arctan2(ymean,xmean))%360.0 - crval2 = fileutil.RADTODEG(np.arctan2(zmean,np.sqrt(xmean*xmean+ymean*ymean))) - - return crval1,crval2 - -def undistortWCS(wcsobj): - """ - Creates an undistorted linear WCS by applying the IDCTAB distortion model - to a 3-point square. The new ORIENTAT angle is calculated as well as the - plate scale in the undistorted frame. - """ - assert isinstance(wcsobj, pywcs.WCS) - import coeff_converter - - cx, cy = coeff_converter.sip2idc(wcsobj) - # cx, cy can be None because either there is no model available - # or updatewcs was not run. - if cx == None or cy == None: - if foundIDCTAB(wcsobj.idctab): - m = """IDCTAB is present but distortion model is missing. - Run updatewcs() to update the headers or - pass 'undistort=False' keyword to output_wcs().\n - """ - raise RuntimeError, m - else: - print 'Distortion model is not available, using input reference image for output WCS.\n' - return wcsobj.copy() - crpix1 = wcsobj.wcs.crpix[0] - crpix2 = wcsobj.wcs.crpix[1] - xy = np.array([(crpix1,crpix2),(crpix1+1.,crpix2),(crpix1,crpix2+1.)],dtype=np.double) - offsets = np.array([wcsobj.ltv1, wcsobj.ltv2]) - px = xy + offsets - #order = wcsobj.sip.a_order - pscale = wcsobj.idcscale - #pixref = np.array([wcsobj.sip.SIPREF1, wcsobj.sip.SIPREF2]) - - tan_pix = apply_idc(px, cx, cy, wcsobj.wcs.crpix, pscale, order=1) - xc = tan_pix[:,0] - yc = tan_pix[:,1] - am = xc[1] - xc[0] - bm = xc[2] - xc[0] - cm = yc[1] - yc[0] - dm = yc[2] - yc[0] - cd_mat = np.array([[am,bm],[cm,dm]],dtype=np.double) - - # Check the determinant for singularity - _det = (am * dm) - (bm * cm) - if ( _det == 0.0): - print 'Singular matrix in updateWCS, aborting ...' - return - - lin_wcsobj = pywcs.WCS() - cd_inv = np.linalg.inv(cd_mat) - cd = np.dot(wcsobj.wcs.cd, cd_inv).astype(np.float64) - lin_wcsobj.wcs.cd = cd - lin_wcsobj.wcs.set() - lin_wcsobj.orientat = arctan2(lin_wcsobj.wcs.cd[0,1],lin_wcsobj.wcs.cd[1,1]) * 180./np.pi - lin_wcsobj.pscale = sqrt(lin_wcsobj.wcs.cd[0,0]**2 + lin_wcsobj.wcs.cd[1,0]**2)*3600. - lin_wcsobj.wcs.crval = np.array([0.,0.]) - lin_wcsobj.wcs.crpix = np.array([0.,0.]) - lin_wcsobj.wcs.ctype = ['RA---TAN', 'DEC--TAN'] - lin_wcsobj.wcs.set() - return lin_wcsobj - -def apply_idc(pixpos, cx, cy, pixref, pscale= None, order=None): - """ - Apply the IDCTAB polynomial distortion model to pixel positions. - pixpos must be already corrected for ltv1/2. - - Parameters - ---------- - pixpos: a 2D numpy array of (x,y) pixel positions to be distortion corrected - cx, cy: IDC model distortion coefficients - pixref: reference opixel position - - """ - if cx == None: - return pixpos - - if order is None: - print 'Unknown order of distortion model \n' - return pixpos - if pscale is None: - print 'Unknown model plate scale\n' - return pixpos - - # Apply in the same way that 'drizzle' would... - _cx = cx/pscale - _cy = cy/ pscale - _p = pixpos - - # Do NOT include any zero-point terms in CX,CY here - # as they should not be scaled by plate-scale like rest - # of coeffs... This makes the computations consistent - # with 'drizzle'. WJH 17-Feb-2004 - _cx[0,0] = 0. - _cy[0,0] = 0. - - dxy = _p - pixref - # Apply coefficients from distortion model here... - - c = _p * 0. - for i in range(order+1): - for j in range(i+1): - c[:,0] = c[:,0] + _cx[i][j] * pow(dxy[:,0],j) * pow(dxy[:,1],(i-j)) - c[:,1] = c[:,1] + _cy[i][j] * pow(dxy[:,0],j) * pow(dxy[:,1],(i-j)) - - return c - -def foundIDCTAB(idctab): - idctab_found = True - try: - idctab = fileutil.osfn(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 - |