pydrizzle.exposure

42 """ 43 This class will provide the basic functionality for keeping 44 track of an exposure's parameters, including distortion model, 45 WCS information, and metachip shape. 46 """ 47

48 - def __init__(self,expname, handle=None, dqname=None, idckey=None, 49 new=no,wcs=None,mask=None,pa_key=None, parity=None, 50 idcdir=None, rot=None, extver=1, exptime=None, 51 ref_pscale=1.0, binned=1, mt_wcs=None, group_indx = None):

52 53 # This name should be formatted for use in image I/O 54 self.name = fileutil.osfn(expname) 55 56 # osfn() will expand '.' unnecessarily, potentially 57 # creating a string-length problem for 'drizzle', which 58 # is limited to strings of 80 chars. 59 _path,_name = os.path.split(self.name) 60 # if path for this filename is the same as the current dir, 61 # then there is no need to pass along the path. 62 if _path == os.getcwd(): self.name = _name 63 64 65 # Keep track of any associated mask file created for 66 # this exposure from its DQ file, or other mask file. 67 _fname,_extn = fileutil.parseFilename(expname) 68 _open = False 69 70 # Make sure we have an open file handle to use for getting the 71 # header and data arrays. 72 if not handle and not new: 73 handle = fileutil.openImage(expname) 74 _open = True 75 76 # If no extension was specified, try to interrogate the file 77 # to find whether the SCI array is in the Primary 78 # (as in Simple FITS) or first extension (as in MEF). 79 if handle and _extn == None: 80 if handle[0].data == None: 81 # Primary extension specified and no data present. 82 # Try looking for data in next extension. 83 if len(handle) > 1 and handle[1].data != None: 84 _extn = 1 85 expname += '[1]' 86 else: 87 raise IOError, "No valid image data in %s.\n"%expname 88 else: 89 _extn = 0 90 91 self.dgeoname = None 92 self.xgeoim = "" 93 self.ygeoim = "" 94 self.exptime = exptime 95 self.group_indx = group_indx 96 if not new: 97 # Read in a copy of the header for this exposure/group/extension 98 _header = fileutil.getHeader(expname,handle=handle) 99 _chip = drutil.getChipId(_header) 100 self.chip = str(_chip) 101 # Keep track of any distortion correction images provided 102 # for this chip 103 self.dgeoname = fileutil.getKeyword(expname,'DGEOFILE',handle=handle) 104 self.xgeoim,self.ygeoim = self.getDGEOExtn() 105 if self.exptime == None: 106 self.exptime = float(_header['EXPTIME']) 107 if self.exptime == 0.: self.exptime = 1.0 108 # 109 # Extract photometric transformation keywords 110 # If they do not exist, use default values of 0 and 1 111 # 112 self.plam = float(fileutil.getKeyword(expname,'PHOTPLAM',handle=handle)) / 10. 113 if self.plam == None: 114 # Setup a default value in case this keyword does not exist 115 self.plam = 555. 116 self.photzpt = float(fileutil.getKeyword(expname,'PHOTZPT',handle=handle)) 117 if self.photzpt == None: self.photzpt = 0.0 118 self.photflam = float(fileutil.getKeyword(expname,'PHOTFLAM',handle=handle)) 119 if self.photflam == None: self.photflam = 1.0 120 121 # Read in date-obs from primary header 122 if _header: 123 if _header.has_key('date-obs'): 124 self.dateobs = _header['date-obs'] 125 elif _header.has_key('date_obs'): 126 self.dateobs = _header['date_obs'] 127 else: 128 self.dateobs = None 129 else: 130 self.dateobs = None 131 132 # Initialize the value of BUNIT based on the header information, if 133 # the header has the keyword 134 if _header.has_key('BUNIT') and _header['BUNIT'].find('ergs') < 0: 135 self.bunit = _header['BUNIT'] 136 else: 137 self.bunit = 'ELECTRONS' 138 139 else: 140 _chip = 1 141 _header = None 142 self.chip = str(_chip) 143 # Set a default value for pivot wavelength 144 self.plam = 555. 145 self.photzpt = 0.0 146 self.photflam = 1.0 147 self.dateobs = None 148 if self.exptime == None: 149 self.exptime = 1. 150 151 self.parity = parity 152 self.header = _header 153 self.extver = extver 154 155 # Create a pointer to the mask file's data array 156 # and the name of the original input DQ file 157 self.maskname = None 158 self.singlemaskname = None 159 self.masklist = None 160 if mask != None: 161 # Specifies filenames to be used if created. 162 self.maskname = mask[0] 163 self.singlemaskname = mask[1] 164 self.masklist = mask[2] 165 166 self.dqname = dqname 167 168 # Remember the name of the coeffs file generated for this chip 169 self.coeffs = self.buildCoeffsName() 170 171 # Read the name of idcfile from image header if not explicitly 172 # provided by user. 173 if idckey != None and idckey.lower() != 'wcs': 174 _indx = expname.find('[') 175 if _indx > -1: 176 _idc_fname = expname[:_indx]+'[0]' 177 else: _idc_fname = expname+'[0]' 178 179 idcfile, idctype = drutil.getIDCFile(self.header,keyword=idckey, 180 directory=idcdir) 181 else: 182 idcfile = None 183 idctype = None 184 185 if (idckey != None) and (idckey.lower() == 'header'): 186 idckey = idctype 187 188 # Get distortion model and WCS info. 189 self.geometry = ObsGeometry(expname, idcfile, idckey=idckey, 190 chip=_chip, new=new, header=self.header, 191 pa_key=pa_key, rot=rot, date=self.dateobs, 192 ref_pscale=ref_pscale, binned=binned, mt_wcs=mt_wcs) 193 194 # Remember the name and type of the IDC file used... 195 self.idcfile = idcfile 196 self.idctype = idctype 197 198 # Remember the names of the filters used for the exposure 199 self.filters = self.geometry.filter1+','+self.geometry.filter2 200 201 # Define shape here... 202 # nx,ny,pixel scale 203 # 204 if wcs != None: 205 # We have been passed a WCS to use 206 self.geometry.wcs = wcs 207 self.geometry.model.pscale = wcs.pscale 208 if expname != None: 209 self.geometry.wcs.rootname = expname 210 211 self.naxis1 = self.geometry.wcs.naxis1 212 self.naxis2 = self.geometry.wcs.naxis2 213 self.pscale = self.geometry.wcs.pscale 214 self.shape = (self.naxis1,self.naxis2,self.pscale) 215 216 # Keep track of the positions of the corners of the exposure 217 # both for the RAW image and the 218 # distortion-corrected, unscaled, unrotated image 219 self.corners = {'raw':np.zeros((4,2),dtype=np.float64),'corrected':np.zeros((4,2),dtype=np.float64)} 220 self.setCorners() 221 222 # Generate BLOT output name specific to this Exposure 223 _blot_extn = '_sci'+repr(extver)+'_blt.fits' 224 self.outblot = fileutil.buildNewRootname(self.name,extn=_blot_extn) 225 226 # Keep track of undistorted frame's position relative to metachip 227 # Zero-point offset for chip relative to meta-chip product 228 # These values get computed using 'setSingleOffsets' from 'writeCoeffs' 229 # to insure that the final XDELTA/YDELTA values have been computed. 230 self.product_wcs = self.geometry.wcslin 231 self.xzero = 0. 232 self.yzero = 0. 233 self.chip_shape = (0.,0.) 234 self.xsh2 = 0. 235 self.ysh2 = 0. 236 237 if _open: 238 handle.close() 239 del handle

