pydrizzle.drutil

1 """ 2 Utility functions for PyDrizzle that rely on PyRAF's interface to IRAF 3 tasks. 4 """ 5 # 6 # Revision History: 7 # Nov 2001: Added function (getLTVOffsets) for reading subarray offsets 8 # from extension header's LTV keywords. WJH 9 # Mar 2002: Restructured to only contain IRAF based functions. 10 # 20-Sept-2002: Replaced all calls to 'keypar' with calls to 'hselect'. 11 # This eliminates any parameter writing, making it safer for 12 # pipeline/multi-process use. 13 # 24-Sept-2003: Replaced all calls to hselect with calls to 14 # fileutil.getKeyword() to remove dependencies on IRAF. 15 # Also, replaced 'yes' and 'no' with Python Bool vars. 16 # 17 18 import string,os,types 19 from math import ceil,floor 20 21 import numpy as N 22 from numpy import linalg 23 24 from pytools import fileutil 25 from pytools.fileutil import buildRotMatrix 26 27 # Convenience definitions 28 DEGTORAD = fileutil.DEGTORAD 29 30 no = False 31 yes = True 32 33 # Constants 34 IDCTAB = 1 35 DRIZZLE = 2 36 TRAUGER = 3 37 38 try: 39 DEFAULT_IDCDIR = fileutil.osfn('stsdas$pkg/analysis/dither/drizzle/coeffs/') 40 except: 41 DEFAULT_IDCDIR = os.getcwd() 42 43

44 -def factorial(n):

45 """ Compute a factorial for integer n. """ 46 m = 1 47 for i in range(int(n)): 48 m = m * (i+1) 49 return m

50

51 -def combin(j,n):

52 """ Return the combinatorial factor for j in n.""" 53 return (factorial(j) / (factorial(n) * factorial( (j-n) ) ) )

54 55 56 ################# 57 # 58 # 59 # Utility Functions based on IRAF 60 # 61 # 62 #################

63 -def findNumExt(filename):

64 # Open the file given by 'rootname' and return the 65 # number of extensions written out for the observation. 66 # It will be up to the Exposure classes to determine how 67 # many IMSETS they have to work with. 68 # 69 # Only look in the primary extension or first group. 70 #_s = iraf.hselect(filename,'NEXTEND',expr='yes',Stdout=1) 71 _s = fileutil.getKeyword(filename,keyword='NEXTEND') 72 if not _s: 73 _s = fileutil.getKeyword(filename,keyword='GCOUNT') 74 # This may need to be changed to support simple image without 75 # extensions (such as simple FITS images). 76 if _s == '': 77 raise ValueError,"There are NO extensions to be read in this image!" 78 79 return _s

80 81 # Utility function to extract subarray offsets 82 # from LTV keywords. Could be expanded to use 83 # different keywords for some cases. 84 # Added: 7-Nov-2001 WJH

85 -def getLTVOffsets(rootname,header=None):

86 87 _ltv1 = None 88 _ltv2 = None 89 if header: 90 if header.has_key('LTV1'): 91 _ltv1 = header['LTV1'] 92 if header.has_key('LTV2'): 93 _ltv2 = header['LTV2'] 94 else: 95 _ltv1 = fileutil.getKeyword(rootname,'LTV1') 96 _ltv2 = fileutil.getKeyword(rootname,'LTV2') 97 98 if _ltv1 == None: _ltv1 = 0. 99 if _ltv2 == None: _ltv2 = 0. 100 101 return _ltv1,_ltv2

102

103 -def getChipId(header):

104 105 if header.has_key('CCDCHIP'): 106 chip = int(header['CCDCHIP']) 107 elif header.has_key('DETECTOR') and str(header['DETECTOR']).isdigit(): 108 chip = int(header['DETECTOR']) 109 elif header.has_key('CAMERA') and str(header['CAMERA']).isdigit(): 110 chip = int(header['CAMERA']) 111 else: 112 chip = 1 113 114 return chip

