diff options
| -rw-r--r-- | distortion/__init__.py | 0 | ||||
| -rw-r--r-- | distortion/models.py | 360 | ||||
| -rw-r--r-- | distortion/mutil.py | 605 | ||||
| -rw-r--r-- | lib/__init__.py | 47 | ||||
| -rw-r--r-- | lib/mappings.py | 38 | ||||
| -rw-r--r-- | setup.py | 23 | ||||
| -rw-r--r-- | setup.py.old | 89 | ||||
| -rw-r--r-- | updatewcs/__init__.py | 262 | ||||
| -rw-r--r-- | updatewcs/corrections.py | 164 | ||||
| -rw-r--r-- | updatewcs/dgeo.py | 206 | ||||
| -rw-r--r-- | updatewcs/makewcs.py | 260 | ||||
| -rw-r--r-- | wcsutil/__init__.py | 166 | ||||
| -rw-r--r-- | wcsutil/instruments.py | 117 |
13 files changed, 2337 insertions, 0 deletions
diff --git a/distortion/__init__.py b/distortion/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/distortion/__init__.py diff --git a/distortion/models.py b/distortion/models.py new file mode 100644 index 0000000..cff9a46 --- /dev/null +++ b/distortion/models.py @@ -0,0 +1,360 @@ +import types + +# Import PyDrizzle utility modules +import mutil +import numpy as N +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 = N.zeros(shape=cx.shape,dtype=cx.dtype) + _cys = N.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 = N.array(_p,dtype=N.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 new file mode 100644 index 0000000..c3f83f2 --- /dev/null +++ b/distortion/mutil.py @@ -0,0 +1,605 @@ +from pytools import fileutil +import numpy as N +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 + + # 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 = N.zeros(shape=(order+1,order+1),dtype=N.float64) + fy = N.zeros(shape=(order+1,order+1),dtype=N.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 + 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: '+filter1+','+filter2+'\n' + ftab.close() + del ftab + raise LookupError,err_str + else: + print '- IDCTAB: Distortion model from row',str(row+1),'for chip',detchip,':',filter1.strip(),'and',filter2.strip() + + # 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 + + # 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 + +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 = N.zeros(shape=(order+1,order+1),dtype=N.float64) + fy = N.zeros(shape=(order+1,order+1),dtype=N.float64) + + # Read in CD matrix + _cd11 = header['cd1_1'] + _cd12 = header['cd1_2'] + _cd21 = header['cd2_1'] + _cd22 = header['cd2_2'] + _cdmat = N.array([[_cd11,_cd12],[_cd21,_cd22]]) + _theta = N.arctan2(-_cd12,_cd22) + _rotmat = N.array([[N.cos(_theta),N.sin(_theta)], + [-N.sin(_theta),N.cos(_theta)]]) + _rCD = N.dot(_rotmat,_cdmat) + _skew = N.arcsin(-_rCD[1][0] / _rCD[0][0]) + _scale = _rCD[0][0] * N.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 = N.zeros(shape=(order+1,order+1),dtype=N.float64) + fy = N.zeros(shape=(order+1,order+1),dtype=N.float64) + # Assign the coefficients to their array positions + fx[0,0] = 0. + fx[1] = N.array([a_coeffs[2],a_coeffs[1],0.,0.],dtype=N.float64) + fx[2] = N.array([a_coeffs[5],a_coeffs[4],a_coeffs[3],0.],dtype=N.float64) + fx[3] = N.array([a_coeffs[9],a_coeffs[8],a_coeffs[7],a_coeffs[6]],dtype=N.float64) + fy[0,0] = 0. + fy[1] = N.array([b_coeffs[2],b_coeffs[1],0.,0.],dtype=N.float64) + fy[2] = N.array([b_coeffs[5],b_coeffs[4],b_coeffs[3],0.],dtype=N.float64) + fy[3] = N.array([b_coeffs[9],b_coeffs[8],b_coeffs[7],b_coeffs[6]],dtype=N.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 = N.zeros(shape=(order+1,order+1),dtype=N.float64) + fy = N.zeros(shape=(order+1,order+1),dtype=N.float64) + # Assign the coefficients to their array positions + fx[0,0] = 0. + fx[1] = N.array([a_coeffs[2],a_coeffs[1],0.,0.],dtype=N.float64) + fx[2] = N.array([a_coeffs[5],a_coeffs[4],a_coeffs[3],0.],dtype=N.float64) + fx[3] = N.array([a_coeffs[9],a_coeffs[8],a_coeffs[7],a_coeffs[6]],dtype=N.float64) + fy[0,0] = 0. + fy[1] = N.array([b_coeffs[2],b_coeffs[1],0.,0.],dtype=N.float64) + fy[2] = N.array([b_coeffs[5],b_coeffs[4],b_coeffs[3],0.],dtype=N.float64) + fy[3] = N.array([b_coeffs[9],b_coeffs[8],b_coeffs[7],b_coeffs[6]],dtype=N.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 = N.zeros(shape=(order+1,order+1),dtype=N.float64) + fy = N.zeros(shape=(order+1,order+1),dtype=N.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 N.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/lib/__init__.py b/lib/__init__.py new file mode 100644 index 0000000..08e421b --- /dev/null +++ b/lib/__init__.py @@ -0,0 +1,47 @@ +import wcsutil +from wcsutil import HSTWCS +from pytools import fileutil, parseinput +import pyfits + +#to avoid relative imports +import mappings +from mappings import basic_wcs +import distortion +import pywcs + +def restoreWCS(fnames): + """ + Given a list of fits file names restore the original basic WCS kw + and write out the files. This overwrites the original files. + """ + files = parseinput.parseinput(fnames)[0] + for f in files: + isfits, ftype = fileutil.isFits(f) + if not isfits or (isfits and ftype == 'waiver'): + print "RestoreWCS works only on multiextension fits files." + return + else: + fobj = pyfits.open(f, mode='update') + hdr0 = fobj[0].header + for ext in fobj: + try: + extname = ext.header['EXTNAME'].lower() + except KeyError: + continue + if extname in ['sci', 'err', 'sdq']: + hdr = ext.header + owcs = HSTWCS(hdr0, hdr) + #Changed restore to update the attributes ??? + # this may need to be changed + backup = owcs.restore() + if not backup: + print '\No archived keywords found.\n' + continue + else: + for kw in basic_wcs: + nkw = ('O'+kw)[:7] + if backup.has_key(kw): + hdr.update(kw, hdr[nkw]) + fobj.close() + + diff --git a/lib/mappings.py b/lib/mappings.py new file mode 100644 index 0000000..fed9b60 --- /dev/null +++ b/lib/mappings.py @@ -0,0 +1,38 @@ +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'], + 'ACS': ['TDDCorr', 'MakeWCS','CompSIP', 'VACorr'] + } + +# 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'] + +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/setup.py b/setup.py new file mode 100644 index 0000000..cc150f1 --- /dev/null +++ b/setup.py @@ -0,0 +1,23 @@ +from distutils.core import setup +from os.path import join +# PyFITS +try: + import pyfits +except ImportError: + print "WARNING: PyFITS must be installed to use hstwcs." + print " Since this is not a build-time dependency, the build will proceed." + +# PyWCS +try: + import pywcs +except ImportError: + print "WARNING: PyWCS must be installed to use hstwcs." + print " Since this is not a build-time dependency, the build will proceed." + +setup(name="hstwcs", + version="0.1", + description="HST WCS Corrections", + packages=['hstwcs', 'hstwcs/updatewcs', 'hstwcs/wcsutil', 'hstwcs/distortion'], + package_dir={'hstwcs':'lib', 'hstwcs/updatewcs': 'updatewcs', + 'hstwcs/wcsutil': 'wcsutil', 'hstwcs/distortion': 'distortion'} + ) diff --git a/setup.py.old b/setup.py.old new file mode 100644 index 0000000..a4a57a2 --- /dev/null +++ b/setup.py.old @@ -0,0 +1,89 @@ +from distutils.core import setup, Extension +from os.path import join + +###################################################################### +# NUMPY +try: + import numpy +except ImportError: + print "numpy must be installed to build pywcs." + print "ABORTING." + sys.exit(1) + +try: + numpy_include = numpy.get_include() +except AttributeError: + numpy_include = numpy.get_numpy_include() + +###################################################################### +# PyFITS +try: + import pyfits +except ImportError: + print "WARNING: PyFITS must be installed to use pywcs." + print " Since this is not a build-time dependency, the build will proceed." + +###################################################################### +###################################################################### +# WCSLIB +WCSVERSION = "4.3" +WCSLIB = "wcslib-%s" % WCSVERSION # Path to wcslib +WCSLIBC = join('pywcs', WCSLIB, "C") # Path to wcslib source files +WCSFILES = [ # List of wcslib files to compile + 'flexed/wcsulex.c', + 'flexed/wcsutrn.c', + 'cel.c', + 'lin.c', + 'log.c', + 'prj.c', + 'spc.c', + 'sph.c', + 'spx.c', + 'tab.c', + 'wcs.c', + 'wcsfix.c', + 'wcshdr.c', + 'wcspih.c', + 'wcstrig.c', + 'wcsunits.c', + 'wcsutil.c'] +WCSFILES = [join(WCSLIBC, x) for x in WCSFILES] +###################################################################### + +###################################################################### +# GENERATE DOCSTRINGS IN C +docstrings = {} +execfile("pywcs/doc/docstrings.py", docstrings) +keys = docstrings.keys() +keys.sort() +fd = open("pywcs/src/docstrings.h", "w") +fd.write('/* This file is autogenerated by setup.py. To edit its contents\n') +fd.write(' edit doc/docstrings.py\n') +fd.write('*/\n') +for key in keys: + if key.startswith('__'): + continue + val = docstrings[key].lstrip().encode("string_escape").replace('"', '\\"') + fd.write("static const char doc_%s[] = \"%s\";\n\n" % (key, val)) +fd.close() + +###################################################################### + +setup(name="hstwcs", + version="0.1", + description="HST WCS Corrections", + packages=['hstwcs', 'hstwcs/pywcs', 'hstwcs/updatewcs', 'hstwcs/wcsutil', 'hstwcs/distortion'], + package_dir={'hstwcs':'lib', 'hstwcs/pywcs': 'pywcs/pywcs', 'hstwcs/updatewcs': 'updatewcs', + 'hstwcs/wcsutil': 'wcsutil', 'hstwcs/distortion': 'distortion'}, + ext_modules=[ + Extension('hstwcs.pywcs._pywcs', + ['pywcs/src/pywcs.c'] + + WCSFILES, + include_dirs=[ + numpy_include, + WCSLIBC, + "src" + ] + ) + ] + ) diff --git a/updatewcs/__init__.py b/updatewcs/__init__.py new file mode 100644 index 0000000..db89d13 --- /dev/null +++ b/updatewcs/__init__.py @@ -0,0 +1,262 @@ +import os +import pyfits +#from .. wcsutil import HSTWCS +from hstwcs.wcsutil import HSTWCS +from hstwcs.mappings import allowed_corrections +#from .. mappings import allowed_corrections +import corrections, makewcs +import dgeo +import time +from pytools import parseinput, fileutil + +#NB! the order of corrections matters + +__docformat__ = 'restructuredtext' + + +def updatewcs(input, vacorr=True, tddcorr=True, checkfiles=True): + """ + Purpose + ======= + Applies corrections to the WCS keywords. + + Example + ======= + >>>from hstwcs import updatewcs + >>>updatewcs.updatewcs(filename) + + Dependencies + ============ + `pytools` + `pyfits` + `pywcs` + `numpy` + + :Parameters: + `input`: a python list of file names or a string (wild card characters allowed) + input files may be in fits, geis or waiver fits format + `vacorr`: boolean + If True, vecocity aberration correction will be applied + `tddcorr`: boolean + If True, time dependen correction will be applied to the distortion model + `checkfiles`: boolean + If True, the format of the input files will be checked, + geis and waiver fits files will be converted to MEF format. + Default value is True for standalone mode. + """ + + files = parseinput.parseinput(input)[0] + if checkfiles: + files = checkFiles(files) + if not files: + print 'No valid input, quitting ...\n' + return + for f in files: + instr = pyfits.getval(f, 'INSTRUME') + try: + acorr = setCorrections(instr,vacorr=vacorr, tddcorr=tddcorr) + except KeyError: + print 'Unsupported instrument %s ' %instr + print 'Removing %s from list of processed files\n' % f + files.remove(f) + continue + + makecorr(f, acorr) + return files + +def makecorr(fname, acorr): + """ + Purpose + ======= + Applies corrections to a single file + + :Parameters: + `fname`: string + file name + `acorr`: list + list of corrections to be applied + + """ + f = pyfits.open(fname, mode='update') + nrefchip, nrefext = getNrefchip(f) + primhdr = f[0].header + + + for extn in f: + refwcs = HSTWCS(primhdr, f[nrefext].header) + refwcs.archive_kw() + refwcs.readModel() + if extn.header.has_key('extname') and extn.header['extname'].lower() == 'sci': + hdr = extn.header + owcs = HSTWCS(primhdr, hdr) + owcs.archive_kw() + owcs.readModel() + for c in acorr: + owcs.__setattr__('DO'+c, 'PERFORM') + corr_klass = corrections.__getattribute__(c) + corr_klass(owcs, refwcs) + + + + #always do dgeo correction + if applyDgeoCorr(fname): + dgeo.DGEO(f) + f.close() + + +def setCorrections(instrument, vacorr=True, tddcorr=True): + """ + Purpose + ======= + 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. + """ + 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') + + return acorr + + + +def applyDgeoCorr(fname): + """ + Purpose + ======= + Adds dgeo extensions to files based on the DGEOFILE keyword in the primary + header. This is a default correction and will always run in the pipeline. + The file used to generate the extensions is + recorded in the DGEOFILE keyword in each science extension. + If 'DGEOFILE' in the primary header is different from 'DGEOFILE' in the + extension header and the file exists on disk and is a 'new type' dgeofile, + then the dgeo extensions will be updated. + """ + applyDGEOCorr = True + try: + # get DGEOFILE kw from primary header + fdgeo0 = pyfits.getval(fname, 'DGEOFILE') + fdgeo0 = fileutil.osfn(fdgeo0) + if not fileutil.findFile(fdgeo0): + print 'Kw DGEOFILE exists in primary header but file %s not found\n' % fdgeo0 + print 'DGEO correction will not be applied\n' + applyDGEOCorr = False + return applyDGEOCorr + try: + # get DGEOFILE kw from first extension header + fdgeo1 = pyfits.getval(fname, 'DGEOFILE', ext=1) + fdgeo1 = fileutil.osfn(fdgeo1) + if fdgeo1 and fileutil.findFile(fdgeo1): + if fdgeo0 != fdgeo1: + applyDGEOCorr = True + else: + applyDGEOCorr = False + else: + # dgeo file defined in first extension may not be found + # but if a valid kw exists in the primary header, dgeo should be applied. + applyDGEOCorr = True + except KeyError: + # the case of DGEOFILE kw present in primary header but missing + # in first extension header + applyDGEOCorr = True + except KeyError: + + print 'DGEOFILE keyword not found in primary header' + applyDGEOCorr = False + + if isOldStyleDGEO(fname, fdgeo0): + applyDGEOCorr = False + + return applyDGEOCorr + +def isOldStyleDGEO(fname, dgname): + # checks if the file defined in a DGEOFILE kw is a full size + # (old style) image + + sci_naxis1 = pyfits.getval(fname, 'NAXIS1', ext=1) + sci_naxis2 = pyfits.getval(fname, 'NAXIS2', ext=1) + dg_naxis1 = pyfits.getval(dgname, 'NAXIS1', ext=1) + dg_naxis2 = pyfits.getval(dgname, 'NAXIS2', ext=1) + if sci_naxis1 <= dg_naxis1 or sci_naxis2 <= dg_naxis2: + print 'Only full size (old style) XY file was found.' + print 'DGEO correction will not be applied.\n' + return True + else: + return False + +def getNrefchip(fobj): + """ + This handles the fact that WFPC2 subarray observations + may not include chip 3 which is the default reference chip for + full observations. Also for subarrays chip 3 may not be the third + extension in a MEF file. + """ + Nrefext = 1 + instrument = fobj[0].header['INSTRUME'] + if instrument == 'WFPC2': + detectors = [img.header['DETECTOR'] for img in fobj[1:]] + + if 3 not in detectors: + Nrefchip=detectors[0] + Nrefext = 1 + else: + Nrefchip = 3 + Nrefext = detectors.index(3) + 1 + elif instrument == 'ACS': + detector = fobj[0].header['DETECTOR'] + if detector == 'WCS': + Nrefchip =2 + else: + Nrefchip = 1 + elif instrument == 'NICMOS': + Nrefchip = fobj[0].header['CAMERA'] + return Nrefchip, Nrefext + +def checkFiles(input): + """ + Purpose + ======= + Checks that input files are in the correct format. + Converts geis and waiver fits files to multietension fits. + """ + from pytools.check_files import geis2mef, waiver2mef + removed_files = [] + newfiles = [] + for file in input: + try: + imgfits,imgtype = fileutil.isFits(file) + except IOError: + print "Warning: File %s could not be found\n" %file + print "Removing file %s from input list" %file + removed_files.append(file) + continue + # Check for existence of waiver FITS input, and quit if found. + # Or should we print a warning and continue but not use that file + if imgfits: + if imgtype == 'waiver': + newfilename = waiver2mef(file, convert_dq=True) + if newfilename == None: + print "Removing file %s from input list - could not convert waiver to mef" %file + removed_files.append(file) + else: + newfiles.append(newfilename) + else: + newfiles.append(file) + + # If a GEIS image is provided as input, create a new MEF file with + # a name generated using 'buildFITSName()' + # Convert the corresponding data quality file if present + if not imgfits: + newfilename = geis2mef(file, convert_dq=True) + if newfilename == None: + print "Removing file %s from input list - could not convert geis to mef" %file + removed_files.append(file) + else: + newfiles.append(newfilename) + if removed_files: + print 'The following files will be removed from the list of files to be processed :\n' + for f in removed_files: + print f + return newfiles + diff --git a/updatewcs/corrections.py b/updatewcs/corrections.py new file mode 100644 index 0000000..800c1c6 --- /dev/null +++ b/updatewcs/corrections.py @@ -0,0 +1,164 @@ +import datetime +import numpy +from numpy import linalg +#from pytools import wcsutil + +from makewcs import MakeWCS +from dgeo import DGEO + +class TDDCorr(object): + """ + Purpose + ======= + Apply time dependent distortion correction to SIP coefficients and basic + WCS keywords. Atthi stime it is applicable only to ACS/WFC data. + + Ref: Jay Anderson, ACS ISR 2007-08, Variation of the Distortion + Solution of the WFC + + """ + def __init__(self, owcs=None, refwcs=None): + """ + :Parameters: + `owcs`: HSTWCS object + An extension HSTWCS object to be modified + `refwcs`: HSTWCS object + A reference HSTWCS object + """ + self.wcsobj = owcs + if self.wcsobj.DOTDDCorr == 'PERFORM': + self.updateWCS() + self.wcsobj.DOTDDCorr = 'COMPLETE' + else: + pass + + def updateWCS(self): + """ + - Calculates alpha and beta for ACS/WFC data. + - Writes 2 new kw to the extension header: TDDALPHA and TDDBETA + """ + if self.wcsobj.instrument == 'ACS' and self.wcsobj.detector == 'WFC': + newkw = ['TDDALPHA', 'TDDBETA'] + self.set_alpha_beta() + self.updatehdr(newkeywords=newkw) + + else: + pass + + def set_alpha_beta(self): + # Compute the time dependent distortion skew terms + # default date of 2004.5 = 2004-7-1 + + datedefault = datetime.datetime(2004,7,1) + year,month,day = self.wcsobj.date_obs.split('-') + rdate = datetime.datetime(int(year),int(month),int(day)) + self.TDDALPHA = 0.095+0.090*((rdate-datedefault).days/365.25)/2.5 + self.TDDBETA = -0.029-0.030*((rdate-datedefault).days/365.25)/2.5 + self.TDDALPHA = 0.0 + self.TDDBETA = 0.0 + + def updatehdr(self, newkeywords=None): + + for kw in newkeywords: + self.wcsobj.hdr.update(kw, self.__getattribute__(kw)) + + +class VACorr(object): + """ + Purpose + ======= + Apply velocity aberation correction to WCS keywords. + Modifies the CD matrix and CRVAL1/2 + """ + def __init__(self, owcs=None, refwcs=None): + self.wcsobj = owcs + self.refwcs = refwcs + if self.wcsobj.DOVACorr == 'PERFORM': + self.updateWCS() + self.wcsobj.DOVACorr == 'COMPLETE' + else: + pass + + def updateWCS(self): + if self.wcsobj.vafactor != 1: + self.wcsobj.cd = self.wcsobj.cd * self.wcsobj.vafactor + + #self.wcsobj.crval[0] = self.wcsobj.ra_targ + self.wcsobj.vafactor*(self.wcsobj.crval[0] - self.wcsobj.ra_targ) + #self.wcsobj.crval[1] = self.wcsobj.dec_targ + self.wcsobj.vafactor*(self.wcsobj.crval[1] - self.wcsobj.dec_targ) + ref_backup = self.refwcs.restore() + crval1 = ref_backup['CRVAL1'] + crval2 = ref_backup['CRVAL2'] + self.wcsobj.crval[0] = crval1 + self.wcsobj.vafactor*diff_angles(self.wcsobj.crval[0], crval1) + self.wcsobj.crval[1] = crval2 + self.wcsobj.vafactor*diff_angles(self.wcsobj.crval[1], crval2) + + kw2update={'CD1_1': self.wcsobj.cd[0,0], 'CD1_2':self.wcsobj.cd[0,1], + 'CD2_1':self.wcsobj.cd[1,0], 'CD2_2':self.wcsobj.cd[1,1], + 'CRVAL1':self.wcsobj.crval[0], 'CRVAL2':self.wcsobj.crval[1]} + self.updatehdr(newkeywords=kw2update) + else: + pass + + def updatehdr(self, newkeywords=None): + for kw in newkeywords: + self.wcsobj.hdr.update(kw, newkeywords[kw]) + +class CompSIP(object): + """ + Purpose + ======= + Compute SIP coefficients from idc table coefficients. + + """ + def __init__(self, owcs, refwcs): + self.wcsobj = owcs + self.updateWCS() + self.DOCOMPSIP = 'COMPLETE' + + def updateWCS(self): + kw2update = {} + order = self.wcsobj.idcmodel.norder + kw2update['A_ORDER'] = order + kw2update['B_ORDER'] = order + pscale = self.wcsobj.idcmodel.refpix['PSCALE'] + + cx = self.wcsobj.idcmodel.cx + cy = self.wcsobj.idcmodel.cy + + matr = numpy.array([[cx[1,1],cx[1,0]], [cy[1,1],cy[1,0]]], dtype=numpy.float) + imatr = linalg.inv(matr) + akeys1 = numpy.zeros((order+1,order+1), dtype=numpy.float) + bkeys1 = numpy.zeros((order+1,order+1), dtype=numpy.float) + for n in range(order+1): + for m in range(order+1): + if n >= m and n>=2: + idcval = numpy.array([[cx[n,m]],[cy[n,m]]]) + sipval = numpy.dot(imatr, idcval) + akeys1[m,n-m] = sipval[0] + bkeys1[m,n-m] = sipval[1] + Akey="A_%d_%d" % (m,n-m) + Bkey="B_%d_%d" % (m,n-m) + kw2update[Akey] = sipval[0,0] + kw2update[Bkey] = sipval[1,0] + self.updatehdr(newkw=kw2update) + + def updatehdr(self, newkw=None): + if not newkw: return + for kw in newkw.keys(): + self.wcsobj.hdr.update(kw, newkw[kw]) + + +def diff_angles(a,b): + """ + Perform angle subtraction a-b taking into account + small-angle differences across 360degree line. + """ + + diff = a - b + + if diff > 180.0: + diff -= 360.0 + + if diff < -180.0: + diff += 360.0 + + return diff diff --git a/updatewcs/dgeo.py b/updatewcs/dgeo.py new file mode 100644 index 0000000..32e1a38 --- /dev/null +++ b/updatewcs/dgeo.py @@ -0,0 +1,206 @@ +import pyfits +from pytools import fileutil +from hstwcs.mappings import dgeo_vals + +class DGEO(object): + """ + Purpose + ======= + Defines a Lookup table prior distortion correction as per WCS paper IV. + It uses a reference file defined by the DGEOFILE keyword in the primary header. + + Algorithm + ========= + - 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 + science extension header + - Add a keyword 'DGEOFILE' to the science extension header, whose + value is the reference file used to create the WCSVARR extension + + """ + def __init__(self, fobj): + """ + :Parameters: + `fobj`: pyfits object + Science file, for which a distortion correc tion in a DGEOFILE is available + + """ + assert isinstance(fobj, pyfits.NP_pyfits.HDUList) + self.fobj = fobj + self.applyDgeoCorr() + + + def applyDgeoCorr(self): + """ + 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(self.fobj[0].header['DGEOFILE']) + wcsdvarr_ind = self.getwcsindex() + for ext in self.fobj: + try: + extname = ext.header['EXTNAME'].lower() + except KeyError: + continue + if extname == 'sci': + extversion = ext.header['EXTVER'] + for ename in ['DX', 'DY']: + dgeover +=1 + + self.addSciExtKw(ext.header, extver=dgeover, ename=ename) + hdu = self.createDgeoHDU(dgeofile=dgfile, dgeover=dgeover,ename=ename, extver=extversion) + if wcsdvarr_ind: + self.fobj[wcsdvarr_ind[dgeover]] = hdu + else: + self.fobj.append(hdu) + self.updateDGEOkw(ext.header) + + def getwcsindex(self): + """ + Returns the index of a WCSDVARR extension in a pyfits file object if it exists. + If it exists subsequent updates will overwrite it. If not, it will be + added to the file object. + """ + wcsd = {} + for e in range(len(self.fobj)): + try: + ename = self.fobj[e].header['EXTNAME'] + except KeyError: + continue + if ename == 'WCSDVARR': + wcsd[self.fobj[e].header['EXTVER']] = e + return wcsd + + return self.fobj.index_of(('WCSDVARR', dgeover)) + + def addSciExtKw(self, hdr, extver=None, ename=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) + """ + if ename =='DX': + j=1 + else: + j=2 + + cperror = 'CPERROR%s' %j + cpdis = 'CPDIS%s' %j + dpext = 'DP%s.' %j + 'EXTVER' + dpnaxes = 'DP%s.' %j +'NAXES' + 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, + dpaxis1: 1, dpaxis2: 2} + + comments = {cperror: 'Maximum error of dgeo correction for axis %s' % (extver/2), + cpdis: 'Prior distortion funcion type', + dpext: 'Version number of WCSDVARR extension containing lookup distortion table', + dpnaxes: 'Number of independent variables in distortion function', + dpaxis1: 'Axis number of the jth independent variable in a distortion function', + dpaxis2: 'Axis number of the jth independent variable in a distortion function' + } + + for key in keys: + hdr.update(key=key, value=values[key], comment=comments[key], before='HISTORY') + + dgfile = self.fobj[0].header['DGEOFILE'] + hdr.update(key='DGEOFILE', value=dgfile, comment='DGEOFILE used to create WCSDVARR extensions', before='HISTORY') + + + def getDgeoData(self, dgfile=None, ename=None, extver=1): + """ + 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)) + + def createDgeoHDU(self, dgeofile=None, dgeover=1, ename=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 = self.createDgeoHdr(**dgeokw) + data = self.getDgeoData(dgfile=dgeofile, ename=ename, extver=extver) + hdu=pyfits.ImageHDU(header=hdr, data=data) + return hdu + + def createDgeoHdr(self, **kw): + """ + Creates a header for the dgeo extension based on dgeo 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'] + + keys = ['XTENSION','BITPIX','NAXIS','NAXIS1','NAXIS2', + 'EXTNAME','EXTVER','PCOUNT','GCOUNT','CRPIX1', + 'CDELT1','CRVAL1','CRPIX2','CDELT2','CRVAL2'] + + comments = {'XTENSION': 'Image extension', + 'BITPIX': 'IEEE floating point', + 'NAXIS': 'Number of image columns', + 'NAXIS1': 'Number of image columns', + 'NAXIS2': 'Number of image rows', + 'EXTNAME': 'WCS distortion array', + 'EXTVER': 'Distortion array version number', + 'PCOUNT': 'Special data area of size 0', + 'GCOUNT': 'One data group', + 'CRPIX1': 'Distortion array reference pixel', + 'CDELT1': 'Grid step size in first coordinate', + 'CRVAL1': 'Image array pixel coordinate', + 'CRPIX2': 'Distortion array reference pixel', + 'CDELT2': 'Grid step size in second coordinate', + 'CRVAL2': 'Image array pixel coordinate'} + + values = {'XTENSION': 'IMAGE', + 'BITPIX': -32, + 'NAXIS': 2, + 'NAXIS1': naxis1, + 'NAXIS2': naxis2, + 'EXTNAME': 'WCSDVARR', + 'EXTVER': extver, + 'PCOUNT': 0, + 'GCOUNT': 1, + 'CRPIX1': crpix1, + 'CDELT1': cdelt1, + 'CRVAL1': crval1, + 'CRPIX2': crpix1, + 'CDELT2': cdelt2, + 'CRVAL2': crval2 + } + + + cdl = pyfits.CardList() + for c in keys: + cdl.append(pyfits.Card(key=c, value=values[c], comment=comments[c])) + + hdr = pyfits.Header(cards=cdl) + return hdr + + def updateDGEOkw(self, hdr): + dgfile = self.fobj[0].header['DGEOFILE'] + hdr.update(key='DGEOFILE', value=dgfile, comment='DGEOFILE used to create WCSDVARR extensions', before='HISTORY') +
\ No newline at end of file diff --git a/updatewcs/makewcs.py b/updatewcs/makewcs.py new file mode 100644 index 0000000..451da13 --- /dev/null +++ b/updatewcs/makewcs.py @@ -0,0 +1,260 @@ +#from .. mappings import DEGTORAD, RADTODEG +from hstwcs.mappings import DEGTORAD, RADTODEG +import numpy +from numpy import math +from math import sin, sqrt, pow, cos, asin, atan2,pi + +class MakeWCS(object): + """ + Purpose + ======= + Recompute WCS keywords based on PA_V3 and distortion model. + + Algorithm + ========= + - update reference chip wcs + - update extension wcs + - update extension header + + """ + + def __init__(self, owcs, refwcs): + self.wcsobj = owcs + self.refwcs = refwcs + self.updateWCS() + self.DOMakeWCS = 'COMPLETE' + + def updateWCS(self): + """ + recomputes the basic WCS kw + """ + backup_wcs = self.wcsobj.restore() + ltvoff, offshift = self.getOffsets(backup_wcs) + + self.uprefwcs(ltvoff, offshift) + self.upextwcs(ltvoff, offshift) + self.updatehdr() + + def upextwcs(self, loff, offsh): + """ + updates an extension wcs + """ + ltvoffx, ltvoffy = loff[0], loff[1] + offshiftx, offshifty = offsh[0], offsh[1] + backup_wcs = self.wcsobj.restore() + ltv1 = self.wcsobj.ltv1 + ltv2 = self.wcsobj.ltv2 + + + + if ltv1 != 0. or ltv2 != 0.: + offsetx = backup_wcs['CRPIX1'] - ltv1 - self.wcsobj.idcmodel.refpix['XREF'] + offsety = backup_wcs['CRPIX2'] - ltv2 - self.wcsobj.idcmodel.refpix['YREF'] + fx,fy = self.wcsobj.idcmodel.shift(self.wcsobj.idcmodel.cx,self.wcsobj.idcmodel.cy,offsetx,offsety) + else: + fx, fy = self.wcsobj.idcmodel.cx, self.wcsobj.idcmodel.cy + + ridcmodel = self.refwcs.idcmodel + v2 = self.wcsobj.idcmodel.refpix['V2REF'] + v3 = self.wcsobj.idcmodel.refpix['V3REF'] + v2ref = self.refwcs.idcmodel.refpix['V2REF'] + + v3ref = self.refwcs.idcmodel.refpix['V3REF'] + R_scale = self.refwcs.idcmodel.refpix['PSCALE']/3600.0 + off = sqrt((v2-v2ref)**2 + (v3-v3ref)**2)/(R_scale*3600.0) + + # Here we must include the PARITY + if v3 == v3ref: + theta=0.0 + else: + theta = atan2(self.wcsobj.parity[0][0]*(v2-v2ref),self.wcsobj.parity[1][1]*(v3-v3ref)) + + if self.refwcs.idcmodel.refpix['THETA']: theta += self.refwcs.idcmodel.refpix['THETA']*pi/180.0 + + dX=(off*sin(theta)) + offshiftx + dY=(off*cos(theta)) + offshifty + px = numpy.array([[dX,dY]]) + newcrval = self.refwcs.pixel2world_fits(px)[0] + newcrpix = numpy.array([self.wcsobj.idcmodel.refpix['XREF'] + ltvoffx, + self.wcsobj.idcmodel.refpix['YREF'] + ltvoffy]) + + self.wcsobj.crval = newcrval + self.wcsobj.crpix = newcrpix + + # Account for subarray offset + # Angle of chip relative to chip + if self.wcsobj.idcmodel.refpix['THETA']: + dtheta = self.wcsobj.idcmodel.refpix['THETA'] - self.refwcs.idcmodel.refpix['THETA'] + else: + dtheta = 0.0 + + + # Create a small vector, in reference image pixel scale + # There is no parity effect here ??? + + delXX=fx[1,1]/R_scale/3600. + delYX=fy[1,1]/R_scale/3600. + delXY=fx[1,0]/R_scale/3600. + delYY=fy[1,0]/R_scale/3600. + + + # Convert to radians + rr=dtheta*pi/180.0 + + # Rotate the vectors + dXX = cos(rr)*delXX - sin(rr)*delYX + dYX = sin(rr)*delXX + cos(rr)*delYX + + dXY = cos(rr)*delXY - sin(rr)*delYY + dYY = sin(rr)*delXY + cos(rr)*delYY + px = numpy.array([[dX+dXX,dY+dYX]]) + + # Transform to sky coordinates + + wc = self.refwcs.pixel2world_fits(px) + + a = wc[0,0] + b = wc[0,1] + px = numpy.array([[dX+dXY,dY+dYY]]) + + wc = self.refwcs.pixel2world_fits(px) + c = wc[0,0] + d = wc[0,1] + + # Calculate the new CDs and convert to degrees + cd11 = diff_angles(a,newcrval[0])*cos(newcrval[1]*pi/180.0) + cd12 = diff_angles(c,newcrval[0])*cos(newcrval[1]*pi/180.0) + cd21 = diff_angles(b,newcrval[1]) + cd22 = diff_angles(d,newcrval[1]) + + cd = numpy.array([[cd11, cd12], [cd21, cd22]]) + self.wcsobj.cd = cd + + def uprefwcs(self, loff, offsh): + """ + Updates the reference chip + """ + ltvoffx, ltvoffy = loff[0], loff[1] + offshift = offsh + backup_refwcs = self.refwcs.restore() + dec = backup_refwcs['CRVAL2'] + + # Get an approximate reference position on the sky + rref = numpy.array([[self.refwcs.idcmodel.refpix['XREF']+ltvoffx, + self.refwcs.idcmodel.refpix['YREF']+ltvoffy]]) + + + crval = self.refwcs.pixel2world_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(self.wcsobj.pav3,dec,self.refwcs.idcmodel.refpix['V2REF'], + self.refwcs.idcmodel.refpix['V3REF']) + # Add the chip rotation angle + if self.refwcs.idcmodel.refpix['THETA']: + pv += self.refwcs.idcmodel.refpix['THETA'] + + + # Set values for the rest of the reference WCS + self.refwcs.crval = crval[0] + self.refwcs.crpix = numpy.array([0.0,0.0])+offsh + parity = self.refwcs.parity + R_scale = self.refwcs.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]]) + self.refwcs.cd = rcd + self.refwcs.set() + + + + def updatehdr(self): + """ + Keywords to update: + CD1_1, CD1_2, CD2_1, CD2_2, CRPIX1, CRPIX2, CRVAL1, CRVAL2 + """ + + self.wcsobj.hdr.update('CD1_1', self.wcsobj.cd[0,0]) + self.wcsobj.hdr.update('CD1_2', self.wcsobj.cd[0,1]) + self.wcsobj.hdr.update('CD2_1', self.wcsobj.cd[1,0]) + self.wcsobj.hdr.update('CD2_2', self.wcsobj.cd[1,1]) + self.wcsobj.hdr.update('CRVAL1', self.wcsobj.crval[0]) + self.wcsobj.hdr.update('CRVAL2', self.wcsobj.crval[1]) + self.wcsobj.hdr.update('CRPIX1', self.wcsobj.crpix[0]) + self.wcsobj.hdr.update('CRPIX2', self.wcsobj.crpix[1]) + self.wcsobj.hdr.update('ORIENTAT', self.wcsobj.orientat) + + def getOffsets(self, backup_wcs): + ltv1 = self.wcsobj.ltv1 + ltv2 = self.wcsobj.ltv2 + + offsetx = backup_wcs['CRPIX1'] - ltv1 - self.wcsobj.idcmodel.refpix['XREF'] + offsety = backup_wcs['CRPIX2'] - ltv2 - self.wcsobj.idcmodel.refpix['YREF'] + + shiftx = self.wcsobj.idcmodel.refpix['XREF'] + ltv1 + shifty = self.wcsobj.idcmodel.refpix['YREF'] + ltv2 + if ltv1 != 0. or ltv2 != 0.: + ltvoffx = ltv1 + offsetx + ltvoffy = ltv2 + offsety + offshiftx = offsetx + shiftx + offshifty = offsety + shifty + else: + ltvoffx = 0. + ltvoffy = 0. + offshiftx = 0. + offshifty = 0. + + ltvoff = numpy.array([ltvoffx, ltvoffy]) + offshift = numpy.array([offshiftx, offshifty]) + return ltvoff, offshift + +def diff_angles(a,b): + """ + Perform angle subtraction a-b taking into account + small-angle differences across 360degree line. + """ + + diff = a - b + + if diff > 180.0: + diff -= 360.0 + + if diff < -180.0: + diff += 360.0 + + return diff + +def troll(roll, dec, v2, v3): + """ Computes the roll angle at the target position based on: + the roll angle at the V1 axis(roll), + the dec of the target(dec), and + the V2/V3 position of the aperture (v2,v3) in arcseconds. + + Based on the algorithm provided by Colin Cox that is used in + Generic Conversion at STScI. + """ + # Convert all angles to radians + _roll = DEGTORAD(roll) + _dec = DEGTORAD(dec) + _v2 = DEGTORAD(v2 / 3600.) + _v3 = DEGTORAD(v3 / 3600.) + + # compute components + sin_rho = sqrt((pow(sin(_v2),2)+pow(sin(_v3),2)) - (pow(sin(_v2),2)*pow(sin(_v3),2))) + rho = asin(sin_rho) + beta = asin(sin(_v3)/sin_rho) + if _v2 < 0: beta = pi - beta + gamma = asin(sin(_v2)/sin_rho) + if _v3 < 0: gamma = pi - gamma + A = pi/2. + _roll - beta + B = atan2( sin(A)*cos(_dec), (sin(_dec)*sin_rho - cos(_dec)*cos(rho)*cos(A))) + + # compute final value + troll = RADTODEG(pi - (gamma+B)) + + return troll + diff --git a/wcsutil/__init__.py b/wcsutil/__init__.py new file mode 100644 index 0000000..46ad42e --- /dev/null +++ b/wcsutil/__init__.py @@ -0,0 +1,166 @@ +#from .. pywcs.sip import SIP +from pywcs.sip import SIP +import pyfits +import instruments +#from .. distortion import models +from hstwcs.distortion import models +import numpy as N + +#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 + +__docformat__ = 'restructuredtext' + +class HSTWCS(SIP): + """ + Purpose + ======= + Create a wcs object based on the instrument. + It has all basic WCS kw as attribbutes (set by pywcs). + It also uses the primary and extension header to define + instrument specific attributes needed by the correction classes. + """ + + def __init__(self, hdr0, hdr): + """ + :Parameters: + `hdr0`: Pyfits Header + primary header + `hdr`: Pyfits Header + extension header + """ + self.hdr = hdr + self.hdr0 = hdr0 + SIP.__init__(self, self.hdr) + self.setHDR0kw(hdr0) + self.detector = self.setDetector() + self.inst_kw = ins_spec_kw + self.setInstrSpecKw() + self.pscale = self.setPscale() + self.orientat = self.setOrient() + + + + def setHDR0kw(self, primhdr): + #sets attributes from kw defined in the primary header + 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) + + + def setDetector(self): + #sets detector attribute for instuments which have more than one detector + if self.instrument in ['ACS', 'WFC3']: + return self.hdr0.get('DETECTOR', None) + + + def setInstrSpecKw(self): + #Based on the instrument kw creates an instalnce of an instrument WCS class + #and sets attributes from instrument specific kw + if self.instrument in inst_mappings.keys(): + inst_kl = inst_mappings[self.instrument] + inst_kl = instruments.__dict__[inst_kl] + insobj = inst_kl(self.hdr0, self.hdr) + for key in self.inst_kw: + self.__setattr__(key, insobj.__getattribute__(key)) + else: + raise KeyError, "Unsupported instrument - %s" %self.instrument + + def setPscale(self): + #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] + return N.sqrt(N.power(cd11,2)+N.power(cd21,2)) * 3600. + + def setOrient(self): + #same as setPscale + cd12 = self.cd[0][1] + cd22 = self.cd[1][1] + return RADTODEG(N.arctan2(cd12,cd22)) + + def verifyPa_V3(self): + """make sure PA_V3 is populated""" + + def verifyKw(self): + """verify that all required kw have meaningful values""" + + def readModel(self): + """ + Purpose + ======= + Read distortion model from idc table. + Save some of the information as kw needed later by pydrizzle + + list of kw - to be revised later + """ + self.idcmodel = models.IDCModel(self.idctab, + chip=self.chip, direction='forward', date=self.date_obs, + filter1=self.filter1, filter2=self.filter2, + offtab=self.offtab, binned=self.binned) + + + self.updatehdr() + + + def updatehdr(self, newkeywords=None): + #kw2add : OCX10, OCX11, OCY10, OCY11 + # record the model in the header for use by pydrizzle + self.hdr.update('OCX10', self.idcmodel.cx[1,0]) + self.hdr.update('OCX11', self.idcmodel.cx[1,1]) + self.hdr.update('OCY10', self.idcmodel.cy[1,0]) + self.hdr.update('OCY11', self.idcmodel.cy[1,1]) + self.hdr.update('IDCSCALE', self.idcmodel.refpix['PSCALE']) + self.hdr.update('IDCTHETA', self.idcmodel.refpix['THETA']) + self.hdr.update('IDCXREF', self.idcmodel.refpix['XREF']) + self.hdr.update('IDCYREF', self.idcmodel.refpix['YREF']) + self.hdr.update('IDCV2REF', self.idcmodel.refpix['V2REF']) + self.hdr.update('IDCV3REF', self.idcmodel.refpix['V3REF']) + self.hdr.update('IDCXSIZE', self.idcmodel.refpix['XSIZE']) + self.hdr.update('IDCYSIZE', self.idcmodel.refpix['YSIZE']) + self.hdr.update('IDCXDELT', self.idcmodel.refpix['XDELTA']) + self.hdr.update('IDCYDELT', self.idcmodel.refpix['YDELTA']) + self.hdr.update('CENTERED', self.idcmodel.refpix['centered']) + + + def restore(self): + """ + restores basic wcs keywords from archive + this should be modified to always return a populated + dictionary, although the values may be None. + """ + + backup = {} + for k in basic_wcs: + try: + nkw = ('O'+k)[:7] + #backup[k] = self.hdr.__getitem__('O'+k) + backup[k] = self.hdr[nkw] + #self.__setattr__(k, self.hdr.__getitem__('O'+k)) + self.__setattr__(k, self.hdr[nkw]) + except KeyError: + print 'Keyword %s not found' % k + + return backup + + def archive_kw(self): + """ + archive original WCS kw before recalculating them. + """ + backup_kw = self.restore() + if backup_kw != {}: + return + else: + # keep the archived kw 8 characters long + cmt = 'archived value' + for kw in basic_wcs: + nkw = 'O'+kw + self.hdr.update(nkw[:7], self.hdr[kw], comment=cmt) + + diff --git a/wcsutil/instruments.py b/wcsutil/instruments.py new file mode 100644 index 0000000..2790d8a --- /dev/null +++ b/wcsutil/instruments.py @@ -0,0 +1,117 @@ +import pyfits +import numpy as N +#from .. mappings import ins_spec_kw +from hstwcs.mappings import ins_spec_kw, prim_hdr_kw + +class InstrWCS(object): + """ + A base class for instrument specific keyword definition. + It prvides a default implementation (modeled by ACS) for + all set_kw methods. + """ + def __init__(self, hdr0, hdr): + self.exthdr = hdr + self.primhdr = hdr0 + + def set_ins_spec_kw(self): + """ + This method MUST call all set_kw methods. + There should be a set_kw method for all kw listed in + mappings.ins_spec_kw + """ + self.set_filter1() + self.set_filter2() + self.set_vafactor() + self.set_naxis1() + self.set_naxis2() + self.set_ltv1() + self.set_ltv2() + self.set_binned() + self.set_chip() + self.set_parity() + + def set_filter1(self): + self.filter1 = self.primhdr.get('FILTER1', None) + + def set_filter2(self): + self.filter2 = self.primhdr.get('FILTER2', None) + + def set_vafactor(self): + self.vafactor = self.exthdr.get('vafactor', 1) + + def set_naxis1(self): + self.naxis1 = self.exthdr.get('naxis1', None) + + def set_naxis2(self): + self.naxis2 = self.exthdr.get('naxis2', None) + + def set_ltv1(self): + self.ltv1 = self.exthdr.get('ltv1', 0.0) + + def set_ltv2(self): + self.ltv2 = self.exthdr.get('ltv2', 0.0) + + def set_binned(self): + self.binned = self.exthdr.get('BINAXIS1', 1) + + def set_chip(self): + self.chip = self.exthdr.get('CCDCHIP', 1) + + def set_parity(self): + self.parity = [[1.0,0.0],[0.0,-1.0]] + +class ACSWCS(InstrWCS): + """ + get instrument specific kw + """ + + def __init__(self, hdr0, hdr): + self.primhdr = hdr0 + self.exthdr = hdr + InstrWCS.__init__(self,hdr0, hdr) + self.set_ins_spec_kw() + + + def set_parity(self): + parity = {'WFC':[[1.0,0.0],[0.0,-1.0]], + 'HRC':[[-1.0,0.0],[0.0,1.0]], + 'SBC':[[-1.0,0.0],[0.0,1.0]]} + detector = self.primhdr.get('detector', None) + print 'detector', detector + self.parity = parity[detector] + + +class WFPC2WCS(InstrWCS): + + + def __init__(self, hdr0, hdr): + self.primhdr = hdr0 + self.exthdr = hdr + InstrWCS.__init__(self,hdr0, hdr) + self.set_ins_spec_kw() + + def set_filter1(self): + self.filter1 = self.primhdr.get('FILTNAM1', None) + if self.filter1 == " " or self.filter1 == None: + self.filter1 = 'CLEAR1' + + def set_filter2(self): + self.filter2 = self.primhdr.get('FILTNAM2', None) + if self.filter2 == " " or self.filter2 == None: + self.filter2 = 'CLEAR1' + + + def set_binned(self): + mode = self.primhdr.get('MODE', 1) + if mode == 'FULL': + self.binned = 1 + elif mode == 'AREA': + self.binned = 2 + + def set_chip(self): + self.chip = self.exthdr.get('DETECTOR', 1) + + def set_parity(self): + self.parity = [[-1.0,0.],[0.,1.0]] + + |