pydrizzle

1 from pytools import numerixenv 2 numerixenv.check() 3 4 import string,os,types,sys 5 import shutil 6 import datetime 7 8 #from pyraf import iraf 9 #from pyraf.iraf import stsdas,analysis,dither 10 # Import drizzle/blot Python callable module 11 import arrdriz 12 13 # Import PyDrizzle utility modules 14 import buildmask, drutil 15 import outputimage, imtype 16 from exposure import Exposure 17 18 from pytools import fileutil, wcsutil 19 import numpy as N 20 import pyfits 21 22 #Add buildasn/updateasn to namespace for use by other programs 23 import buildasn, updateasn 24 import traits102 25 26 yes = True # 1 27 no = False # 0 28 29 # List of supported instruments/detectors 30 INSTRUMENT = ["ACS","WFPC2","STIS","NICMOS"] 31 DEXPTIME = 'EXPTIME' 32 DEFAULT_PARITY = [[1.0,0.0],[0.0,1.0]] 33 34 from drutil import DEFAULT_IDCDIR 35 from pytools.fileutil import RADTODEG, DEGTORAD 36 from math import * 37 38 39 # Version 40 __version__ = "6.1 (25-Jan-2007" 41 42 # For History of changes and updates, see 'History' 43

44 -def _ptime():

45 import time 46 47 # Format time values for keywords IRAF-TLM, and DATE 48 _ltime = time.localtime(time.time()) 49 tlm_str = time.strftime('%H:%M:%S (%d/%m/%Y)',_ltime) 50 #date_str = time.strftime('%Y-%m-%dT%H:%M:%S',_ltime) 51 return tlm_str

52 53 ################# 54 # 55 # 56 # Custom Parameter Dictionary Class 57 # 58 # 59 ################# 60

61 -class ParDict(dict):

62 - def __init__(self):

63 dict.__init__(self)

64

65 - def __str__(self):

66 if self['xsh'] == None: _xsh = repr(self['xsh']) 67 else: _xsh = '%0.4f'%self['xsh'] 68 if self['ysh'] == None: _ysh = repr(self['ysh']) 69 else: _ysh = '%0.4f'%self['ysh'] 70 71 # Insure that valid strings are available for 72 # 'driz_mask' and 'single_driz_mask' 73 if self['driz_mask'] == None: _driz_mask = '' 74 else: _driz_mask = self['driz_mask'] 75 if self['single_driz_mask'] == None: _sdriz_mask = '' 76 else: _sdriz_mask = self['single_driz_mask'] 77 78 79 _str = '' 80 _str += 'Parameters for input chip: '+self['data']+'\n' 81 _str += ' Shifts: %s %s\n'%(_xsh,_ysh) 82 _str += ' Output size: %d %d \n'%(self['outnx'],self['outny']) 83 _str += ' Rotation: %0.4g deg. Scale: %0.4f \n'%(self['rot'],self['scale']) 84 _str += ' pixfrac: %0.4f Kernel: %s Units: %s\n'%(self['pixfrac'],self['kernel'],self['units']) 85 _str += ' ORIENTAT: %0.4g deg. Pixel Scale: %0.4f arcsec/pix\n'%(self['orient'],self['outscl']) 86 _str += ' Center at RA: %0.9g Dec: %0.9f \n'%(self['racen'],self['deccen']) 87 _str += ' Output product: %s\n'%self['output'] 88 _str += ' Coeffs: '+self['coeffs']+' Geomode: '+self['geomode']+'\n' 89 _str += ' XGeoImage : '+self['xgeoim']+'\n' 90 _str += ' YGeoImage : '+self['ygeoim']+'\n' 91 _str += ' Single mask file: '+_sdriz_mask+'\n' 92 _str += ' Final mask file : '+_driz_mask+'\n' 93 _str += ' Output science : '+self['outdata']+'\n' 94 _str += ' Output weight : '+self['outweight']+'\n' 95 _str += ' Output context : '+self['outcontext']+'\n' 96 _str += ' Exptime -- total: %0.4f single: %0.4f\n'%(self['texptime'],self['exptime']) 97 _str += ' start: %s end: %s\n'%(repr(self['expstart']),repr(self['expend'])) 98 _str += ' Single image products-- output: %s\n'%self['outsingle'] 99 _str += ' weight: %s \n '%self['outsweight'] 100 _str += ' Blot output: %s \n'%self['outblot'] 101 _str += ' Size of original image to blot: %d %d \n'%(self['blotnx'],self['blotny']) 102 103 _str += '\n' 104 105 return _str

106

107 -class Pattern:

108 """ 109 Set default values for these to be overridden by 110 instrument specific class variables as necessary. 111 """ 112 IDCKEY = 'IDCTAB' 113 PARITY = {'detector':DEFAULT_PARITY} 114 REFDATA = {'detector':[[1.,1.],[1.,1.]]} 115 116 NUM_IMSET = 3 # Number of extensions in an IMSET 117 PA_KEY = 'PA_V3' 118 DETECTOR_NAME = 'detector' 119 COPY_SUFFIX = '.orig' # suffix to use for filename of copy 120

121 - def __init__(self, filename, output=None, pars=None):

122 # Set these up for use... 123 self.members = [] 124 self.pars = pars 125 self.output = output 126 self.name = filename 127 128 # Set default value, to be reset as needed by sub-class 129 self.nmembers = 1 130 self.nimages = 1 131 132 # Extract bit values to be used for this instrument 133 self.bitvalue = self.pars['bits'] 134 135 # Set IDCKEY, if specified by user... 136 if self.pars['idckey'] == '': 137 self.idckey = self.IDCKEY 138 else: 139 self.idckey = self.pars['idckey'] 140 self.idcdir = self.pars['idcdir'] 141 142 # Keyword which provides PA_V3 value; ie., orientation of 143 # telescope at center axis relative to North. 144 # Each instrument has their own pre-defined keyword; mostly, PA_V3 145 self.pa_key = self.PA_KEY 146 147 self.exptime = None 148 149 self.binned = 1 150 151 # These attributes are used for keeping track of the reference 152 # image used for computing the shifts, and the shifts computed 153 # for each observation, respectively. 154 self.refimage = None 155 self.offsets = None 156 self.v2com = None 157 self.v3com = None 158 self.alpha = 0 159 self.beta = 0 160 161 # Read in Primary Header to reduce the overhead of getting 162 # keyword values and setup PyFITS object 163 # as self.header and self.image_handle 164 # 165 image_handle = self.getHeaderHandle() 166 167 if self.pars['section'] != None: 168 # If a section was specified, check the length of the list 169 # for the number of groups specified... 170 self.nmembers = len(self.pars['section']) 171 172 # Determine type of input image and syntax needed to read data 173 self.imtype = imtype.Imtype(filename,handle=self.image_handle, 174 dqsuffix=self.pars['dqsuffix']) 175 176 # Keep file I/O localized to same method/routine 177 if image_handle: 178 image_handle.close() 179 180 if self.header and self.header.has_key(self.DETECTOR_NAME): 181 self.detector = self.header[self.DETECTOR_NAME] 182 else: 183 self.detector = 'detector'

184 185

186 - def getHeaderHandle(self):

187 """ Sets up the PyFITS image handle and Primary header 188 as self.image_handle and self.header. 189 190 When Pattern being used for output product, filename will be 191 set to None and this returns None for header and image_handle. 192 """ 193 194 _numsci = 0 195 if self.name: 196 _handle = fileutil.openImage(self.name,mode='readonly',memmap=self.pars['memmap']) 197 _fname,_extn = fileutil.parseFilename(self.name) 198 _hdr = _handle['PRIMARY'].header.copy() 199 # Count number of SCI extensions 200 for _fext in _handle: 201 if _fext.header.has_key('extname') and _fext.header['extname'] == 'SCI': 202 _numsci += 1 203 204 if _extn > 0: 205 # Append correct extension/chip/group header to PRIMARY... 206 for _card in fileutil.getExtn(_handle,_extn).header.ascardlist(): 207 _hdr.ascard.append(_card) 208 else: 209 # Default to None 210 _handle = None 211 _hdr = None 212 213 # Set attribute to point to these products 214 self.image_handle = None 215 self.header = _hdr 216 self.nmembers = _numsci 217 218 return _handle

219

220 - def closeHandle(self):

221 """ Closes image_handle. """ 222 if self.image_handle: 223 self.image_handle.close() 224 #print 'Closing file handle for image: ',self.name 225 self.image_handle = None

226 227

228 - def addMembers(self,filename):

