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 - def applyAsnShifts(self,output):

1633 """ 1634 Apply any shifts read in from the ASN table /shiftfile to 1635 the WCS of each input's product. 1636 1637 In the case there are no shifts to apply, copy the product 1638 WCS into WCSLIN as the default case. 1639 1640 This method should (eventually?) support updating the shifts 1641 without starting from the beginning with these datasets, similar 1642 to the 'resetPars' functionality. 1643 1644 """ 1645 # If there were any shifts to be applied, update input 1646 # image's WCS with the shifts, then re-build the final 1647 # product Metachip using corrected product WCS values. 1648 if self.pars['abshift'] or self.pars['dshift']: 1649 for prod in self.members: 1650 prod.applyAsnShifts() 1651 1652 #If we have shifts of any sort, 1653 # recompute the DitherProduct's meta-chip with corrected values 1654 output_wcs = self.buildMetachip() 1655 self.size = (output_wcs.naxis1,output_wcs.naxis2) 1656 self.product = Exposure(output,wcs=output_wcs,new=yes) 1657 1658 self.product.geometry.wcslin = self.product.geometry.wcs.copy()

1659

1660 - def computeOffsets(self):

1661 """ 1662 This method will rely on final product's 'rd2xy' method 1663 to compute offsets between the different input chips. 1664 """ 1665 ref_wcs = self.product.geometry.wcs 1666 ref_crv= (ref_wcs.crval1,ref_wcs.crval2) 1667 1668 _tot_ref = {'pix_shift':(0.,0.),'ra_shift':(0.,0.), 1669 'val':999999999.,'name':''} 1670 1671 _member_num = 0 1672 for member in self.members: 1673 in_wcs = member.product.geometry.wcs 1674 1675 x,y = ref_wcs.rd2xy((in_wcs.crval1,in_wcs.crval2)) 1676 xoff = x - ref_wcs.crpix1 1677 yoff = y - ref_wcs.crpix2 1678 1679 raoff = (in_wcs.crval1 - ref_wcs.crval1,in_wcs.crval2 - ref_wcs.crval2) 1680 1681 _delta_rot = in_wcs.orient - ref_wcs.orient 1682 # Determine which image has the smallest offset 1683 _tot = N.sqrt(N.power(xoff,2)+N.power(yoff,2)) 1684 1685 # Use first image as reference 1686 if _member_num == 0: 1687 _tot_ref['val'] = _tot 1688 _tot_ref['pix_shift'] = (xoff,yoff) 1689 _tot_ref['name'] = member.rootname 1690 _tot_ref['ra_shift'] = raoff 1691 # This will be used primarily for manual verification 1692 _tot_ref['wcs'] = in_wcs 1693 1694 _member_num += 1 1695 1696 # Keep track of the results for later use/comparison 1697 member.offsets = {'pixels':(xoff,yoff),'arcseconds':raoff} 1698 1699 for member in self.members: 1700 member.refimage = _tot_ref

1701 1702

1703 - def getExptime(self):

1704 """ 1705 Add up the exposure time for all the members in 1706 the pattern, since 'drizzle' doesn't have the necessary 1707 information to correctly set this itself. 1708 """ 1709 _exptime = 0. 1710 _start = [] 1711 _end = [] 1712 for member in self.members: 1713 _memexp = member.product.exptime 1714 _exptime = _exptime + _memexp[0] 1715 _start.append(_memexp[1]) 1716 _end.append(_memexp[2]) 1717 _expstart = min(_start) 1718 _expend = max(_end) 1719 1720 return (_exptime,_expstart,_expend)

1721

1722 - def addMembers(self,prodlist,pars,output):

1723 """ 1724 For each entry in prodlist, append the appropriate type 1725 of Observation to the members list. 1726 1727 If it's a moving target observation 1728 apply the wcs of the first observation to 1729 all other observations. 1730 """ 1731 1732 member = prodlist['order'][0] 1733 filename = fileutil.buildRootname(member) 1734 mtflag = fileutil.getKeyword(filename, 'MTFLAG') 1735 1736 if mtflag == 'T': 1737 mt_member = selectInstrument(filename,output, pars=pars) 1738 mt_wcs = {} 1739 for member in mt_member.members: 1740 mt_wcs[member.chip] = member.geometry.wcs 1741 1742 pars['mt_wcs'] = mt_wcs 1743 del mt_member 1744 1745 for memname in prodlist['order']: 1746 pardict = self.pars[memname] 1747 pardict.update(pars) 1748 filename = fileutil.buildRootname(memname) 1749 if filename: 1750 self.members.append(selectInstrument(filename,output,pars=pardict)) 1751 else: 1752 print 'No recognizable input! Not building parameters for ',memname

1753

1754 - def getMember(self,memname):

1755 """ Return the class instance for the member with name memname.""" 1756 member = None 1757 for mem in self.members: 1758 member = mem.getMember(memname) 1759 if member != None: 1760 break 1761 return member

1762

1763 - def getMemberNames(self):

1764 """ Returns a dictionary with one key for each member. 1765 The value for each key is a list of all the extension/chip 1766 names that make up each member. 1767 Output: {member1:[extname1,extname2,...],...} 1768 """ 1769 memlist = [] 1770 for member in self.members: 1771 memlist.extend(member.getMemberNames()) 1772 1773 return memlist

1774

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

1776 1777 # If an output field is specified, update product WCS to match 1778 if ref != None: 1779 _field = ref 1780 _field.mergeWCS(self.product.geometry.wcslin) 1781 #_field.exptime = self.exptime 1782 1783 # Update product WCS based on input 1784 self.transformMetachip(_field) 1785 #_field.mergeWCS(self.product.geometry.wcs) 1786 else: 1787 _field = SkyField() 1788 _field.mergeWCS(self.product.geometry.wcslin) 1789 # Reset Product WCS to reflect use of default WCS 1790 self.product.geometry.wcs = self.product.geometry.wcslin.copy() 1791 _field.exptime = self.exptime 1792 # Specify that this product comes represents a dither product 1793 # not just a single observation product. 1794 _field.dither = yes 1795 1796 # Copy the new reference WCS into the product WCS 1797 self.product.geometry.wcs = _field.wcs.copy() 1798 1799 parlist = [] 1800 for member in self.members: 1801 parlist = parlist + member.buildPars(ref=_field) 1802 1803 # Set value for number of images combined by PyDrizzle 1804 # based on DitherPattern attribute. 1805 for pl in parlist: 1806 pl['nimages'] = self.nimages 1807 1808 return parlist

1809