115 116

117 -def getIDCFile(image,keyword=None,directory=None):

118 # Open the primary header of the file and read the name of 119 # the IDCTAB. 120 # Parameters: 121 # header - primary and extension header object to read IDCTAB info 122 # 123 # keyword(optional) - header keyword with name of IDCTAB 124 # --OR-- 'HEADER' if need to build IDCTAB name from scratch 125 # (default value: 'IDCTAB') 126 # directory(optional) - directory with default drizzle coeffs tables 127 # (default value: 'drizzle$coeffs') 128 # 129 # This function needs to be generalized to support 130 # keyword='HEADER', as would be the case for WFPC2 data. 131 # 132 133 134 if isinstance(image, types.StringType): 135 # We were provided an image name, so read in the header... 136 header = fileutil.getHeader(image) 137 else: 138 # otherwise, we were provided an image header we can work with directly 139 header = image 140 141 if keyword.lower() == 'header': 142 idcfile,idctype = __getIDCTAB(header) 143 if (idcfile == None): 144 idcfile,idctype = __buildIDCTAB(header,directory) 145 146 elif string.lower(keyword) == 'idctab': 147 # keyword specifies header keyword with IDCTAB name 148 idcfile,idctype = __getIDCTAB(header) 149 150 elif keyword == '': 151 idcfile = None 152 idctype = None 153 else: 154 # Need to build IDCTAB filename from scratch 155 idcfile,idctype = __buildIDCTAB(header,directory,kw = keyword) 156 157 # Account for possible absence of IDCTAB name in header 158 if idcfile == 'N/A': 159 idcfile = None 160 161 if idcfile != None and idcfile != '': 162 # Now we need to recursively expand any IRAF symbols to full paths... 163 idcfile = fileutil.osfn(idcfile) 164 165 if idcfile == None: 166 print 'WARNING: No valid distortion coefficients available!' 167 print 'Using default unshifted, unscaled, unrotated model.' 168 169 return idcfile,idctype

170 171

172 -def __buildIDCTAB(header, directory, kw = None):

173 # Need to build IDCTAB filename from scratch 174 instrument = header['INSTRUME'] 175 if instrument != 'NICMOS': 176 detector = header['DETECTOR'] 177 else: 178 detector = str(header['CAMERA']) 179 180 # Default non-IDCTAB distortion model 181 if (kw == None): 182 keyword = 'cubic' 183 else : 184 keyword = kw 185 186 if not directory: 187 default_dir = DEFAULT_IDCDIR 188 else: 189 default_dir = directory 190 191 if instrument == 'WFPC2': 192 if detector == 1: 193 detname = 'pc' 194 else: 195 detname = 'wf' 196 idcfile = default_dir+detname+str(detector)+'-'+string.lower(keyword) 197 198 elif instrument == 'STIS': 199 idcfile = default_dir+'stis-'+string.lower(detector) 200 201 elif instrument == 'NICMOS': 202 if detector != None: 203 idcfile = default_dir+'nic-'+detector 204 else: 205 idcfile = None 206 else: 207 idcfile = None 208 209 idctype = getIDCFileType(fileutil.osfn(idcfile)) 210 211 return idcfile,idctype

212

213 -def __getIDCTAB(header):

214 # keyword specifies header keyword with IDCTAB name 215 try: 216 idcfile = header['idctab'] 217 except: 218 print 'Warning: No IDCTAB specified in header!' 219 idcfile = None 220 221 return idcfile,'idctab'

222

223 -def getIDCFileType(idcfile):