229 """ Build rootname for each SCI extension, and 230 create the mask image from the DQ extension. 231 It would then append a new Exposure object to 'members' 232 list for each extension. 233 """ 234 235 self.detector = detector = str(self.header[self.DETECTOR_NAME]) 236 237 if self.pars['section'] == None: 238 self.pars['section'] = [None]*self.nmembers 239 # Build rootname here for each SCI extension... 240 for i in range(self.nmembers): 241 _sciname = self.imtype.makeSciName(i+1,section=self.pars['section'][i]) 242 _dqname = self.imtype.makeDQName(i+1) 243 _extname = self.imtype.dq_extname 244 245 # Build mask files based on input 'bits' parameter values 246 _masklist = [] 247 _masknames = [] 248 # 249 # If we have a valid bits value... 250 # Creat the name of the output mask file 251 _maskname = buildmask.buildMaskName(_dqname,i+1) 252 _masknames.append(_maskname) 253 # Create the actual mask file now... 254 outmask = buildmask.buildMaskImage(_dqname,self.bitvalue[0],_maskname,extname=_extname,extver=i+1) 255 _masklist.append(outmask) 256 257 # 258 # Check to see if a bits value was provided for single drizzling... 259 # Different bits value specified for single drizzle step 260 # create new filename for single_drizzle mask file 261 _maskname = _maskname.replace('final_mask','single_mask') 262 _masknames.append(_maskname) 263 264 # Create new mask file with different bit value. 265 outmask = buildmask.buildMaskImage(_dqname,self.bitvalue[1],_maskname,extname=_extname,extver=i+1) 266 # Add new name to list for single drizzle step 267 _masklist.append(outmask) 268 _masklist.append(_masknames) 269 270 self.members.append(Exposure(_sciname, idckey=self.idckey, dqname=_dqname, 271 mask=_masklist, pa_key=self.pa_key, parity=self.PARITY[detector], 272 idcdir=self.pars['idcdir'], group_indx = i+1, 273 handle=self.image_handle,extver=i+1,exptime=self.exptime[0], mt_wcs=self.pars['mt_wcs']))

274

275 - def applyAsnShifts(self):

276 """ Apply ASN Shifts to each member and the observations product. """ 277 _prod_geo = self.product.geometry 278 _geo = self.product.geometry.wcslin 279 _geo0 = self.product.geometry.wcslin.copy() 280 281 _crval0 = (_geo.crval1,_geo.crval2) 282 _rot0 = _geo.orient 283 self._applyShifts(_geo) 284 285 _dcrval = (_crval0[0] - _geo.crval1,_crval0[1] - _geo.crval2) 286 287 _delta_rot = _rot0 - _geo.orient 288 _crpixlin = _geo.rd2xy(_crval0) 289 _delta_crpix = (_geo.crpix1 - _crpixlin[0],_geo.crpix2 - _crpixlin[1]) 290 291 # 292 # Do NOT apply AsnShifts to product.geometry.wcs as it serves 293 # as the default WCS for computing the output frame, allowing 294 # the images to be shifted relative to original WCS. This will 295 # result in output image where inputs may both be shifted off-center 296 # but presumably those shifts correct errors in WCS. WJH 28-Apr-03 297 #self._applyShifts(self.product.geometry.wcs) 298 # 299 # Update each members WCS with the new shifts, to allow for 300 # accurate RA/Dec measurements using internal methods. 301 for member in self.members: 302 # Compute chip specific shifts now... 303 _mem_geo = member.geometry 304 _mem_wcs = member.geometry.wcs 305 _wcs_copy = _mem_wcs.copy() 306 _mem_orient = member.geometry.wcs.orient 307 308 # Remember secondary shifts from ASN table 309 member.geometry.gpar_xsh = _delta_crpix[0] 310 member.geometry.gpar_ysh = _delta_crpix[1] 311 member.geometry.gpar_rot = _delta_rot 312 313 # Implement 'WBACKWCS' to update input WCS with shift computed 314 # in the output/reference WCS frame. 315 # 316 # Compute the pixel position of the input reference image in 317 # reference WCS 318 _oxy1 = _mem_geo.wtraxy((_mem_wcs.crpix1, _mem_wcs.crpix2 ),_geo0) 319 320 # Compute new CRPIX values for original WCS model to use for 321 # removing change in orientation from this shift. WJH 322 _oxy1o = _mem_geo.wtraxy((_mem_wcs.crpix1, _mem_wcs.crpix2),_geo) 323 324 # Use 'xy2rd' to convert these pixel positions to RA/Dec 325 ra1,dec1 = _geo.xy2rd((_oxy1[0], _oxy1[1])) 326 327 # Convert these numbers into a new WCS 328 _mem_wcs.crval1 = ra1 329 _mem_wcs.crval2 = dec1 330 331 """ 332 Now update CD matrix itself 333 We will have a lever-arm effect in place given the shift in 334 reference position. 335 336 """ 337 # 338 # Compute change in orientation due to change in CRPIX. 339 # This delta gets introduced due to the use of the updated 340 # undistorted WCS with the distorted WCS's reference positions. 341 # It must be taken out in order to get proper alignment for 342 # observations with large delta shifts. WJH 5-May-2004 343 # 344 _ref_wcs = _mem_wcs.copy() 345 _orig_orient = _ref_wcs.get_orient() 346 _ref_wcs.crpix1 += _oxy1o[0] - _oxy1[0] 347 _ref_wcs.crpix2 += _oxy1o[1] - _oxy1[1] 348 _ref_wcs.recenter() 349 350 _new_orient = _ref_wcs.get_orient() 351 _delta_orient = _new_orient - _orig_orient 352 _final_orient = _new_orient - _delta_rot 353 354 _mem_wcs.rotateCD(_final_orient) 355 _mem_wcs.orient = _mem_orient - _delta_orient 356 357 """ 358 Finish updating CD matrix 359 """ 360 361 # Make sure updates get translated back to undistorted WCS's 362 _mem_geo.undistortWCS() 363 364 # Update Member corner positions 365 member.corners['corrected'] += _delta_crpix 366 367 # Now update product corners 368 self.getProductCorners()

369

370 - def _applyShifts(self,in_wcs):