240

241 - def set_bunit(self,value=None):

242 """Sets the value of bunit to user specified value, if one is provided. 243 """ 244 if value is not None and self.bunit is not None: 245 self.bunit = value

246

247 - def setCorners(self):

248 """ Initializes corners for the raw image. """ 249 self.corners['raw'] = np.array([(1.,1.),(1.,self.naxis2),(self.naxis1,1.),(self.naxis1,self.naxis2)],dtype=np.float64)

250

251 - def setSingleOffsets(self):

252 """ Computes the zero-point offset and shape of undistorted single chip relative 253 to the full final output product metachip. 254 """ 255 _wcs = self.geometry.wcs 256 _corners = np.array([(1.,1.),(1.,_wcs.naxis2),(_wcs.naxis1,1.),(_wcs.naxis1,_wcs.naxis2)],dtype=np.int64) 257 _wc = self.geometry.wtraxy(_corners,self.product_wcs) 258 259 _xrange = (_wc[:,0].min(),_wc[:,0].max()) 260 _yrange = (_wc[:,1].min(),_wc[:,1].max()) 261 262 self.xzero = int(_xrange[0] - 1) 263 self.yzero = int(_yrange[0] - 1) 264 if self.xzero < 0: self.xzero = 0 265 if self.yzero < 0: self.yzero = 0 266 267 _out_naxis1 = int(ceil(_xrange[1]) - floor(_xrange[0])) 268 _out_naxis2 = int(ceil(_yrange[1]) - floor(_yrange[0])) 269 _max_x = _out_naxis1 + self.xzero 270 _max_y = _out_naxis2 + self.yzero 271 if _max_x > self.product_wcs.naxis1: _out_naxis1 -= (_max_x - self.product_wcs.naxis1) 272 if _max_y > self.product_wcs.naxis2: _out_naxis2 -= (_max_y - self.product_wcs.naxis2) 273 274 self.chip_shape = (_out_naxis1,_out_naxis2) 275 self.xsh2 = int(ceil((self.product_wcs.naxis1 - _out_naxis1)/2.)) - self.xzero 276 self.ysh2 = int(ceil((self.product_wcs.naxis2 - _out_naxis2)/2.)) - self.yzero