224 """ Open ASCII IDCFILE to determine the type: cubic,trauger,... """ 225 if idcfile == None: 226 return None 227 228 ifile = open(idcfile,'r') 229 # Search for the first line of the coefficients 230 _line = fileutil.rAsciiLine(ifile) 231 232 # Search for first non-commented line... 233 while _line[0] == '#': 234 _line = fileutil.rAsciiLine(ifile) 235 236 _type = string.rstrip(_line.lower()) 237 238 if _type in ['cubic','quartic','quintic'] or _type.find('poly') > -1: 239 _type = 'cubic' 240 elif _type == 'trauger': 241 _type = 'trauger' 242 else: 243 _type = None 244 245 ifile.close() 246 del ifile 247 248 return _type

249 250 # 251 # Function to read in coefficients from ASCII Cubic Drizzle 252 # Coefficients table. 253 # 254

255 -def readCubicTable(idcfile):

256 # Assumption: this will only be used for cubic file... 257 order = 3 258 # Also, this function does NOT perform any scaling on 259 # the coefficients, it simply passes along what is found 260 # in the file as is... 261 262 # Return a default geometry model if no coefficients filename 263 # is given. This model will not distort the data in any way. 264 if idcfile == None: 265 return fileutil.defaultModel() 266 267 ifile = open(idcfile,'r') 268 # Search for the first line of the coefficients 269 _line = fileutil.rAsciiLine(ifile) 270 271 _found = no 272 while _found == no: 273 if _line[:7] in ['cubic','quartic','quintic'] or _line[:4] == 'poly': 274 found = yes 275 break 276 _line = fileutil.rAsciiLine(ifile) 277 278 # Read in each row of coefficients, without line breaks or newlines 279 # split them into their values, and create a list for A coefficients 280 # and another list for the B coefficients 281 _line = fileutil.rAsciiLine(ifile) 282 a_coeffs = string.split(_line) 283 284 x0 = float(a_coeffs[0]) 285 _line = fileutil.rAsciiLine(ifile) 286 a_coeffs[len(a_coeffs):] = string.split(_line) 287 # Scale coefficients for use within PyDrizzle 288 for i in range(len(a_coeffs)): 289 a_coeffs[i] = float(a_coeffs[i]) 290 291 _line = fileutil.rAsciiLine(ifile) 292 b_coeffs = string.split(_line) 293 y0 = float(b_coeffs[0]) 294 _line = fileutil.rAsciiLine(ifile) 295 b_coeffs[len(b_coeffs):] = string.split(_line) 296 # Scale coefficients for use within PyDrizzle 297 for i in range(len(b_coeffs)): 298 b_coeffs[i] = float(b_coeffs[i]) 299 300 ifile.close() 301 del ifile 302 # Now, convert the coefficients into a Numeric array 303 # with the right coefficients in the right place. 304 # Populate output values now... 305 fx = N.zeros(shape=(order+1,order+1),dtype=N.float64) 306 fy = N.zeros(shape=(order+1,order+1),dtype=N.float64) 307 # Assign the coefficients to their array positions 308 fx[0,0] = 0. 309 fx[1] = N.array([a_coeffs[2],a_coeffs[1],0.,0.],dtype=N.float64) 310 fx[2] = N.array([a_coeffs[5],a_coeffs[4],a_coeffs[3],0.],dtype=N.float64) 311 fx[3] = N.array([a_coeffs[9],a_coeffs[8],a_coeffs[7],a_coeffs[6]],dtype=N.float64) 312 fy[0,0] = 0. 313 fy[1] = N.array([b_coeffs[2],b_coeffs[1],0.,0.],dtype=N.float64) 314 fy[2] = N.array([b_coeffs[5],b_coeffs[4],b_coeffs[3],0.],dtype=N.float64) 315 fy[3] = N.array([b_coeffs[9],b_coeffs[8],b_coeffs[7],b_coeffs[6]],dtype=N.float64) 316 317 # Used in Pattern.computeOffsets() 318 refpix = {} 319 refpix['XREF'] = None 320 refpix['YREF'] = None 321 refpix['V2REF'] = x0 322 refpix['V3REF'] = y0 323 refpix['XDELTA'] = 0. 324 refpix['YDELTA'] = 0. 325 refpix['PSCALE'] = None 326 refpix['DEFAULT_SCALE'] = no 327 refpix['centered'] = yes 328 329 return fx,fy,refpix,order