371 """ This method updates each member's WCS to reflect any 372 offsets/corrections specified in the ASN table. 373 374 This method converts shifts given in output pixels 375 into the input frame by using a reference image's 376 WCS. This reference image must exist and have an 377 header with a valid WCS; specifically, one which can 378 be read using WCSObject. 379 380 """ 381 # Check to see if there are any offsets given for this member... 382 # pars['shift'] will be yes if absolute shifts are given, or 383 # if rotations and/or scale changes are provided with/without shifts 384 # 385 386 # First, are there any shifts at all.. 387 if not self.pars.has_key('xshift'): 388 return 389 390 # Next, are they non-zero... 391 if not self.pars['abshift'] and not self.pars['dshift']: 392 return 393 # 394 # We have shifts, so apply them... 395 # 396 # Directly apply the shifts provided by the user in the ASN 397 # table to the original CRVAL values. 398 # 399 # Start by converting pixel shifts to units of arcseconds 400 # in RA and Dec. 401 # 402 #in_wcs = self.product.geometry.wcs 403 _crval = (in_wcs.crval1,in_wcs.crval2) 404 _crpix = (in_wcs.crpix1,in_wcs.crpix2) 405 # 406 # Setup tuples for containing the final delta's in input units 407 #converted from the given shifts regardless of the frame or type. 408 # 409 _dcrval = (0.,0.) 410 _dcrpix = (0.,0.) 411 412 _refimage = False 413 # If we are working with shifts provided in output frame... 414 if self.pars['refimage'] != '' and self.pars['refimage'] != None: 415 # Extract the reference image's WCS for use in converting 416 # the shifts into input pixels 417 # Only necessary for 'pixel' shifts... 418 # 419 _out_wcs = wcsutil.WCSObject(self.pars['refimage']) 420 _out_wcs.recenter() 421 _refimage = True 422 423 """ 424 Each product now has 'refimage' and 'offsets' attributes 425 refimage = {'pix_shift':(),'ra_shift':(),'name':'','val':0.} 426 offsets = {'pixels':(),'arcseconds':()} 427 428 PyDrizzle measures ALL shifts relative to the center of the final 429 output frame. The user's shifts will be measured relative to a 430 reference image's reference point. The difference between the 431 two frames must be accounted for when comparing shifts measured 432 in the two frames. 433 434 This offset needs to be subtracted from the shift provided by 435 the user in order to apply it to the PyDrizzle shifts. 436 """ 437 _ra_offset = (self.offsets['arcseconds'][0] - self.refimage['ra_shift'][0], 438 self.offsets['arcseconds'][1] - self.refimage['ra_shift'][1]) 439 440 _rotmat = fileutil.buildRotMatrix(self.pars['delta_rot']) 441 442 if self.pars['abshift']: 443 # Run 'computeOffsets()' to determine default shifts 444 # Extract shifts into numarray object and shift to input pixels 445 # Compute delta between abshift and 'offsets' as _delta_x,_delta_y 446 _delta_x = self.pars['xshift'] 447 _delta_y = self.pars['yshift'] 448 449 # Insure that 'pixel' shifts are in 'input' units 450 if self.pars['shift_units'] == 'pixels': 451 if _refimage: 452 # Compute delta RA,Dec for output frame ref point 453 # using given pixel shifts 454 _dcrpix = N.dot((-_delta_x, -_delta_y),_rotmat) 455 456 # CRVAL of output reference frame user shifts were measured from 457 _refcrval = _out_wcs.xy2rd((_dcrpix[0]+_out_wcs.crpix1,_dcrpix[1]+_out_wcs.crpix2)) 458 # Total arcsecond shift for this image 459 _ra_shift = ( (_refcrval[0] - _out_wcs.crval1), 460 (_refcrval[1] - _out_wcs.crval2)) 461 462 else: 463 # Tested against WFPC2 CADC ASNs... 464 _dcrpix = (_crpix[0] - _delta_x,_crpix[1] - _delta_y) 465 466 _refcrval = in_wcs.xy2rd(_dcrpix) 467 # Total arcsecond shift for this image 468 _ra_shift = ((_refcrval[0] - in_wcs.crval1),(_refcrval[1] - in_wcs.crval2)) 469 else: 470 # 471 # Arcsecond shifts 472 # 473 # Directly apply shift in arcseconds of RA/Dec to CRVALs 474 _ra_shift = (_delta_x/3600.,_delta_y/3600.) 475 476 # Delta RA/Dec based on output pixel shift for this image 477 _dcrval = (_ra_shift[0] -_ra_offset[0], _ra_shift[1] -_ra_offset[1]) 478 479 else: 480 # 481 # Delta Shifts 482 # Work with deltas from ASN table directly here 483 _delta_x = self.pars['delta_x'] 484 _delta_y = self.pars['delta_y'] 485 486 if self.pars['shift_units'] == 'pixels': 487 if _refimage: 488 # Compute delta RA,Dec for output frame ref point 489 # using given pixel shifts 490 #_dcrpix = N.dot((_out_wcs.crpix1-_delta_x, _out_wcs.crpix2-_delta_y),_rotmat) 491 _dcrpix = N.dot((-_delta_x, -_delta_y),_rotmat) 492 493 _refcrval = _out_wcs.xy2rd((_dcrpix[0]+_out_wcs.crpix1,_dcrpix[1]+_out_wcs.crpix2)) 494 _dcrval = ((_out_wcs.crval1 - _refcrval[0] ), 495 (_out_wcs.crval2 - _refcrval[1] )) 496 else: 497 _dcrpix = (_crpix[0] - _delta_x,_crpix[1] - _delta_y) 498 _refcrval = in_wcs.xy2rd(_dcrpix) 499 _dcrval = (_crval[0] - _refcrval[0],_crval[1] - _refcrval[1]) 500 else: 501 # 502 # Arcsecond Shifts 503 # 504 _dcrval = (_delta_x/3600.,_delta_y/3600.) 505 506 # Now, apply computed delta shift, _dcrval, to input WCS 507 # and recenter WCS... 508 # Update current CRVAL with delta we have computed 509 in_wcs.crval1 -= _dcrval[0] 510 in_wcs.crval2 -= _dcrval[1] 511 512 # Also, update orientation and plate scale 513 _orient = None 514 if self.pars['delta_rot'] != 0.: 515 _orient = in_wcs.orient + self.pars['delta_rot'] 516 517 _scale = None 518 if self.pars['delta_scale'] != 0.: 519 _scale = in_wcs.pscale * self.pars['delta_scale'] 520 521 # Compute new CRVAL for current CRPIX position 522 in_wcs.crval1,in_wcs.crval2 = in_wcs.xy2rd(_crpix) 523 in_wcs.crpix1 = _crpix[0] 524 in_wcs.crpix2 = _crpix[1] 525 526 #if not _refimage: 527 # Update product WCS with the new values 528 if _orient != None or _scale != None: 529 in_wcs.updateWCS(orient=_orient,pixel_scale=_scale)

530

531 - def getProductCorners(self):

532 """ Compute the product's corner positions based on input exposure's 533 corner positions. 534 """ 535 _prodcorners = [] 536 for member in self.members: 537 _prodcorners += member.corners['corrected'].tolist() 538 539 self.product.corners['corrected'] = N.array(_prodcorners,dtype=N.float64)

540

541 - def buildProduct(self,filename,output):

542 """ 543 Create Exposure object for meta-chip product after applying 544 distortion model to input members. 545 """ 546 #output_wcs = self._buildMetachip(psize=_psize,orient=_rot) 547 output_wcs = self.buildMetachip() 548 549 # For each member, update WCSLIN to account for chip-to-chip 550 # offsets, [X/Y]DELTA in refpix. 551 552 553 # Update the final output size with the meta-chip size 554 self.size = (output_wcs.naxis1,output_wcs.naxis2) 555 556 # Need to compute shape for dither product here based on output_shapes 557 # and reference pixels for each observation. 558 self.product = Exposure(self.rootname,wcs=output_wcs, new=yes) 559 560 self.product.exptime = self.exptime 561 562 # Preserve the original WCS as WCSLIN... 563 self.product.geometry.wcslin = self.product.geometry.wcs.copy() 564 565 # Combine all corner positions for each input into product's corner 566 self.getProductCorners()

567

568 - def buildMetachip(self,update=yes):

569 """ Build up the new metashape based on the 570 corrected size and position for each Exposure. 571 (Pattern method) 572 """ 573 # Build an initial value for size of meta-chip 574 _geo_ref = self.members[0].geometry 575 _wcs_ref = _geo_ref.wcslin 576 _model_ref = _geo_ref.model 577 578 # Get pscale from linearized WCS 579 pscale = _wcs_ref.pscale 580 581 _shape = (_wcs_ref.naxis1,_wcs_ref.naxis2,pscale) 582 # Build new WCS for output metachip here 583 # It will be based on the undistorted version of members[0] WCS 584 meta_wcs = _wcs_ref.copy() 585 586 # Now, check to see if members have subarray offsets, but was 587 # taken as a full image... 588 if len(self.members) > 1 and _geo_ref.wcs.subarray == yes: 589 for member in self.members: 590 member.geometry.wcslin.subarray = no 591 member.geometry.model.refpix['centered'] = no 592 593 #Determine range of pixel values for corrected image 594 # Using verbose=yes will return additional info on range calculation 595 meta_range = drutil.getRange(self.members,meta_wcs,verbose=no) 596 597 # Update WCS based on new size 598 xsize = int(meta_range['xsize']) 599 ysize = int(meta_range['ysize']) 600 meta_wcs.naxis1 = xsize 601 meta_wcs.naxis2 = ysize 602 cen = ((xsize/2.),(ysize/2.)) 603 meta_wcs.recenter() 604 605 _nref = meta_range['nref'] 606 for member in self.members: 607 member.corners['corrected'] -= (_nref[0]/2.,_nref[1]/2.) 608 609 if update: 610 # Shifts position of CRPIX to reflect new size 611 # Instead of being centered on (0.,0.) like the original guess. 612 # CRVAL for this position remains the same, though, as it is still 613 # the same point in the sky/image. 614 # 615 # We need this in order to correctly handle WFPC2 data where 616 # the XDELTA/YDELTA values are computed relative to the image center 617 # already, so we don't need this correction. 618 if not _model_ref.refpix['centered']: 619 _nref = meta_range['nref'] 620 621 for member in self.members: 622 _refpix = member.geometry.model.refpix 623 # Update XDELTA,YDELTA (zero-point of coefficients) to adjust for 624 # uncentered output 625 _refpix['XDELTA'] -= _nref[0]/2. 626 _refpix['YDELTA'] -= _nref[1]/2. 627 # 628 # TROLL computation not needed, as this get corrected for both 629 # in 'recenter()' and in 'wcsfit'... 630 # WJH 19-Feb-2004 631 # 632 # Do a full fit between the input WCS and meta_wcs now... 633 # 634 # Rigorously compute the orientation changes from WCS 635 # information using algorithm provided by R. Hook from WDRIZZLE. 636 abxt,cdyt = drutil.wcsfit(self.members[0].geometry, meta_wcs) 637 #Compute the rotation between input and reference from fit coeffs. 638 _delta_rot = RADTODEG(N.arctan2(abxt[1],cdyt[0])) 639 _crpix = (meta_wcs.crpix1 + abxt[2], meta_wcs.crpix2 + cdyt[2]) 640 641 meta_wcs.crval1,meta_wcs.crval2 = meta_wcs.xy2rd(_crpix) 642 # Insure output WCS has exactly orthogonal CD matrix 643 #meta_wcs.rotateCD(meta_wcs.orient+_delta_rot) 644 meta_wcs.updateWCS(orient=meta_wcs.orient+_delta_rot) 645 646 return meta_wcs

