source: trunk/python/scantable.py@ 938

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

fixed self.stm to self._math

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