330 331 # Function to compute the index of refraction for MgF2 at 332 # the specified wavelength for use with Trauger coefficients

333 -def _MgF2(lam):

334 _sig = pow((1.0e7/lam),2) 335 return N.sqrt(1.0 + 2.590355e10/(5.312993e10-_sig) + 336 4.4543708e9/(11.17083e9-_sig) + 4.0838897e5/(1.766361e5-_sig))

337 338 # 339 # Function to read Trauger ASCII file and return cubic coefficients 340 #

341 -def readTraugerTable(idcfile,wavelength):

342 343 # Return a default geometry model if no coefficients filename 344 # is given. This model will not distort the data in any way. 345 if idcfile == None: 346 return fileutil.defaultModel() 347 348 # Trauger coefficients only result in a cubic file... 349 order = 3 350 numco = 10 351 a_coeffs = [0] * numco 352 b_coeffs = [0] * numco 353 indx = _MgF2(wavelength) 354 355 ifile = open(idcfile,'r') 356 # Search for the first line of the coefficients 357 _line = fileutil.rAsciiLine(ifile) 358 while string.lower(_line[:7]) != 'trauger': 359 _line = fileutil.rAsciiLine(ifile) 360 # Read in each row of coefficients,split them into their values, 361 # and convert them into cubic coefficients based on 362 # index of refraction value for the given wavelength 363 # Build X coefficients from first 10 rows of Trauger coefficients 364 j = 0 365 while j < 20: 366 _line = fileutil.rAsciiLine(ifile) 367 if _line == '': continue 368 _lc = string.split(_line) 369 if j < 10: 370 a_coeffs[j] = float(_lc[0])+float(_lc[1])*(indx-1.5)+float(_lc[2])*(indx-1.5)**2 371 else: 372 b_coeffs[j-10] = float(_lc[0])+float(_lc[1])*(indx-1.5)+float(_lc[2])*(indx-1.5)**2 373 j = j + 1 374 375 ifile.close() 376 del ifile 377 378 # Now, convert the coefficients into a Numeric array 379 # with the right coefficients in the right place. 380 # Populate output values now... 381 fx = N.zeros(shape=(order+1,order+1),dtype=N.float64) 382 fy = N.zeros(shape=(order+1,order+1),dtype=N.float64) 383 # Assign the coefficients to their array positions 384 fx[0,0] = 0. 385 fx[1] = N.array([a_coeffs[2],a_coeffs[1],0.,0.],dtype=N.float64) 386 fx[2] = N.array([a_coeffs[5],a_coeffs[4],a_coeffs[3],0.],dtype=N.float64) 387 fx[3] = N.array([a_coeffs[9],a_coeffs[8],a_coeffs[7],a_coeffs[6]],dtype=N.float64) 388 fy[0,0] = 0. 389 fy[1] = N.array([b_coeffs[2],b_coeffs[1],0.,0.],dtype=N.float64) 390 fy[2] = N.array([b_coeffs[5],b_coeffs[4],b_coeffs[3],0.],dtype=N.float64) 391 fy[3] = N.array([b_coeffs[9],b_coeffs[8],b_coeffs[7],b_coeffs[6]],dtype=N.float64) 392 393 # Used in Pattern.computeOffsets() 394 refpix = {} 395 refpix['XREF'] = None 396 refpix['YREF'] = None 397 refpix['V2REF'] = None 398 refpix['V3REF'] = None 399 refpix['XDELTA'] = 0. 400 refpix['YDELTA'] = 0. 401 refpix['PSCALE'] = None 402 refpix['DEFAULT_SCALE'] = no 403 refpix['centered'] = yes 404 405 return fx,fy,refpix,order