647

648 - def transformMetachip(self,ref):

649 """ 650 This method transforms this Exposure's WCS to be consistent 651 with the provided reference WCS 'ref'. This method only 652 operates on the product MetaChip, with the original WCS 653 being preserved as 'wcslin'. 654 655 Primarily, this transformation involves scaling and rotating 656 the chip to match the reference frame values. Also, any specified 657 size for the output frame would replace the default rotated/scaled 658 size. All rotations would be about the center, and the reference 659 pixel position gets shifted to accomodate this rotation. 660 661 """ 662 # Start by getting product WCS 663 # Use values from 'wcslin' as they will keep track of the 664 # default situation 665 _in_wcs = self.product.geometry.wcslin 666 _out_wcs = self.product.geometry.wcs 667 668 # Make sure that wcslin has all the WCS information 669 if _in_wcs.rootname == 'New': 670 # Copy the original WCS data into 'wcslin' to preserve the values 671 _in_wcs = self.product.geometry.wcs.copy() 672 673 #_dcrpix = (_in_wcs.naxis1/2.- _in_wcs.crpix1,_in_wcs.naxis2/2.- _in_wcs.crpix2) 674 675 # Check to see if there is any rotation needed 676 if ref.orient != None: 677 _angle = _in_wcs.orient - ref.wcs.orient 678 _ref_orient = ref.wcs.orient 679 else: 680 _angle = 0. 681 _ref_orient = _in_wcs.orient 682 683 # Apply any pixel scale changes to delta_crpix and size of axes 684 if ref.psize != None: 685 _scale = _in_wcs.pscale / ref.wcs.pscale 686 _ref_pscale = ref.wcs.pscale 687 else: 688 _scale = 1.0 689 _ref_pscale = _in_wcs.pscale 690 691 if ref.wcs.naxis1 != 0 and ref.wcs.naxis2 != 0: 692 _naxis1 = ref.wcs.naxis1 693 _naxis2 = ref.wcs.naxis2 694 # Delta between same-scaled frames 695 _delta_cens = (ref.wcs.naxis1/_scale - _in_wcs.naxis1,ref.wcs.naxis2/_scale - _in_wcs.naxis2) 696 else: 697 _delta_cens = (0.,0.) 698 _naxis1 = _in_wcs.naxis1 699 _naxis2 = _in_wcs.naxis2 700 701 _delta_cen = (_naxis1 - _in_wcs.naxis1, _naxis2 - _in_wcs.naxis2) 702 703 # Rotate axes as necessary 704 if _angle != 0.: 705 # Rotate axes to find default rotated size and new center 706 _xrange,_yrange = drutil.getRotatedSize(self.product.corners['corrected'],_angle) 707 _range = [_xrange[1] - _xrange[0] + _delta_cens[0]/2.,_yrange[1] - _yrange[0] + _delta_cens[1]/2.] 708 #_dcrpix = ((_xrange[0] + _xrange[1])/2.,(_yrange[0] + _yrange[1])/2.) 709 else: 710 _range = [_naxis1,_naxis2] 711 712 # Now update CD matrix and reference position 713 _out_wcs.naxis1 = int(_range[0]) 714 _out_wcs.naxis2 = int(_range[1]) 715 716 _crpix = (_in_wcs.crpix1 + _delta_cen[0]/2., 717 _in_wcs.crpix2 + _delta_cen[1]/2.) 718 719 if _scale != 0.: 720 _out_wcs.updateWCS(orient=_ref_orient,pixel_scale=_ref_pscale) 721 _delta_crpix = (_naxis1 - _out_wcs.naxis1, _naxis2 - _out_wcs.naxis2) 722 else: 723 _out_wcs.rotateCD(_ref_orient) 724 _delta_crpix = (0.,0.) 725 726 _out_wcs.crpix1 = _crpix[0] - _delta_crpix[0]/2. 727 _out_wcs.crpix2 = _crpix[1] - _delta_crpix[1]/2. 728 729 _out_wcs.recenter() 730 """ 731 # Update the size and rotated position of reference pixel 732 _cen = (_out_wcs.naxis1/2.,_out_wcs.naxis2/2.) 733 734 _out_wcs.crval1,_out_wcs.crval2 = _out_wcs.xy2rd(_cen) 735 _out_wcs.crpix1 = _cen[0] 736 _out_wcs.crpix2 = _cen[1] 737 """

738 # This method would use information from the product class and exposure class 739 # to build the complete parameter list for running the 'drizzle' task

740 - def buildPars(self,ref=None):