1810 - def buildMetachip(self):

1811 """ 1812 This method combines the results of the member's buildMetachip() 1813 methods into a composite WCS. 1814 (DitherProduct method) 1815 """ 1816 prodlist = [] 1817 1818 _ref = self.members[0].product.geometry 1819 1820 # Get pscale from model 1821 pscale = _ref.model.pscale 1822 1823 for member in self.members: 1824 # merge the corrected shapes into a corrected meta-chip here 1825 # Start by computing the corected positions for reference points 1826 prodlist.append(member.product) 1827 1828 # Use first input image WCS as initial reference WCS 1829 meta_wcs = _ref.wcs.copy() 1830 1831 ref_crval = (meta_wcs.crval1,meta_wcs.crval2) 1832 ref_crpix = (meta_wcs.crpix1,meta_wcs.crpix2) 1833 1834 # Specify the output orientation angle 1835 ref_orient = meta_wcs.orient 1836 1837 _delta_x = _delta_y = 0.0 1838 _xr,_yr,_dpx,_dpy = [],[],[],[] 1839 # Compute offsets for each input relative to this reference WCS 1840 for member in prodlist: 1841 _wcs = member.geometry.wcs 1842 1843 # Rigorously compute the orientation changes from WCS 1844 # information using algorithm provided by R. Hook from WDRIZZLE. 1845 abxt,cdyt = drutil.wcsfit(member.geometry, meta_wcs) 1846 1847 # Compute the rotation between input and reference from fit coeffs. 1848 _angle = RADTODEG(N.arctan2(abxt[1],cdyt[0])) 1849 _dpos = (abxt[2],cdyt[2]) 1850 1851 _delta_x += _dpos[0] 1852 _delta_y += _dpos[1] 1853 _dpx.append(_dpos[0]) 1854 _dpy.append(_dpos[1]) 1855 1856 # Compute the range of pixels each input spans in the 1857 # output meta-frame coordinates 1858 # Each product in this list could have a different plate scale. 1859 # apply 1860 _scale = meta_wcs.pscale / _wcs.pscale 1861 1862 # Now, account for any rotation difference between 1863 # input and output frames 1864 #_angle = meta_wcs.orient - _wcs.orient 1865 _xrange,_yrange = drutil.getRotatedSize(member.corners['corrected'],_angle) 1866 1867 #_range = [(_xrange[1] - _xrange[0]),(_yrange[1] - _yrange[0])] 1868 #_range[0] *= _scale 1869 #_range[1] *= _scale 1870 1871 _xr.append(_dpos[0] + _xrange[0]*_scale) 1872 _xr.append(_dpos[0] + _xrange[1]*_scale) 1873 _yr.append(_dpos[1] + _yrange[0]*_scale) 1874 _yr.append(_dpos[1] + _yrange[1]*_scale) 1875 1876 # Determine the full size of the metachip 1877 _xmin = N.minimum.reduce(_xr) 1878 _ymin = N.minimum.reduce(_yr) 1879 _xmax = N.maximum.reduce(_xr) 1880 _ymax = N.maximum.reduce(_yr) 1881 1882 _dxmin = N.minimum.reduce(_dpx) 1883 _dymin = N.minimum.reduce(_dpy) 1884 _dxmax = N.maximum.reduce(_dpx) 1885 _dymax = N.maximum.reduce(_dpy) 1886 1887 _nimg = len(prodlist) 1888 _delta_x /= _nimg 1889 _delta_y /= _nimg 1890 1891 # Computes the offset from center of the overall set of shifts 1892 # as applied to the pixels. This insures that the visible pixels 1893 # are centered in the output. 1894 _delta_dx = ((_xmax - _delta_x) + (_xmin - _delta_x))/2. 1895 _delta_dy = ((_ymax - _delta_y) + (_ymin - _delta_y))/2. 1896 1897 if _xmin > 0.: _xmin = 0. 1898 if _ymin > 0.: _ymin = 0. 1899 1900 # Need to account for overall offset in images. 1901 # by adding in average offset induced by shifts to 1902 # output size: delta_dx/dy. 1903 # This accounts for relative offsets that are not centered 1904 # on the final output image. 1905 #xsize = int(_xmax - _xmin + 2.0 + abs(_delta_dx) ) 1906 #ysize = int(_ymax - _ymin + 2.0 + abs(_delta_dy) ) 1907 xsize = int(_xmax - _xmin) 1908 ysize = int(_ymax - _ymin) 1909 1910 # Compute difference between new size and 1911 # original size of meta_wcs frame. 1912 _dxsize = xsize - meta_wcs.naxis1 1913 _dysize = ysize - meta_wcs.naxis2 1914 1915 # Update reference WCS for new size 1916 meta_wcs.naxis1 = xsize 1917 meta_wcs.naxis2 = ysize 1918 _cen = (float(xsize)/2.,float(ysize)/2.) 1919 1920 # Determine offset of centers from center of image. 1921 # d[x/y]size : difference between input and output frame sizes 1922 # delta_[x/y] : average of all shifts applied to input images 1923 # delta_d[x/y]: amount off-center of all shifts 1924 # 1925 # Need to take into account overall offset in shifts, such that 1926 # minimum, not maximum, shift is always 0.0. This is needed to allow 1927 # subarrays to align properly with full frame observations without 1928 # introducing an overall offset to the output. WJH 13 Sept 2004. 1929 # 1930 _nref = ( _dxsize/2. - _delta_x - _delta_dx, _dysize/2. - _delta_y - _delta_dy) 1931 1932 # Have to adjust the CRPIX by how much the center needs to shift 1933 # to fit into the reference frame. 1934 meta_wcs.crpix1 = meta_wcs.crpix1 + _nref[0] 1935 meta_wcs.crpix2 = meta_wcs.crpix2 + _nref[1] 1936 meta_wcs.recenter() 1937 1938 return meta_wcs

1939 1940 # 1941 # 1942 # Set up default pars dictionary for external calls 1943 _default_pars = {'psize':None,'rot':None,'idckey':None} 1944

1945 -def selectInstrument(filename,output,pars=_default_pars):