406

407 -def rotateCubic(fxy,theta):

408 # This function transforms cubic coefficients so that 409 # they calculate pixel positions oriented by theta (the same 410 # orientation as the PC). 411 # Parameters: fxy - cubic-coefficients Numeric array 412 # theta - angle to rotate coefficients 413 # Returns new array with same order as 'fxy' 414 # 415 # Set up some simplifications 416 newf = fxy * 0. 417 cost = N.cos(DEGTORAD(theta)) 418 sint = N.sin(DEGTORAD(theta)) 419 cos2t = pow(cost,2) 420 sin2t = pow(sint,2) 421 cos3t = pow(cost,3) 422 sin3t = pow(sint,3) 423 424 # Now compute the new coefficients 425 newf[1][1] = fxy[1][1] * cost - fxy[1][0] * sint 426 newf[1][0] = fxy[1][1] * sint + fxy[1][0] * cost 427 428 newf[2][2] = fxy[2][2] * cos2t - fxy[2][1] * cost * sint + fxy[2][0] * sin2t 429 newf[2][1] = fxy[2][2] * 2 * cost * sint + fxy[2][1] * (cos2t - sin2t) + fxy[2][0] * 2 * sint * cost 430 newf[2][0] = fxy[2][2] * sin2t + fxy[2][1] * cost * sint + fxy[2][0] * cos2t 431 432 newf[3][3] = fxy[3][3] * cos3t - fxy[3][2] * sint * cos2t + fxy[3][1] * sin2t * cost - fxy[3][0] * sin3t 433 newf[3][2] = fxy[3][3] * 3. * cos2t * sint + fxy[3][2]* (cos3t - 2. * sin2t * cost) +fxy[3][1] * (sin3t + 2 * sint * cos2t) - fxy[3][0] * sin2t * cost 434 newf[3][1] = fxy[3][3] * 3. * cost * sin2t + fxy[3][2] *(2.*cos2t*sint - sin3t) + fxy[3][1] * (2 * sin2t * cost + cos3t) + fxy[3][0] * sint * cos2t 435 newf[3][0] = fxy[3][3] * sin3t + fxy[3][2] * sin2t * cost + fxy[3][1] * sint * cos2t + fxy[3][0] * cos3t 436 437 return newf

438 439

440 -def rotatePos(pos, theta,offset=None,scale=None):

441 442 if scale == None: 443 scale = 1. 444 445 if offset == None: 446 offset = N.array([0.,0.],dtype=N.float64) 447 mrot = buildRotMatrix(theta) 448 xr = ((pos[0] * mrot[0][0]) + (pos[1]*mrot[0][1]) )/ scale + offset[0] 449 yr = ((pos[0] * mrot[1][0]) + (pos[1]*mrot[1][1]) )/ scale + offset[1] 450 451 return xr,yr

452 453 # Function for determining the positions of the image 454 # corners after applying the geometry model. 455 # Returns a dictionary with the limits and size of the 456 # image.

457 -def getRange(members,ref_wcs,verbose=None):