741 """This method would build a list of parameters to run 'drizzle' 742 one a single input image. 743 The reference image info 'ref' will be passed as a SkyField object. 744 The default output reference frame will be passed as 'def_wcs' 745 for comparison to the user's selected object 'ref'. 746 """ 747 # pars contains the drizzle parameters for each input(in order): 748 # data,outdata,outnx,outny,scale,xsh,ysh,rot 749 parlist = [] 750 751 # Now define the default WCS for this product 752 def_wcs = self.product.geometry.wcslin.copy() 753 754 if ref != None: 755 756 # Extract the total exptime for this output object 757 if ref.exptime == None: 758 _texptime = self.exptime 759 else: 760 _texptime = ref.exptime 761 # 762 # Initialize the SkyField Object with the 763 # user settings. 764 _field = ref 765 766 # Transform self.product.geometry.wcs to match ref.wcs 767 #self.transformMetachip(_field) 768 769 _out_wcs = self.product.geometry.wcs.copy() 770 771 #if not ref.dither or ref.dither == None: 772 # Check to make sure we have a complete WCS 773 # if not, fill in using the product's default WCS 774 _field.mergeWCS(_out_wcs) 775 _field.wcs.rootname=def_wcs.rootname 776 777 self.product.geometry.wcslin = _out_wcs 778 self.product.geometry.wcs = _field.wcs.copy() 779 780 # Set reference for computing shifts to be transformed 781 # product's WCS 782 ref_wcs = _field.wcs 783 784 else: 785 # 786 # Single observation case... 787 # 788 # Define a default WCS 789 ref_wcs = self.product.geometry.wcslin.copy() 790 791 # Update product WCS to reflect default value again 792 self.product.geometry.wcs = ref_wcs.copy() 793 794 # Extract the total exposure time 795 _texptime = self.product.exptime 796 797 # Insure that reference WCS is properly centered in order to 798 # correctly fit to each input WCS 799 ref_wcs.recenter() 800 801 for member in self.members: 802 in_wcs = member.geometry.wcslin 803 in_wcs_orig = member.geometry.wcs 804 805 #_delta_rot = in_wcs.orient - ref_wcs.orient 806 #_scale = ref_wcs.pscale / in_wcs.pscale 807 808 # Rigorously compute the orientation changes from WCS 809 # information using algorithm provided by R. Hook from WDRIZZLE. 810 abxt,cdyt = drutil.wcsfit(member.geometry, ref_wcs) 811 812 # Compute the rotation between input and reference from fit coeffs. 813 _delta_roty = _delta_rot = RADTODEG(N.arctan2(abxt[1],cdyt[0])) 814 _delta_rotx = RADTODEG(N.arctan2(abxt[0],cdyt[1])) 815 # Compute scale from fit to allow WFPC2 (and similar) data to be handled correctly 816 #_scale = 1./((N.sqrt(abxt[0]**2 + abxt[1]**2)+N.sqrt(cdyt[0]**2+cdyt[1]**2))/2.) 817 _scale = 1./N.sqrt(abxt[0]**2 + abxt[1]**2) 818 819 # Correct for additional shifts from shiftfile now 820 #_delta_x = abxt[2] + member.geometry.gpar_xsh 821 # _delta_y = cdyt[2] + member.geometry.gpar_ysh 822 _delta_x = abxt[2] 823 _delta_y = cdyt[2] 824 825 # Start building parameter dictionary for this chip 826 parameters = ParDict() 827 parameters['data'] = member.name 828 parameters['output'] = self.output 829 parameters['exposure'] = member 830 parameters['group'] = member.group_indx 831 832 parameters['instrument'] = self.instrument 833 parameters['detector'] = self.detector 834 835 parameters['driz_mask'] = member.maskname 836 parameters['single_driz_mask'] = member.singlemaskname 837 838 # Setup parameters for special cases here... 839 parameters['outsingle'] = self.outsingle 840 parameters['outsweight'] = self.outsweight 841 parameters['outscontext'] = self.outscontext 842 parameters['outblot'] = member.outblot 843 parameters['blotnx'] = member.naxis1 844 parameters['blotny'] = member.naxis2 845 846 #Setup parameters for normal operations 847 parameters['outdata'] = self.outdata 848 parameters['outweight'] = self.outweight 849 parameters['outcontext'] = self.outcontext 850 851 parameters['outnx'] = ref_wcs.naxis1 852 parameters['outny'] = ref_wcs.naxis2 853 854 parameters['xsh'] = _delta_x 855 parameters['ysh'] = _delta_y 856 857 parameters['alpha'] = self.alpha 858 parameters['beta'] = self.beta 859 860 # Calculate any rotation relative to the orientation 861 # AFTER applying ONLY the distortion coefficients without 862 # applying any additional rotation... 863 parameters['rot'] = _delta_rot 864 865 # Keep track of both the exposure information and 866 # the combined product exposure time information. 867 # Single exposure information will be used for 'single=yes' 868 # header updates... 869 parameters['exptime'] = self.exptime[0] 870 parameters['expstart'] = self.exptime[1] 871 parameters['expend'] = self.exptime[2] 872 parameters['texptime'] = _texptime[0] 873 parameters['texpstart'] = _texptime[1] 874 parameters['texpend'] = _texptime[2] 875 # Need to revise how this gets computed... 876 # The pixel scale of the product corresponds to the 877 # desired output pixel scale, and the model pscale for 878 # the member represents the un-scaled pixel size for the input. 879 parameters['scale'] = _scale 880 881 # Parameters only used by 'wdrizzle' 882 parameters['geomode'] = 'wcs' 883 parameters['racen'] = ref_wcs.crval1 884 parameters['deccen'] = ref_wcs.crval2 885 parameters['orient'] = ref_wcs.orient 886 parameters['outscl'] = ref_wcs.pscale 887 parameters['gpar_xsh'] = member.geometry.gpar_xsh 888 parameters['gpar_ysh'] = member.geometry.gpar_ysh 889 parameters['gpar_rot'] = member.geometry.gpar_rot 890 891 # Insure that the product WCS applied to each exposure gets set 892 member.product_wcs = ref_wcs 893 # Set up the idcfile for use by 'drizzle' 894 member.writeCoeffs() 895 parameters['coeffs'] = member.coeffs 896 897 parameters['plam'] = member.plam 898 899 # Set up the distortion image names as parameters 900 parameters['xgeoim'] = member.xgeoim 901 parameters['ygeoim'] = member.ygeoim 902 903 # Now pass along the remainder of the user specified parameters 904 if self.pars['units'] != None: 905 parameters['units'] = self.pars['units'] 906 else: 907 parameters['units'] = 'cps' 908 909 if self.pars['in_units'] != None: 910 parameters['in_units'] = self.pars['in_units'] 911 else: 912 parameters['in_units'] = 'counts' 913 914 if self.pars['pixfrac'] != None: 915 parameters['pixfrac'] = self.pars['pixfrac'] 916 else: 917 parameters['pixfrac'] = 1.0 918 919 if self.pars['kernel'] != None: 920 parameters['kernel'] = self.pars['kernel'] 921 else: 922 parameters['kernel'] = 'square' 923 924 if self.pars['wt_scl'] != None: 925 parameters['wt_scl'] = self.pars['wt_scl'] 926 else: 927 parameters['wt_scl'] = 'exptime' 928 929 if self.pars['fillval'] != None: 930 parameters['fillval'] = str(self.pars['fillval']) 931 else: 932 parameters['fillval'] = 'INDEF' 933 934 # Parameters useful for header keywords 935 parameters['version'] = 'PyDrizzle Version '+__version__ 936 parameters['driz_version'] = '' 937 parameters['nimages'] = self.nimages 938 939 parlist.append(parameters) 940 941 # Now, combine them for a complete set of pars for all exposures 942 # in this pattern/observation. 943 # 944 return parlist

945

946 - def computeCubicCoeffs(self):

947 """ 948 Method for converting cubic and Trauger coefficients tables 949 into a usable form. It also replaces 'computeOffsets' for 950 those tables as well. 951 """ 952 # For each chip in the observation... 953 _pscale1 = None 954 for img in self.members: 955 _chip = img.chip 956 _detector = str(img.header[self.DETECTOR_NAME]) 957 # scale all chips to first chip plate scale... 958 if _pscale1 == None or img.chip == '1': 959 _pscale1 = self.REFDATA[_detector]['psize'] 960 _reftheta = self.REFDATA[_detector]['theta'] 961 962 _v2ref = 0. 963 _v3ref = 0. 964 _nmembers = 0 965 for img in self.members: 966 # ... get the model and type of coefficients table used... 967 _model = img.geometry.model 968 _ikey = img.geometry.ikey 969 _chip = img.chip 970 _detector = str(img.header[self.DETECTOR_NAME]) 971 _refdata = self.REFDATA[_detector] 972 973 # ... determine the plate scale and scaling factor between chips... 974 if img.chip == '1': 975 _pscale = _refdata['psize'] 976 _ratio = 1.0 977 else: 978 _pscale = _refdata['psize'] 979 _ratio = _refdata['psize'] / _pscale1 980 # Record the plate scale for each chip's model that was used 981 # to compute the coefficients, not the plate scale from the 982 # image header. 983 _model.pscale = _pscale 984 _model.refpix['PSCALE'] = _pscale 985 986 if _ikey == 'trauger': 987 _ratio = 1. 988 989 # V2REF/V3REF was not available in idc file, so we 990 # must use our own data... 991 _model.refpix['V2REF'] = _refdata['xoff'] 992 _model.refpix['V3REF'] = _refdata['yoff'] 993 994 else: 995 _model.refpix['V2REF'] = _refdata['xoff'] * _pscale 996 _model.refpix['V3REF'] = _refdata['yoff'] * _pscale 997 998 # Correct the coefficients for the differences in plate scales 999 _model.cx = _model.cx * N.array([_model.pscale/_ratio],dtype=N.float64) 1000 _model.cy = _model.cy * N.array([_model.pscale/_ratio],dtype=N.float64) 1001 _model.refpix['XREF'] = self.REFPIX['x'] 1002 _model.refpix['YREF'] = self.REFPIX['y'] 1003 1004 # Correct the offsets for the plate scales as well... 1005 _model.refpix['XDELTA'] = _refdata['xoff'] 1006 _model.refpix['YDELTA'] = _refdata['yoff'] 1007 1008 _model.refpix['centered'] = yes 1009 1010 if _ikey != 'trauger': 1011 _model.refpix['V2REF'] = _model.refpix['V2REF'] / _ratio 1012 _model.refpix['V3REF'] = _model.refpix['V3REF'] / _ratio 1013 _v2ref += _model.refpix['V2REF'] 1014 _v3ref += _model.refpix['V3REF'] 1015 _nmembers += 1

1016 1017

1018 - def computeOffsets(self,parity=None,refchip=None):