277

278 - def getShape(self):

279 """ 280 This method gets the shape after opening and reading 281 the input image. This version will be the default 282 way of doing this, but each instrument may override 283 this method with one specific to their data format. 284 """ 285 return self.shape

286

287 - def setShape(self,size,pscale):

288 """ 289 This method will set the shape for a new file if 290 there is no image header information. 291 292 Size will be defined as (nx,ny) and pixel size 293 """ 294 self.shape = (size[0],size[1],pscale)

295 296

297 - def buildCoeffsName(self):

298 """ Define the name of the coeffs file used for this chip. """ 299 indx = self.name.rfind('.') 300 return self.name[:indx]+'_coeffs'+self.chip+'.dat'

301

302 - def writeCoeffs(self):

303 """ Write out coeffs file for this chip. """ 304 # Pass along the reference position assumed by Drizzle 305 # based on 'align=center' according to the conventions 306 # described in the help page for 'drizzle'. 26Mar03 WJH 307 # 308 # coord for image array reference pixel in full chip coords 309 # 310 if not self.geometry.model.refpix.has_key('empty_model'): 311 _xref = self.geometry.wcs.chip_xref 312 _yref = self.geometry.wcs.chip_yref 313 else: 314 _xref = None 315 _yref = None 316 317 # If we have a subarray, pass along the offset reference position 318 _delta = not (self.geometry.wcslin.subarray) 319 320 # Set up the idcfile for use by 'drizzle' 321 self.geometry.model.convert(self.coeffs,xref=_xref,yref=_yref,delta=_delta) 322 323 # set exposure zero-point, now that all values are computed... 324 self.setSingleOffsets()

325

326 - def getDGEOExtn(self):

