pydrizzle.makewcs

1 """ 2 MAKEWCS.PY - Updated the WCS in an image header so that 3 it matches the geometric distortion defined in an IDC table 4 which is referenced in the image header. 5 6 License: http://www.stsci.edu/resources/software_hardware/pyraf/LICENSE 7 8 This version tries to implement a full updating of the WCS based on 9 information about the V2/V3 plane which is obtained from th IDCTAB and, 10 in the case of WFPC2, the OFFTAB. 11 12 The only parameters from the original WCS which are retained are 13 the CRVALs of the reference chip. 14 15 The original WCS are first copied to MCD1_1 etc before being updated. 16 17 :UPINCD History: 18 First try, Richard Hook, ST-ECF/STScI, August 2002. 19 Version 0.0.1 (WJH) - Obtain IDCTAB using PyDrizzle function. 20 Version 0.1 (WJH) - Added support for processing image lists. 21 Revised to base CD matrix on ORIENTAT, instead of PA_V3 22 Supports subarrays by shifting coefficients as needed. 23 Version 0.2 (WJH) - Implemented orientation computation based on PA_V3 using 24 Troll function from Colin to compute new ORIENTAT value. 25 Version 0.3 (WJH) - Supported filter dependent distortion models in IDCTAB 26 fixed bugs in applying Troll function to WCS. 27 Version 0.4 (WJH) - Updated to support use of 'defaultModel' for generic 28 cases: XREF/YREF defaults to image center and idctab 29 name defaults to None. 30 Version 0.5 (WJH) - Added support for WFPC2 OFFTAB updates, which updates 31 the CRVALs. However, for WFPC2 data, the creation of 32 the backup values does not currently work. 33 :MAKEWCS History: 34 MAKEWCS V0.0 (RNH) - Created new version to implement more complete 35 WCS creation based on a reference tangent plane. 36 37 V0.1 (RNH) - First working version for tests. May 20th 2004. 38 V0.11 (RNH) - changed reference chip for ACS/WFC. May 26th 2004. 39 V0.2 (WJH) - Removed all dependencies from IRAF and use new WCSObject 40 class for all WCS operations. 41 V0.4 (WJH/CJH) - Corrected logic for looping of extension in FITS image. 42 V0.5 (RNH) - Chip to chip CRVAL shifting logic change. 43 V0.6 (CJH/WJH) - Added support for non-associated STIS data. 44 V0.6.2 (WJH) - Added support for NICMOS data. This required 45 new versions of wcsutil and fileutil in PyDrizzle. 46 V0.6.3 (WJH) - Modified to support new version of WCSUtil which correctly 47 sets up and uses archived WCS keywords. 48 V0.7.0 (WJH) - Revised algorithm to work properly with subarray images. 49 Also, simplified keyword access using PyFITS object. 50 V0.8.0 (CJH) - Modified to work with either numarray or numpy through 51 the use of the numerix interface layer. 52 53 """ 54 from __future__ import division # confidence high 55 56 from stsci.tools import numerixenv 57 numerixenv.check() 58 59 #import iraf 60 from math import * 61 import os.path 62 import pyfits 63 64 import drutil 65 from distortion import models, mutil 66 from stsci.tools import fileutil, wcsutil, parseinput 67 import numpy as N 68 69 yes = True 70 no = False 71 72 # Define parity matrices for supported detectors. 73 # These provide conversion from XY to V2/V3 coordinate systems. 74 # Ideally, this information could be included in IDCTAB... 75 PARITY = {'WFC':[[1.0,0.0],[0.0,-1.0]],'HRC':[[-1.0,0.0],[0.0,1.0]], 76 'SBC':[[-1.0,0.0],[0.0,1.0]],'default':[[1.0,0.0],[0.0,1.0]], 77 'WFPC2':[[-1.0,0.],[0.,1.0]],'STIS':[[-1.0,0.],[0.,1.0]], 78 'NICMOS':[[-1.0,0.],[0.,1.0]], 'UVIS':[[-1.0,0.0],[0.0,1.0]], 79 'IR':[[-1.0,0.0],[0.0,1.0]] } 80 81 NUM_PER_EXTN = {'ACS':3,'WFPC2':1,'STIS':3,'NICMOS':5, 'WFC3':3} 82 83 __version__ = '1.1.7 (6 Jul 2010)'

84 -def run(input,quiet=yes,restore=no,prepend='O', tddcorr=True):