1019 """ 1020 This version of 'computeOffsets' calculates the zero-point 1021 shifts to be included in the distortion coefficients table 1022 used by 'drizzle'. 1023 It REQUIRES a parity matrix to convert from 1024 V2/V3 coordinates into detector image X/Y coordinates. This 1025 matrix will be specific to each detector. 1026 """ 1027 vref = [] 1028 1029 # Check to see if any chip-to-chip offsets need to be computed at all 1030 if len(self.members) == 1: 1031 refp = self.members[0].geometry.model.refpix 1032 refp['XDELTA'] = 0. 1033 refp['YDELTA'] = 0. 1034 #refp['centered'] = yes 1035 return 1036 1037 # Set up the parity matrix here for a SINGLE chip 1038 if parity == None: 1039 # Use class defined dictionary as default 1040 parity = self.PARITY 1041 # Get reference chip information 1042 ref_model=None 1043 for member in self.members: 1044 if not refchip or refchip == int(member.chip): 1045 ref_model = member.geometry.model 1046 ref_scale = ref_model.refpix['PSCALE'] 1047 ref_v2v3 = N.array([ref_model.refpix['V2REF'],ref_model.refpix['V3REF']]) 1048 ref_theta = ref_model.refpix['THETA'] 1049 if ref_theta == None: ref_theta = 0.0 1050 ref_pmat = N.dot(fileutil.buildRotMatrix(ref_theta), member.parity) 1051 ref_xy = (ref_model.refpix['XREF'],ref_model.refpix['YREF']) 1052 break 1053 1054 if not ref_model: 1055 ref_scale = 1.0 1056 ref_theta = 0.0 1057 ref_v2v3 = [0.,0.] 1058 ref_xy = [0.,0.] 1059 ref_pmat = N.array([[1.,0.],[0.,1.0]]) 1060 1061 # Now compute the offset for each chip 1062 # Compute position of each chip's common point relative 1063 # to the output chip's reference position. 1064 for member in self.members: 1065 in_model = member.geometry.model 1066 refp = in_model.refpix 1067 pscale = in_model.pscale 1068 memwcs = member.geometry.wcs 1069 1070 v2v3 = N.array([in_model.refpix['V2REF'],in_model.refpix['V3REF']]) 1071 scale = refp['PSCALE'] 1072 theta = refp['THETA'] 1073 if theta == None: theta = 0.0 1074 1075 chipcen = ( (memwcs.naxis1/2.) + memwcs.offset_x, 1076 (memwcs.naxis2/2.) + memwcs.offset_y) 1077 xypos = N.dot(ref_pmat,v2v3-ref_v2v3) / scale + ref_xy 1078 chiprot = fileutil.buildRotMatrix(theta - ref_theta) 1079 1080 offcen = ((refp['XREF'] - chipcen[0]), (refp['YREF'] - chipcen[1])) 1081 1082 # Update member's geometry model with computed 1083 # reference position... 1084 #refp['XDELTA'] = vref[i][0] - v2com + chip.geometry.delta_x 1085 #refp['YDELTA'] = vref[i][1] - v3com + chip.geometry.delta_y 1086 offset_xy = N.dot(chiprot,xypos-offcen)*scale/ref_scale 1087 refp['XDELTA'] = offset_xy[0] 1088 refp['YDELTA'] = offset_xy[1] 1089 1090 # Only set centered to yes for full exposures... 1091 if member.geometry.wcs.subarray != yes: 1092 refp['centered'] = no 1093 else: 1094 refp['centered'] = yes

1095 1096

1097 - def setNames(self,filename,output):

1098 """ 1099 Define standard name attibutes: 1100 outname - Default final output name 1101 outdata - Name for drizzle science output 1102 outsingle - Name for output for single image 1103 """ 1104 self.rootname = filename 1105 self.outname = output 1106 1107 # Define FITS output filenames for intermediate products 1108 # to be used when 'build=no' 1109 self.outdata = fileutil.buildNewRootname(output,extn='_sci.fits') 1110 self.outweight = fileutil.buildNewRootname(output,extn='_weight.fits') 1111 self.outcontext = fileutil.buildNewRootname(output,extn='_context.fits') 1112 1113 # Define output file names for separate output for each input 1114 self.outsingle = fileutil.buildNewRootname(filename,extn='_single_sci.fits') 1115 self.outsweight = fileutil.buildNewRootname(filename,extn='_single_wht.fits') 1116 self.outscontext = None

1117 #self.outscontext = fileutil.buildNewRootname(filename,extn='_single_ctx.fits') 1118 1119 ######## 1120 # 1121 # User interface methods 1122 # 1123 ########

1124 - def getWCS(self):

1125 return self.members[0].getWCS()

1126

1127 - def getMember(self,memname):

1128 """ Return the class instance for the member with name memname.""" 1129 member = None 1130 for mem in self.members: 1131 if mem.name == memname: 1132 member = mem 1133 return member

1134

1135 - def getMemberNames(self):

1136 """ Return the names of all members for this Class. 1137 Output: [{self.name:[list of member names]}] 1138 """ 1139 memlist = [] 1140 for member in self.members: 1141 memlist.append(member.name) 1142 return [{self.name:memlist}]

1143

1144 - def getExptime(self):

1145 _exptime = float(self.header['EXPTIME']) 1146 if _exptime == 0.: _exptime = 1.0 1147 1148 if self.header.has_key('EXPSTART'): 1149 _expstart = float(self.header['EXPSTART']) 1150 _expend = float(self.header['EXPEND']) 1151 else: 1152 _expstart = 0. 1153 _expend = _exptime 1154 1155 return (_exptime,_expstart,_expend)

1156

1157 - def DeltaXYtoOffset(self,delta):

1158 """ 1159 Converts provided delta(x,y) pixel offset into a 1160 delta(RA,Dec) offset in arcseconds. 1161 """ 1162 _wcs = self.product.getWCS() 1163 _geom = self.product.geometry 1164 1165 new_rd = _geom.XYtoSky((_wcs.crpix1 - delta[0],_wcs.crpix2 - delta[1])) 1166 delta_ra = (_wcs.crval1 - new_rd[0]) * 3600. 1167 delta_dec = (_wcs.crval2 - new_rd[1]) * 3600. 1168 1169 return (delta_ra,delta_dec)

1170

1171 -class GenericObservation(Pattern):

1172 """ 1173 This class defines an observation stored in a Simple FITS format; 1174 i.e., only a Primary header and image without extensions. 1175 """ 1176 1177 REFPIX = {'x':512.,'y':512.} 1178 DETECTOR_NAME = 'INSTRUME' 1179

1180 - def __init__(self, filename, output, pars=None):

1181 1182 # Now initialize Pattern with all member exposures... 1183 Pattern.__init__(self, filename, output=output, pars=pars) 1184 1185 # Determine the instrument... 1186 if self.header.has_key('INSTRUME'): 1187 _instrument = self.header['INSTRUME'] 1188 else: 1189 _instrument = self.DETECTOR_NAME 1190 self.instrument = _instrument 1191 1192 if self.header.has_key('crpix1'): 1193 self.REFPIX['x'] = self.header['crpix1'] 1194 self.REFPIX['y'] = self.header['crpix2'] 1195 else: 1196 self.REFPIX['x'] = self.header['naxis1'] /2. 1197 self.REFPIX['y'] = self.header['naxis2'] /2. 1198 1199 # build output rootnames here... 1200 self.setNames(filename,output) 1201 1202 # Set EXPTIME for exposure 1203 self.exptime = self.getExptime() 1204 1205 self.nmembers = 1 1206 1207 # Now, build list of members and initialize them 1208 self.addMembers(filename) 1209 1210 _ikey = self.members[0].geometry.ikey 1211 if _ikey != 'idctab' and _ikey != 'wcs' : 1212 # Correct distortion coefficients to match output pixel scale 1213 self.computeCubicCoeffs() 1214 else: 1215 self.computeOffsets() 1216 1217 # Set up the input members and create the product meta-chip 1218 self.buildProduct(filename, output)

1219

1220 -class ACSObservation(Pattern):

1221 """This class defines an observation with information specific 1222 to ACS WFC exposures, including knowledge of how to mosaic both 1223 chips.""" 1224 1225 # Define a class variable for the gap between the chips 1226 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]]} 1227

1228 - def __init__(self, filename, output, pars=None):

1229 1230 # Now, initialize Pattern with all member Exposures... 1231 Pattern.__init__(self, filename, output=output, pars=pars) 1232 1233 self.instrument = 'ACS' 1234 1235 # build required name attributes: 1236 # outname, output, outsingle 1237 self.setNames(filename,output) 1238 1239 # Set EXPTIME for exposure 1240 self.exptime = self.getExptime() 1241 1242 # Build up list of chips in observation 1243 self.addMembers(filename) 1244 1245 # Only need to worry about IDC tables for ACS... 1246 self.computeOffsets() 1247 1248 # Set up the input members and create the product meta-chip 1249 self.buildProduct(filename, output) 1250 1251 # Compute the time dependent distrotion skew terms 1252 1253 # default date of 2004.5 = 2004-7-1 1254 #if self.header['DETECTOR'] == 'WFC': 1255 #datedefault = datetime.datetime(2004,7,1) 1256 #year,month,day = self.header['date-obs'].split('-') 1257 #rdate = datetime.datetime(int(year),int(month),int(day)) 1258 #self.alpha = 0.095+0.090*((rdate-datedefault).days/365.25)/2.5 1259 #self.beta = -0.029-0.030*((rdate-datedefault).days/365.25)/2.5 1260 #else: 1261 self.alpha = 0 1262 self.beta = 0

1263 1264

1265 -class STISObservation(Pattern):