327 """ Builds filename with extension to access distortion 328 correction image extension appropriate to each chip. 329 """ 330 # If no DGEOFILE has been given, then simply return blanks 331 # and 'drizzle' will not use any. 332 if not self.dgeoname or self.dgeoname == 'N/A': 333 return '','' 334 335 # Open file for introspection. 336 fimg = fileutil.openImage(self.dgeoname) 337 dx_extver = None 338 dy_extver = None 339 # Find which extensions match this chip ID 340 # We need to identify both a DX and DY EXTNAME extension 341 for hdu in fimg: 342 hdr = hdu.header 343 if not hdr.has_key('CCDCHIP'): 344 _chip = 1 345 else: 346 _chip = int(hdr['CCDCHIP']) 347 348 if hdr.has_key('EXTNAME'): 349 _extname = hdr['EXTNAME'].lower() 350 351 if _chip == int(self.chip): 352 if _extname == 'dx': 353 dx_extver = hdr['EXTVER'] 354 355 if _extname == 'dy': 356 dy_extver = hdr['EXTVER'] 357 358 fimg.close() 359 del fimg 360 361 # Set the name for each extension here... 362 _dxgeo = self.dgeoname+'[DX,'+str(dx_extver)+']' 363 _dygeo = self.dgeoname+'[DY,'+str(dy_extver)+']' 364 365 return _dxgeo,_dygeo

366

367 - def getDGEOArrays(self):

368 """ Return numpy objects for the distortion correction 369 image arrays. 370 371 If no DGEOFILE is specified, it will return 372 empty 2x2 arrays. 373 """ 374 375 # Instantiate array objects for distortion correction image arrays 376 if self.xgeoim == '': 377 # No distortion image specified. 378 # Defaulting to empty 2x2 array. 379 xgdim = ygdim = 2 380 _pxg = np.zeros((ygdim,xgdim),dtype=np.float32) 381 _pyg = np.zeros((ygdim,xgdim),dtype=np.float32) 382 else: 383 # Open distortion correction FITS file 384 _xgfile = fileutil.openImage(self.xgeoim) 385 #_xgfile.info() 386 387 # Access the extensions which correspond to this Exposure 388 _xgname,_xgext = fileutil.parseFilename(self.xgeoim) 389 _ygname,_ygext = fileutil.parseFilename(self.ygeoim) 390 391 _pxgext = fileutil.getExtn(_xgfile,extn=_xgext) 392 _pygext = fileutil.getExtn(_xgfile,extn=_ygext) 393 394 # Copy out the numpy objects for output 395 _ltv1 = int(self.geometry.wcs.offset_x) 396 _ltv2 = int(self.geometry.wcs.offset_y) 397 if _ltv1 != 0. or _ltv2 != 0.: 398 # subarray section only 399 _pxg = _pxgext.data[_ltv2:_ltv2+self.naxis2,_ltv1:_ltv1+self.naxis1].copy() 400 _pyg = _pygext.data[_ltv2:_ltv2+self.naxis2,_ltv1:_ltv1+self.naxis1].copy() 401 else: 402 # full array 403 _pxg = _pxgext.data.copy() 404 _pyg = _pygext.data.copy() 405 406 # Close file handles now... 407 _xgfile.close() 408 del _xgfile 409 410 return _pxg,_pyg

411 412

413 - def applyDeltaWCS(self,dcrval=None,dcrpix=None,drot=None,dscale=None):

414 """ 415 Apply shifts to the WCS of this exposure. 416 417 Shifts are always relative to the current value and in the 418 same reference frame. 419 """ 420 in_wcs = self.geometry.wcs 421 in_wcslin = self.geometry.wcslin 422 423 if dcrpix: 424 _crval = in_wcs.xy2rd((in_wcs.crpix1-dcrpix[0],in_wcs.crpix2-dcrpix[1])) 425 elif dcrval: 426 _crval = (in_wcs.crval1 - dcrval[0],in_wcs.crval2 - dcrval[1]) 427 else: 428 _crval = None 429 430 _orient = None 431 _scale = None 432 if drot: 433 _orient = in_wcs.orient + drot 434 if dscale: 435 _scale = in_wcs.pscale * dscale 436 437 _crpix = (in_wcs.crpix1,in_wcs.crpix2) 438 439 in_wcs.updateWCS(pixel_scale=_scale,orient=_orient,refval=_crval,refpos=_crpix) 440 in_wcslin.updateWCS(pixel_scale=_scale,orient=_orient,refval=_crval,refpos=_crpix)