85 86 print "+ MAKEWCS Version %s" % __version__ 87 88 _prepend = prepend 89 90 files = parseinput.parseinput(input)[0] 91 newfiles = [] 92 if files == []: 93 print "No valid input files found.\n" 94 raise IOError 95 96 for image in files: 97 #find out what the input is 98 imgfits,imgtype = fileutil.isFits(image) 99 100 # Check for existence of waiver FITS input, and quit if found. 101 if imgfits and imgtype == 'waiver': 102 """ 103 errormsg = '\n\nPyDrizzle does not support waiver fits format.\n' 104 errormsg += 'Convert the input files to GEIS or multiextension FITS.\n\n' 105 raise ValueError, errormsg 106 """ 107 newfilename = fileutil.buildNewRootname(image, extn='_c0h.fits') 108 # Convert GEIS image to MEF file 109 newimage = fileutil.openImage(image,writefits=True,fitsname=newfilename,clobber=True) 110 del newimage 111 # Work with new file 112 image = newfilename 113 newfiles.append(image) 114 # If a GEIS image is provided as input, create a new MEF file with 115 # a name generated using 'buildFITSName()' and update that new MEF file. 116 if not imgfits: 117 # Create standardized name for MEF file 118 newfilename = fileutil.buildFITSName(image) 119 # Convert GEIS image to MEF file 120 newimage = fileutil.openImage(image,writefits=True,fitsname=newfilename,clobber=True) 121 del newimage 122 # Work with new file 123 image = newfilename 124 newfiles.append(image) 125 126 if not quiet: 127 print "Input files: ",files 128 129 # First get the name of the IDC table 130 #idctab = drutil.getIDCFile(_files[0][0],keyword='idctab')[0] 131 idctab = drutil.getIDCFile(image,keyword='idctab')[0] 132 _found = fileutil.findFile(idctab) 133 if idctab == None or idctab == '': 134 print '#\n No IDCTAB specified. No correction can be done for file %s.Quitting makewcs\n' %image 135 #raise ValueError 136 continue 137 elif not _found: 138 print '#\n IDCTAB: ',idctab,' could not be found. \n' 139 print 'WCS keywords for file %s will not be updated.\n' %image 140 #raise IOError 141 continue 142 143 _phdu = image + '[0]' 144 _instrument = fileutil.getKeyword(_phdu,keyword='INSTRUME') 145 if _instrument == 'WFPC2': 146 Nrefchip, Nrefext = getNrefchip(image) 147 else: 148 Nrefchip = None 149 Nrefext = None 150 if not NUM_PER_EXTN.has_key(_instrument): 151 152 raise ValueError("Instrument %s not supported yet. Exiting..." \ 153 %_instrument) 154 155 _detector = fileutil.getKeyword(_phdu, keyword='DETECTOR') 156 _nimsets = get_numsci(image) 157 158 for i in xrange(_nimsets): 159 if image.find('.fits') > 0: 160 _img = image+'[sci,'+repr(i+1)+']' 161 else: 162 _img = image+'['+repr(i+1)+']' 163 if not restore: 164 if not quiet: 165 print 'Updating image: ', _img 166 167 _update(_img,idctab, _nimsets, apply_tdd=False, 168 quiet=quiet,instrument=_instrument,prepend=_prepend, 169 nrchip=Nrefchip, nrext = Nrefext) 170 if _instrument == 'ACS' and _detector == 'WFC': 171 tddswitch = fileutil.getKeyword(_phdu,keyword='TDDCORR') 172 # This logic requires that TDDCORR be in the primary header 173 # and set to PERFORM in order to turn this on at all. It can 174 # be turned off by setting either tddcorr=False or setting 175 # the keyword to anything but PERFORM or by deleting the 176 # keyword altogether. PyDrizzle will rely simply on the 177 # values of alpha and beta as computed here to apply the 178 # correction to the coefficients. 179 if (tddcorr and tddswitch != 'OMIT'): 180 print 'Applying time-dependent distortion corrections...' 181 _update(_img,idctab, _nimsets, apply_tdd=True, \ 182 quiet=quiet,instrument=_instrument,prepend=_prepend, nrchip=Nrefchip, nrext = Nrefext) 183 else: 184 if not quiet: 185 print 'Restoring original WCS values for',_img 186 restoreCD(_img,_prepend) 187 188 #fimg = fileutil.openImage(image,mode='update') 189 #if fimg[0].header.has_key('TDDCORR') and fimg[0].header['TDDCORR'] == 'PERFORM': 190 # fimg[0].header['TDDCORR'] = 'COMPLETE' 191 #fimg.close() 192 193 if newfiles == []: 194 return files 195 else: 196 return newfiles