458 xma,yma = [],[] 459 xmi,ymi = [],[] 460 461 # Compute corrected positions of each chip's common point 462 crpix = (ref_wcs.crpix1,ref_wcs.crpix2) 463 ref_rot = ref_wcs.orient 464 _rot = ref_wcs.orient - members[0].geometry.wcslin.orient 465 466 for member in members: 467 # Need to make sure this is populated with proper defaults 468 # for ALL exposures! 469 _model = member.geometry.model 470 _wcs = member.geometry.wcs 471 _wcslin = member.geometry.wcslin 472 _theta = _wcslin.orient - ref_rot 473 474 # Need to scale corner positions to final output scale 475 _scale =_wcslin.pscale/ ref_wcs.pscale 476 477 # Compute the corrected,scaled corner positions for each chip 478 #xyedge = member.calcNewEdges(pscale=ref_wcs.pscale) 479 xypos = member.geometry.calcNewCorners() * _scale 480 481 if _theta != 0.0: 482 #rotate coordinates to match output orientation 483 # Now, rotate new coord 484 _mrot = buildRotMatrix(_theta) 485 xypos = N.dot(xypos,_mrot) 486 487 _oxmax = N.maximum.reduce(xypos[:,0]) 488 _oymax = N.maximum.reduce(xypos[:,1]) 489 _oxmin = N.minimum.reduce(xypos[:,0]) 490 _oymin = N.minimum.reduce(xypos[:,1]) 491 492 # Update the corners attribute of the member with the 493 # positions of the computed, distortion-corrected corners 494 #member.corners['corrected'] = N.array([(_oxmin,_oymin),(_oxmin,_oymax),(_oxmax,_oymin),(_oxmax,_oymax)],dtype=N.float64) 495 member.corners['corrected'] = xypos 496 497 xma.append(_oxmax) 498 yma.append(_oymax) 499 xmi.append(_oxmin) 500 ymi.append(_oymin) 501 502 # Determine the full size of the metachip 503 xmax = N.maximum.reduce(xma) 504 ymax = N.maximum.reduce(yma) 505 ymin = N.minimum.reduce(ymi) 506 xmin = N.minimum.reduce(xmi) 507 508 # Compute offset from center that distortion correction shifts the image. 509 # This accounts for the fact that the output is no longer symmetric around 510 # the reference position... 511 # Scale by ratio of plate-scales so that DELTAs are always in input frame 512 # 513 _ratio = ref_wcs.pscale / _wcslin.pscale 514 515 nref = ( (xmin + xmax)*_ratio, (ymin + ymax)*_ratio ) 516 if _rot != 0.: 517 _mrot = buildRotMatrix(_rot) 518 nref = N.dot(nref,_mrot) 519 520 # Now, compute overall size based on range of pixels and offset from center. 521 #xsize = int(xmax - xmin + nref[0]) 522 #ysize = int(ymax - ymin + nref[1]) 523 # Add '2' to each dimension to allow for fractional pixels at the 524 # edge of the image. Also, 'drizzle' centers the output, so 525 # adding 2 only expands image by 1 row on each edge. 526 # An additional two is added to accomodate floating point errors in drizzle. 527 xsize = int(ceil(xmax)) - int(floor(xmin)) 528 ysize = int(ceil(ymax)) - int(floor(ymin)) 529 530 meta_range = {} 531 meta_range = {'xmin':xmin,'xmax':xmax,'ymin':ymin,'ymax':ymax,'nref':nref} 532 meta_range['xsize'] = xsize 533 meta_range['ysize'] = ysize 534 535 if verbose: 536 print 'Meta_WCS:' 537 print ' NREF :',nref 538 print ' X range :',xmin,xmax 539 print ' Y range :',ymin,ymax 540 print ' Computed Size: ',xsize,ysize 541 542 return meta_range

543

544 -def computeRange(corners):

545 """ Determine the range spanned by an array of pixel positions. """ 546 _xrange = (N.minimum.reduce(corners[:,0]),N.maximum.reduce(corners[:,0])) 547 _yrange = (N.minimum.reduce(corners[:,1]),N.maximum.reduce(corners[:,1])) 548 return _xrange,_yrange

549

550 -def convertWCS(inwcs,drizwcs):

551 """ Copy WCSObject WCS into Drizzle compatible array.""" 552 drizwcs[0] = inwcs.crpix1 553 drizwcs[1] = inwcs.crval1 554 drizwcs[2] = inwcs.crpix2 555 drizwcs[3] = inwcs.crval2 556 drizwcs[4] = inwcs.cd11 557 drizwcs[5] = inwcs.cd21 558 drizwcs[6] = inwcs.cd12 559 drizwcs[7] = inwcs.cd22 560 561 return drizwcs