441 442

443 - def calcNewEdges(self,pscale=None):

444 """ 445 This method will compute arrays for all the pixels around 446 the edge of an image AFTER applying the geometry model. 447 448 Parameter: delta - offset from nominal crpix center position 449 450 Output: xsides - array which contains the new positions for 451 all pixels along the LEFT and RIGHT edges 452 ysides - array which contains the new positions for 453 all pixels along the TOP and BOTTOM edges 454 The new position for each pixel is calculated by calling 455 self.geometry.apply() on each position. 456 """ 457 # build up arrays for pixel positions for the edges 458 # These arrays need to be: array([(x,y),(x1,y1),...]) 459 numpix = self.naxis1*2 + self.naxis2 * 2 460 border = np.zeros(shape=(numpix,2),dtype=np.float64) 461 462 # Now determine the appropriate values for this array 463 # We also need to account for any subarray offsets 464 xmin = 1. 465 xmax = self.naxis1 466 ymin = 1. 467 ymax = self.naxis2 468 469 # Build range of pixel values for each side 470 # Add 1 to make them consistent with pixel numbering in IRAF 471 # Also include the LTV offsets to represent position in full chip 472 # since the model works relative to full chip positions. 473 xside = np.arange(self.naxis1) + xmin 474 yside = np.arange(self.naxis2) + ymin 475 476 #Now apply them to the array to generate the appropriate tuples 477 #bottom 478 _range0 = 0 479 _range1 = self.naxis1 480 border[_range0:_range1,0] = xside 481 border[_range0:_range1,1] = ymin 482 #top 483 _range0 = _range1 484 _range1 = _range0 + self.naxis1 485 border[_range0:_range1,0] = xside 486 border[_range0:_range1,1] = ymax 487 #left 488 _range0 = _range1 489 _range1 = _range0 + self.naxis2 490 border[_range0:_range1,0] = xmin 491 border[_range0:_range1,1] = yside 492 #right 493 _range0 = _range1 494 _range1 = _range0 + self.naxis2 495 border[_range0:_range1,0] = xmax 496 border[_range0:_range1,1] = yside 497 498 # calculate new edge positions 499 border[:,0],border[:,1] = self.geometry.apply(border,pscale=pscale) 500 501 #print 'Calculated corrected border positions at ',_ptime() 502 #Now apply any chip-to-chip offset from REFPIX 503 _refpix = self.geometry.model.refpix 504 if _refpix != None: 505 _ratio = pscale / _refpix['PSCALE'] 506 _delta = (_refpix['XDELTA']/_ratio, _refpix['YDELTA']/_ratio) 507 else: 508 _delta = (0.,0.) 509 510 return border + _delta

511

512 - def getWCS(self):

513 return self.geometry.wcs

514 - def showWCS(self):

515 print self.geometry.wcs

516

517 - def runDriz(self,pixfrac=1.0,kernel='turbo',fillval='INDEF'):