197

198 -def restoreCD(image,prepend):

199 200 _prepend = prepend 201 try: 202 _wcs = wcsutil.WCSObject(image) 203 _wcs.restoreWCS(prepend=_prepend) 204 del _wcs 205 except: 206 print 'ERROR: Could not restore WCS keywords for %s.'%image

207

208 -def _update(image,idctab,nimsets,apply_tdd=False, 209 quiet=None,instrument=None,prepend=None,nrchip=None, nrext=None):

210 211 tdd_xyref = {1: [2048, 3072], 2:[2048, 1024]} 212 _prepend = prepend 213 _dqname = None 214 # Make a copy of the header for keyword access 215 # This copy includes both Primary header and 216 # extension header 217 hdr = fileutil.getHeader(image) 218 219 # Try to get the instrument if we don't have it already 220 instrument = readKeyword(hdr,'INSTRUME') 221 222 binned = 1 223 # Read in any specified OFFTAB, if present (WFPC2) 224 offtab = readKeyword(hdr,'OFFTAB') 225 dateobs = readKeyword(hdr,'DATE-OBS') 226 if not quiet: 227 print "OFFTAB, DATE-OBS: ",offtab,dateobs 228 229 print "-Updating image ",image 230 231 if not quiet: 232 print "-Reading IDCTAB file ",idctab 233 234 # Get telescope orientation from image header 235 # If PA_V# is not present of header, try to get it from the spt file 236 pvt = readKeyword(hdr,'PA_V3') 237 if pvt == None: 238 sptfile = fileutil.buildNewRootname(image, extn='_spt.fits') 239 if os.path.exists(sptfile): 240 spthdr = fileutil.getHeader(sptfile) 241 pvt = readKeyword(spthdr,'PA_V3') 242 if pvt != None: 243 pvt = float(pvt) 244 else: 245 print 'PA_V3 keyword not found, WCS cannot be updated. Quitting ...' 246 raise ValueError 247 248 # Find out about instrument, detector & filters 249 detector = readKeyword(hdr,'DETECTOR') 250 251 Nrefchip=1 252 if instrument == 'WFPC2': 253 filter1 = readKeyword(hdr,'FILTNAM1') 254 filter2 = readKeyword(hdr,'FILTNAM2') 255 mode = readKeyword(hdr,'MODE') 256 if os.path.exists(fileutil.buildNewRootname(image, extn='_c1h.fits')): 257 _dqname = fileutil.buildNewRootname(image, extn='_c1h.fits') 258 dqhdr = pyfits.getheader(_dqname,1) 259 dqext = readKeyword(dqhdr, 'EXTNAME') 260 if mode == 'AREA': 261 binned = 2 262 Nrefchip=nrchip 263 elif instrument == 'NICMOS': 264 filter1 = readKeyword(hdr,'FILTER') 265 filter2 = None 266 elif instrument == 'WFC3': 267 filter1 = readKeyword(hdr,'FILTER') 268 filter2 = None 269 # use value of 'BINAXIS' keyword to set binning value for WFC3 data 270 binned = readKeyword(hdr,'BINAXIS1') 271 else: 272 filter1 = readKeyword(hdr,'FILTER1') 273 filter2 = readKeyword(hdr,'FILTER2') 274 275 if filter1 == None or filter1.strip() == '': filter1 = 'CLEAR' 276 else: filter1 = filter1.strip() 277 if filter2 == None or filter2.strip() == '': filter2 = 'CLEAR' 278 else: filter2 = filter2.strip() 279 280 if filter1.find('CLEAR') == 0: filter1 = 'CLEAR' 281 if filter2.find('CLEAR') == 0: filter2 = 'CLEAR' 282 283 # Set up parity matrix for chip 284 if instrument == 'WFPC2' or instrument =='STIS' or instrument == 'NICMOS': 285 parity = PARITY[instrument] 286 elif PARITY.has_key(detector): 287 parity = PARITY[detector] 288 else: 289 raise ValueError('Detector ',detector, 290 ' Not supported at this time. Exiting...') 291 292 # Get the VAFACTOR keyword if it exists, otherwise set to 1.0 293 # we also need the reference pointing position of the target 294 # as this is where 295 _va_key = readKeyword(hdr,'VAFACTOR') 296 if _va_key != None: 297 VA_fac = float(_va_key) 298 else: 299 VA_fac=1.0 300 301 if not quiet: 302 print 'VA factor: ',VA_fac 303 304 #ra_targ = float(readKeyword(hdr,'RA_TARG')) 305 #dec_targ = float(readKeyword(hdr,'DEC_TARG')) 306 307 # Get the chip number 308 _c = readKeyword(hdr,'CAMERA') 309 _s = readKeyword(hdr,'CCDCHIP') 310 _d = readKeyword(hdr,'DETECTOR') 311 if _c != None and str(_c).isdigit(): 312 chip = int(_c) 313 elif _s == None and _d == None: 314 chip = 1 315 else: 316 if _s: 317 chip = int(_s) 318 elif str(_d).isdigit(): 319 chip = int(_d) 320 else: 321 chip = 1 322 # For the ACS/WFC case the chip number doesn't match the image 323 # extension 324 nr = 1 325 if (instrument == 'ACS' and detector == 'WFC') or (instrument == 'WFC3' and detector == 'UVIS'): 326 if nimsets > 1: 327 Nrefchip = 2 328 else: 329 Nrefchip = chip 330 elif instrument == 'NICMOS': 331 Nrefchip = readKeyword(hdr,'CAMERA') 332 elif instrument == 'WFPC2': 333 nr = nrext 334 else: 335 if nimsets > 1: 336 nr = Nrefchip 337 338 if not quiet: 339 print "-PA_V3 : ",pvt," CHIP #",chip 340 341 342 # Extract the appropriate information from the IDCTAB 343 #fx,fy,refpix,order=fileutil.readIDCtab(idctab,chip=chip,direction='forward', 344 # filter1=filter1,filter2=filter2,offtab=offtab,date=dateobs) 345 idcmodel = models.IDCModel(idctab, 346 chip=chip, direction='forward', date=dateobs, 347 filter1=filter1, filter2=filter2, offtab=offtab, binned=binned, 348 tddcorr=apply_tdd) 349 fx = idcmodel.cx 350 fy = idcmodel.cy 351 refpix = idcmodel.refpix 352 order = idcmodel.norder 353 354 # Determine whether to perform time-dependent correction 355 # Construct matrices neded to correct the zero points for TDD 356 if apply_tdd: 357 #alpha,beta = mutil.compute_wfc_tdd_coeffs(dateobs,skew_coeffs) 358 alpha = refpix['TDDALPHA'] 359 beta = refpix['TDDBETA'] 360 tdd = N.array([[beta, alpha], [alpha, -beta]]) 361 mrotp = fileutil.buildRotMatrix(2.234529)/2048. 362 363 else: 364 alpha = 0.0 365 beta = 0.0 366 367 # Get the original image WCS 368 Old=wcsutil.WCSObject(image,prefix=_prepend) 369 370 # Reset the WCS keywords to original archived values. 371 Old.restore() 372 373 # 374 # Look for any subarray offset 375 # 376 ltv1,ltv2 = drutil.getLTVOffsets(image) 377 # 378 # If reference point is not centered on distortion model 379 # shift coefficients to be applied relative to observation 380 # reference position 381 # 382 offsetx = Old.crpix1 - ltv1 - refpix['XREF'] 383 offsety = Old.crpix2 - ltv2 - refpix['YREF'] 384 shiftx = refpix['XREF'] + ltv1 385 shifty = refpix['YREF'] + ltv2 386 if ltv1 != 0. or ltv2 != 0.: 387 ltvoffx = ltv1 + offsetx 388 ltvoffy = ltv2 + offsety 389 offshiftx = offsetx + shiftx 390 offshifty = offsety + shifty 391 else: 392 ltvoffx = 0. 393 ltvoffy = 0. 394 offshiftx = 0. 395 offshifty = 0. 396 397 if ltv1 != 0. or ltv2 != 0.: 398 fx,fy = idcmodel.shift(idcmodel.cx,idcmodel.cy,offsetx,offsety) 399 400 # Extract the appropriate information for reference chip 401 402 ridcmodel = models.IDCModel(idctab, 403 chip=Nrefchip, direction='forward', date=dateobs, 404 filter1=filter1, filter2=filter2, offtab=offtab, binned=binned, 405 tddcorr=apply_tdd) 406 rfx = ridcmodel.cx 407 rfy = ridcmodel.cy 408 rrefpix = ridcmodel.refpix 409 rorder = ridcmodel.norder 410 """ 411 rfx,rfy,rrefpix,rorder=mutil.readIDCtab(idctab,chip=Nrefchip, 412 direction='forward', filter1=filter1,filter2=filter2,offtab=offtab, 413 date=dateobs,tddcorr=apply_tdd) 414 """ 415 # Create the reference image name 416 rimage = image.split('[')[0]+"[sci,%d]" % nr 417 if not quiet: 418 print "Reference image: ",rimage 419 420 # Create the tangent plane WCS on which the images are defined 421 # This is close to that of the reference chip 422 R=wcsutil.WCSObject(rimage) 423 R.write_archive(rimage) 424 R.restore() 425 426 # Reacd in declination of target (for computing orientation at aperture) 427 # Note that this is from the reference image 428 #dec = float(fileutil.getKeyword(rimage,'CRVAL2')) 429 #crval1 = float(fileutil.getKeyword(rimage,'CRVAL1')) 430 #crval1 = float(R.crval1) 431 #crval2 = dec 432 dec = float(R.crval2) 433 434 # Get an approximate reference position on the sky 435 rref = (rrefpix['XREF']+ltvoffx, rrefpix['YREF']+ltvoffy) 436 437 crval1,crval2=R.xy2rd(rref) 438 439 if apply_tdd: 440 # Correct zero points for TDD 441 tddscale = (R.pscale/fx[1][1]) 442 rxy0 = N.array([[tdd_xyref[Nrefchip][0]-2048.],[ tdd_xyref[Nrefchip][1]-2048.]]) 443 xy0 = N.array([[tdd_xyref[chip][0]-2048.], [tdd_xyref[chip][1]-2048.]]) 444 rv23_corr = N.dot(mrotp,N.dot(tdd,rxy0))*tddscale 445 v23_corr = N.dot(mrotp,N.dot(tdd,xy0))*tddscale 446 else: 447 rv23_corr = N.array([[0],[0]]) 448 v23_corr = N.array([[0],[0]]) 449 450 # Convert the PA_V3 orientation to the orientation at the aperture 451 # This is for the reference chip only - we use this for the 452 # reference tangent plane definition 453 # It has the same orientation as the reference chip 454 v2ref = rrefpix['V2REF'] + rv23_corr[0][0]*0.05 455 v3ref = rrefpix['V3REF'] - rv23_corr[1][0]*0.05 456 v2 = refpix['V2REF'] + v23_corr[0][0]*0.05 457 v3 = refpix['V3REF'] - v23_corr[1][0] *0.05 458 459 pv = wcsutil.troll(pvt,dec,v2ref,v3ref) 460 461 # Add the chip rotation angle 462 if rrefpix['THETA']: 463 pv += rrefpix['THETA'] 464 465 466 # Set values for the rest of the reference WCS 467 R.crval1=crval1 468 R.crval2=crval2 469 R.crpix1=0.0 + offshiftx 470 R.crpix2=0.0 + offshifty 471 472 R_scale=rrefpix['PSCALE']/3600.0 473 R.cd11=parity[0][0] * cos(pv*pi/180.0)*R_scale 474 R.cd12=parity[0][0] * -sin(pv*pi/180.0)*R_scale 475 R.cd21=parity[1][1] * sin(pv*pi/180.0)*R_scale 476 R.cd22=parity[1][1] * cos(pv*pi/180.0)*R_scale 477 478 ##print R 479 R_cdmat = N.array([[R.cd11,R.cd12],[R.cd21,R.cd22]]) 480 481 if not quiet: 482 print " Reference Chip Scale (arcsec/pix): ",rrefpix['PSCALE'] 483 484 # Offset and angle in V2/V3 from reference chip to 485 # new chip(s) - converted to reference image pixels 486 487 off = sqrt((v2-v2ref)**2 + (v3-v3ref)**2)/(R_scale*3600.0) 488 489 # Here we must include the PARITY 490 if v3 == v3ref: 491 theta=0.0 492 else: 493 theta = atan2(parity[0][0]*(v2-v2ref),parity[1][1]*(v3-v3ref)) 494 495 if rrefpix['THETA']: theta += rrefpix['THETA']*pi/180.0 496 497 dX=(off*sin(theta)) + offshiftx 498 dY=(off*cos(theta)) + offshifty 499 500 # Check to see whether we are working with GEIS or FITS input 501 _fname,_iextn = fileutil.parseFilename(image) 502 503 if _fname.find('.fits') < 0: 504 # Input image is NOT a FITS file, so 505 # build a FITS name for it's copy. 506 _fitsname = fileutil.buildFITSName(_fname) 507 else: 508 _fitsname = None 509 # Create a new instance of a WCS 510 if _fitsname == None: 511 _new_name = image 512 else: 513 _new_name = _fitsname+'['+str(_iextn)+']' 514 515 #New=wcsutil.WCSObject(_new_name,new=yes) 516 New = Old.copy() 517 518 # Calculate new CRVALs and CRPIXs 519 New.crval1,New.crval2=R.xy2rd((dX,dY)) 520 New.crpix1=refpix['XREF'] + ltvoffx 521 New.crpix2=refpix['YREF'] + ltvoffy 522 523 # Account for subarray offset 524 # Angle of chip relative to chip 525 if refpix['THETA']: 526 dtheta = refpix['THETA'] - rrefpix['THETA'] 527 else: 528 dtheta = 0.0 529 530 # Create a small vector, in reference image pixel scale 531 # There is no parity effect here ??? 532 delXX=fx[1,1]/R_scale/3600. 533 delYX=fy[1,1]/R_scale/3600. 534 delXY=fx[1,0]/R_scale/3600. 535 delYY=fy[1,0]/R_scale/3600. 536 537 # Convert to radians 538 rr=dtheta*pi/180.0 539 540 # Rotate the vectors 541 dXX= cos(rr)*delXX - sin(rr)*delYX 542 dYX= sin(rr)*delXX + cos(rr)*delYX 543 544 dXY= cos(rr)*delXY - sin(rr)*delYY 545 dYY= sin(rr)*delXY + cos(rr)*delYY 546 547 # Transform to sky coordinates 548 a,b=R.xy2rd((dX+dXX,dY+dYX)) 549 c,d=R.xy2rd((dX+dXY,dY+dYY)) 550 551 # Calculate the new CDs and convert to degrees 552 New.cd11=diff_angles(a,New.crval1)*cos(New.crval2*pi/180.0) 553 New.cd12=diff_angles(c,New.crval1)*cos(New.crval2*pi/180.0) 554 New.cd21=diff_angles(b,New.crval2) 555 New.cd22=diff_angles(d,New.crval2) 556 557 # Apply the velocity aberration effect if applicable 558 if VA_fac != 1.0: 559 560 # First shift the CRVALs apart 561 # New.crval1 = ra_targ + VA_fac*(New.crval1 - ra_targ) 562 # New.crval2 = dec_targ + VA_fac*(New.crval2 - dec_targ) 563 # First shift the CRVALs apart 564 # This is now relative to the reference chip, not the 565 # target position. 566 New.crval1 = R.crval1 + VA_fac*diff_angles(New.crval1, R.crval1) 567 New.crval2 = R.crval2 + VA_fac*diff_angles(New.crval2, R.crval2) 568 569 # and scale the CDs 570 New.cd11 = New.cd11*VA_fac 571 New.cd12 = New.cd12*VA_fac 572 New.cd21 = New.cd21*VA_fac 573 New.cd22 = New.cd22*VA_fac 574 575 New_cdmat = N.array([[New.cd11,New.cd12],[New.cd21,New.cd22]]) 576 577 # Store new one 578 # archive=yes specifies to also write out archived WCS keywords 579 # overwrite=no specifies do not overwrite any pre-existing archived keywords 580 581 New.write(fitsname=_new_name,overwrite=no,quiet=quiet,archive=yes) 582 if _dqname: 583 _dq_iextn = _iextn.replace('sci', dqext.lower()) 584 _new_dqname = _dqname +'['+_dq_iextn+']' 585 dqwcs = wcsutil.WCSObject(_new_dqname) 586 dqwcs.write(fitsname=_new_dqname, wcs=New,overwrite=no,quiet=quiet, archive=yes) 587 588 """ Convert distortion coefficients into SIP style 589 values and write out to image (assumed to be FITS). 590 """ 591 #First the CD matrix: 592 f = refpix['PSCALE']/3600.0 593 a = fx[1,1]/3600.0 594 b = fx[1,0]/3600.0 595 c = fy[1,1]/3600.0 596 d = fy[1,0]/3600.0 597 det = (a*d - b*c)*refpix['PSCALE'] 598 599 # Write to header 600 fimg = fileutil.openImage(_new_name,mode='update') 601 _new_root,_nextn = fileutil.parseFilename(_new_name) 602 _new_extn = fileutil.getExtn(fimg,_nextn) 603 604 605 # Transform the higher-order coefficients 606 for n in range(order+1): 607 for m in range(order+1): 608 if n >= m and n>=2: 609 610 # Form SIP-style keyword names 611 Akey="A_%d_%d" % (m,n-m) 612 Bkey="B_%d_%d" % (m,n-m) 613 614 # Assign them values 615 Aval= f*(d*fx[n,m]-b*fy[n,m])/det 616 Bval= f*(a*fy[n,m]-c*fx[n,m])/det 617 618 _new_extn.header.update(Akey,Aval) 619 _new_extn.header.update(Bkey,Bval) 620 621 # Update the SIP flag keywords as well 622 #iraf.hedit(image,"CTYPE1","RA---TAN-SIP",verify=no,show=no) 623 #iraf.hedit(image,"CTYPE2","DEC--TAN-SIP",verify=no,show=no) 624 _new_extn.header.update("CTYPE1","RA---TAN-SIP") 625 _new_extn.header.update("CTYPE2","DEC--TAN-SIP") 626 627 # Finally we also need the order 628 #iraf.hedit(image,"A_ORDER","%d" % order,add=yes,verify=no,show=no) 629 #iraf.hedit(image,"B_ORDER","%d" % order,add=yes,verify=no,show=no) 630 _new_extn.header.update("A_ORDER",order) 631 _new_extn.header.update("B_ORDER",order) 632 633 # Update header with additional keywords required for proper 634 # interpretation of SIP coefficients by PyDrizzle. 635 636 _new_extn.header.update("IDCSCALE",refpix['PSCALE']) 637 _new_extn.header.update("IDCV2REF",refpix['V2REF']) 638 _new_extn.header.update("IDCV3REF",refpix['V3REF']) 639 _new_extn.header.update("IDCTHETA",refpix['THETA']) 640 _new_extn.header.update("OCX10",fx[1][0]) 641 _new_extn.header.update("OCX11",fx[1][1]) 642 _new_extn.header.update("OCY10",fy[1][0]) 643 _new_extn.header.update("OCY11",fy[1][1]) 644 #_new_extn.header.update("TDDXOFF",rv23_corr[0][0] - v23_corr[0][0]) 645 #_new_extn.header.update("TDDYOFF",-(rv23_corr[1][0] - v23_corr[1][0])) 646 647 # Report time-dependent coeffs, if computed 648 if instrument == 'ACS' and detector == 'WFC': 649 _new_extn.header.update("TDDALPHA",alpha) 650 _new_extn.header.update("TDDBETA",beta) 651 652 653 # Close image now 654 fimg.close() 655 del fimg