562

563 -def updateWCS(drizwcs,inwcs):

564 """ Copy output WCS array from Drizzle into WCSObject.""" 565 inwcs.crpix1 = drizwcs[0] 566 inwcs.crval1 = drizwcs[1] 567 inwcs.crpix2 = drizwcs[2] 568 inwcs.crval2 = drizwcs[3] 569 inwcs.cd11 = drizwcs[4] 570 inwcs.cd21 = drizwcs[5] 571 inwcs.cd12 = drizwcs[6] 572 inwcs.cd22 = drizwcs[7] 573 inwcs.pscale = N.sqrt(N.power(inwcs.cd11,2)+N.power(inwcs.cd21,2)) * 3600. 574 inwcs.orient = N.arctan2(inwcs.cd12,inwcs.cd22) * 180./N.pi

575

576 -def wcsfit(img_geom, ref):

577 """ 578 Perform a linear fit between 2 WCS for shift, rotation and scale. 579 Based on 'WCSLIN' from 'drutil.f'(Drizzle V2.9) and modified to 580 allow for differences in reference positions assumed by PyDrizzle's 581 distortion model and the coeffs used by 'drizzle'. 582 583 Parameters: 584 img - ObsGeometry instance for input image 585 ref_wcs - Undistorted WCSObject instance for output frame 586 """ 587 # Define objects that we need to use for the fit... 588 img_wcs = img_geom.wcs 589 in_refpix = img_geom.model.refpix 590 591 # Only work on a copy to avoid unexpected behavior in the 592 # call routine... 593 ref_wcs = ref.copy() 594 595 # Convert the RA/Dec positions back to X/Y in output product image 596 _cpix_xyref = N.zeros((4,2),dtype=N.float64) 597 598 # Start by setting up an array of points +/-0.5 pixels around CRVAL1,2 599 # However, we must shift these positions by 1.0pix to match what 600 # drizzle will use as its reference position for 'align=center'. 601 _cpix = (img_wcs.crpix1,img_wcs.crpix2) 602 _cpix_arr = N.array([_cpix,(_cpix[0],_cpix[1]+1.), 603 (_cpix[0]+1.,_cpix[1]+1.),(_cpix[0]+1.,_cpix[1])], dtype=N.float64) 604 # Convert these positions to RA/Dec 605 _cpix_rd = img_wcs.xy2rd(_cpix_arr) 606 for pix in xrange(len(_cpix_rd[0])): 607 _cpix_xyref[pix,0],_cpix_xyref[pix,1] = ref_wcs.rd2xy((_cpix_rd[0][pix],_cpix_rd[1][pix])) 608 609 # needed to handle correctly subarrays and wfpc2 data 610 if img_wcs.delta_refx == 0.0 and img_wcs.delta_refy == 0.0: 611 offx, offy = (0.0,0.0) 612 else: 613 offx, offy = (1.0, 1.0) 614 615 # Now, apply distortion model to input image XY positions 616 _cpix_xyc = N.zeros((4,2),dtype=N.float64) 617 _cpix_xyc[:,0],_cpix_xyc[:,1] = img_geom.apply(_cpix_arr - (offx, offy), order=1) 618 619 if in_refpix: 620 _cpix_xyc += (in_refpix['XDELTA'], in_refpix['YDELTA']) 621 622 # Perform a fit between: 623 # - undistorted, input positions: _cpix_xyc 624 # - X/Y positions in reference frame: _cpix_xyref 625 abxt,cdyt = fitlin(_cpix_xyc,_cpix_xyref) 626 627 # This correction affects the final fit when you are fitting 628 # a WCS to itself (no distortion coeffs), so it needs to be 629 # taken out in the coeffs file by modifying the zero-point value. 630 # WJH 17-Mar-2005 631 abxt[2] -= ref_wcs.crpix1 + offx 632 cdyt[2] -= ref_wcs.crpix2 + offy 633 634 return abxt,cdyt

