| Home | Trees | Indices | Help |
|
|---|
|
|
1 from __future__ import division # confidence medium
2 from drutil import DEFAULT_IDCDIR
3 from stsci.tools import fileutil, wcsutil, asnutil
4
5 import string, os, types, sys
6 import shutil
7 import arrdriz
8 import outputimage
9 from pattern import *
10 from observations import *
11
12 # Hackish, but not invalid
13 from __init__ import __version__
14 #import buildasn
15 #import p_m_input
16
17 import pyfits
18 import numpy as np
19
20 yes = True # 1
21 no = False # 0
22 # default_pars is passed as a default parameter dictionary in selectInstrument.
23 # Looks like this is not needed any more. 'default_rot' was originally 'rot'.
24 _default_pars = {'psize':None,'default_rot':None,'idckey':None}
25
26 INSTRUMENT = ["ACS","WFPC2","STIS","NICMOS","WFC3"]
27
28
30 """
31 Program to process and/or dither-combine image(s) using (t)drizzle.
32 To create an object named 'test' that corresponds to a drizzle product:
33 --> test = pydrizzle.PyDrizzle(input)
34 where input is the FULL filename of an ACS observation or ASN table.
35 This computes all the parameters necessary for running drizzle on all
36 the input images. Once this object is created, you can run drizzle using:
37 --> test.run()
38
39 The 'clean()' method can be used to remove files which would interfere with
40 running Drizzle again using the 'run()' method after a product has already
41 been created.
42
43 Optional parameters:
44 output User-specified name for output products
45 field User-specified parameters for output image
46 includes: psize, orient, ra, dec, shape
47 units Units for final product: 'counts' or 'cps'(DEFAULT)
48 section Extension/group to be drizzled: list or single FITS extension(s)
49 or group(s) syntax ('1' or 'sci,1') or None (DEFAULT: Use all chips).
50 kernel Specify which kernel to use in TDRIZZLE
51 'square'(default),'point','gaussian','turbo','tophat'
52 pixfrac drizzle pixfrac value (Default: 1.0)
53 idckey User-specified keyword for determining IDCTAB filename
54 'IDCTAB'(ACS default),'TRAUGER'(WFPC2),'CUBIC'(WFPC2)
55 idcdir User-specified directory for finding coeffs files:
56 'drizzle$coeffs' (default)
57 bits_single Specify DQ values to be considered good when
58 drizzling with 'single=yes'. (Default: 0)
59 bits_final Specify DQ values to be considered good when
60 drizzling with 'single=no'. (Default: 0)
61 updatewcs Flag whether to run makewcs on the input (Default: True)
62 Bits parameters will be interpreted as:
63 None - No DQ information to be used, no mask created
64 Int - sum of DQ values to be considered good
65
66 Optional Parameters for '.run()':
67 build create multi-extension output: yes (Default) or no
68 save keeps the individual inputs from drizzle: yes or no (Default)
69 single drizzle to separate output for each input: yes or no (Default)
70 blot run blot on drizzled products: yes or no (Default)
71 clean remove coeffs and static mask files: yes or no (Default)
72
73 Optional Parameters for '.clean()':
74 coeffs Removes coeffs and static mask files: yes or no (Default)
75 final Removes final product: yes or no (Default)
76
77 Usage of optional parameters:
78 --> test = pydrizzle.PyDrizzle('test_asn.fits',units='counts')
79 To keep the individual 'drizzle' output products:
80 --> test.run(save=yes)
81
82 Output frame parameters can be modified 'on-the-fly' using 'resetPars'.
83 Given an already drizzled image 'refimg_drz.fits' as a reference,
84 reset drizzle parameters using:
85 --> wcsref = wcsutil.WCSObject('refimg_drz.fits[sci,1]')
86 --> f = pydrizzle.SkyField(wcs=wcsref)
87 Use either:
88 --> test.resetPars(wcsref)
89 Or:
90 --> test.resetPars(f)
91 Return to default parameters using no parameters at all:
92 --> test.resetPars()
93 More help on SkyField objects and their parameters can be obtained using:
94 --> f.help()
95 """
96
97 - def __init__(self, input, output=None, field=None, units=None, section=None,
98 kernel=None,pixfrac=None,bits_final=0,bits_single=0,
99 wt_scl='exptime', fillval=0.,idckey='', in_units='counts',
100 idcdir=DEFAULT_IDCDIR,memmap=0,dqsuffix=None):
101
102 if idcdir == None: idcdir = DEFAULT_IDCDIR
103
104
105 print 'Starting PyDrizzle Version ',__version__,' at ', _ptime()
106
107
108 self.input = input
109 self.output = output
110 asndict = input
111
112 # Insure that the section parameter always becomes a list
113 if isinstance(section,list) == False and section != None:
114 section = [section]
115
116 # Set the default value for 'build'
117 self.build = yes
118
119 # Extract user-specified parameters, if any have been set...
120 # 'field' will be a SkyField object...
121 if field != None:
122 psize = field.psize
123 orient = field.orient
124 orient_rot = orient
125 else:
126 psize = None
127 orient = None
128 orient_rot = 0.
129
130 #self.pars['orient_rot'] was originaly self.pars['rot']. It is set based on
131 # orientat but it looks like it is not used. In SkyField.mergeWCS() rot is
132 #defined again based on orientat.
133
134 # These can also be set by the user.
135 # Minimum set needed: psize, rot, and idckey
136 self.pars = {'psize':psize,'units':units,'kernel':kernel,'orient_rot':orient_rot,
137 'pixfrac':pixfrac,'idckey':idckey,'wt_scl':wt_scl,
138 'fillval':fillval,'section':section, 'idcdir':idcdir+os.sep,
139 'memmap':memmap,'dqsuffix':dqsuffix, 'in_units':in_units,
140 'bits':[bits_final,bits_single], 'mt_wcs': None}
141
142 # Watch out for any errors.
143 # If they arise, all open files need to be closed...
144 self.observation = None
145 self.parlist = None
146
147 if len(asndict['order']) > 1:
148 self.observation = DitherProduct(asndict,pars=self.pars)
149 else:
150 inroot = asndict['order'][0]
151 pardict = asndict['members'][inroot]
152 infile = fileutil.buildRootname(inroot)
153 if infile == None:
154 raise IOError,'No product found for association table.'
155 # Append user specified parameters to shifts dictionary
156 pardict.update(self.pars)
157 #self.observation = selectInstrument(infile,output,pars=pardict)
158 self.observation = selectInstrument(infile,self.output,pars=pardict)
159 """
160 if self.output == None:
161 self.output = fileutil.buildNewRootname(asndict['output'],extn='_drz.fits')
162 print 'Setting up output name: ',output
163 else:
164 self.output = output
165 """
166 # This call puts together the parameters for the input image
167 # with those for the output to create a final parameter list
168 # for running 'drizzle'.
169 # It relies on the buildPars() methods for each exposure to
170 # generate a complete set of parameters for all inputs
171 #
172 self.translateShifts(asndict)
173 self.parlist = self.observation.buildPars(ref=field)
174
175 # Something went wrong, so we need to make sure all file
176 # handles get closed...
177 #self.close()
178 #print 'PyDrizzle could not initialize due to errors.'
179 self.observation.closeHandle()
180
181 # Let the user know parameters have been successfully calculated
182 print 'Drizzle parameters have been calculated. Ready to .run()...'
183 print 'Finished calculating parameters at ',_ptime()
184
186 """
187 Translate the shifts specified in the ASNDICT (as read in from the
188 shiftfile) into offsets in the sky, so they can be translated back
189 into the WCS of the PyDrizzle output product.
190
191 NOTE: Only works with 'delta' shifts now, and
192 requires that a 'refimage' be specified.
193 """
194
195 # for each set of shifts, translate them into delta(ra,dec) based on refwcs
196 for img in asndict['order']:
197
198 xsh = asndict['members'][img]['xshift']
199 ysh = asndict['members'][img]['yshift']
200
201 if xsh == 0.0 and ysh == 0.0:
202 delta_ra = 0.0
203 delta_dec = 0.0
204 else:
205 #check the units for the shifts...
206 if asndict['members'][img]['shift_units'] == 'pixels':
207 # Initialize the reference WCS for use in translation
208 # NOTE: This assumes that a 'refimage' has been specified for
209 # every set of shifts.
210 refwcs = wcsutil.WCSObject(asndict['members'][img]['refimage'])
211 cp1 = refwcs.crpix1
212 cp2 = refwcs.crpix2
213
214 nra,ndec = refwcs.xy2rd((cp1+xsh,cp2+ysh))
215
216 delta_ra = refwcs.crval1-nra
217 delta_dec = refwcs.crval2-ndec
218 else:
219 # Shifts already in units of RA/Dec (decimal degrees)
220 # No conversion necessary
221 delta_ra = xsh
222 delta_dec = ysh
223
224 asndict['members'][img]['delta_ra'] = delta_ra
225 asndict['members'][img]['delta_dec'] = delta_dec
226
228 """ Removes intermediate products from disk. """
229 for img in self.parlist:
230 # If build == no, then we do not want to delete the only
231 # products created by PyDrizzle; namely,
232 # outdata, outcontext, outweight
233 if self.build == yes:
234 fileutil.removeFile([img['outdata'],img['outcontext'],img['outweight']])
235
236 fileutil.removeFile([img['outsingle'],img['outsweight']])
237 #fileutil.removeFile([img['outsingle'],img['outsweight'],img['outscontext']])
238 fileutil.removeFile(img['outblot'])
239 if coeffs:
240 os.remove(img['coeffs'])
241 if img['driz_mask'] != None:
242 fileutil.removeFile(img['driz_mask'])
243 if img['single_driz_mask'] != None:
244 fileutil.removeFile(img['single_driz_mask'])
245 if final:
246 fileutil.removeFile(self.output)
247
248
249 # Run 'drizzle' here...
250 #
251 - def run(self,save=no,build=yes,blot=no,single=no,clean=no,interp='linear',sinscl=1.0, debug=no):
252 """Perform drizzle operation on input to create output.
253 This method would rely on the buildPars() method for
254 the output product to produce correct parameters
255 based on the inputs. The output for buildPars() is a list
256 of dictionaries, one for each input, that matches the
257 primary parameters for an IRAF drizzle task.
258
259 This method would then loop over all the entries in the
260 list and run 'drizzle' for each entry. """
261
262 # Store the value of build set by the user for use, if desired,
263 # in the 'clean()' method.
264 self.build = build
265 self.debug = debug
266
267 print 'PyDrizzle: drizzle task started at ',_ptime()
268 _memmap = self.pars['memmap']
269
270 # Check for existance of output file.
271 if single == no and build == yes and fileutil.findFile(self.output):
272 print 'Removing previous output product...'
273 os.remove(self.output)
274 #
275 # Setup the versions info dictionary for output to PRIMARY header
276 # The keys will be used as the name reported in the header, as-is
277 #
278 _versions = {'PyDrizzle':__version__,'PyFITS':pyfits.__version__,'Numpy':np.__version__}
279
280 # Set parameters for each input and run drizzle on it here.
281
282 if blot:
283 #
284 # Run blot on data...
285 #
286
287 _hdrlist = []
288
289 for plist in self.parlist:
290
291 _insci = np.zeros((plist['outny'],plist['outnx']),dtype=np.float32)
292 _outsci = np.zeros((plist['blotny'],plist['blotnx']),dtype=np.float32)
293 _hdrlist.append(plist)
294 # Open input image as PyFITS object
295 if plist['outsingle'] != plist['outdata']:
296 _data = plist['outsingle']
297 else:
298 _data = plist['outdata']
299
300 # PyFITS can be used here as it will always operate on
301 # output from PyDrizzle (which will always be a FITS file)
302 # Open the input science file
303 _fname,_sciextn = fileutil.parseFilename(_data)
304 _inimg = fileutil.openImage(_fname)
305
306 # Return the PyFITS HDU corresponding to the named extension
307 _scihdu = fileutil.getExtn(_inimg,_sciextn)
308 _insci = _scihdu.data.copy()
309
310 # Read in the distortion correction arrays, if specified
311 _pxg,_pyg = plist['exposure'].getDGEOArrays()
312
313 # Now pass numpy objects to callable version of Blot...
314 #runBlot(plist)
315 build=no
316 misval = 0.0
317 kscale = 1.0
318
319 xmin = 1
320 xmax = plist['outnx']
321 ymin = 1
322 ymax = plist['outny']
323 #
324 # Convert shifts to input units
325 #
326 xsh = plist['xsh'] * plist['scale']
327 ysh = plist['ysh'] * plist['scale']
328 # ARRDRIZ.TBLOT needs to be updated to support 'poly5' interpolation,
329 # and exptime scaling of output image.
330 #
331 """
332 #
333 # This call to 'arrdriz.tdriz' uses the F2PY syntax
334 #
335 arrdriz.tblot(_insci, _outsci,xmin,xmax,ymin,ymax,
336 plist['xsh'],plist['ysh'],
337 plist['rot'],plist['scale'], kscale, _pxg, _pyg,
338 'center',interp, plist['coeffs'], plist['exptime'],
339 misval, sinscl, 1)
340 #
341 # End of F2PY syntax
342 #
343 """
344 #
345 # This call to 'arrdriz.tdriz' uses the F2C syntax
346 #
347 if (_insci.dtype > np.float32):
348 #WARNING: Input array recast as a float32 array
349 _insci = _insci.astype(np.float32)
350 t = arrdriz.tblot(_insci, _outsci,xmin,xmax,ymin,ymax,
351 xsh,ysh, plist['rot'],plist['scale'], kscale,
352 0.0,0.0, 1.0,1.0, 0.0, 'output',
353 _pxg, _pyg,
354 'center',interp, plist['coeffs'], plist['exptime'],
355 misval, sinscl, 1,0.0,0.0)
356
357 #
358 # End of F2C syntax
359 #
360
361 # Write output Numpy objects to a PyFITS file
362 # Blotting only occurs from a drizzled SCI extension
363 # to a blotted SCI extension...
364 _header = fileutil.getHeader(plist['data'])
365 _wcs = wcsutil.WCSObject(plist['data'],header=_header)
366
367 _outimg = outputimage.OutputImage(_hdrlist, build=no, wcs=_wcs, blot=yes)
368 _outimg.outweight = None
369 _outimg.outcontext = None
370 _outimg.writeFITS(plist['data'],_outsci,None,versions=_versions)
371
372 #_buildOutputFits(_outsci,None,plist['outblot'])
373 _insci *= 0.
374 _outsci *= 0.
375 _inimg.close()
376 del _inimg
377 _hdrlist = []
378
379 del _pxg,_pyg
380
381 del _insci,_outsci
382 del _outimg
383
384 else:
385 #
386 # Perform drizzling...
387 #
388 # Only work on a copy of the product WCS, so that when
389 # this gets updated for the output image, it does not
390 # modify the original WCS computed by PyDrizzle
391 _wcs = self.observation.product.geometry.wcs.copy()
392
393 _numctx = {'all':len(self.parlist)}
394 # if single:
395 # Determine how many chips make up each single image
396 for plist in self.parlist:
397 plsingle = plist['outsingle']
398 if _numctx.has_key(plsingle): _numctx[plsingle] += 1
399 else: _numctx[plsingle] = 1
400 #
401 # A image buffer needs to be setup for converting the input
402 # arrays (sci and wht) from FITS format to native format
403 # with respect to byteorder and byteswapping.
404 # This buffer should be reused for each input.
405 #
406 plist = self.parlist[0]
407 _outsci = np.zeros((plist['outny'],plist['outnx']),dtype=np.float32)
408 _outwht = np.zeros((plist['outny'],plist['outnx']),dtype=np.float32)
409 _inwcs = np.zeros([8],dtype=np.float64)
410
411 # Compute how many planes will be needed for the context image.
412 _nplanes = int((_numctx['all']-1) / 32) + 1
413 # For single drizzling or when context is turned off,
414 # minimize to 1 plane only...
415 if single or self.parlist[0]['outcontext'] == '':
416 _nplanes = 1
417
418 # Always initialize context images to a 3-D array
419 # and only pass the appropriate plane to drizzle as needed
420 _outctx = np.zeros((_nplanes,plist['outny'],plist['outnx']),dtype=np.int32)
421
422 # Keep track of how many chips have been processed
423 # For single case, this will determine when to close
424 # one product and open the next.
425 _numchips = 0
426 _nimg = 0
427 _hdrlist = []
428
429 for plist in self.parlist:
430 # Read in the distortion correction arrays, if specifij8cw08n4q_raw.fitsed
431 _pxg,_pyg = plist['exposure'].getDGEOArrays()
432
433 _hdrlist.append(plist)
434 # Open the SCI image
435 _expname = plist['data']
436 _handle = fileutil.openImage(_expname,mode='readonly',memmap=0)
437 _fname,_extn = fileutil.parseFilename(_expname)
438 _sciext = fileutil.getExtn(_handle,extn=_extn)
439
440 # Make a local copy of SCI array and WCS info
441 #_insci = _sciext.data.copy()
442 _inwcs = drutil.convertWCS(wcsutil.WCSObject(_fname,header=_sciext.header),_inwcs)
443
444 # Determine output value of BUNITS
445 # and make sure it is not specified as 'ergs/cm...'
446 # we want to use the actual value from the input image header directly
447 # when possible in order to account for any unit conversions that may
448 # be applied to the input image between initialization of PyDrizzle
449 # and the calling of this run() method.
450 if _sciext.header.has_key('bunit') and _sciext.header['bunit'] not in ['','N/A']:
451 _bunit = _sciext.header['bunit']
452 else:
453 # default based on instrument-specific logic
454 _bunit = plist['bunit']
455
456 _bindx = _bunit.find('/')
457 if plist['units'] == 'cps':
458 # If BUNIT value does not specify count rate already...
459 if _bindx < 1:
460 # ... append '/SEC' to value
461 _bunit += '/S'
462 else:
463 # reset _bunit here to None so it does not
464 # overwrite what is already in header
465 _bunit = None
466 else:
467 if _bindx > 0:
468 # remove '/S'
469 _bunit = _bunit[:_bindx]
470 else:
471 # reset _bunit here to None so it does not
472 # overwrite what is already in header
473 _bunit = None
474
475 # Compute what plane of the context image this input would
476 # correspond to:
477 _planeid = int(_numchips /32)
478
479 # Select which mask needs to be read in for drizzling
480 if single:
481 _mask = plist['single_driz_mask']
482 else:
483 _mask = plist['driz_mask']
484
485 # Check to see whether there is a mask_array at all to use...
486 if isinstance(_mask,types.StringType):
487 if _mask != None and _mask != '':
488 _wht_handle = fileutil.openImage(_mask,mode='readonly',memmap=0)
489 _inwht = _wht_handle[0].data.astype(np.float32)
490 _wht_handle.close()
491 del _wht_handle
492 else:
493 print 'No weight or mask file specified! Assuming all pixels are good.'
494 _inwht = np.ones((plist['blotny'],plist['blotnx']),dtype=np.float32)
495 elif _mask != None:
496 _inwht = _mask.astype(np.float32)
497 else:
498 print 'No weight or mask file specified! Assuming all pixels are good.'
499 _inwht = np.ones((plist['blotny'],plist['blotnx']),dtype=np.float32)
500
501 if plist['wt_scl'] != None:
502 if isinstance(plist['wt_scl'],types.StringType):
503 if plist['wt_scl'].isdigit() == False :
504 # String passed in as value, check for 'exptime' or 'expsq'
505 _wtscl_float = None
506 try:
507 _wtscl_float = float(plist['wt_scl'])
508 except ValueError:
509 _wtscl_float = None
510 if _wtscl_float != None:
511 _wtscl = _wtscl_float
512 elif plist['wt_scl'] == 'expsq':
513 _wtscl = plist['exptime']*plist['exptime']
514 else:
515 # Default to the case of 'exptime', if
516 # not explicitly specified as 'expsq'
517 _wtscl = plist['exptime']
518 else:
519 # int value passed in as a string, convert to float
520 _wtscl = float(plist['wt_scl'])
521 else:
522 # We have a non-string value passed in...
523 _wtscl = float(plist['wt_scl'])
524 else:
525 # Default case: wt_scl = exptime
526 _wtscl = plist['exptime']
527
528 #print 'WT_SCL: ',plist['wt_scl'],' _wtscl: ',_wtscl
529 # Set additional parameters needed by 'drizzle'
530 _in_units = plist['in_units']
531 if _in_units == 'cps':
532 _expin = 1.0
533 else:
534 _expin = plist['exptime']
535 _shift_fr = 'output'
536 _shift_un = 'output'
537 _uniqid = _numchips + 1
538 ystart = 0
539 nmiss = 0
540 nskip = 0
541 _vers = plist['driz_version']
542
543 _con = yes
544 _imgctx = _numctx['all']
545 if single:
546 _imgctx = _numctx[plist['outsingle']]
547 #if single or (plist['outcontext'] == '' and single == yes):
548 if _nplanes == 1:
549 _con = no
550 # We need to reset what gets passed to TDRIZ
551 # when only 1 context image plane gets generated
552 # to prevent overflow problems with trying to access
553 # planes that weren't created for large numbers of inputs.
554 _planeid = 0
555 _uniqid = ((_uniqid-1) % 32) + 1
556
557 """
558 #
559 # This call to 'arrdriz.tdriz' uses the F2PY syntax
560 #
561 #_dny = plist['blotny']
562 # Call 'drizzle' to perform image combination
563 tdriz,nmiss,nskip,_vers = arrdriz.tdriz(
564 _sciext.data,_inwht, _outsci, _outwht,
565 _outctx[_planeid], _con, _uniqid, ystart, 1, 1,
566 plist['xsh'],plist['ysh'], 'output','output',
567 plist['rot'],plist['scale'], _pxg,_pyg,
568 'center', plist['pixfrac'], plist['kernel'],
569 plist['coeffs'], 'counts', _expin,_wtscl,
570 plist['fillval'], _inwcs, 1, nmiss, nskip,_vers)
571 #
572 # End of F2PY syntax
573 #
574 """
575 #
576 # This call to 'arrdriz.tdriz' uses the F2C syntax
577 #
578 _dny = plist['blotny']
579 # Call 'drizzle' to perform image combination
580 if (_sciext.data.dtype > np.float32):
581 #WARNING: Input array recast as a float32 array
582 _sciext.data = _sciext.data.astype(np.float32)
583
584 _vers,nmiss,nskip = arrdriz.tdriz(_sciext.data,_inwht, _outsci, _outwht,
585 _outctx[_planeid], _uniqid, ystart, 1, 1, _dny,
586 plist['xsh'],plist['ysh'], 'output','output',
587 plist['rot'],plist['scale'],
588 0.0,0.0, 1.0,1.0,0.0,'output',
589 _pxg,_pyg, 'center', plist['pixfrac'], plist['kernel'],
590 plist['coeffs'], _in_units, _expin,_wtscl,
591 plist['fillval'], _inwcs, nmiss, nskip, 1,0.0,0.0)
592 """
593 _vers,nmiss,nskip = arrdriz.tdriz(_sciext.data,_inwht, _outsci, _outwht,
594 _outctx[_planeid], _uniqid, ystart, 1, 1, _dny,
595 plist['xsh'],plist['ysh'], 'output','output',
596 plist['rot'],plist['scale'],
597 _pxg,_pyg, 'center', plist['pixfrac'], plist['kernel'],
598 plist['coeffs'], 'counts', _expin,_wtscl,
599 plist['fillval'], _inwcs, nmiss, nskip, 1)
600 """
601 #
602 # End of F2C syntax
603 #
604
605 plist['driz_version'] = _vers
606
607 if nmiss > 0:
608 print '! Warning, ',nmiss,' points were outside the output image.'
609 if nskip > 0:
610 print '! Note, ',nskip,' input lines were skipped completely.'
611 # Close image handle
612 _handle.close()
613 del _handle,_fname,_extn,_sciext
614 del _inwht
615
616 del _pxg,_pyg
617
618 # Remember the name of the first image that goes into
619 # this particular product
620 # This will insure that the header reports the proper
621 # values for the start of the exposure time used to make
622 # this product; in particular, TIME-OBS and DATE-OBS.
623 if _numchips == 0:
624 _template = plist['data']
625
626 if _nimg == 0 and self.debug == yes:
627 # Only update the WCS from drizzling the
628 # first image in the list, just like IRAF DRIZZLE
629 drutil.updateWCS(_inwcs,_wcs)
630 #print '[run] Updated WCS now:'
631 #print _wcs
632
633 # Increment number of chips processed for single output
634 _numchips += 1
635 if _numchips == _imgctx:
636 ###########################
637 #
638 # IMPLEMENTATION REQUIREMENT:
639 #
640 # Need to implement scaling of the output image
641 # from 'cps' to 'counts' in the case where 'units'
642 # was set to 'counts'... 21-Mar-2005
643 #
644 ###########################
645 # Start by determining what exposure time needs to be used
646 # to rescale the product.
647 if single:
648 _expscale = plist['exptime']
649 else:
650 _expscale = plist['texptime']
651
652 #If output units were set to 'counts', rescale the array in-place
653 if plist['units'] == 'counts':
654 np.multiply(_outsci, _expscale, _outsci)
655
656 #
657 # Write output arrays to FITS file(s) and reset chip counter
658 #
659 _outimg = outputimage.OutputImage(_hdrlist, build=build, wcs=_wcs, single=single)
660 _outimg.set_bunit(_bunit)
661 _outimg.set_units(plist['units'])
662
663 _outimg.writeFITS(_template,_outsci,_outwht,ctxarr=_outctx,versions=_versions)
664 del _outimg
665 #
666 # Reset chip counter for next output image...
667 #
668 _numchips = 0
669 _nimg = 0
670 np.multiply(_outsci,0.,_outsci)
671 np.multiply(_outwht,0.,_outwht)
672 np.multiply(_outctx,0,_outctx)
673
674 _hdrlist = []
675 else:
676 _nimg += 1
677
678 del _outsci,_outwht,_inwcs,_outctx, _hdrlist
679
680 # Remove temp files, if desired
681 # Files to be removed are:
682 # parlist['coeffs']
683 # parlist['driz_mask']
684 if not save and clean:
685 for img in self.parlist:
686 fileutil.removeFile(img['coeffs'])
687 if img['driz_mask'] != None:
688 fileutil.removeFile(img['driz_mask'])
689 if img['single_driz_mask'] != None:
690 fileutil.removeFile(img['single_driz_mask'])
691
692 print 'PyDrizzle drizzling completed at ',_ptime()
693
694
696 """
697 Recompute the output parameters based on a new
698 SkyField or WCSObject object.
699 """
700 if field and not isinstance(field, SkyField):
701 if isinstance(field, wcsutil.WCSObject):
702 _ref = SkyField(wcs=field)
703 else:
704 raise TypeError, 'No valid WCSObject or SkyField object entered...'
705 else:
706 _ref = field
707 # Create new version of the parlist with the new values of the
708 # parameters based on the reference image.
709 new_parlist = self.observation.buildPars(ref=_ref)
710
711
712 # Copy the parameters from the new parlist into the existing
713 # parlist (self.parlist) to preserve any changes/updates made
714 # to the parlist externally to PyDrizzle.
715 for i in xrange(len(self.parlist)):
716 for key in new_parlist[i]:
717 self.parlist[i][key] = new_parlist[i][key]
718 del new_parlist
719
720 if pixfrac or kernel or units:
721 #print 'resetting additional parameter(s)...'
722 for p in self.parlist:
723 if kernel:
724 p['kernel'] = kernel
725 if pixfrac:
726 p['pixfrac'] = pixfrac
727 if units:
728 if units == 'counts' or units == 'cps':
729 p['units'] = units
730 else:
731 print 'Units ',units,' not valid! Parameter not reset.'
732 print 'Please use either "cps" or "counts".'
733
739
741 """ Prints common parameters for review. """
742 if format:
743 _title = pars.replace(',',' ')
744 print _title
745 print '-'*72
746
747 _pars = pars.split(',')
748 for pl in self.parlist:
749 for _p in _pars: print pl[_p],
750 print ''
751
753 """ Returns the class instance for the specified member name."""
754 return self.observation.getMember(memname)
755
756
757
758
760 import time
761
762 # Format time values for keywords IRAF-TLM, and DATE
763 _ltime = time.localtime(time.time())
764 tlm_str = time.strftime('%H:%M:%S (%d/%m/%Y)',_ltime)
765 #date_str = time.strftime('%Y-%m-%dT%H:%M:%S',_ltime)
766 return tlm_str
767
768 #################
769
770
771
773 """
774 An class for specifying the parameters and building a WCS object
775 for a user-specified drizzle product.
776 The user may optionally modify the values for:
777 psize - size of image's pixels in arc-seconds
778 orient - value of ORIENTAT for the field
779 shape - tuple containing the sizes of the field's x/y axes
780 ra,dec - position of the center of the field
781 decimal (124.5678) or
782 sexagesimal string _in quotes_ ('hh:mm:ss.ss')
783 crpix - tuple for pixel position of reference point in
784 output image
785 Usage:
786 To specify a new field with a fixed output size of 1024x1024:
787 --> field = pydrizzle.SkyField(shape=(1024,1024))
788
789 The 'set()' method modifies one of the parameters listed above
790 without affecting the remainder of the parameters.
791 --> field.set(psize=0.1,orient=0.0)
792 --> field.set(ra=123.45678,dec=0.1000,crpix=(521,576))
793
794 View the WCS or user-specified values for this object:
795 --> print field
796
797 """
799
800 self.shape = shape
801 self.ra = None
802 self.dec = None
803 self.orient = None
804 self.psize = psize
805
806 self.crpix = None
807
808 # Set up proper shape tuple for WCSObject
809 if shape != None and psize != None:
810 wshape = (shape[0],shape[1],psize)
811 else:
812 wshape = None
813 if wcs:
814 self.crpix = (wcs.crpix1,wcs.crpix2)
815
816 # Specify 'new=yes' with given rootname to unambiguously create
817 # a new WCS from scratch.
818 if wcs == None:
819 self.wcs = wcsutil.WCSObject("New",new=yes)
820 else:
821 self.wcs = wcs.copy()
822 # Need to adjust WCS to match the 'drizzle' task
823 # align=center conventions. WJH 29-Aug-2005
824 #self.wcs.crpix1 -= 1.0
825 #self.wcs.crpix2 -= 1.0
826 self.wcs.recenter()
827
828 self.wcs.updateWCS(size=shape,pixel_scale=psize)
829
830 # Set this to keep track of total exposure time for
831 # output frame... Not user settable.
832 self.exptime = None
833 # Keep track of whether this is for a Dither product or not
834 self.dither = None
835
837 """ Sets up the WCS for this object based on another WCS.
838 This method will NOT update object attributes other
839 than WCS, as all other attributes reflect user-settings.
840 """
841 #
842 # Start by making a copy of the input WCS...
843 #
844 if self.wcs.rootname == 'New':
845 self.wcs = wcs.copy()
846 else:
847 return
848 self.wcs.recenter()
849
850 if self.ra == None:
851 _crval = None
852 else:
853 _crval = (self.ra,self.dec)
854
855 if self.psize == None:
856 _ratio = 1.0
857 _psize = None
858 # Need to resize the WCS for any changes in pscale
859 else:
860 _ratio = wcs.pscale / self.psize
861 _psize = self.psize
862
863 if self.orient == None:
864 _orient = None
865 _delta_rot = 0.
866 else:
867 _orient = self.orient
868 _delta_rot = wcs.orient - self.orient
869
870 _mrot = fileutil.buildRotMatrix(_delta_rot)
871
872 if self.shape == None:
873 _corners = np.array([[0.,0.],[wcs.naxis1,0.],[0.,wcs.naxis2],[wcs.naxis1,wcs.naxis2]])
874 _corners -= (wcs.naxis1/2.,wcs.naxis2/2.)
875 _range = drutil.getRotatedSize(_corners,_delta_rot)
876 shape = ((_range[0][1] - _range[0][0])*_ratio,(_range[1][1]-_range[1][0])*_ratio)
877 old_shape = (wcs.naxis1*_ratio,wcs.naxis2*_ratio)
878
879 _cen = (shape[0]/2., shape[1]/2.)
880
881 #if _delta_rot == 0.:
882 # _crpix = (self.wcs.crpix1,self.wcs.crpix2)
883 #else:
884 # Rotate original scaled crpix position to new orientation
885 #_crpix = N.dot((wcs.crpix1*_ratio - _cen[0],wcs.crpix2*_ratio -_cen[1]),_mrot)+_cen
886 _crpix = _cen
887 else:
888 shape = self.shape
889 if self.crpix == None:
890 _crpix = (self.shape[0]/2.,self.shape[1]/2.)
891 else:
892 _crpix = self.crpix
893
894 # Set up the new WCS based on values from old one.
895 self.wcs.updateWCS(pixel_scale=_psize,orient=_orient,refpos=_crpix,refval=_crval)
896 self.wcs.naxis1 = int(shape[0])
897 self.wcs.naxis2 = int(shape[1])
898
901 """
902 Modifies the attributes of the SkyField and
903 updates it's WCS when appropriate.
904 """
905 # Converts(if necessary), then updates the RA and Dec.
906 _ra,_dec = None,None
907 if ra != None:
908 if string.find(repr(ra),':') > 0:
909 _hms = string.split(repr(ra)[1:-1],':')
910 if _hms[0][0] == '-': _sign = -1
911 else: _sign = 1
912
913 for i in range(len(_hms)): _hms[i] = float(_hms[i])
914 _ra = _sign * (_hms[0] + ((_hms[1] + _hms[2]/60.) / 60.)) * 15.
915 else:
916 _ra = float(ra)
917 self.ra = _ra
918
919 if dec != None:
920 if string.find(repr(dec),':') > 0:
921 _dms = string.split(repr(dec)[1:-1],':')
922 if _dms[0][0] == '-': _sign = -1
923 else: _sign = 1
924
925 for i in range(len(_dms)): _dms[i] = float(_dms[i])
926 _dec = _sign * (_dms[0] + ((_dms[1] + _dms[2]/60.) / 60.))
927 else:
928 _dec = float(dec)
929 self.dec = _dec
930
931 if self.ra != None and self.dec != None:
932 _crval = (self.ra,self.dec)
933 else:
934 _crval = None
935
936 # Updates the shape, and reference position,
937 # only if a new value is specified.
938 _crpix = None
939 if crpix == None:
940 if shape != None:
941 self.shape = shape
942 _crpix = (self.shape[0]/2.,self.shape[1]/2.)
943 else:
944 _crpix = crpix
945
946 self.crpix=_crpix
947
948 if psize != None:
949 self.psize = psize
950
951 if orient != None:
952 self.orient = orient
953
954 # Updates the WCS with all the changes, if there is enough info.
955 self.wcs.updateWCS(pixel_scale=psize,orient=orient,refpos=_crpix,
956 refval=_crval,size=self.shape)
957
959 """ Prints the WCS information set for this object.
960 """
961 if self.psize != None and self.orient != None:
962 block = self.wcs.__str__()
963 else:
964 block = 'User parameters for SkyField object: \n'
965 block = block + ' ra = '+repr(self.ra)+' \n'
966 block = block + ' dec = '+repr(self.dec)+' \n'
967 block = block + ' shape = '+repr(self.shape)+' \n'
968 block = block + ' crpix = '+repr(self.crpix)+' \n'
969 block = block + ' psize = '+repr(self.psize)+' \n'
970 block = block + ' orient = '+repr(self.orient)+' \n'
971 block = block + ' No WCS.\n'
972
973 return block
974
979
981 """
982 Builds an object for a set of dithered inputs, each of which
983 will be one of the Observation objects.
984 """
986 # Build temporary output drizzle product name
987 if prodlist['output'].find('.fits') < 0:
988 if prodlist['output'].rfind('_drz') < 0:
989 output = fileutil.buildNewRootname(prodlist['output'],extn='_drz.fits')
990 else:
991 output = prodlist['output']+'.fits'
992 else:
993 output = prodlist['output']
994
995 # Setup a default exposure to contain the results
996 Pattern.__init__(self, None, output=output, pars=pars)
997 self.pars = prodlist['members']
998
999 self.nmembers = self.nimages = len(prodlist['members'])
1000 self.offsets = None
1001
1002 self.addMembers(prodlist,pars,output)
1003
1004 if len(self.members) == 0:
1005 print 'No suitable inputs from ASN table. Quitting PyDrizzle...'
1006 raise Exception
1007
1008 self.exptime = self.getExptime()
1009
1010 self.buildProduct(output)
1011
1013 """ Close image handle for each member."""
1014 for member in self.members:
1015 member.closeHandle()
1016
1017
1018
1020 # Build default Metachip based on unmodified header WCS values
1021 output_wcs = self.buildMetachip()
1022
1023 self.size = (output_wcs.naxis1,output_wcs.naxis2)
1024 self.product = Exposure(output,wcs=output_wcs,new=yes)
1025 self.product.geometry.wcslin = self.product.geometry.wcs.copy()
1026
1027 # Compute default offsets between all inputs...
1028 # This will be used when applying relative shifts from ASN table
1029 # or shiftfile.
1030 self.computeOffsets()
1031 # Preserve default DitherProduct Metachip WCS as wcslin
1032 self.product.exptime = self.exptime
1033 # Update the corners arrays for the product now...
1034 self.product.setCorners()
1035 #self.product.corners['corrected'] = self.product.corners['raw']
1036 _prodcorners = []
1037 for prod in self.members:
1038 _prodcorners += prod.product.corners['corrected'].tolist()
1039 self.product.corners['corrected'] = np.array(_prodcorners,dtype=np.float64)
1040
1042 """
1043 This method will rely on final product's 'rd2xy' method
1044 to compute offsets between the different input chips.
1045 """
1046 ref_wcs = self.product.geometry.wcs
1047 ref_crv= (ref_wcs.crval1,ref_wcs.crval2)
1048
1049 _tot_ref = {'pix_shift':(0.,0.),'ra_shift':(0.,0.),
1050 'val':999999999.,'name':''}
1051
1052 _member_num = 0
1053 for member in self.members:
1054 in_wcs = member.product.geometry.wcs
1055
1056 x,y = ref_wcs.rd2xy((in_wcs.crval1,in_wcs.crval2))
1057 xoff = x - ref_wcs.crpix1
1058 yoff = y - ref_wcs.crpix2
1059
1060 raoff = (in_wcs.crval1 - ref_wcs.crval1,in_wcs.crval2 - ref_wcs.crval2)
1061
1062 _delta_rot = in_wcs.orient - ref_wcs.orient
1063 # Determine which image has the smallest offset
1064 _tot = np.sqrt(np.power(xoff,2)+np.power(yoff,2))
1065
1066 """
1067 # Use first image as reference
1068 if _member_num == 0:
1069 _tot_ref['val'] = _tot
1070 _tot_ref['pix_shift'] = (xoff,yoff)
1071 _tot_ref['name'] = member.rootname
1072 _tot_ref['ra_shift'] = raoff
1073 # This will be used primarily for manual verification
1074 _tot_ref['wcs'] = in_wcs
1075 """
1076
1077 _member_num += 1
1078
1079 # Keep track of the results for later use/comparison
1080 member.offsets = {'pixels':(xoff,yoff),'arcseconds':raoff}
1081 """
1082 for member in self.members:
1083 member.refimage = _tot_ref
1084 """
1085
1087 """
1088 Add up the exposure time for all the members in
1089 the pattern, since 'drizzle' doesn't have the necessary
1090 information to correctly set this itself.
1091 """
1092 _exptime = 0.
1093 _start = []
1094 _end = []
1095 for member in self.members:
1096 _memexp = member.product.exptime
1097 _exptime = _exptime + _memexp[0]
1098 _start.append(_memexp[1])
1099 _end.append(_memexp[2])
1100 _expstart = min(_start)
1101 _expend = max(_end)
1102
1103 return (_exptime,_expstart,_expend)
1104
1106 """
1107 For each entry in prodlist, append the appropriate type
1108 of Observation to the members list.
1109
1110 If it's a moving target observation apply the wcs of the first
1111 observation to all other observations.
1112 """
1113
1114 member = prodlist['order'][0]
1115 filename = fileutil.buildRootname(member)
1116 mtflag = fileutil.getKeyword(filename, 'MTFLAG')
1117
1118 if mtflag == 'T':
1119 mtpars = pars.copy() # necessary in order to avoid problems with use of pars in next loop
1120 mt_member = selectInstrument(filename,output, pars=mtpars)
1121 mt_wcs = {}
1122 for member in mt_member.members:
1123 mt_wcs[member.chip] = member.geometry.wcs
1124
1125 pars['mt_wcs'] = mt_wcs
1126 del mt_member
1127
1128 for memname in prodlist['order']:
1129 pardict = self.pars[memname]
1130 pardict.update(pars)
1131 filename = fileutil.buildRootname(memname)
1132 if filename:
1133 self.members.append(selectInstrument(filename,output,pars=pardict))
1134 else:
1135 print 'No recognizable input! Not building parameters for ',memname
1136
1138 """ Return the class instance for the member with name memname."""
1139 member = None
1140 for mem in self.members:
1141 member = mem.getMember(memname)
1142 if member != None:
1143 break
1144 return member
1145
1147 """ Returns a dictionary with one key for each member.
1148 The value for each key is a list of all the extension/chip
1149 names that make up each member.
1150 Output: {member1:[extname1,extname2,...],...}
1151 """
1152 memlist = []
1153 for member in self.members:
1154 memlist.extend(member.getMemberNames())
1155
1156 return memlist
1157
1159
1160 # If an output field is specified, update product WCS to match
1161 if ref != None:
1162 _field = ref
1163 _field.mergeWCS(self.product.geometry.wcslin)
1164 #_field.exptime = self.exptime
1165
1166 # Update product WCS based on input
1167 self.transformMetachip(_field)
1168 #_field.mergeWCS(self.product.geometry.wcs)
1169 else:
1170 _field = SkyField()
1171 _field.mergeWCS(self.product.geometry.wcslin)
1172 # Reset Product WCS to reflect use of default WCS
1173 self.product.geometry.wcs = self.product.geometry.wcslin.copy()
1174 _field.exptime = self.exptime
1175 # Specify that this product comes represents a dither product
1176 # not just a single observation product.
1177 _field.dither = yes
1178
1179 # Copy the new reference WCS into the product WCS
1180 self.product.geometry.wcs = _field.wcs.copy()
1181 parlist = []
1182 for member in self.members:
1183 parlist = parlist + member.buildPars(ref=_field)
1184
1185 # Set value for number of images combined by PyDrizzle
1186 # based on DitherPattern attribute.
1187 for pl in parlist:
1188 pl['nimages'] = self.nimages
1189
1190 return parlist
1191
1193 """
1194 This method combines the results of the member's buildMetachip()
1195 methods into a composite WCS.
1196 (DitherProduct method)
1197 """
1198 prodlist = []
1199
1200 _ref = self.members[0].product.geometry
1201
1202 # Get pscale from model
1203 pscale = _ref.model.pscale
1204
1205 for member in self.members:
1206 # merge the corrected shapes into a corrected meta-chip here
1207 # Start by computing the corected positions for reference points
1208 prodlist.append(member.product)
1209
1210 # Use first input image WCS as initial reference WCS
1211 meta_wcs = _ref.wcs.copy()
1212
1213 ref_crval = (meta_wcs.crval1,meta_wcs.crval2)
1214 ref_crpix = (meta_wcs.crpix1,meta_wcs.crpix2)
1215
1216 # Specify the output orientation angle
1217 ref_orient = meta_wcs.orient
1218
1219 _delta_x = _delta_y = 0.0
1220 _xr,_yr,_dpx,_dpy = [],[],[],[]
1221 # Compute offsets for each input relative to this reference WCS
1222 for member in prodlist:
1223 _wcs = member.geometry.wcs
1224
1225 # Rigorously compute the orientation changes from WCS
1226 # information using algorithm provided by R. Hook from WDRIZZLE.
1227 abxt,cdyt = drutil.wcsfit(member.geometry, meta_wcs)
1228
1229 # Compute the rotation between input and reference from fit coeffs.
1230 _angle = RADTODEG(np.arctan2(abxt[1],cdyt[0]))
1231 _dpos = (abxt[2],cdyt[2])
1232
1233 _delta_x += _dpos[0]
1234 _delta_y += _dpos[1]
1235 _dpx.append(_dpos[0])
1236 _dpy.append(_dpos[1])
1237
1238 # Compute the range of pixels each input spans in the
1239 # output meta-frame coordinates
1240 # Each product in this list could have a different plate scale.
1241 # apply
1242 _scale = _wcs.pscale / meta_wcs.pscale
1243
1244 # Now, account for any rotation difference between
1245 # input and output frames
1246 #_angle = meta_wcs.orient - _wcs.orient
1247 _xrange,_yrange = drutil.getRotatedSize(member.corners['corrected'],_angle)
1248
1249 #_range = [(_xrange[1] - _xrange[0]),(_yrange[1] - _yrange[0])]
1250 #_range[0] *= _scale
1251 #_range[1] *= _scale
1252
1253 _xr.append(_dpos[0] + _xrange[0]*_scale)
1254 _xr.append(_dpos[0] + _xrange[1]*_scale)
1255 _yr.append(_dpos[1] + _yrange[0]*_scale)
1256 _yr.append(_dpos[1] + _yrange[1]*_scale)
1257
1258 # Determine the full size of the metachip
1259 _xmin = np.minimum.reduce(_xr)
1260 _ymin = np.minimum.reduce(_yr)
1261 _xmax = np.maximum.reduce(_xr)
1262 _ymax = np.maximum.reduce(_yr)
1263
1264 _dxmin = np.minimum.reduce(_dpx)
1265 _dymin = np.minimum.reduce(_dpy)
1266 _dxmax = np.maximum.reduce(_dpx)
1267 _dymax = np.maximum.reduce(_dpy)
1268
1269 _nimg = len(prodlist)
1270 _delta_x /= _nimg
1271 _delta_y /= _nimg
1272
1273 # Computes the offset from center of the overall set of shifts
1274 # as applied to the pixels. This insures that the visible pixels
1275 # are centered in the output.
1276 _delta_dx = ((_xmax - _delta_x) + (_xmin - _delta_x))/2.
1277 _delta_dy = ((_ymax - _delta_y) + (_ymin - _delta_y))/2.
1278
1279 if _xmin > 0.: _xmin = 0.
1280 if _ymin > 0.: _ymin = 0.
1281
1282 # Need to account for overall offset in images.
1283 # by adding in average offset induced by shifts to
1284 # output size: delta_dx/dy.
1285 # This accounts for relative offsets that are not centered
1286 # on the final output image.
1287 #xsize = int(_xmax - _xmin + 2.0 + abs(_delta_dx) )
1288 #ysize = int(_ymax - _ymin + 2.0 + abs(_delta_dy) )
1289 xsize = int(_xmax - _xmin)
1290 ysize = int(_ymax - _ymin)
1291
1292 # Compute difference between new size and
1293 # original size of meta_wcs frame.
1294 _dxsize = xsize - meta_wcs.naxis1
1295 _dysize = ysize - meta_wcs.naxis2
1296
1297 # Update reference WCS for new size
1298 meta_wcs.naxis1 = xsize
1299 meta_wcs.naxis2 = ysize
1300 _cen = (float(xsize)/2.,float(ysize)/2.)
1301
1302 # Determine offset of centers from center of image.
1303 # d[x/y]size : difference between input and output frame sizes
1304 # delta_[x/y] : average of all shifts applied to input images
1305 # delta_d[x/y]: amount off-center of all shifts
1306 #
1307 # Need to take into account overall offset in shifts, such that
1308 # minimum, not maximum, shift is always 0.0. This is needed to allow
1309 # subarrays to align properly with full frame observations without
1310 # introducing an overall offset to the output. WJH 13 Sept 2004.
1311 #
1312 _nref = ( _dxsize/2. - _delta_x - _delta_dx, _dysize/2. - _delta_y - _delta_dy)
1313
1314 # Have to adjust the CRPIX by how much the center needs to shift
1315 # to fit into the reference frame.
1316 meta_wcs.crpix1 = meta_wcs.crpix1 + _nref[0]
1317 meta_wcs.crpix2 = meta_wcs.crpix2 + _nref[1]
1318 meta_wcs.recenter()
1319
1320 return meta_wcs
1321
1322 #
1323 #
1324 # Set up default pars dictionary for external calls
1325
1327 """
1328 Method which encapsulates the logic for determining
1329 which class to instantiate for each file.
1330
1331 """
1332 # Determine the instrument...
1333 instrument = fileutil.getKeyword(filename+'[0]','INSTRUME')
1334
1335 #try:
1336 # ... then create an appropriate object.
1337 if instrument == INSTRUMENT[0]:
1338 member = ACSObservation(filename,output,pars=pars)
1339 elif instrument == INSTRUMENT[1]:
1340 member = WFPCObservation(filename,output,pars=pars)
1341 elif instrument == INSTRUMENT[2]:
1342 member = STISObservation(filename,output,pars=pars)
1343 elif instrument == INSTRUMENT[3]:
1344 member = NICMOSObservation(filename,output,pars=pars)
1345 elif instrument == INSTRUMENT[4]:
1346 member = WFC3Observation(filename,output,pars=pars)
1347 else:
1348 #raise AttributeError, "Instrument '%s' not supported now."%instrument
1349 member = GenericObservation(filename,output,pars=pars)
1350 #except:
1351 # raise IOError,"Image %s could not be processed."%filename
1352
1353 return member
1354
| Home | Trees | Indices | Help |
|
|---|
| Generated by Epydoc 3.0.1 on Mon Aug 22 14:36:03 2011 | http://epydoc.sourceforge.net |