1946 """ 1947 Method which encapsulates the logic for determining 1948 which class to instantiate for each file. 1949 1950 """ 1951 # Determine the instrument... 1952 instrument = fileutil.getKeyword(filename+'[0]','INSTRUME') 1953 1954 #try: 1955 # ... then create an appropriate object. 1956 if instrument == INSTRUMENT[0]: 1957 member = ACSObservation(filename,output,pars=pars) 1958 elif instrument == INSTRUMENT[1]: 1959 member = WFPCObservation(filename,output,pars=pars) 1960 elif instrument == INSTRUMENT[2]: 1961 member = STISObservation(filename,output,pars=pars) 1962 elif instrument == INSTRUMENT[3]: 1963 member = NICMOSObservation(filename,output,pars=pars) 1964 else: 1965 #raise AttributeError, "Instrument '%s' not supported now."%instrument 1966 member = GenericObservation(filename,output,pars=pars) 1967 #except: 1968 # raise IOError,"Image %s could not be processed."%filename 1969 1970 return member

1971

1972 -class SkyField:

1973 """ 1974 An class for specifying the parameters and building a WCS object 1975 for a user-specified drizzle product. 1976 The user may optionally modify the values for: 1977 psize - size of image's pixels in arc-seconds 1978 orient - value of ORIENTAT for the field 1979 shape - tuple containing the sizes of the field's x/y axes 1980 ra,dec - position of the center of the field 1981 decimal (124.5678) or 1982 sexagesimal string _in quotes_ ('hh:mm:ss.ss') 1983 crpix - tuple for pixel position of reference point in 1984 output image 1985 Usage: 1986 To specify a new field with a fixed output size of 1024x1024: 1987 --> field = pydrizzle.SkyField(shape=(1024,1024)) 1988 1989 The 'set()' method modifies one of the parameters listed above 1990 without affecting the remainder of the parameters. 1991 --> field.set(psize=0.1,orient=0.0) 1992 --> field.set(ra=123.45678,dec=0.1000,crpix=(521,576)) 1993 1994 View the WCS or user-specified values for this object: 1995 --> print field 1996 1997 """

1998 - def __init__(self,shape=None,psize=None,wcs=None):

1999 2000 self.shape = shape 2001 self.ra = None 2002 self.dec = None 2003 self.orient = None 2004 self.psize = psize 2005 2006 self.crpix = None 2007 2008 # Set up proper shape tuple for WCSObject 2009 if shape != None and psize != None: 2010 wshape = (shape[0],shape[1],psize) 2011 else: 2012 wshape = None 2013 if wcs: 2014 self.crpix = (wcs.crpix1,wcs.crpix2) 2015 2016 # Specify 'new=yes' with given rootname to unambiguously create 2017 # a new WCS from scratch. 2018 if wcs == None: 2019 self.wcs = wcsutil.WCSObject("New",new=yes) 2020 else: 2021 self.wcs = wcs.copy() 2022 # Need to adjust WCS to match the 'drizzle' task 2023 # align=center conventions. WJH 29-Aug-2005 2024 #self.wcs.crpix1 -= 1.0 2025 #self.wcs.crpix2 -= 1.0 2026 self.wcs.recenter() 2027 2028 self.wcs.updateWCS(size=shape,pixel_scale=psize) 2029 2030 # Set this to keep track of total exposure time for 2031 # output frame... Not user settable. 2032 self.exptime = None 2033 # Keep track of whether this is for a Dither product or not 2034 self.dither = None

2035

2036 - def mergeWCS(self,wcs,overwrite=yes):

2037 """ Sets up the WCS for this object based on another WCS. 2038 This method will NOT update object attributes other 2039 than WCS, as all other attributes reflect user-settings. 2040 """ 2041 # 2042 # Start by making a copy of the input WCS... 2043 # 2044 if self.wcs.rootname == 'New': 2045 self.wcs = wcs.copy() 2046 else: 2047 return 2048 self.wcs.recenter() 2049 2050 if self.ra == None: 2051 _crval = None 2052 else: 2053 _crval = (self.ra,self.dec) 2054 2055 if self.psize == None: 2056 _ratio = 1.0 2057 _psize = None 2058 # Need to resize the WCS for any changes in pscale 2059 else: 2060 _ratio = wcs.pscale / self.psize 2061 _psize = self.psize 2062 2063 if self.orient == None: 2064 _orient = None 2065 _delta_rot = 0. 2066 else: 2067 _orient = self.orient 2068 _delta_rot = wcs.orient - self.orient 2069 2070 _mrot = fileutil.buildRotMatrix(_delta_rot) 2071 2072 if self.shape == None: 2073 _corners = N.array([[0.,0.],[wcs.naxis1,0.],[0.,wcs.naxis2],[wcs.naxis1,wcs.naxis2]]) 2074 _corners -= (wcs.naxis1/2.,wcs.naxis2/2.) 2075 _range = drutil.getRotatedSize(_corners,_delta_rot) 2076 shape = ((_range[0][1] - _range[0][0])*_ratio,(_range[1][1]-_range[1][0])*_ratio) 2077 old_shape = (wcs.naxis1*_ratio,wcs.naxis2*_ratio) 2078 2079 _cen = (shape[0]/2., shape[1]/2.) 2080 2081 #if _delta_rot == 0.: 2082 # _crpix = (self.wcs.crpix1,self.wcs.crpix2) 2083 #else: 2084 # Rotate original scaled crpix position to new orientation 2085 #_crpix = N.dot((wcs.crpix1*_ratio - _cen[0],wcs.crpix2*_ratio -_cen[1]),_mrot)+_cen 2086 _crpix = _cen 2087 else: 2088 shape = self.shape 2089 if self.crpix == None: 2090 _crpix = (self.shape[0]/2.,self.shape[1]/2.) 2091 else: 2092 _crpix = self.crpix 2093 2094 # Set up the new WCS based on values from old one. 2095 self.wcs.updateWCS(pixel_scale=_psize,orient=_orient,refpos=_crpix,refval=_crval) 2096 self.wcs.naxis1 = int(shape[0]) 2097 self.wcs.naxis2 = int(shape[1])

2098

2099 - def set(self,psize=None,orient=None,ra=None,dec=None, 2100 shape=None,crpix=None):