635 636

637 -def fitlin(imgarr,refarr):

638 """ Compute the least-squares fit between two arrays. 639 A Python translation of 'FITLIN' from 'drutil.f' (Drizzle V2.9). 640 """ 641 # Initialize variables 642 _mat = N.zeros((3,3),dtype=N.float64) 643 _xorg = imgarr[0][0] 644 _yorg = imgarr[0][1] 645 _xoorg = refarr[0][0] 646 _yoorg = refarr[0][1] 647 _sigxox = 0. 648 _sigxoy = 0. 649 _sigxo = 0. 650 _sigyox = 0. 651 _sigyoy = 0. 652 _sigyo = 0. 653 654 _npos = len(imgarr) 655 # Populate matrices 656 for i in xrange(_npos): 657 _mat[0][0] += N.power((imgarr[i][0] - _xorg),2) 658 _mat[0][1] += (imgarr[i][0] - _xorg) * (imgarr[i][1] - _yorg) 659 _mat[0][2] += (imgarr[i][0] - _xorg) 660 _mat[1][1] += N.power((imgarr[i][1] - _yorg),2) 661 _mat[1][2] += imgarr[i][1] - _yorg 662 663 _sigxox += (refarr[i][0] - _xoorg)*(imgarr[i][0] - _xorg) 664 _sigxoy += (refarr[i][0] - _xoorg)*(imgarr[i][1] - _yorg) 665 _sigxo += refarr[i][0] - _xoorg 666 _sigyox += (refarr[i][1] - _yoorg)*(imgarr[i][0] -_xorg) 667 _sigyoy += (refarr[i][1] - _yoorg)*(imgarr[i][1] - _yorg) 668 _sigyo += refarr[i][1] - _yoorg 669 670 _mat[2][2] = _npos 671 _mat[1][0] = _mat[0][1] 672 _mat[2][0] = _mat[0][2] 673 _mat[2][1] = _mat[1][2] 674 675 # Now invert this matrix 676 _mat = linalg.inv(_mat) 677 678 _a = _sigxox*_mat[0][0]+_sigxoy*_mat[0][1]+_sigxo*_mat[0][2] 679 _b = -1*(_sigxox*_mat[1][0]+_sigxoy*_mat[1][1]+_sigxo*_mat[1][2]) 680 #_x0 = _sigxox*_mat[2][0]+_sigxoy*_mat[2][1]+_sigxo*_mat[2][2] 681 682 _c = _sigyox*_mat[1][0]+_sigyoy*_mat[1][1]+_sigyo*_mat[1][2] 683 _d = _sigyox*_mat[0][0]+_sigyoy*_mat[0][1]+_sigyo*_mat[0][2] 684 #_y0 = _sigyox*_mat[2][0]+_sigyoy*_mat[2][1]+_sigyo*_mat[2][2] 685 686 _xt = _xoorg - _a*_xorg+_b*_yorg 687 _yt = _yoorg - _d*_xorg-_c*_yorg 688 689 return [_a,_b,_xt],[_c,_d,_yt]

690

691 -def getRotatedSize(corners,angle):

692 """ Determine the size of a rotated (meta)image.""" 693 # If there is no rotation, simply return original values 694 if angle == 0.: 695 _corners = corners 696 else: 697 # Find center 698 #_xr,_yr = computeRange(corners) 699 #_cen = ( ((_xr[1] - _xr[0])/2.)+_xr[0],((_yr[1]-_yr[0])/2.)+_yr[0]) 700 _rotm = buildRotMatrix(angle) 701 # Rotate about the center 702 #_corners = N.dot(corners - _cen,_rotm) 703 _corners = N.dot(corners,_rotm) 704 705 return computeRange(_corners)

706

Source Code for Module pydrizzle.drutil