518 """ Runs the 'drizzle' algorithm on this specific chip to create 519 a numpy object of the undistorted image. 520 521 The resulting drizzled image gets returned as a numpy object. 522 """ 523 # 524 # Perform drizzling... 525 # 526 _wcs = self.geometry.wcs 527 _wcsout = self.product_wcs 528 529 # Rigorously compute the orientation changes from WCS 530 # information using algorithm provided by R. Hook from WDRIZZLE. 531 abxt,cdyt = drutil.wcsfit(self.geometry, self.product_wcs) 532 533 # Compute the rotation and shifts between input and reference WCS. 534 xsh = abxt[2] 535 ysh = cdyt[2] 536 rot = fileutil.RADTODEG(np.arctan2(abxt[1],cdyt[0])) 537 scale = self.product_wcs.pscale / self.geometry.wcslin.pscale 538 539 # Now, trim the final output to only match this chip 540 _out_naxis1,_out_naxis2 = self.chip_shape 541 542 # 543 # Insure that the coeffs file was created 544 # 545 if not os.path.exists(self.coeffs): 546 self.writeCoeffs() 547 548 # A image buffer needs to be setup for converting the input 549 # arrays (sci and wht) from FITS format to native format 550 # with respect to byteorder and byteswapping. 551 # This buffer should be reused for each input. 552 # 553 _outsci = np.zeros((_out_naxis2,_out_naxis1),dtype=np.float32) 554 _outwht = np.zeros((_out_naxis2,_out_naxis1),dtype=np.float32) 555 _inwcs = np.zeros([8],dtype=np.float64) 556 557 # Only one chip will ever be drizzled using this method, so 558 # the context image will only ever contain 1 bit-plane 559 _outctx = np.zeros((_out_naxis2,_out_naxis1),dtype=np.int32) 560 561 # Read in the distortion correction arrays, if specifij8cw08n4q_raw.fitsed 562 _pxg,_pyg = self.getDGEOArrays() 563 564 # Open the SCI image 565 _expname = self.name 566 _handle = fileutil.openImage(_expname,mode='readonly',memmap=0) 567 _fname,_extn = fileutil.parseFilename(_expname) 568 _sciext = fileutil.getExtn(_handle,extn=_extn) 569 570 # Make a local copy of SCI array and WCS info 571 #_insci = _sciext.data.copy() 572 _inwcs = drutil.convertWCS(wcsutil.WCSObject(_fname,header=_sciext.header),_inwcs) 573 574 # Compute what plane of the context image this input would 575 # correspond to: 576 _planeid = 1 577 578 # Select which mask needs to be read in for drizzling 579 _inwht = np.ones((self.naxis2,self.naxis1),dtype=np.float32) 580 581 # Default case: wt_scl = exptime 582 _wtscl = self.exptime 583 584 # Set additional parameters needed by 'drizzle' 585 _expin = self.exptime 586 #_in_un = 'counts' 587 #_shift_fr = 'output' 588 #_shift_un = 'output' 589 _uniqid = 1 590 ystart = 0 591 nmiss = 0 592 nskip = 0 593 _vers = '' 594 595 # 596 # This call to 'arrdriz.tdriz' uses the F2C syntax 597 # 598 _dny = self.naxis2 599 # Call 'drizzle' to perform image combination 600 _vers,nmiss,nskip = arrdriz.tdriz(_sciext.data.copy(),_inwht, 601 _outsci, _outwht, _outctx, 602 _uniqid, ystart, 1, 1, _dny, 603 xsh,ysh, 'output','output', rot,scale, 604 self.xsh2,self.ysh2, 1.0, 1.0, 0.0, 'output', 605 _pxg,_pyg, 606 'center', pixfrac, kernel, 607 self.coeffs, 'counts', _expin,_wtscl, 608 fillval, _inwcs, nmiss, nskip, 1,0.,0.) 609 # 610 # End of F2C syntax 611 # 612 613 if nmiss > 0: 614 print '! Warning, ',nmiss,' points were outside the output image.' 615 if nskip > 0: 616 print '! Note, ',nskip,' input lines were skipped completely.' 617 618 # Close image handle 619 _handle.close() 620 del _handle,_fname,_extn,_sciext 621 del _inwht 622 623 del _pxg,_pyg 624 625 del _outwht,_outctx 626 627 return _outsci

Source Code for Module pydrizzle.exposure