2101 """ 2102 Modifies the attributes of the SkyField and 2103 updates it's WCS when appropriate. 2104 """ 2105 # Converts(if necessary), then updates the RA and Dec. 2106 _ra,_dec = None,None 2107 if ra != None: 2108 if string.find(repr(ra),':') > 0: 2109 _hms = string.split(repr(ra)[1:-1],':') 2110 if _hms[0][0] == '-': _sign = -1 2111 else: _sign = 1 2112 2113 for i in range(len(_hms)): _hms[i] = float(_hms[i]) 2114 _ra = _sign * (_hms[0] + ((_hms[1] + _hms[2]/60.) / 60.)) * 15. 2115 else: 2116 _ra = float(ra) 2117 self.ra = _ra 2118 2119 if dec != None: 2120 if string.find(repr(dec),':') > 0: 2121 _dms = string.split(repr(dec)[1:-1],':') 2122 if _dms[0][0] == '-': _sign = -1 2123 else: _sign = 1 2124 2125 for i in range(len(_dms)): _dms[i] = float(_dms[i]) 2126 _dec = _sign * (_dms[0] + ((_dms[1] + _dms[2]/60.) / 60.)) 2127 else: 2128 _dec = float(dec) 2129 self.dec = _dec 2130 2131 if self.ra != None and self.dec != None: 2132 _crval = (self.ra,self.dec) 2133 else: 2134 _crval = None 2135 2136 # Updates the shape, and reference position, 2137 # only if a new value is specified. 2138 _crpix = None 2139 if crpix == None: 2140 if shape != None: 2141 self.shape = shape 2142 _crpix = (self.shape[0]/2.,self.shape[1]/2.) 2143 else: 2144 _crpix = crpix 2145 2146 self.crpix=_crpix 2147 2148 if psize != None: 2149 self.psize = psize 2150 2151 if orient != None: 2152 self.orient = orient 2153 2154 # Updates the WCS with all the changes, if there is enough info. 2155 self.wcs.updateWCS(pixel_scale=psize,orient=orient,refpos=_crpix, 2156 refval=_crval,size=self.shape)

2157

2158 - def __str__(self):

2159 """ Prints the WCS information set for this object. 2160 """ 2161 if self.psize != None and self.orient != None: 2162 block = self.wcs.__str__() 2163 else: 2164 block = 'User parameters for SkyField object: \n' 2165 block = block + ' ra = '+repr(self.ra)+' \n' 2166 block = block + ' dec = '+repr(self.dec)+' \n' 2167 block = block + ' shape = '+repr(self.shape)+' \n' 2168 block = block + ' crpix = '+repr(self.crpix)+' \n' 2169 block = block + ' psize = '+repr(self.psize)+' \n' 2170 block = block + ' orient = '+repr(self.orient)+' \n' 2171 block = block + ' No WCS.\n' 2172 2173 return block

2174

2175 - def help(self):

2176 """ Creates and prints usage information for this class. 2177 """ 2178 print self.__doc__

2179

2180 -class PyDrizzle:

2181 """ 2182 Program to process and/or dither-combine image(s) using (t)drizzle. 2183 To create an object named 'test' that corresponds to a drizzle product: 2184 --> test = pydrizzle.PyDrizzle(input) 2185 where input is the FULL filename of an ACS observation or ASN table. 2186 This computes all the parameters necessary for running drizzle on all 2187 the input images. Once this object is created, you can run drizzle using: 2188 --> test.run() 2189 2190 The 'clean()' method can be used to remove files which would interfere with 2191 running Drizzle again using the 'run()' method after a product has already 2192 been created. 2193 2194 Optional parameters: 2195 output User-specified name for output products 2196 field User-specified parameters for output image 2197 includes: psize, orient, ra, dec, shape 2198 units Units for final product: 'counts' or 'cps'(DEFAULT) 2199 section Extension/group to be drizzled: list or single FITS extension(s) 2200 or group(s) syntax ('1' or 'sci,1') or None (DEFAULT: Use all chips). 2201 kernel Specify which kernel to use in TDRIZZLE 2202 'square'(default),'point','gaussian','turbo','tophat' 2203 pixfrac drizzle pixfrac value (Default: 1.0) 2204 idckey User-specified keyword for determining IDCTAB filename 2205 'IDCTAB'(ACS default),'TRAUGER'(WFPC2),'CUBIC'(WFPC2) 2206 idcdir User-specified directory for finding coeffs files: 2207 'drizzle$coeffs' (default) 2208 bits_single Specify DQ values to be considered good when 2209 drizzling with 'single=yes'. (Default: 0) 2210 bits_final Specify DQ values to be considered good when 2211 drizzling with 'single=no'. (Default: 0) 2212 Bits parameters will be interpreted as: 2213 None - No DQ information to be used, no mask created 2214 Int - sum of DQ values to be considered good 2215 2216 Optional Parameters for '.run()': 2217 build create multi-extension output: yes (Default) or no 2218 save keeps the individual inputs from drizzle: yes or no (Default) 2219 single drizzle to separate output for each input: yes or no (Default) 2220 blot run blot on drizzled products: yes or no (Default) 2221 clean remove coeffs and static mask files: yes or no (Default) 2222 2223 Optional Parameters for '.clean()': 2224 coeffs Removes coeffs and static mask files: yes or no (Default) 2225 final Removes final product: yes or no (Default) 2226 2227 Usage of optional parameters: 2228 --> test = pydrizzle.PyDrizzle('test_asn.fits',units='counts') 2229 To keep the individual 'drizzle' output products: 2230 --> test.run(save=yes) 2231 2232 Output frame parameters can be modified 'on-the-fly' using 'resetPars'. 2233 Given an already drizzled image 'refimg_drz.fits' as a reference, 2234 reset drizzle parameters using: 2235 --> wcsref = wcsutil.WCSObject('refimg_drz.fits[sci,1]') 2236 --> f = pydrizzle.SkyField(wcs=wcsref) 2237 Use either: 2238 --> test.resetPars(wcsref) 2239 Or: 2240 --> test.resetPars(f) 2241 Return to default parameters using no parameters at all: 2242 --> test.resetPars() 2243 More help on SkyField objects and their parameters can be obtained using: 2244 --> f.help() 2245 """

2246 - def __init__(self, input, output=None, field=None, units=None, section=None, 2247 kernel=None,pixfrac=None,bits_final=0,bits_single=0, 2248 wt_scl='exptime', fillval=0.,idckey='', in_units='counts', 2249 idcdir=DEFAULT_IDCDIR,memmap=0,dqsuffix=None,prodonly=yes,shiftfile=None):