1266 """This class defines an observation with information specific 1267 to STIS exposures. 1268 """ 1269 1270 # Default coefficients table to use for this instrument 1271 IDCKEY = 'cubic' 1272 1273 __theta = 0.0 1274 __parity = fileutil.buildRotMatrix(__theta) * N.array([[-1.,1.],[-1.,1.]]) 1275 PARITY = {'CCD':__parity,'NUV-MAMA':__parity,'FUV-MAMA':__parity} 1276 1277 # The dictionaries 'REFDATA' and 'REFPIX' are required for use with 1278 # cubic and Trauger coefficients tables in 'computeCubicCoeffs'. 1279 # 1280 # This information provides the absolute relationship between the chips 1281 # Latest plate scales: 0.05071, 0.0246 1282 #REFDATA = {'CCD':{'psize':0.05,'xoff':0.0,'yoff':0.0,'v2':-213.999,'v3':-224.897,'theta':__theta}, 1283 # 'NUV-MAMA':{'psize':0.024,'xoff':0.0,'yoff':0.0,'v2':-213.999,'v3':-224.897,'theta':__theta}, 1284 # 'FUV-MAMA':{'psize':0.024,'xoff':0.0,'yoff':0.0,'v2':-213.999,'v3':-224.897,'theta':__theta}} 1285 #REFPIX = {'x':512.,'y':512.} 1286

1287 - def __init__(self, filename, output, pars=None):

1288 1289 # Now initialize Pattern with all member exposures... 1290 Pattern.__init__(self, filename, output=output, pars=pars) 1291 1292 self.instrument = 'STIS' 1293 self.__theta = 0.0 1294 self.REFDATA = {'CCD':{'psize':0.05,'xoff':0.0,'yoff':0.0,'v2':-213.999,'v3':-224.897,'theta':self.__theta}, 1295 'NUV-MAMA':{'psize':0.024,'xoff':0.0,'yoff':0.0,'v2':-213.999,'v3':-224.897,'theta':self.__theta}, 1296 'FUV-MAMA':{'psize':0.024,'xoff':0.0,'yoff':0.0,'v2':-213.999,'v3':-224.897,'theta':self.__theta}} 1297 self.REFPIX = {'x':512.,'y':512.} 1298 1299 # build output rootnames here... 1300 self.setNames(filename,output) 1301 1302 # Set EXPTIME for exposure 1303 self.exptime = self.getExptime() 1304 1305 # Now, build list of members and initialize them 1306 self.addMembers(filename) 1307 1308 if self.members[0].geometry.ikey != 'idctab': 1309 # Correct distortion coefficients to match output pixel scale 1310 self.computeCubicCoeffs() 1311 else: 1312 self.computeOffsets() 1313 1314 # Set up the input members and create the product meta-chip 1315 self.buildProduct(filename, output)

1316

1317 - def getExptime(self):

1318 1319 header = fileutil.getHeader(self.name+'[sci,1]') 1320 _exptime = float(header['EXPTIME']) 1321 if _exptime == 0.: _exptime = 1.0 1322 1323 if header.has_key('EXPSTART'): 1324 _expstart = float(header['EXPSTART']) 1325 _expend = float(header['EXPEND']) 1326 else: 1327 _expstart = 0. 1328 _expend = _exptime 1329 1330 return (_exptime,_expstart,_expend)

1331 1332 1333 1334

1335 -class NICMOSObservation(Pattern):

1336 """This class defines an observation with information specific 1337 to NICMOS exposures. 1338 """ 1339 1340 # Default coefficients table to use for this instrument 1341 IDCKEY = 'cubic' 1342 1343 DETECTOR_NAME = 'camera' 1344 NUM_IMSET = 5 1345 1346 __theta = 0.0 1347 __parity = fileutil.buildRotMatrix(__theta) * N.array([[-1.,1.],[-1.,1.]]) 1348 PARITY = {'1':__parity,'2':__parity,'3':__parity} 1349 1350 # The dictionaries 'REFDATA' and 'REFPIX' are required for use with 1351 # cubic and Trauger coefficients tables in 'computeCubicCoeffs'. 1352 # 1353 # This information provides the absolute relationship between the chips 1354 # Latest plate scales: 0.05071, 0.0246 1355 #REFDATA = {'1':{'psize':0.0432,'xoff':0.0,'yoff':0.0,'v2':-296.9228,'v3':290.1827,'theta':__theta}, 1356 # '2':{'psize':0.076,'xoff':0.0,'yoff':0.0,'v2':-319.9464,'v3':311.8579,'theta':__theta}, 1357 # '3':{'psize':0.0203758,'xoff':0.0,'yoff':0.0,'v2':-249.8170,'v3':235.2371,'theta':__theta}} 1358 #REFPIX = {'x':128.,'y':128.} 1359

1360 - def __init__(self, filename, output, pars=None):

1361 1362 # Now initialize Pattern with all member exposures... 1363 Pattern.__init__(self, filename, output=output, pars=pars) 1364 1365 self.instrument = 'NICMOS' 1366 self.__theta = 0.0 1367 self.REFDATA = {'1':{'psize':0.0432,'xoff':0.0,'yoff':0.0,'v2':-296.9228,'v3':290.1827,'theta':self.__theta}, 1368 '2':{'psize':0.076,'xoff':0.0,'yoff':0.0,'v2':-319.9464,'v3':311.8579,'theta':self.__theta}, 1369 '3':{'psize':0.0203758,'xoff':0.0,'yoff':0.0,'v2':-249.8170,'v3':235.2371,'theta':self.__theta}} 1370 self.REFPIX = {'x':128.,'y':128.} 1371 # build output rootnames here... 1372 self.setNames(filename,output) 1373 1374 # Set EXPTIME for exposure 1375 self.exptime = self.getExptime() 1376 1377 # Now, build list of members and initialize them 1378 self.addMembers(filename) 1379 1380 if self.members[0].geometry.ikey != 'idctab': 1381 # Correct distortion coefficients to match output pixel scale 1382 self.computeCubicCoeffs() 1383 else: 1384 self.computeOffsets() 1385 1386 # Set up the input members and create the product meta-chip 1387 self.buildProduct(filename, output)

1388 1389

1390 -class WFPCObservation(Pattern):

1391 """This class defines an observation with information specific 1392 to WFPC2 exposures, including knowledge of how to mosaic the 1393 chips.""" 1394 1395 # Default coefficients table to use for this instrument 1396 IDCKEY = 'idctab' 1397 1398 # This parity is the multiplication of PC1 rotation matrix with 1399 # a flip in X for output image. 1400 #__theta = 44.67 1401 __pmat = N.array([[-1.,0.],[0.,1.]]) 1402 __refchip = 3 1403 PARITY = {'1':__pmat,'2':__pmat,'3':__pmat,'4':__pmat,'WFPC':__pmat} 1404 1405 NUM_IMSET = 1 1406 1407 # The dictionaries 'REFDATA' and 'REFPIX' are required for use with 1408 # cubic and Trauger coefficients tables in 'computeCubicCoeffs'. 1409 # 1410 # This information provides the absolute relationship between the chips 1411 1412 #REFDATA ={'1':{'psize':0.04554,'xoff':354.356,'yoff':343.646,'v2':2.374,'v3':-30.268,'theta':224.8480}, 1413 # '2':{'psize':0.0996,'xoff':345.7481,'yoff':375.28818,'v2':-51.368,'v3':-5.698,'theta':314.3520}, 1414 # '3':{'psize':0.0996,'xoff':366.56876,'yoff':354.79435,'v2':0.064,'v3':48.692,'theta':44.67}, 1415 # '4':{'psize':0.0996,'xoff':355.85016,'yoff':351.29183,'v2':55.044,'v3':-6.098,'theta':135.2210}} 1416 #REFPIX = {'x':400.,'y':400.} 1417

1418 - def __init__(self, filename, output, pars=None):

