diff options
author | Nadia Dencheva <nadia.astropy@gmail.com> | 2016-08-11 14:42:28 -0400 |
---|---|---|
committer | GitHub <noreply@github.com> | 2016-08-11 14:42:28 -0400 |
commit | 892edca9d767cf2dfea48a60c9cc0eb97ef2e4e0 (patch) | |
tree | d81f29b7f54ff7132ac1073a2b9c2186da34aef9 /stwcs | |
parent | 6ee1b08a2bc2fea4e61fb05d6c3d9250c15a1a75 (diff) | |
parent | 9932a03e2abe7c871195af0dc21f43415c12027d (diff) | |
download | stwcs_hcf-892edca9d767cf2dfea48a60c9cc0eb97ef2e4e0.tar.gz |
Merge pull request #12 from nden/pep8-pylint
Pep8 pylint
Diffstat (limited to 'stwcs')
35 files changed, 1474 insertions, 1427 deletions
diff --git a/stwcs/distortion/coeff_converter.py b/stwcs/distortion/coeff_converter.py index 415b512..99cf74d 100644 --- a/stwcs/distortion/coeff_converter.py +++ b/stwcs/distortion/coeff_converter.py @@ -1,9 +1,10 @@ -from __future__ import division, print_function # confidence high +from __future__ import (absolute_import, unicode_literals, division, print_function) import numpy as np from astropy.io import fits from astropy import wcs as pywcs + def sip2idc(wcs): """ Converts SIP style coefficients to IDCTAB coefficients. @@ -49,28 +50,28 @@ def sip2idc(wcs): 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 None, None 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]]]) + 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[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 + 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. @@ -84,15 +85,15 @@ def _read_sip_kw(header): "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 = 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 = 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 "B_ORDER" in header: raise ValueError( @@ -102,7 +103,7 @@ def _read_sip_kw(header): a = None b = None - return a , b + return a, b """ diff --git a/stwcs/distortion/models.py b/stwcs/distortion/models.py index 231a9f1..136c6cf 100644 --- a/stwcs/distortion/models.py +++ b/stwcs/distortion/models.py @@ -1,14 +1,9 @@ -from __future__ import absolute_import, division, print_function # confidence high - +from __future__ import absolute_import, division, print_function import numpy as np - -# Import PyDrizzle utility modules from . import mutil from .mutil import combin -yes = True -no = False ################# # @@ -72,8 +67,8 @@ class GeometryModel: to the reference position of the chip. """ - _cxs = np.zeros(shape=self.cx.shape,dtype=self.cx.dtype) - _cys = np.zeros(shape=self.cy.shape,dtype=self.cy.dtype) + _cxs = np.zeros(shape=self.cx.shape, dtype=self.cx.dtype) + _cys = np.zeros(shape=self.cy.shape, dtype=self.cy.dtype) _k = self.norder + 1 # loop over each input coefficient for m in range(_k): @@ -83,15 +78,17 @@ class GeometryModel: _ilist = list(range(m, _k)) # sum from m to k for i in _ilist: - _jlist = list(range(n, i - (m-n)+1)) + _jlist = list(range(n, i - (m - n) + 1)) # sum from n to i-(m-n) for j in _jlist: - _cxs[m,n] += self.cx[i,j]*combin(j,n)*combin((i-j),(m-n))*pow(xs,(j-n))*pow(ys,((i-j)-(m-n))) - _cys[m,n] += self.cy[i,j]*combin(j,n)*combin((i-j),(m-n))*pow(xs,(j-n))*pow(ys,((i-j)-(m-n))) + _cxs[m, n] += self.cx[i, j] * combin(j, n) * combin((i - j), (m - n)) * \ + pow(xs, (j - n)) * pow(ys, ((i - j) - (m - n))) + _cys[m, n] += self.cy[i, j] * combin(j, n) * combin((i - j), (m - n)) * \ + pow(xs, (j - n)) * pow(ys, ((i - j) - (m - n))) self.cx = _cxs.copy() self.cy = _cys.copy() - def convert(self, tmpname, xref=None,yref=None,delta=yes): + def convert(self, tmpname, xref=None, yref=None, delta=True): """ Open up an ASCII file, output coefficients in drizzle format after converting them as necessary. @@ -104,89 +101,62 @@ class GeometryModel: 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'] + 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['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 + self.x0 = x0 + self.y0 = y0 # 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('# 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: + elif self.norder == 4: lines.append('quartic\n') - elif self.norder==5: + 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]) + 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]) + 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]) + 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]) + 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]) + 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]) + 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]) + 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]) + 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 = open(tmpname, 'w') output.writelines(lines) output.close() - - def apply(self, pixpos,scale=1.0,order=None): + 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. @@ -195,13 +165,13 @@ class GeometryModel: 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 self.cx is None: + return pixpos[:, 0], pixpos[:, 1] if order is None: order = self.norder @@ -209,29 +179,29 @@ class GeometryModel: # Apply in the same way that 'drizzle' would... _cx = self.cx / (self.pscale * scale) _cy = self.cy / (self.pscale * scale) - _convert = no + _convert = False _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. + _cx[0, 0] = 0. + _cy[0, 0] = 0. if isinstance(_p, list) or isinstance(_p, tuple): - _p = np.array(_p,dtype=np.float64) - _convert = yes + _p = np.array(_p, dtype=np.float64) + _convert = True - dxy = _p - (self.refpix['XREF'],self.refpix['YREF']) + 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] + 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: @@ -242,11 +212,11 @@ class GeometryModel: xc = xc[0] yc = yc[0] - return xc,yc + return xc, yc - def setPScaleCoeffs(self,pscale): - self.cx[1,1] = pscale - self.cy[1,0] = pscale + def setPScaleCoeffs(self, pscale): + self.cx[1, 1] = pscale + self.cy[1, 0] = pscale self.refpix['PSCALE'] = pscale self.pscale = pscale @@ -259,7 +229,7 @@ class IDCModel(GeometryModel): 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): + filter1='CLEAR1', filter2='CLEAR2', offtab=None, binned=1): GeometryModel.__init__(self) # # Norder must be derived from the coeffs file itself, @@ -269,9 +239,9 @@ class IDCModel(GeometryModel): # 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) + self.cx, self.cy, self.refpix, self.norder = mutil.readIDCtab(idcfile, + chip=chip, direction=direction, filter1=filter1, filter2=filter2, + date=date, offtab=offtab) if 'empty_model' in self.refpix and self.refpix['empty_model']: pass @@ -285,61 +255,54 @@ class IDCModel(GeometryModel): 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): + def __init__(self, header, rootname): GeometryModel.__init__(self) - if 'rootname' in header: 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() + # 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.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): + 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) + self.cx, self.cy, self.refpix, self.norder = mutil.readCubicTable(idcfile) - # scale is the ratio wcs.pscale/model.pscale. + # 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: + + if scale is not 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): """ @@ -349,14 +312,8 @@ class TraugerModel(GeometryModel): """ NORDER = 3 - def __init__(self, idcfile,lam): + def __init__(self, idcfile, lam): GeometryModel.__init__(self) self.name = idcfile - self.cx,self.cy,self.refpix,self.norder = mutil.readTraugerTable(idcfile,lam) + 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/stwcs/distortion/mutil.py b/stwcs/distortion/mutil.py index ed6a1ea..0e193c8 100644 --- a/stwcs/distortion/mutil.py +++ b/stwcs/distortion/mutil.py @@ -1,25 +1,22 @@ -from __future__ import division, print_function # confidence high +from __future__ import division, print_function from stsci.tools import fileutil import numpy as np 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. +# 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): + +def readIDCtab(tabname, chip=1, date=None, direction='forward', + filter1=None, filter2=None, offtab=None): """ Read IDCTAB, and optional OFFTAB if sepcified, and generate @@ -30,17 +27,17 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward', """ - # Return a default geometry model if no IDCTAB filename + # Return a default geometry model if no IDCTAB filename # is given. This model will not distort the data in any way. - if tabname == None: + if tabname is 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 or filter1.strip() == '': + if filter1 is None or filter1.find('CLEAR') == 0 or filter1.strip() == '': filter1 = 'CLEAR' - if filter2 == None or filter2.find('CLEAR') == 0 or filter2.strip() == '': + if filter2 is None or filter2.find('CLEAR') == 0 or filter2.strip() == '': filter2 = 'CLEAR' # Insure that tabname is full filename with fully expanded @@ -51,7 +48,7 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward', try: ftab = fileutil.openImage(tabname) except: - err_str = "------------------------------------------------------------------------ \n" + 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" @@ -63,7 +60,7 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward', err_str += "------------------------------------------------------------------------ \n" raise IOError(err_str) - #First thing we need, is to read in the coefficients from the IDC + # First thing we need, is to read in the coefficients from the IDC # table and populate the Fx and Fy matrices. if 'DETECTOR' in ftab['PRIMARY'].header: @@ -96,10 +93,10 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward', 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) + 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. + # 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 @@ -111,14 +108,14 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward', try: # Match FILTER combo to appropriate row, - #if there is a filter column in the IDCTAB... + # 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] + if filt1.find('CLEAR') > -1: filt1 = filt1[: 5] filt2 = ftab[1].data.field('FILTER2')[i] - if filt2.find('CLEAR') > -1: filt2 = filt2[:5] + if filt2.find('CLEAR') > -1: filt2 = filt2[: 5] else: if 'OPT_ELEM' in colnames: filt1 = ftab[1].data.field('OPT_ELEM') @@ -170,16 +167,17 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward', joinstr = '' else: f2str = filter2.strip() - filtstr = (joinstr.join([f1str,f2str])).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: '+filtstr+'\n' + err_str += ' CHIP: ' + str(detchip) + '\n' + err_str += ' FILTERS: ' + filtstr + '\n' ftab.close() del ftab raise LookupError(err_str) else: - print('- IDCTAB: Distortion model from row',str(row+1),'for chip',detchip,':',filtstr) + 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) @@ -190,12 +188,12 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward', else: # Read V2REF/V3REF from offset table (OFFTAB) if offtab: - v2ref,v3ref,theta = readOfftab(offtab, date, chip=detchip) + v2ref, v3ref, theta = readOfftab(offtab, date, chip=detchip) else: v2ref = 0.0 v3ref = 0.0 - if theta == None: + if theta is None: if 'THETA' in colnames: theta = ftab[1].data.field('THETA')[row] else: @@ -206,14 +204,14 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward', 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['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['DEFAULT_SCALE'] = True + refpix['centered'] = False refpix['skew_coeffs'] = skew_coeffs # Now that we know which row to look at, read coefficients into the # numeric arrays we have set up... @@ -226,13 +224,13 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward', cxstr = 'A' cystr = 'B' - for i in range(norder+1): + for i in range(norder + 1): if i > 0: - for j in range(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] + for j in range(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 @@ -240,21 +238,21 @@ def readIDCtab (tabname, chip=1, date=None, direction='forward', # 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']: + 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) -# -# + return fx, fy, refpix, order + + def read_tdd_coeffs(phdr, chip=1): - ''' Read in the TDD related keywords from the PRIMARY header of the IDCTAB - ''' + """ + Time-dependent skew correction coefficients (only ACS/WFC). + + Read in the TDD related keywords from the PRIMARY header of the IDCTAB + """ # Insure we have an integer form of chip ic = int(chip) @@ -269,7 +267,7 @@ def read_tdd_coeffs(phdr, chip=1): skew_coeffs['TDD_CX_ALPHA'] = None # Skew-based TDD coefficients - skew_terms = ['TDD_CTB','TDD_CTA','TDD_CYA','TDD_CYB','TDD_CXA','TDD_CXB'] + skew_terms = ['TDD_CTB', 'TDD_CTA', 'TDD_CYA', 'TDD_CYB', 'TDD_CXA', 'TDD_CXB'] for s in skew_terms: skew_coeffs[s] = None @@ -280,7 +278,7 @@ def read_tdd_coeffs(phdr, chip=1): print("Using 2015-calibrated VAFACTOR-corrected TDD correction...") skew_coeffs['TDD_DATE'] = phdr['TDD_DATE'] for s in skew_terms: - skew_coeffs[s] = phdr.get('{0}{1}'.format(s,ic),None) + skew_coeffs[s] = phdr.get('{0}{1}'.format(s, ic), None) elif "TDD_CYB1" in phdr: # We have 2014-calibrated TDD correction to apply, not J.A.-derived values @@ -288,17 +286,17 @@ def read_tdd_coeffs(phdr, chip=1): skew_coeffs['TDD_DATE'] = phdr['TDD_DATE'] # Read coefficients for TDD Y coefficient cyb_kw = 'TDD_CYB{0}'.format(int(chip)) - skew_coeffs['TDD_CY_BETA'] = phdr.get(cyb_kw,None) + skew_coeffs['TDD_CY_BETA'] = phdr.get(cyb_kw, None) cya_kw = 'TDD_CYA{0}'.format(int(chip)) - tdd_cya = phdr.get(cya_kw,None) + tdd_cya = phdr.get(cya_kw, None) if tdd_cya == 0 or tdd_cya == 'N/A': tdd_cya = None skew_coeffs['TDD_CY_ALPHA'] = tdd_cya # Read coefficients for TDD X coefficient cxb_kw = 'TDD_CXB{0}'.format(int(chip)) - skew_coeffs['TDD_CX_BETA'] = phdr.get(cxb_kw,None) + skew_coeffs['TDD_CX_BETA'] = phdr.get(cxb_kw, None) cxa_kw = 'TDD_CXA{0}'.format(int(chip)) - tdd_cxa = phdr.get(cxa_kw,None) + tdd_cxa = phdr.get(cxa_kw, None) if tdd_cxa == 0 or tdd_cxa == 'N/A': tdd_cxa = None skew_coeffs['TDD_CX_ALPHA'] = tdd_cxa @@ -309,12 +307,12 @@ def read_tdd_coeffs(phdr, chip=1): 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 = 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(): + 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) @@ -325,15 +323,15 @@ def read_tdd_coeffs(phdr, chip=1): 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. +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 is None: + return 0., 0. # Provide a default value for chip if chip: @@ -347,7 +345,7 @@ def readOfftab(offtab, date, chip=None): except: raise IOError("Offset table '%s' not valid as specified!" % offtab) - #Determine row from which to get the coefficients. + # 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 @@ -374,7 +372,7 @@ def readOfftab(offtab, date, chip=None): obsdate = convertDate(ftab[1].data.field('OBSDATE')[i]) # If the row is appropriate for the chip... - # Interpolate between dates + # Interpolate between dates if int(detchip) == int(chip) or int(detchip) == -999: if num_date <= obsdate: date_end = obsdate @@ -384,7 +382,7 @@ def readOfftab(offtab, date, chip=None): row_end = i continue - if row_end == None and (num_date > obsdate): + if row_end is None and (num_date > obsdate): date_end = obsdate v2end = ftab[1].data.field('V2REF')[i] v3end = ftab[1].data.field('V3REF')[i] @@ -403,16 +401,17 @@ def readOfftab(offtab, date, chip=None): ftab.close() del ftab - if row_start == None and row_end == None: - print('Row corresponding to DETCHIP of ',detchip,' was not found!') + if row_start is None and row_end is 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)) + elif row_start is 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)) + 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: + if row_start is None or row_end is row_start: # We are processing an observation taken after the last calibration date_start = date_end v2start = v2end @@ -426,31 +425,32 @@ def readOfftab(offtab, date, chip=None): v3ref = _fraction * (v3end - v3start) + v3start theta = _fraction * (theta_end - theta_start) + theta_start - return v2ref,v3ref,theta + return v2ref, v3ref, theta -def readWCSCoeffs(header): - - #Read distortion coeffs from WCS header keywords and - #populate distortion coeffs arrays. +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) + 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) + 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) + _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. @@ -467,40 +467,40 @@ def readWCSCoeffs(header): refpix['THETA'] = np.rad2deg(_theta) refpix['XDELTA'] = 0.0 refpix['YDELTA'] = 0.0 - refpix['DEFAULT_SCALE'] = yes - refpix['centered'] = yes - + refpix['DEFAULT_SCALE'] = True + refpix['centered'] = True # Set up template for coeffs keyword names cxstr = 'A_' cystr = 'B_' # Read coeffs into their own matrix - for i in range(_xorder+1): - for j in range(i+1): - xcname = cxstr+str(j)+'_'+str(i-j) + for i in range(_xorder + 1): + for j in range(i + 1): + xcname = cxstr + str(j) + '_' + str(i - j) if xcname in header: - fx[i,j] = header[xcname] + fx[i, j] = header[xcname] # Extract Y coeffs separately as a different order may # have been used to fit it. - for i in range(_yorder+1): - for j in range(i+1): - ycname = cystr+str(j)+'_'+str(i-j) + for i in range(_yorder + 1): + for j in range(i + 1): + ycname = cystr + str(j) + '_' + str(i - j) if ycname in header: - fy[i,j] = header[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 + 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: +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 is None: return fileutil.defaultModel() # Trauger coefficients only result in a cubic file... @@ -510,10 +510,10 @@ def readTraugerTable(idcfile,wavelength): b_coeffs = [0] * numco indx = _MgF2(wavelength) - ifile = open(idcfile,'r') + ifile = open(idcfile, 'r') # Search for the first line of the coefficients _line = fileutil.rAsciiLine(ifile) - while _line[:7].lower() != 'trauger': + while _line[: 7].lower() != 'trauger': _line = fileutil.rAsciiLine(ifile) # Read in each row of coefficients,split them into their values, # and convert them into cubic coefficients based on @@ -525,9 +525,11 @@ def readTraugerTable(idcfile,wavelength): if _line == '': continue _lc = _line.split() if j < 10: - a_coeffs[j] = float(_lc[0])+float(_lc[1])*(indx-1.5)+float(_lc[2])*(indx-1.5)**2 + 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 + 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() @@ -536,17 +538,17 @@ def readTraugerTable(idcfile,wavelength): # 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) + 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) + 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 = {} @@ -557,10 +559,10 @@ def readTraugerTable(idcfile,wavelength): refpix['XDELTA'] = 0. refpix['YDELTA'] = 0. refpix['PSCALE'] = None - refpix['DEFAULT_SCALE'] = no - refpix['centered'] = yes + refpix['DEFAULT_SCALE'] = False + refpix['centered'] = True - return fx,fy,refpix,order + return fx, fy, refpix, order def readCubicTable(idcfile): @@ -572,17 +574,17 @@ def readCubicTable(idcfile): # Return a default geometry model if no coefficients filename # is given. This model will not distort the data in any way. - if idcfile == None: + if idcfile is None: return fileutil.defaultModel() - ifile = open(idcfile,'r') + 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 + _found = False + while not _found: + if _line[:7] in ['cubic', 'quartic', 'quintic'] or _line[: 4] == 'poly': + found = True break _line = fileutil.rAsciiLine(ifile) @@ -613,17 +615,17 @@ def readCubicTable(idcfile): # 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) + 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) + 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 = {} @@ -634,21 +636,23 @@ def readCubicTable(idcfile): refpix['XDELTA'] = 0. refpix['YDELTA'] = 0. refpix['PSCALE'] = None - refpix['DEFAULT_SCALE'] = no - refpix['centered'] = yes + refpix['DEFAULT_SCALE'] = False + refpix['centered'] = True + + return fx, fy, refpix, order - 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) + m = m * (i + 1) return m -def combin(j,n): + +def combin(j, n): """ Return the combinatorial factor for j in n.""" - return (factorial(j) / (factorial(n) * factorial( (j-n) ) ) ) + return (factorial(j) / (factorial(n) * factorial((j - n)))) def defaultModel(): @@ -657,15 +661,15 @@ def defaultModel(): """ 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 = 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. + fx[1, 1] = 1. + fy[1, 0] = 1. # Used in Pattern.computeOffsets() refpix = {} - refpix['empty_model'] = yes + refpix['empty_model'] = True refpix['XREF'] = None refpix['YREF'] = None refpix['V2REF'] = 0. @@ -675,17 +679,20 @@ def defaultModel(): refpix['XDELTA'] = 0. refpix['YDELTA'] = 0. refpix['PSCALE'] = None - refpix['DEFAULT_SCALE'] = no + refpix['DEFAULT_SCALE'] = False refpix['THETA'] = 0. - refpix['centered'] = yes - return fx,fy,refpix,order + refpix['centered'] = True + 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)) + """ + Function to compute the index of refraction for MgF2 at + the specified wavelength for use with Trauger coefficients + """ + _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): diff --git a/stwcs/distortion/utils.py b/stwcs/distortion/utils.py index 4228f62..1baad3f 100644 --- a/stwcs/distortion/utils.py +++ b/stwcs/distortion/utils.py @@ -1,4 +1,4 @@ -from __future__ import division, print_function # confidence high +from __future__ import absolute_import, division, print_function import os @@ -6,11 +6,11 @@ import numpy as np from numpy import linalg from astropy import wcs as pywcs -from stwcs import wcsutil -from stwcs import updatewcs +from .. import updatewcs from numpy import sqrt, arctan2 from stsci.tools import fileutil + def output_wcs(list_of_wcsobj, ref_wcs=None, owcs=None, undistort=True): """ Create an output WCS. @@ -33,63 +33,65 @@ def output_wcs(list_of_wcsobj, ref_wcs=None, owcs=None, undistort=True): # 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) + crval1, crval2 = computeFootprintCenter(fra_dec) - crval = np.array([crval1,crval2], dtype=np.float64) # this value is now zero-based + crval = np.array([crval1, crval2], dtype=np.float64) # this value is now zero-based if owcs is None: if ref_wcs is None: ref_wcs = list_of_wcsobj[0].deepcopy() if undistort: - #outwcs = undistortWCS(ref_wcs) + # outwcs = undistortWCS(ref_wcs) outwcs = make_orthogonal_cd(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 + 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 + 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._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()]) + 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() - outwcs.wcs.name = wcsname # keep track of label for this solution + outwcs.wcs.name = wcsname # keep track of label for this solution 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 = np.deg2rad(edges[:,0]) - dec = np.deg2rad(edges[:,1]) + alpha = np.deg2rad(edges[:, 0]) + dec = np.deg2rad(edges[:, 1]) - xmean = np.mean(np.cos(dec)*np.cos(alpha)) - ymean = np.mean(np.cos(dec)*np.sin(alpha)) + 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 = np.rad2deg(np.arctan2(ymean,xmean))%360.0 - crval2 = np.rad2deg(np.arctan2(zmean,np.sqrt(xmean*xmean+ymean*ymean))) + crval1 = np.rad2deg(np.arctan2(ymean, xmean)) % 360.0 + crval2 = np.rad2deg(np.arctan2(zmean, np.sqrt(xmean * xmean + ymean * ymean))) + + return crval1, crval2 - return crval1,crval2 def make_orthogonal_cd(wcs): """ Create a perfect (square, orthogonal, undistorted) CD matrix from the @@ -98,21 +100,20 @@ def make_orthogonal_cd(wcs): # get determinant of the CD matrix: cd = wcs.celestial.pixel_scale_matrix - if hasattr(wcs, 'idcv2ref') and wcs.idcv2ref is not None: # 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 = updatewcs.makewcs.troll(wcs.pav3,wcs.wcs.crval[1],wcs.idcv2ref,wcs.idcv3ref) + pv = updatewcs.makewcs.troll(wcs.pav3, wcs.wcs.crval[1], wcs.idcv2ref, wcs.idcv3ref) # Add the chip rotation angle if wcs.idctheta: pv += wcs.idctheta cs = np.cos(np.deg2rad(pv)) sn = np.sin(np.deg2rad(pv)) - pvmat = np.dot(np.array([[cs,sn],[-sn,cs]]),wcs.parity) - rot = np.arctan2(pvmat[0,1],pvmat[1,1]) - scale = wcs.idcscale/3600. + pvmat = np.dot(np.array([[cs, sn], [-sn, cs]]), wcs.parity) + rot = np.arctan2(pvmat[0, 1], pvmat[1, 1]) + scale = wcs.idcscale / 3600. det = linalg.det(wcs.parity) @@ -122,36 +123,37 @@ def make_orthogonal_cd(wcs): # find pixel scale: if hasattr(wcs, 'idcscale'): - scale = (wcs.idcscale) / 3600. # HST pixel scale provided + scale = (wcs.idcscale) / 3600. # HST pixel scale provided else: - scale = np.sqrt(np.abs(det)) # find as sqrt(pixel area) + scale = np.sqrt(np.abs(det)) # find as sqrt(pixel area) # find Y-axis orientation: if hasattr(wcs, 'orientat') and not ignoreHST: - rot = np.deg2rad(wcs.orientat) # use HST ORIENTAT + rot = np.deg2rad(wcs.orientat) # use HST ORIENTAT else: - rot = np.arctan2(wcs.wcs.cd[0,1], wcs.wcs.cd[1,1]) # angle of the Y-axis + rot = np.arctan2(wcs.wcs.cd[0, 1], wcs.wcs.cd[1, 1]) # angle of the Y-axis par = -1 if det < 0.0 else 1 # create a perfectly square, orthogonal WCS sn = np.sin(rot) cs = np.cos(rot) - orthogonal_cd = scale * np.array([[par*cs, sn], [-par*sn, cs]]) + orthogonal_cd = scale * np.array([[par * cs, sn], [-par * sn, cs]]) lin_wcsobj = pywcs.WCS() lin_wcsobj.wcs.cd = orthogonal_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.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 undistortWCS(wcsobj): + +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 @@ -175,25 +177,26 @@ def undistortWCS(wcsobj): 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) + 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 + # order = wcsobj.sip.a_order pscale = wcsobj.idcscale - #pixref = np.array([wcsobj.sip.SIPREF1, wcsobj.sip.SIPREF2]) + # 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] + 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) + 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): + if (_det == 0.0): print('Singular matrix in updateWCS, aborting ...') return @@ -202,15 +205,16 @@ def undistortWCS(wcsobj): 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.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): + +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. @@ -233,27 +237,28 @@ def apply_idc(pixpos, cx, cy, pixref, pscale= None, order=None): return pixpos # Apply in the same way that 'drizzle' would... - _cx = cx/pscale - _cy = cy/ pscale + _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. + _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)) + 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 - return c def foundIDCTAB(idctab): idctab_found = True diff --git a/stwcs/gui/apply_headerlet.py b/stwcs/gui/apply_headerlet.py index d517e9f..0cef612 100644 --- a/stwcs/gui/apply_headerlet.py +++ b/stwcs/gui/apply_headerlet.py @@ -1,31 +1,34 @@ +from __future__ import absolute_import, division, print_function import os from stsci.tools import teal, parseinput -import stwcs -from stwcs.wcsutil import headerlet +from .. import __version__ +from ..wcsutil import headerlet -__taskname__ = __name__.split('.')[-1] # needed for help string +__taskname__ = __name__.split('.')[-1] # needed for help string __package__ = headerlet.__name__ -__version__ = stwcs.__version__ -# -#### Interfaces used by TEAL -# +# __version__ = stwcs.__version__ + + +############### Interfaces used by TEAL ############### + def getHelpAsString(docstring=False): """ return useful help from a file in the script directory called __taskname__.help """ install_dir = os.path.dirname(__file__) - htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html') - helpfile = os.path.join(install_dir,__taskname__+'.help') + htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html') + helpfile = os.path.join(install_dir, __taskname__ + '.help') if docstring or (not docstring and not os.path.exists(htmlfile)): - helpString = __taskname__+' Version '+__version__+'\n\n' + helpString = __taskname__ + ' Version ' + __version__ + '\n\n' if os.path.exists(helpfile): - helpString += teal.getHelpFileAsString(__taskname__,__file__) + helpString += teal.getHelpFileAsString(__taskname__, __file__) else: - helpString = 'file://'+htmlfile + helpString = 'file://' + htmlfile return helpString + def run(configObj=None): # start by interpreting filename and hdrlet inputs @@ -37,18 +40,18 @@ def run(configObj=None): # Syntax: apply_headerlet_as_primary(filename, hdrlet, attach=True, # archive=True, force=False, verbose=False) headerlet.apply_headerlet_as_primary(filename, - hdrlet,attach=configObj['attach'], - archive=configObj['archive'],force=configObj['force'], - logging=configObj['logging']) + hdrlet, attach=configObj['attach'], + archive=configObj['archive'], force=configObj['force'], + logging=configObj['logging']) else: wcsname = configObj['wcsname'] - if wcsname in ['',' ','INDEF']: wcsname = None + if wcsname in ['', ' ', 'INDEF']: wcsname = None wcskey = configObj['wcskey'] if wcskey == '': wcskey = None # Call function with properly interpreted input parameters # apply_headerlet_as_alternate(filename, hdrlet, attach=True, # wcskey=None, wcsname=None, verbose=False) headerlet.apply_headerlet_as_alternate(filename, - hdrlet, attach=configObj['attach'], - wcsname=wcsname, wcskey=wcskey, - logging=configObj['logging']) + hdrlet, attach=configObj['attach'], + wcsname=wcsname, wcskey=wcskey, + logging=configObj['logging']) diff --git a/stwcs/gui/archive_headerlet.py b/stwcs/gui/archive_headerlet.py index 7ad3d4d..577a348 100644 --- a/stwcs/gui/archive_headerlet.py +++ b/stwcs/gui/archive_headerlet.py @@ -1,49 +1,49 @@ -from __future__ import print_function +from __future__ import absolute_import, division, print_function import os -from astropy.io import fits from stsci.tools import teal +from .. import __version__ +from ..wcsutil import headerlet -import stwcs -from stwcs.wcsutil import headerlet - -__taskname__ = __name__.split('.')[-1] # needed for help string +__taskname__ = __name__.split('.')[-1] # needed for help string __package__ = headerlet.__name__ -__version__ = stwcs.__version__ +# __version__ = stwcs.__version__ + # #### Interfaces used by TEAL # def getHelpAsString(docstring=False): """ - return useful help from a file in the script directory called __taskname__.help + Return useful help from a file in the script directory called __taskname__.help """ install_dir = os.path.dirname(__file__) - htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html') - helpfile = os.path.join(install_dir,__taskname__+'.help') + htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html') + helpfile = os.path.join(install_dir, __taskname__ + '.help') if docstring or (not docstring and not os.path.exists(htmlfile)): - helpString = __taskname__+' Version '+__version__+'\n\n' + helpString = __taskname__ + ' Version ' + __version__ + '\n\n' if os.path.exists(helpfile): - helpString += teal.getHelpFileAsString(__taskname__,__file__) + helpString += teal.getHelpFileAsString(__taskname__, __file__) else: helpString += headerlet.archive_as_headerlet.__doc__ else: - helpString = 'file://'+htmlfile + helpString = 'file://' + htmlfile return helpString + def run(configObj=None): - if configObj['hdrname'] in ['',' ','INDEF']: - print('='*60) + if configObj['hdrname'] in ['', ' ', 'INDEF']: + print('=' * 60) print('ERROR:') print(' No valid "hdrname" parameter value provided!') print(' Please restart this task and provide a value for this parameter.') - print('='*60) + print('=' * 60) return - str_kw = ['wcsname','destim','sipname','npolfile','d2imfile', - 'descrip','history','author'] + str_kw = ['wcsname', 'destim', 'sipname', 'npolfile', 'd2imfile', + 'descrip', 'history', 'author'] # create dictionary of remaining parameters, deleting extraneous ones # such as those above @@ -66,4 +66,4 @@ def run(configObj=None): # author=None, descrip=None, history=None, # hdrlet=None, clobber=False) headerlet.archive_as_headerlet(configObj['filename'], configObj['hdrname'], - **cdict) + **cdict) diff --git a/stwcs/gui/attach_headerlet.py b/stwcs/gui/attach_headerlet.py index 873c549..f3ab690 100644 --- a/stwcs/gui/attach_headerlet.py +++ b/stwcs/gui/attach_headerlet.py @@ -1,13 +1,13 @@ +from __future__ import absolute_import, division, print_function import os from stsci.tools import teal +from .. import __version +from ..wcsutil import headerlet -import stwcs -from stwcs.wcsutil import headerlet - -__taskname__ = __name__.split('.')[-1] # needed for help string +__taskname__ = __name__.split('.')[-1] # needed for help string __package__ = headerlet.__name__ -__version__ = stwcs.__version__ +#__version__ = stwcs.__version__ # #### Interfaces used by TEAL # @@ -16,21 +16,21 @@ def getHelpAsString(docstring=False): return useful help from a file in the script directory called __taskname__.help """ install_dir = os.path.dirname(__file__) - htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html') - helpfile = os.path.join(install_dir,__taskname__+'.help') + htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html') + helpfile = os.path.join(install_dir, __taskname__ + '.help') if docstring or (not docstring and not os.path.exists(htmlfile)): - helpString = __taskname__+' Version '+__version__+'\n\n' + helpString = __taskname__ + ' Version ' + __version__ + '\n\n' if os.path.exists(helpfile): - helpString += teal.getHelpFileAsString(__taskname__,__file__) + helpString += teal.getHelpFileAsString(__taskname__, __file__) else: - helpString += eval('.'.join([__package__,__taskname__,'__doc__'])) + helpString += eval('.'.join([__package__, __taskname__, '__doc__'])) else: - helpString = 'file://'+htmlfile + helpString = 'file://' + htmlfile return helpString + def run(configObj=None): - headerlet.attach_headerlet(configObj['filename'],configObj['hdrlet'], + headerlet.attach_headerlet(configObj['filename'], configObj['hdrlet'], configObj['logging']) - diff --git a/stwcs/gui/delete_headerlet.py b/stwcs/gui/delete_headerlet.py index b3df5a7..29937f7 100644 --- a/stwcs/gui/delete_headerlet.py +++ b/stwcs/gui/delete_headerlet.py @@ -1,15 +1,15 @@ -from __future__ import print_function -import os +from __future__ import absolute_import, division, print_function +import os from stsci.tools import teal from stsci.tools import parseinput -import stwcs -from stwcs.wcsutil import headerlet +from .. import __version__ +from ..wcsutil import headerlet -__taskname__ = __name__.split('.')[-1] # needed for help string +__taskname__ = __name__.split('.')[-1] # needed for help string __package__ = headerlet.__name__ -__version__ = stwcs.__version__ +# __version__ = stwcs.__version__ # #### Interfaces used by TEAL # @@ -18,35 +18,34 @@ def getHelpAsString(docstring=False): return useful help from a file in the script directory called __taskname__.help """ install_dir = os.path.dirname(__file__) - htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html') - helpfile = os.path.join(install_dir,__taskname__+'.help') + htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html') + helpfile = os.path.join(install_dir, __taskname__ + '.help') if docstring or (not docstring and not os.path.exists(htmlfile)): - helpString = __taskname__+' Version '+__version__+'\n\n' + helpString = __taskname__ + ' Version ' + __version__ + '\n\n' if os.path.exists(helpfile): - helpString += teal.getHelpFileAsString(__taskname__,__file__) + helpString += teal.getHelpFileAsString(__taskname__, __file__) else: - helpString += eval('.'.join([__package__,__taskname__,'__doc__'])) + helpString += eval('.'.join([__package__, __taskname__, '__doc__'])) else: - helpString = 'file://'+htmlfile + helpString = 'file://' + htmlfile return helpString def run(configObj=None): if configObj['hdrname'] == '' and configObj['hdrext'] is None and \ - configObj['distname'] == '': - print('='*60) + configObj['distname'] == '': + print('=' * 60) print('ERROR:') print(' No valid "hdrname", "hdrext" or "distname" parameter value provided!') print(' Please restart this task and provide a value for one of these parameters.') - print('='*60) + print('=' * 60) return filename = parseinput.parseinput(configObj['filename'])[0] # Call function with properly interpreted input parameters # Syntax: delete_headerlet(filename, hdrname=None, hdrext=None, distname=None) headerlet.delete_headerlet(filename, - hdrname = configObj['hdrname'], - hdrext = configObj['hdrext'], - distname = configObj['distname'], - logging = configObj['logging']) - + hdrname=configObj['hdrname'], + hdrext=configObj['hdrext'], + distname=configObj['distname'], + logging=configObj['logging']) diff --git a/stwcs/gui/extract_headerlet.py b/stwcs/gui/extract_headerlet.py index 02ecd7a..1e88221 100644 --- a/stwcs/gui/extract_headerlet.py +++ b/stwcs/gui/extract_headerlet.py @@ -1,14 +1,14 @@ -from __future__ import print_function +from __future__ import absolute_import, division, print_function + import os from stsci.tools import teal - -import stwcs -from stwcs.wcsutil import headerlet +from .. import __version__ +from ..wcsutil import headerlet __taskname__ = __name__.split('.')[-1] # needed for help string __package__ = headerlet.__name__ -__version__ = stwcs.__version__ +# __version__ = stwcs.__version__ # #### Interfaces used by TEAL # @@ -17,28 +17,29 @@ def getHelpAsString(docstring=False): return useful help from a file in the script directory called __taskname__.help """ install_dir = os.path.dirname(__file__) - htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html') - helpfile = os.path.join(install_dir,__taskname__+'.help') + htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html') + helpfile = os.path.join(install_dir, __taskname__ + '.help') if docstring or (not docstring and not os.path.exists(htmlfile)): - helpString = __taskname__+' Version '+__version__+'\n\n' + helpString = __taskname__ + ' Version ' + __version__ + '\n\n' if os.path.exists(helpfile): - helpString += teal.getHelpFileAsString(__taskname__,__file__) + helpString += teal.getHelpFileAsString(__taskname__, __file__) else: - helpString += eval('.'.join([__package__,__taskname__,'__doc__'])) + helpString += eval('.'.join([__package__, __taskname__, '__doc__'])) else: - helpString = 'file://'+htmlfile + helpString = 'file://' + htmlfile return helpString + def run(configObj=None): - if configObj['output'] in ['',' ',None]: - print('='*60) + if configObj['output'] in ['', ' ', None]: + print('=' * 60) print('ERROR:') print(' No valid "output" parameter value provided!') print(' Please restart this task and provide a value for this parameter.') - print('='*60) + print('=' * 60) return # create dictionary of remaining parameters, deleting extraneous ones @@ -55,4 +56,3 @@ def run(configObj=None): # clobber=False, verbose=100) headerlet.extract_headerlet(configObj['filename'], configObj['output'], **cdict) - diff --git a/stwcs/gui/headerlet_summary.py b/stwcs/gui/headerlet_summary.py index 82a3e0c..eac59ee 100644 --- a/stwcs/gui/headerlet_summary.py +++ b/stwcs/gui/headerlet_summary.py @@ -1,12 +1,13 @@ +from __future__ import absolute_import, division, print_function import os from stsci.tools import teal -import stwcs -from stwcs.wcsutil import headerlet +from .. import __version__ +from ..wcsutil import headerlet __taskname__ = __name__.split('.')[-1] # needed for help string __package__ = headerlet.__name__ -__version__ = stwcs.__version__ +# __version__ = stwcs.__version__ # #### Interfaces used by TEAL # @@ -15,19 +16,20 @@ def getHelpAsString(docstring=False): return useful help from a file in the script directory called __taskname__.help """ install_dir = os.path.dirname(__file__) - htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html') - helpfile = os.path.join(install_dir,__taskname__+'.help') + htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html') + helpfile = os.path.join(install_dir, __taskname__ + '.help') if docstring or (not docstring and not os.path.exists(htmlfile)): - helpString = __taskname__+' Version '+__version__+'\n\n' + helpString = __taskname__ + ' Version ' + __version__ + '\n\n' if os.path.exists(helpfile): - helpString += teal.getHelpFileAsString(__taskname__,__file__) + helpString += teal.getHelpFileAsString(__taskname__, __file__) else: - helpString += eval('.'.join([__package__,__taskname__,'__doc__'])) + helpString += eval('.'.join([__package__, __taskname__, '__doc__'])) else: - helpString = 'file://'+htmlfile + helpString = 'file://' + htmlfile return helpString + def run(configObj=None): # create dictionary of remaining parameters, deleting extraneous ones @@ -43,5 +45,4 @@ def run(configObj=None): # Syntax: headerlet_summary(filename,columns=None,pad=2,maxwidth=None, # output=None,clobber=True,quiet=False) - headerlet.headerlet_summary(configObj['filename'],**cdict) - + headerlet.headerlet_summary(configObj['filename'], **cdict) diff --git a/stwcs/gui/restore_headerlet.py b/stwcs/gui/restore_headerlet.py index 7570d76..790c239 100644 --- a/stwcs/gui/restore_headerlet.py +++ b/stwcs/gui/restore_headerlet.py @@ -1,13 +1,14 @@ +from __future__ import absolute_import, division, print_function import os from stsci.tools import teal -import stwcs -from stwcs.wcsutil import headerlet +from .. import __version__ +from ..wcsutil import headerlet __taskname__ = __name__.split('.')[-1] # needed for help string __package__ = headerlet.__name__ -__version__ = stwcs.__version__ +#__version__ = stwcs.__version__ # #### Interfaces used by TEAL # @@ -16,33 +17,33 @@ def getHelpAsString(docstring=False): return useful help from a file in the script directory called __taskname__.help """ install_dir = os.path.dirname(__file__) - htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html') - helpfile = os.path.join(install_dir,__taskname__+'.help') + htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html') + helpfile = os.path.join(install_dir, __taskname__ + '.help') if docstring or (not docstring and not os.path.exists(htmlfile)): - helpString = __taskname__+' Version '+__version__+'\n\n' + helpString = __taskname__ + ' Version ' + __version__ + '\n\n' if os.path.exists(helpfile): - helpString += teal.getHelpFileAsString(__taskname__,__file__) + helpString += teal.getHelpFileAsString(__taskname__, __file__) else: - helpString = 'file://'+htmlfile + helpString = 'file://' + htmlfile return helpString + def run(configObj=None): - if configObj['distname'] not in ['',' ','INDEF']: + if configObj['distname'] not in ['', ' ', 'INDEF']: # Call function with properly interpreted input parameters # Syntax: restore_all_with_distname(filename, distname, primary, # archive=True, sciext='SCI', verbose=False) headerlet.restore_all_with_distname(configObj['filename'], - configObj['distname'],configObj['primary'], - archive=configObj['archive'],sciext=configObj['sciext'], - logging=configObj['logging']) + configObj['distname'], configObj['primary'], + archive=configObj['archive'], sciext=configObj['sciext'], + logging=configObj['logging']) else: # Call function with properly interpreted input parameters # restore_from_headerlet(filename, hdrname=None, hdrext=None, # archive=True, force=False) headerlet.restore_from_headerlet(configObj['filename'], - hdrname=configObj['hdrname'],hdrext=configObj['hdrext'], - archive=configObj['archive'], force=configObj['force'], - logging=configObj['logging']) - + hdrname=configObj['hdrname'], hdrext=configObj['hdrext'], + archive=configObj['archive'], force=configObj['force'], + logging=configObj['logging']) diff --git a/stwcs/gui/updatewcs.py b/stwcs/gui/updatewcs.py index 3dacb67..b5dc700 100644 --- a/stwcs/gui/updatewcs.py +++ b/stwcs/gui/updatewcs.py @@ -1,19 +1,22 @@ -from __future__ import print_function +from __future__ import absolute_import, division, print_function + import os from astropy.io import fits from stsci.tools import parseinput from stsci.tools import fileutil from stsci.tools import teal -import stwcs -from stwcs import updatewcs -from stwcs.wcsutil import convertwcs +from .. import __version__ +from .. import updatewcs + + +allowed_corr_dict = {'vacorr': 'VACorr', 'tddcorr': 'TDDCorr', 'npolcorr': 'NPOLCorr', + 'd2imcorr': 'DET2IMCorr'} -allowed_corr_dict = {'vacorr':'VACorr','tddcorr':'TDDCorr','npolcorr':'NPOLCorr','d2imcorr':'DET2IMCorr'} __taskname__ = __name__.split('.')[-1] # needed for help string __package__ = updatewcs.__name__ -__version__ = stwcs.__version__ +#__version__ = stwcs.__version__ # #### Interfaces used by TEAL @@ -23,24 +26,24 @@ def getHelpAsString(docstring=False): return useful help from a file in the script directory called __taskname__.help """ install_dir = os.path.dirname(__file__) - htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html') - helpfile = os.path.join(install_dir,__taskname__+'.help') + htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html') + helpfile = os.path.join(install_dir, __taskname__ + '.help') if docstring or (not docstring and not os.path.exists(htmlfile)): - helpString = __taskname__+' Version '+__version__+'\n\n' + helpString = __taskname__ + ' Version ' + __version__ + '\n\n' if os.path.exists(helpfile): - helpString += teal.getHelpFileAsString(__taskname__,__file__) + helpString += teal.getHelpFileAsString(__taskname__, __file__) else: - helpString += eval('.'.join([__package__,__taskname__,'__doc__'])) + helpString += eval('.'.join([__package__, __taskname__, '__doc__'])) else: - helpString = 'file://'+htmlfile + helpString = 'file://' + htmlfile return helpString + def run(configObj=None): # Interpret primary parameters from configObj instance - extname = configObj['extname'] input = configObj['input'] # create dictionary of remaining parameters, deleting extraneous ones @@ -53,7 +56,7 @@ def run(configObj=None): del cdict['extname'] # parse input - input,altfiles = parseinput.parseinput(configObj['input']) + input, altfiles = parseinput.parseinput(configObj['input']) # Insure that all input files have a correctly archived # set of OPUS WCS keywords @@ -69,7 +72,7 @@ def run(configObj=None): # Check to insure that there is a valid reference file to be used idctab = fits.getval(file, 'idctab') if not os.path.exists(fileutil.osfn(idctab)): - print('No valid distortion reference file ',idctab,' found in ',file,'!') + print('No valid distortion reference file ', idctab, ' found in ', file, '!') raise ValueError # Re-define 'cdict' to only have switches for steps supported by that instrument @@ -78,7 +81,7 @@ def run(configObj=None): # for file in input: # get instrument name from input file - instr = fits.getval(file,'INSTRUME') + instr = fits.getval(file, 'INSTRUME') # make copy of input parameters dict for this file fdict = cdict.copy() # Remove any parameter that is not part of this instrument's allowed corrections @@ -86,5 +89,4 @@ def run(configObj=None): if allowed_corr_dict[step] not in updatewcs.apply_corrections.allowed_corrections[instr]: fdict[step] # Call 'updatewcs' on correctly archived file - updatewcs.updatewcs(file,**fdict) - + updatewcs.updatewcs(file, **fdict) diff --git a/stwcs/gui/write_headerlet.py b/stwcs/gui/write_headerlet.py index e18bed8..56520b3 100644 --- a/stwcs/gui/write_headerlet.py +++ b/stwcs/gui/write_headerlet.py @@ -1,15 +1,15 @@ -from __future__ import print_function +from __future__ import absolute_import, division, print_function import os from stsci.tools import teal from stsci.tools import parseinput -import stwcs -from stwcs.wcsutil import headerlet +from .. import __version__ +from ..wcsutil import headerlet __taskname__ = __name__.split('.')[-1] # needed for help string __package__ = headerlet.__name__ -__version__ = stwcs.__version__ +#__version__ = stwcs.__version__ # #### Interfaces used by TEAL # @@ -18,43 +18,44 @@ def getHelpAsString(docstring=False): return useful help from a file in the script directory called __taskname__.help """ install_dir = os.path.dirname(__file__) - htmlfile = os.path.join(install_dir,'htmlhelp',__taskname__+'.html') - helpfile = os.path.join(install_dir,__taskname__+'.help') + htmlfile = os.path.join(install_dir, 'htmlhelp', __taskname__ + '.html') + helpfile = os.path.join(install_dir, __taskname__ + '.help') if docstring or (not docstring and not os.path.exists(htmlfile)): - helpString = __taskname__+' Version '+__version__+'\n\n' + helpString = __taskname__ + ' Version ' + __version__ + '\n\n' if os.path.exists(helpfile): - helpString += teal.getHelpFileAsString(__taskname__,__file__) + helpString += teal.getHelpFileAsString(__taskname__, __file__) else: helpString += headerlet.write_headerlet.__doc__ else: - helpString = 'file://'+htmlfile + helpString = 'file://' + htmlfile return helpString + def run(configObj=None): - flist,oname = parseinput.parseinput(configObj['filename']) + flist, oname = parseinput.parseinput(configObj['filename']) if len(flist) == 0: - print('='*60) + print('=' * 60) print('ERROR:') print(' No valid "filename" parameter value provided!') print(' Please check the working directory and restart this task.') - print('='*60) + print('=' * 60) return - if configObj['hdrname'] in ['',' ','INDEF']: - print('='*60) + if configObj['hdrname'] in ['', ' ', 'INDEF']: + print('=' * 60) print('ERROR:') print(' No valid "hdrname" parameter value provided!') print(' Please restart this task and provide a value for this parameter.') - print('='*60) + print('=' * 60) return - if configObj['output'] in ['',' ','INDEF']: + if configObj['output'] in ['', ' ', 'INDEF']: configObj['output'] = None - str_kw = ['wcsname','destim','sipname','npolfile','d2imfile', - 'descrip','history','author','output','catalog'] + str_kw = ['wcsname', 'destim', 'sipname', 'npolfile', 'd2imfile', + 'descrip', 'history', 'author', 'output', 'catalog'] # create dictionary of remaining parameters, deleting extraneous ones # such as those above diff --git a/stwcs/tests/test_altwcs.py b/stwcs/tests/test_altwcs.py index 86d100f..d646106 100644 --- a/stwcs/tests/test_altwcs.py +++ b/stwcs/tests/test_altwcs.py @@ -1,11 +1,14 @@ +from __future__ import absolute_import, division, print_function + import shutil import os from astropy.io import fits as pyfits -from stwcs.wcsutil import altwcs -from stwcs import updatewcs -from stwcs.wcsutil import HSTWCS +from ..wcsutil import altwcs +from .. import updatewcs +from .. import wcsutil import numpy as np from numpy.testing import utils +import pytest from . import data data_path = os.path.split(os.path.abspath(data.__file__))[0] @@ -70,25 +73,25 @@ class TestAltWCS(object): updatewcs.updatewcs(acs_file) self.acs_file = acs_file self.simplefits = simple_file - self.ww = HSTWCS(self.acs_file, ext=1) + self.ww = wcsutil.HSTWCS(self.acs_file, ext=1) def test_archive(self): altwcs.archiveWCS(self.acs_file, ext=1, wcskey='Z', wcsname='ZTEST', reusekey=False) - w1 = HSTWCS(self.acs_file, ext=1) - w1z = HSTWCS(self.acs_file, ext=1, wcskey='Z') + w1 = wcsutil.HSTWCS(self.acs_file, ext=1) + w1z = wcsutil.HSTWCS(self.acs_file, ext=1, wcskey='Z') compare_wcs(w1, w1z) def test_archive_clobber(self): altwcs.archiveWCS(self.acs_file, ext=1, wcskey='Z', wcsname='ZTEST', reusekey=True) - w1 = HSTWCS(self.acs_file, ext=1) - w1z = HSTWCS(self.acs_file, ext=1, wcskey='Z') + w1 = wcsutil.HSTWCS(self.acs_file, ext=1) + w1z = wcsutil.HSTWCS(self.acs_file, ext=1, wcskey='Z') compare_wcs(w1, w1z) def test_restore_wcs(self): # test restore on a file altwcs.restoreWCS(self.acs_file, ext=1, wcskey='O') - w1o = HSTWCS(self.acs_file, ext=1, wcskey='O') - w1 = HSTWCS(self.acs_file, ext=1) + w1o = wcsutil.HSTWCS(self.acs_file, ext=1, wcskey='O') + w1 = wcsutil.HSTWCS(self.acs_file, ext=1) compare_wcs(w1, w1o, exclude_keywords=['ctype']) def test_restore_wcs_mem(self): @@ -99,8 +102,8 @@ class TestAltWCS(object): f = pyfits.open(self.acs_file, mode='update') altwcs.restoreWCS(f, ext=1, wcskey='T') f.close() - w1o = HSTWCS(self.acs_file, ext=1, wcskey='T') - w1 = HSTWCS(self.acs_file, ext=1) + w1o = wcsutil.HSTWCS(self.acs_file, ext=1, wcskey='T') + w1 = wcsutil.HSTWCS(self.acs_file, ext=1) compare_wcs(w1, w1o) def test_restore_simple(self): @@ -108,8 +111,8 @@ class TestAltWCS(object): altwcs.archiveWCS(self.simplefits, ext=0, wcskey='R') pyfits.setval(self.simplefits, ext=0, keyword='CRVAL1R', value=1) altwcs.restoreWCS(self.simplefits, ext=0, wcskey='R') - wo = HSTWCS(self.simplefits, ext=0, wcskey='R') - ws = HSTWCS(self.simplefits, ext=0) + wo = wcsutil.HSTWCS(self.simplefits, ext=0, wcskey='R') + ws = wcsutil.HSTWCS(self.simplefits, ext=0) compare_wcs(ws, wo) def test_restore_wcs_from_to(self): @@ -121,25 +124,27 @@ class TestAltWCS(object): altwcs.restore_from_to(f, fromext='SCI', toext=['SCI', 'ERR', 'DQ'], wcskey='T') f.close() - w1o = HSTWCS(self.acs_file, ext=('SCI', 1), wcskey='T') - w1 = HSTWCS(self.acs_file, ext=('SCI', 1)) + w1o = wcsutil.HSTWCS(self.acs_file, ext=('SCI', 1), wcskey='T') + w1 = wcsutil.HSTWCS(self.acs_file, ext=('SCI', 1)) compare_wcs(w1, w1o) - w2 = HSTWCS(self.acs_file, ext=('ERR', 1)) + w2 = wcsutil.HSTWCS(self.acs_file, ext=('ERR', 1)) compare_wcs(w2, w1o, exclude_keywords=['ctype']) - w3 = HSTWCS(self.acs_file, ext=('DQ', 1)) + w3 = wcsutil.HSTWCS(self.acs_file, ext=('DQ', 1)) compare_wcs(w3, w1o, exclude_keywords=['ctype']) - w4o = HSTWCS(self.acs_file, ext=4, wcskey='T') - w4 = HSTWCS(self.acs_file, ext=('SCI', 2)) + w4o = wcsutil.HSTWCS(self.acs_file, ext=4, wcskey='T') + w4 = wcsutil.HSTWCS(self.acs_file, ext=('SCI', 2)) compare_wcs(w4, w4o) - w5 = HSTWCS(self.acs_file, ext=('ERR', 2)) + w5 = wcsutil.HSTWCS(self.acs_file, ext=('ERR', 2)) compare_wcs(w5, w4o, exclude_keywords=['ctype']) - w6 = HSTWCS(self.acs_file, ext=('DQ', 2)) + w6 = wcsutil.HSTWCS(self.acs_file, ext=('DQ', 2)) compare_wcs(w3, w1o, exclude_keywords=['ctype']) def test_delete_wcs(self): #altwcs.archiveWCS(self.acs_file, ext=1, wcskey='Z') altwcs.deleteWCS(self.acs_file, ext=1, wcskey='Z') - utils.assert_raises(KeyError, HSTWCS, self.acs_file, ext=1, wcskey='Z') + #utils.assert_raises(KeyError, wcsutil.HSTWCS, self.acs_file, ext=1, wcskey='Z') + with pytest.raises(KeyError): + wcsutil.HSTWCS(self.acs_file, ext=1, wcskey='Z') def test_pars_file_mode1(self): assert(not altwcs._parpasscheck(self.acs_file, ext=1, wcskey='Z')) diff --git a/stwcs/tests/test_headerlet.py b/stwcs/tests/test_headerlet.py index 48fb470..f4fa050 100644 --- a/stwcs/tests/test_headerlet.py +++ b/stwcs/tests/test_headerlet.py @@ -1,9 +1,11 @@ +from __future__ import absolute_import, division, print_function + import shutil import os from astropy.io import fits -from stwcs import updatewcs -from stwcs.wcsutil import headerlet, wcsdiff -from stwcs.wcsutil import HSTWCS +from .. import updatewcs +from ..wcsutil import headerlet, wcsdiff +from ..wcsutil import HSTWCS import numpy as np from numpy.testing import utils from nose.tools import * diff --git a/stwcs/tests/test_updatewcs.py b/stwcs/tests/test_updatewcs.py index dc5a34f..b7e75c8 100644 --- a/stwcs/tests/test_updatewcs.py +++ b/stwcs/tests/test_updatewcs.py @@ -1,12 +1,14 @@ +from __future__ import absolute_import, division, print_function + import shutil import os from astropy import wcs from astropy.io import fits -from stwcs import updatewcs -from stwcs.updatewcs import apply_corrections -from stwcs.distortion import utils as dutils -from stwcs.wcsutil import HSTWCS +from .. import updatewcs +from ..updatewcs import apply_corrections +from ..distortion import utils as dutils +from ..wcsutil import HSTWCS import numpy as np from numpy.testing import utils import pytest @@ -111,7 +113,7 @@ class TestStwcs(object): #print('outwcs.wcs.crval = {0}'.format(outwcs.wcs.crval)) utils.assert_allclose( - outwcs.wcs.crval, np.array([5.65109952, -72.0674181]), rtol=1e-7) + outwcs.wcs.crval, np.array([5.65109952, -72.0674181]), atol=1e-7) utils.assert_almost_equal(outwcs.wcs.crpix, np.array([2107.0, 2118.5])) utils.assert_almost_equal( diff --git a/stwcs/updatewcs/__init__.py b/stwcs/updatewcs/__init__.py index eb5e850..a59f106 100644 --- a/stwcs/updatewcs/__init__.py +++ b/stwcs/updatewcs/__init__.py @@ -1,14 +1,16 @@ -from __future__ import absolute_import, division, print_function # confidence high +from __future__ import absolute_import, division, print_function +import atexit from astropy.io import fits -from stwcs import wcsutil -from stwcs.wcsutil import HSTWCS -import stwcs +from .. import wcsutil +#from ..wcsutil.hwstwcs import HSTWCS + +from .. import __version__ from astropy import wcs as pywcs import astropy -from . import utils, corrections, makewcs +from . import utils, corrections from . import npol, det2im from stsci.tools import parseinput, fileutil from . import apply_corrections @@ -17,10 +19,10 @@ import time import logging logger = logging.getLogger('stwcs.updatewcs') -import atexit atexit.register(logging.shutdown) -#Note: The order of corrections is important +# Note: The order of corrections is important + def updatewcs(input, vacorr=True, tddcorr=True, npolcorr=True, d2imcorr=True, checkfiles=True, verbose=False): @@ -62,7 +64,7 @@ def updatewcs(input, vacorr=True, tddcorr=True, npolcorr=True, d2imcorr=True, geis and waiver fits files will be converted to MEF format. Default value is True for standalone mode. """ - if verbose == False: + if not verbose: logger.setLevel(100) else: formatter = logging.Formatter("%(asctime)s - %(name)s - %(levelname)s - %(message)s") @@ -74,12 +76,12 @@ def updatewcs(input, vacorr=True, tddcorr=True, npolcorr=True, d2imcorr=True, logger.setLevel(verbose) args = "vacorr=%s, tddcorr=%s, npolcorr=%s, d2imcorr=%s, checkfiles=%s, \ " % (str(vacorr), str(tddcorr), str(npolcorr), - str(d2imcorr), str(checkfiles)) + str(d2imcorr), str(checkfiles)) logger.info('\n\tStarting UPDATEWCS: %s', time.asctime()) files = parseinput.parseinput(input)[0] logger.info("\n\tInput files: %s, " % [i for i in files]) - logger.info("\n\tInput arguments: %s" %args) + logger.info("\n\tInput arguments: %s" % args) if checkfiles: files = checkFiles(files) if not files: @@ -87,8 +89,8 @@ def updatewcs(input, vacorr=True, tddcorr=True, npolcorr=True, d2imcorr=True, return for f in files: - acorr = apply_corrections.setCorrections(f, vacorr=vacorr, \ - tddcorr=tddcorr,npolcorr=npolcorr, d2imcorr=d2imcorr) + acorr = apply_corrections.setCorrections(f, vacorr=vacorr, tddcorr=tddcorr, + npolcorr=npolcorr, d2imcorr=d2imcorr) if 'MakeWCS' in acorr and newIDCTAB(f): logger.warning("\n\tNew IDCTAB file detected. All current WCSs will be deleted") cleanWCS(f) @@ -97,6 +99,7 @@ def updatewcs(input, vacorr=True, tddcorr=True, npolcorr=True, d2imcorr=True, return files + def makecorr(fname, allowed_corr): """ Purpose @@ -111,11 +114,11 @@ def makecorr(fname, allowed_corr): """ logger.info("Allowed corrections: {0}".format(allowed_corr)) f = fits.open(fname, mode='update') - #Determine the reference chip and create the reference HSTWCS object + # Determine the reference chip and create the reference HSTWCS object nrefchip, nrefext = getNrefchip(f) wcsutil.restoreWCS(f, nrefext, wcskey='O') - rwcs = HSTWCS(fobj=f, ext=nrefext) - rwcs.readModel(update=True,header=f[nrefext].header) + rwcs = wcsutil.HSTWCS(fobj=f, ext=nrefext) + rwcs.readModel(update=True, header=f[nrefext].header) if 'DET2IMCorr' in allowed_corr: kw2update = det2im.DET2IMCorr.updateWCS(f) @@ -127,16 +130,16 @@ def makecorr(fname, allowed_corr): if 'extname' in extn.header: extname = extn.header['extname'].lower() - if extname == 'sci': + if extname == 'sci': wcsutil.restoreWCS(f, ext=i, wcskey='O') sciextver = extn.header['extver'] ref_wcs = rwcs.deepcopy() hdr = extn.header - ext_wcs = HSTWCS(fobj=f, ext=i) - ### check if it exists first!!! + ext_wcs = wcsutil.HSTWCS(fobj=f, ext=i) + # check if it exists first!!! # 'O ' can be safely archived again because it has been restored first. wcsutil.archiveWCS(f, ext=i, wcskey="O", wcsname="OPUS", reusekey=True) - ext_wcs.readModel(update=True,header=hdr) + ext_wcs.readModel(update=True, header=hdr) for c in allowed_corr: if c != 'NPOLCorr' and c != 'DET2IMCorr': corr_klass = corrections.__getattribute__(c) @@ -147,14 +150,14 @@ def makecorr(fname, allowed_corr): idcname = f[0].header.get('IDCTAB', " ") if idcname.strip() and 'idc.fits' in idcname: wname = ''.join(['IDC_', - utils.extract_rootname(idcname,suffix='_idc')]) + utils.extract_rootname(idcname, suffix='_idc')]) else: wname = " " hdr['WCSNAME'] = wname elif extname in ['err', 'dq', 'sdq', 'samp', 'time']: cextver = extn.header['extver'] if cextver == sciextver: - hdr = f[('SCI',sciextver)].header + hdr = f[('SCI', sciextver)].header w = pywcs.WCS(hdr, f) copyWCS(w, extn.header) @@ -167,14 +170,14 @@ def makecorr(fname, allowed_corr): f[1].header[kw] = kw2update[kw] # Finally record the version of the software which updated the WCS if 'HISTORY' in f[0].header: - f[0].header.set('UPWCSVER', value=stwcs.__version__, + f[0].header.set('UPWCSVER', value=__version__, comment="Version of STWCS used to updated the WCS", before='HISTORY') f[0].header.set('PYWCSVER', value=astropy.__version__, comment="Version of PYWCS used to updated the WCS", before='HISTORY') elif 'ASN_MTYP' in f[0].header: - f[0].header.set('UPWCSVER', value=stwcs.__version__, + f[0].header.set('UPWCSVER', value=__version__, comment="Version of STWCS used to updated the WCS", after='ASN_MTYP') f[0].header.set('PYWCSVER', value=astropy.__version__, @@ -185,20 +188,21 @@ def makecorr(fname, allowed_corr): for i in range(len(f[0].header) - 1, 0, -1): if f[0].header[i].strip() != '': break - f[0].header.set('UPWCSVER', stwcs.__version__, + f[0].header.set('UPWCSVER', __version__, "Version of STWCS used to updated the WCS", after=i) f[0].header.set('PYWCSVER', astropy.__version__, "Version of PYWCS used to updated the WCS", after=i) # add additional keywords to be used by headerlets - distdict = utils.construct_distname(f,rwcs) + distdict = utils.construct_distname(f, rwcs) f[0].header['DISTNAME'] = distdict['DISTNAME'] f[0].header['SIPNAME'] = distdict['SIPNAME'] # Make sure NEXTEND keyword remains accurate - f[0].header['NEXTEND'] = len(f)-1 + f[0].header['NEXTEND'] = len(f) - 1 f.close() + def copyWCS(w, ehdr): """ This is a convenience function to copy a WCS object @@ -214,6 +218,7 @@ def copyWCS(w, ehdr): key = k[:7] ehdr[key] = hwcs[k] + def getNrefchip(fobj): """ @@ -235,15 +240,15 @@ def getNrefchip(fobj): if instrument == 'WFPC2': chipkw = 'DETECTOR' - extvers = [("SCI",img.header['EXTVER']) for img in - fobj[1:] if img.header['EXTNAME'].lower()=='sci'] + extvers = [("SCI", img.header['EXTVER']) for img in + fobj[1:] if img.header['EXTNAME'].lower() == 'sci'] detectors = [img.header[chipkw] for img in fobj[1:] if - img.header['EXTNAME'].lower()=='sci'] + img.header['EXTNAME'].lower() == 'sci'] fitsext = [i for i in range(len(fobj))[1:] if - fobj[i].header['EXTNAME'].lower()=='sci'] - det2ext=dict(list(zip(detectors, extvers))) - ext2det=dict(list(zip(extvers, detectors))) - ext2fitsext=dict(list(zip(extvers, fitsext))) + fobj[i].header['EXTNAME'].lower() == 'sci'] + det2ext = dict(list(zip(detectors, extvers))) + ext2det = dict(list(zip(extvers, detectors))) + ext2fitsext = dict(list(zip(extvers, fitsext))) if 3 not in detectors: nrefchip = ext2det.pop(extvers[0]) @@ -256,15 +261,15 @@ def getNrefchip(fobj): elif (instrument == 'ACS' and fobj[0].header['DETECTOR'] == 'WFC') or \ (instrument == 'WFC3' and fobj[0].header['DETECTOR'] == 'UVIS'): chipkw = 'CCDCHIP' - extvers = [("SCI",img.header['EXTVER']) for img in - fobj[1:] if img.header['EXTNAME'].lower()=='sci'] + extvers = [("SCI", img.header['EXTVER']) for img in + fobj[1:] if img.header['EXTNAME'].lower() == 'sci'] detectors = [img.header[chipkw] for img in fobj[1:] if - img.header['EXTNAME'].lower()=='sci'] + img.header['EXTNAME'].lower() == 'sci'] fitsext = [i for i in range(len(fobj))[1:] if - fobj[i].header['EXTNAME'].lower()=='sci'] - det2ext=dict(list(zip(detectors, extvers))) - ext2det=dict(list(zip(extvers, detectors))) - ext2fitsext=dict(list(zip(extvers, fitsext))) + fobj[i].header['EXTNAME'].lower() == 'sci'] + det2ext = dict(list(zip(detectors, extvers))) + ext2det = dict(list(zip(extvers, detectors))) + ext2fitsext = dict(list(zip(extvers, fitsext))) if 2 not in detectors: nrefchip = ext2det.pop(extvers[0]) @@ -276,12 +281,13 @@ def getNrefchip(fobj): else: for i in range(len(fobj)): extname = fobj[i].header.get('EXTNAME', None) - if extname != None and extname.lower == 'sci': + if extname is not None and extname.lower == 'sci': nrefext = i break return nrefchip, nrefext + def checkFiles(input): """ Checks that input files are in the correct format. @@ -292,13 +298,14 @@ def checkFiles(input): removed_files = [] newfiles = [] if not isinstance(input, list): - input=[input] + input = [input] for file in input: try: - imgfits,imgtype = fileutil.isFits(file) + imgfits, imgtype = fileutil.isFits(file) except IOError: - logger.warning( "\n\tFile %s could not be found, removing it from the input list.\n" %file) + logger.warning("File {0} could not be found, removing it from the" + "input list.".format(file)) removed_files.append(file) continue # Check for existence of waiver FITS input, and quit if found. @@ -306,8 +313,9 @@ def checkFiles(input): if imgfits: if imgtype == 'waiver': newfilename = waiver2mef(file, convert_dq=True) - if newfilename == None: - logger.warning("\n\tRemoving file %s from input list - could not convert waiver to mef" %file) + if newfilename is None: + logger.warning("Could not convert waiver to mef." + "Removing file {0} from input list".format(file)) removed_files.append(file) else: newfiles.append(newfilename) @@ -319,24 +327,26 @@ def checkFiles(input): # Convert the corresponding data quality file if present if not imgfits: newfilename = geis2mef(file, convert_dq=True) - if newfilename == None: - logger.warning("\n\tRemoving file %s from input list - could not convert geis to mef" %file) + if newfilename is None: + logger.warning("Could not convert file {0} from geis to mef, removing it from input list".format(file)) removed_files.append(file) else: newfiles.append(newfilename) if removed_files: - logger.warning('\n\tThe following files will be removed from the list of files to be processed %s' % removed_files) + logger.warning('\n\tThe following files will be removed from the list of files', + 'to be processed %s' % removed_files) newfiles = checkFiles(newfiles)[0] logger.info("\n\tThese files passed the input check and will be processed: %s" % newfiles) return newfiles + def newIDCTAB(fname): - #When this is called we know there's a kw IDCTAB in the header + # When this is called we know there's a kw IDCTAB in the header hdul = fits.open(fname) idctab = fileutil.osfn(hdul[0].header['IDCTAB']) try: - #check for the presence of IDCTAB in the first extension + # check for the presence of IDCTAB in the first extension oldidctab = fileutil.osfn(hdul[1].header['IDCTAB']) except KeyError: return False @@ -345,6 +355,7 @@ def newIDCTAB(fname): else: return True + def cleanWCS(fname): # A new IDCTAB means all previously computed WCS's are invalid # We are deleting all of them except the original OPUS WCS. @@ -363,6 +374,7 @@ def cleanWCS(fname): # Some extensions don't have the alternate (or any) WCS keywords continue + def getCorrections(instrument): """ Print corrections available for an instrument @@ -373,4 +385,4 @@ def getCorrections(instrument): acorr = apply_corrections.allowed_corrections[instrument] print("The following corrections will be performed for instrument %s\n" % instrument) - for c in acorr: print(c,': ' , apply_corrections.cnames[c]) + for c in acorr: print(c, ': ', apply_corrections.cnames[c]) diff --git a/stwcs/updatewcs/apply_corrections.py b/stwcs/updatewcs/apply_corrections.py index 86b623a..4a8e5ac 100644 --- a/stwcs/updatewcs/apply_corrections.py +++ b/stwcs/updatewcs/apply_corrections.py @@ -1,37 +1,35 @@ -from __future__ import division, print_function # confidence high +from __future__ import absolute_import, division, print_function -import os +import os.path from astropy.io import fits -import time from stsci.tools import fileutil -import os.path -from stwcs.wcsutil import altwcs from . import utils from . import wfpc2_dgeo import logging logger = logging.getLogger("stwcs.updatewcs.apply_corrections") -#Note: The order of corrections is important +# Note: The order of corrections is important # A dictionary which lists the allowed corrections for each instrument. # These are the default corrections applied also in the pipeline. -allowed_corrections={'WFPC2': ['DET2IMCorr', 'MakeWCS','CompSIP', 'VACorr'], - 'ACS': ['DET2IMCorr', 'TDDCorr', 'MakeWCS', 'CompSIP','VACorr', 'NPOLCorr'], - 'STIS': ['MakeWCS', 'CompSIP','VACorr'], - 'NICMOS': ['MakeWCS', 'CompSIP','VACorr'], - 'WFC3': ['DET2IMCorr', 'MakeWCS', 'CompSIP','VACorr', 'NPOLCorr'], - } +allowed_corrections = {'WFPC2': ['DET2IMCorr', 'MakeWCS', 'CompSIP', 'VACorr'], + 'ACS': ['DET2IMCorr', 'TDDCorr', 'MakeWCS', 'CompSIP', 'VACorr', 'NPOLCorr'], + 'STIS': ['MakeWCS', 'CompSIP', 'VACorr'], + 'NICMOS': ['MakeWCS', 'CompSIP', 'VACorr'], + 'WFC3': ['DET2IMCorr', 'MakeWCS', 'CompSIP', 'VACorr', 'NPOLCorr'], + } cnames = {'DET2IMCorr': 'Detector to Image Correction', - 'TDDCorr': 'Time Dependent Distortion Correction', - 'MakeWCS': 'Recalculate basic WCS keywords based on the distortion model', - 'CompSIP': 'Given IDCTAB distortion model calculate the SIP coefficients', - 'VACorr': 'Velocity Aberration Correction', - 'NPOLCorr': 'Lookup Table Distortion' - } + 'TDDCorr': 'Time Dependent Distortion Correction', + 'MakeWCS': 'Recalculate basic WCS keywords based on the distortion model', + 'CompSIP': 'Given IDCTAB distortion model calculate the SIP coefficients', + 'VACorr': 'Velocity Aberration Correction', + 'NPOLCorr': 'Lookup Table Distortion' + } + def setCorrections(fname, vacorr=True, tddcorr=True, npolcorr=True, d2imcorr=True): """ @@ -40,7 +38,7 @@ def setCorrections(fname, vacorr=True, tddcorr=True, npolcorr=True, d2imcorr=Tru for the instrument. """ instrument = fits.getval(fname, 'INSTRUME') - # make a copy of this list ! + # make a copy of this list acorr = allowed_corrections[instrument][:] # For WFPC2 images, the old-style DGEOFILE needs to be @@ -56,18 +54,22 @@ def setCorrections(fname, vacorr=True, tddcorr=True, npolcorr=True, d2imcorr=Tru if 'MakeWCS' in acorr: acorr.remove('MakeWCS') if 'CompSIP' in acorr: acorr.remove('CompSIP') - if 'VACorr' in acorr and vacorr==False: acorr.remove('VACorr') + if 'VACorr' in acorr and not vacorr: + acorr.remove('VACorr') if 'TDDCorr' in acorr: tddcorr = applyTDDCorr(fname, tddcorr) - if tddcorr == False: acorr.remove('TDDCorr') + if not tddcorr: + acorr.remove('TDDCorr') if 'NPOLCorr' in acorr: npolcorr = applyNpolCorr(fname, npolcorr) - if npolcorr == False: acorr.remove('NPOLCorr') + if not npolcorr: + acorr.remove('NPOLCorr') if 'DET2IMCorr' in acorr: d2imcorr = applyD2ImCorr(fname, d2imcorr) - if d2imcorr == False: acorr.remove('DET2IMCorr') - logger.info("\n\tCorrections to be applied to %s: %s " % (fname, str(acorr))) + if not d2imcorr: + acorr.remove('DET2IMCorr') + logger.info("Corrections to be applied to {0} {1}".format(fname, acorr)) return acorr @@ -119,22 +121,21 @@ def applyTDDCorr(fname, utddcorr): except KeyError: tddswitch = 'PERFORM' - if instrument == 'ACS' and detector == 'WFC' and utddcorr == True and tddswitch == 'PERFORM': + if instrument == 'ACS' and detector == 'WFC' and utddcorr and tddswitch == 'PERFORM': tddcorr = True try: idctab = phdr['IDCTAB'] except KeyError: tddcorr = False - #print "***IDCTAB keyword not found - not applying TDD correction***\n" if os.path.exists(fileutil.osfn(idctab)): tddcorr = True else: tddcorr = False - #print "***IDCTAB file not found - not applying TDD correction***\n" else: tddcorr = False return tddcorr + def applyNpolCorr(fname, unpolcorr): """ Determines whether non-polynomial distortion lookup tables should be added @@ -156,9 +157,8 @@ def applyNpolCorr(fname, unpolcorr): return False fnpol0 = fileutil.osfn(fnpol0) if not fileutil.findFile(fnpol0): - msg = """\n\tKw "NPOLFILE" exists in primary header but file %s not found - Non-polynomial distortion correction will not be applied\n - """ % fnpol0 + msg = '"NPOLFILE" exists in primary header but file {0} not found.' + 'Non-polynomial distortion correction will not be applied.'.format(file) logger.critical(msg) raise IOError("NPOLFILE {0} not found".format(fnpol0)) try: @@ -176,7 +176,7 @@ def applyNpolCorr(fname, unpolcorr): else: # npl file defined in first extension may not be found # but if a valid kw exists in the primary header, non-polynomial - #distortion correction should be applied. + # distortion correction should be applied. applyNPOLCorr = True except KeyError: # the case of "NPOLFILE" kw present in primary header but "NPOLEXT" missing @@ -191,6 +191,7 @@ def applyNpolCorr(fname, unpolcorr): applyNPOLCorr = False return (applyNPOLCorr and unpolcorr) + def isOldStyleDGEO(fname, dgname): # checks if the file defined in a NPOLFILE kw is a full size # (old style) image @@ -209,6 +210,7 @@ def isOldStyleDGEO(fname, dgname): else: return False + def applyD2ImCorr(fname, d2imcorr): applyD2IMCorr = True try: diff --git a/stwcs/updatewcs/corrections.py b/stwcs/updatewcs/corrections.py index d3641eb..975a806 100644 --- a/stwcs/updatewcs/corrections.py +++ b/stwcs/updatewcs/corrections.py @@ -1,8 +1,9 @@ -from __future__ import division, print_function # confidence high +from __future__ import absolute_import, division, print_function import copy import datetime -import logging, time +import logging +import time import numpy as np from numpy import linalg from stsci.tools import fileutil @@ -11,11 +12,12 @@ from . import npol from . import makewcs from .utils import diff_angles -logger=logging.getLogger('stwcs.updatewcs.corrections') +logger = logging.getLogger('stwcs.updatewcs.corrections') MakeWCS = makewcs.MakeWCS NPOLCorr = npol.NPOLCorr + class TDDCorr(object): """ Apply time dependent distortion correction to distortion coefficients and basic @@ -70,39 +72,39 @@ class TDDCorr(object): ext_wcs.idcmodel.ocx = copy.deepcopy(ext_wcs.idcmodel.cx) ext_wcs.idcmodel.ocy = copy.deepcopy(ext_wcs.idcmodel.cy) - newkw = {'TDDALPHA': None, 'TDDBETA':None, - 'OCX10':ext_wcs.idcmodel.ocx[1,0], - 'OCX11':ext_wcs.idcmodel.ocx[1,1], - 'OCY10':ext_wcs.idcmodel.ocy[1,0], - 'OCY11':ext_wcs.idcmodel.ocy[1,1], - 'TDD_CTA':None, 'TDD_CTB':None, - 'TDD_CYA':None, 'TDD_CYB':None, - 'TDD_CXA':None, 'TDD_CXB':None} - - if ext_wcs.idcmodel.refpix['skew_coeffs']is not None and \ - ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CTB'] is not None: + newkw = {'TDDALPHA': None, 'TDDBETA': None, + 'OCX10': ext_wcs.idcmodel.ocx[1, 0], + 'OCX11': ext_wcs.idcmodel.ocx[1, 1], + 'OCY10': ext_wcs.idcmodel.ocy[1, 0], + 'OCY11': ext_wcs.idcmodel.ocy[1, 1], + 'TDD_CTA': None, 'TDD_CTB': None, + 'TDD_CYA': None, 'TDD_CYB': None, + 'TDD_CXA': None, 'TDD_CXB': None + } + + if ext_wcs.idcmodel.refpix['skew_coeffs'] is not None and \ + ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CTB'] is not None: cls.apply_tdd2idc2015(ref_wcs) cls.apply_tdd2idc2015(ext_wcs) - newkw.update({ - 'TDD_CTA':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CTA'], - 'TDD_CYA':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CYA'], - 'TDD_CXA':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CXA'], - 'TDD_CTB':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CTB'], - 'TDD_CYB':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CYB'], - 'TDD_CXB':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CXB']}) - - elif ext_wcs.idcmodel.refpix['skew_coeffs']is not None and \ - ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CY_BETA'] is not None: - logger.info("\n\t Applying 2014-calibrated TDD: %s" % time.asctime()) + newkw.update({'TDD_CTA': ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CTA'], + 'TDD_CYA': ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CYA'], + 'TDD_CXA': ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CXA'], + 'TDD_CTB': ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CTB'], + 'TDD_CYB': ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CYB'], + 'TDD_CXB': ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CXB'] + }) + + elif ext_wcs.idcmodel.refpix['skew_coeffs'] is not None and \ + ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CY_BETA'] is not None: + logger.info("Applying 2014-calibrated TDD: {0}".format(time.asctime())) # We have 2014-calibrated TDD, not J.A.-style TDD cls.apply_tdd2idc2(ref_wcs) cls.apply_tdd2idc2(ext_wcs) - newkw.update({'TDD_CYA':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CY_ALPHA'], - 'TDD_CYB':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CY_BETA'], - 'TDD_CXA':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CX_ALPHA'], - 'TDD_CXB':ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CX_BETA']}) - + newkw.update({'TDD_CYA': ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CY_ALPHA'], + 'TDD_CYB': ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CY_BETA'], + 'TDD_CXA': ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CX_ALPHA'], + 'TDD_CXB': ext_wcs.idcmodel.refpix['skew_coeffs']['TDD_CX_BETA']}) else: alpha, beta = cls.compute_alpha_beta(ext_wcs) cls.apply_tdd2idc(ref_wcs, alpha, beta) @@ -111,8 +113,7 @@ class TDDCorr(object): ext_wcs.idcmodel.refpix['TDDBETA'] = beta ref_wcs.idcmodel.refpix['TDDALPHA'] = alpha ref_wcs.idcmodel.refpix['TDDBETA'] = beta - - newkw.update( {'TDDALPHA': alpha, 'TDDBETA':beta} ) + newkw.update({'TDDALPHA': alpha, 'TDDBETA': beta} ) return newkw updateWCS = classmethod(updateWCS) @@ -121,10 +122,10 @@ class TDDCorr(object): """ Applies 2015-calibrated TDD correction to a couple of IDCTAB coefficients for ACS/WFC observations. """ - if not isinstance(hwcs.date_obs,float): - year,month,day = hwcs.date_obs.split('-') - rdate = datetime.datetime(int(year),int(month),int(day)) - rday = float(rdate.strftime("%j"))/365.25 + rdate.year + if not isinstance(hwcs.date_obs, float): + year, month, day = hwcs.date_obs.split('-') + rdate = datetime.datetime(int(year), int(month), int(day)) + rday = float(rdate.strftime("%j")) / 365.25 + rdate.year else: rday = hwcs.date_obs @@ -132,12 +133,11 @@ class TDDCorr(object): delta_date = rday - skew_coeffs['TDD_DATE'] if skew_coeffs['TDD_CXB'] is not None: - hwcs.idcmodel.cx[1,1] += skew_coeffs['TDD_CXB']*delta_date + hwcs.idcmodel.cx[1, 1] += skew_coeffs['TDD_CXB'] * delta_date if skew_coeffs['TDD_CTB'] is not None: - hwcs.idcmodel.cy[1,1] += skew_coeffs['TDD_CTB']*delta_date + hwcs.idcmodel.cy[1, 1] += skew_coeffs['TDD_CTB'] * delta_date if skew_coeffs['TDD_CYB'] is not None: - hwcs.idcmodel.cy[1,0] += skew_coeffs['TDD_CYB']*delta_date - #print("CX[1,1]_TDD={}, CY[1,1]_TDD={}, CY[1,0]_TDD={}".format(hwcs.idcmodel.cx[1,1],hwcs.idcmodel.cy[1,1],hwcs.idcmodel.cy[1,0])) + hwcs.idcmodel.cy[1, 0] += skew_coeffs['TDD_CYB'] * delta_date apply_tdd2idc2015 = classmethod(apply_tdd2idc2015) @@ -145,10 +145,10 @@ class TDDCorr(object): """ Applies 2014-calibrated TDD correction to single IDCTAB coefficient of an ACS/WFC observation. """ - if not isinstance(hwcs.date_obs,float): - year,month,day = hwcs.date_obs.split('-') - rdate = datetime.datetime(int(year),int(month),int(day)) - rday = float(rdate.strftime("%j"))/365.25 + rdate.year + if not isinstance(hwcs.date_obs, float): + year, month, day = hwcs.date_obs.split('-') + rdate = datetime.datetime(int(year), int(month), int(day)) + rday = float(rdate.strftime("%j")) / 365.25 + rdate.year else: rday = hwcs.date_obs @@ -156,21 +156,24 @@ class TDDCorr(object): cy_beta = skew_coeffs['TDD_CY_BETA'] cy_alpha = skew_coeffs['TDD_CY_ALPHA'] delta_date = rday - skew_coeffs['TDD_DATE'] - print("DELTA_DATE: {0} based on rday: {1}, TDD_DATE: {2}".format(delta_date,rday,skew_coeffs['TDD_DATE'])) + logger.info("DELTA_DATE: {0} based on rday: {1}, TDD_DATE: {2}".format(delta_date, rday, + skew_coeffs['TDD_DATE'])) if cy_alpha is None: - hwcs.idcmodel.cy[1,1] += cy_beta*delta_date + hwcs.idcmodel.cy[1, 1] += cy_beta * delta_date else: - new_beta = cy_alpha + cy_beta*delta_date - hwcs.idcmodel.cy[1,1] = new_beta - print("CY11: {0} based on alpha: {1}, beta: {2}".format(hwcs.idcmodel.cy[1,1],cy_alpha,cy_beta)) - + new_beta = cy_alpha + cy_beta * delta_date + hwcs.idcmodel.cy[1, 1] = new_beta + logger.info("CY11: {0} based on alpha: {1}, beta: {2}".format(hwcs.idcmodel.cy[1, 1], + cy_alpha, cy_beta)) + cx_beta = skew_coeffs['TDD_CX_BETA'] cx_alpha = skew_coeffs['TDD_CX_ALPHA'] if cx_alpha is not None: - new_beta = cx_alpha + cx_beta*delta_date - hwcs.idcmodel.cx[1,1] = new_beta - print("CX11: {0} based on alpha: {1}, beta: {2}".format(new_beta,cx_alpha,cx_beta)) + new_beta = cx_alpha + cx_beta * delta_date + hwcs.idcmodel.cx[1, 1] = new_beta + logger.info("CX11: {0} based on alpha: {1}, beta: {2}".format(new_beta, + cx_alpha, cx_beta)) apply_tdd2idc2 = classmethod(apply_tdd2idc2) @@ -182,11 +185,12 @@ class TDDCorr(object): theta_v2v3 = 2.234529 mrotp = fileutil.buildRotMatrix(theta_v2v3) mrotn = fileutil.buildRotMatrix(-theta_v2v3) - tdd_mat = np.array([[1+(beta/2048.), alpha/2048.],[alpha/2048.,1-(beta/2048.)]],np.float64) + tdd_mat = np.array([[1 + (beta / 2048.), alpha / 2048.], + [alpha / 2048., 1 - (beta / 2048.)]], np.float64) abmat1 = np.dot(tdd_mat, mrotn) - abmat2 = np.dot(mrotp,abmat1) + abmat2 = np.dot(mrotp, abmat1) xshape, yshape = hwcs.idcmodel.cx.shape, hwcs.idcmodel.cy.shape - icxy = np.dot(abmat2,[hwcs.idcmodel.cx.ravel(), hwcs.idcmodel.cy.ravel()]) + icxy = np.dot(abmat2, [hwcs.idcmodel.cx.ravel(), hwcs.idcmodel.cy.ravel()]) hwcs.idcmodel.cx = icxy[0] hwcs.idcmodel.cy = icxy[1] hwcs.idcmodel.cx.shape = xshape @@ -209,10 +213,10 @@ class TDDCorr(object): alpha = 0.095 + 0.090*(rday-dday)/2.5 beta = -0.029 - 0.030*(rday-dday)/2.5 """ - if not isinstance(ext_wcs.date_obs,float): - year,month,day = ext_wcs.date_obs.split('-') - rdate = datetime.datetime(int(year),int(month),int(day)) - rday = float(rdate.strftime("%j"))/365.25 + rdate.year + if not isinstance(ext_wcs.date_obs, float): + year, month, day = ext_wcs.date_obs.split('-') + rdate = datetime.datetime(int(year), int(month), int(day)) + rday = float(rdate.strftime("%j")) / 365.25 + rdate.year else: rday = ext_wcs.date_obs @@ -220,26 +224,27 @@ class TDDCorr(object): if skew_coeffs is None: # Only print out warning for post-SM4 data where this may matter if rday > 2009.0: - err_str = "------------------------------------------------------------------------ \n" + err_str = "------------------------------------------------------------------------ \n" err_str += "WARNING: the IDCTAB geometric distortion file specified in the image \n" err_str += " header did not have the time-dependent distortion coefficients. \n" err_str += " The pre-SM4 time-dependent skew solution will be used by default.\n" err_str += " Please update IDCTAB with new reference file from HST archive. \n" - err_str += "------------------------------------------------------------------------ \n" + err_str += "------------------------------------------------------------------------ \n" print(err_str) # Using default pre-SM4 coefficients - skew_coeffs = {'TDD_A':[0.095,0.090/2.5], - 'TDD_B':[-0.029,-0.030/2.5], - 'TDD_DATE':2004.5,'TDDORDER':1} + skew_coeffs = {'TDD_A': [0.095, 0.090 / 2.5], + 'TDD_B': [-0.029, -0.030 / 2.5], + 'TDD_DATE': 2004.5, + 'TDDORDER': 1} alpha = 0 beta = 0 # Compute skew terms, allowing for non-linear coefficients as well - for c in range(skew_coeffs['TDDORDER']+1): - alpha += skew_coeffs['TDD_A'][c]* np.power((rday-skew_coeffs['TDD_DATE']),c) - beta += skew_coeffs['TDD_B'][c]*np.power((rday-skew_coeffs['TDD_DATE']),c) + for c in range(skew_coeffs['TDDORDER'] + 1): + alpha += skew_coeffs['TDD_A'][c] * np.power((rday - skew_coeffs['TDD_DATE']), c) + beta += skew_coeffs['TDD_B'][c] * np.power((rday - skew_coeffs['TDD_DATE']), c) - return alpha,beta + return alpha, beta compute_alpha_beta = classmethod(compute_alpha_beta) @@ -254,21 +259,21 @@ class VACorr(object): """ def updateWCS(cls, ext_wcs, ref_wcs): - logger.info("\n\tStarting VACorr: %s" % time.asctime()) + logger.info("Starting VACorr: %s" % time.asctime()) if ext_wcs.vafactor != 1: ext_wcs.wcs.cd = ext_wcs.wcs.cd * ext_wcs.vafactor - crval0 = ref_wcs.wcs.crval[0] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[0], - ref_wcs.wcs.crval[0]) - crval1 = ref_wcs.wcs.crval[1] + ext_wcs.vafactor*diff_angles(ext_wcs.wcs.crval[1], - ref_wcs.wcs.crval[1]) - crval = np.array([crval0,crval1]) + crval0 = ref_wcs.wcs.crval[0] + ext_wcs.vafactor * diff_angles(ext_wcs.wcs.crval[0], + ref_wcs.wcs.crval[0]) + crval1 = ref_wcs.wcs.crval[1] + ext_wcs.vafactor * diff_angles(ext_wcs.wcs.crval[1], + ref_wcs.wcs.crval[1]) + crval = np.array([crval0, crval1]) ext_wcs.wcs.crval = crval ext_wcs.wcs.set() else: pass - kw2update={'CD1_1': ext_wcs.wcs.cd[0,0], 'CD1_2':ext_wcs.wcs.cd[0,1], - 'CD2_1':ext_wcs.wcs.cd[1,0], 'CD2_2':ext_wcs.wcs.cd[1,1], - 'CRVAL1':ext_wcs.wcs.crval[0], 'CRVAL2':ext_wcs.wcs.crval[1]} + kw2update = {'CD1_1': ext_wcs.wcs.cd[0, 0], 'CD1_2': ext_wcs.wcs.cd[0, 1], + 'CD2_1': ext_wcs.wcs.cd[1, 0], 'CD2_2': ext_wcs.wcs.cd[1, 1], + 'CRVAL1': ext_wcs.wcs.crval[0], 'CRVAL2': ext_wcs.wcs.crval[1]} return kw2update updateWCS = classmethod(updateWCS) @@ -291,34 +296,34 @@ class CompSIP(object): """ def updateWCS(cls, ext_wcs, ref_wcs): - logger.info("\n\tStarting CompSIP: %s" %time.asctime()) + logger.info("Starting CompSIP: {0}".format(time.asctime())) kw2update = {} if not ext_wcs.idcmodel: - logger.info("IDC model not found, SIP coefficient will not be computed") + logger.info("IDC model not found, SIP coefficient will not be computed.") return kw2update order = ext_wcs.idcmodel.norder kw2update['A_ORDER'] = order kw2update['B_ORDER'] = order - pscale = ext_wcs.idcmodel.refpix['PSCALE'] + # pscale = ext_wcs.idcmodel.refpix['PSCALE'] cx = ext_wcs.idcmodel.cx cy = ext_wcs.idcmodel.cy - matr = np.array([[cx[1,1],cx[1,0]], [cy[1,1],cy[1,0]]], dtype=np.float64) + matr = np.array([[cx[1, 1], cx[1, 0]], [cy[1, 1], cy[1, 0]]], dtype=np.float64) imatr = linalg.inv(matr) - akeys1 = np.zeros((order+1,order+1), dtype=np.float64) - bkeys1 = np.zeros((order+1,order+1), dtype=np.float64) - for n in range(order+1): - for m in range(order+1): - if n >= m and n>=2: - idcval = np.array([[cx[n,m]],[cy[n,m]]]) + akeys1 = np.zeros((order + 1, order + 1), dtype=np.float64) + bkeys1 = np.zeros((order + 1, order + 1), dtype=np.float64) + for n in range(order + 1): + for m in range(order + 1): + if n >= m and n >= 2: + idcval = np.array([[cx[n, m]], [cy[n, m]]]) sipval = np.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] * ext_wcs.binned - kw2update[Bkey] = sipval[1,0] * ext_wcs.binned + 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] * ext_wcs.binned + kw2update[Bkey] = sipval[1, 0] * ext_wcs.binned kw2update['CTYPE1'] = 'RA---TAN-SIP' kw2update['CTYPE2'] = 'DEC--TAN-SIP' return kw2update diff --git a/stwcs/updatewcs/det2im.py b/stwcs/updatewcs/det2im.py index bddb37b..bb4cd9c 100644 --- a/stwcs/updatewcs/det2im.py +++ b/stwcs/updatewcs/det2im.py @@ -1,14 +1,13 @@ -from __future__ import absolute_import, division # confidence high +from __future__ import absolute_import, division, print_function -import numpy as np from astropy.io import fits from stsci.tools import fileutil -from . import utils - -import logging, time +import logging +import time logger = logging.getLogger('stwcs.updatewcs.d2im') + class DET2IMCorr(object): """ Defines a Lookup table prior distortion correction as per WCS paper IV. @@ -38,11 +37,11 @@ class DET2IMCorr(object): Science file, for which a distortion correction in a NPOLFILE is available """ - logger.info("\n\tStarting DET2IM: %s" %time.asctime()) + logger.info("Starting DET2IM: {0}".format(time.asctime())) try: assert isinstance(fobj, fits.HDUList) except AssertionError: - logger.exception('\n\tInput must be a fits.HDUList object') + logger.exception('Input must be a fits.HDUList object') raise cls.applyDet2ImCorr(fobj) @@ -84,7 +83,7 @@ class DET2IMCorr(object): hdu = cls.createD2ImHDU(header, d2imfile=d2imfile, wdvarr_ver=d2im_num_ext, d2im_extname=ename[0], - data=ename[1],ccdchip=ccdchip) + data=ename[1], ccdchip=ccdchip) if wcsdvarr_ind and d2im_num_ext in wcsdvarr_ind: fobj[wcsdvarr_ind[d2im_num_ext]] = hdu else: @@ -108,7 +107,7 @@ class DET2IMCorr(object): continue if ename == 'D2IMARR': wcsd[fobj[e].header['EXTVER']] = e - logger.debug("A map of D2IMARR extensions %s" % wcsd) + logger.debug("A map of D2IMARR extensions {0}".format(wcsd)) return wcsd getWCSIndex = classmethod(getWCSIndex) @@ -118,17 +117,17 @@ class DET2IMCorr(object): Adds kw to sci extension to define WCSDVARR lookup table extensions """ - if d2im_extname =='DX': - j=1 + if d2im_extname == 'DX': + j = 1 else: - j=2 - - d2imerror = 'D2IMERR%s' %j - d2imdis = 'D2IMDIS%s' %j - d2imext = 'D2IM%s.' %j + 'EXTVER' - d2imnaxes = 'D2IM%s.' %j +'NAXES' - d2imaxis1 = 'D2IM%s.' %j+'AXIS.1' - d2imaxis2 = 'D2IM%s.' %j+'AXIS.2' + j = 2 + + d2imerror = 'D2IMERR%s' % j + d2imdis = 'D2IMDIS%s' % j + d2imext = 'D2IM%s.' % j + 'EXTVER' + d2imnaxes = 'D2IM%s.' % j + 'NAXES' + d2imaxis1 = 'D2IM%s.' % j + 'AXIS.1' + d2imaxis2 = 'D2IM%s.' % j + 'AXIS.2' keys = [d2imerror, d2imdis, d2imext, d2imnaxes, d2imaxis1, d2imaxis2] values = {d2imerror: error_val, d2imdis: 'Lookup', @@ -154,7 +153,7 @@ class DET2IMCorr(object): addSciExtKw = classmethod(addSciExtKw) - def getData(cls,d2imfile, ccdchip): + def getData(cls, d2imfile, ccdchip): """ Get the data arrays from the reference D2I files Make sure 'CCDCHIP' in the npolfile matches "CCDCHIP' in the science file. @@ -162,8 +161,8 @@ class DET2IMCorr(object): xdata, ydata = (None, None) d2im = fits.open(d2imfile) for ext in d2im: - d2imextname = ext.header.get('EXTNAME',"") - d2imccdchip = ext.header.get('CCDCHIP',1) + d2imextname = ext.header.get('EXTNAME', "") + d2imccdchip = ext.header.get('CCDCHIP', 1) if d2imextname == 'DX' and d2imccdchip == ccdchip: xdata = ext.data.copy() continue @@ -177,7 +176,7 @@ class DET2IMCorr(object): getData = classmethod(getData) def createD2ImHDU(cls, sciheader, d2imfile=None, wdvarr_ver=1, - d2im_extname=None,data = None, ccdchip=1): + d2im_extname=None, data=None, ccdchip=1): """ Creates an HDU to be added to the file object. """ @@ -216,26 +215,26 @@ class DET2IMCorr(object): naxis = d2im[1].header['NAXIS'] ccdchip = d2imextname - kw = { 'NAXIS': 'Size of the axis', - 'CDELT': 'Coordinate increment along axis', - 'CRPIX': 'Coordinate system reference pixel', - 'CRVAL': 'Coordinate system value at reference pixel', - } + kw = {'NAXIS': 'Size of the axis', + 'CDELT': 'Coordinate increment along axis', + 'CRPIX': 'Coordinate system reference pixel', + 'CRVAL': 'Coordinate system value at reference pixel', + } kw_comm1 = {} kw_val1 = {} for key in kw.keys(): - for i in range(1, naxis+1): + for i in range(1, naxis + 1): si = str(i) - kw_comm1[key+si] = kw[key] + kw_comm1[key + si] = kw[key] - for i in range(1, naxis+1): + for i in range(1, naxis + 1): si = str(i) - kw_val1['NAXIS'+si] = d2im_header.get('NAXIS'+si) - kw_val1['CDELT'+si] = d2im_header.get('CDELT'+si, 1.0) * sciheader.get('LTM'+si+'_'+si, 1) - kw_val1['CRPIX'+si] = d2im_header.get('CRPIX'+si, 0.0) - kw_val1['CRVAL'+si] = (d2im_header.get('CRVAL'+si, 0.0) - \ - sciheader.get('LTV'+str(i), 0)) + kw_val1['NAXIS' + si] = d2im_header.get('NAXIS' + si) + kw_val1['CDELT' + si] = d2im_header.get('CDELT' + si, 1.0) * sciheader.get('LTM' + si + '_' + si, 1) + kw_val1['CRPIX' + si] = d2im_header.get('CRPIX' + si, 0.0) + kw_val1['CRVAL' + si] = (d2im_header.get('CRVAL' + si, 0.0) - + sciheader.get('LTV' + str(i), 0)) kw_comm0 = {'XTENSION': 'Image extension', 'BITPIX': 'IEEE floating point', 'NAXIS': 'Number of axes', @@ -245,15 +244,15 @@ class DET2IMCorr(object): 'GCOUNT': 'One data group', } - kw_val0 = { 'XTENSION': 'IMAGE', - 'BITPIX': -32, - 'NAXIS': naxis, - 'EXTNAME': 'D2IMARR', - 'EXTVER': wdvarr_ver, - 'PCOUNT': 0, - 'GCOUNT': 1, - 'CCDCHIP': ccdchip, - } + kw_val0 = {'XTENSION': 'IMAGE', + 'BITPIX': -32, + 'NAXIS': naxis, + 'EXTNAME': 'D2IMARR', + 'EXTVER': wdvarr_ver, + 'PCOUNT': 0, + 'GCOUNT': 1, + 'CCDCHIP': ccdchip, + } cdl = [] for key in kw_comm0.keys(): cdl.append((key, kw_val0[key], kw_comm0[key])) @@ -266,8 +265,8 @@ class DET2IMCorr(object): for i, c in enumerate(d2im_phdr): if c == 'FILENAME': start_indx = i - if c == '': # remove blanks from end of header - end_indx = i+1 + if c == '': # remove blanks from end of header + end_indx = i + 1 break if start_indx >= 0: for card in d2im_phdr.cards[start_indx:end_indx]: @@ -287,7 +286,7 @@ class DET2IMCorr(object): if fobj[0].header['INSTRUME'] == 'ACS' and fobj[0].header['DETECTOR'] == 'WFC': ccdchip = fobj[extname, extver].header['CCDCHIP'] elif fobj[0].header['INSTRUME'] == 'WFC3' and fobj[0].header['DETECTOR'] == 'UVIS': - ccdchip = fobj[extname, extver].header['CCDCHIP'] + ccdchip = fobj[extname, extver].header['CCDCHIP'] elif fobj[0].header['INSTRUME'] == 'WFPC2': ccdchip = fobj[extname, extver].header['DETECTOR'] elif fobj[0].header['INSTRUME'] == 'STIS': diff --git a/stwcs/updatewcs/makewcs.py b/stwcs/updatewcs/makewcs.py index 06c6f9c..4c6b3df 100644 --- a/stwcs/updatewcs/makewcs.py +++ b/stwcs/updatewcs/makewcs.py @@ -1,16 +1,16 @@ -from __future__ import absolute_import, division # confidence high - -import datetime +from __future__ import absolute_import, division, print_function import numpy as np -import logging, time -from math import sin, sqrt, pow, cos, asin, atan2,pi +import logging +import time +from math import sin, sqrt, pow, cos, asin, atan2, pi from stsci.tools import fileutil from . import utils logger = logging.getLogger(__name__) + class MakeWCS(object): """ Recompute basic WCS keywords based on PA_V3 and distortion model. @@ -36,7 +36,8 @@ class MakeWCS(object): - Time dependent distortion correction is applied to both chips' V2/V3 values. """ - tdd_xyref = {1: [2048, 3072], 2:[2048, 1024]} + tdd_xyref = {1: [2048, 3072], 2: [2048, 1024]} + def updateWCS(cls, ext_wcs, ref_wcs): """ recomputes the basic WCS kw @@ -53,20 +54,20 @@ class MakeWCS(object): cls.uprefwcs(ext_wcs, ref_wcs, rv23_corr, ltvoff, offshift) cls.upextwcs(ext_wcs, ref_wcs, v23_corr, rv23_corr, ltvoff, offshift) - kw2update = {'CD1_1': ext_wcs.wcs.cd[0,0], - 'CD1_2': ext_wcs.wcs.cd[0,1], - 'CD2_1': ext_wcs.wcs.cd[1,0], - 'CD2_2': ext_wcs.wcs.cd[1,1], - 'CRVAL1': ext_wcs.wcs.crval[0], - 'CRVAL2': ext_wcs.wcs.crval[1], - 'CRPIX1': ext_wcs.wcs.crpix[0], - 'CRPIX2': ext_wcs.wcs.crpix[1], - 'IDCTAB': ext_wcs.idctab, - 'OCX10' : ext_wcs.idcmodel.cx[1,0], - 'OCX11' : ext_wcs.idcmodel.cx[1,1], - 'OCY10' : ext_wcs.idcmodel.cy[1,0], - 'OCY11' : ext_wcs.idcmodel.cy[1,1] - } + kw2update = {'CD1_1': ext_wcs.wcs.cd[0, 0], + 'CD1_2': ext_wcs.wcs.cd[0, 1], + 'CD2_1': ext_wcs.wcs.cd[1, 0], + 'CD2_2': ext_wcs.wcs.cd[1, 1], + 'CRVAL1': ext_wcs.wcs.crval[0], + 'CRVAL2': ext_wcs.wcs.crval[1], + 'CRPIX1': ext_wcs.wcs.crpix[0], + 'CRPIX2': ext_wcs.wcs.crpix[1], + 'IDCTAB': ext_wcs.idctab, + 'OCX10': ext_wcs.idcmodel.cx[1, 0], + 'OCX11': ext_wcs.idcmodel.cx[1, 1], + 'OCY10': ext_wcs.idcmodel.cy[1, 0], + 'OCY11': ext_wcs.idcmodel.cy[1, 1] + } return kw2update updateWCS = classmethod(updateWCS) @@ -82,41 +83,42 @@ class MakeWCS(object): if ltv1 != 0. or ltv2 != 0.: offsetx = ext_wcs.wcs.crpix[0] - ltv1 - ext_wcs.idcmodel.refpix['XREF'] offsety = ext_wcs.wcs.crpix[1] - ltv2 - ext_wcs.idcmodel.refpix['YREF'] - ext_wcs.idcmodel.shift(offsetx,offsety) + ext_wcs.idcmodel.shift(offsetx, offsety) fx, fy = ext_wcs.idcmodel.cx, ext_wcs.idcmodel.cy ext_wcs.setPscale() - tddscale = (ref_wcs.pscale/fx[1,1]) - v2 = ext_wcs.idcmodel.refpix['V2REF'] + v23_corr[0,0] * tddscale - v3 = ext_wcs.idcmodel.refpix['V3REF'] - v23_corr[1,0] * tddscale - v2ref = ref_wcs.idcmodel.refpix['V2REF'] + rv23_corr[0,0] * tddscale - v3ref = ref_wcs.idcmodel.refpix['V3REF'] - rv23_corr[1,0] * tddscale + tddscale = (ref_wcs.pscale / fx[1, 1]) + v2 = ext_wcs.idcmodel.refpix['V2REF'] + v23_corr[0, 0] * tddscale + v3 = ext_wcs.idcmodel.refpix['V3REF'] - v23_corr[1, 0] * tddscale + v2ref = ref_wcs.idcmodel.refpix['V2REF'] + rv23_corr[0, 0] * tddscale + v3ref = ref_wcs.idcmodel.refpix['V3REF'] - rv23_corr[1, 0] * tddscale - R_scale = ref_wcs.idcmodel.refpix['PSCALE']/3600.0 - off = sqrt((v2-v2ref)**2 + (v3-v3ref)**2)/(R_scale*3600.0) + R_scale = ref_wcs.idcmodel.refpix['PSCALE'] / 3600.0 + off = sqrt((v2 - v2ref) ** 2 + (v3 - v3ref) ** 2) / (R_scale * 3600.0) if v3 == v3ref: - theta=0.0 + theta = 0.0 else: - theta = atan2(ext_wcs.parity[0][0]*(v2-v2ref), ext_wcs.parity[1][1]*(v3-v3ref)) + theta = atan2(ext_wcs.parity[0][0] * (v2 - v2ref), + ext_wcs.parity[1][1] * (v3 - v3ref)) - if ref_wcs.idcmodel.refpix['THETA']: theta += ref_wcs.idcmodel.refpix['THETA']*pi/180.0 + if ref_wcs.idcmodel.refpix['THETA']: theta += ref_wcs.idcmodel.refpix['THETA'] * pi / 180.0 - dX=(off*sin(theta)) + offshiftx - dY=(off*cos(theta)) + offshifty + dX = (off * sin(theta)) + offshiftx + dY = (off * cos(theta)) + offshifty - px = np.array([[dX,dY]]) + px = np.array([[dX, dY]]) newcrval = ref_wcs.wcs.p2s(px, 1)['world'][0] newcrpix = np.array([ext_wcs.idcmodel.refpix['XREF'] + ltvoffx, - ext_wcs.idcmodel.refpix['YREF'] + ltvoffy]) + ext_wcs.idcmodel.refpix['YREF'] + ltvoffy]) ext_wcs.wcs.crval = newcrval ext_wcs.wcs.crpix = newcrpix ext_wcs.wcs.set() # Create a small vector, in reference image pixel scale - delmat = np.array([[fx[1,1], fy[1,1]], \ - [fx[1,0], fy[1,0]]]) / R_scale/3600. + delmat = np.array([[fx[1, 1], fy[1, 1]], + [fx[1, 0], fy[1, 0]]]) / R_scale / 3600. # Account for subarray offset # Angle of chip relative to chip @@ -131,10 +133,10 @@ class MakeWCS(object): wc = ref_wcs.wcs.p2s((px + dxy), 1)['world'] # Calculate the new CDs and convert to degrees - cd11 = utils.diff_angles(wc[0,0],newcrval[0])*cos(newcrval[1]*pi/180.0) - cd12 = utils.diff_angles(wc[1,0],newcrval[0])*cos(newcrval[1]*pi/180.0) - cd21 = utils.diff_angles(wc[0,1],newcrval[1]) - cd22 = utils.diff_angles(wc[1,1],newcrval[1]) + cd11 = utils.diff_angles(wc[0, 0], newcrval[0]) * cos(newcrval[1] * pi / 180.0) + cd12 = utils.diff_angles(wc[1, 0], newcrval[0]) * cos(newcrval[1] * pi / 180.0) + cd21 = utils.diff_angles(wc[0, 1], newcrval[1]) + cd22 = utils.diff_angles(wc[1, 1], newcrval[1]) cd = np.array([[cd11, cd12], [cd21, cd22]]) ext_wcs.wcs.cd = cd ext_wcs.wcs.set() @@ -151,39 +153,38 @@ class MakeWCS(object): if ref_wcs.ltv1 != 0. or ref_wcs.ltv2 != 0.: offsetx = ref_wcs.wcs.crpix[0] - ltv1 - ref_wcs.idcmodel.refpix['XREF'] offsety = ref_wcs.wcs.crpix[1] - ltv2 - ref_wcs.idcmodel.refpix['YREF'] - ref_wcs.idcmodel.shift(offsetx,offsety) + ref_wcs.idcmodel.shift(offsetx, offsety) - rfx, rfy = ref_wcs.idcmodel.cx, ref_wcs.idcmodel.cy + # rfx, rfy = ref_wcs.idcmodel.cx, ref_wcs.idcmodel.cy - offshift = offsh + # offshift = offsh dec = ref_wcs.wcs.crval[1] - tddscale = (ref_wcs.pscale/ref_wcs.idcmodel.cx[1,1]) - rv23 = [ref_wcs.idcmodel.refpix['V2REF'] + (rv23_corr_tdd[0,0] *tddscale), - ref_wcs.idcmodel.refpix['V3REF'] - (rv23_corr_tdd[1,0] * tddscale)] + tddscale = (ref_wcs.pscale / ref_wcs.idcmodel.cx[1, 1]) + rv23 = [ref_wcs.idcmodel.refpix['V2REF'] + (rv23_corr_tdd[0, 0] * tddscale), + ref_wcs.idcmodel.refpix['V3REF'] - (rv23_corr_tdd[1, 0] * tddscale)] # Get an approximate reference position on the sky - rref = np.array([[ref_wcs.idcmodel.refpix['XREF']+ltvoffx , - ref_wcs.idcmodel.refpix['YREF']+ltvoffy]]) + rref = np.array([[ref_wcs.idcmodel.refpix['XREF'] + ltvoffx, + ref_wcs.idcmodel.refpix['YREF'] + ltvoffy]]) crval = ref_wcs.wcs.p2s(rref, 1)['world'][0] # 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(ext_wcs.pav3,dec,rv23[0], rv23[1]) + pv = troll(ext_wcs.pav3, dec, rv23[0], rv23[1]) # Add the chip rotation angle if ref_wcs.idcmodel.refpix['THETA']: pv += ref_wcs.idcmodel.refpix['THETA'] - # Set values for the rest of the reference WCS ref_wcs.wcs.crval = crval - ref_wcs.wcs.crpix = np.array([0.0,0.0])+offsh + ref_wcs.wcs.crpix = np.array([0.0, 0.0]) + offsh parity = ref_wcs.parity - R_scale = ref_wcs.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 + R_scale = ref_wcs.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 = np.array([[cd11, cd12], [cd21, cd22]]) ref_wcs.wcs.cd = rcd @@ -191,9 +192,9 @@ class MakeWCS(object): uprefwcs = classmethod(uprefwcs) - def zero_point_corr(cls,hwcs): + def zero_point_corr(cls, hwcs): if hwcs.idcmodel.refpix['skew_coeffs'] is not None and 'TDD_CY_BETA' in hwcs.idcmodel.refpix['skew_coeffs']: - v23_corr = np.array([[0.],[0.]]) + v23_corr = np.array([[0.], [0.]]) return v23_corr else: try: @@ -202,15 +203,17 @@ class MakeWCS(object): except KeyError: alpha = 0.0 beta = 0.0 - v23_corr = np.array([[0.],[0.]]) - logger.debug("\n\tTDD Zero point correction for chip %s defaulted to: %s" % (hwcs.chip, v23_corr)) + v23_corr = np.array([[0.], [0.]]) + logger.debug("TDD Zero point correction for chip" + "{0} defaulted to: {1}".format(hwcs.chip, v23_corr)) return v23_corr tdd = np.array([[beta, alpha], [alpha, -beta]]) - mrotp = fileutil.buildRotMatrix(2.234529)/2048. - xy0 = np.array([[cls.tdd_xyref[hwcs.chip][0]-2048.], [cls.tdd_xyref[hwcs.chip][1]-2048.]]) - v23_corr = np.dot(mrotp,np.dot(tdd,xy0)) * 0.05 - logger.debug("\n\tTDD Zero point correction for chip %s: %s" % (hwcs.chip, v23_corr)) + mrotp = fileutil.buildRotMatrix(2.234529) / 2048. + xy0 = np.array([[cls.tdd_xyref[hwcs.chip][0] - 2048.], + [cls.tdd_xyref[hwcs.chip][1] - 2048.]]) + v23_corr = np.dot(mrotp, np.dot(tdd, xy0)) * 0.05 + logger.debug("TDD Zero point correction for chip {0}: {1}".format(hwcs.chip, v23_corr)) return v23_corr zero_point_corr = classmethod(zero_point_corr) @@ -258,16 +261,16 @@ def troll(roll, dec, v2, v3): _v3 = np.deg2rad(v3 / 3600.) # compute components - sin_rho = sqrt((pow(sin(_v2),2)+pow(sin(_v3),2)) - (pow(sin(_v2),2)*pow(sin(_v3),2))) + 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) + beta = asin(sin(_v3) / sin_rho) if _v2 < 0: beta = pi - beta - gamma = asin(sin(_v2)/sin_rho) + 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))) + 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 = np.rad2deg(pi - (gamma+B)) + troll = np.rad2deg(pi - (gamma + B)) return troll diff --git a/stwcs/updatewcs/npol.py b/stwcs/updatewcs/npol.py index 33579f0..9c88150 100644 --- a/stwcs/updatewcs/npol.py +++ b/stwcs/updatewcs/npol.py @@ -1,15 +1,16 @@ -from __future__ import absolute_import, division # confidence high +from __future__ import absolute_import, division, print_function -import logging, time +import logging +import time import numpy as np from astropy.io import fits from stsci.tools import fileutil -from . import utils logger = logging.getLogger('stwcs.updatewcs.npol') + class NPOLCorr(object): """ Defines a Lookup table prior distortion correction as per WCS paper IV. @@ -43,7 +44,7 @@ class NPOLCorr(object): Science file, for which a distortion correction in a NPOLFILE is available """ - logger.info("\n\tStarting NPOL: %s" %time.asctime()) + logger.info("\n\tStarting NPOL: %s" % time.asctime()) try: assert isinstance(fobj, fits.HDUList) except AssertionError: @@ -79,28 +80,28 @@ class NPOLCorr(object): header = ext.header # get the data arrays from the reference file and transform # them for use with SIP - dx,dy = cls.getData(nplfile, ccdchip) + dx, dy = cls.getData(nplfile, ccdchip) idccoeffs = cls.getIDCCoeffs(header) if idccoeffs is not None: - dx, dy = cls.transformData(dx,dy, idccoeffs) + dx, dy = cls.transformData(dx, dy, idccoeffs) # Determine EXTVER for the WCSDVARR extension from the # NPL file (EXTNAME, EXTVER) kw. # This is used to populate DPj.EXTVER kw - wcsdvarr_x_version = 2 * extversion -1 + wcsdvarr_x_version = 2 * extversion - 1 wcsdvarr_y_version = 2 * extversion - for ename in zip(['DX', 'DY'], [wcsdvarr_x_version,wcsdvarr_y_version],[dx, dy]): + for ename in zip(['DX', 'DY'], [wcsdvarr_x_version, wcsdvarr_y_version], [dx, dy]): error_val = ename[2].max() cls.addSciExtKw(header, wdvarr_ver=ename[1], npol_extname=ename[0], error_val=error_val) - hdu = cls.createNpolHDU(header, npolfile=nplfile, \ - wdvarr_ver=ename[1], npl_extname=ename[0], data=ename[2],ccdchip=ccdchip) + hdu = cls.createNpolHDU(header, npolfile=nplfile, + wdvarr_ver=ename[1], npl_extname=ename[0], + data=ename[2], ccdchip=ccdchip) if wcsdvarr_ind: fobj[wcsdvarr_ind[ename[1]]] = hdu else: fobj.append(hdu) - applyNPOLCorr = classmethod(applyNPOLCorr) def getWCSIndex(cls, fobj): @@ -129,17 +130,17 @@ class NPOLCorr(object): Adds kw to sci extension to define WCSDVARR lookup table extensions """ - if npol_extname =='DX': - j=1 + if npol_extname == 'DX': + j = 1 else: - j=2 - - cperror = 'CPERR%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' + j = 2 + + cperror = 'CPERR%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: error_val, cpdis: 'Lookup', @@ -165,15 +166,15 @@ class NPOLCorr(object): addSciExtKw = classmethod(addSciExtKw) - def getData(cls,nplfile, ccdchip): + def getData(cls, nplfile, ccdchip): """ Get the data arrays from the reference NPOL files Make sure 'CCDCHIP' in the npolfile matches "CCDCHIP' in the science file. """ npl = fits.open(nplfile) for ext in npl: - nplextname = ext.header.get('EXTNAME',"") - nplccdchip = ext.header.get('CCDCHIP',1) + nplextname = ext.header.get('EXTNAME', "") + nplccdchip = ext.header.get('CCDCHIP', 1) if nplextname == 'DX' and nplccdchip == ccdchip: xdata = ext.data.copy() continue @@ -192,7 +193,7 @@ class NPOLCorr(object): """ ndx, ndy = np.dot(coeffs, [dx.ravel(), dy.ravel()]).astype(np.float32) ndx.shape = dx.shape - ndy.shape=dy.shape + ndy.shape = dy.shape return ndx, ndy transformData = classmethod(transformData) @@ -206,7 +207,7 @@ class NPOLCorr(object): ocx11 = header['OCX11'] ocy10 = header['OCY10'] ocy11 = header['OCY11'] - coeffs = np.array([[ocx11, ocx10], [ocy11,ocy10]], dtype=np.float32) + coeffs = np.array([[ocx11, ocx10], [ocy11, ocy10]], dtype=np.float32) except KeyError: logger.exception('\n\tFirst order IDCTAB coefficients are not available. \n\ Cannot convert SIP to IDC coefficients.') @@ -217,15 +218,17 @@ class NPOLCorr(object): logger.exception("IDCSCALE not found in header - setting it to 1.") idcscale = 1 - return np.linalg.inv(coeffs/idcscale) + return np.linalg.inv(coeffs / idcscale) getIDCCoeffs = classmethod(getIDCCoeffs) - def createNpolHDU(cls, sciheader, npolfile=None, wdvarr_ver=1, npl_extname=None,data = None, ccdchip=1): + def createNpolHDU(cls, sciheader, npolfile=None, wdvarr_ver=1, npl_extname=None, + data=None, ccdchip=1): """ Creates an HDU to be added to the file object. """ - hdr = cls.createNpolHdr(sciheader, npolfile=npolfile, wdvarr_ver=wdvarr_ver, npl_extname=npl_extname, ccdchip=ccdchip) + hdr = cls.createNpolHdr(sciheader, npolfile=npolfile, wdvarr_ver=wdvarr_ver, + npl_extname=npl_extname, ccdchip=ccdchip) hdu = fits.ImageHDU(header=hdr, data=data) return hdu @@ -256,29 +259,29 @@ class NPOLCorr(object): npl.close() naxis = npl[1].header['NAXIS'] - ccdchip = nplextname #npol_header['CCDCHIP'] + ccdchip = nplextname # npol_header['CCDCHIP'] - kw = { 'NAXIS': 'Size of the axis', - 'CDELT': 'Coordinate increment along axis', - 'CRPIX': 'Coordinate system reference pixel', - 'CRVAL': 'Coordinate system value at reference pixel', - } + kw = {'NAXIS': 'Size of the axis', + 'CDELT': 'Coordinate increment along axis', + 'CRPIX': 'Coordinate system reference pixel', + 'CRVAL': 'Coordinate system value at reference pixel', + } kw_comm1 = {} kw_val1 = {} for key in kw.keys(): - for i in range(1, naxis+1): + for i in range(1, naxis + 1): si = str(i) - kw_comm1[key+si] = kw[key] + kw_comm1[key + si] = kw[key] - for i in range(1, naxis+1): + for i in range(1, naxis + 1): si = str(i) - kw_val1['NAXIS'+si] = npol_header.get('NAXIS'+si) - kw_val1['CDELT'+si] = npol_header.get('CDELT'+si, 1.0) * \ - sciheader.get('LTM'+si+'_'+si, 1) - kw_val1['CRPIX'+si] = npol_header.get('CRPIX'+si, 0.0) - kw_val1['CRVAL'+si] = (npol_header.get('CRVAL'+si, 0.0) - \ - sciheader.get('LTV'+str(i), 0)) + kw_val1['NAXIS' + si] = npol_header.get('NAXIS' + si) + kw_val1['CDELT' + si] = npol_header.get('CDELT' + si, 1.0) * \ + sciheader.get('LTM' + si + '_' + si, 1) + kw_val1['CRPIX' + si] = npol_header.get('CRPIX' + si, 0.0) + kw_val1['CRVAL' + si] = (npol_header.get('CRVAL' + si, 0.0) - + sciheader.get('LTV' + str(i), 0)) kw_comm0 = {'XTENSION': 'Image extension', 'BITPIX': 'IEEE floating point', @@ -289,15 +292,15 @@ class NPOLCorr(object): 'GCOUNT': 'One data group', } - kw_val0 = { 'XTENSION': 'IMAGE', - 'BITPIX': -32, - 'NAXIS': naxis, - 'EXTNAME': 'WCSDVARR', - 'EXTVER': wdvarr_ver, - 'PCOUNT': 0, - 'GCOUNT': 1, - 'CCDCHIP': ccdchip, - } + kw_val0 = {'XTENSION': 'IMAGE', + 'BITPIX': -32, + 'NAXIS': naxis, + 'EXTNAME': 'WCSDVARR', + 'EXTVER': wdvarr_ver, + 'PCOUNT': 0, + 'GCOUNT': 1, + 'CCDCHIP': ccdchip, + } cdl = [] for key in kw_comm0.keys(): cdl.append((key, kw_val0[key], kw_comm0[key])) @@ -310,11 +313,11 @@ class NPOLCorr(object): for i, c in enumerate(npol_phdr): if c == 'FILENAME': start_indx = i - if c == '': # remove blanks from end of header - end_indx = i+1 + if c == '': # remove blanks from end of header + end_indx = i + 1 break if start_indx >= 0: - for card in npol_phdr.cards[start_indx:end_indx]: + for card in npol_phdr.cards[start_indx: end_indx]: cdl.append(card) hdr = fits.Header(cards=cdl) @@ -331,7 +334,7 @@ class NPOLCorr(object): if fobj[0].header['INSTRUME'] == 'ACS' and fobj[0].header['DETECTOR'] == 'WFC': ccdchip = fobj[extname, extver].header['CCDCHIP'] elif fobj[0].header['INSTRUME'] == 'WFC3' and fobj[0].header['DETECTOR'] == 'UVIS': - ccdchip = fobj[extname, extver].header['CCDCHIP'] + ccdchip = fobj[extname, extver].header['CCDCHIP'] elif fobj[0].header['INSTRUME'] == 'WFPC2': ccdchip = fobj[extname, extver].header['DETECTOR'] elif fobj[0].header['INSTRUME'] == 'STIS': diff --git a/stwcs/updatewcs/utils.py b/stwcs/updatewcs/utils.py index ebfc701..8c8d416 100644 --- a/stwcs/updatewcs/utils.py +++ b/stwcs/updatewcs/utils.py @@ -1,4 +1,4 @@ -from __future__ import division # confidence high +from __future__ import absolute_import, division, print_function import os from astropy.io import fits from stsci.tools import fileutil @@ -7,7 +7,7 @@ import logging logger = logging.getLogger("stwcs.updatewcs.utils") -def diff_angles(a,b): +def diff_angles(a, b): """ Perform angle subtraction a-b taking into account small-angle differences across 360degree line. @@ -23,6 +23,7 @@ def diff_angles(a,b): return diff + def getBinning(fobj, extver=1): # Return the binning factor binned = 1 @@ -30,9 +31,10 @@ def getBinning(fobj, extver=1): mode = fobj[0].header.get('MODE', "") if mode == 'AREA': binned = 2 else: - binned = fobj['SCI', extver].header.get('BINAXIS',1) + binned = fobj['SCI', extver].header.get('BINAXIS', 1) return binned + def updateNEXTENDKw(fobj): """ Updates PRIMARY header with correct value for NEXTEND, if present. @@ -44,9 +46,10 @@ def updateNEXTENDKw(fobj): """ if 'nextend' in fobj[0].header: - fobj[0].header['nextend'] = len(fobj) -1 + fobj[0].header['nextend'] = len(fobj) - 1 + -def extract_rootname(kwvalue,suffix=""): +def extract_rootname(kwvalue, suffix=""): """ Returns the rootname from a full reference filename If a non-valid value (any of ['','N/A','NONE','INDEF',None]) is input, @@ -55,13 +58,13 @@ def extract_rootname(kwvalue,suffix=""): This function will also replace any 'suffix' specified with a blank. """ # check to see whether a valid kwvalue has been provided as input - if kwvalue.strip() in ['','N/A','NONE','INDEF',None]: - return 'NONE' # no valid value, so return 'NONE' + if kwvalue.strip() in ['', 'N/A', 'NONE', 'INDEF', None]: + return 'NONE' # no valid value, so return 'NONE' # for a valid kwvalue, parse out the rootname # strip off any environment variable from input filename, if any are given if '$' in kwvalue: - fullval = kwvalue[kwvalue.find('$')+1:] + fullval = kwvalue[kwvalue.find('$') + 1: ] else: fullval = kwvalue # Extract filename without path from kwvalue @@ -71,10 +74,11 @@ def extract_rootname(kwvalue,suffix=""): rootname = fileutil.buildNewRootname(fname) # Now, remove any known suffix from rootname - rootname = rootname.replace(suffix,'') + rootname = rootname.replace(suffix, '') return rootname.strip() -def construct_distname(fobj,wcsobj): + +def construct_distname(fobj, wcsobj): """ This function constructs the value for the keyword 'DISTNAME'. It relies on the reference files specified by the keywords 'IDCTAB', @@ -85,8 +89,8 @@ def construct_distname(fobj,wcsobj): and have a value of 'NONE' if no reference files are specified. """ idcname = extract_rootname(fobj[0].header.get('IDCTAB', "NONE"), - suffix='_idc') - if (idcname is None or idcname=='NONE') and wcsobj.sip is not None: + suffix='_idc') + if (idcname is None or idcname == 'NONE') and wcsobj.sip is not None: idcname = 'UNKNOWN' npolname, npolfile = build_npolname(fobj) @@ -99,29 +103,30 @@ def construct_distname(fobj,wcsobj): sipname, idctab = build_sipname(fobj) - distname = build_distname(sipname,npolname,d2imname) - return {'DISTNAME':distname,'SIPNAME':sipname} + distname = build_distname(sipname, npolname, d2imname) + return {'DISTNAME': distname, 'SIPNAME': sipname} + -def build_distname(sipname,npolname,d2imname): +def build_distname(sipname, npolname, d2imname): """ Core function to build DISTNAME keyword value without the HSTWCS input. """ - distname = sipname.strip() if npolname != 'NONE' or d2imname != 'NONE': if d2imname != 'NONE': - distname+= '-'+npolname.strip() + '-'+d2imname.strip() + distname += '-' + npolname.strip() + '-' + d2imname.strip() else: - distname+='-'+npolname.strip() + distname += '-' + npolname.strip() return distname -def build_default_wcsname(idctab): - idcname = extract_rootname(idctab,suffix='_idc') +def build_default_wcsname(idctab): + idcname = extract_rootname(idctab, suffix='_idc') wcsname = 'IDC_' + idcname return wcsname + def build_sipname(fobj, fname=None, sipname=None): """ Build a SIPNAME from IDCTAB @@ -141,10 +146,10 @@ def build_sipname(fobj, fname=None, sipname=None): """ try: idctab = fobj[0].header['IDCTAB'] - idcname = extract_rootname(idctab,suffix='_idc') + idcname = extract_rootname(idctab, suffix='_idc') except KeyError: idctab = 'N/A' - idcname= 'N/A' + idcname = 'N/A' if not fname: try: fname = fobj.filename() @@ -161,7 +166,7 @@ def build_sipname(fobj, fname=None, sipname=None): rootname = fobj[0].header['rootname'] except KeyError: rootname = fname - sipname = rootname +'_'+ idcname + sipname = rootname + '_' + idcname else: if 'A_ORDER' in fobj[1].header or 'B_ORDER' in fobj[1].header: sipname = 'UNKNOWN' @@ -170,6 +175,7 @@ def build_sipname(fobj, fname=None, sipname=None): return sipname, idctab + def build_npolname(fobj, npolfile=None): """ Build a NPOLNAME from NPOLFILE @@ -203,6 +209,7 @@ def build_npolname(fobj, npolfile=None): npolname = 'NOMODEL' return npolname, npolfile + def build_d2imname(fobj, d2imfile=None): """ Build a D2IMNAME from D2IMFILE @@ -227,11 +234,11 @@ def build_d2imname(fobj, d2imfile=None): d2imname = 'UNKNOWN' else: d2imname = 'NOMODEL' - d2imname = extract_rootname(d2imfile,suffix='_d2i') + d2imname = extract_rootname(d2imfile, suffix='_d2i') if d2imname == 'NONE': d2imname = 'NOMODEL' else: - d2imname = extract_rootname(d2imfile,suffix='_d2i') + d2imname = extract_rootname(d2imfile, suffix='_d2i') if d2imname == 'NONE': d2imname = 'NOMODEL' return d2imname, d2imfile diff --git a/stwcs/updatewcs/wfpc2_dgeo.py b/stwcs/updatewcs/wfpc2_dgeo.py index e57bb5c..b76198a 100644 --- a/stwcs/updatewcs/wfpc2_dgeo.py +++ b/stwcs/updatewcs/wfpc2_dgeo.py @@ -1,6 +1,7 @@ """ wfpc2_dgeo - Functions to convert WFPC2 DGEOFILE into D2IMFILE """ +from __future__ import absolute_import, division, print_function import os import datetime @@ -13,6 +14,7 @@ from stsci.tools import fileutil import logging logger = logging.getLogger("stwcs.updatewcs.apply_corrections") + def update_wfpc2_d2geofile(filename, fhdu=None): """ Creates a D2IMFILE from the DGEOFILE for a WFPC2 image (input), and @@ -35,7 +37,7 @@ def update_wfpc2_d2geofile(filename, fhdu=None): image header will be updated/added to point to this newly created file. """ - + close_fhdu = False if fhdu is None: fhdu = fileutil.openImage(filename, mode='update') @@ -48,7 +50,7 @@ def update_wfpc2_d2geofile(filename, fhdu=None): dgeofile = fhdu['PRIMARY'].header.get('ODGEOFIL', None) logger.info('Converting DGEOFILE %s into D2IMFILE...' % dgeofile) rootname = filename[:filename.find('.fits')] - d2imfile = convert_dgeo_to_d2im(dgeofile,rootname) + d2imfile = convert_dgeo_to_d2im(dgeofile, rootname) fhdu['PRIMARY'].header['ODGEOFIL'] = dgeofile fhdu['PRIMARY'].header['DGEOFILE'] = 'N/A' fhdu['PRIMARY'].header['D2IMFILE'] = d2imfile @@ -67,25 +69,26 @@ def update_wfpc2_d2geofile(filename, fhdu=None): # (multidrizzle clean=True mode of operation) return d2imfile -def convert_dgeo_to_d2im(dgeofile,output,clobber=True): + +def convert_dgeo_to_d2im(dgeofile, output, clobber=True): """ Routine that converts the WFPC2 DGEOFILE into a D2IMFILE. """ dgeo = fileutil.openImage(dgeofile) - outname = output+'_d2im.fits' + outname = output + '_d2im.fits' removeFileSafely(outname) - data = np.array([dgeo['dy',1].data[:,0]]) + data = np.array([dgeo['dy', 1].data[:, 0]]) scihdu = fits.ImageHDU(data=data) dgeo.close() # add required keywords for D2IM header scihdu.header['EXTNAME'] = ('DY', 'Extension name') scihdu.header['EXTVER'] = (1, 'Extension version') - fits_str = 'PYFITS Version '+str(astropy.__version__) + fits_str = 'PYFITS Version ' + str(astropy.__version__) scihdu.header['ORIGIN'] = (fits_str, 'FITS file originator') scihdu.header['INHERIT'] = (False, 'Inherits global header') dnow = datetime.datetime.now() - scihdu.header['DATE'] = (str(dnow).replace(' ','T'), + scihdu.header['DATE'] = (str(dnow).replace(' ', 'T'), 'Date FITS file was generated') scihdu.header['CRPIX1'] = (0, 'Distortion array reference pixel') @@ -116,7 +119,7 @@ def convert_dgeo_to_d2im(dgeofile,output,clobber=True): return outname -def removeFileSafely(filename,clobber=True): +def removeFileSafely(filename, clobber=True): """ Delete the file specified, but only if it exists and clobber is True. """ if filename is not None and filename.strip() != '': diff --git a/stwcs/wcsutil/__init__.py b/stwcs/wcsutil/__init__.py index 65280be..d9d21c5 100644 --- a/stwcs/wcsutil/__init__.py +++ b/stwcs/wcsutil/__init__.py @@ -1,4 +1,4 @@ -from __future__ import absolute_import, print_function # confidence high +from __future__ import absolute_import, print_function from .altwcs import * from .hstwcs import HSTWCS diff --git a/stwcs/wcsutil/altwcs.py b/stwcs/wcsutil/altwcs.py index 5dc7dd6..64cdc87 100644 --- a/stwcs/wcsutil/altwcs.py +++ b/stwcs/wcsutil/altwcs.py @@ -1,7 +1,6 @@ -from __future__ import division, print_function # confidence high -import os -import string +from __future__ import absolute_import, division, print_function +import string import numpy as np from astropy import wcs as pywcs from astropy.io import fits @@ -9,9 +8,11 @@ from stsci.tools import fileutil as fu altwcskw = ['WCSAXES', 'CRVAL', 'CRPIX', 'PC', 'CDELT', 'CD', 'CTYPE', 'CUNIT', 'PV', 'PS'] -altwcskw_extra = ['LATPOLE','LONPOLE','RESTWAV','RESTFRQ'] +altwcskw_extra = ['LATPOLE', 'LONPOLE', 'RESTWAV', 'RESTFRQ'] # file operations + + def archiveWCS(fname, ext, wcskey=" ", wcsname=" ", reusekey=False): """ Copy the primary WCS to the header as an alternate WCS @@ -77,7 +78,7 @@ def archiveWCS(fname, ext, wcskey=" ", wcsname=" ", reusekey=False): closefobj(fname, f) raise KeyError("Wcskey %s is aready used. \ Run archiveWCS() with reusekey=True to overwrite this alternate WCS. \ - Alternatively choose another wcskey with altwcs.available_wcskeys()." %wcskey) + Alternatively choose another wcskey with altwcs.available_wcskeys()." % wcskey) elif wcskey == " ": # wcsname exists, overwrite it if reuse is True or get the next key if wcsname.strip() in wcsnames(f[wcsext].header).values(): @@ -89,13 +90,13 @@ def archiveWCS(fname, ext, wcskey=" ", wcsname=" ", reusekey=False): wkey = next_wcskey(f[wcsext].header) elif wkey is None: closefobj(fname, f) - raise KeyError("Could not get a valid wcskey from wcsname %s" %wcsname) + raise KeyError("Could not get a valid wcskey from wcsname %s" % wcsname) else: closefobj(fname, f) raise KeyError("Wcsname %s is aready used. \ Run archiveWCS() with reusekey=True to overwrite this alternate WCS. \ Alternatively choose another wcskey with altwcs.available_wcskeys() or\ - choose another wcsname." %wcsname) + choose another wcsname." % wcsname) else: wkey = next_wcskey(f[wcsext].header) if wcsname.strip(): @@ -103,7 +104,7 @@ def archiveWCS(fname, ext, wcskey=" ", wcsname=" ", reusekey=False): else: # determine which WCSNAME needs to be replicated in archived WCS wnames = wcsnames(f[wcsext].header) - if 'O' in wnames: del wnames['O'] # we don't want OPUS/original + if 'O' in wnames: del wnames['O'] # we don't want OPUS/original if len(wnames) > 0: if ' ' in wnames: wname = wnames[' '] @@ -139,15 +140,16 @@ def archiveWCS(fname, ext, wcskey=" ", wcsname=" ", reusekey=False): f[e].header[wcsnamekey] = wname try: - old_wcsname=hwcs.pop('WCSNAME') + old_wcsname = hwcs.pop('WCSNAME') except: pass for k in hwcs.keys(): - key = k[:7] + wkey + key = k[: 7] + wkey f[e].header[key] = hwcs[k] closefobj(fname, f) + def restore_from_to(f, fromext=None, toext=None, wcskey=" ", wcsname=" "): """ Copy an alternate WCS from one extension as a primary WCS of another extension @@ -188,14 +190,14 @@ def restore_from_to(f, fromext=None, toext=None, wcskey=" ", wcsname=" "): raise ValueError("Input parameters problem") # Interpret input 'ext' value to get list of extensions to process - #ext = _buildExtlist(fobj, ext) + # ext = _buildExtlist(fobj, ext) if isinstance(toext, str): toext = [toext] # the case of an HDUList object in memory without an associated file - #if fobj.filename() is not None: + # if fobj.filename() is not None: # name = fobj.filename() simplefits = fu.isFits(fobj)[1] is 'simple' @@ -221,14 +223,15 @@ def restore_from_to(f, fromext=None, toext=None, wcskey=" ", wcsname=" "): if not countext: raise KeyError("File does not have extension with extname %s", fromext) else: - for i in range(1, countext+1): + for i in range(1, countext + 1): for toe in toext: _restore(fobj, fromextnum=i, fromextnam=fromext, toextnum=i, toextnam=toe, ukey=wkey) if fobj.filename() is not None: - #fobj.writeto(name) + # fobj.writeto(name) closefobj(f, fobj) + def restoreWCS(f, ext, wcskey=" ", wcsname=" "): """ Copy a WCS with key "WCSKEY" to the primary WCS @@ -272,9 +275,6 @@ def restoreWCS(f, ext, wcskey=" ", wcsname=" "): # the case of an HDUList object in memory without an associated file - #if fobj.filename() is not None: - # name = fobj.filename() - simplefits = fu.isFits(fobj)[1] is 'simple' if simplefits: wcskeyext = 0 @@ -297,6 +297,7 @@ def restoreWCS(f, ext, wcskey=" ", wcsname=" "): if fobj.filename() is not None: closefobj(f, fobj) + def deleteWCS(fname, ext, wcskey=" ", wcsname=" "): """ Delete an alternate WCS defined with wcskey. @@ -351,7 +352,7 @@ def deleteWCS(fname, ext, wcskey=" ", wcsname=" "): prexts = [] for i in ext: hdr = fobj[i].header - hwcs = readAltWCS(fobj,i,wcskey=wkey) + hwcs = readAltWCS(fobj, i, wcskey=wkey) if hwcs is None: continue for k in hwcs[::-1]: @@ -363,6 +364,7 @@ def deleteWCS(fname, ext, wcskey=" ", wcsname=" "): print("Did not find WCS with key %s in any of the extensions" % wkey) closefobj(fname, fobj) + def _buildExtlist(fobj, ext): """ Utility function to interpret the provided value of 'ext' and return a list @@ -378,8 +380,8 @@ def _buildExtlist(fobj, ext): If a string is provided, it should specify the EXTNAME of extensions with WCSs to be archived """ - if not isinstance(ext,list): - if isinstance(ext,str): + if not isinstance(ext, list): + if isinstance(ext, str): extstr = ext ext = [] for extn in range(1, len(fobj)): @@ -391,6 +393,7 @@ def _buildExtlist(fobj, ext): raise KeyError("Valid extensions in 'ext' parameter need to be specified.") return ext + def _restore(fobj, ukey, fromextnum, toextnum=None, fromextnam=None, toextnam=None, verbose=True): """ @@ -415,7 +418,7 @@ def _restore(fobj, ukey, fromextnum, if toextnam: toextension = (toextnam, toextnum) else: - toextension =toextnum + toextension = toextnum else: toextension = fromextension @@ -445,34 +448,38 @@ def _restore(fobj, ukey, fromextnum, fobj[toextension].header['TDDALPHA'] = 0.0 fobj[toextension].header['TDDBETA'] = 0.0 if 'ORIENTAT' in fobj[toextension].header: - norient = np.rad2deg(np.arctan2(hwcs['CD1_2'+'%s' %ukey], hwcs['CD2_2'+'%s' %ukey])) + norient = np.rad2deg(np.arctan2(hwcs['CD1_2' + '%s' % ukey], hwcs['CD2_2' + '%s' % ukey])) fobj[toextension].header['ORIENTAT'] = norient # Reset 2014 TDD keywords prior to computing new values (if any are computed) - for kw in ['TDD_CYA','TDD_CYB','TDD_CXA','TDD_CXB']: + for kw in ['TDD_CYA', 'TDD_CYB', 'TDD_CXA', 'TDD_CXB']: if kw in fobj[toextension].header: fobj[toextension].header[kw] = 0.0 -#header operations + +# header operations + + def _check_headerpars(fobj, ext): if not isinstance(fobj, fits.Header) and not isinstance(fobj, fits.HDUList) \ - and not isinstance(fobj, str): + and not isinstance(fobj, str): raise ValueError("Expected a file name, a file object or a header\n") if not isinstance(fobj, fits.Header): - #raise ValueError("Expected a valid ext parameter when input is a file") if not isinstance(ext, int) and not isinstance(ext, tuple): raise ValueError("Expected ext to be a number or a tuple, e.g. ('SCI', 1)\n") + def _getheader(fobj, ext): if isinstance(fobj, str): - hdr = fits.getheader(fobj,ext) + hdr = fits.getheader(fobj, ext) elif isinstance(fobj, fits.Header): hdr = fobj else: hdr = fobj[ext].header return hdr -def readAltWCS(fobj, ext, wcskey=' ',verbose=False): + +def readAltWCS(fobj, ext, wcskey=' ', verbose=False): """ Reads in alternate primary WCS from specified extension. @@ -495,12 +502,12 @@ def readAltWCS(fobj, ext, wcskey=' ',verbose=False): if isinstance(fobj, str): fobj = fits.open(fobj) - hdr = _getheader(fobj,ext) + hdr = _getheader(fobj, ext) try: nwcs = pywcs.WCS(hdr, fobj=fobj, key=wcskey) except KeyError: if verbose: - print('readAltWCS: Could not read WCS with key %s' %wcskey) + print('readAltWCS: Could not read WCS with key %s' % wcskey) print(' Skipping %s[%s]' % (fobj.filename(), str(ext))) return None hwcs = nwcs.to_header() @@ -510,7 +517,8 @@ def readAltWCS(fobj, ext, wcskey=' ',verbose=False): return hwcs -def convertAltWCS(fobj,ext,oldkey=" ",newkey=' '): + +def convertAltWCS(fobj, ext, oldkey=" ", newkey=' '): """ Translates the alternate/primary WCS with one key to an alternate/primary WCS with another key. @@ -534,7 +542,7 @@ def convertAltWCS(fobj,ext,oldkey=" ",newkey=' '): hdr: `astropy.io.fits.Header` header object with keywords renamed from oldkey to newkey """ - hdr = readAltWCS(fobj,ext,wcskey=oldkey) + hdr = readAltWCS(fobj, ext, wcskey=oldkey) if hdr is None: return None # Converting WCS to new key @@ -543,10 +551,11 @@ def convertAltWCS(fobj,ext,oldkey=" ",newkey=' '): cname = card else: cname = card.rstrip(oldkey) - hdr.rename_key(card,cname+newkey,force=True) + hdr.rename_key(card, cname + newkey, force=True) return hdr + def wcskeys(fobj, ext=None): """ Returns a list of characters used in the header for alternate @@ -565,11 +574,12 @@ def wcskeys(fobj, ext=None): names = hdr["WCSNAME*"] d = [] for key in names: - wkey = key.replace('WCSNAME','') + wkey = key.replace('WCSNAME', '') if wkey == '': wkey = ' ' d.append(wkey) return d + def wcsnames(fobj, ext=None): """ Returns a dictionary of wcskey: WCSNAME pairs @@ -588,11 +598,12 @@ def wcsnames(fobj, ext=None): names = hdr["WCSNAME*"] d = {} for keyword, value in names.items(): - wkey = keyword.replace('WCSNAME','') + wkey = keyword.replace('WCSNAME', '') if wkey == '': wkey = ' ' d[wkey] = value return d + def available_wcskeys(fobj, ext=None): """ Returns a list of characters which are not used in the header @@ -618,6 +629,7 @@ def available_wcskeys(fobj, ext=None): [all_keys.remove(key) for key in used_keys] return all_keys + def next_wcskey(fobj, ext=None): """ Returns next available character to be used for an alternate WCS @@ -638,6 +650,7 @@ def next_wcskey(fobj, ext=None): else: return None + def getKeyFromName(header, wcsname): """ If WCSNAME is found in header, return its key, else return @@ -663,6 +676,7 @@ def getKeyFromName(header, wcsname): wkey = None return wkey + def pc2cd(hdr, key=' '): """ Convert a CD matrix to a CD matrix. @@ -690,6 +704,7 @@ def pc2cd(hdr, key=' '): hdr['CD{0}{1}'.format(c, key)] = val return hdr + def _parpasscheck(fobj, ext, wcskey, fromext=None, toext=None, reusekey=False): """ Check input parameters to altwcs functions @@ -722,7 +737,7 @@ def _parpasscheck(fobj, ext, wcskey, fromext=None, toext=None, reusekey=False): return False if not isinstance(ext, int) and not isinstance(ext, tuple) \ - and not isinstance(ext,str) \ + and not isinstance(ext, str) \ and not isinstance(ext, list) and ext is not None: print("Ext must be integer, tuple, string,a list of int extension numbers, \n\ or a list of tuples representing a fits extension, for example ('sci', 1).") @@ -733,7 +748,7 @@ def _parpasscheck(fobj, ext, wcskey, fromext=None, toext=None, reusekey=False): return False if not isinstance(toext, list) and not isinstance(toext, str) and \ - toext is not None : + toext is not None: print("toext must be a string or a list of strings representing extname") return False @@ -743,6 +758,7 @@ def _parpasscheck(fobj, ext, wcskey, fromext=None, toext=None, reusekey=False): return True + def closefobj(fname, f): """ Functions in this module accept as input a file name or a file object. @@ -752,6 +768,7 @@ def closefobj(fname, f): if isinstance(fname, str): f.close() + def mapFitsExt2HDUListInd(fname, extname): """ Map FITS extensions with 'EXTNAME' to HDUList indexes. diff --git a/stwcs/wcsutil/convertwcs.py b/stwcs/wcsutil/convertwcs.py index a384eb1..e47a829 100644 --- a/stwcs/wcsutil/convertwcs.py +++ b/stwcs/wcsutil/convertwcs.py @@ -1,4 +1,5 @@ -from __future__ import print_function +from __future__ import absolute_import, division, print_function +from astropy.io import fits try: import stwcs from stwcs import wcsutil @@ -7,12 +8,12 @@ except: from stsci.tools import fileutil -OPUS_WCSKEYS = ['OCRVAL1','OCRVAL2','OCRPIX1','OCRPIX2', - 'OCD1_1','OCD1_2','OCD2_1','OCD2_2', - 'OCTYPE1','OCTYPE2'] +OPUS_WCSKEYS = ['OCRVAL1', 'OCRVAL2', 'OCRPIX1', 'OCRPIX2', + 'OCD1_1', 'OCD1_2', 'OCD2_1', 'OCD2_2', + 'OCTYPE1', 'OCTYPE2'] -def archive_prefix_OPUS_WCS(fobj,extname='SCI'): +def archive_prefix_OPUS_WCS(fobj, extname='SCI'): """ Identifies WCS keywords which were generated by OPUS and archived using a prefix of 'O' for all 'SCI' extensions in the file @@ -28,43 +29,42 @@ def archive_prefix_OPUS_WCS(fobj,extname='SCI'): print('=====================') raise ImportError - closefits = False - if isinstance(fobj,str): + if isinstance(fobj, str): # A filename was provided as input - fobj = fits.open(fobj,mode='update') - closefits=True + fobj = fits.open(fobj, mode='update') + closefits = True # Define the header - ext = ('sci',1) + ext = ('sci', 1) hdr = fobj[ext].header numextn = fileutil.countExtn(fobj) extlist = [] - for e in range(1,numextn+1): - extlist.append(('sci',e)) + for e in range(1, numextn + 1): + extlist.append(('sci', e)) # Insure that the 'O' alternate WCS is present if 'O' not in wcsutil.wcskeys(hdr): # if not, archive the Primary WCS as the default OPUS WCS - wcsutil.archiveWCS(fobj,extlist, wcskey='O', wcsname='OPUS') + wcsutil.archiveWCS(fobj, extlist, wcskey='O', wcsname='OPUS') # find out how many SCI extensions are in the image - numextn = fileutil.countExtn(fobj,extname=extname) + numextn = fileutil.countExtn(fobj, extname=extname) if numextn == 0: extname = 'PRIMARY' # create HSTWCS object from PRIMARY WCS - wcsobj = wcsutil.HSTWCS(fobj,ext=ext,wcskey='O') + wcsobj = wcsutil.HSTWCS(fobj, ext=ext, wcskey='O') # get list of WCS keywords wcskeys = list(wcsobj.wcs2header().keys()) # For each SCI extension... - for e in range(1,numextn+1): + for e in range(1, numextn + 1): # Now, look for any WCS keywords with a prefix of 'O' for key in wcskeys: - okey = 'O'+key[:7] - hdr = fobj[(extname,e)].header + okey = 'O' + key[: 7] + hdr = fobj[(extname, e)].header if okey in hdr: # Update alternate WCS keyword with prefix-O OPUS keyword value hdr[key] = hdr[okey] @@ -72,7 +72,8 @@ def archive_prefix_OPUS_WCS(fobj,extname='SCI'): if closefits: fobj.close() -def create_prefix_OPUS_WCS(fobj,extname='SCI'): + +def create_prefix_OPUS_WCS(fobj, extname='SCI'): """ Creates alternate WCS with a prefix of 'O' for OPUS generated WCS values to work with old MultiDrizzle. @@ -91,10 +92,10 @@ def create_prefix_OPUS_WCS(fobj,extname='SCI'): owcskeys = OPUS_WCSKEYS closefits = False - if isinstance(fobj,str): + if isinstance(fobj, str): # A filename was provided as input fobj = fits.open(fobj, mode='update') - closefits=True + closefits = True else: # check to make sure this FITS obj has been opened in update mode if fobj.fileinfo(0)['filemode'] != 'update': @@ -102,16 +103,16 @@ def create_prefix_OPUS_WCS(fobj,extname='SCI'): raise IOError # check for existance of O-prefix WCS - if owcskeys[0] not in fobj['sci',1].header: + if owcskeys[0] not in fobj['sci', 1].header: # find out how many SCI extensions are in the image - numextn = fileutil.countExtn(fobj,extname=extname) + numextn = fileutil.countExtn(fobj, extname=extname) if numextn == 0: extname = '' - for extn in range(1,numextn+1): - hdr = fobj[(extname,extn)].header + for extn in range(1, numextn + 1): + hdr = fobj[(extname, extn)].header for okey in owcskeys: - hdr[okey] = hdr[okey[1:]+'O'] + hdr[okey] = hdr[okey[1: ] + 'O'] # Close FITS image if we had to open it... if closefits: diff --git a/stwcs/wcsutil/getinput.py b/stwcs/wcsutil/getinput.py index 8ee1123..c8b3b1b 100644 --- a/stwcs/wcsutil/getinput.py +++ b/stwcs/wcsutil/getinput.py @@ -1,23 +1,27 @@ +from __future__ import absolute_import, division, print_function + from astropy.io import fits from stsci.tools import irafglob, fileutil, parseinput +#from . import HSTWCS + def parseSingleInput(f=None, ext=None): if isinstance(f, str): # create an HSTWCS object from a filename - if ext != None: + if ext is not None: filename = f - if isinstance(ext,tuple): + if isinstance(ext, tuple): if ext[0] == '': - extnum = ext[1] # handle ext=('',1) + extnum = ext[1] # handle ext=('',1) else: extnum = ext else: extnum = int(ext) - elif ext == None: + elif ext is None: filename, ext = fileutil.parseFilename(f) ext = fileutil.parseExtn(ext) if ext[0] == '': - extnum = int(ext[1]) #handle ext=('',extnum) + extnum = int(ext[1]) # handle ext=('',extnum) else: extnum = ext phdu = fits.open(filename) @@ -29,7 +33,7 @@ def parseSingleInput(f=None, ext=None): elif isinstance(f, fits.HDUList): phdu = f - if ext == None: + if ext is None: extnum = 0 else: extnum = ext @@ -54,9 +58,9 @@ def parseMultipleInput(input): filelist, output = parseinput.parseinput(input) except IOError: raise elif isinstance(input, list): - if isinstance(input[0], wcsutil.HSTWCS): - # a list of HSTWCS objects - return input - else: - filelist = input[:] + #if isinstance(input[0], HSTWCS): + ## a list of HSTWCS objects + #return input + #else: + filelist = input[:] return filelist diff --git a/stwcs/wcsutil/headerlet.py b/stwcs/wcsutil/headerlet.py index 5a980aa..03861df 100644 --- a/stwcs/wcsutil/headerlet.py +++ b/stwcs/wcsutil/headerlet.py @@ -34,7 +34,9 @@ from . import wcscorr from .hstwcs import HSTWCS from .mappings import basic_wcs -#### Logging support functions +# Logging support functions + + class FuncNameLoggingFormatter(logging.Formatter): def __init__(self, fmt=None, datefmt=None): if '%(funcName)s' not in fmt: @@ -59,8 +61,8 @@ logger.addHandler(ch) logger.setLevel(logging.DEBUG) FITS_STD_KW = ['XTENSION', 'BITPIX', 'NAXIS', 'PCOUNT', - 'GCOUNT', 'EXTNAME', 'EXTVER', 'ORIGIN', - 'INHERIT', 'DATE', 'IRAF-TLM'] + 'GCOUNT', 'EXTNAME', 'EXTVER', 'ORIGIN', + 'INHERIT', 'DATE', 'IRAF-TLM'] DEFAULT_SUMMARY_COLS = ['HDRNAME', 'WCSNAME', 'DISTNAME', 'AUTHOR', 'DATE', 'SIPNAME', 'NPOLFILE', 'D2IMFILE', 'DESCRIP'] @@ -120,10 +122,13 @@ def with_logging(func): return func(*args, **kw) return wrapped -#### Utility functions +# Utility functions + + def is_par_blank(par): return par in ['', ' ', 'INDEF', "None", None] + def parse_filename(fname, mode='readonly'): """ Interprets the input as either a filename of a file that needs to be opened @@ -174,6 +179,7 @@ def parse_filename(fname, mode='readonly'): fname = '' return fobj, fname, close_fobj + def get_headerlet_kw_names(fobj, kw='HDRNAME'): """ Returns a list of specified keywords from all HeaderletHDU @@ -198,6 +204,7 @@ def get_headerlet_kw_names(fobj, kw='HDRNAME'): return hdrnames + def get_header_kw_vals(hdr, kwname, kwval, default=0): if kwval is None: if kwname in hdr: @@ -206,6 +213,7 @@ def get_header_kw_vals(hdr, kwname, kwval, default=0): kwval = default return kwval + @with_logging def find_headerlet_HDUs(fobj, hdrext=None, hdrname=None, distname=None, strict=True, logging=False, logmode='w'): @@ -259,9 +267,9 @@ def find_headerlet_HDUs(fobj, hdrext=None, hdrname=None, distname=None, hdrlets = [] if hdrext is not None and isinstance(hdrext, int): - if hdrext in range(len(fobj)): # insure specified hdrext is in fobj + if hdrext in range(len(fobj)): # insure specified hdrext is in fobj if isinstance(fobj[hdrext], fits.hdu.base.NonstandardExtHDU) and \ - fobj[hdrext].header['EXTNAME'] == 'HDRLET': + fobj[hdrext].header['EXTNAME'] == 'HDRLET': hdrlets.append(hdrext) else: for ext in fobj: @@ -280,9 +288,9 @@ def find_headerlet_HDUs(fobj, hdrext=None, hdrname=None, distname=None, (hdrextnum == ext.header['EXTVER']) and (hdrextname == ext.header['EXTNAME'])) hdrname_match = ((hdrname is not None) and - (hdrname == ext.header['HDRNAME'])) + (hdrname == ext.header['HDRNAME'])) distname_match = ((distname is not None) and - (distname == ext.header['DISTNAME'])) + (distname == ext.header['DISTNAME'])) if hdrext_match or hdrname_match or distname_match: hdrlets.append(fobj.index(ext)) @@ -310,6 +318,7 @@ def find_headerlet_HDUs(fobj, hdrext=None, hdrname=None, distname=None, return hdrlets + def verify_hdrname_is_unique(fobj, hdrname): """ Verifies that no other HeaderletHDU extension has the specified hdrname. @@ -331,6 +340,7 @@ def verify_hdrname_is_unique(fobj, hdrname): return unique + def update_versions(sourcehdr, desthdr): """ Update keywords which store version numbers @@ -344,6 +354,7 @@ def update_versions(sourcehdr, desthdr): except KeyError: desthdr[key] = (" ", phdukw[key]) + def update_ref_files(source, dest): """ Update the reference files name in the primary header of 'dest' @@ -372,8 +383,9 @@ def update_ref_files(source, dest): phdukw[key] = False return phdukw + def print_summary(summary_cols, summary_dict, pad=2, maxwidth=None, idcol=None, - output=None, clobber=True, quiet=False ): + output=None, clobber=True, quiet=False ): """ Print out summary dictionary to STDOUT, and possibly an output file @@ -404,9 +416,9 @@ def print_summary(summary_cols, summary_dict, pad=2, maxwidth=None, idcol=None, for row in range(nrows): if idcol: outstr += COLUMN_FMT.format(idcol['vals'][row], - width=idcol['width']+pad) + width=idcol['width'] + pad) for kw in summary_cols: - val = summary_dict[kw]['vals'][row][:(column_widths[kw]-pad)] + val = summary_dict[kw]['vals'][row][:(column_widths[kw] - pad)] outstr += COLUMN_FMT.format(val, width=column_widths[kw]) outstr += '\n' if not quiet: @@ -415,7 +427,7 @@ def print_summary(summary_cols, summary_dict, pad=2, maxwidth=None, idcol=None, # If specified, write info to separate text file write_file = False if output: - output = fu.osfn(output) # Expand any environment variables in filename + output = fu.osfn(output) # Expand any environment variables in filename write_file = True if os.path.exists(output): if clobber: @@ -430,11 +442,13 @@ def print_summary(summary_cols, summary_dict, pad=2, maxwidth=None, idcol=None, fout.write(outstr) fout.close() -#### Private utility functions +# Private utility functions + + def _create_primary_HDU(fobj, fname, wcsext, destim, hdrname, wcsname, - sipname, npolfile, d2imfile, - nmatch,catalog, wcskey, - author, descrip, history): + sipname, npolfile, d2imfile, + nmatch, catalog, wcskey, + author, descrip, history): # convert input values into valid FITS kw values if author is None: author = '' @@ -447,7 +461,7 @@ def _create_primary_HDU(fobj, fname, wcsext, destim, hdrname, wcsname, npolname, npolfile = utils.build_npolname(fobj, npolfile) logger.info("Setting npolfile value to %s" % npolname) - d2imname, d2imfile = utils.build_d2imname(fobj,d2imfile) + d2imname, d2imfile = utils.build_d2imname(fobj, d2imfile) logger.info("Setting d2imfile value to %s" % d2imname) distname = utils.build_distname(sipname, npolname, d2imname) @@ -461,23 +475,14 @@ def _create_primary_HDU(fobj, fname, wcsext, destim, hdrname, wcsname, else: history = '' - rms_ra = fobj[wcsext].header.get("CRDER1"+wcskey, 0) - rms_dec = fobj[wcsext].header.get("CRDER2"+wcskey, 0) + rms_ra = fobj[wcsext].header.get("CRDER1" + wcskey, 0) + rms_dec = fobj[wcsext].header.get("CRDER2" + wcskey, 0) if not nmatch: - nmatch = fobj[wcsext].header.get("NMATCH"+wcskey, 0) + nmatch = fobj[wcsext].header.get("NMATCH" + wcskey, 0) if not catalog: - catalog = fobj[wcsext].header.get('CATALOG'+wcskey, "") + catalog = fobj[wcsext].header.get('CATALOG' + wcskey, "") # get the version of STWCS used to create the WCS of the science file. - #try: - #upwcsver = fobj[0].header.cards[fobj[0].header.index('UPWCSVER')] - #except KeyError: - #upwcsver = pyfits.Card("UPWCSVER", " ", - #"Version of STWCS used to update the WCS") - #try: - #pywcsver = fobj[0].header.cards[fobj[0].header.index('PYWCSVER')] - #except KeyError: - #pywcsver = pyfits.Card("PYWCSVER", " ", - #"Version of PYWCS used to update the WCS") + upwcsver = fobj[0].header.get('UPWCSVER', "") pywcsver = fobj[0].header.get('PYWCSVER', "") # build Primary HDU @@ -495,7 +500,7 @@ def _create_primary_HDU(fobj, fname, wcsext, destim, hdrname, wcsname, phdu.header['D2IMFILE'] = (d2imfile, 'origin of detector to image correction') phdu.header['IDCTAB'] = (idctab, - 'origin of Polynomial Distortion') + 'origin of Polynomial Distortion') phdu.header['AUTHOR'] = (author, 'headerlet created by this user') phdu.header['DESCRIP'] = (descrip, 'Short description of headerlet solution') @@ -526,7 +531,9 @@ def _create_primary_HDU(fobj, fname, wcsext, destim, hdrname, wcsname, return phdu -#### Public Interface functions +# Public Interface functions + + @with_logging def extract_headerlet(filename, output, extnum=None, hdrname=None, clobber=False, logging=True): @@ -611,11 +618,11 @@ def extract_headerlet(filename, output, extnum=None, hdrname=None, @with_logging def write_headerlet(filename, hdrname, output=None, sciext='SCI', - wcsname=None, wcskey=None, destim=None, - sipname=None, npolfile=None, d2imfile=None, - author=None, descrip=None, history=None, - nmatch=None, catalog=None, - attach=True, clobber=False, logging=False): + wcsname=None, wcskey=None, destim=None, + sipname=None, npolfile=None, d2imfile=None, + author=None, descrip=None, history=None, + nmatch=None, catalog=None, + attach=True, clobber=False, logging=False): """ Save a WCS as a headerlet FITS file. @@ -733,22 +740,22 @@ def write_headerlet(filename, hdrname, output=None, sciext='SCI', # Interpret sciext input for this file if isinstance(sciext, int): - sciextlist = [sciext] # allow for specification of simple FITS header + sciextlist = [sciext] # allow for specification of simple FITS header elif isinstance(sciext, str): numsciext = countExtn(fobj, sciext) if numsciext > 0: - sciextlist = [tuple((sciext,i)) for i in range(1, numsciext+1)] + sciextlist = [tuple((sciext, i)) for i in range(1, numsciext + 1)] else: sciextlist = [0] elif isinstance(sciext, list): sciextlist = sciext else: - errstr = "Expected sciext to be a list of FITS extensions with science data\n"+\ - " a valid EXTNAME string, or an integer." + errstr = "Expected sciext to be a list of FITS extensions with science data\n" + \ + " a valid EXTNAME string, or an integer." logger.critical(errstr) raise ValueError - wnames = altwcs.wcsnames(fobj,ext=sciextlist[0]) + wnames = altwcs.wcsnames(fobj, ext=sciextlist[0]) # Insure that WCSCORR table has been created with all original # WCS's recorded prior to adding the headerlet WCS @@ -756,7 +763,7 @@ def write_headerlet(filename, hdrname, output=None, sciext='SCI', if wcsname is None: scihdr = fobj[sciextlist[0]].header - wname = scihdr['wcsname'+wcskey] + wname = scihdr['wcsname' + wcskey] else: wname = wcsname if hdrname in [None, ' ', '']: @@ -764,17 +771,17 @@ def write_headerlet(filename, hdrname, output=None, sciext='SCI', logger.critical('Creating the headerlet from image %s' % fname) hdrletobj = create_headerlet(fobj, sciext=sciextlist, - wcsname=wname, wcskey=wcskey, - hdrname=hdrname, - sipname=sipname, npolfile=npolfile, - d2imfile=d2imfile, author=author, - descrip=descrip, history=history, - nmatch=nmatch, catalog=catalog, - logging=False) + wcsname=wname, wcskey=wcskey, + hdrname=hdrname, + sipname=sipname, npolfile=npolfile, + d2imfile=d2imfile, author=author, + descrip=descrip, history=history, + nmatch=nmatch, catalog=catalog, + logging=False) if attach: # Check to see whether or not a HeaderletHDU with - #this hdrname already exists + # this hdrname already exists hdrnames = get_headerlet_kw_names(fobj) if hdrname not in hdrnames: hdrlet_hdu = HeaderletHDU.fromheaderlet(hdrletobj) @@ -810,14 +817,15 @@ def write_headerlet(filename, hdrname, output=None, sciext='SCI', outname = output if not outname.endswith('.fits'): - outname = '{0}_{1}_hlet.fits'.format(frootname,outname) + outname = '{0}_{1}_hlet.fits'.format(frootname, outname) # If user specifies an output filename for headerlet, write it out hdrletobj.tofile(outname, clobber=clobber) - logger.critical( 'Created Headerlet file %s ' % outname) + logger.critical('Created Headerlet file %s ' % outname) del hdrletobj + @with_logging def create_headerlet(filename, sciext='SCI', hdrname=None, destim=None, wcskey=" ", wcsname=None, @@ -856,7 +864,8 @@ def create_headerlet(filename, sciext='SCI', hdrname=None, destim=None, if " ", use the primary (default) if None use wcsname wcsname: string or None - if wcskey is None use wcsname specified here to choose an alternate WCS for the headerlet + if wcskey is None use wcsname specified here to choose an alternate WCS + for the headerlet sipname: string or None (default) Name of unique file where the polynomial distortion coefficients were read from. If None, the behavior is: @@ -935,7 +944,7 @@ def create_headerlet(filename, sciext='SCI', hdrname=None, destim=None, if not wcsname: # User did not specify a value for 'wcsname' if wcsnamekw in fobj[wcsext].header: - #check if there's a WCSNAME for this wcskey in the header + # check if there's a WCSNAME for this wcskey in the header wcsname = fobj[wcsext].header[wcsnamekw] logger.info("Setting wcsname from header[%s] to %s" % (wcsnamekw, wcsname)) else: @@ -973,8 +982,9 @@ def create_headerlet(filename, sciext='SCI', hdrname=None, destim=None, wkeys = altwcs.wcskeys(fobj, ext=wcsext) if wcskey != ' ': if wcskey not in wkeys: - logger.critical('No WCS with wcskey=%s found in extension %s. Skipping...' % (wcskey, str(wcsext))) - raise ValueError("No WCS with wcskey=%s found in extension %s. Skipping...' % (wcskey, str(wcsext))") + mess = "Skipping extension {0} - no WCS with wcskey={1} found.".format(wcsext, wcskey) + logger.critical(mess) + raise ValueError(mess) # get remaining required keywords if destim is None: @@ -1005,13 +1015,11 @@ def create_headerlet(filename, sciext='SCI', hdrname=None, destim=None, logger.critical(message) raise KeyError - - hdul = [] phdu = _create_primary_HDU(fobj, fname, wcsext, destim, hdrname, wcsname, - sipname, npolfile, d2imfile, - nmatch, catalog, wcskey, - author, descrip, history) + sipname, npolfile, d2imfile, + nmatch, catalog, wcskey, + author, descrip, history) hdul.append(phdu) wcsdvarr_extns = [] """ @@ -1100,14 +1108,6 @@ def create_headerlet(filename, sciext='SCI', hdrname=None, destim=None, whdu.update_ext_version(ihdu.header['D2IM2.EXTVER']) hdul.append(whdu) - - #if hwcs.det2im1 or hwcs.det2im2: - #try: - #darr = hdul[('D2IMARR', 1)] - #except KeyError: - #whdu = whdul[('D2IMARR')] - #whdu.update_ext_version(1) - #hdul.append(whdu) if close_file: fobj.close() @@ -1115,9 +1115,10 @@ def create_headerlet(filename, sciext='SCI', hdrname=None, destim=None, hlet.init_attrs() return hlet + @with_logging def apply_headerlet_as_primary(filename, hdrlet, attach=True, archive=True, - force=False, logging=False, logmode='a'): + force=False, logging=False, logmode='a'): """ Apply headerlet 'hdrfile' to a science observation 'destfile' as the primary WCS @@ -1146,18 +1147,19 @@ def apply_headerlet_as_primary(filename, hdrlet, attach=True, archive=True, hdrlet = [hdrlet] if len(hdrlet) != len(filename): logger.critical("Filenames must have matching headerlets. " - "{0:d} filenames and {1:d} headerlets specified".format(len(filename),len(hdrlet))) + "{0:d} filenames and {1:d} headerlets specified".format(len(filename), + len(hdrlet))) - for fname,h in zip(filename,hdrlet): - print("Applying {0} as Primary WCS to {1}".format(h,fname)) + for fname, h in zip(filename, hdrlet): + print("Applying {0} as Primary WCS to {1}".format(h, fname)) hlet = Headerlet.fromfile(h, logging=logging, logmode=logmode) hlet.apply_as_primary(fname, attach=attach, archive=archive, - force=force) + force=force) @with_logging def apply_headerlet_as_alternate(filename, hdrlet, attach=True, wcskey=None, - wcsname=None, logging=False, logmode='w'): + wcsname=None, logging=False, logmode='w'): """ Apply headerlet to a science observation as an alternate WCS @@ -1188,13 +1190,14 @@ def apply_headerlet_as_alternate(filename, hdrlet, attach=True, wcskey=None, hdrlet = [hdrlet] if len(hdrlet) != len(filename): logger.critical("Filenames must have matching headerlets. " - "{0:d} filenames and {1:d} headerlets specified".format(len(filename),len(hdrlet))) + "{0:d} filenames and {1:d} headerlets specified".format(len(filename), + len(hdrlet))) - for fname,h in zip(filename,hdrlet): - print('Applying {0} as an alternate WCS to {1}'.format(h,fname)) + for fname, h in zip(filename, hdrlet): + print('Applying {0} as an alternate WCS to {1}'.format(h, fname)) hlet = Headerlet.fromfile(h, logging=logging, logmode=logmode) hlet.apply_as_alternate(fname, attach=attach, - wcsname=wcsname, wcskey=wcskey) + wcsname=wcsname, wcskey=wcskey) @with_logging @@ -1218,12 +1221,13 @@ def attach_headerlet(filename, hdrlet, logging=False, logmode='a'): hdrlet = [hdrlet] if len(hdrlet) != len(filename): logger.critical("Filenames must have matching headerlets. " - "{0:d} filenames and {1:d} headerlets specified".format(len(filename),len(hdrlet))) + "{0:d} filenames and {1:d} headerlets specified".format(len(filename), + len(hdrlet))) - for fname,h in zip(filename,hdrlet): - print('Attaching {0} as Headerlet extension to {1}'.format(h,fname)) + for fname, h in zip(filename, hdrlet): + print('Attaching {0} as Headerlet extension to {1}'.format(h, fname)) hlet = Headerlet.fromfile(h, logging=logging, logmode=logmode) - hlet.attach_to_file(fname,archive=True) + hlet.attach_to_file(fname, archive=True) @with_logging @@ -1262,12 +1266,13 @@ def delete_headerlet(filename, hdrname=None, hdrext=None, distname=None, filename = [filename] for f in filename: - print("Deleting Headerlet from ",f) + print("Deleting Headerlet from ", f) _delete_single_headerlet(f, hdrname=hdrname, hdrext=hdrext, - distname=distname, logging=logging, logmode='a') + distname=distname, logging=logging, logmode='a') + def _delete_single_headerlet(filename, hdrname=None, hdrext=None, distname=None, - logging=False, logmode='w'): + logging=False, logmode='w'): """ Deletes HeaderletHDU(s) from a SINGLE science file @@ -1297,7 +1302,7 @@ def _delete_single_headerlet(filename, hdrname=None, hdrext=None, distname=None, logmode: 'a' or 'w' """ hdrlet_ind = find_headerlet_HDUs(filename, hdrname=hdrname, hdrext=hdrext, - distname=distname, logging=logging, logmode='a') + distname=distname, logging=logging, logmode='a') if len(hdrlet_ind) == 0: message = """ No HDUs deleted... No Headerlet HDUs found with ' @@ -1404,8 +1409,8 @@ def headerlet_summary(filename, columns=None, pad=2, maxwidth=None, # Print out the summary dictionary print_summary(summary_cols, summary_dict, pad=pad, maxwidth=maxwidth, - idcol=extnums_col, output=output, - clobber=clobber, quiet=quiet) + idcol=extnums_col, output=output, + clobber=clobber, quiet=quiet) @with_logging @@ -1452,7 +1457,7 @@ def restore_from_headerlet(filename, hdrname=None, hdrext=None, archive=True, message = """ Multiple Headerlet extensions found with the same name. %d Headerlets with "%s" = %s found in %s. - """% (len(hdrlet_ind), kwerr, kwval, fname) + """ % (len(hdrlet_ind), kwerr, kwval, fname) if close_fobj: fobj.close() logger.critical(message) @@ -1464,7 +1469,7 @@ def restore_from_headerlet(filename, hdrname=None, hdrext=None, archive=True, if hasattr(fobj[hdrlet_ind[0]], 'hdulist'): hdrlet = fobj[hdrlet_indx].hdulist else: - hdrlet = fobj[hdrlet_indx].headerlet # older convention in PyFITS + hdrlet = fobj[hdrlet_indx].headerlet # older convention in PyFITS # read in the names of the extensions which HeaderletHDU updates extlist = [] @@ -1503,7 +1508,7 @@ def restore_from_headerlet(filename, hdrname=None, hdrext=None, archive=True, else: if 'idctab' in scihdr: priwcs_hdrname = ''.join(['IDC_', - utils.extract_rootname(scihdr['idctab'], suffix='_idc')]) + utils.extract_rootname(scihdr['idctab'], suffix='_idc')]) else: priwcs_hdrname = 'UNKNOWN' priwcs_name = priwcs_hdrname @@ -1513,7 +1518,7 @@ def restore_from_headerlet(filename, hdrname=None, hdrext=None, archive=True, if archive and priwcs_unique: if priwcs_unique: newhdrlet = create_headerlet(fobj, sciext=scihdr['extname'], - hdrname=priwcs_hdrname) + hdrname=priwcs_hdrname) newhdrlet.attach_to_file(fobj) # # copy hdrlet as a primary @@ -1598,7 +1603,7 @@ def restore_all_with_distname(filename, distname, primary, archive=True, if hasattr(fobj[primary_ind], 'hdulist'): primary_hdrlet = fobj[primary_ind].hdulist else: - primary_hdrlet = fobj[primary_ind].headerlet # older convention in PyFITS + primary_hdrlet = fobj[primary_ind].headerlet # older convention in PyFITS pri_distname = primary_hdrlet[0].header['distname'] if pri_distname != distname: if close_fobj: @@ -1625,7 +1630,7 @@ def restore_all_with_distname(filename, distname, primary, archive=True, if hasattr(fobj[hlet], 'hdulist'): hdrlet = fobj[hlet].hdulist else: - hdrlet = fobj[hlet].headerlet # older convention in PyFITS + hdrlet = fobj[hlet].headerlet # older convention in PyFITS if hlet == primary_ind: hdrlet.apply_as_primary(fobj, attach=False, archive=archive, force=True) @@ -1641,11 +1646,11 @@ def restore_all_with_distname(filename, distname, primary, archive=True, @with_logging def archive_as_headerlet(filename, hdrname, sciext='SCI', - wcsname=None, wcskey=None, destim=None, - sipname=None, npolfile=None, d2imfile=None, - author=None, descrip=None, history=None, - nmatch=None, catalog=None, - logging=False, logmode='w'): + wcsname=None, wcskey=None, destim=None, + sipname=None, npolfile=None, d2imfile=None, + author=None, descrip=None, history=None, + nmatch=None, catalog=None, + logging=False, logmode='w'): """ Save a WCS as a headerlet extension and write it out to a file. @@ -1735,7 +1740,7 @@ def archive_as_headerlet(filename, hdrname, sciext='SCI', if wcsname is None: scihdr = fobj[sciext, 1].header - wcsname = scihdr['wcsname'+wcskey] + wcsname = scihdr['wcsname' + wcskey] if hdrname in [None, ' ', '']: hdrname = wcsname @@ -1745,13 +1750,13 @@ def archive_as_headerlet(filename, hdrname, sciext='SCI', hdrnames = get_headerlet_kw_names(fobj) if hdrname not in hdrnames: hdrletobj = create_headerlet(fobj, sciext=sciext, - wcsname=wcsname, wcskey=wcskey, - hdrname=hdrname, - sipname=sipname, npolfile=npolfile, - d2imfile=d2imfile, author=author, - descrip=descrip, history=history, - nmatch=nmatch, catalog=catalog, - logging=False) + wcsname=wcsname, wcskey=wcskey, + hdrname=hdrname, + sipname=sipname, npolfile=npolfile, + d2imfile=d2imfile, author=author, + descrip=descrip, history=history, + nmatch=nmatch, catalog=catalog, + logging=False) hlt_hdu = HeaderletHDU.fromheaderlet(hdrletobj) if destim is not None: @@ -1771,7 +1776,10 @@ def archive_as_headerlet(filename, hdrname, sciext='SCI', if close_fobj: fobj.close() -#### Headerlet Class definitions + +# Headerlet Class definitions + + class Headerlet(fits.HDUList): """ A Headerlet class @@ -1811,9 +1819,9 @@ class Headerlet(fits.HDUList): self.distname = self[0].header["DISTNAME"] try: - self.vafactor = self[("SIPWCS", 1)].header.get("VAFACTOR", 1) #None instead of 1? + self.vafactor = self[("SIPWCS", 1)].header.get("VAFACTOR", 1) # None instead of 1? except (IndexError, KeyError): - self.vafactor = self[0].header.get("VAFACTOR", 1) #None instead of 1? + self.vafactor = self[0].header.get("VAFACTOR", 1) # None instead of 1? self.author = self[0].header["AUTHOR"] self.descrip = self[0].header["DESCRIP"] @@ -1846,7 +1854,7 @@ class Headerlet(fits.HDUList): init_logging('class Headerlet', level=logging, mode=logmode) return hlet - def apply_as_primary(self, fobj, attach=True, archive=True, force=False): + def apply_as_primary(self, fobj, attach=True, archive=True, force=False): """ Copy this headerlet as a primary WCS to fobj @@ -1877,7 +1885,8 @@ class Headerlet(fits.HDUList): if close_dest: fobj.close() raise ValueError("Destination name does not match headerlet" - "Observation %s cannot be updated with headerlet %s" % (fname, self.hdrname)) + "Observation {0} cannot be updated with" + "headerlet {1}".format((fname, self.hdrname))) # Check to see whether the distortion model in the destination # matches the distortion model in the headerlet being applied @@ -1886,7 +1895,7 @@ class Headerlet(fits.HDUList): dist_models_equal = self.equal_distmodel(dname) if not dist_models_equal and not force: raise ValueError("Distortion models do not match" - " To overwrite the distortion model, set force=True") + " To overwrite the distortion model, set force=True") orig_hlt_hdu = None numhlt = countExtn(fobj, 'HDRLET') @@ -1896,15 +1905,14 @@ class Headerlet(fits.HDUList): # WCS's recorded prior to adding the headerlet WCS wcscorr.init_wcscorr(fobj) - - ### start archive + # start archive # If archive has been specified - # regardless of whether or not the distortion models are equal... + # regardless of whether or not the distortion models are equal... numsip = countExtn(self, 'SIPWCS') sciext_list = [] alt_hlethdu = [] - for i in range(1, numsip+1): + for i in range(1, numsip + 1): sipheader = self[('SIPWCS', i)].header sciext_list.append((sipheader['TG_ENAME'], sipheader['TG_EVER'])) target_ext = sciext_list[0] @@ -1920,10 +1928,10 @@ class Headerlet(fits.HDUList): # Create a headerlet for the original Primary WCS data in the file, # create an HDU from the original headerlet, and append it to # the file - orig_hlt = create_headerlet(fobj, sciext=sciext_list, #[target_ext], - wcsname=wcsname, - hdrname=hdrname, - logging=self.logging) + orig_hlt = create_headerlet(fobj, sciext=sciext_list, # [target_ext], + wcsname=wcsname, + hdrname=hdrname, + logging=self.logging) orig_hlt_hdu = HeaderletHDU.fromheaderlet(orig_hlt) numhlt += 1 orig_hlt_hdu.header['EXTVER'] = numhlt @@ -1931,7 +1939,6 @@ class Headerlet(fits.HDUList): else: logger.info("Headerlet with name %s is already attached" % hdrname) - if dist_models_equal: # Use the WCSNAME to determine whether or not to archive # Primary WCS as altwcs @@ -1945,8 +1952,8 @@ class Headerlet(fits.HDUList): else: if 'idctab' in scihdr: priwcs_name = ''.join(['IDC_', - utils.extract_rootname(scihdr['idctab'], - suffix='_idc')]) + utils.extract_rootname(scihdr['idctab'], + suffix='_idc')]) else: priwcs_name = 'UNKNOWN' nextkey = altwcs.next_wcskey(fobj, ext=target_ext) @@ -1958,11 +1965,11 @@ class Headerlet(fits.HDUList): if hname != 'OPUS' and hname not in hdrlet_extnames: # get HeaderletHDU for alternate WCS as well alt_hlet = create_headerlet(fobj, sciext=sciext_list, - wcsname=hname, wcskey=wcskey, - hdrname=hname, sipname=None, - npolfile=None, d2imfile=None, - author=None, descrip=None, history=None, - logging=self.logging) + wcsname=hname, wcskey=wcskey, + hdrname=hname, sipname=None, + npolfile=None, d2imfile=None, + author=None, descrip=None, history=None, + logging=self.logging) numhlt += 1 alt_hlet_hdu = HeaderletHDU.fromheaderlet(alt_hlet) alt_hlet_hdu.header['EXTVER'] = numhlt @@ -1970,8 +1977,8 @@ class Headerlet(fits.HDUList): hdrlet_extnames.append(hname) self._del_dest_WCS_ext(fobj) - for i in range(1, numsip+1): - target_ext = sciext_list[i-1] + for i in range(1, numsip + 1): + target_ext = sciext_list[i - 1] self._del_dest_WCS(fobj, target_ext) sipwcs = HSTWCS(self, ('SIPWCS', i)) idckw = sipwcs._idc2hdr() @@ -2035,19 +2042,19 @@ class Headerlet(fits.HDUList): fobj[target_ext].header.extend(priwcs[0].header) if sipwcs.cpdis1: - whdu = priwcs[('WCSDVARR', (i-1)*numnpol+1)].copy() + whdu = priwcs[('WCSDVARR', (i - 1) * numnpol + 1)].copy() whdu.update_ext_version(self[('SIPWCS', i)].header['DP1.EXTVER']) fobj.append(whdu) if sipwcs.cpdis2: - whdu = priwcs[('WCSDVARR', i*numnpol)].copy() + whdu = priwcs[('WCSDVARR', i * numnpol)].copy() whdu.update_ext_version(self[('SIPWCS', i)].header['DP2.EXTVER']) fobj.append(whdu) - if sipwcs.det2im1: #or sipwcs.det2im2: - whdu = priwcs[('D2IMARR', (i-1)*numd2im+1)].copy() + if sipwcs.det2im1: # or sipwcs.det2im2: + whdu = priwcs[('D2IMARR', (i - 1) * numd2im + 1)].copy() whdu.update_ext_version(self[('SIPWCS', i)].header['D2IM1.EXTVER']) fobj.append(whdu) if sipwcs.det2im2: - whdu = priwcs[('D2IMARR', i*numd2im)].copy() + whdu = priwcs[('D2IMARR', i * numd2im)].copy() whdu.update_ext_version(self[('SIPWCS', i)].header['D2IM2.EXTVER']) fobj.append(whdu) @@ -2070,7 +2077,6 @@ class Headerlet(fits.HDUList): if close_dest: fobj.close() - def apply_as_alternate(self, fobj, attach=True, wcskey=None, wcsname=None): """ Copy this headerlet as an alternate WCS to fobj @@ -2098,17 +2104,19 @@ class Headerlet(fits.HDUList): if close_dest: fobj.close() raise ValueError("Destination name does not match headerlet" - "Observation %s cannot be updated with headerlet %s" % (fname, self.hdrname)) + "Observation %s cannot be updated with" + "headerlet %s" % (fname, self.hdrname)) # Verify whether this headerlet has the same distortion - #found in the image being updated + # found in the image being updated dname = self.get_destination_model(fobj) dist_models_equal = self.equal_distmodel(dname) if not dist_models_equal: raise ValueError("Distortion models do not match \n" "Headerlet: %s \n" "Destination file: %s\n" - "attach_to_file() can be used to append this headerlet" %(self.distname, dname)) + "attach_to_file() can be used to append this headerlet" % + (self.distname, dname)) # Insure that WCSCORR table has been created with all original # WCS's recorded prior to adding the headerlet WCS @@ -2148,12 +2156,12 @@ class Headerlet(fits.HDUList): hwcs_header = altwcs.pc2cd(hwcs_header, key=wkey) for ax in range(1, hwcs.naxis + 1): hwcs_header['CTYPE{0}{1}'.format(ax, wkey)] = \ - self[('SIPWCS', 1)].header['CTYPE{0}'.format(ax)] + self[('SIPWCS', 1)].header['CTYPE{0}'.format(ax)] fhdr.extend(hwcs_header) fhdr['WCSNAME' + wkey] = wname # also update with HDRNAME (a non-WCS-standard kw) for kw in self.fit_kws: - #fhdr.insert(wind, pyfits.Card(kw + wkey, + # fhdr.insert(wind, pyfits.Card(kw + wkey, # self[0].header[kw])) fhdr.append(fits.Card(kw + wkey, self[0].header[kw])) # Update the WCSCORR table with new rows from the headerlet's WCSs @@ -2210,7 +2218,7 @@ class Headerlet(fits.HDUList): fobj.close() def info(self, columns=None, pad=2, maxwidth=None, - output=None, clobber=True, quiet=False): + output=None, clobber=True, quiet=False): """ Prints a summary of this headerlet The summary includes: @@ -2240,7 +2248,7 @@ class Headerlet(fits.HDUList): """ summary_cols, summary_dict = self.summary(columns=columns) print_summary(summary_cols, summary_dict, pad=pad, maxwidth=maxwidth, - idcol=None, output=output, clobber=clobber, quiet=quiet) + idcol=None, output=output, clobber=clobber, quiet=quiet) def summary(self, columns=None): """ @@ -2301,7 +2309,7 @@ class Headerlet(fits.HDUList): HDRNAME. """ unique = verify_hdrname_is_unique(dest, self.hdrname) - logger.debug("verify_hdrname() returned %s"%unique) + logger.debug("verify_hdrname() returned %s" % unique) return unique def get_destination_model(self, dest): @@ -2318,7 +2326,7 @@ class Headerlet(fits.HDUList): else: destim = dest dname = destim[0].header['DISTNAME'] if 'distname' in destim[0].header \ - else self.build_distname(dest) + else self.build_distname(dest) if destim_opened: destim.close() return dname @@ -2357,7 +2365,8 @@ class Headerlet(fits.HDUList): else: logger.debug("verify_destim() returned False") logger.critical("Destination name does not match headerlet. " - "Observation %s cannot be updated with headerlet %s" % (fname, self.hdrname)) + "Observation %s cannot be updated with" + "headerlet %s" % (fname, self.hdrname)) return False def build_distname(self, dest): @@ -2377,7 +2386,7 @@ class Headerlet(fits.HDUList): sipname, idctab = utils.build_sipname(dest, dest, None) npolname, npolfile = utils.build_npolname(dest, npolfile) d2imname, d2imfile = utils.build_d2imname(dest, d2imfile) - dname = utils.build_distname(sipname,npolname,d2imname) + dname = utils.build_distname(sipname, npolname, d2imname) return dname def tofile(self, fname, destim=None, hdrname=None, clobber=False): @@ -2418,8 +2427,7 @@ class Headerlet(fits.HDUList): else: for idx in range(numext): # Only delete WCS from extensions which may have WCS keywords - if ('XTENSION' in dest[idx].header and - dest[idx].header['XTENSION'] == 'IMAGE'): + if ('XTENSION' in dest[idx].header and dest[idx].header['XTENSION'] == 'IMAGE'): self._remove_d2im(dest[idx]) self._remove_sip(dest[idx]) self._remove_lut(dest[idx]) @@ -2427,12 +2435,7 @@ class Headerlet(fits.HDUList): self._remove_idc_coeffs(dest[idx]) self._remove_fit_values(dest[idx]) self._remove_ref_files(dest[0]) - """ - if not ext: - self._remove_alt_WCS(dest, ext=range(numext)) - else: - self._remove_alt_WCS(dest, ext=ext) - """ + def _del_dest_WCS_ext(self, dest): numwdvarr = countExtn(dest, 'WCSDVARR') numd2im = countExtn(dest, 'D2IMARR') @@ -2459,13 +2462,13 @@ class Headerlet(fits.HDUList): Remove the any existing astrometric fit values from a FITS extension """ - logger.debug("Removing astrometric fit values from (%s, %s)"% + logger.debug("Removing astrometric fit values from (%s, %s)" % (ext.name, ext.ver)) dkeys = altwcs.wcskeys(ext.header) - if 'O' in dkeys: dkeys.remove('O') # Do not remove wcskey='O' values + if 'O' in dkeys: dkeys.remove('O') # Do not remove wcskey='O' values for fitkw in ['NMATCH', 'CATALOG']: for k in dkeys: - fkw = (fitkw+k).rstrip() + fkw = (fitkw + k).rstrip() if fkw in ext.header: del ext.header[fkw] @@ -2544,7 +2547,7 @@ class Headerlet(fits.HDUList): dkeys = altwcs.wcskeys(dest[('SCI', 1)].header) for val in ['O', '', ' ']: if val in dkeys: - dkeys.remove(val) # Never delete WCS with wcskey='O' + dkeys.remove(val) # Never delete WCS with wcskey='O' logger.debug("Removing alternate WCSs with keys %s from %s" % (dkeys, dest.filename())) @@ -2588,6 +2591,7 @@ class Headerlet(fits.HDUList): except KeyError: pass + @with_logging def _idc2hdr(fromhdr, tohdr, towkey=' '): """ @@ -2598,7 +2602,7 @@ def _idc2hdr(fromhdr, tohdr, towkey=' '): coeffs = ['OCX10', 'OCX11', 'OCY10', 'OCY11', 'IDCSCALE'] for c in coeffs: try: - tohdr[c+towkey] = fromhdr[c] + tohdr[c + towkey] = fromhdr[c] logger.debug("Copied %s to header") except KeyError: continue @@ -2663,8 +2667,8 @@ def get_extname_extver_list(fobj, sciext): else: extlist = sciext[:] else: - errstr = "Expected sciext to be a list of FITS extensions with science data\n"+\ - " a valid EXTNAME string, or an integer." + errstr = "Expected sciext to be a list of FITS extensions with science data\n" + \ + " a valid EXTNAME string, or an integer." logger.critical(errstr) raise ValueError return extlist diff --git a/stwcs/wcsutil/hstwcs.py b/stwcs/wcsutil/hstwcs.py index bfebcfc..27f6467 100644 --- a/stwcs/wcsutil/hstwcs.py +++ b/stwcs/wcsutil/hstwcs.py @@ -1,25 +1,21 @@ -from __future__ import absolute_import, division, print_function # confidence high +from __future__ import absolute_import, division, print_function import os from astropy.wcs import WCS from astropy.io import fits -from stwcs.distortion import models, coeff_converter +from ..distortion import models, coeff_converter import numpy as np from stsci.tools import fileutil -from . import altwcs +from . import pc2cd from . import getinput -from . import mappings from . import instruments from .mappings import inst_mappings, ins_spec_kw -from .mappings import basic_wcs -__docformat__ = 'restructuredtext' +__all__ = ['HSTWCS'] -# -#### Utility functions copied from 'updatewcs.utils' to avoid circular imports -# -def extract_rootname(kwvalue,suffix=""): + +def extract_rootname(kwvalue, suffix=""): """ Returns the rootname from a full reference filename If a non-valid value (any of ['','N/A','NONE','INDEF',None]) is input, @@ -28,13 +24,13 @@ def extract_rootname(kwvalue,suffix=""): This function will also replace any 'suffix' specified with a blank. """ # check to see whether a valid kwvalue has been provided as input - if kwvalue.strip() in ['','N/A','NONE','INDEF',None]: - return 'NONE' # no valid value, so return 'NONE' + if kwvalue.strip() in ['', 'N/A', 'NONE', 'INDEF', None]: + return 'NONE' # no valid value, so return 'NONE' # for a valid kwvalue, parse out the rootname # strip off any environment variable from input filename, if any are given if '$' in kwvalue: - fullval = kwvalue[kwvalue.find('$')+1:] + fullval = kwvalue[kwvalue.find('$') + 1:] else: fullval = kwvalue # Extract filename without path from kwvalue @@ -44,12 +40,13 @@ def extract_rootname(kwvalue,suffix=""): rootname = fileutil.buildNewRootname(fname) # Now, remove any known suffix from rootname - rootname = rootname.replace(suffix,'') + rootname = rootname.replace(suffix, '') return rootname.strip() + def build_default_wcsname(idctab): - idcname = extract_rootname(idctab,suffix='_idc') + idcname = extract_rootname(idctab, suffix='_idc') wcsname = 'IDC_' + idcname return wcsname @@ -93,12 +90,9 @@ class NoConvergence(Exception): self.accuracy = kwargs.pop('accuracy', None) self.niter = kwargs.pop('niter', None) self.divergent = kwargs.pop('divergent', None) - self.failed2converge= kwargs.pop('failed2converge', None) + self.failed2converge = kwargs.pop('failed2converge', None) -# -#### HSTWCS Class definition -# class HSTWCS(WCS): def __init__(self, fobj=None, ext=None, minerr=0.0, wcskey=" "): @@ -138,7 +132,6 @@ class HSTWCS(WCS): self.filename = filename instrument_name = hdr0.get('INSTRUME', 'DEFAULT') if instrument_name == 'DEFAULT' or instrument_name not in list(inst_mappings.keys()): - #['IRAF/ARTDATA','',' ','N/A']: self.instrument = 'DEFAULT' else: self.instrument = instrument_name @@ -189,7 +182,7 @@ class HSTWCS(WCS): Reads in first order IDCTAB coefficients if present in the header """ coeffs = ['ocx10', 'ocx11', 'ocy10', 'ocy11', 'idcscale', - 'idcv2ref','idcv3ref', 'idctheta'] + 'idcv2ref', 'idcv3ref', 'idctheta'] for c in coeffs: self.__setattr__(c, header.get(c, None)) @@ -225,7 +218,7 @@ class HSTWCS(WCS): pass else: - raise KeyError("Unsupported instrument - %s" %self.instrument) + raise KeyError("Unsupported instrument - %s" % self.instrument) def setPscale(self): """ @@ -234,7 +227,7 @@ class HSTWCS(WCS): try: cd11 = self.wcs.cd[0][0] cd21 = self.wcs.cd[1][0] - self.pscale = np.sqrt(np.power(cd11,2)+np.power(cd21,2)) * 3600. + self.pscale = np.sqrt(np.power(cd11, 2) + np.power(cd21, 2)) * 3600. except AttributeError: if self.wcs.has_cd(): print("This file has a PC matrix. You may want to convert it \n \ @@ -249,7 +242,7 @@ class HSTWCS(WCS): try: cd12 = self.wcs.cd[0][1] cd22 = self.wcs.cd[1][1] - self.orientat = np.rad2deg(np.arctan2(cd12,cd22)) + self.orientat = np.rad2deg(np.arctan2(cd12, cd22)) except AttributeError: if self.wcs.has_cd(): print("This file has a PC matrix. You may want to convert it \n \ @@ -261,7 +254,7 @@ class HSTWCS(WCS): """ Updates the CD matrix with a new plate scale """ - self.wcs.cd = self.wcs.cd/self.pscale*scale + self.wcs.cd = self.wcs.cd / self.pscale * scale self.setPscale() def readModel(self, update=False, header=None): @@ -282,8 +275,8 @@ class HSTWCS(WCS): CX10, CX11, CY10, CY11, IDCSCALE, IDCTHETA, IDCXREF, IDCYREF, IDCV2REF, IDCV3REF """ - if self.idctab in [None, '', ' ','N/A']: - #Keyword idctab is not present in header - check for sip coefficients + if self.idctab in [None, '', ' ', 'N/A']: + # Keyword idctab is not present in header - check for sip coefficients if header is not None and 'IDCSCALE' in header: self._readModelFromHeader(header) else: @@ -319,7 +312,6 @@ class HSTWCS(WCS): self.idcmodel = model - def readModelFromIDCTAB(self, header=None, update=False): """ Read distortion model from idc table. @@ -334,25 +326,26 @@ class HSTWCS(WCS): IDCV2REF, IDCV3REF """ - if self.date_obs == None: + if self.date_obs is None: print('date_obs not available\n') self.idcmodel = None return - if self.filter1 == None and self.filter2 == None: + if self.filter1 is None and self.filter2 is None: 'No filter information available\n' self.idcmodel = None return 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) + chip=self.chip, direction='forward', + date=self.date_obs, + filter1=self.filter1, filter2=self.filter2, + offtab=self.offtab, binned=self.binned) if self.ltv1 != 0. or self.ltv2 != 0.: self.resetLTV() if update: - if header==None: + if header is None: print('Update header with IDC model kw requested but header was not provided\n.') else: self._updatehdr(header) @@ -393,7 +386,7 @@ class HSTWCS(WCS): if not wcskey: wcskey = self.wcs.alt if self.wcs.has_cd(): - h = altwcs.pc2cd(h, key=wcskey) + h = pc2cd(h, key=wcskey) if 'wcsname' not in h: if self.idctab is not None: @@ -440,13 +433,13 @@ class HSTWCS(WCS): k - one of 'a', 'b', 'ap', 'bp' """ - cards = [] #fits.CardList() - korder = self.sip.__getattribute__(k+'_order') - cards.append(fits.Card(keyword=k.upper()+'_ORDER', value=korder)) + cards = [] + korder = self.sip.__getattribute__(k + '_order') + cards.append(fits.Card(keyword=k.upper() + '_ORDER', value=korder)) coeffs = self.sip.__getattribute__(k) ind = coeffs.nonzero() for i in range(len(ind[0])): - card = fits.Card(keyword=k.upper()+'_'+str(ind[0][i])+'_'+str(ind[1][i]), + card = fits.Card(keyword=k.upper() + '_' + str(ind[0][i]) + '_' + str(ind[1][i]), value=coeffs[ind[0][i], ind[1][i]]) cards.append(card) return cards @@ -454,7 +447,7 @@ class HSTWCS(WCS): def _idc2hdr(self): # save some of the idc coefficients coeffs = ['ocx10', 'ocx11', 'ocy10', 'ocy11', 'idcscale'] - cards = [] #fits.CardList() + cards = [] for c in coeffs: try: val = self.__getattribute__(c) @@ -711,27 +704,25 @@ adaptive=False, detect_divergence=False, quiet=False) try: ra = np.asarray(args[0], dtype=np.float64) dec = np.asarray(args[1], dtype=np.float64) - #assert( len(ra.shape) == 1 and len(dec.shape) == 1 ) + # assert( len(ra.shape) == 1 and len(dec.shape) == 1 ) origin = int(args[2]) vect1D = True except: - raise TypeError("When providing three arguments, they must " \ - "be (Ra, Dec, origin) where Ra and Dec are " \ - "Nx1 vectors.") + raise TypeError("When providing three arguments, they must " + "be (Ra, Dec, origin) where Ra and Dec are " + "Nx1 vectors.") elif nargs == 2: try: rd = np.asarray(args[0], dtype=np.float64) - #assert( rd.shape[1] == 2 ) - ra = rd[:,0] - dec = rd[:,1] + ra = rd[:, 0] + dec = rd[:, 1] origin = int(args[1]) vect1D = False except: - raise TypeError("When providing two arguments, they must " \ - "be (RaDec, origin) where RaDec is a Nx2 array.") + raise TypeError("When providing two arguments, they must " + "be (RaDec, origin) where RaDec is a Nx2 array.") else: - raise TypeError("Expected 2 or 3 arguments, {:d} given." \ - .format(nargs)) + raise TypeError("Expected 2 or 3 arguments, {:d} given.".format(nargs)) # process optional arguments: accuracy = kwargs.pop('accuracy', 1.0e-4) @@ -743,8 +734,8 @@ adaptive=False, detect_divergence=False, quiet=False) ##################################################################### ## INITIALIZE ITERATIVE PROCESS: ## ##################################################################### - x0, y0 = self.wcs_world2pix(ra, dec, origin) # <-- initial approximation - # (WCS based only) + x0, y0 = self.wcs_world2pix(ra, dec, origin) # <-- initial approximation + # (WCS based only) # see if an iterative solution is required (when any of the # non-CD-matrix corrections are present). If not required @@ -757,17 +748,17 @@ adaptive=False, detect_divergence=False, quiet=False) if vect1D: return [x0, y0] else: - return np.dstack([x0,y0])[0] + return np.dstack([x0, y0])[0] - x = x0.copy() # 0-order solution - y = y0.copy() # 0-order solution + x = x0.copy() # 0-order solution + y = y0.copy() # 0-order solution # initial correction: dx, dy = self.pix2foc(x, y, origin) # If pix2foc does not apply all the required distortion # corrections then replace the above line with: - #r0, d0 = self.all_pix2world(x, y, origin) - #dx, dy = self.wcs_world2pix(r0, d0, origin ) + # r0, d0 = self.all_pix2world(x, y, origin) + # dx, dy = self.wcs_world2pix(r0, d0, origin ) dx -= x0 dy -= y0 @@ -776,21 +767,21 @@ adaptive=False, detect_divergence=False, quiet=False) y -= dy # norn (L2) squared of the correction: - dn2prev = dx**2+dy**2 - dn2 = dn2prev + dn2prev = dx ** 2 + dy ** 2 + dn2 = dn2prev # prepare for iterative process - iterlist = list(range(1, maxiter+1)) - accuracy2 = accuracy**2 - ind = None - inddiv = None + iterlist = list(range(1, maxiter + 1)) + accuracy2 = accuracy ** 2 + ind = None + inddiv = None - npts = x.shape[0] + npts = x.shape[0] # turn off numpy runtime warnings for 'invalid' and 'over': old_invalid = np.geterr()['invalid'] - old_over = np.geterr()['over'] - np.seterr(invalid = 'ignore', over = 'ignore') + old_over = np.geterr()['over'] + np.seterr(invalid='ignore', over='ignore') ##################################################################### ## NON-ADAPTIVE ITERATIONS: ## @@ -805,13 +796,13 @@ adaptive=False, detect_divergence=False, quiet=False) dx, dy = self.pix2foc(x, y, origin) # If pix2foc does not apply all the required distortion # corrections then replace the above line with: - #r0, d0 = self.all_pix2world(x, y, origin) - #dx, dy = self.wcs_world2pix(r0, d0, origin ) + # r0, d0 = self.all_pix2world(x, y, origin) + # dx, dy = self.wcs_world2pix(r0, d0, origin ) dx -= x0 dy -= y0 # update norn (L2) squared of the correction: - dn2 = dx**2+dy**2 + dn2 = dx ** 2 + dy ** 2 # check for divergence (we do this in two stages # to optimize performance for the most common @@ -826,12 +817,12 @@ adaptive=False, detect_divergence=False, quiet=False) x[ind] -= dx[ind] y[ind] -= dy[ind] # switch to adaptive iterations: - ind, = np.where((dn2 >= accuracy2) & \ - (dn2 <= dn2prev) & np.isfinite(dn2)) + ind, = np.where((dn2 >= accuracy2) & + (dn2 <= dn2prev) & np.isfinite(dn2)) iterlist = iterlist[k:] adaptive = True break - #dn2prev[ind] = dn2[ind] + # dn2prev[ind] = dn2[ind] dn2prev = dn2 # apply correction: @@ -854,22 +845,22 @@ adaptive=False, detect_divergence=False, quiet=False) dx[ind], dy[ind] = self.pix2foc(x[ind], y[ind], origin) # If pix2foc does not apply all the required distortion # corrections then replace the above line with: - #r0[ind], d0[ind] = self.all_pix2world(x[ind], y[ind], origin) - #dx[ind], dy[ind] = self.wcs_world2pix(r0[ind], d0[ind], origin) + # r0[ind], d0[ind] = self.all_pix2world(x[ind], y[ind], origin) + # dx[ind], dy[ind] = self.wcs_world2pix(r0[ind], d0[ind], origin) dx[ind] -= x0[ind] dy[ind] -= y0[ind] # update norn (L2) squared of the correction: - dn2 = dx**2+dy**2 + dn2 = dx ** 2 + dy ** 2 # update indices of elements that still need correction: if detect_divergence: ind, = np.where((dn2 >= accuracy2) & (dn2 <= dn2prev)) - #ind = ind[np.where((dn2[ind] >= accuracy2) & (dn2[ind] <= dn2prev))] + # ind = ind[np.where((dn2[ind] >= accuracy2) & (dn2[ind] <= dn2prev))] dn2prev[ind] = dn2[ind] else: ind, = np.where(dn2 >= accuracy2) - #ind = ind[np.where(dn2[ind] >= accuracy2)] + # ind = ind[np.where(dn2[ind] >= accuracy2)] # apply correction: x[ind] -= dx[ind] @@ -880,9 +871,9 @@ adaptive=False, detect_divergence=False, quiet=False) ## AND FAILED-TO-CONVERGE POINTS ## ##################################################################### # Identify diverging and/or invalid points: - invalid = (((~np.isfinite(y)) | (~np.isfinite(x)) | \ - (~np.isfinite(dn2))) & \ - (np.isfinite(ra)) & (np.isfinite(dec))) + invalid = (((~np.isfinite(y)) | (~np.isfinite(x)) | + (~np.isfinite(dn2))) & + (np.isfinite(ra)) & (np.isfinite(dec))) # When detect_divergence==False, dn2prev is outdated (it is the # norm^2 of the very first correction). Still better than nothing... inddiv, = np.where(((dn2 >= accuracy2) & (dn2 > dn2prev)) | invalid) @@ -891,7 +882,7 @@ adaptive=False, detect_divergence=False, quiet=False) # identify points that did not converge within # 'maxiter' iterations: if k >= maxiter: - ind,= np.where((dn2 >= accuracy2) & (dn2 <= dn2prev) & (~invalid)) + ind, = np.where((dn2 >= accuracy2) & (dn2 <= dn2prev) & (~invalid)) if ind.shape[0] == 0: ind = None else: @@ -907,50 +898,50 @@ adaptive=False, detect_divergence=False, quiet=False) sol = [x, y] err = [np.abs(dx), np.abs(dy)] else: - sol = np.dstack( [x, y] )[0] - err = np.dstack( [np.abs(dx), np.abs(dy)] )[0] + sol = np.dstack([x, y] )[0] + err = np.dstack([np.abs(dx), np.abs(dy)] )[0] # restore previous numpy error settings: - np.seterr(invalid = old_invalid, over = old_over) + np.seterr(invalid=old_invalid, over=old_over) if inddiv is None: - raise NoConvergence("'HSTWCS.all_world2pix' failed to " \ - "converge to the requested accuracy after {:d} " \ - "iterations.".format(k), best_solution = sol, \ - accuracy = err, niter = k, failed2converge = ind, \ - divergent = None) + raise NoConvergence("'HSTWCS.all_world2pix' failed to " + "converge to the requested accuracy after {:d} " + "iterations.".format(k), best_solution=sol, + accuracy=err, niter=k, failed2converge=ind, + divergent=None) else: - raise NoConvergence("'HSTWCS.all_world2pix' failed to " \ - "converge to the requested accuracy.{0:s}" \ - "After {1:d} iterations, the solution is diverging " \ - "at least for one input point." \ - .format(os.linesep, k), best_solution = sol, \ - accuracy = err, niter = k, failed2converge = ind, \ - divergent = inddiv) + raise NoConvergence("'HSTWCS.all_world2pix' failed to " + "converge to the requested accuracy.{0:s}" + "After {1:d} iterations, the solution is diverging " + "at least for one input point." + .format(os.linesep, k), best_solution=sol, + accuracy=err, niter=k, failed2converge=ind, + divergent=inddiv) ##################################################################### ## FINALIZE AND FORMAT DATA FOR RETURN: ## ##################################################################### # restore previous numpy error settings: - np.seterr(invalid = old_invalid, over = old_over) + np.seterr(invalid=old_invalid, over=old_over) if vect1D: return [x, y] else: - return np.dstack( [x, y] )[0] + return np.dstack([x, y] )[0] def _updatehdr(self, ext_hdr): - #kw2add : OCX10, OCX11, OCY10, OCY11 + # kw2add : OCX10, OCX11, OCY10, OCY11 # record the model in the header for use by pydrizzle - ext_hdr['OCX10'] = self.idcmodel.cx[1,0] - ext_hdr['OCX11'] = self.idcmodel.cx[1,1] - ext_hdr['OCY10'] = self.idcmodel.cy[1,0] - ext_hdr['OCY11'] = self.idcmodel.cy[1,1] + ext_hdr['OCX10'] = self.idcmodel.cx[1, 0] + ext_hdr['OCX11'] = self.idcmodel.cx[1, 1] + ext_hdr['OCY10'] = self.idcmodel.cy[1, 0] + ext_hdr['OCY11'] = self.idcmodel.cy[1, 1] ext_hdr['IDCSCALE'] = self.idcmodel.refpix['PSCALE'] ext_hdr['IDCTHETA'] = self.idcmodel.refpix['THETA'] ext_hdr['IDCXREF'] = self.idcmodel.refpix['XREF'] ext_hdr['IDCYREF'] = self.idcmodel.refpix['YREF'] - ext_hdr['IDCV2REF'] = self.idcmodel.refpix['V2REF'] + ext_hdr['IDCV2REF'] = self.idcmodel.refpix['V2REF'] ext_hdr['IDCV3REF'] = self.idcmodel.refpix['V3REF'] def printwcs(self): @@ -958,8 +949,8 @@ adaptive=False, detect_divergence=False, quiet=False) Print the basic WCS keywords. """ print('WCS Keywords\n') - print('CD_11 CD_12: %r %r' % (self.wcs.cd[0,0], self.wcs.cd[0,1])) - print('CD_21 CD_22: %r %r' % (self.wcs.cd[1,0], self.wcs.cd[1,1])) + print('CD_11 CD_12: %r %r' % (self.wcs.cd[0, 0], self.wcs.cd[0, 1])) + print('CD_21 CD_22: %r %r' % (self.wcs.cd[1, 0], self.wcs.cd[1, 1])) print('CRVAL : %r %r' % (self.wcs.crval[0], self.wcs.crval[1])) print('CRPIX : %r %r' % (self.wcs.crpix[0], self.wcs.crpix[1])) print('NAXIS : %d %d' % (self.naxis1, self.naxis2)) diff --git a/stwcs/wcsutil/instruments.py b/stwcs/wcsutil/instruments.py index f662513..83fb1bd 100644 --- a/stwcs/wcsutil/instruments.py +++ b/stwcs/wcsutil/instruments.py @@ -1,6 +1,5 @@ -from __future__ import absolute_import, division, print_function # confidence high +from __future__ import absolute_import, division, print_function -from .mappings import ins_spec_kw class InstrWCS(object): """ @@ -135,13 +134,14 @@ class InstrWCS(object): self.chip = 1 def set_parity(self): - self.parity = [[1.0,0.0],[0.0,-1.0]] + self.parity = [[1.0, 0.0], [0.0, -1.0]] def set_detector(self): # each instrument has a different kw for detector and it can be # in a different header, so this is to be handled by the instrument classes self.detector = 'DEFAULT' + class ACSWCS(InstrWCS): """ get instrument specific kw @@ -150,7 +150,7 @@ class ACSWCS(InstrWCS): def __init__(self, hdr0, hdr): self.primhdr = hdr0 self.exthdr = hdr - InstrWCS.__init__(self,hdr0, hdr) + InstrWCS.__init__(self, hdr0, hdr) self.set_ins_spec_kw() def set_detector(self): @@ -161,9 +161,9 @@ class ACSWCS(InstrWCS): raise 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]]} + 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]]} if self.detector not in list(parity.keys()): parity = InstrWCS.set_parity(self) @@ -173,24 +173,22 @@ class ACSWCS(InstrWCS): class WFPC2WCS(InstrWCS): - def __init__(self, hdr0, hdr): self.primhdr = hdr0 self.exthdr = hdr - InstrWCS.__init__(self,hdr0, 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: + if self.filter1 == " " or self.filter1 is None: self.filter1 = 'CLEAR1' def set_filter2(self): self.filter2 = self.primhdr.get('FILTNAM2', None) - if self.filter2 == " " or self.filter2 == None: + if self.filter2 == " " or self.filter2 is None: self.filter2 = 'CLEAR2' - def set_binned(self): mode = self.primhdr.get('MODE', 1) if mode == 'FULL': @@ -202,7 +200,7 @@ class WFPC2WCS(InstrWCS): self.chip = self.exthdr.get('DETECTOR', 1) def set_parity(self): - self.parity = [[-1.0,0.],[0.,1.0]] + self.parity = [[-1.0, 0.], [0., 1.0]] def set_detector(self): try: @@ -220,7 +218,7 @@ class WFC3WCS(InstrWCS): def __init__(self, hdr0, hdr): self.primhdr = hdr0 self.exthdr = hdr - InstrWCS.__init__(self,hdr0, hdr) + InstrWCS.__init__(self, hdr0, hdr) self.set_ins_spec_kw() def set_detector(self): @@ -232,22 +230,23 @@ class WFC3WCS(InstrWCS): def set_filter1(self): self.filter1 = self.primhdr.get('FILTER', None) - if self.filter1 == " " or self.filter1 == None: + if self.filter1 == " " or self.filter1 is None: self.filter1 = 'CLEAR' def set_filter2(self): - #Nicmos idc tables do not allow 2 filters. + # Nicmos idc tables do not allow 2 filters. self.filter2 = 'CLEAR' def set_parity(self): - parity = {'UVIS':[[-1.0,0.0],[0.0,1.0]], - 'IR':[[-1.0,0.0],[0.0,1.0]]} + parity = {'UVIS': [[-1.0, 0.0], [0.0, 1.0]], + 'IR': [[-1.0, 0.0], [0.0, 1.0]]} if self.detector not in list(parity.keys()): parity = InstrWCS.set_parity(self) else: self.parity = parity[self.detector] + class NICMOSWCS(InstrWCS): """ Create a NICMOS specific class @@ -256,19 +255,19 @@ class NICMOSWCS(InstrWCS): def __init__(self, hdr0, hdr): self.primhdr = hdr0 self.exthdr = hdr - InstrWCS.__init__(self,hdr0, hdr) + InstrWCS.__init__(self, hdr0, hdr) self.set_ins_spec_kw() def set_parity(self): - self.parity = [[-1.0,0.],[0.,1.0]] + self.parity = [[-1.0, 0.], [0., 1.0]] def set_filter1(self): self.filter1 = self.primhdr.get('FILTER', None) - if self.filter1 == " " or self.filter1 == None: + if self.filter1 == " " or self.filter1 is None: self.filter1 = 'CLEAR' def set_filter2(self): - #Nicmos idc tables do not allow 2 filters. + # Nicmos idc tables do not allow 2 filters. self.filter2 = 'CLEAR' def set_chip(self): @@ -281,6 +280,7 @@ class NICMOSWCS(InstrWCS): print('ERROR: Detector kw not found.\n') raise + class STISWCS(InstrWCS): """ A STIS specific class @@ -289,20 +289,20 @@ class STISWCS(InstrWCS): def __init__(self, hdr0, hdr): self.primhdr = hdr0 self.exthdr = hdr - InstrWCS.__init__(self,hdr0, hdr) + InstrWCS.__init__(self, hdr0, hdr) self.set_ins_spec_kw() def set_parity(self): - self.parity = [[-1.0,0.],[0.,1.0]] + self.parity = [[-1.0, 0.], [0., 1.0]] def set_filter1(self): self.filter1 = self.exthdr.get('OPT_ELEM', None) - if self.filter1 == " " or self.filter1 == None: + if self.filter1 == " " or self.filter1 is None: self.filter1 = 'CLEAR1' def set_filter2(self): self.filter2 = self.exthdr.get('FILTER', None) - if self.filter2 == " " or self.filter2 == None: + if self.filter2 == " " or self.filter2 is None: self.filter2 = 'CLEAR2' def set_detector(self): diff --git a/stwcs/wcsutil/mappings.py b/stwcs/wcsutil/mappings.py index 24038bf..a85df13 100644 --- a/stwcs/wcsutil/mappings.py +++ b/stwcs/wcsutil/mappings.py @@ -1,29 +1,29 @@ -from __future__ import division # confidence high +from __future__ import division # This dictionary maps an instrument into an instrument class # The instrument class handles instrument specific keywords -inst_mappings={'WFPC2': 'WFPC2WCS', - 'ACS': 'ACSWCS', - 'NICMOS': 'NICMOSWCS', - 'STIS': 'STISWCS', - 'WFC3': 'WFC3WCS', - 'DEFAULT': 'InstrWCS' - } +inst_mappings = {'WFPC2': 'WFPC2WCS', + 'ACS': 'ACSWCS', + 'NICMOS': 'NICMOSWCS', + 'STIS': 'STISWCS', + 'WFC3': 'WFC3WCS', + 'DEFAULT': 'InstrWCS' + } # A list of instrument specific keywords # Every instrument class must have methods which define each of these # as class attributes. -ins_spec_kw = [ 'idctab', 'offtab', 'date_obs', 'ra_targ', 'dec_targ', 'pav3', \ - 'detector', 'ltv1', 'ltv2', 'parity', 'binned','vafactor', \ - 'chip', 'naxis1', 'naxis2', 'filter1', 'filter2'] +ins_spec_kw = ['idctab', 'offtab', 'date_obs', 'ra_targ', 'dec_targ', + 'pav3', 'detector', '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'] +# 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_', 'CD2_', 'CRVAL', 'CTYPE', 'CRPIX', 'CTYPE', 'CDELT', 'CUNIT'] - +basic_wcs = ['CD1_', 'CD2_', 'CRVAL', 'CTYPE', 'CRPIX', 'CTYPE', + 'CDELT', 'CUNIT'] diff --git a/stwcs/wcsutil/mosaic.py b/stwcs/wcsutil/mosaic.py index 9d2d0a3..a757e0c 100644 --- a/stwcs/wcsutil/mosaic.py +++ b/stwcs/wcsutil/mosaic.py @@ -1,15 +1,17 @@ -from __future__ import division, print_function +from __future__ import absolute_import, division, print_function import numpy as np from matplotlib import pyplot as plt from astropy.io import fits import string from stsci.tools import parseinput, irafglob -from stwcs.distortion import utils -from stwcs import updatewcs, wcsutil -from stwcs.wcsutil import altwcs +from ..distortion import utils +from .. import wcsutil +from ..wcsutil import altwcs -def vmosaic(fnames, outwcs=None, ref_wcs=None, ext=None, extname=None, undistort=True, wkey='V', wname='VirtualMosaic', plot=False, clobber=False): + +def vmosaic(fnames, outwcs=None, ref_wcs=None, ext=None, extname=None, undistort=True, + wkey='V', wname='VirtualMosaic', plot=False, clobber=False): """ Create a virtual mosaic using the WCS of the input images. @@ -61,27 +63,29 @@ def vmosaic(fnames, outwcs=None, ref_wcs=None, ext=None, extname=None, undistort tangent plane and the virtual WCS is recorded in the header. """ wcsobjects = readWCS(fnames, ext, extname) - if outwcs != None: + if outwcs is not None: outwcs = outwcs.deepcopy() else: - if ref_wcs != None: + if ref_wcs is not None: outwcs = utils.output_wcs(wcsobjects, ref_wcs=ref_wcs, undistort=undistort) else: outwcs = utils.output_wcs(wcsobjects, undistort=undistort) if plot: - outc=np.array([[0.,0], [outwcs._naxis1, 0], - [outwcs._naxis1, outwcs._naxis2], - [0, outwcs._naxis2], [0, 0]]) - plt.plot(outc[:,0], outc[:,1]) + outc = np.array([[0., 0], [outwcs._naxis1, 0], + [outwcs._naxis1, outwcs._naxis2], + [0, outwcs._naxis2], [0, 0]]) + plt.plot(outc[:, 0], outc[:, 1]) for wobj in wcsobjects: - outcorners = outwcs.wcs_world2pix(wobj.calc_footprint(),1) + outcorners = outwcs.wcs_world2pix(wobj.calc_footprint(), 1) if plot: - plt.plot(outcorners[:,0], outcorners[:,1]) + plt.plot(outcorners[:, 0], outcorners[:, 1]) objwcs = outwcs.deepcopy() objwcs.wcs.crpix = objwcs.wcs.crpix - (outcorners[0]) - updatehdr(wobj.filename, objwcs,wkey=wkey, wcsname=wname, ext=wobj.extname, clobber=clobber) + updatehdr(wobj.filename, objwcs, wkey=wkey, wcsname=wname, ext=wobj.extname, + clobber=clobber) return outwcs + def updatehdr(fname, wcsobj, wkey, wcsname, ext=1, clobber=False): hdr = fits.getheader(fname, ext=ext) all_keys = list(string.ascii_uppercase) @@ -89,7 +93,9 @@ def updatehdr(fname, wcsobj, wkey, wcsname, ext=1, clobber=False): raise KeyError("wkey must be one character: A-Z") if wkey not in altwcs.available_wcskeys(hdr): if not clobber: - raise ValueError("wkey %s is already in use. Use clobber=True to overwrite it or specify a different key." %wkey) + raise ValueError("wkey {0} is already in use." + "Use clobber=True to overwrite it or" + "specify a different key.".format(wkey)) else: altwcs.deleteWCS(fname, ext=ext, wcskey='V') f = fits.open(fname, mode='update') @@ -98,24 +104,25 @@ def updatehdr(fname, wcsobj, wkey, wcsname, ext=1, clobber=False): wcsnamekey = 'WCSNAME' + wkey f[ext].header[wcsnamekey] = wcsname for k in hwcs: - f[ext].header[k[:7]+wkey] = hwcs[k] + f[ext].header[k[: 7] + wkey] = hwcs[k] f.close() -def wcs2header(wcsobj): +def wcs2header(wcsobj): h = wcsobj.to_header() if wcsobj.wcs.has_cd(): altwcs.pc2cd(h) h['CTYPE1'] = 'RA---TAN' h['CTYPE2'] = 'DEC--TAN' - norient = np.rad2deg(np.arctan2(h['CD1_2'],h['CD2_2'])) - #okey = 'ORIENT%s' % wkey + norient = np.rad2deg(np.arctan2(h['CD1_2'], h['CD2_2'])) + # okey = 'ORIENT%s' % wkey okey = 'ORIENT' h[okey] = norient return h + def readWCS(input, exts=None, extname=None): if isinstance(input, str): if input[0] == '@': @@ -134,15 +141,15 @@ def readWCS(input, exts=None, extname=None): wcso = [] fomited = [] # figure out which FITS extension(s) to use - if exts == None and extname == None: - #Assume it's simple FITS and the data is in the primary HDU + if exts is None and extname is None: + # Assume it's simple FITS and the data is in the primary HDU for f in filelist: try: wcso = wcsutil.HSTWCS(f) except AttributeError: fomited.append(f) continue - elif exts != None and validateExt(exts): + elif exts is not None and validateExt(exts): exts = [exts] for f in filelist: try: @@ -150,7 +157,7 @@ def readWCS(input, exts=None, extname=None): except KeyError: fomited.append(f) continue - elif extname != None: + elif extname is not None: for f in filelist: fobj = fits.open(f) for i in range(len(fobj)): @@ -159,7 +166,7 @@ def readWCS(input, exts=None, extname=None): except KeyError: continue if ename.lower() == extname.lower(): - wcso.append(wcsutil.HSTWCS(f,ext=i)) + wcso.append(wcsutil.HSTWCS(f, ext=i)) else: continue fobj.close() @@ -179,5 +186,3 @@ def validateExt(ext): else: return True - - diff --git a/stwcs/wcsutil/wcscorr.py b/stwcs/wcsutil/wcscorr.py index 3f9b7d5..48a1978 100644 --- a/stwcs/wcsutil/wcscorr.py +++ b/stwcs/wcsutil/wcscorr.py @@ -1,21 +1,20 @@ from __future__ import absolute_import, division, print_function -import os,copy +import copy import numpy as np from astropy.io import fits import stwcs -from stwcs.wcsutil import altwcs -from stwcs.updatewcs import utils +from . import altwcs +from ..updatewcs import utils from stsci.tools import fileutil -from . import convertwcs -DEFAULT_WCS_KEYS = ['CRVAL1','CRVAL2','CRPIX1','CRPIX2', - 'CD1_1','CD1_2','CD2_1','CD2_2', - 'CTYPE1','CTYPE2','ORIENTAT'] -DEFAULT_PRI_KEYS = ['HDRNAME','SIPNAME','NPOLNAME','D2IMNAME','DESCRIP'] -COL_FITSKW_DICT = {'RMS_RA':'sci.crder1','RMS_DEC':'sci.crder2', - 'NMatch':'sci.nmatch','Catalog':'sci.catalog'} +DEFAULT_WCS_KEYS = ['CRVAL1', 'CRVAL2', 'CRPIX1', 'CRPIX2', + 'CD1_1', 'CD1_2', 'CD2_1', 'CD2_2', + 'CTYPE1', 'CTYPE2', 'ORIENTAT'] +DEFAULT_PRI_KEYS = ['HDRNAME', 'SIPNAME', 'NPOLNAME', 'D2IMNAME', 'DESCRIP'] +COL_FITSKW_DICT = {'RMS_RA': 'sci.crder1', 'RMS_DEC': 'sci.crder2', + 'NMatch': 'sci.nmatch', 'Catalog': 'sci.catalog'} ### ### WCSEXT table related keyword archive functions @@ -53,7 +52,7 @@ def init_wcscorr(input, force=False): return else: del fimg['wcscorr'] - print('Initializing new WCSCORR table for ',fimg.filename()) + print('Initializing new WCSCORR table for ', fimg.filename()) used_wcskeys = altwcs.wcskeys(fimg['SCI', 1].header) @@ -64,14 +63,14 @@ def init_wcscorr(input, force=False): # create new table with more rows than needed initially to make it easier to # add new rows later - wcsext = create_wcscorr(descrip=True,numrows=numsci, padding=(numsci*numwcs) + numsci * 4) + wcsext = create_wcscorr(descrip=True, numrows=numsci, padding=(numsci * numwcs) + numsci * 4) # Assign the correct EXTNAME value to this table extension wcsext.header['TROWS'] = (numsci * 2, 'Number of updated rows in table') wcsext.header['EXTNAME'] = ('WCSCORR', 'Table with WCS Update history') wcsext.header['EXTVER'] = 1 # define set of WCS keywords which need to be managed and copied to the table - wcs1 = stwcs.wcsutil.HSTWCS(fimg,ext=('SCI',1)) + wcs1 = stwcs.wcsutil.HSTWCS(fimg, ext=('SCI', 1)) idc2header = True if wcs1.idcscale is None: idc2header = False @@ -79,18 +78,18 @@ def init_wcscorr(input, force=False): prihdr = fimg[0].header prihdr_keys = DEFAULT_PRI_KEYS - pri_funcs = {'SIPNAME':stwcs.updatewcs.utils.build_sipname, - 'NPOLNAME':stwcs.updatewcs.utils.build_npolname, - 'D2IMNAME':stwcs.updatewcs.utils.build_d2imname} + pri_funcs = {'SIPNAME': stwcs.updatewcs.utils.build_sipname, + 'NPOLNAME': stwcs.updatewcs.utils.build_npolname, + 'D2IMNAME': stwcs.updatewcs.utils.build_d2imname} # Now copy original OPUS values into table for extver in range(1, numsci + 1): rowind = find_wcscorr_row(wcsext.data, {'WCS_ID': 'OPUS', 'EXTVER': extver, - 'WCS_key':'O'}) + 'WCS_key': 'O'}) # There should only EVER be a single row for each extension with OPUS values rownum = np.where(rowind)[0][0] - #print 'Archiving OPUS WCS in row number ',rownum,' in WCSCORR table for SCI,',extver + # print 'Archiving OPUS WCS in row number ',rownum,' in WCSCORR table for SCI,',extver hdr = fimg['SCI', extver].header # define set of WCS keywords which need to be managed and copied to the table @@ -100,7 +99,7 @@ def init_wcscorr(input, force=False): # if so, get its values directly, otherwise, archive the PRIMARY WCS # as the OPUS values in the WCSCORR table if 'O' not in used_wcskeys: - altwcs.archiveWCS(fimg,('SCI',extver),wcskey='O', wcsname='OPUS') + altwcs.archiveWCS(fimg, ('SCI', extver), wcskey='O', wcsname='OPUS') wkey = 'O' wcs = stwcs.wcsutil.HSTWCS(fimg, ext=('SCI', extver), wcskey=wkey) @@ -113,7 +112,7 @@ def init_wcscorr(input, force=False): break for key in wcs_keywords: if key in wcsext.data.names: - wcsext.data.field(key)[rownum] = wcshdr[(key+wkey)[:8]] + wcsext.data.field(key)[rownum] = wcshdr[(key + wkey)[:8]] # Now get any keywords from PRIMARY header needed for WCS updates for key in prihdr_keys: if key in prihdr: @@ -143,7 +142,7 @@ def init_wcscorr(input, force=False): # identify next empty row rowind = find_wcscorr_row(wcsext.data, - selections={'wcs_id':['','0.0']}) + selections={'wcs_id': ['', '0.0']}) rows = np.where(rowind) if len(rows[0]) > 0: rownum = np.where(rowind)[0][0] @@ -196,26 +195,26 @@ def find_wcscorr_row(wcstab, selections): mask = None for i in selections: icol = wcstab.field(i) - if isinstance(icol,np.chararray): icol = icol.rstrip() + if isinstance(icol, np.chararray): icol = icol.rstrip() selecti = selections[i] - if not isinstance(selecti,list): - if isinstance(selecti,str): + if not isinstance(selecti, list): + if isinstance(selecti, str): selecti = selecti.rstrip() bmask = (icol == selecti) if mask is None: mask = bmask.copy() else: - mask = np.logical_and(mask,bmask) + mask = np.logical_and(mask, bmask) del bmask else: for si in selecti: - if isinstance(si,str): + if isinstance(si, str): si = si.rstrip() bmask = (icol == si) if mask is None: mask = bmask.copy() else: - mask = np.logical_or(mask,bmask) + mask = np.logical_or(mask, bmask) del bmask return mask @@ -268,7 +267,7 @@ def update_wcscorr(dest, source=None, extname='SCI', wcs_id=None, active=True): headerlet """ if not isinstance(dest, fits.HDUList): - dest = fits.open(dest,mode='update') + dest = fits.open(dest, mode='update') fname = dest.filename() if source is None: @@ -280,7 +279,7 @@ def update_wcscorr(dest, source=None, extname='SCI', wcs_id=None, active=True): numext = fileutil.countExtn(source, extname) if numext == 0: raise ValueError('No %s extensions found in the source HDU list.' - % extname) + % extname) # Initialize the WCSCORR table extension in dest if not already present init_wcscorr(dest) try: @@ -291,7 +290,7 @@ def update_wcscorr(dest, source=None, extname='SCI', wcs_id=None, active=True): # check to see whether or not this is an up-to-date table # replace with newly initialized table with current format old_table = dest['WCSCORR'] - wcscorr_cols = ['WCS_ID','EXTVER', 'SIPNAME', + wcscorr_cols = ['WCS_ID', 'EXTVER', 'SIPNAME', 'HDRNAME', 'NPOLNAME', 'D2IMNAME'] for colname in wcscorr_cols: @@ -308,7 +307,7 @@ def update_wcscorr(dest, source=None, extname='SCI', wcs_id=None, active=True): # modified... wcs_keys = altwcs.wcskeys(source[(extname, 1)].header) wcs_keys = [kk for kk in wcs_keys if kk] - if ' ' not in wcs_keys: wcs_keys.append(' ') # Insure that primary WCS gets used + if ' ' not in wcs_keys: wcs_keys.append(' ') # Insure that primary WCS gets used # apply logic for only updating WCSCORR table with specified keywords # corresponding to the WCS with WCSNAME=wcs_id if wcs_id is not None: @@ -324,10 +323,10 @@ def update_wcscorr(dest, source=None, extname='SCI', wcs_id=None, active=True): wcs_keywords = list(wcshdr.keys()) if 'O' in wcs_keys: - wcs_keys.remove('O') # 'O' is reserved for original OPUS WCS + wcs_keys.remove('O') # 'O' is reserved for original OPUS WCS # create new table for hdr and populate it with the newly updated values - new_table = create_wcscorr(descrip=True,numrows=0, padding=len(wcs_keys)*numext) + new_table = create_wcscorr(descrip=True, numrows=0, padding=len(wcs_keys) * numext) prihdr = source[0].header # Get headerlet related keywords here @@ -344,18 +343,18 @@ def update_wcscorr(dest, source=None, extname='SCI', wcs_id=None, active=True): for extver in range(1, numext + 1): extn = (extname, extver) if 'SIPWCS' in extname and not active: - tab_extver = 0 # Since it has not been added to the SCI header yet + tab_extver = 0 # Since it has not been added to the SCI header yet else: tab_extver = extver hdr = source[extn].header - if 'WCSNAME'+wcs_key in hdr: + if 'WCSNAME' + wcs_key in hdr: wcsname = hdr['WCSNAME' + wcs_key] else: wcsname = utils.build_default_wcsname(hdr['idctab']) selection = {'WCS_ID': wcsname, 'EXTVER': tab_extver, - 'SIPNAME':sipname, 'HDRNAME': hdrname, - 'NPOLNAME': npolname, 'D2IMNAME':d2imname + 'SIPNAME': sipname, 'HDRNAME': hdrname, + 'NPOLNAME': npolname, 'D2IMNAME': d2imname } # Ensure that an entry for this WCS is not already in the dest @@ -387,15 +386,14 @@ def update_wcscorr(dest, source=None, extname='SCI', wcs_id=None, active=True): # Interpret any 'pri.hdrname' or # 'sci.crpix1' formatted keyword names if '.' in fitkw: - srchdr,fitkw = fitkw.split('.') + srchdr, fitkw = fitkw.split('.') if 'pri' in srchdr.lower(): srchdr = prihdr else: srchdr = source[extn].header else: srchdr = source[extn].header - if fitkw+wcs_key in srchdr: - new_table.data.field(key)[idx] = srchdr[fitkw+wcs_key] - + if fitkw + wcs_key in srchdr: + new_table.data.field(key)[idx] = srchdr[fitkw + wcs_key] # If idx was never incremented, no rows were added, so there's nothing else # to do... @@ -403,16 +401,16 @@ def update_wcscorr(dest, source=None, extname='SCI', wcs_id=None, active=True): return # Now, we need to merge this into the existing table - rowind = find_wcscorr_row(old_table.data, {'wcs_id':['','0.0']}) + rowind = find_wcscorr_row(old_table.data, {'wcs_id': ['', '0.0']}) old_nrows = np.where(rowind)[0][0] new_nrows = new_table.data.shape[0] # check to see if there is room for the new row - if (old_nrows + new_nrows) > old_table.data.shape[0]-1: + if (old_nrows + new_nrows) > old_table.data.shape[0] - 1: pad_rows = 2 * new_nrows # if not, create a new table with 'pad_rows' new empty rows - upd_table = fits.new_table(old_table.columns,header=old_table.header, - nrows=old_table.data.shape[0]+pad_rows) + upd_table = fits.new_table(old_table.columns, header=old_table.header, + nrows=old_table.data.shape[0] + pad_rows) else: upd_table = old_table pad_rows = 0 @@ -421,10 +419,10 @@ def update_wcscorr(dest, source=None, extname='SCI', wcs_id=None, active=True): if name in new_table.data.names: # reset the default values to ones specific to the row definitions for i in range(pad_rows): - upd_table.data.field(name)[old_nrows+i] = old_table.data.field(name)[-1] + upd_table.data.field(name)[old_nrows + i] = old_table.data.field(name)[-1] # Now populate with values from new table upd_table.data.field(name)[old_nrows:old_nrows + new_nrows] = \ - new_table.data.field(name) + new_table.data.field(name) upd_table.header['TROWS'] = old_nrows + new_nrows # replace old extension with newly updated table extension @@ -447,14 +445,14 @@ def restore_file_from_wcscorr(image, id='OPUS', wcskey=''): wcs_table = fimg['WCSCORR'] orig_rows = (wcs_table.data.field('WCS_ID') == 'OPUS') # create an HSTWCS object to figure out what WCS keywords need to be updated - wcsobj = stwcs.wcsutil.HSTWCS(fimg,ext=('sci',1)) + wcsobj = stwcs.wcsutil.HSTWCS(fimg, ext=('sci', 1)) wcshdr = wcsobj.wcs2header() - for extn in range(1,numsci+1): + for extn in range(1, numsci + 1): # find corresponding row from table ext_rows = (wcs_table.data.field('EXTVER') == extn) - erow = np.where(np.logical_and(ext_rows,orig_rows))[0][0] + erow = np.where(np.logical_and(ext_rows, orig_rows))[0][0] for key in wcshdr: - if key in wcs_table.data.names: # insure that keyword is column in table + if key in wcs_table.data.names: # insure that keyword is column in table tkey = key if 'orient' in key.lower(): @@ -462,14 +460,14 @@ def restore_file_from_wcscorr(image, id='OPUS', wcskey=''): if wcskey == '': skey = key else: - skey = key[:7]+wcskey - fimg['sci',extn].header[skey] = wcs_table.data.field(tkey)[erow] + skey = key[:7] + wcskey + fimg['sci', extn].header[skey] = wcs_table.data.field(tkey)[erow] for key in DEFAULT_PRI_KEYS: if key in wcs_table.data.names: if wcskey == '': pkey = key else: - pkey = key[:7]+wcskey + pkey = key[: 7] + wcskey fimg[0].header[pkey] = wcs_table.data.field(key)[erow] utils.updateNEXTENDKw(fimg) @@ -498,25 +496,25 @@ def create_wcscorr(descrip=False, numrows=1, padding=0): def_float64_zeros = np.array([0.0] * trows, dtype=np.float64) def_float64_ones = def_float64_zeros + 1.0 def_float_col = {'format': 'D', 'array': def_float64_zeros.copy()} - def_float1_col = {'format': 'D', 'array':def_float64_ones.copy()} + def_float1_col = {'format': 'D', 'array': def_float64_ones.copy()} def_str40_col = {'format': '40A', 'array': np.array([''] * trows, dtype='S40')} def_str24_col = {'format': '24A', 'array': np.array([''] * trows, dtype='S24')} def_int32_col = {'format': 'J', - 'array': np.array([0]*trows,dtype=np.int32)} + 'array': np.array([0] * trows, dtype=np.int32)} # If more columns are needed, simply add their definitions to this list - col_names = [('HDRNAME', def_str24_col), ('SIPNAME', def_str24_col), - ('NPOLNAME', def_str24_col), ('D2IMNAME', def_str24_col), - ('CRVAL1', def_float_col), ('CRVAL2', def_float_col), - ('CRPIX1', def_float_col), ('CRPIX2', def_float_col), - ('CD1_1', def_float_col), ('CD1_2', def_float_col), - ('CD2_1', def_float_col), ('CD2_2', def_float_col), - ('CTYPE1', def_str24_col), ('CTYPE2', def_str24_col), - ('ORIENTAT', def_float_col), ('PA_V3', def_float_col), - ('RMS_RA', def_float_col), ('RMS_Dec', def_float_col), - ('NMatch', def_int32_col), ('Catalog', def_str40_col)] + col_names = [('HDRNAME', def_str24_col), ('SIPNAME', def_str24_col), + ('NPOLNAME', def_str24_col), ('D2IMNAME', def_str24_col), + ('CRVAL1', def_float_col), ('CRVAL2', def_float_col), + ('CRPIX1', def_float_col), ('CRPIX2', def_float_col), + ('CD1_1', def_float_col), ('CD1_2', def_float_col), + ('CD2_1', def_float_col), ('CD2_2', def_float_col), + ('CTYPE1', def_str24_col), ('CTYPE2', def_str24_col), + ('ORIENTAT', def_float_col), ('PA_V3', def_float_col), + ('RMS_RA', def_float_col), ('RMS_Dec', def_float_col), + ('NMatch', def_int32_col), ('Catalog', def_str40_col)] # Define selector columns id_col = fits.Column(name='WCS_ID', format='40A', @@ -529,7 +527,7 @@ def create_wcscorr(descrip=False, numrows=1, padding=0): array=np.array(['O'] * numrows + [''] * padding, dtype='S')) # create list of remaining columns to be added to table - col_list = [id_col, extver_col, wcskey_col] # start with selector columns + col_list = [id_col, extver_col, wcskey_col] # start with selector columns for c in col_names: cdef = copy.deepcopy(c[1]) @@ -552,7 +550,8 @@ def create_wcscorr(descrip=False, numrows=1, padding=0): return newtab -def delete_wcscorr_row(wcstab,selections=None,rows=None): + +def delete_wcscorr_row(wcstab, selections=None, rows=None): """ Sets all values in a specified row or set of rows to default values @@ -588,7 +587,7 @@ def delete_wcscorr_row(wcstab,selections=None,rows=None): rowind = find_wcscorr_row(wcstab, selections=selections) delete_rows = np.where(rowind)[0].tolist() else: - if not isinstance(rows,list): + if not isinstance(rows, list): rows = [rows] delete_rows = rows @@ -601,13 +600,14 @@ def delete_wcscorr_row(wcstab,selections=None,rows=None): return # identify next empty row - rowind = find_wcscorr_row(wcstab, selections={'wcs_id':['','0.0']}) + rowind = find_wcscorr_row(wcstab, selections={'wcs_id': ['', '0.0']}) last_blank_row = np.where(rowind)[0][-1] # copy values from blank row into user-specified rows for colname in wcstab.names: wcstab[colname][delete_rows] = wcstab[colname][last_blank_row] + def update_wcscorr_column(wcstab, column, values, selections=None, rows=None): """ Update the values in 'column' with 'values' for selected rows @@ -645,7 +645,7 @@ def update_wcscorr_column(wcstab, column, values, selections=None, rows=None): rowind = find_wcscorr_row(wcstab, selections=selections) update_rows = np.where(rowind)[0].tolist() else: - if not isinstance(rows,list): + if not isinstance(rows, list): rows = [rows] update_rows = rows @@ -654,8 +654,8 @@ def update_wcscorr_column(wcstab, column, values, selections=None, rows=None): # Expand single input value to apply to all selected rows if len(values) > 1 and len(values) < len(update_rows): - print('ERROR: Number of new values',len(values)) - print(' does not match number of rows',len(update_rows),' to be updated!') + print('ERROR: Number of new values', len(values)) + print(' does not match number of rows', len(update_rows), ' to be updated!') print(' Please enter either 1 value or the same number of values') print(' as there are rows to be updated.') print(' Table will not be updated...') diff --git a/stwcs/wcsutil/wcsdiff.py b/stwcs/wcsutil/wcsdiff.py index cfc2d66..9f00f1a 100644 --- a/stwcs/wcsutil/wcsdiff.py +++ b/stwcs/wcsutil/wcsdiff.py @@ -1,10 +1,12 @@ -from __future__ import print_function +from __future__ import absolute_import, print_function + from astropy import wcs as pywcs from collections import OrderedDict from astropy.io import fits from .headerlet import parse_filename import numpy as np + def is_wcs_identical(scifile, file2, sciextlist, fextlist, scikey=" ", file2key=" ", verbose=False): """ @@ -48,28 +50,29 @@ def is_wcs_identical(scifile, file2, sciextlist, fextlist, scikey=" ", sciobj, sciname, close_scifile = parse_filename(scifile) diff['file_names'] = [scifile, file2] if get_rootname(scifile) != get_rootname(file2): - #logger.info('Rootnames do not match.') - diff['rootname'] = ("%s: %s", "%s: %s") % (sciname, get_rootname(scifile), file2, get_rootname(file2)) + # logger.info('Rootnames do not match.') + diff['rootname'] = ("%s: %s", "%s: %s") % (sciname, get_rootname(scifile), file2, + get_rootname(file2)) result = False for i, j in zip(sciextlist, fextlist): w1 = pywcs.WCS(sciobj[i].header, sciobj, key=scikey) w2 = pywcs.WCS(fobj[j].header, fobj, key=file2key) diff['extension'] = [get_extname_extnum(sciobj[i]), get_extname_extnum(fobj[j])] if not np.allclose(w1.wcs.crval, w2.wcs.crval, rtol=10**(-7)): - #logger.info('CRVALs do not match') + # logger.info('CRVALs do not match') diff['CRVAL'] = w1.wcs.crval, w2.wcs.crval result = False if not np.allclose(w1.wcs.crpix, w2.wcs.crpix, rtol=10**(-7)): - #logger.info('CRPIX do not match') - diff ['CRPIX'] = w1.wcs.crpix, w2.wcs.crpix + # logger.info('CRPIX do not match') + diff['CRPIX'] = w1.wcs.crpix, w2.wcs.crpix result = False if not np.allclose(w1.wcs.cd, w2.wcs.cd, rtol=10**(-7)): - #logger.info('CDs do not match') - diff ['CD'] = w1.wcs.cd, w2.wcs.cd + # logger.info('CDs do not match') + diff['CD'] = w1.wcs.cd, w2.wcs.cd result = False if not (np.array(w1.wcs.ctype) == np.array(w2.wcs.ctype)).all(): - #logger.info('CTYPEs do not match') - diff ['CTYPE'] = w1.wcs.ctype, w2.wcs.ctype + # logger.info('CTYPEs do not match') + diff['CTYPE'] = w1.wcs.ctype, w2.wcs.ctype result = False if w1.sip or w2.sip: if (w2.sip and not w1.sip) or (w1.sip and not w2.sip): @@ -79,41 +82,41 @@ def is_wcs_identical(scifile, file2, sciextlist, fextlist, scikey=" ", diff['SIP_A'] = 'SIP_A differ' result = False if not np.allclose(w1.sip.b, w2.sip.b, rtol=10**(-7)): - #logger.info('SIP coefficients do not match') - diff ['SIP_B'] = (w1.sip.b, w2.sip.b) + # logger.info('SIP coefficients do not match') + diff['SIP_B'] = (w1.sip.b, w2.sip.b) result = False if w1.cpdis1 or w2.cpdis1: if w1.cpdis1 and not w2.cpdis1 or w2.cpdis1 and not w1.cpdis1: diff['CPDIS1'] = "CPDIS1 missing" - result=False + result = False if w1.cpdis2 and not w2.cpdis2 or w2.cpdis2 and not w1.cpdis2: diff['CPDIS2'] = "CPDIS2 missing" result = False if not np.allclose(w1.cpdis1.data, w2.cpdis1.data, rtol=10**(-7)): - #logger.info('NPOL distortions do not match') - diff ['CPDIS1_data'] = (w1.cpdis1.data, w2.cpdis1.data) + # logger.info('NPOL distortions do not match') + diff['CPDIS1_data'] = (w1.cpdis1.data, w2.cpdis1.data) result = False if not np.allclose(w1.cpdis2.data, w2.cpdis2.data, rtol=10**(-7)): - #logger.info('NPOL distortions do not match') - diff ['CPDIS2_data'] = (w1.cpdis2.data, w2.cpdis2.data) + # logger.info('NPOL distortions do not match') + diff['CPDIS2_data'] = (w1.cpdis2.data, w2.cpdis2.data) result = False if w1.det2im1 or w2.det2im1: if w1.det2im1 and not w2.det2im1 or \ - w2.det2im1 and not w1.det2im1: + w2.det2im1 and not w1.det2im1: diff['DET2IM'] = "Det2im1 missing" result = False if not np.allclose(w1.det2im1.data, w2.det2im1.data, rtol=10**(-7)): - #logger.info('Det2Im corrections do not match') - diff ['D2IM1_data'] = (w1.det2im1.data, w2.det2im1.data) - result = False + # logger.info('Det2Im corrections do not match') + diff['D2IM1_data'] = (w1.det2im1.data, w2.det2im1.data) + result = False if w1.det2im2 or w2.det2im2: if w1.det2im2 and not w2.det2im2 or \ w2.det2im2 and not w1.det2im2: diff['DET2IM2'] = "Det2im2 missing" result = False if not np.allclose(w1.det2im2.data, w2.det2im2.data, rtol=10**(-7)): - #logger.info('Det2Im corrections do not match') - diff ['D2IM2_data'] = (w1.det2im2.data, w2.det2im2.data) + # logger.info('Det2Im corrections do not match') + diff['D2IM2_data'] = (w1.det2im2.data, w2.det2im2.data) result = False if not result and verbose: for key in diff: @@ -124,6 +127,7 @@ def is_wcs_identical(scifile, file2, sciextlist, fextlist, scikey=" ", sciobj.close() return result, diff + def get_rootname(fname): """ Returns the value of ROOTNAME or DESTIM @@ -139,12 +143,13 @@ def get_rootname(fname): rootname = fname return rootname + def get_extname_extnum(ext): """ Return (EXTNAME, EXTNUM) of a FITS extension """ extname = "" - extnum=1 + extnum = 1 extname = ext.header.get('EXTNAME', extname) extnum = ext.header.get('EXTVER', extnum) return (extname, extnum) |