2250 2251 if idcdir == None: idcdir = DEFAULT_IDCDIR 2252 2253 print 'Starting PyDrizzle Version ',__version__,' at ', _ptime() 2254 2255 # Do a quick sanity check on the input filename 2256 # Can it be found in the current directory? 2257 if isinstance(input,list) == False: 2258 _input = input 2259 else: 2260 _input = input[0] 2261 2262 # Insure that the section parameter always becomes a list 2263 if isinstance(section,list) == False and section != None: 2264 section = [section] 2265 2266 if not fileutil.findFile(_input): 2267 raise IOError,'Can not find input file in current directory!' 2268 2269 # Does it have a recognizable extension 2270 if string.find(_input,'.') < 0: 2271 raise ValueError,"Please specify extension (i.e., .fits) for input '%s' "%input 2272 2273 # These parameters are needed for buildPars() 2274 self.input = input 2275 2276 # Set the default value for 'build' 2277 self.build = yes 2278 2279 # Extract user-specified parameters, if any have been set... 2280 # 'field' will be a SkyField object... 2281 if field != None: 2282 psize = field.psize 2283 orient = field.orient 2284 else: 2285 psize = None 2286 orient = None 2287 2288 # These can also be set by the user. 2289 # Minimum set needed: psize, rot, and idckey 2290 self.pars = {'psize':psize,'units':units,'kernel':kernel,'rot':orient, 2291 'pixfrac':pixfrac,'idckey':idckey,'wt_scl':wt_scl, 2292 'fillval':fillval,'section':section, 'idcdir':idcdir+os.sep, 2293 'memmap':memmap,'dqsuffix':dqsuffix, 'in_units':in_units, 2294 'bits':[bits_final,bits_single], 'mt_wcs': None} 2295 2296 # Check to see if user-supplied output name is complete 2297 # Append .FITS suffix to output name if necessary 2298 self.output = output 2299 #if output != None: 2300 # _indx = string.rfind(output,'.') 2301 # if _indx < 0: output = output + '.fits' 2302 2303 # Watch out for any errors. 2304 # If they arise, all open files need to be closed... 2305 self.observation = None 2306 self.parlist = None 2307 #try: 2308 2309 # Decipher input to determine whether we are working with 2310 # an ASN table or single image... 2311 # 2312 # Should this check for ASN table and 2313 if (isinstance(input, list) == False) and \ 2314 (input.find('_asn') < 0 and input.find('_asc') < 0) : 2315 # Start with single image case... 2316 # Check to see if an output name was provided. 2317 if output == None: 2318 # We need to build a default output name... 2319 output = fileutil.buildNewRootname(input,extn='_drz.fits') 2320 print 'Setting up default output name: ',output 2321 elif output.rfind('.fits') < 0: 2322 if output.rfind('_drz') < 0: 2323 output += '_drz' 2324 2325 output += '.fits' 2326 2327 self.observation = selectInstrument(input,output,pars=self.pars) 2328 else: 2329 # We are dealing with an ASN table... 2330 # 'input' - full filename for ASN table 2331 if isinstance(input, list): 2332 asndict = fileutil.buildAsnDict(input, output=output,shiftfile=shiftfile) 2333 else: 2334 asndict = fileutil.readAsnTable(input,output=output,prodonly=prodonly) 2335 2336 # Build output filename 2337 if output == None: 2338 output = fileutil.buildNewRootname(asndict['output'],extn='_drz.fits') 2339 print 'Setting up output name: ',output 2340 2341 if len(asndict['order']) > 1: 2342 self.observation = DitherProduct(asndict,pars=self.pars) 2343 else: 2344 inroot = asndict['order'][0] 2345 pardict = asndict['members'][inroot] 2346 infile = fileutil.buildRootname(inroot) 2347 if infile == None: 2348 raise IOError,'No product found for association table.' 2349 # Append user specified parameters to shifts dictionary 2350 pardict.update(self.pars) 2351 self.observation = selectInstrument(infile,output,pars=pardict) 2352 self.output = output 2353 2354 # This call puts together the parameters for the input image 2355 # with those for the output to create a final parameter list 2356 # for running 'drizzle'. 2357 # It relies on the buildPars() methods for each exposure to 2358 # generate a complete set of parameters for all inputs 2359 # 2360 self.parlist = self.observation.buildPars(ref=field) 2361 2362 #finally: 2363 # Something went wrong, so we need to make sure all file 2364 # handles get closed... 2365 #self.close() 2366 #print 'PyDrizzle could not initialize due to errors.' 2367 self.observation.closeHandle() 2368 2369 # Let the user know parameters have been successfully calculated 2370 print 'Drizzle parameters have been calculated. Ready to .run()...' 2371 print 'Finished calculating parameters at ',_ptime()

2372

2373 - def clean(self,coeffs=no,final=no):

2374 """ Removes intermediate products from disk. """ 2375 for img in self.parlist: 2376 # If build == no, then we do not want to delete the only 2377 # products created by PyDrizzle; namely, 2378 # outdata, outcontext, outweight 2379 if self.build == yes: 2380 fileutil.removeFile([img['outdata'],img['outcontext'],img['outweight']]) 2381 2382 fileutil.removeFile([img['outsingle'],img['outsweight']]) 2383 #fileutil.removeFile([img['outsingle'],img['outsweight'],img['outscontext']]) 2384 fileutil.removeFile(img['outblot']) 2385 if coeffs: 2386 os.remove(img['coeffs']) 2387 if img['driz_mask'] != None: 2388 fileutil.removeFile(img['driz_mask']) 2389 if img['single_driz_mask'] != None: 2390 fileutil.removeFile(img['single_driz_mask']) 2391 if final: 2392 fileutil.removeFile(self.output)

2393 2394 2395 # Run 'drizzle' here... 2396 #

2397 - def run(self,save=no,build=yes,blot=no,single=no,clean=no,interp='linear',sinscl=1.0, debug=no):

