source: trunk/python/scantable.py@ 964

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

re-enabled average_pol; made use of selector class

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