1419 1420 # Now initialize Pattern with all member exposures... 1421 Pattern.__init__(self, filename, output=output, pars=pars) 1422 1423 self.instrument = 'WFPC2' 1424 self.REFDATA ={'1':{'psize':0.04554,'xoff':354.356,'yoff':343.646,'v2':2.374,'v3':-30.268,'theta':224.8480}, 1425 '2':{'psize':0.0996,'xoff':345.7481,'yoff':375.28818,'v2':-51.368,'v3':-5.698,'theta':314.3520}, 1426 '3':{'psize':0.0996,'xoff':366.56876,'yoff':354.79435,'v2':0.064,'v3':48.692,'theta':44.67}, 1427 '4':{'psize':0.0996,'xoff':355.85016,'yoff':351.29183,'v2':55.044,'v3':-6.098,'theta':135.2210}} 1428 1429 self.REFPIX = {'x':400.,'y':400.} 1430 gcount = None 1431 1432 # build output rootnames here... 1433 self.setNames(filename,output) 1434 1435 # Set EXPTIME for exposure 1436 self.exptime = self.getExptime() 1437 1438 _mode = fileutil.getKeyword(filename, 'MODE') 1439 if _mode == 'AREA': 1440 self.binned = 2 1441 if self.idckey == 'cubic': 1442 for l in self.REFPIX.keys(): self.REFPIX[l]= self.REFPIX[l] / self.binned 1443 for l in self.REFDATA.keys(): 1444 self.REFDATA[l]['psize'] = self.REFDATA[l]['psize'] * self.binned 1445 self.REFDATA[l]['xoff'] = self.REFDATA[l]['xoff'] / self.binned 1446 self.REFDATA[l]['yoff'] = self.REFDATA[l]['yoff'] / self.binned 1447 1448 # Now, build list of members and initialize them 1449 self.addMembers(filename) 1450 1451 if self.members[0].geometry.ikey != 'idctab': 1452 # Correct distortion coefficients to match output pixel scale 1453 self.computeCubicCoeffs() 1454 else: 1455 self.computeOffsets(refchip=self.__refchip) 1456 1457 # Determine desired orientation of product 1458 self.setOrient() 1459 1460 # Set up the input members and create the product meta-chip 1461 self.buildProduct(filename, output)

1462 1463 # Correct the crpix position of the metachip product 1464 # in order to account for 'align=center' convention. 1465 #self.product.geometry.wcs.crpix1 -= 1.0 1466 #self.product.geometry.wcs.crpix2 -= 1.0 1467 1468

1469 - def addMembers(self,filename):

1470 1471 # The PC chip defines the orientation of the metachip, so use 1472 # it for the PARITY as well. 1473 self.detector = 'WFPC' 1474 1475 _chip1_rot = None 1476 # Build rootname here for each SCI extension... 1477 if self.pars['section'] == None: 1478 self.pars['section'] = [None]*self.nmembers 1479 1480 for i in range(self.nmembers): 1481 _extname = self.imtype.makeSciName(i+1,section=self.pars['section'][i]) 1482 1483 _detnum = fileutil.getKeyword(_extname,self.DETECTOR_NAME) 1484 1485 # Start by looking for the corresponding WFPC2 'c1h' files 1486 _dqfile, _dqextn = self._findDQFile() 1487 1488 # Reset dqfile name in ImType class to point to new file 1489 self.imtype.dqfile = _dqfile 1490 self.imtype.dq_extn = _dqextn 1491 1492 # Build mask file for this member chip 1493 _dqname = self.imtype.makeDQName(extver=_detnum) 1494 _masklist = [] 1495 _masknames = [] 1496 1497 if _dqname != None: 1498 _maskname = buildmask.buildMaskName(fileutil.buildNewRootname(_dqname),_detnum) 1499 else: 1500 _maskname = None 1501 _masknames.append(_maskname) 1502 1503 outmask = buildmask.buildShadowMaskImage(_dqname,_detnum,_maskname, bitvalue=self.bitvalue[0], binned=self.binned) 1504 _masklist.append(outmask) 1505 1506 _maskname = _maskname.replace('final_mask','single_mask') 1507 _masknames.append(_maskname) 1508 outmask = buildmask.buildShadowMaskImage(_dqname,_detnum,_maskname, bitvalue=self.bitvalue[1], binned=self.binned) 1509 _masklist.append(outmask) 1510 _masklist.append(_masknames) 1511 1512 1513 self.members.append(Exposure(_extname, idckey=self.idckey, dqname=_dqname, 1514 mask=_masklist, parity=self.PARITY[str(i+1)], 1515 idcdir=self.pars['idcdir'], group_indx = i+1, 1516 rot=_chip1_rot, handle=self.image_handle, extver=_detnum, 1517 exptime=self.exptime[0], ref_pscale=self.REFDATA['1']['psize'], binned=self.binned)) 1518 1519 if self.idckey != 'idctab': 1520 _chip1_rot = self.members[0].geometry.def_rot

1521

1522 - def _findDQFile(self):

1523 """ Find the DQ file which corresponds to the input WFPC2 image. """ 1524 if self.name.find('.fits') < 0: 1525 # Working with a GEIS image... 1526 dqfile = self.name[:-2]+'1h' 1527 dqextn = '[sdq,1]' 1528 else: 1529 # Looking for c1f FITS DQ file... 1530 # In the WFPC2 pipeline the extensions are called 'c0f' and 'c1f' 1531 # and EXTNAME is 'sci'. Files are in MEF format. 1532 # Readgeis creates output with extensions 'c0h' and 'c1h' 1533 # and EXPNAME is 'sdq' 1534 # Hence the code below ... 1535 if 'c0h.fits' in self.name: 1536 dqfile = self.name.replace('0h.fits','1h.fits') 1537 dqextn = '[sdq,1]' 1538 elif 'c0f.fits' in self.name: 1539 dqfile = self.name.replace('0f.fits','1f.fits') 1540 dqextn = '[sci,1]' 1541 1542 return dqfile, dqextn

1543

1544 - def setOrient(self):

1545 1546 """ Determine desired orientation of product.""" 1547 meta_orient = None 1548 for exp in self.members: 1549 if int(exp.chip) == 1: 1550 meta_orient = exp.geometry.wcslin.orient 1551 1552 if meta_orient == None: 1553 meta_orient = self.members[0].geometry.wcslin.orient 1554 1555 # Set orient for all groups 1556 # Dither coefficients rotate all chips to chip 1 orientation 1557 # if not using Trauger coefficients 1558 for exp in self.members: 1559 exp.geometry.wcs.orient = meta_orient

1560 # We need to correct each chip's full-frame reference 1561 # position to account for 'align=center' convention. 1562 #exp.geometry.wcs.chip_xref += 1.0 1563 #exp.geometry.wcs.chip_yref += 1.0 1564 1565

1566 -class DitherProduct(Pattern):

1567 """ 1568 Builds an object for a set of dithered inputs, each of which 1569 will be one of the Observation objects. 1570 """

1571 - def __init__(self, prodlist, pars=None):

1572 1573 # Build temporary output drizzle product name 1574 if prodlist['output'].find('.fits') < 0: 1575 if prodlist['output'].rfind('_drz') < 0: 1576 output = fileutil.buildNewRootname(prodlist['output'],extn='_drz.fits') 1577 else: 1578 output = prodlist['output']+'.fits' 1579 else: 1580 output = prodlist['output'] 1581 1582 # Setup a default exposure to contain the results 1583 Pattern.__init__(self, None, output=output, pars=pars) 1584 1585 self.pars = prodlist['members'] 1586 1587 # Subtract 2 from len(members) to account for 'abshift' and 'dshift' 1588 # entries in prodlist 1589 self.nmembers = self.nimages = len(prodlist['members']) - 2 1590 self.offsets = None 1591 1592 self.addMembers(prodlist,pars,output) 1593 1594 if len(self.members) == 0: 1595 print 'No suitable inputs from ASN table. Quitting PyDrizzle...' 1596 raise Exception 1597 1598 self.exptime = self.getExptime() 1599 1600 self.buildProduct(output)

1601

1602 - def closeHandle(self):

1603 """ Close image handle for each member.""" 1604 for member in self.members: 1605 member.closeHandle()

1606

1607 - def buildProduct(self,output):

1608 # Build default Metachip based on unmodified header WCS values 1609 output_wcs = self.buildMetachip() 1610 1611 self.size = (output_wcs.naxis1,output_wcs.naxis2) 1612 self.product = Exposure(output,wcs=output_wcs,new=yes) 1613 1614 # Compute default offsets between all inputs... 1615 # This will be used when applying relative shifts from ASN table 1616 # or shiftfile. 1617 self.computeOffsets() 1618 1619 # Apply any additional shifts from ASN table/shiftfile 1620 self.applyAsnShifts(output) 1621 1622 # Preserve default DitherProduct Metachip WCS as wcslin 1623 self.product.exptime = self.exptime 1624 # Update the corners arrays for the product now... 1625 self.product.setCorners() 1626 #self.product.corners['corrected'] = self.product.corners['raw'] 1627 _prodcorners = [] 1628 for prod in self.members: 1629 _prodcorners += prod.product.corners['corrected'].tolist() 1630 self.product.corners['corrected'] = N.array(_prodcorners,dtype=N.float64)

1631

1632

Source Code for Package pydrizzle