2398 """Perform drizzle operation on input to create output. 2399 This method would rely on the buildPars() method for 2400 the output product to produce correct parameters 2401 based on the inputs. The output for buildPars() is a list 2402 of dictionaries, one for each input, that matches the 2403 primary parameters for an IRAF drizzle task. 2404 2405 This method would then loop over all the entries in the 2406 list and run 'drizzle' for each entry. """ 2407 2408 # Store the value of build set by the user for use, if desired, 2409 # in the 'clean()' method. 2410 self.build = build 2411 self.debug = debug 2412 2413 print 'PyDrizzle: drizzle task started at ',_ptime() 2414 _memmap = self.pars['memmap'] 2415 2416 # Check for existance of output file. 2417 if single == no and build == yes and fileutil.findFile(self.output): 2418 print 'Removing previous output product...' 2419 os.remove(self.output) 2420 # 2421 # Setup the versions info dictionary for output to PRIMARY header 2422 # The keys will be used as the name reported in the header, as-is 2423 # 2424 _versions = {'PyDrizzle':__version__,'PyFITS':pyfits.__version__,'Numpy':N.__version__} 2425 2426 # Set parameters for each input and run drizzle on it here. 2427 2428 if blot: 2429 # 2430 # Run blot on data... 2431 # 2432 2433 _hdrlist = [] 2434 2435 for plist in self.parlist: 2436 2437 _insci = N.zeros((plist['outny'],plist['outnx']),dtype=N.float32) 2438 _outsci = N.zeros((plist['blotny'],plist['blotnx']),dtype=N.float32) 2439 _hdrlist.append(plist) 2440 # Open input image as PyFITS object 2441 if plist['outsingle'] != plist['outdata']: 2442 _data = plist['outsingle'] 2443 else: 2444 _data = plist['outdata'] 2445 2446 # PyFITS can be used here as it will always operate on 2447 # output from PyDrizzle (which will always be a FITS file) 2448 # Open the input science file 2449 _fname,_sciextn = fileutil.parseFilename(_data) 2450 _inimg = fileutil.openImage(_fname) 2451 2452 # Return the PyFITS HDU corresponding to the named extension 2453 _scihdu = fileutil.getExtn(_inimg,_sciextn) 2454 _insci = _scihdu.data.copy() 2455 2456 # Read in the distortion correction arrays, if specified 2457 _pxg,_pyg = plist['exposure'].getDGEOArrays() 2458 2459 # Now pass numarray objects to callable version of Blot... 2460 #runBlot(plist) 2461 build=no 2462 misval = 0.0 2463 kscale = 1.0 2464 2465 xmin = 1 2466 xmax = plist['outnx'] 2467 ymin = 1 2468 ymax = plist['outny'] 2469 # 2470 # Convert shifts to input units 2471 # 2472 xsh = plist['xsh'] * plist['scale'] 2473 ysh = plist['ysh'] * plist['scale'] 2474 # ARRDRIZ.TBLOT needs to be updated to support 'poly5' interpolation, 2475 # and exptime scaling of output image. 2476 # 2477 """ 2478 # 2479 # This call to 'arrdriz.tdriz' uses the F2PY syntax 2480 # 2481 arrdriz.tblot(_insci, _outsci,xmin,xmax,ymin,ymax, 2482 plist['xsh'],plist['ysh'], 2483 plist['rot'],plist['scale'], kscale, _pxg, _pyg, 2484 'center',interp, plist['coeffs'], plist['exptime'], 2485 misval, sinscl, 1) 2486 # 2487 # End of F2PY syntax 2488 # 2489 """ 2490 # 2491 # This call to 'arrdriz.tdriz' uses the F2C syntax 2492 # 2493 if (_insci.dtype > N.float32): 2494 _insci = _insci.astype(N.float32) 2495 t = arrdriz.tblot(_insci, _outsci,xmin,xmax,ymin,ymax, 2496 xsh,ysh, plist['rot'],plist['scale'], kscale, 2497 0.0,0.0, 1.0,1.0, 0.0, 'output', 2498 _pxg, _pyg, 2499 'center',interp, plist['coeffs'], plist['exptime'], 2500 misval, sinscl, 1,plist['alpha'],plist['beta']) 2501 2502 # 2503 # End of F2C syntax 2504 # 2505 2506 # Write output Numarray objects to a PyFITS file 2507 # Blotting only occurs from a drizzled SCI extension 2508 # to a blotted SCI extension... 2509 _header = fileutil.getHeader(plist['data']) 2510 _wcs = wcsutil.WCSObject(plist['data'],header=_header) 2511 2512 _outimg = outputimage.OutputImage(_hdrlist, build=no, wcs=_wcs, blot=yes) 2513 _outimg.outweight = None 2514 _outimg.outcontext = None 2515 _outimg.writeFITS(plist['data'],_outsci,None,versions=_versions) 2516 2517 #_buildOutputFits(_outsci,None,plist['outblot']) 2518 _insci *= 0. 2519 _outsci *= 0. 2520 _inimg.close() 2521 del _inimg 2522 _hdrlist = [] 2523 2524 del _pxg,_pyg 2525 2526 del _insci,_outsci 2527 del _outimg 2528 2529 else: 2530 # 2531 # Perform drizzling... 2532 # 2533 # Only work on a copy of the product WCS, so that when 2534 # this gets updated for the output image, it does not 2535 # modify the original WCS computed by PyDrizzle 2536 _wcs = self.observation.product.geometry.wcs.copy() 2537 2538 _numctx = {'all':len(self.parlist)} 2539 # if single: 2540 # Determine how many chips make up each single image 2541 for plist in self.parlist: 2542 plsingle = plist['outsingle'] 2543 if _numctx.has_key(plsingle): _numctx[plsingle] += 1 2544 else: _numctx[plsingle] = 1 2545 # 2546 # A image buffer needs to be setup for converting the input 2547 # arrays (sci and wht) from FITS format to native format 2548 # with respect to byteorder and byteswapping. 2549 # This buffer should be reused for each input. 2550 # 2551 plist = self.parlist[0] 2552 _outsci = N.zeros((plist['outny'],plist['outnx']),dtype=N.float32) 2553 _outwht = N.zeros((plist['outny'],plist['outnx']),dtype=N.float32) 2554 _inwcs = N.zeros([8],dtype=N.float64) 2555 2556 # Compute how many planes will be needed for the context image. 2557 _nplanes = int((_numctx['all']-1) / 32) + 1 2558 # For single drizzling or when context is turned off, 2559 # minimize to 1 plane only... 2560 if single or self.parlist[0]['outcontext'] == '': 2561 _nplanes = 1 2562 2563 # Always initialize context images to a 3-D array 2564 # and only pass the appropriate plane to drizzle as needed 2565 _outctx = N.zeros((_nplanes,plist['outny'],plist['outnx']),dtype=N.int32) 2566 2567 # Keep track of how many chips have been processed 2568 # For single case, this will determine when to close 2569 # one product and open the next. 2570 _numchips = 0 2571 _nimg = 0 2572 _hdrlist = [] 2573 2574 for plist in self.parlist: 2575 # Read in the distortion correction arrays, if specifij8cw08n4q_raw.fitsed 2576 _pxg,_pyg = plist['exposure'].getDGEOArrays() 2577 2578 _hdrlist.append(plist) 2579 # Open the SCI image 2580 _expname = plist['data'] 2581 _handle = fileutil.openImage(_expname,mode='readonly',memmap=0) 2582 _fname,_extn = fileutil.parseFilename(_expname) 2583 _sciext = fileutil.getExtn(_handle,extn=_extn) 2584 2585 # Make a local copy of SCI array and WCS info 2586 #_insci = _sciext.data.copy() 2587 _inwcs = drutil.convertWCS(wcsutil.WCSObject(_fname,header=_sciext.header),_inwcs) 2588 2589 # Determine output value of BUNITS 2590 # and make sure it is not specified as 'ergs/cm...' 2591 _bunit = None 2592 if _sciext.header.has_key('BUNIT') and _sciext.header['BUNIT'].find('ergs') < 0: 2593 _bunit = _sciext.header['BUNIT'] 2594 2595 if _bunit is not None: 2596 _bindx = _bunit.find('/') 2597 if plist['units'] == 'cps': 2598 # If BUNIT value does not specify count rate already... 2599 if _bindx < 1: 2600 # ... append '/SEC' to value 2601 _bunit += '/SEC' 2602 else: 2603 # reset _bunit here to None so it does not 2604 # overwrite what is already in header 2605 _bunit = None 2606 else: 2607 if _bindx > 0: 2608 # remove '/SEC' 2609 _bunit = _bunit[:_bindx] 2610 else: 2611 # reset _bunit here to None so it does not 2612 # overwrite what is already in header 2613 _bunit = None 2614 2615 # Compute what plane of the context image this input would 2616 # correspond to: 2617 _planeid = int(_numchips /32) 2618 2619 # Select which mask needs to be read in for drizzling 2620 if single: 2621 _mask = plist['single_driz_mask'] 2622 else: 2623 _mask = plist['driz_mask'] 2624 2625 # Check to see whether there is a mask_array at all to use... 2626 if isinstance(_mask,types.StringType): 2627 if _mask != None and _mask != '': 2628 _wht_handle = fileutil.openImage(_mask,mode='readonly',memmap=0) 2629 _inwht = _wht_handle[0].data.astype(N.float32) 2630 _wht_handle.close() 2631 del _wht_handle 2632 else: 2633 print 'No weight or mask file specified! Assuming all pixels are good.' 2634 _inwht = N.ones((plist['blotny'],plist['blotnx']),dtype=N.float32) 2635 elif _mask != None: 2636 _inwht = _mask.astype(N.float32) 2637 else: 2638 print 'No weight or mask file specified! Assuming all pixels are good.' 2639 _inwht = N.ones((plist['blotny'],plist['blotnx']),dtype=N.float32) 2640 2641 if plist['wt_scl'] != None: 2642 if isinstance(plist['wt_scl'],types.StringType): 2643 if plist['wt_scl'].isdigit() == False : 2644 # String passed in as value, check for 'exptime' or 'expsq' 2645 _wtscl_float = None 2646 try: 2647 _wtscl_float = float(plist['wt_scl']) 2648 except ValueError: 2649 _wtscl_float = None 2650 if _wtscl_float != None: 2651 _wtscl = _wtscl_float 2652 elif plist['wt_scl'] == 'expsq': 2653 _wtscl = plist['exptime']*plist['exptime'] 2654 else: 2655 # Default to the case of 'exptime', if 2656 # not explicitly specified as 'expsq' 2657 _wtscl = plist['exptime'] 2658 else: 2659 # int value passed in as a string, convert to float 2660 _wtscl = float(plist['wt_scl']) 2661 else: 2662 # We have a non-string value passed in... 2663 _wtscl = float(plist['wt_scl']) 2664 else: 2665 # Default case: wt_scl = exptime 2666 _wtscl = plist['exptime'] 2667 2668 #print 'WT_SCL: ',plist['wt_scl'],' _wtscl: ',_wtscl 2669 # Set additional parameters needed by 'drizzle' 2670 _in_units = plist['in_units'] 2671 if _in_units == 'cps': 2672 _expin = 1.0 2673 else: 2674 _expin = plist['exptime'] 2675 _shift_fr = 'output' 2676 _shift_un = 'output' 2677 _uniqid = _numchips + 1 2678 ystart = 0 2679 nmiss = 0 2680 nskip = 0 2681 _vers = plist['driz_version'] 2682 2683 _con = yes 2684 _imgctx = _numctx['all'] 2685 if single: 2686 _imgctx = _numctx[plist['outsingle']] 2687 #if single or (plist['outcontext'] == '' and single == yes): 2688 if _nplanes == 1: 2689 _con = no 2690 # We need to reset what gets passed to TDRIZ 2691 # when only 1 context image plane gets generated 2692 # to prevent overflow problems with trying to access 2693 # planes that weren't created for large numbers of inputs. 2694 _planeid = 0 2695 _uniqid = ((_uniqid-1) % 32) + 1 2696 2697 """ 2698 # 2699 # This call to 'arrdriz.tdriz' uses the F2PY syntax 2700 # 2701 #_dny = plist['blotny'] 2702 # Call 'drizzle' to perform image combination 2703 tdriz,nmiss,nskip,_vers = arrdriz.tdriz( 2704 _sciext.data,_inwht, _outsci, _outwht, 2705 _outctx[_planeid], _con, _uniqid, ystart, 1, 1, 2706 plist['xsh'],plist['ysh'], 'output','output', 2707 plist['rot'],plist['scale'], _pxg,_pyg, 2708 'center', plist['pixfrac'], plist['kernel'], 2709 plist['coeffs'], 'counts', _expin,_wtscl, 2710 plist['fillval'], _inwcs, 1, nmiss, nskip,_vers) 2711 # 2712 # End of F2PY syntax 2713 # 2714 """ 2715 # 2716 # This call to 'arrdriz.tdriz' uses the F2C syntax 2717 # 2718 _dny = plist['blotny'] 2719 # Call 'drizzle' to perform image combination 2720 if (_sciext.data.dtype > N.float32): 2721 _sciext.data = _sciext.data.astype(N.float32) 2722 _vers,nmiss,nskip = arrdriz.tdriz(_sciext.data,_inwht, _outsci, _outwht, 2723 _outctx[_planeid], _uniqid, ystart, 1, 1, _dny, 2724 plist['xsh'],plist['ysh'], 'output','output', 2725 plist['rot'],plist['scale'], 2726 0.0,0.0, 1.0,1.0,0.0,'output', 2727 _pxg,_pyg, 'center', plist['pixfrac'], plist['kernel'], 2728 plist['coeffs'], _in_units, _expin,_wtscl, 2729 plist['fillval'], _inwcs, nmiss, nskip, 1,plist['alpha'],plist['beta']) 2730 """ 2731 _vers,nmiss,nskip = arrdriz.tdriz(_sciext.data,_inwht, _outsci, _outwht, 2732 _outctx[_planeid], _uniqid, ystart, 1, 1, _dny, 2733 plist['xsh'],plist['ysh'], 'output','output', 2734 plist['rot'],plist['scale'], 2735 _pxg,_pyg, 'center', plist['pixfrac'], plist['kernel'], 2736 plist['coeffs'], 'counts', _expin,_wtscl, 2737 plist['fillval'], _inwcs, nmiss, nskip, 1) 2738 """ 2739 # 2740 # End of F2C syntax 2741 # 2742 2743 plist['driz_version'] = _vers 2744 2745 if nmiss > 0: 2746 print '! Warning, ',nmiss,' points were outside the output image.' 2747 if nskip > 0: 2748 print '! Note, ',nskip,' input lines were skipped completely.' 2749 # Close image handle 2750 _handle.close() 2751 del _handle,_fname,_extn,_sciext 2752 del _inwht 2753 2754 del _pxg,_pyg 2755 2756 2757 if _nimg == 0 and self.debug == yes: 2758 # Only update the WCS from drizzling the 2759 # first image in the list, just like IRAF DRIZZLE 2760 drutil.updateWCS(_inwcs,_wcs) 2761 #print '[run] Updated WCS now:' 2762 #print _wcs 2763 2764 # Increment number of chips processed for single output 2765 _numchips += 1 2766 if _numchips == _imgctx: 2767 ########################### 2768 # 2769 # IMPLEMENTATION REQUIREMENT: 2770 # 2771 # Need to implement scaling of the output image 2772 # from 'cps' to 'counts' in the case where 'units' 2773 # was set to 'counts'... 21-Mar-2005 2774 # 2775 ########################### 2776 # Start by determining what exposure time needs to be used 2777 # to rescale the product. 2778 if single: 2779 _expscale = plist['exptime'] 2780 else: 2781 _expscale = plist['texptime'] 2782 2783 #If output units were set to 'counts', rescale the array in-place 2784 if plist['units'] == 'counts': 2785 N.multiply(_outsci, _expscale, _outsci) 2786 2787 # 2788 # Write output arrays to FITS file(s) and reset chip counter 2789 # 2790 _outimg = outputimage.OutputImage(_hdrlist, build=build, wcs=_wcs, single=single) 2791 _outimg.bunit = _bunit 2792 _outimg.writeFITS(plist['data'],_outsci,_outwht,ctxarr=_outctx,versions=_versions) 2793 del _outimg 2794 # 2795 # Reset chip counter for next output image... 2796 # 2797 _numchips = 0 2798 _nimg = 0 2799 N.multiply(_outsci,0.,_outsci) 2800 N.multiply(_outwht,0.,_outwht) 2801 N.multiply(_outctx,0,_outctx) 2802 2803 _hdrlist = [] 2804 else: 2805 _nimg += 1 2806 2807 del _outsci,_outwht,_inwcs,_outctx, _hdrlist 2808 2809 # Remove temp files, if desired 2810 # Files to be removed are: 2811 # parlist['coeffs'] 2812 # parlist['driz_mask'] 2813 if not save and clean: 2814 for img in self.parlist: 2815 fileutil.removeFile(img['coeffs']) 2816 if img['driz_mask'] != None: 2817 fileutil.removeFile(img['driz_mask']) 2818 if img['single_driz_mask'] != None: 2819 fileutil.removeFile(img['single_driz_mask']) 2820 2821 print 'PyDrizzle drizzling completed at ',_ptime()