656 657

658 -def diff_angles(a,b):

659 """ Perform angle subtraction a-b taking into account 660 small-angle differences across 360degree line. """ 661 662 diff = a - b 663 664 if diff > 180.0: 665 diff -= 360.0 666 667 if diff < -180.0: 668 diff += 360.0 669 670 return diff

671

672 -def readKeyword(hdr,keyword):

673 674 try: 675 value = hdr[keyword] 676 except KeyError: 677 value = None 678 679 # NOTE: Need to clean up the keyword.. Occasionally the keyword value 680 # goes right up to the "/" FITS delimiter, and iraf.keypar is incapable 681 # of realizing this, so it incorporates "/" along with the keyword value. 682 # For example, after running "pydrizzle" on the image "j8e601bkq_flt.fits", 683 # the CD keywords look like this: 684 # 685 # CD1_1 = 9.221627430999639E-06/ partial of first axis coordinate w.r.t. x 686 # CD1_2 = -1.0346992614799E-05 / partial of first axis coordinate w.r.t. y 687 # 688 # so for CD1_1, iraf.keypar returns: 689 # "9.221627430999639E-06/" 690 # 691 # So, the following piece of code CHECKS for this and FIXES the string, 692 # very simply by removing the last character if it is a "/". 693 # This fix courtesy of Anton Koekemoer, 2002. 694 if isinstance(value, basestring): 695 if value[-1:] == '/': 696 value = value[:-1] 697 698 return value

