source: trunk/python/scantable.py@ 989

Last change on this file since 989 was 989, checked in by mar637, 19 years ago

reactivated scantable.save as in Ticket #4.
added scantable.set_dirframe as in ticket #13

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 54.3 KB
Line 
1from asap._asap import Scantable
2from asap import rcParams
3from asap import print_log, asaplog
4from asap import selector
5from numarray import ones,zeros
6import sys
7
8class scantable(Scantable):
9 """
10 The ASAP container for scans
11 """
12
13 def __init__(self, filename, average=None, unit=None):
14 """
15 Create a scantable from a saved one or make a reference
16 Parameters:
17 filename: the name of an asap table on disk
18 or
19 the name of a rpfits/sdfits/ms file
20 (integrations within scans are auto averaged
21 and the whole file is read)
22 or
23 [advanced] a reference to an existing
24 scantable
25 average: average all integrations withinb a scan on read.
26 The default (True) is taken from .asaprc.
27 unit: brightness unit; must be consistent with K or Jy.
28 Over-rides the default selected by the reader
29 (input rpfits/sdfits/ms) or replaces the value
30 in existing scantables
31 """
32 if average is None:
33 average = rcParams['scantable.autoaverage']
34 varlist = vars()
35 from asap._asap import stmath
36 self._math = stmath()
37 if isinstance(filename, Scantable):
38 Scantable.__init__(self, filename)
39 else:
40 if isinstance(filename,str):
41 import os.path
42 filename = os.path.expandvars(filename)
43 filename = os.path.expanduser(filename)
44 if not os.path.exists(filename):
45 s = "File '%s' not found." % (filename)
46 if rcParams['verbose']:
47 asaplog.push(s)
48 print asaplog.pop().strip()
49 return
50 raise IOError(s)
51 if os.path.isdir(filename):
52 # crude check if asap table
53 if os.path.exists(filename+'/table.info'):
54 Scantable.__init__(self, filename, "memory")
55 if unit is not None:
56 self.set_fluxunit(unit)
57 self.set_freqframe(rcParams['scantable.freqframe'])
58 else:
59 msg = "The given file '%s'is not a valid asap table." % (filename)
60 if rcParams['verbose']:
61 print msg
62 return
63 else:
64 raise IOError(msg)
65 else:
66 self._fill([filename],unit, average)
67 elif (isinstance(filename,list) or isinstance(filename,tuple)) \
68 and isinstance(filename[-1], str):
69 self._fill(filename, unit, average)
70 print_log()
71
72 def save(self, name=None, format=None, overwrite=False):
73 """
74 Store the scantable on disk. This can be an asap (aips++) Table, SDFITS,
75 Image FITS or MS2 format.
76 Parameters:
77 name: the name of the outputfile. For format="FITS" this
78 is the directory file name into which all the files
79 will be written (default is 'asap_FITS'). For format
80 "ASCII" this is the root file name (data in 'name'.txt
81 and header in 'name'_header.txt)
82 format: an optional file format. Default is ASAP.
83 Allowed are - 'ASAP' (save as ASAP [aips++] Table),
84 'SDFITS' (save as SDFITS file)
85 'FITS' (saves each row as a FITS Image)
86 'ASCII' (saves as ascii text file)
87 'MS2' (saves as an aips++
88 MeasurementSet V2)
89 overwrite: If the file should be overwritten if it exists.
90 The default False is to return with warning
91 without writing the output. USE WITH CARE.
92 Example:
93 scan.save('myscan.asap')
94 scan.save('myscan.sdfits','SDFITS')
95 """
96 from os import path
97 if format is None: format = rcParams['scantable.save']
98 suffix = '.'+format.lower()
99 if name is None or name =="":
100 name = 'scantable'+suffix
101 from asap import asaplog
102 msg = "No filename given. Using default name %s..." % name
103 asaplog.push(msg)
104 name = path.expandvars(name)
105 if path.isfile(name) or path.isdir(name):
106 if not overwrite:
107 msg = "File %s exists." % name
108 if rcParams['verbose']:
109 print msg
110 return
111 else:
112 raise IOError(msg)
113 format2 = format.upper()
114 if format2 == 'ASAP':
115 self._save(name)
116 else:
117 from asap._asap import stwriter as stw
118 w = stw(format2)
119 w.write(self, name)
120 print_log()
121 return
122
123 def copy(self):
124 """
125 Return a copy of this scantable.
126 Parameters:
127 none
128 Example:
129 copiedscan = scan.copy()
130 """
131 sd = scantable(Scantable._copy(self))
132 return sd
133
134 def get_scan(self, scanid=None):
135 """
136 Return a specific scan (by scanno) or collection of scans (by
137 source name) in a new scantable.
138 Parameters:
139 scanid: a (list of) scanno or a source name, unix-style
140 patterns are accepted for source name matching, e.g.
141 '*_R' gets all 'ref scans
142 Example:
143 # get all scans containing the source '323p459'
144 newscan = scan.get_scan('323p459')
145 # get all 'off' scans
146 refscans = scan.get_scan('*_R')
147 # get a susbset of scans by scanno (as listed in scan.summary())
148 newscan = scan.get_scan([0,2,7,10])
149 """
150 if scanid is None:
151 if rcParams['verbose']:
152 print "Please specify a scan no or name to retrieve from the scantable"
153 return
154 else:
155 raise RuntimeError("No scan given")
156
157 try:
158 bsel = self.get_selection()
159 sel = selector()
160 if type(scanid) is str:
161 sel.set_name(scanid)
162 self.set_selection(bsel+sel)
163 scopy = self._copy()
164 self.set_selection(bsel)
165 return scantable(scopy)
166 elif type(scanid) is int:
167 sel.set_scans([scanid])
168 self.set_selection(bsel+sel)
169 scopy = self._copy()
170 self.set_selection(bsel)
171 return scantable(scopy)
172 elif type(scanid) is list:
173 sel.set_scans(scanid)
174 self.set_selection(sel)
175 scopy = self._copy()
176 self.set_selection(bsel)
177 return scantable(scopy)
178 else:
179 msg = "Illegal scanid type, use 'int' or 'list' if ints."
180 if rcParams['verbose']:
181 print msg
182 else:
183 raise TypeError(msg)
184 except RuntimeError:
185 if rcParams['verbose']: print "Couldn't find any match."
186 else: raise
187
188 def __str__(self):
189 return Scantable._summary(self,True)
190
191 def summary(self, filename=None):
192 """
193 Print a summary of the contents of this scantable.
194 Parameters:
195 filename: the name of a file to write the putput to
196 Default - no file output
197 verbose: print extra info such as the frequency table
198 The default (False) is taken from .asaprc
199 """
200 info = Scantable._summary(self, True)
201 #if verbose is None: verbose = rcParams['scantable.verbosesummary']
202 if filename is not None:
203 if filename is "":
204 filename = 'scantable_summary.txt'
205 from os.path import expandvars, isdir
206 filename = expandvars(filename)
207 if not isdir(filename):
208 data = open(filename, 'w')
209 data.write(info)
210 data.close()
211 else:
212 msg = "Illegal file name '%s'." % (filename)
213 if rcParams['verbose']:
214 print msg
215 else:
216 raise IOError(msg)
217 if rcParams['verbose']:
218 try:
219 from IPython.genutils import page as pager
220 except ImportError:
221 from pydoc import pager
222 pager(info)
223 else:
224 return info
225
226
227 def get_selection(self):
228 """
229 """
230 return selector(self._getselection())
231
232 def set_selection(self, selection):
233 """
234 """
235 self._setselection(selection)
236
237 def set_cursor(self, beam=0, IF=0, pol=0):
238 """
239 Set the spectrum for individual operations.
240 Parameters:
241 beam, IF, pol: a number
242 Example:
243 scan.set_cursor(0,0,1)
244 pol1sig = scan.stats(all=False) # returns std dev for beam=0
245 # if=0, pol=1
246 """
247 print "DEPRECATED"
248 varlist = vars()
249 sel = asap._asap.Selector()
250 sel._setbeams([beam])
251 sel._setpols([pol])
252 sel._setifs([IF])
253 self._add_history("set_cursor", varlist)
254 return
255
256 def get_cursor(self):
257 """
258 Return/print a the current 'cursor' into the Beam/IF/Pol cube.
259 Parameters:
260 none
261 Returns:
262 a list of values (currentBeam,currentIF,currentPol)
263 Example:
264 none
265 """
266 print "DEPRECATED"
267 sel = self._getselection()
268 i = sel.getbeams()[0]
269 j = sel.getifs()[0]
270 k = sel.getpols()[0]
271 from asap import asaplog
272 out = "--------------------------------------------------\n"
273 out += " Cursor position\n"
274 out += "--------------------------------------------------\n"
275 out += 'Beam=%d IF=%d Pol=%d ' % (i,j,k)
276 asaplog.push(out)
277 print_log()
278 return i,j,k
279
280 def stats(self, stat='stddev', mask=None):
281 """
282 Determine the specified statistic of the current beam/if/pol
283 Takes a 'mask' as an optional parameter to specify which
284 channels should be excluded.
285 Parameters:
286 stat: 'min', 'max', 'sumsq', 'sum', 'mean'
287 'var', 'stddev', 'avdev', 'rms', 'median'
288 mask: an optional mask specifying where the statistic
289 should be determined.
290 Example:
291 scan.set_unit('channel')
292 msk = scan.create_mask([100,200],[500,600])
293 scan.stats(stat='mean', mask=m)
294 """
295 from numarray import array,zeros,Float
296 if mask == None:
297 mask = []
298 axes = ['Beam','IF','Pol','Time']
299 if not self._check_ifs():
300 raise ValueError("Cannot apply mask as the IFs have different number of channels"
301 "Please use setselection() to select individual IFs")
302
303 statvals = self._math._stats(self, mask, stat)
304 out = ''
305 axes = []
306 for i in range(self.nrow()):
307 axis = []
308 axis.append(self.getscan(i))
309 axis.append(self.getbeam(i))
310 axis.append(self.getif(i))
311 axis.append(self.getpol(i))
312 axis.append(self.getcycle(i))
313 axes.append(axis)
314 tm = self._gettime(i)
315 src = self._getsourcename(i)
316 out += 'Scan[%d] (%s) ' % (axis[0], src)
317 out += 'Time[%s]:\n' % (tm)
318 if self.nbeam(-1) > 1: out += ' Beam[%d] ' % (axis[1])
319 if self.nif(-1) > 1: out += ' IF[%d] ' % (axis[2])
320 if self.npol(-1) > 1: out += ' Pol[%d] ' % (axis[3])
321 out += '= %3.3f\n' % (statvals[i])
322 out += "--------------------------------------------------\n"
323
324 if rcParams['verbose']:
325 print "--------------------------------------------------"
326 print " ",stat
327 print "--------------------------------------------------"
328 print out
329 retval = { 'axesnames': ['scanno','beamno','ifno','polno','cycleno'],
330 'axes' : axes,
331 'data': statvals}
332 return retval
333
334 def stddev(self,mask=None):
335 """
336 Determine the standard deviation of the current beam/if/pol
337 Takes a 'mask' as an optional parameter to specify which
338 channels should be excluded.
339 Parameters:
340 mask: an optional mask specifying where the standard
341 deviation should be determined.
342
343 Example:
344 scan.set_unit('channel')
345 msk = scan.create_mask([100,200],[500,600])
346 scan.stddev(mask=m)
347 """
348 return self.stats(stat='stddev',mask=mask);
349
350 def get_tsys(self):
351 """
352 Return the System temperatures.
353 Parameters:
354
355 Returns:
356 a list of Tsys values for the current selection
357 """
358
359 return self._row_callback(self._gettsys, "Tsys")
360
361 def _row_callback(self, callback, label):
362 axes = []
363 axesnames = ['scanno','beamno','ifno','polno','cycleno']
364 out = ""
365 outvec =[]
366 for i in range(self.nrow()):
367 axis = []
368 axis.append(self.getscan(i))
369 axis.append(self.getbeam(i))
370 axis.append(self.getif(i))
371 axis.append(self.getpol(i))
372 axis.append(self.getcycle(i))
373 axes.append(axis)
374 tm = self._gettime(i)
375 src = self._getsourcename(i)
376 out += 'Scan[%d] (%s) ' % (axis[0], src)
377 out += 'Time[%s]:\n' % (tm)
378 if self.nbeam(-1) > 1: out += ' Beam[%d] ' % (axis[1])
379 if self.nif(-1) > 1: out += ' IF[%d] ' % (axis[2])
380 if self.npol(-1) > 1: out += ' Pol[%d] ' % (axis[3])
381 outvec.append(callback(i))
382 out += '= %3.3f\n' % (outvec[i])
383 out += "--------------------------------------------------\n"
384 if rcParams['verbose']:
385 print "--------------------------------------------------"
386 print " %s" % (label)
387 print "--------------------------------------------------"
388 print out
389 retval = {'axesnames': axesnames, 'axes': axes, 'data': outvec}
390 return retval
391
392
393 def get_time(self, row=-1):
394 """
395 Get a list of time stamps for the observations.
396 Return a string for each integration in the scantable.
397 Parameters:
398 row: row no of integration. Default -1 return all rows
399 Example:
400 none
401 """
402 out = []
403 if row == -1:
404 for i in range(self.nrow()):
405 out.append(self._gettime(i))
406 return out
407 else:
408 if row < self.nrow():
409 return self._gettime(row)
410
411 def get_sourcename(self, row=-1):
412 """
413 Get a list source names for the observations.
414 Return a string for each integration in the scantable.
415 Parameters:
416 row: row no of integration. Default -1 return all rows
417 Example:
418 none
419 """
420 out = []
421 if row == -1:
422 return [self._getsourcename(i) for i in range(self.nrow())]
423 else:
424 if 0 <= row < self.nrow():
425 return self._getsourcename(row)
426
427 def get_elevation(self, row=-1):
428 """
429 Get a list of elevations for the observations.
430 Return a float for each integration in the scantable.
431 Parameters:
432 row: row no of integration. Default -1 return all rows
433 Example:
434 none
435 """
436 out = []
437 if row == -1:
438 return [self._getelevation(i) for i in range(self.nrow())]
439 else:
440 if 0 <= row < self.nrow():
441 return self._getelevation(row)
442
443 def get_azimuth(self, row=-1):
444 """
445 Get a list of azimuths for the observations.
446 Return a float for each integration in the scantable.
447 Parameters:
448 row: row no of integration. Default -1 return all rows
449 Example:
450 none
451 """
452 out = []
453 if row == -1:
454 return [self._getazimuth(i) for i in range(self.nrow())]
455 else:
456 if 0 <= row < self.nrow():
457 return self._getazimuth(row)
458
459 def get_parangle(self, row=-1):
460 """
461 Get a list of parallactic angles for the observations.
462 Return a float for each integration in the scantable.
463 Parameters:
464 row: row no of integration. Default -1 return all rows
465 Example:
466 none
467 """
468 out = []
469 if row == -1:
470 return [self._getparangle(i) for i in range(self.nrow())]
471 else:
472 if 0 <= row < self.nrow():
473 return self._getparangle(row)
474
475 def set_unit(self, unit='channel'):
476 """
477 Set the unit for all following operations on this scantable
478 Parameters:
479 unit: optional unit, default is 'channel'
480 one of '*Hz','km/s','channel', ''
481 """
482 varlist = vars()
483 if unit in ['','pixel', 'channel']:
484 unit = ''
485 inf = list(self._getcoordinfo())
486 inf[0] = unit
487 self._setcoordinfo(inf)
488 self._add_history("set_unit",varlist)
489
490 def set_instrument(self, instr):
491 """
492 Set the instrument for subsequent processing
493 Parameters:
494 instr: Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
495 'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
496 """
497 self._setInstrument(instr)
498 self._add_history("set_instument",vars())
499 print_log()
500
501 def set_doppler(self, doppler='RADIO'):
502 """
503 Set the doppler for all following operations on this scantable.
504 Parameters:
505 doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
506 """
507 varlist = vars()
508 inf = list(self._getcoordinfo())
509 inf[2] = doppler
510 self._setcoordinfo(inf)
511 self._add_history("set_doppler",vars())
512 print_log()
513
514 def set_freqframe(self, frame=None):
515 """
516 Set the frame type of the Spectral Axis.
517 Parameters:
518 frame: an optional frame type, default 'LSRK'. Valid frames are:
519 'REST','TOPO','LSRD','LSRK','BARY',
520 'GEO','GALACTO','LGROUP','CMB'
521 Examples:
522 scan.set_freqframe('BARY')
523 """
524 if frame is None: frame = rcParams['scantable.freqframe']
525 varlist = vars()
526 valid = ['REST','TOPO','LSRD','LSRK','BARY', \
527 'GEO','GALACTO','LGROUP','CMB']
528
529 if frame in valid:
530 inf = list(self._getcoordinfo())
531 inf[1] = frame
532 self._setcoordinfo(inf)
533 self._add_history("set_freqframe",varlist)
534 else:
535 msg = "Please specify a valid freq type. Valid types are:\n",valid
536 if rcParams['verbose']:
537 print msg
538 else:
539 raise TypeError(msg)
540 print_log()
541
542 def set_dirframe(self, frame=""):
543 """
544 Set the frame type of the Direction on the sky.
545 Parameters:
546 frame: an optional frame type, default ''. Valid frames are:
547 'J2000', 'B1950', 'GALACTIC'
548 Examples:
549 scan.set_dirframe('GALACTIC')
550 """
551 varlist = vars()
552 try:
553 Scantable.set_dirframe(self, frame)
554 except RuntimeError,msg:
555 if rcParams['verbose']:
556 print msg
557 else:
558 raise
559 self._add_history("set_dirframe",varlist)
560
561 def get_unit(self):
562 """
563 Get the default unit set in this scantable
564 Parameters:
565 Returns:
566 A unit string
567 """
568 inf = self._getcoordinfo()
569 unit = inf[0]
570 if unit == '': unit = 'channel'
571 return unit
572
573 def get_abcissa(self, rowno=0):
574 """
575 Get the abcissa in the current coordinate setup for the currently
576 selected Beam/IF/Pol
577 Parameters:
578 rowno: an optional row number in the scantable. Default is the
579 first row, i.e. rowno=0
580 Returns:
581 The abcissa values and it's format string (as a dictionary)
582 """
583 abc = self._getabcissa(rowno)
584 lbl = self._getabcissalabel(rowno)
585 print_log()
586 return abc, lbl
587
588 def create_mask(self, *args, **kwargs):
589 """
590 Compute and return a mask based on [min,max] windows.
591 The specified windows are to be INCLUDED, when the mask is
592 applied.
593 Parameters:
594 [min,max],[min2,max2],...
595 Pairs of start/end points specifying the regions
596 to be masked
597 invert: optional argument. If specified as True,
598 return an inverted mask, i.e. the regions
599 specified are EXCLUDED
600 row: create the mask using the specified row for
601 unit conversions, default is row=0
602 only necessary if frequency varies over rows.
603 Example:
604 scan.set_unit('channel')
605
606 a)
607 msk = scan.set_mask([400,500],[800,900])
608 # masks everything outside 400 and 500
609 # and 800 and 900 in the unit 'channel'
610
611 b)
612 msk = scan.set_mask([400,500],[800,900], invert=True)
613 # masks the regions between 400 and 500
614 # and 800 and 900 in the unit 'channel'
615
616 """
617 row = 0
618 if kwargs.has_key("row"):
619 row = kwargs.get("row")
620 data = self._getabcissa(row)
621 u = self._getcoordinfo()[0]
622 if rcParams['verbose']:
623 if u == "": u = "channel"
624 from asap import asaplog
625 msg = "The current mask window unit is %s" % u
626 if not self._check_ifs():
627 msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
628 asaplog.push(msg)
629 n = self.nchan()
630 msk = zeros(n)
631 # test if args is a 'list' or a 'normal *args - UGLY!!!
632
633 ws = (isinstance(args[-1][-1],int) or isinstance(args[-1][-1],float)) and args or args[0]
634 for window in ws:
635 if (len(window) != 2 or window[0] > window[1] ):
636 raise TypeError("A window needs to be defined as [min,max]")
637 for i in range(n):
638 if data[i] >= window[0] and data[i] < window[1]:
639 msk[i] = 1
640 if kwargs.has_key('invert'):
641 if kwargs.get('invert'):
642 from numarray import logical_not
643 msk = logical_not(msk)
644 print_log()
645 return msk
646
647 def get_restfreqs(self):
648 """
649 Get the restfrequency(s) stored in this scantable.
650 The return value(s) are always of unit 'Hz'
651 Parameters:
652 none
653 Returns:
654 a list of doubles
655 """
656 return list(self._getrestfreqs())
657
658
659 def set_restfreqs(self, freqs=None, unit='Hz'):
660 """
661 Set or replace the restfrequency specified and
662 If the 'freqs' argument holds a scalar,
663 then that rest frequency will be applied to all the selected
664 data. If the 'freqs' argument holds
665 a vector, then it MUST be of equal or smaller length than
666 the number of IFs (and the available restfrequencies will be
667 replaced by this vector). In this case, *all* data have
668 the restfrequency set per IF according
669 to the corresponding value you give in the 'freqs' vector.
670 E.g. 'freqs=[1e9,2e9]' would mean IF 0 gets restfreq 1e9 and
671 IF 1 gets restfreq 2e9.
672 You can also specify the frequencies via known line names
673 from the built-in Lovas table.
674 Parameters:
675 freqs: list of rest frequency values or string idenitfiers
676 unit: unit for rest frequency (default 'Hz')
677
678 Example:
679 # set the given restfrequency for the whole table
680 scan.set_restfreqs(freqs=1.4e9)
681 # If thee number of IFs in the data is >= 2 the IF0 gets the first
682 # value IF1 the second...
683 scan.set_restfreqs(freqs=[1.4e9,1.67e9])
684 #set the given restfrequency for the whole table (by name)
685 scan.set_restfreqs(freqs="OH1667")
686
687 Note:
688 To do more sophisticate Restfrequency setting, e.g. on a
689 source and IF basis, use scantable.set_selection() before using
690 this function.
691 # provide your scantable is call scan
692 selection = selector()
693 selection.set_name("ORION*")
694 selection.set_ifs([1])
695 scan.set_selection(selection)
696 scan.set_restfreqs(freqs=86.6e9)
697
698 """
699 varlist = vars()
700
701 t = type(freqs)
702 if isinstance(freqs, int) or isinstance(freqs,float):
703 self._setrestfreqs(freqs, unit)
704 elif isinstance(freqs, list) or isinstance(freqs,tuple):
705 if isinstance(freqs[-1], int) or isinstance(freqs[-1],float):
706 sel = selector()
707 savesel = self._getselection()
708 for i in xrange(len(freqs)):
709 sel.set_ifs([i])
710 self._setselection(sel)
711 self._setrestfreqs(freqs[i], unit)
712 self._setselection(savesel)
713 elif isinstance(freqs[-1], str):
714 # not yet implemented
715 pass
716 else:
717 return
718 self._add_history("set_restfreqs", varlist)
719
720
721
722 def history(self):
723 hist = list(self._gethistory())
724 out = "-"*80
725 for h in hist:
726 if h.startswith("---"):
727 out += "\n"+h
728 else:
729 items = h.split("##")
730 date = items[0]
731 func = items[1]
732 items = items[2:]
733 out += "\n"+date+"\n"
734 out += "Function: %s\n Parameters:" % (func)
735 for i in items:
736 s = i.split("=")
737 out += "\n %s = %s" % (s[0],s[1])
738 out += "\n"+"-"*80
739 try:
740 from IPython.genutils import page as pager
741 except ImportError:
742 from pydoc import pager
743 pager(out)
744 return
745
746 #
747 # Maths business
748 #
749
750 def average_time(self, mask=None, scanav=False, weight='tint', align=False):
751 """
752 Return the (time) average of a scan, or apply it 'insitu'.
753 Note:
754 in channels only
755 The cursor of the output scan is set to 0.
756 Parameters:
757 one scan or comma separated scans
758 mask: an optional mask (only used for 'var' and 'tsys'
759 weighting)
760 scanav: True averages each scan separately
761 False (default) averages all scans together,
762 weight: Weighting scheme. 'none', 'var' (1/var(spec)
763 weighted), 'tsys' (1/Tsys**2 weighted), 'tint'
764 (integration time weighted) or 'tintsys' (Tint/Tsys**2).
765 The default is 'tint'
766 align: align the spectra in velocity before averaging. It takes
767 the time of the first spectrum as reference time.
768 Example:
769 # time average the scantable without using a mask
770 newscan = scan.average_time()
771 """
772 varlist = vars()
773 if weight is None: weight = 'TINT'
774 if mask is None: mask = ()
775 if scanav:
776 scanav = "SCAN"
777 else:
778 scanav = "NONE"
779 scan = (self,)
780 try:
781 if align:
782 scan = (self.freq_align(insitu=False),)
783 s = scantable(self._math._average(scan, mask, weight.upper(),
784 scanav))
785 except RuntimeError,msg:
786 if rcParams['verbose']:
787 print msg
788 return
789 else: raise
790 s._add_history("average_time",varlist)
791 print_log()
792 return s
793
794 def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None):
795 """
796 Return a scan where all spectra are converted to either
797 Jansky or Kelvin depending upon the flux units of the scan table.
798 By default the function tries to look the values up internally.
799 If it can't find them (or if you want to over-ride), you must
800 specify EITHER jyperk OR eta (and D which it will try to look up
801 also if you don't set it). jyperk takes precedence if you set both.
802 Parameters:
803 jyperk: the Jy / K conversion factor
804 eta: the aperture efficiency
805 d: the geomtric diameter (metres)
806 insitu: if False a new scantable is returned.
807 Otherwise, the scaling is done in-situ
808 The default is taken from .asaprc (False)
809 allaxes: if True apply to all spectra. Otherwise
810 apply only to the selected (beam/pol/if)spectra only
811 The default is taken from .asaprc (True if none)
812 """
813 if insitu is None: insitu = rcParams['insitu']
814 self._math._setinsitu(insitu)
815 varlist = vars()
816 if jyperk is None: jyperk = -1.0
817 if d is None: d = -1.0
818 if eta is None: eta = -1.0
819 s = scantable(self._math._convertflux(self, d, eta, jyperk))
820 s._add_history("convert_flux", varlist)
821 print_log()
822 if insitu: self._assign(s)
823 else: return s
824
825 def gain_el(self, poly=None, filename="", method="linear", insitu=None):
826 """
827 Return a scan after applying a gain-elevation correction.
828 The correction can be made via either a polynomial or a
829 table-based interpolation (and extrapolation if necessary).
830 You specify polynomial coefficients, an ascii table or neither.
831 If you specify neither, then a polynomial correction will be made
832 with built in coefficients known for certain telescopes (an error
833 will occur if the instrument is not known).
834 The data and Tsys are *divided* by the scaling factors.
835 Parameters:
836 poly: Polynomial coefficients (default None) to compute a
837 gain-elevation correction as a function of
838 elevation (in degrees).
839 filename: The name of an ascii file holding correction factors.
840 The first row of the ascii file must give the column
841 names and these MUST include columns
842 "ELEVATION" (degrees) and "FACTOR" (multiply data
843 by this) somewhere.
844 The second row must give the data type of the
845 column. Use 'R' for Real and 'I' for Integer.
846 An example file would be
847 (actual factors are arbitrary) :
848
849 TIME ELEVATION FACTOR
850 R R R
851 0.1 0 0.8
852 0.2 20 0.85
853 0.3 40 0.9
854 0.4 60 0.85
855 0.5 80 0.8
856 0.6 90 0.75
857 method: Interpolation method when correcting from a table.
858 Values are "nearest", "linear" (default), "cubic"
859 and "spline"
860 insitu: if False a new scantable is returned.
861 Otherwise, the scaling is done in-situ
862 The default is taken from .asaprc (False)
863 """
864
865 if insitu is None: insitu = rcParams['insitu']
866 self._math._setinsitu(insitu)
867 varlist = vars()
868 if poly is None:
869 poly = ()
870 from os.path import expandvars
871 filename = expandvars(filename)
872 s = scantable(self._math._gainel(self, poly, filename, method))
873 s._add_history("gain_el", varlist)
874 print_log()
875 if insitu: self._assign(s)
876 else: return s
877
878 def freq_align(self, reftime=None, method='cubic', insitu=None):
879 """
880 Return a scan where all rows have been aligned in frequency/velocity.
881 The alignment frequency frame (e.g. LSRK) is that set by function
882 set_freqframe.
883 Parameters:
884 reftime: reference time to align at. By default, the time of
885 the first row of data is used.
886 method: Interpolation method for regridding the spectra.
887 Choose from "nearest", "linear", "cubic" (default)
888 and "spline"
889 insitu: if False a new scantable is returned.
890 Otherwise, the scaling is done in-situ
891 The default is taken from .asaprc (False)
892 """
893 if insitu is None: insitu = rcParams["insitu"]
894 self._math._setinsitu(insitu)
895 varlist = vars()
896 if reftime is None: reftime = ""
897 s = scantable(self._math._freq_align(self, reftime, method))
898 s._add_history("freq_align", varlist)
899 print_log()
900 if insitu: self._assign(s)
901 else: return s
902
903 def opacity(self, tau, insitu=None):
904 """
905 Apply an opacity correction. The data
906 and Tsys are multiplied by the correction factor.
907 Parameters:
908 tau: Opacity from which the correction factor is
909 exp(tau*ZD)
910 where ZD is the zenith-distance
911 insitu: if False a new scantable is returned.
912 Otherwise, the scaling is done in-situ
913 The default is taken from .asaprc (False)
914 """
915 if insitu is None: insitu = rcParams['insitu']
916 self._math._setinsitu(insitu)
917 varlist = vars()
918 s = scantable(self._math._opacity(self, tau))
919 s._add_history("opacity", varlist)
920 print_log()
921 if insitu: self._assign(s)
922 else: return s
923
924 def bin(self, width=5, insitu=None):
925 """
926 Return a scan where all spectra have been binned up.
927 width: The bin width (default=5) in pixels
928 insitu: if False a new scantable is returned.
929 Otherwise, the scaling is done in-situ
930 The default is taken from .asaprc (False)
931 """
932 if insitu is None: insitu = rcParams['insitu']
933 self._math._setinsitu(insitu)
934 varlist = vars()
935 s = scantable(self._math._bin(self, width))
936 s._add_history("bin",varlist)
937 print_log()
938 if insitu: self._assign(s)
939 else: return s
940
941
942 def resample(self, width=5, method='cubic', insitu=None):
943 """
944 Return a scan where all spectra have been binned up
945 width: The bin width (default=5) in pixels
946 method: Interpolation method when correcting from a table.
947 Values are "nearest", "linear", "cubic" (default)
948 and "spline"
949 insitu: if False a new scantable is returned.
950 Otherwise, the scaling is done in-situ
951 The default is taken from .asaprc (False)
952 """
953 if insitu is None: insitu = rcParams['insitu']
954 self._math._setinsitu(insitu)
955 varlist = vars()
956 s = scantable(self._math._resample(self, method, width))
957 s._add_history("resample",varlist)
958 print_log()
959 if insitu: self._assign(s)
960 else: return s
961
962
963 def average_pol(self, mask=None, weight='none'):
964 """
965 Average the Polarisations together.
966 Parameters:
967 mask: An optional mask defining the region, where the
968 averaging will be applied. The output will have all
969 specified points masked.
970 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
971 weighted), or 'tsys' (1/Tsys**2 weighted)
972 """
973 varlist = vars()
974 if mask is None:
975 mask = ()
976 s = self._math._averagepol(self, mask, weight)
977 s._add_history("average_pol",varlist)
978 print_log()
979 return scantable(s)
980
981 def smooth(self, kernel="hanning", width=5.0, insitu=None):
982 """
983 Smooth the spectrum by the specified kernel (conserving flux).
984 Parameters:
985 scan: The input scan
986 kernel: The type of smoothing kernel. Select from
987 'hanning' (default), 'gaussian' and 'boxcar'.
988 The first three characters are sufficient.
989 width: The width of the kernel in pixels. For hanning this is
990 ignored otherwise it defauls to 5 pixels.
991 For 'gaussian' it is the Full Width Half
992 Maximum. For 'boxcar' it is the full width.
993 insitu: if False a new scantable is returned.
994 Otherwise, the scaling is done in-situ
995 The default is taken from .asaprc (False)
996 Example:
997 none
998 """
999 if insitu is None: insitu = rcParams['insitu']
1000 self._math._setinsitu(insitu)
1001 varlist = vars()
1002 s = scantable(self._math._smooth(self,kernel,width))
1003 s._add_history("smooth", varlist)
1004 print_log()
1005 if insitu: self._assign(s)
1006 else: return s
1007
1008
1009 def poly_baseline(self, mask=None, order=0, insitu=None):
1010 """
1011 Return a scan which has been baselined (all rows) by a polynomial.
1012 Parameters:
1013 scan: a scantable
1014 mask: an optional mask
1015 order: the order of the polynomial (default is 0)
1016 insitu: if False a new scantable is returned.
1017 Otherwise, the scaling is done in-situ
1018 The default is taken from .asaprc (False)
1019 allaxes: If True (default) apply to all spectra. Otherwise
1020 apply only to the selected (beam/pol/if)spectra only
1021 The default is taken from .asaprc (True if none)
1022 Example:
1023 # return a scan baselined by a third order polynomial,
1024 # not using a mask
1025 bscan = scan.poly_baseline(order=3)
1026 """
1027 if insitu is None: insitu = rcParams['insitu']
1028 varlist = vars()
1029 if mask is None:
1030 from numarray import ones
1031 mask = list(ones(self.nchan(-1)))
1032 from asap.asapfitter import fitter
1033 f = fitter()
1034 f.set_scan(self, mask)
1035 f.set_function(poly=order)
1036 s = f.auto_fit(insitu)
1037 s._add_history("poly_baseline", varlist)
1038 print_log()
1039 if insitu: self._assign(s)
1040 else: return s
1041
1042 def auto_poly_baseline(self, mask=[], edge=(0,0), order=0,
1043 threshold=3, insitu=None):
1044 """
1045 Return a scan which has been baselined (all rows) by a polynomial.
1046 Spectral lines are detected first using linefinder and masked out
1047 to avoid them affecting the baseline solution.
1048
1049 Parameters:
1050 mask: an optional mask retreived from scantable
1051 edge: an optional number of channel to drop at
1052 the edge of spectrum. If only one value is
1053 specified, the same number will be dropped from
1054 both sides of the spectrum. Default is to keep
1055 all channels. Nested tuples represent individual
1056 edge selection for different IFs (a number of spectral
1057 channels can be different)
1058 order: the order of the polynomial (default is 0)
1059 threshold: the threshold used by line finder. It is better to
1060 keep it large as only strong lines affect the
1061 baseline solution.
1062 insitu: if False a new scantable is returned.
1063 Otherwise, the scaling is done in-situ
1064 The default is taken from .asaprc (False)
1065
1066 Example:
1067 scan2=scan.auto_poly_baseline(order=7)
1068 """
1069 if insitu is None: insitu = rcParams['insitu']
1070 varlist = vars()
1071 from asap.asapfitter import fitter
1072 from asap.asaplinefind import linefinder
1073 from asap import _is_sequence_or_number as _is_valid
1074
1075 # check whether edge is set up for each IF individually
1076 individualEdge = False;
1077 if len(edge)>1:
1078 if isinstance(edge[0],list) or isinstance(edge[0],tuple):
1079 individualEdge = True;
1080
1081 if not _is_valid(edge, int) and not individualEdge:
1082 raise ValueError, "Parameter 'edge' has to be an integer or a \
1083 pair of integers specified as a tuple. Nested tuples are allowed \
1084 to make individual selection for different IFs."
1085
1086 curedge = (0,0)
1087 if individualEdge:
1088 for edge_par in edge:
1089 if not _is_valid(edge,int):
1090 raise ValueError, "Each element of the 'edge' tuple has \
1091 to be a pair of integers or an integer."
1092 else:
1093 curedge = edge;
1094
1095 # setup fitter
1096 f = fitter()
1097 f.set_function(poly=order)
1098
1099 # setup line finder
1100 fl=linefinder()
1101 fl.set_options(threshold=threshold)
1102
1103 if not insitu:
1104 workscan=self.copy()
1105 else:
1106 workscan=self
1107
1108 fl.set_scan(workscan)
1109
1110 rows=range(workscan.nrow())
1111 from asap import asaplog
1112 asaplog.push("Processing:")
1113 for r in rows:
1114 msg = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % (workscan.getscan(r),workscan.getbeam(r),workscan.getif(r),workscan.getpol(r), workscan.getcycle(r))
1115 asaplog.push(msg, False)
1116
1117 # figure out edge parameter
1118 if individualEdge:
1119 if len(edge)>=workscan.getif(r):
1120 raise RuntimeError, "Number of edge elements appear to be less than the number of IFs"
1121 curedge = edge[workscan.getif(r)]
1122
1123 # setup line finder
1124 fl.find_lines(r,mask,curedge)
1125 f.set_scan(workscan, fl.get_mask())
1126 f.x = workscan._getabcissa(r)
1127 f.y = workscan._getspectrum(r)
1128 f.data = None
1129 f.fit()
1130 x = f.get_parameters()
1131 workscan._setspectrum(f.fitter.getresidual(), r)
1132 workscan._add_history("poly_baseline", varlist)
1133 if insitu:
1134 self._assign(workscan)
1135 else:
1136 return workscan
1137
1138 def rotate_linpolphase(self, angle):
1139 """
1140 Rotate the phase of the complex polarization O=Q+iU correlation.
1141 This is always done in situ in the raw data. So if you call this
1142 function more than once then each call rotates the phase further.
1143 Parameters:
1144 angle: The angle (degrees) to rotate (add) by.
1145 Examples:
1146 scan.rotate_linpolphase(2.3)
1147 """
1148 varlist = vars()
1149 self._math._rotate_linpolphase(self, angle)
1150 self._add_history("rotate_linpolphase", varlist)
1151 print_log()
1152 return
1153
1154
1155 def rotate_xyphase(self, angle):
1156 """
1157 Rotate the phase of the XY correlation. This is always done in situ
1158 in the data. So if you call this function more than once
1159 then each call rotates the phase further.
1160 Parameters:
1161 angle: The angle (degrees) to rotate (add) by.
1162 Examples:
1163 scan.rotate_xyphase(2.3)
1164 """
1165 varlist = vars()
1166 self._math._rotate_xyphase(self, angle)
1167 self._add_history("rotate_xyphase", varlist)
1168 print_log()
1169 return
1170
1171 def swap_linears(self):
1172 """
1173 Swap the linear polarisations XX and YY
1174 """
1175 varlist = vars()
1176 self._math._swap_linears(self)
1177 self._add_history("swap_linears", varlist)
1178 print_log()
1179 return
1180
1181 def invert_phase(self):
1182 """
1183 Invert the phase of the complex polarisation
1184 """
1185 varlist = vars()
1186 self._math._invert_phase(self)
1187 self._add_history("invert_phase", varlist)
1188 print_log()
1189 return
1190
1191 def add(self, offset, insitu=None):
1192 """
1193 Return a scan where all spectra have the offset added
1194 Parameters:
1195 offset: the offset
1196 insitu: if False a new scantable is returned.
1197 Otherwise, the scaling is done in-situ
1198 The default is taken from .asaprc (False)
1199 """
1200 if insitu is None: insitu = rcParams['insitu']
1201 self._math._setinsitu(insitu)
1202 varlist = vars()
1203 s = scantable(self._math._unaryop(self, offset, "ADD", False))
1204 s._add_history("add",varlist)
1205 print_log()
1206 if insitu:
1207 self._assign(s)
1208 else:
1209 return s
1210
1211 def scale(self, factor, tsys=True, insitu=None,):
1212 """
1213 Return a scan where all spectra are scaled by the give 'factor'
1214 Parameters:
1215 factor: the scaling factor
1216 insitu: if False a new scantable is returned.
1217 Otherwise, the scaling is done in-situ
1218 The default is taken from .asaprc (False)
1219 tsys: if True (default) then apply the operation to Tsys
1220 as well as the data
1221 """
1222 if insitu is None: insitu = rcParams['insitu']
1223 self._math._setinsitu(insitu)
1224 varlist = vars()
1225 s = scantable(self._math._unaryop(self, factor, "MUL", tsys))
1226 s._add_history("scale",varlist)
1227 print_log()
1228 if insitu:
1229 self._assign(s)
1230 else:
1231 return s
1232
1233 def auto_quotient(self, mode='time', preserve=True):
1234 """
1235 This function allows to build quotients automatically.
1236 It assumes the observation to have the same numer of
1237 "ons" and "offs"
1238 It will support "closest off in time" in the future
1239 Parameters:
1240 mode: the on/off detection mode; 'suffix' (default)
1241 'suffix' identifies 'off' scans by the
1242 trailing '_R' (Mopra/Parkes) or
1243 '_e'/'_w' (Tid)
1244 preserve: you can preserve (default) the continuum or
1245 remove it. The equations used are
1246 preserve: Output = Toff * (on/off) - Toff
1247 remove: Output = Tref * (on/off) - Ton
1248 """
1249 modes = ["time"]
1250 if not mode in modes:
1251 msg = "please provide valid mode. Valid modes are %s" % (modes)
1252 raise ValueError(msg)
1253 varlist = vars()
1254 s = scantable(self._math._quotient(self, mode, preserve))
1255 s._add_history("auto_quotient",varlist)
1256 print_log()
1257 return s
1258
1259
1260
1261
1262 def freq_switch(self, insitu=None):
1263 """
1264 Apply frequency switching to the data.
1265 Parameters:
1266 insitu: if False a new scantable is returned.
1267 Otherwise, the swictching is done in-situ
1268 The default is taken from .asaprc (False)
1269 Example:
1270 none
1271 """
1272 if insitu is None: insitu = rcParams['insitu']
1273 self._math._setinsitu(insitu)
1274 varlist = vars()
1275 s = scantable(self._math._freqswitch(self))
1276 s._add_history("freq_switch",varlist)
1277 print_log()
1278 if insitu: self._assign(s)
1279 else: return s
1280
1281 def recalc_azel(self):
1282 """
1283 Recalculate the azimuth and elevation for each position.
1284 Parameters:
1285 none
1286 Example:
1287 """
1288 varlist = vars()
1289 self._recalcazel()
1290 self._add_history("recalc_azel", varlist)
1291 print_log()
1292 return
1293
1294 def __add__(self, other):
1295 varlist = vars()
1296 s = None
1297 if isinstance(other, scantable):
1298 print "scantable + scantable NYI"
1299 return
1300 elif isinstance(other, float):
1301 s = scantable(self._math._unaryop(self, other, "ADD", False))
1302 else:
1303 raise TypeError("Other input is not a scantable or float value")
1304 s._add_history("operator +", varlist)
1305 print_log()
1306 return s
1307
1308 def __sub__(self, other):
1309 """
1310 implicit on all axes and on Tsys
1311 """
1312 varlist = vars()
1313 s = None
1314 if isinstance(other, scantable):
1315 print "scantable - scantable NYI"
1316 return
1317 elif isinstance(other, float):
1318 s = scantable(self._math._unaryop(self, other, "SUB", False))
1319 else:
1320 raise TypeError("Other input is not a scantable or float value")
1321 s._add_history("operator -", varlist)
1322 print_log()
1323 return s
1324
1325 def __mul__(self, other):
1326 """
1327 implicit on all axes and on Tsys
1328 """
1329 varlist = vars()
1330 s = None
1331 if isinstance(other, scantable):
1332 print "scantable * scantable NYI"
1333 return
1334 elif isinstance(other, float):
1335 s = scantable(self._math._unaryop(self, other, "MUL", False))
1336 else:
1337 raise TypeError("Other input is not a scantable or float value")
1338 s._add_history("operator *", varlist)
1339 print_log()
1340 return s
1341
1342
1343 def __div__(self, other):
1344 """
1345 implicit on all axes and on Tsys
1346 """
1347 varlist = vars()
1348 s = None
1349 if isinstance(other, scantable):
1350 print "scantable / scantable NYI"
1351 return
1352 elif isinstance(other, float):
1353 if other == 0.0:
1354 raise ZeroDivisionError("Dividing by zero is not recommended")
1355 s = scantable(self._math._unaryop(self, other, "DIV", False))
1356 else:
1357 raise TypeError("Other input is not a scantable or float value")
1358 s._add_history("operator /", varlist)
1359 print_log()
1360 return s
1361
1362 def get_fit(self, row=0):
1363 """
1364 Print or return the stored fits for a row in the scantable
1365 Parameters:
1366 row: the row which the fit has been applied to.
1367 """
1368 if row > self.nrow():
1369 return
1370 from asap.asapfit import asapfit
1371 fit = asapfit(self._getfit(row))
1372 if rcParams['verbose']:
1373 print fit
1374 return
1375 else:
1376 return fit.as_dict()
1377
1378 def _add_history(self, funcname, parameters):
1379 # create date
1380 sep = "##"
1381 from datetime import datetime
1382 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
1383 hist = dstr+sep
1384 hist += funcname+sep#cdate+sep
1385 if parameters.has_key('self'): del parameters['self']
1386 for k,v in parameters.iteritems():
1387 if type(v) is dict:
1388 for k2,v2 in v.iteritems():
1389 hist += k2
1390 hist += "="
1391 if isinstance(v2,scantable):
1392 hist += 'scantable'
1393 elif k2 == 'mask':
1394 if isinstance(v2,list) or isinstance(v2,tuple):
1395 hist += str(self._zip_mask(v2))
1396 else:
1397 hist += str(v2)
1398 else:
1399 hist += str(v2)
1400 else:
1401 hist += k
1402 hist += "="
1403 if isinstance(v,scantable):
1404 hist += 'scantable'
1405 elif k == 'mask':
1406 if isinstance(v,list) or isinstance(v,tuple):
1407 hist += str(self._zip_mask(v))
1408 else:
1409 hist += str(v)
1410 else:
1411 hist += str(v)
1412 hist += sep
1413 hist = hist[:-2] # remove trailing '##'
1414 self._addhistory(hist)
1415
1416
1417 def _zip_mask(self, mask):
1418 mask = list(mask)
1419 i = 0
1420 segments = []
1421 while mask[i:].count(1):
1422 i += mask[i:].index(1)
1423 if mask[i:].count(0):
1424 j = i + mask[i:].index(0)
1425 else:
1426 j = len(mask)
1427 segments.append([i,j])
1428 i = j
1429 return segments
1430
1431 def _get_ordinate_label(self):
1432 fu = "("+self.get_fluxunit()+")"
1433 import re
1434 lbl = "Intensity"
1435 if re.match(".K.",fu):
1436 lbl = "Brightness Temperature "+ fu
1437 elif re.match(".Jy.",fu):
1438 lbl = "Flux density "+ fu
1439 return lbl
1440
1441 def _check_ifs(self):
1442 nchans = [self.nchan(i) for i in range(self.nif(-1))]
1443 nchans = filter(lambda t: t > 0, nchans)
1444 return (sum(nchans)/len(nchans) == nchans[0])
1445
1446 def _fill(self, names, unit, average):
1447 import os
1448 varlist = vars()
1449 from asap._asap import stfiller
1450 first = True
1451 fullnames = []
1452 for name in names:
1453 name = os.path.expandvars(name)
1454 name = os.path.expanduser(name)
1455 if not os.path.exists(name):
1456 msg = "File '%s' does not exists" % (name)
1457 if rcParams['verbose']:
1458 asaplog.push(msg)
1459 print asaplog.pop().strip()
1460 return
1461 raise IOError(msg)
1462 fullnames.append(name)
1463 if average:
1464 asaplog.push('Auto averaging integrations')
1465 for name in fullnames:
1466 r = stfiller()
1467 msg = "Importing %s..." % (name)
1468 asaplog.push(msg,False)
1469 print_log()
1470 r._open(name,-1,-1)
1471 r._read()
1472 tbl = r._getdata()
1473 if average:
1474 tbl = self._math._average((tbl,),(),'NONE','SCAN')
1475 #tbl = tbl2
1476 if not first:
1477 tbl = self._math._merge([self, tbl])
1478 #tbl = tbl2
1479 Scantable.__init__(self, tbl)
1480 r._close()
1481 del r,tbl
1482 first = False
1483 if unit is not None:
1484 self.set_fluxunit(unit)
1485 self.set_freqframe(rcParams['scantable.freqframe'])
1486 #self._add_history("scantable", varlist)
1487
Note: See TracBrowser for help on using the repository browser.