2822 2823

2824 - def resetPars(self,field=None,pixfrac=None,kernel=None,units=None):

2825 """ 2826 Recompute the output parameters based on a new 2827 SkyField or WCSObject object. 2828 """ 2829 if field and not isinstance(field, SkyField): 2830 if isinstance(field, wcsutil.WCSObject): 2831 _ref = SkyField(wcs=field) 2832 else: 2833 raise TypeError, 'No valid WCSObject or SkyField object entered...' 2834 else: 2835 _ref = field 2836 # Create new version of the parlist with the new values of the 2837 # parameters based on the reference image. 2838 new_parlist = self.observation.buildPars(ref=_ref) 2839 2840 2841 # Copy the parameters from the new parlist into the existing 2842 # parlist (self.parlist) to preserve any changes/updates made 2843 # to the parlist externally to PyDrizzle. 2844 for i in xrange(len(self.parlist)): 2845 for key in new_parlist[i]: 2846 self.parlist[i][key] = new_parlist[i][key] 2847 del new_parlist 2848 2849 if pixfrac or kernel or units: 2850 #print 'resetting additional parameter(s)...' 2851 for p in self.parlist: 2852 if kernel: 2853 p['kernel'] = kernel 2854 if pixfrac: 2855 p['pixfrac'] = pixfrac 2856 if units: 2857 if units == 'counts' or units == 'cps': 2858 p['units'] = units 2859 else: 2860 print 'Units ',units,' not valid! Parameter not reset.' 2861 print 'Please use either "cps" or "counts".'