699

700 -def get_numsci(image):

701 """ Find the number of SCI extensions in the image. 702 Input: 703 image - name of single input image 704 """ 705 handle = fileutil.openImage(image) 706 num_sci = 0 707 for extn in handle: 708 if extn.header.has_key('extname'): 709 if extn.header['extname'].lower() == 'sci': 710 num_sci += 1 711 handle.close() 712 return num_sci

713

714 -def shift_coeffs(cx,cy,xs,ys,norder):

715 """ 716 Shift reference position of coefficients to new center 717 where (xs,ys) = old-reference-position - subarray/image center. 718 This will support creating coeffs files for drizzle which will 719 be applied relative to the center of the image, rather than relative 720 to the reference position of the chip. 721 722 Derived directly from PyDrizzle V3.3d. 723 """ 724 725 _cxs = N.zeros(shape=cx.shape,dtype=cx.dtype.name) 726 _cys = N.zeros(shape=cy.shape,dtype=cy.dtype.name) 727 _k = norder + 1 728 729 # loop over each input coefficient 730 for m in xrange(_k): 731 for n in xrange(_k): 732 if m >= n: 733 # For this coefficient, shift by xs/ys. 734 _ilist = N.array(range(_k - m)) + m 735 # sum from m to k 736 for i in _ilist: 737 _jlist = N.array(range( i - (m-n) - n + 1)) + n 738 # sum from n to i-(m-n) 739 for j in _jlist: 740 _cxs[m,n] = _cxs[m,n] + cx[i,j]*pydrizzle._combin(j,n)*pydrizzle._combin((i-j),(m-n))*pow(xs,(j-n))*pow(ys,((i-j)-(m-n))) 741 _cys[m,n] = _cys[m,n] + cy[i,j]*pydrizzle._combin(j,n)*pydrizzle._combin((i-j),(m-n))*pow(xs,(j-n))*pow(ys,((i-j)-(m-n))) 742 _cxs[0,0] = _cxs[0,0] - xs 743 _cys[0,0] = _cys[0,0] - ys 744 #_cxs[0,0] = 0. 745 #_cys[0,0] = 0. 746 747 return _cxs,_cys