2862

2863 - def help(self):

2864 """ 2865 Run the module level help function to provide syntax information. 2866 """ 2867 help()

2868

2869 - def printPars(self,pars='xsh,ysh,rot,scale,outnx,outny,data',format=no):

2870 """ Prints common parameters for review. """ 2871 if format: 2872 _title = pars.replace(',',' ') 2873 print _title 2874 print '-'*72 2875 2876 _pars = pars.split(',') 2877 for pl in self.parlist: 2878 for _p in _pars: print pl[_p], 2879 print ''

2880

2881 - def getMember(self,memname):

2882 """ Returns the class instance for the specified member name.""" 2883 return self.observation.getMember(memname)

2884 2885

2886 -def _buildOutputFits(sci,wht,fname,ctx=None,extlist=['SCI','ERR','DQ']):

2887 2888 # Remove previous output products, rather than appending to them 2889 if fileutil.findFile(fname): 2890 fileutil.removeFile(fname) 2891 2892 # Setup primary header as an HDU ready for appending to output FITS file 2893 prihdu = pyfits.PrimaryHDU() 2894 scihdu = pyfits.ImageHDU(data=sci,name=extlist[0]) 2895 whthdu = pyfits.ImageHDU(data=wht,name=extlist[1]) 2896 ctxhdu = pyfits.ImageHDU(data=ctx,name=extlist[2]) 2897 2898 fimg = pyfits.open(fname,mode='append') 2899 fimg.append(prihdu) 2900 fimg.append(scihdu) 2901 fimg.append(whthdu) 2902 fimg.append(ctxhdu) 2903 fimg.close() 2904 del fimg

2905 2906 # End of 'runDrizzle'

2907 -def help():

2908 print PyDrizzle.__doc__

2909

Source Code for Package pydrizzle