748

749 -def getNrefchip(image,instrument='WFPC2'):

750 """ 751 This handles the fact that WFPC2 subarray observations 752 may not include chip 3 which is the default reference chip for 753 full observations. Also for subarrays chip 3 may not be the third 754 extension in a MEF file. It is a kludge but this whole module is 755 one big kludge. ND 756 """ 757 hdu = fileutil.openImage(image) 758 if instrument == 'WFPC2': 759 detectors = [img.header['DETECTOR'] for img in hdu[1:]] 760 761 if 3 not in detectors: 762 Nrefchip=detectors[0] 763 Nrefext = 1 764 else: 765 Nrefchip = 3 766 Nrefext = detectors.index(3) + 1 767 768 hdu.close() 769 return Nrefchip, Nrefext

770 771 772 _help_str = """ makewcs - a task for updating an image header WCS to make 773 it consistent with the distortion model and velocity aberration. 774 775 This task will read in a distortion model from the IDCTAB and generate 776 a new WCS matrix based on the value of ORIENTAT. It will support subarrays 777 by shifting the distortion coefficients to image reference position before 778 applying them to create the new WCS, including velocity aberration. 779 Original WCS values will be moved to an O* keywords (OCD1_1,...). 780 Currently, this task will only support ACS and WFPC2 observations. 781 782 Parameters 783 ---------- 784 input: str 785 The filename(s) of image(s) to be updated given either as: 786 * a single image with extension specified, 787 * a substring common to all desired image names, 788 * a wildcarded filename 789 * '@file' where file is a file containing a list of images 790 791 quiet: bool 792 turns off ALL reporting messages: 'yes' or 'no'(default) 793 794 prepend: char 795 This parameter specifies what prefix letter should be used to 796 create a new set of WCS keywords for recording the original values 797 [Default: 'O'] 798 799 restore: bool 800 restore WCS for all input images to defaults if possible: 801 'yes' or 'no'(default) 802 803 tddcorr: bool 804 applies the time-dependent skew terms to the SIP coefficients 805 written out to the header: 'yes' or True or, 'no' or False (default). 806 807 Notes 808 ----- 809 This function can be run using the syntax: 810 makewcs.run(image,quiet=no,prepend='O',restore=no,tddcorr=True) 811 An example of how this can be used is given as:: 812 813 >>> import makewcs 814 >>> makewcs.run('raw') # This will update all _raw files in directory 815 >>> makewcs.run('j8gl03igq_raw.fits[sci,1]') 816 """ 817

818 -def help():

819 print _help_str

820 821 run.__doc__ = _help_str 822

Source Code for Module pydrizzle.makewcs