source: branches/alma/python/scantable.py@ 1613

Last change on this file since 1613 was 1613, checked in by Takeshi Nakazato, 16 years ago

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Temporary file name has been changed.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 80.0 KB
Line 
1from asap._asap import Scantable
2from asap import rcParams
3from asap import print_log
4from asap import asaplog
5from asap import selector
6from asap import linecatalog
7from asap import _n_bools, mask_not, mask_and, mask_or
8from taskinit import *
9
10class scantable(Scantable):
11 """
12 The ASAP container for scans
13 """
14
15 def __init__(self, filename, average=None, unit=None, getpt=None):
16 """
17 Create a scantable from a saved one or make a reference
18 Parameters:
19 filename: the name of an asap table on disk
20 or
21 the name of a rpfits/sdfits/ms file
22 (integrations within scans are auto averaged
23 and the whole file is read)
24 or
25 [advanced] a reference to an existing
26 scantable
27 average: average all integrations withinb a scan on read.
28 The default (True) is taken from .asaprc.
29 unit: brightness unit; must be consistent with K or Jy.
30 Over-rides the default selected by the reader
31 (input rpfits/sdfits/ms) or replaces the value
32 in existing scantables
33 getpt: for MeasurementSet input data only:
34 If True, all pointing data are filled.
35 The deafult is False, which makes time to load
36 the MS data faster in some cases.
37 """
38 if average is None:
39 average = rcParams['scantable.autoaverage']
40 if getpt is None:
41 getpt = False
42 varlist = vars()
43 from asap._asap import stmath
44 self._math = stmath()
45 if isinstance(filename, Scantable):
46 Scantable.__init__(self, filename)
47 else:
48 if isinstance(filename, str):# or \
49# (isinstance(filename, list) or isinstance(filename, tuple)) \
50# and isinstance(filename[-1], str):
51 import os.path
52 filename = os.path.expandvars(filename)
53 filename = os.path.expanduser(filename)
54 if not os.path.exists(filename):
55 s = "File '%s' not found." % (filename)
56 if rcParams['verbose']:
57 asaplog.push(s)
58 #print asaplog.pop().strip()
59 print_log()
60 return
61 raise IOError(s)
62 if os.path.isdir(filename) \
63 and not os.path.exists(filename+'/table.f1'):
64 # crude check if asap table
65 if os.path.exists(filename+'/table.info'):
66 ondisk = rcParams['scantable.storage'] == 'disk'
67 Scantable.__init__(self, filename, ondisk)
68 if unit is not None:
69 self.set_fluxunit(unit)
70 # do not reset to the default freqframe
71 #self.set_freqframe(rcParams['scantable.freqframe'])
72 else:
73 msg = "The given file '%s'is not a valid " \
74 "asap table." % (filename)
75 if rcParams['verbose']:
76 #print msg
77 casalog.post( msg, 'WARN' )
78 return
79 else:
80 raise IOError(msg)
81 else:
82 self._fill([filename], unit, average, getpt)
83 elif (isinstance(filename, list) or isinstance(filename, tuple)) \
84 and isinstance(filename[-1], str):
85 self._fill(filename, unit, average, getpt)
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,
92 SDFITS or MS2 format.
93 Parameters:
94 name: the name of the outputfile. For format "ASCII"
95 this is the root file name (data in 'name'.txt
96 and header in 'name'_header.txt)
97 format: an optional file format. Default is ASAP.
98 Allowed are - 'ASAP' (save as ASAP [aips++] Table),
99 'SDFITS' (save as SDFITS file)
100 'ASCII' (saves as ascii text file)
101 'MS2' (saves as an aips++
102 MeasurementSet V2)
103 'FITS' (save as image FITS - not
104 readable by class)
105 'CLASS' (save as FITS readable by CLASS)
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 msg = "No filename given. Using default name %s..." % name
119 asaplog.push(msg)
120 name = path.expandvars(name)
121 if path.isfile(name) or path.isdir(name):
122 if not overwrite:
123 msg = "File %s exists." % name
124 if rcParams['verbose']:
125 #print msg
126 casalog.post( msg, 'WARN' )
127 return
128 else:
129 raise IOError(msg)
130 format2 = format.upper()
131 if format2 == 'ASAP':
132 self._save(name)
133 else:
134 from asap._asap import stwriter as stw
135 writer = stw(format2)
136 writer.write(self, name)
137 print_log()
138 return
139
140 def copy(self):
141 """
142 Return a copy of this scantable.
143 Note:
144 This makes a full (deep) copy. scan2 = scan1 makes a reference.
145 Parameters:
146 none
147 Example:
148 copiedscan = scan.copy()
149 """
150 sd = scantable(Scantable._copy(self))
151 return sd
152
153 def drop_scan(self, scanid=None):
154 """
155 Return a new scantable where the specified scan number(s) has(have)
156 been dropped.
157 Parameters:
158 scanid: a (list of) scan number(s)
159 """
160 from asap import _is_sequence_or_number as _is_valid
161 from asap import _to_list
162 from asap import unique
163 if not _is_valid(scanid):
164 if rcParams['verbose']:
165 #print "Please specify a scanno to drop from the scantable"
166 casalog.post( "Please specify a scanno to drop from the scantable", 'WARN' )
167 return
168 else:
169 raise RuntimeError("No scan given")
170 try:
171 scanid = _to_list(scanid)
172 allscans = unique([ self.getscan(i) for i in range(self.nrow())])
173 for sid in scanid: allscans.remove(sid)
174 if len(allscans) == 0:
175 raise ValueError("Can't remove all scans")
176 except ValueError:
177 if rcParams['verbose']:
178 #print "Couldn't find any match."
179 casalog.post( "Couldn't find any match.", 'WARN' )
180 return
181 else: raise
182 try:
183 bsel = self.get_selection()
184 sel = selector()
185 sel.set_scans(allscans)
186 self.set_selection(bsel+sel)
187 scopy = self._copy()
188 self.set_selection(bsel)
189 return scantable(scopy)
190 except RuntimeError:
191 if rcParams['verbose']:
192 #print "Couldn't find any match."
193 casalog.post( "Couldn't find any match.", 'WARN' )
194 else:
195 raise
196
197
198 def get_scan(self, scanid=None):
199 """
200 Return a specific scan (by scanno) or collection of scans (by
201 source name) in a new scantable.
202 Note:
203 See scantable.drop_scan() for the inverse operation.
204 Parameters:
205 scanid: a (list of) scanno or a source name, unix-style
206 patterns are accepted for source name matching, e.g.
207 '*_R' gets all 'ref scans
208 Example:
209 # get all scans containing the source '323p459'
210 newscan = scan.get_scan('323p459')
211 # get all 'off' scans
212 refscans = scan.get_scan('*_R')
213 # get a susbset of scans by scanno (as listed in scan.summary())
214 newscan = scan.get_scan([0, 2, 7, 10])
215 """
216 if scanid is None:
217 if rcParams['verbose']:
218 #print "Please specify a scan no or name to " \
219 # "retrieve from the scantable"
220 casalog.post( "Please specify a scan no or name to retrieve from the scantable", 'WARN' )
221 return
222 else:
223 raise RuntimeError("No scan given")
224
225 try:
226 bsel = self.get_selection()
227 sel = selector()
228 if type(scanid) is str:
229 sel.set_name(scanid)
230 self.set_selection(bsel+sel)
231 scopy = self._copy()
232 self.set_selection(bsel)
233 return scantable(scopy)
234 elif type(scanid) is int:
235 sel.set_scans([scanid])
236 self.set_selection(bsel+sel)
237 scopy = self._copy()
238 self.set_selection(bsel)
239 return scantable(scopy)
240 elif type(scanid) is list:
241 sel.set_scans(scanid)
242 self.set_selection(sel)
243 scopy = self._copy()
244 self.set_selection(bsel)
245 return scantable(scopy)
246 else:
247 msg = "Illegal scanid type, use 'int' or 'list' if ints."
248 if rcParams['verbose']:
249 #print msg
250 casalog.post( msg, 'WARN' )
251 else:
252 raise TypeError(msg)
253 except RuntimeError:
254 if rcParams['verbose']:
255 #print "Couldn't find any match."
256 casalog.post( "Couldn't find any match.", 'WARN' )
257 else: raise
258
259 def __str__(self):
260 return Scantable._summary(self, True)
261
262 def summary(self, filename=None):
263 """
264 Print a summary of the contents of this scantable.
265 Parameters:
266 filename: the name of a file to write the putput to
267 Default - no file output
268 verbose: print extra info such as the frequency table
269 The default (False) is taken from .asaprc
270 """
271 info = Scantable._summary(self, True)
272 #if verbose is None: verbose = rcParams['scantable.verbosesummary']
273 if filename is not None:
274 if filename is "":
275 filename = 'scantable_summary.txt'
276 from os.path import expandvars, isdir
277 filename = expandvars(filename)
278 if not isdir(filename):
279 data = open(filename, 'w')
280 data.write(info)
281 data.close()
282 else:
283 msg = "Illegal file name '%s'." % (filename)
284 if rcParams['verbose']:
285 #print msg
286 casalog.post( msg, 'WARN' )
287 else:
288 raise IOError(msg)
289 if rcParams['verbose']:
290 try:
291 from IPython.genutils import page as pager
292 except ImportError:
293 from pydoc import pager
294 pager(info)
295 else:
296 return info
297
298 def get_spectrum(self, rowno):
299 """Return the spectrum for the current row in the scantable as a list.
300 Parameters:
301 rowno: the row number to retrieve the spectrum from
302 """
303 return self._getspectrum(rowno)
304
305 def get_mask(self, rowno):
306 """Return the mask for the current row in the scantable as a list.
307 Parameters:
308 rowno: the row number to retrieve the mask from
309 """
310 return self._getmask(rowno)
311
312 def set_spectrum(self, spec, rowno):
313 """Return the spectrum for the current row in the scantable as a list.
314 Parameters:
315 spec: the spectrum
316 rowno: the row number to set the spectrum for
317 """
318 assert(len(spec) == self.nchan())
319 return self._setspectrum(spec, rowno)
320
321 def get_selection(self):
322 """
323 Get the selection object currently set on this scantable.
324 Parameters:
325 none
326 Example:
327 sel = scan.get_selection()
328 sel.set_ifs(0) # select IF 0
329 scan.set_selection(sel) # apply modified selection
330 """
331 return selector(self._getselection())
332
333 def set_selection(self, selection=selector()):
334 """
335 Select a subset of the data. All following operations on this scantable
336 are only applied to thi selection.
337 Parameters:
338 selection: a selector object (default unset the selection)
339 Examples:
340 sel = selector() # create a selection object
341 self.set_scans([0, 3]) # select SCANNO 0 and 3
342 scan.set_selection(sel) # set the selection
343 scan.summary() # will only print summary of scanno 0 an 3
344 scan.set_selection() # unset the selection
345 """
346 self._setselection(selection)
347
348 def get_row(self, row=0, insitu=None):
349 """
350 Select a row in the scantable.
351 Return a scantable with single row.
352 Parameters:
353 row: row no of integration, default is 0.
354 insitu: if False a new scantable is returned.
355 Otherwise, the scaling is done in-situ
356 The default is taken from .asaprc (False)
357 """
358 if insitu is None: insitu = rcParams['insitu']
359 if not insitu:
360 workscan = self.copy()
361 else:
362 workscan = self
363 # Select a row
364 sel=selector()
365 sel.set_scans([workscan.getscan(row)])
366 sel.set_cycles([workscan.getcycle(row)])
367 sel.set_beams([workscan.getbeam(row)])
368 sel.set_ifs([workscan.getif(row)])
369 sel.set_polarisations([workscan.getpol(row)])
370 sel.set_name(workscan._getsourcename(row))
371 workscan.set_selection(sel)
372 if not workscan.nrow() == 1:
373 msg = "Cloud not identify single row. %d rows selected."%(workscan.nrow())
374 raise RuntimeError(msg)
375 del sel
376 if insitu:
377 self._assign(workscan)
378 else:
379 return workscan
380
381 def stats(self, stat='stddev', mask=None):
382 """
383 Determine the specified statistic of the current beam/if/pol
384 Takes a 'mask' as an optional parameter to specify which
385 channels should be excluded.
386 Parameters:
387 stat: 'min', 'max', 'min_abc', 'max_abc', 'sumsq', 'sum',
388 'mean', 'var', 'stddev', 'avdev', 'rms', 'median'
389 mask: an optional mask specifying where the statistic
390 should be determined.
391 Example:
392 scan.set_unit('channel')
393 msk = scan.create_mask([100, 200], [500, 600])
394 scan.stats(stat='mean', mask=m)
395 """
396 if mask == None:
397 mask = []
398 axes = ['Beam', 'IF', 'Pol', 'Time']
399 if not self._check_ifs():
400 raise ValueError("Cannot apply mask as the IFs have different "
401 "number of channels. Please use setselection() "
402 "to select individual IFs")
403 rtnabc = False
404 if stat.lower().endswith('_abc'): rtnabc = True
405 getchan = False
406 if stat.lower().startswith('min') or stat.lower().startswith('max'):
407 chan = self._math._minmaxchan(self, mask, stat)
408 getchan = True
409 statvals = []
410 if not rtnabc: statvals = self._math._stats(self, mask, stat)
411
412 out = ''
413 axes = []
414 for i in range(self.nrow()):
415 axis = []
416 axis.append(self.getscan(i))
417 axis.append(self.getbeam(i))
418 axis.append(self.getif(i))
419 axis.append(self.getpol(i))
420 axis.append(self.getcycle(i))
421 axes.append(axis)
422 tm = self._gettime(i)
423 src = self._getsourcename(i)
424 refstr = ''
425 statunit= ''
426 if getchan:
427 qx, qy = self.chan2data(rowno=i, chan=chan[i])
428 if rtnabc:
429 statvals.append(qx['value'])
430 refstr = '(value: %3.3f' % (qy['value'])+' ['+qy['unit']+'])'
431 statunit= '['+qx['unit']+']'
432 else:
433 refstr = '(@ %3.3f' % (qx['value'])+' ['+qx['unit']+'])'
434 #statunit= ' ['+qy['unit']+']'
435 out += 'Scan[%d] (%s) ' % (axis[0], src)
436 out += 'Time[%s]:\n' % (tm)
437 if self.nbeam(-1) > 1: out += ' Beam[%d] ' % (axis[1])
438 if self.nif(-1) > 1: out += ' IF[%d] ' % (axis[2])
439 if self.npol(-1) > 1: out += ' Pol[%d] ' % (axis[3])
440 out += '= %3.3f ' % (statvals[i]) +refstr+'\n'
441 out += "--------------------------------------------------\n"
442
443 if rcParams['verbose']:
444 usr=os.environ['USER']
445 tmpfile='/tmp/tmp_'+usr+'_casapy_asap_scantable_stats'
446 f=open(tmpfile,'w')
447 print >> f, "--------------------------------------------------"
448 print >> f, " ", stat, statunit
449 print >> f, "--------------------------------------------------"
450 print >> f, out
451 f.close()
452 f=open(tmpfile,'r')
453 x=f.readlines()
454 f.close()
455 for xx in x:
456 casalog.post( xx )
457 #else:
458 #retval = { 'axesnames': ['scanno', 'beamno', 'ifno', 'polno', 'cycleno'],
459 # 'axes' : axes,
460 # 'data': statvals}
461 return statvals
462
463 def chan2data(self, rowno=0, chan=0):
464 """
465 Returns channel/frequency/velocity and spectral value
466 at an arbitrary row and channel in the scantable.
467 Parameters:
468 rowno: a row number in the scantable. Default is the
469 first row, i.e. rowno=0
470 chan: a channel in the scantable. Default is the first
471 channel, i.e. pos=0
472 """
473 if isinstance(rowno, int) and isinstance(chan, int):
474 qx = {'unit': self.get_unit(),
475 'value': self._getabcissa(rowno)[chan]}
476 qy = {'unit': self.get_fluxunit(),
477 'value': self._getspectrum(rowno)[chan]}
478 return qx, qy
479
480 def stddev(self, mask=None):
481 """
482 Determine the standard deviation of the current beam/if/pol
483 Takes a 'mask' as an optional parameter to specify which
484 channels should be excluded.
485 Parameters:
486 mask: an optional mask specifying where the standard
487 deviation should be determined.
488
489 Example:
490 scan.set_unit('channel')
491 msk = scan.create_mask([100, 200], [500, 600])
492 scan.stddev(mask=m)
493 """
494 return self.stats(stat='stddev', mask=mask);
495
496
497 def get_column_names(self):
498 """
499 Return a list of column names, which can be used for selection.
500 """
501 return list(Scantable.get_column_names(self))
502
503 def get_tsys(self):
504 """
505 Return the System temperatures.
506 Returns:
507 a list of Tsys values for the current selection
508 """
509
510 return self._row_callback(self._gettsys, "Tsys")
511
512 def _row_callback(self, callback, label):
513 axes = []
514 axesnames = ['scanno', 'beamno', 'ifno', 'polno', 'cycleno']
515 out = ""
516 outvec = []
517 for i in range(self.nrow()):
518 axis = []
519 axis.append(self.getscan(i))
520 axis.append(self.getbeam(i))
521 axis.append(self.getif(i))
522 axis.append(self.getpol(i))
523 axis.append(self.getcycle(i))
524 axes.append(axis)
525 tm = self._gettime(i)
526 src = self._getsourcename(i)
527 out += 'Scan[%d] (%s) ' % (axis[0], src)
528 out += 'Time[%s]:\n' % (tm)
529 if self.nbeam(-1) > 1: out += ' Beam[%d] ' % (axis[1])
530 if self.nif(-1) > 1: out += ' IF[%d] ' % (axis[2])
531 if self.npol(-1) > 1: out += ' Pol[%d] ' % (axis[3])
532 outvec.append(callback(i))
533 out += '= %3.3f\n' % (outvec[i])
534 out += "--------------------------------------------------\n"
535 if rcParams['verbose']:
536 #print "--------------------------------------------------"
537 #print " %s" % (label)
538 #print "--------------------------------------------------"
539 #print out
540 casalog.post( "--------------------------------------------------" )
541 casalog.post( " %s" % (label) )
542 casalog.post( "--------------------------------------------------" )
543 casalog.post( out )
544 # disabled because the vector seems more useful
545 #retval = {'axesnames': axesnames, 'axes': axes, 'data': outvec}
546 return outvec
547
548 def _get_column(self, callback, row=-1):
549 """
550 """
551 if row == -1:
552 return [callback(i) for i in range(self.nrow())]
553 else:
554 if 0 <= row < self.nrow():
555 return callback(row)
556
557
558 def get_time(self, row=-1, asdatetime=False):
559 """
560 Get a list of time stamps for the observations.
561 Return a datetime object for each integration time stamp in the scantable.
562 Parameters:
563 row: row no of integration. Default -1 return all rows
564 asdatetime: return values as datetime objects rather than strings
565 Example:
566 none
567 """
568 from time import strptime
569 from datetime import datetime
570 times = self._get_column(self._gettime, row)
571 if not asdatetime:
572 return times
573 format = "%Y/%m/%d/%H:%M:%S"
574 if isinstance(times, list):
575 return [datetime(*strptime(i, format)[:6]) for i in times]
576 else:
577 return datetime(*strptime(times, format)[:6])
578
579
580 def get_inttime(self, row=-1):
581 """
582 Get a list of integration times for the observations.
583 Return a time in seconds for each integration in the scantable.
584 Parameters:
585 row: row no of integration. Default -1 return all rows.
586 Example:
587 none
588 """
589 return self._get_column(self._getinttime, row)
590
591
592 def get_sourcename(self, row=-1):
593 """
594 Get a list source names for the observations.
595 Return a string for each integration in the scantable.
596 Parameters:
597 row: row no of integration. Default -1 return all rows.
598 Example:
599 none
600 """
601 return self._get_column(self._getsourcename, row)
602
603 def get_elevation(self, row=-1):
604 """
605 Get a list of elevations for the observations.
606 Return a float for each integration in the scantable.
607 Parameters:
608 row: row no of integration. Default -1 return all rows.
609 Example:
610 none
611 """
612 return self._get_column(self._getelevation, row)
613
614 def get_azimuth(self, row=-1):
615 """
616 Get a list of azimuths for the observations.
617 Return a float for each integration in the scantable.
618 Parameters:
619 row: row no of integration. Default -1 return all rows.
620 Example:
621 none
622 """
623 return self._get_column(self._getazimuth, row)
624
625 def get_parangle(self, row=-1):
626 """
627 Get a list of parallactic angles for the observations.
628 Return a float for each integration in the scantable.
629 Parameters:
630 row: row no of integration. Default -1 return all rows.
631 Example:
632 none
633 """
634 return self._get_column(self._getparangle, row)
635
636 def get_direction(self, row=-1):
637 """
638 Get a list of Positions on the sky (direction) for the observations.
639 Return a float for each integration in the scantable.
640 Parameters:
641 row: row no of integration. Default -1 return all rows
642 Example:
643 none
644 """
645 return self._get_column(self._getdirection, row)
646
647 def get_directionval(self, row=-1):
648 """
649 Get a list of Positions on the sky (direction) for the observations.
650 Return a float for each integration in the scantable.
651 Parameters:
652 row: row no of integration. Default -1 return all rows
653 Example:
654 none
655 """
656 return self._get_column(self._getdirectionvec, row)
657
658 def set_unit(self, unit='channel'):
659 """
660 Set the unit for all following operations on this scantable
661 Parameters:
662 unit: optional unit, default is 'channel'
663 one of '*Hz', 'km/s', 'channel', ''
664 """
665 varlist = vars()
666 if unit in ['', 'pixel', 'channel']:
667 unit = ''
668 inf = list(self._getcoordinfo())
669 inf[0] = unit
670 self._setcoordinfo(inf)
671 self._add_history("set_unit", varlist)
672
673 def set_instrument(self, instr):
674 """
675 Set the instrument for subsequent processing.
676 Parameters:
677 instr: Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
678 'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
679 """
680 self._setInstrument(instr)
681 self._add_history("set_instument", vars())
682 print_log()
683
684 def set_feedtype(self, feedtype):
685 """
686 Overwrite the feed type, which might not be set correctly.
687 Parameters:
688 feedtype: 'linear' or 'circular'
689 """
690 self._setfeedtype(feedtype)
691 self._add_history("set_feedtype", vars())
692 print_log()
693
694 def set_doppler(self, doppler='RADIO'):
695 """
696 Set the doppler for all following operations on this scantable.
697 Parameters:
698 doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
699 """
700 varlist = vars()
701 inf = list(self._getcoordinfo())
702 inf[2] = doppler
703 self._setcoordinfo(inf)
704 self._add_history("set_doppler", vars())
705 print_log()
706
707 def set_freqframe(self, frame=None):
708 """
709 Set the frame type of the Spectral Axis.
710 Parameters:
711 frame: an optional frame type, default 'LSRK'. Valid frames are:
712 'REST', 'TOPO', 'LSRD', 'LSRK', 'BARY',
713 'GEO', 'GALACTO', 'LGROUP', 'CMB'
714 Examples:
715 scan.set_freqframe('BARY')
716 """
717 if frame is None: frame = rcParams['scantable.freqframe']
718 varlist = vars()
719 valid = ['REST', 'TOPO', 'LSRD', 'LSRK', 'BARY', \
720 'GEO', 'GALACTO', 'LGROUP', 'CMB']
721
722 if frame in valid:
723 inf = list(self._getcoordinfo())
724 inf[1] = frame
725 self._setcoordinfo(inf)
726 self._add_history("set_freqframe", varlist)
727 else:
728 msg = "Please specify a valid freq type. Valid types are:\n", valid
729 if rcParams['verbose']:
730 #print msg
731 casalog.post( msg, 'WARN' )
732 else:
733 raise TypeError(msg)
734 print_log()
735
736 def set_dirframe(self, frame=""):
737 """
738 Set the frame type of the Direction on the sky.
739 Parameters:
740 frame: an optional frame type, default ''. Valid frames are:
741 'J2000', 'B1950', 'GALACTIC'
742 Examples:
743 scan.set_dirframe('GALACTIC')
744 """
745 varlist = vars()
746 try:
747 Scantable.set_dirframe(self, frame)
748 except RuntimeError, msg:
749 if rcParams['verbose']:
750 #print msg
751 casalog.post( msg, 'WARN' )
752 else:
753 raise
754 self._add_history("set_dirframe", varlist)
755
756 def get_unit(self):
757 """
758 Get the default unit set in this scantable
759 Returns:
760 A unit string
761 """
762 inf = self._getcoordinfo()
763 unit = inf[0]
764 if unit == '': unit = 'channel'
765 return unit
766
767 def get_abcissa(self, rowno=0):
768 """
769 Get the abcissa in the current coordinate setup for the currently
770 selected Beam/IF/Pol
771 Parameters:
772 rowno: an optional row number in the scantable. Default is the
773 first row, i.e. rowno=0
774 Returns:
775 The abcissa values and the format string (as a dictionary)
776 """
777 abc = self._getabcissa(rowno)
778 lbl = self._getabcissalabel(rowno)
779 print_log()
780 return abc, lbl
781
782 def flag(self, mask=None, unflag=False):
783 """
784 Flag the selected data using an optional channel mask.
785 Parameters:
786 mask: an optional channel mask, created with create_mask. Default
787 (no mask) is all channels.
788 unflag: if True, unflag the data
789 """
790 varlist = vars()
791 if mask is None:
792 mask = []
793 try:
794 self._flag(mask, unflag)
795 except RuntimeError, msg:
796 if rcParams['verbose']:
797 #print msg
798 casalog.post( msg, 'WARN' )
799 return
800 else: raise
801 self._add_history("flag", varlist)
802
803 def lag_flag(self, frequency, width=0.0, unit="GHz", insitu=None):
804 """
805 Flag the data in 'lag' space by providing a frequency to remove.
806 Flagged data in the scantable gets set to 0.0 before the fft.
807 No taper is applied.
808 Parameters:
809 frequency: the frequency (really a period within the bandwidth)
810 to remove
811 width: the width of the frequency to remove, to remove a
812 range of frequencies around the centre.
813 unit: the frequency unit (default "GHz")
814 Notes:
815 It is recommended to flag edges of the band or strong
816 signals beforehand.
817 """
818 if insitu is None: insitu = rcParams['insitu']
819 self._math._setinsitu(insitu)
820 varlist = vars()
821 base = { "GHz": 1000000000., "MHz": 1000000., "kHz": 1000., "Hz": 1. }
822 if not base.has_key(unit):
823 raise ValueError("%s is not a valid unit." % unit)
824 try:
825 s = scantable(self._math._lag_flag(self, frequency*base[unit],
826 width*base[unit]))
827 except RuntimeError, msg:
828 if rcParams['verbose']:
829 #print msg
830 casalog.post( msg, 'WARN' )
831 return
832 else: raise
833 s._add_history("lag_flag", varlist)
834 print_log()
835 if insitu:
836 self._assign(s)
837 else:
838 return s
839
840
841 def create_mask(self, *args, **kwargs):
842 """
843 Compute and return a mask based on [min, max] windows.
844 The specified windows are to be INCLUDED, when the mask is
845 applied.
846 Parameters:
847 [min, max], [min2, max2], ...
848 Pairs of start/end points (inclusive)specifying the regions
849 to be masked
850 invert: optional argument. If specified as True,
851 return an inverted mask, i.e. the regions
852 specified are EXCLUDED
853 row: create the mask using the specified row for
854 unit conversions, default is row=0
855 only necessary if frequency varies over rows.
856 Example:
857 scan.set_unit('channel')
858 a)
859 msk = scan.create_mask([400, 500], [800, 900])
860 # masks everything outside 400 and 500
861 # and 800 and 900 in the unit 'channel'
862
863 b)
864 msk = scan.create_mask([400, 500], [800, 900], invert=True)
865 # masks the regions between 400 and 500
866 # and 800 and 900 in the unit 'channel'
867 c)
868 mask only channel 400
869 msk = scan.create_mask([400, 400])
870 """
871 row = 0
872 if kwargs.has_key("row"):
873 row = kwargs.get("row")
874 data = self._getabcissa(row)
875 u = self._getcoordinfo()[0]
876 if rcParams['verbose']:
877 if u == "": u = "channel"
878 msg = "The current mask window unit is %s" % u
879 i = self._check_ifs()
880 if not i:
881 msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
882 asaplog.push(msg)
883 n = self.nchan()
884 msk = _n_bools(n, False)
885 # test if args is a 'list' or a 'normal *args - UGLY!!!
886
887 ws = (isinstance(args[-1][-1], int) or isinstance(args[-1][-1], float)) \
888 and args or args[0]
889 for window in ws:
890 if (len(window) != 2 or window[0] > window[1] ):
891 raise TypeError("A window needs to be defined as [min, max]")
892 for i in range(n):
893 if data[i] >= window[0] and data[i] <= window[1]:
894 msk[i] = True
895 if kwargs.has_key('invert'):
896 if kwargs.get('invert'):
897 msk = mask_not(msk)
898 print_log()
899 return msk
900
901 def get_masklist(self, mask=None, row=0):
902 """
903 Compute and return a list of mask windows, [min, max].
904 Parameters:
905 mask: channel mask, created with create_mask.
906 row: calcutate the masklist using the specified row
907 for unit conversions, default is row=0
908 only necessary if frequency varies over rows.
909 Returns:
910 [min, max], [min2, max2], ...
911 Pairs of start/end points (inclusive)specifying
912 the masked regions
913 """
914 if not (isinstance(mask,list) or isinstance(mask, tuple)):
915 raise TypeError("The mask should be list or tuple.")
916 if len(mask) < 2:
917 raise TypeError("The mask elements should be > 1")
918 if self.nchan() != len(mask):
919 msg = "Number of channels in scantable != number of mask elements"
920 raise TypeError(msg)
921 data = self._getabcissa(row)
922 u = self._getcoordinfo()[0]
923 if rcParams['verbose']:
924 if u == "": u = "channel"
925 msg = "The current mask window unit is %s" % u
926 i = self._check_ifs()
927 if not i:
928 msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
929 asaplog.push(msg)
930 masklist=[]
931 ist, ien = None, None
932 ist, ien=self.get_mask_indices(mask)
933 if ist is not None and ien is not None:
934 for i in xrange(len(ist)):
935 range=[data[ist[i]],data[ien[i]]]
936 range.sort()
937 masklist.append([range[0],range[1]])
938 return masklist
939
940 def get_mask_indices(self, mask=None):
941 """
942 Compute and Return lists of mask start indices and mask end indices.
943 Parameters:
944 mask: channel mask, created with create_mask.
945 Returns:
946 List of mask start indices and that of mask end indices,
947 i.e., [istart1,istart2,....], [iend1,iend2,....].
948 """
949 if not (isinstance(mask,list) or isinstance(mask, tuple)):
950 raise TypeError("The mask should be list or tuple.")
951 if len(mask) < 2:
952 raise TypeError("The mask elements should be > 1")
953 istart=[]
954 iend=[]
955 if mask[0]: istart.append(0)
956 for i in range(len(mask)-1):
957 if not mask[i] and mask[i+1]:
958 istart.append(i+1)
959 elif mask[i] and not mask[i+1]:
960 iend.append(i)
961 if mask[len(mask)-1]: iend.append(len(mask)-1)
962 if len(istart) != len(iend):
963 raise RuntimeError("Numbers of mask start != mask end.")
964 for i in range(len(istart)):
965 if istart[i] > iend[i]:
966 raise RuntimeError("Mask start index > mask end index")
967 break
968 return istart,iend
969
970# def get_restfreqs(self):
971# """
972# Get the restfrequency(s) stored in this scantable.
973# The return value(s) are always of unit 'Hz'
974# Parameters:
975# none
976# Returns:
977# a list of doubles
978# """
979# return list(self._getrestfreqs())
980
981 def get_restfreqs(self, ids=None):
982 """
983 Get the restfrequency(s) stored in this scantable.
984 The return value(s) are always of unit 'Hz'
985 Parameters:
986 ids: (optional) a list of MOLECULE_ID for that restfrequency(s) to
987 be retrieved
988 Returns:
989 dictionary containing ids and a list of doubles for each id
990 """
991 if ids is None:
992 rfreqs={}
993 idlist = self.getmolnos()
994 for i in idlist:
995 rfreqs[i]=list(self._getrestfreqs(i))
996 return rfreqs
997 else:
998 if type(ids)==list or type(ids)==tuple:
999 rfreqs={}
1000 for i in ids:
1001 rfreqs[i]=list(self._getrestfreqs(i))
1002 return rfreqs
1003 else:
1004 return list(self._getrestfreqs(ids))
1005 #return list(self._getrestfreqs(ids))
1006
1007 def set_restfreqs(self, freqs=None, unit='Hz'):
1008 """
1009 ********NEED TO BE UPDATED begin************
1010 Set or replace the restfrequency specified and
1011 If the 'freqs' argument holds a scalar,
1012 then that rest frequency will be applied to all the selected
1013 data. If the 'freqs' argument holds
1014 a vector, then it MUST be of equal or smaller length than
1015 the number of IFs (and the available restfrequencies will be
1016 replaced by this vector). In this case, *all* data have
1017 the restfrequency set per IF according
1018 to the corresponding value you give in the 'freqs' vector.
1019 E.g. 'freqs=[1e9, 2e9]' would mean IF 0 gets restfreq 1e9 and
1020 IF 1 gets restfreq 2e9.
1021 ********NEED TO BE UPDATED end************
1022 You can also specify the frequencies via a linecatalog.
1023
1024 Parameters:
1025 freqs: list of rest frequency values or string idenitfiers
1026 unit: unit for rest frequency (default 'Hz')
1027
1028 Example:
1029 # set the given restfrequency for the all currently selected IFs
1030 scan.set_restfreqs(freqs=1.4e9)
1031 # set multiple restfrequencies to all the selected data
1032 scan.set_restfreqs(freqs=[1.4e9, 1.41e9, 1.42e9])
1033 # If the number of IFs in the data is >= 2 the IF0 gets the first
1034 # value IF1 the second... NOTE that freqs needs to be
1035 # specified in list of list (e.g. [[],[],...] ).
1036 scan.set_restfreqs(freqs=[[1.4e9],[1.67e9]])
1037 #set the given restfrequency for the whole table (by name)
1038 scan.set_restfreqs(freqs="OH1667")
1039
1040 Note:
1041 To do more sophisticate Restfrequency setting, e.g. on a
1042 source and IF basis, use scantable.set_selection() before using
1043 this function.
1044 # provide your scantable is call scan
1045 selection = selector()
1046 selection.set_name("ORION*")
1047 selection.set_ifs([1])
1048 scan.set_selection(selection)
1049 scan.set_restfreqs(freqs=86.6e9)
1050
1051 """
1052 varlist = vars()
1053 from asap import linecatalog
1054 # simple value
1055 if isinstance(freqs, int) or isinstance(freqs, float):
1056 # TT mod
1057 #self._setrestfreqs(freqs, "",unit)
1058 self._setrestfreqs([freqs], [""],unit)
1059 # list of values
1060 elif isinstance(freqs, list) or isinstance(freqs, tuple):
1061 # list values are scalars
1062 if isinstance(freqs[-1], int) or isinstance(freqs[-1], float):
1063 self._setrestfreqs(freqs, [""],unit)
1064 # list values are tuples, (value, name)
1065 elif isinstance(freqs[-1], dict):
1066 #sel = selector()
1067 #savesel = self._getselection()
1068 #iflist = self.getifnos()
1069 #for i in xrange(len(freqs)):
1070 # sel.set_ifs(iflist[i])
1071 # self._setselection(sel)
1072 # self._setrestfreqs(freqs[i], "",unit)
1073 #self._setselection(savesel)
1074 self._setrestfreqs(freqs["value"],
1075 freqs["name"], "MHz")
1076 elif isinstance(freqs[-1], list) or isinstance(freqs[-1], tuple):
1077 sel = selector()
1078 savesel = self._getselection()
1079 iflist = self.getifnos()
1080 if len(freqs)>len(iflist):
1081 raise ValueError("number of elements in list of list exeeds the current IF selections")
1082 for i in xrange(len(freqs)):
1083 sel.set_ifs(iflist[i])
1084 self._setselection(sel)
1085 self._setrestfreqs(freqs[i]["value"],
1086 freqs[i]["name"], "MHz")
1087 self._setselection(savesel)
1088 # freqs are to be taken from a linecatalog
1089 elif isinstance(freqs, linecatalog):
1090 sel = selector()
1091 savesel = self._getselection()
1092 for i in xrange(freqs.nrow()):
1093 sel.set_ifs(iflist[i])
1094 self._setselection(sel)
1095 self._setrestfreqs(freqs.get_frequency(i),
1096 freqs.get_name(i), "MHz")
1097 # ensure that we are not iterating past nIF
1098 if i == self.nif()-1: break
1099 self._setselection(savesel)
1100 else:
1101 return
1102 self._add_history("set_restfreqs", varlist)
1103
1104 def shift_refpix(self, delta):
1105 """
1106 Shift the reference pixel of the Spectra Coordinate by an
1107 integer amount.
1108 Parameters:
1109 delta: the amount to shift by
1110 Note:
1111 Be careful using this with broadband data.
1112 """
1113 Scantable.shift(self, delta)
1114
1115 def history(self, filename=None):
1116 """
1117 Print the history. Optionally to a file.
1118 Parameters:
1119 filename: The name of the file to save the history to.
1120 """
1121 hist = list(self._gethistory())
1122 out = "-"*80
1123 for h in hist:
1124 if h.startswith("---"):
1125 out += "\n"+h
1126 else:
1127 items = h.split("##")
1128 date = items[0]
1129 func = items[1]
1130 items = items[2:]
1131 out += "\n"+date+"\n"
1132 out += "Function: %s\n Parameters:" % (func)
1133 for i in items:
1134 s = i.split("=")
1135 out += "\n %s = %s" % (s[0], s[1])
1136 out += "\n"+"-"*80
1137 if filename is not None:
1138 if filename is "":
1139 filename = 'scantable_history.txt'
1140 import os
1141 filename = os.path.expandvars(os.path.expanduser(filename))
1142 if not os.path.isdir(filename):
1143 data = open(filename, 'w')
1144 data.write(out)
1145 data.close()
1146 else:
1147 msg = "Illegal file name '%s'." % (filename)
1148 if rcParams['verbose']:
1149 #print msg
1150 casalog.post( msg, 'WARN' )
1151 else:
1152 raise IOError(msg)
1153 if rcParams['verbose']:
1154 try:
1155 from IPython.genutils import page as pager
1156 except ImportError:
1157 from pydoc import pager
1158 pager(out)
1159 else:
1160 return out
1161 return
1162 #
1163 # Maths business
1164 #
1165
1166 def average_time(self, mask=None, scanav=False, weight='tint', align=False):
1167 """
1168 Return the (time) weighted average of a scan.
1169 Note:
1170 in channels only - align if necessary
1171 Parameters:
1172 mask: an optional mask (only used for 'var' and 'tsys'
1173 weighting)
1174 scanav: True averages each scan separately
1175 False (default) averages all scans together,
1176 weight: Weighting scheme.
1177 'none' (mean no weight)
1178 'var' (1/var(spec) weighted)
1179 'tsys' (1/Tsys**2 weighted)
1180 'tint' (integration time weighted)
1181 'tintsys' (Tint/Tsys**2)
1182 'median' ( median averaging)
1183 The default is 'tint'
1184 align: align the spectra in velocity before averaging. It takes
1185 the time of the first spectrum as reference time.
1186 Example:
1187 # time average the scantable without using a mask
1188 newscan = scan.average_time()
1189 """
1190 varlist = vars()
1191 if weight is None: weight = 'TINT'
1192 if mask is None: mask = ()
1193 if scanav: scanav = "SCAN"
1194 else: scanav = "NONE"
1195 scan = (self, )
1196 try:
1197 if align:
1198 scan = (self.freq_align(insitu=False), )
1199 s = None
1200 if weight.upper() == 'MEDIAN':
1201 s = scantable(self._math._averagechannel(scan[0], 'MEDIAN',
1202 scanav))
1203 else:
1204 s = scantable(self._math._average(scan, mask, weight.upper(),
1205 scanav))
1206 except RuntimeError, msg:
1207 if rcParams['verbose']:
1208 #print msg
1209 casalog.post( msg, 'WARN' )
1210 return
1211 else: raise
1212 s._add_history("average_time", varlist)
1213 print_log()
1214 return s
1215
1216 def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None):
1217 """
1218 Return a scan where all spectra are converted to either
1219 Jansky or Kelvin depending upon the flux units of the scan table.
1220 By default the function tries to look the values up internally.
1221 If it can't find them (or if you want to over-ride), you must
1222 specify EITHER jyperk OR eta (and D which it will try to look up
1223 also if you don't set it). jyperk takes precedence if you set both.
1224 Parameters:
1225 jyperk: the Jy / K conversion factor
1226 eta: the aperture efficiency
1227 d: the geomtric diameter (metres)
1228 insitu: if False a new scantable is returned.
1229 Otherwise, the scaling is done in-situ
1230 The default is taken from .asaprc (False)
1231 """
1232 if insitu is None: insitu = rcParams['insitu']
1233 self._math._setinsitu(insitu)
1234 varlist = vars()
1235 if jyperk is None: jyperk = -1.0
1236 if d is None: d = -1.0
1237 if eta is None: eta = -1.0
1238 s = scantable(self._math._convertflux(self, d, eta, jyperk))
1239 s._add_history("convert_flux", varlist)
1240 print_log()
1241 if insitu: self._assign(s)
1242 else: return s
1243
1244 def gain_el(self, poly=None, filename="", method="linear", insitu=None):
1245 """
1246 Return a scan after applying a gain-elevation correction.
1247 The correction can be made via either a polynomial or a
1248 table-based interpolation (and extrapolation if necessary).
1249 You specify polynomial coefficients, an ascii table or neither.
1250 If you specify neither, then a polynomial correction will be made
1251 with built in coefficients known for certain telescopes (an error
1252 will occur if the instrument is not known).
1253 The data and Tsys are *divided* by the scaling factors.
1254 Parameters:
1255 poly: Polynomial coefficients (default None) to compute a
1256 gain-elevation correction as a function of
1257 elevation (in degrees).
1258 filename: The name of an ascii file holding correction factors.
1259 The first row of the ascii file must give the column
1260 names and these MUST include columns
1261 "ELEVATION" (degrees) and "FACTOR" (multiply data
1262 by this) somewhere.
1263 The second row must give the data type of the
1264 column. Use 'R' for Real and 'I' for Integer.
1265 An example file would be
1266 (actual factors are arbitrary) :
1267
1268 TIME ELEVATION FACTOR
1269 R R R
1270 0.1 0 0.8
1271 0.2 20 0.85
1272 0.3 40 0.9
1273 0.4 60 0.85
1274 0.5 80 0.8
1275 0.6 90 0.75
1276 method: Interpolation method when correcting from a table.
1277 Values are "nearest", "linear" (default), "cubic"
1278 and "spline"
1279 insitu: if False a new scantable is returned.
1280 Otherwise, the scaling is done in-situ
1281 The default is taken from .asaprc (False)
1282 """
1283
1284 if insitu is None: insitu = rcParams['insitu']
1285 self._math._setinsitu(insitu)
1286 varlist = vars()
1287 if poly is None:
1288 poly = ()
1289 from os.path import expandvars
1290 filename = expandvars(filename)
1291 s = scantable(self._math._gainel(self, poly, filename, method))
1292 s._add_history("gain_el", varlist)
1293 print_log()
1294 if insitu: self._assign(s)
1295 else: return s
1296
1297 def freq_align(self, reftime=None, method='cubic', insitu=None):
1298 """
1299 Return a scan where all rows have been aligned in frequency/velocity.
1300 The alignment frequency frame (e.g. LSRK) is that set by function
1301 set_freqframe.
1302 Parameters:
1303 reftime: reference time to align at. By default, the time of
1304 the first row of data is used.
1305 method: Interpolation method for regridding the spectra.
1306 Choose from "nearest", "linear", "cubic" (default)
1307 and "spline"
1308 insitu: if False a new scantable is returned.
1309 Otherwise, the scaling is done in-situ
1310 The default is taken from .asaprc (False)
1311 """
1312 if insitu is None: insitu = rcParams["insitu"]
1313 self._math._setinsitu(insitu)
1314 varlist = vars()
1315 if reftime is None: reftime = ""
1316 s = scantable(self._math._freq_align(self, reftime, method))
1317 s._add_history("freq_align", varlist)
1318 print_log()
1319 if insitu: self._assign(s)
1320 else: return s
1321
1322 def opacity(self, tau, insitu=None):
1323 """
1324 Apply an opacity correction. The data
1325 and Tsys are multiplied by the correction factor.
1326 Parameters:
1327 tau: Opacity from which the correction factor is
1328 exp(tau*ZD)
1329 where ZD is the zenith-distance
1330 insitu: if False a new scantable is returned.
1331 Otherwise, the scaling is done in-situ
1332 The default is taken from .asaprc (False)
1333 """
1334 if insitu is None: insitu = rcParams['insitu']
1335 self._math._setinsitu(insitu)
1336 varlist = vars()
1337 s = scantable(self._math._opacity(self, tau))
1338 s._add_history("opacity", varlist)
1339 print_log()
1340 if insitu: self._assign(s)
1341 else: return s
1342
1343 def bin(self, width=5, insitu=None):
1344 """
1345 Return a scan where all spectra have been binned up.
1346 Parameters:
1347 width: The bin width (default=5) in pixels
1348 insitu: if False a new scantable is returned.
1349 Otherwise, the scaling is done in-situ
1350 The default is taken from .asaprc (False)
1351 """
1352 if insitu is None: insitu = rcParams['insitu']
1353 self._math._setinsitu(insitu)
1354 varlist = vars()
1355 s = scantable(self._math._bin(self, width))
1356 s._add_history("bin", varlist)
1357 print_log()
1358 if insitu: self._assign(s)
1359 else: return s
1360
1361
1362 def resample(self, width=5, method='cubic', insitu=None):
1363 """
1364 Return a scan where all spectra have been binned up.
1365
1366 Parameters:
1367 width: The bin width (default=5) in pixels
1368 method: Interpolation method when correcting from a table.
1369 Values are "nearest", "linear", "cubic" (default)
1370 and "spline"
1371 insitu: if False a new scantable is returned.
1372 Otherwise, the scaling is done in-situ
1373 The default is taken from .asaprc (False)
1374 """
1375 if insitu is None: insitu = rcParams['insitu']
1376 self._math._setinsitu(insitu)
1377 varlist = vars()
1378 s = scantable(self._math._resample(self, method, width))
1379 s._add_history("resample", varlist)
1380 print_log()
1381 if insitu: self._assign(s)
1382 else: return s
1383
1384
1385 def average_pol(self, mask=None, weight='none'):
1386 """
1387 Average the Polarisations together.
1388 Parameters:
1389 mask: An optional mask defining the region, where the
1390 averaging will be applied. The output will have all
1391 specified points masked.
1392 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
1393 weighted), or 'tsys' (1/Tsys**2 weighted)
1394 """
1395 varlist = vars()
1396 if mask is None:
1397 mask = ()
1398 s = scantable(self._math._averagepol(self, mask, weight.upper()))
1399 s._add_history("average_pol", varlist)
1400 print_log()
1401 return s
1402
1403 def average_beam(self, mask=None, weight='none'):
1404 """
1405 Average the Beams together.
1406 Parameters:
1407 mask: An optional mask defining the region, where the
1408 averaging will be applied. The output will have all
1409 specified points masked.
1410 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
1411 weighted), or 'tsys' (1/Tsys**2 weighted)
1412 """
1413 varlist = vars()
1414 if mask is None:
1415 mask = ()
1416 s = scantable(self._math._averagebeams(self, mask, weight.upper()))
1417 s._add_history("average_beam", varlist)
1418 print_log()
1419 return s
1420
1421 def convert_pol(self, poltype=None):
1422 """
1423 Convert the data to a different polarisation type.
1424 Parameters:
1425 poltype: The new polarisation type. Valid types are:
1426 "linear", "stokes" and "circular"
1427 """
1428 varlist = vars()
1429 try:
1430 s = scantable(self._math._convertpol(self, poltype))
1431 except RuntimeError, msg:
1432 if rcParams['verbose']:
1433 #print msg
1434 casalog.post( msg, 'WARN' )
1435 return
1436 else:
1437 raise
1438 s._add_history("convert_pol", varlist)
1439 print_log()
1440 return s
1441
1442 def smooth(self, kernel="hanning", width=5.0, insitu=None):
1443 """
1444 Smooth the spectrum by the specified kernel (conserving flux).
1445 Parameters:
1446 kernel: The type of smoothing kernel. Select from
1447 'hanning' (default), 'gaussian', 'boxcar' and
1448 'rmedian'
1449 width: The width of the kernel in pixels. For hanning this is
1450 ignored otherwise it defauls to 5 pixels.
1451 For 'gaussian' it is the Full Width Half
1452 Maximum. For 'boxcar' it is the full width.
1453 For 'rmedian' it is the half width.
1454 insitu: if False a new scantable is returned.
1455 Otherwise, the scaling is done in-situ
1456 The default is taken from .asaprc (False)
1457 Example:
1458 none
1459 """
1460 if insitu is None: insitu = rcParams['insitu']
1461 self._math._setinsitu(insitu)
1462 varlist = vars()
1463 s = scantable(self._math._smooth(self, kernel.lower(), width))
1464 s._add_history("smooth", varlist)
1465 print_log()
1466 if insitu: self._assign(s)
1467 else: return s
1468
1469
1470 def poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None):
1471 """
1472 Return a scan which has been baselined (all rows) by a polynomial.
1473 Parameters:
1474 mask: an optional mask
1475 order: the order of the polynomial (default is 0)
1476 plot: plot the fit and the residual. In this each
1477 indivual fit has to be approved, by typing 'y'
1478 or 'n'
1479 uselin: use linear polynomial fit
1480 insitu: if False a new scantable is returned.
1481 Otherwise, the scaling is done in-situ
1482 The default is taken from .asaprc (False)
1483 Example:
1484 # return a scan baselined by a third order polynomial,
1485 # not using a mask
1486 bscan = scan.poly_baseline(order=3)
1487 """
1488 if insitu is None: insitu = rcParams['insitu']
1489 varlist = vars()
1490 if mask is None:
1491 mask = [True for i in xrange(self.nchan(-1))]
1492 from asap.asapfitter import fitter
1493 try:
1494 f = fitter()
1495 f.set_scan(self, mask)
1496 if uselin:
1497 f.set_function(lpoly=order)
1498 else:
1499 f.set_function(poly=order)
1500 s = f.auto_fit(insitu, plot=plot)
1501 # Save parameters of baseline fits as a class attribute.
1502 # NOTICE: It does not reflect changes in scantable!
1503 self.blpars = f.blpars
1504 s._add_history("poly_baseline", varlist)
1505 print_log()
1506 if insitu: self._assign(s)
1507 else: return s
1508 except RuntimeError:
1509 msg = "The fit failed, possibly because it didn't converge."
1510 if rcParams['verbose']:
1511 #print msg
1512 casalog.post( msg, 'WARN' )
1513 return
1514 else:
1515 raise RuntimeError(msg)
1516
1517
1518 def auto_poly_baseline(self, mask=[], edge=(0, 0), order=0,
1519 threshold=3, chan_avg_limit=1, plot=False,
1520 insitu=None):
1521 """
1522 Return a scan which has been baselined (all rows) by a polynomial.
1523 Spectral lines are detected first using linefinder and masked out
1524 to avoid them affecting the baseline solution.
1525
1526 Parameters:
1527 mask: an optional mask retreived from scantable
1528 edge: an optional number of channel to drop at
1529 the edge of spectrum. If only one value is
1530 specified, the same number will be dropped from
1531 both sides of the spectrum. Default is to keep
1532 all channels. Nested tuples represent individual
1533 edge selection for different IFs (a number of spectral
1534 channels can be different)
1535 order: the order of the polynomial (default is 0)
1536 threshold: the threshold used by line finder. It is better to
1537 keep it large as only strong lines affect the
1538 baseline solution.
1539 chan_avg_limit:
1540 a maximum number of consequtive spectral channels to
1541 average during the search of weak and broad lines.
1542 The default is no averaging (and no search for weak
1543 lines). If such lines can affect the fitted baseline
1544 (e.g. a high order polynomial is fitted), increase this
1545 parameter (usually values up to 8 are reasonable). Most
1546 users of this method should find the default value
1547 sufficient.
1548 plot: plot the fit and the residual. In this each
1549 indivual fit has to be approved, by typing 'y'
1550 or 'n'
1551 insitu: if False a new scantable is returned.
1552 Otherwise, the scaling is done in-situ
1553 The default is taken from .asaprc (False)
1554
1555 Example:
1556 scan2=scan.auto_poly_baseline(order=7)
1557 """
1558 if insitu is None: insitu = rcParams['insitu']
1559 varlist = vars()
1560 from asap.asapfitter import fitter
1561 from asap.asaplinefind import linefinder
1562 from asap import _is_sequence_or_number as _is_valid
1563
1564 # check whether edge is set up for each IF individually
1565 individualedge = False;
1566 if len(edge) > 1:
1567 if isinstance(edge[0], list) or isinstance(edge[0], tuple):
1568 individualedge = True;
1569
1570 if not _is_valid(edge, int) and not individualedge:
1571 raise ValueError, "Parameter 'edge' has to be an integer or a \
1572 pair of integers specified as a tuple. Nested tuples are allowed \
1573 to make individual selection for different IFs."
1574
1575 curedge = (0, 0)
1576 if individualedge:
1577 for edgepar in edge:
1578 if not _is_valid(edgepar, int):
1579 raise ValueError, "Each element of the 'edge' tuple has \
1580 to be a pair of integers or an integer."
1581 else:
1582 curedge = edge;
1583
1584 # setup fitter
1585 f = fitter()
1586 f.set_function(poly=order)
1587
1588 # setup line finder
1589 fl = linefinder()
1590 fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
1591
1592 if not insitu:
1593 workscan = self.copy()
1594 else:
1595 workscan = self
1596
1597 fl.set_scan(workscan)
1598
1599 rows = range(workscan.nrow())
1600 # Save parameters of baseline fits & masklists as a class attribute.
1601 # NOTICE: It does not reflect changes in scantable!
1602 if len(rows) > 0:
1603 self.blpars=[]
1604 self.masklists=[]
1605 asaplog.push("Processing:")
1606 for r in rows:
1607 msg = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % \
1608 (workscan.getscan(r), workscan.getbeam(r), workscan.getif(r), \
1609 workscan.getpol(r), workscan.getcycle(r))
1610 asaplog.push(msg, False)
1611
1612 # figure out edge parameter
1613 if individualedge:
1614 if len(edge) >= workscan.getif(r):
1615 raise RuntimeError, "Number of edge elements appear to " \
1616 "be less than the number of IFs"
1617 curedge = edge[workscan.getif(r)]
1618
1619 # setup line finder
1620 fl.find_lines(r, mask, curedge)
1621 outmask=fl.get_mask()
1622 f.set_scan(workscan, fl.get_mask())
1623 f.x = workscan._getabcissa(r)
1624 f.y = workscan._getspectrum(r)
1625 f.data = None
1626 f.fit()
1627
1628 # Show mask list
1629 masklist=workscan.get_masklist(fl.get_mask(),row=r)
1630 msg = "mask range: "+str(masklist)
1631 asaplog.push(msg, False)
1632
1633 fpar = f.get_parameters()
1634 if plot:
1635 f.plot(residual=True)
1636 x = raw_input("Accept fit ( [y]/n ): ")
1637 if x.upper() == 'N':
1638 self.blpars.append(None)
1639 self.masklists.append(None)
1640 continue
1641 workscan._setspectrum(f.fitter.getresidual(), r)
1642 self.blpars.append(fpar)
1643 self.masklists.append(masklist)
1644 if plot:
1645 f._p.unmap()
1646 f._p = None
1647 workscan._add_history("auto_poly_baseline", varlist)
1648 if insitu:
1649 self._assign(workscan)
1650 else:
1651 return workscan
1652
1653 def rotate_linpolphase(self, angle):
1654 """
1655 Rotate the phase of the complex polarization O=Q+iU correlation.
1656 This is always done in situ in the raw data. So if you call this
1657 function more than once then each call rotates the phase further.
1658 Parameters:
1659 angle: The angle (degrees) to rotate (add) by.
1660 Examples:
1661 scan.rotate_linpolphase(2.3)
1662 """
1663 varlist = vars()
1664 self._math._rotate_linpolphase(self, angle)
1665 self._add_history("rotate_linpolphase", varlist)
1666 print_log()
1667 return
1668
1669
1670 def rotate_xyphase(self, angle):
1671 """
1672 Rotate the phase of the XY correlation. This is always done in situ
1673 in the data. So if you call this function more than once
1674 then each call rotates the phase further.
1675 Parameters:
1676 angle: The angle (degrees) to rotate (add) by.
1677 Examples:
1678 scan.rotate_xyphase(2.3)
1679 """
1680 varlist = vars()
1681 self._math._rotate_xyphase(self, angle)
1682 self._add_history("rotate_xyphase", varlist)
1683 print_log()
1684 return
1685
1686 def swap_linears(self):
1687 """
1688 Swap the linear polarisations XX and YY, or better the first two
1689 polarisations as this also works for ciculars.
1690 """
1691 varlist = vars()
1692 self._math._swap_linears(self)
1693 self._add_history("swap_linears", varlist)
1694 print_log()
1695 return
1696
1697 def invert_phase(self):
1698 """
1699 Invert the phase of the complex polarisation
1700 """
1701 varlist = vars()
1702 self._math._invert_phase(self)
1703 self._add_history("invert_phase", varlist)
1704 print_log()
1705 return
1706
1707 def add(self, offset, insitu=None):
1708 """
1709 Return a scan where all spectra have the offset added
1710 Parameters:
1711 offset: the offset
1712 insitu: if False a new scantable is returned.
1713 Otherwise, the scaling is done in-situ
1714 The default is taken from .asaprc (False)
1715 """
1716 if insitu is None: insitu = rcParams['insitu']
1717 self._math._setinsitu(insitu)
1718 varlist = vars()
1719 s = scantable(self._math._unaryop(self, offset, "ADD", False))
1720 s._add_history("add", varlist)
1721 print_log()
1722 if insitu:
1723 self._assign(s)
1724 else:
1725 return s
1726
1727 def scale(self, factor, tsys=True, insitu=None):
1728 """
1729 Return a scan where all spectra are scaled by the give 'factor'
1730 Parameters:
1731 factor: the scaling factor
1732 insitu: if False a new scantable is returned.
1733 Otherwise, the scaling is done in-situ
1734 The default is taken from .asaprc (False)
1735 tsys: if True (default) then apply the operation to Tsys
1736 as well as the data
1737 """
1738 if insitu is None: insitu = rcParams['insitu']
1739 self._math._setinsitu(insitu)
1740 varlist = vars()
1741 s = scantable(self._math._unaryop(self, factor, "MUL", tsys))
1742 s._add_history("scale", varlist)
1743 print_log()
1744 if insitu:
1745 self._assign(s)
1746 else:
1747 return s
1748
1749 def set_sourcetype(self, match, matchtype="pattern",
1750 sourcetype="reference"):
1751 """
1752 Set the type of the source to be an source or reference scan
1753 using the provided pattern:
1754 Parameters:
1755 match: a Unix style pattern, regular expression or selector
1756 matchtype: 'pattern' (default) UNIX style pattern or
1757 'regex' regular expression
1758 sourcetype: the type of the source to use (source/reference)
1759 """
1760 varlist = vars()
1761 basesel = self.get_selection()
1762 stype = -1
1763 if sourcetype.lower().startswith("r"):
1764 stype = 1
1765 elif sourcetype.lower().startswith("s"):
1766 stype = 0
1767 else:
1768 raise ValueError("Illegal sourcetype use s(ource) or r(eference)")
1769 if matchtype.lower().startswith("p"):
1770 matchtype = "pattern"
1771 elif matchtype.lower().startswith("r"):
1772 matchtype = "regex"
1773 else:
1774 raise ValueError("Illegal matchtype, use p(attern) or r(egex)")
1775 sel = selector()
1776 if isinstance(match, selector):
1777 sel = match
1778 else:
1779 sel.set_query("SRCNAME == %s('%s')" % (matchtype, match))
1780 self.set_selection(basesel+sel)
1781 self._setsourcetype(stype)
1782 self.set_selection(basesel)
1783 s._add_history("set_sourcetype", varlist)
1784
1785 def auto_quotient(self, preserve=True, mode='paired'):
1786 """
1787 This function allows to build quotients automatically.
1788 It assumes the observation to have the same numer of
1789 "ons" and "offs"
1790 Parameters:
1791 preserve: you can preserve (default) the continuum or
1792 remove it. The equations used are
1793 preserve: Output = Toff * (on/off) - Toff
1794 remove: Output = Toff * (on/off) - Ton
1795 mode: the on/off detection mode
1796 'paired' (default)
1797 identifies 'off' scans by the
1798 trailing '_R' (Mopra/Parkes) or
1799 '_e'/'_w' (Tid) and matches
1800 on/off pairs from the observing pattern
1801 'time'
1802 finds the closest off in time
1803
1804 """
1805 modes = ["time", "paired"]
1806 if not mode in modes:
1807 msg = "please provide valid mode. Valid modes are %s" % (modes)
1808 raise ValueError(msg)
1809 varlist = vars()
1810 s = None
1811 if mode.lower() == "paired":
1812 basesel = self.get_selection()
1813 sel = selector()+basesel
1814 sel.set_query("SRCTYPE==1")
1815 self.set_selection(sel)
1816 offs = self.copy()
1817 sel.set_query("SRCTYPE==0")
1818 self.set_selection(sel)
1819 ons = self.copy()
1820 s = scantable(self._math._quotient(ons, offs, preserve))
1821 self.set_selection(basesel)
1822 elif mode.lower() == "time":
1823 s = scantable(self._math._auto_quotient(self, mode, preserve))
1824 s._add_history("auto_quotient", varlist)
1825 print_log()
1826 return s
1827
1828 def mx_quotient(self, mask = None, weight='median', preserve=True):
1829 """
1830 Form a quotient using "off" beams when observing in "MX" mode.
1831 Parameters:
1832 mask: an optional mask to be used when weight == 'stddev'
1833 weight: How to average the off beams. Default is 'median'.
1834 preserve: you can preserve (default) the continuum or
1835 remove it. The equations used are
1836 preserve: Output = Toff * (on/off) - Toff
1837 remove: Output = Toff * (on/off) - Ton
1838 """
1839 if mask is None: mask = ()
1840 varlist = vars()
1841 on = scantable(self._math._mx_extract(self, 'on'))
1842 preoff = scantable(self._math._mx_extract(self, 'off'))
1843 off = preoff.average_time(mask=mask, weight=weight, scanav=False)
1844 from asapmath import quotient
1845 q = quotient(on, off, preserve)
1846 q._add_history("mx_quotient", varlist)
1847 print_log()
1848 return q
1849
1850 def freq_switch(self, insitu=None):
1851 """
1852 Apply frequency switching to the data.
1853 Parameters:
1854 insitu: if False a new scantable is returned.
1855 Otherwise, the swictching is done in-situ
1856 The default is taken from .asaprc (False)
1857 Example:
1858 none
1859 """
1860 if insitu is None: insitu = rcParams['insitu']
1861 self._math._setinsitu(insitu)
1862 varlist = vars()
1863 s = scantable(self._math._freqswitch(self))
1864 s._add_history("freq_switch", varlist)
1865 print_log()
1866 if insitu: self._assign(s)
1867 else: return s
1868
1869 def recalc_azel(self):
1870 """
1871 Recalculate the azimuth and elevation for each position.
1872 Parameters:
1873 none
1874 Example:
1875 """
1876 varlist = vars()
1877 self._recalcazel()
1878 self._add_history("recalc_azel", varlist)
1879 print_log()
1880 return
1881
1882 def __add__(self, other):
1883 varlist = vars()
1884 s = None
1885 if isinstance(other, scantable):
1886 s = scantable(self._math._binaryop(self, other, "ADD"))
1887 elif isinstance(other, float):
1888 s = scantable(self._math._unaryop(self, other, "ADD", False))
1889 else:
1890 raise TypeError("Other input is not a scantable or float value")
1891 s._add_history("operator +", varlist)
1892 print_log()
1893 return s
1894
1895 def __sub__(self, other):
1896 """
1897 implicit on all axes and on Tsys
1898 """
1899 varlist = vars()
1900 s = None
1901 if isinstance(other, scantable):
1902 s = scantable(self._math._binaryop(self, other, "SUB"))
1903 elif isinstance(other, float):
1904 s = scantable(self._math._unaryop(self, other, "SUB", False))
1905 else:
1906 raise TypeError("Other input is not a scantable or float value")
1907 s._add_history("operator -", varlist)
1908 print_log()
1909 return s
1910
1911 def __mul__(self, other):
1912 """
1913 implicit on all axes and on Tsys
1914 """
1915 varlist = vars()
1916 s = None
1917 if isinstance(other, scantable):
1918 s = scantable(self._math._binaryop(self, other, "MUL"))
1919 elif isinstance(other, float):
1920 s = scantable(self._math._unaryop(self, other, "MUL", False))
1921 else:
1922 raise TypeError("Other input is not a scantable or float value")
1923 s._add_history("operator *", varlist)
1924 print_log()
1925 return s
1926
1927
1928 def __div__(self, other):
1929 """
1930 implicit on all axes and on Tsys
1931 """
1932 varlist = vars()
1933 s = None
1934 if isinstance(other, scantable):
1935 s = scantable(self._math._binaryop(self, other, "DIV"))
1936 elif isinstance(other, float):
1937 if other == 0.0:
1938 raise ZeroDivisionError("Dividing by zero is not recommended")
1939 s = scantable(self._math._unaryop(self, other, "DIV", False))
1940 else:
1941 raise TypeError("Other input is not a scantable or float value")
1942 s._add_history("operator /", varlist)
1943 print_log()
1944 return s
1945
1946 def get_fit(self, row=0):
1947 """
1948 Print or return the stored fits for a row in the scantable
1949 Parameters:
1950 row: the row which the fit has been applied to.
1951 """
1952 if row > self.nrow():
1953 return
1954 from asap.asapfit import asapfit
1955 fit = asapfit(self._getfit(row))
1956 if rcParams['verbose']:
1957 #print fit
1958 casalog.post( '%s' %(fit) )
1959 return
1960 else:
1961 return fit.as_dict()
1962
1963 def flag_nans(self):
1964 """
1965 Utility function to flag NaN values in the scantable.
1966 """
1967 import numpy
1968 basesel = self.get_selection()
1969 for i in range(self.nrow()):
1970 sel = selector()+basesel
1971 sel.set_scans(self.getscan(i))
1972 sel.set_beams(self.getbeam(i))
1973 sel.set_ifs(self.getif(i))
1974 sel.set_polarisations(self.getpol(i))
1975 self.set_selection(sel)
1976 nans = numpy.isnan(self._getspectrum(0))
1977 if numpy.any(nans):
1978 bnans = [ bool(v) for v in nans]
1979 self.flag(bnans)
1980 self.set_selection(basesel)
1981
1982
1983 def _add_history(self, funcname, parameters):
1984 if not rcParams['scantable.history']:
1985 return
1986 # create date
1987 sep = "##"
1988 from datetime import datetime
1989 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
1990 hist = dstr+sep
1991 hist += funcname+sep#cdate+sep
1992 if parameters.has_key('self'): del parameters['self']
1993 for k, v in parameters.iteritems():
1994 if type(v) is dict:
1995 for k2, v2 in v.iteritems():
1996 hist += k2
1997 hist += "="
1998 if isinstance(v2, scantable):
1999 hist += 'scantable'
2000 elif k2 == 'mask':
2001 if isinstance(v2, list) or isinstance(v2, tuple):
2002 hist += str(self._zip_mask(v2))
2003 else:
2004 hist += str(v2)
2005 else:
2006 hist += str(v2)
2007 else:
2008 hist += k
2009 hist += "="
2010 if isinstance(v, scantable):
2011 hist += 'scantable'
2012 elif k == 'mask':
2013 if isinstance(v, list) or isinstance(v, tuple):
2014 hist += str(self._zip_mask(v))
2015 else:
2016 hist += str(v)
2017 else:
2018 hist += str(v)
2019 hist += sep
2020 hist = hist[:-2] # remove trailing '##'
2021 self._addhistory(hist)
2022
2023
2024 def _zip_mask(self, mask):
2025 mask = list(mask)
2026 i = 0
2027 segments = []
2028 while mask[i:].count(1):
2029 i += mask[i:].index(1)
2030 if mask[i:].count(0):
2031 j = i + mask[i:].index(0)
2032 else:
2033 j = len(mask)
2034 segments.append([i, j])
2035 i = j
2036 return segments
2037
2038 def _get_ordinate_label(self):
2039 fu = "("+self.get_fluxunit()+")"
2040 import re
2041 lbl = "Intensity"
2042 if re.match(".K.", fu):
2043 lbl = "Brightness Temperature "+ fu
2044 elif re.match(".Jy.", fu):
2045 lbl = "Flux density "+ fu
2046 return lbl
2047
2048 def _check_ifs(self):
2049 nchans = [self.nchan(i) for i in range(self.nif(-1))]
2050 nchans = filter(lambda t: t > 0, nchans)
2051 return (sum(nchans)/len(nchans) == nchans[0])
2052
2053 def _fill(self, names, unit, average, getpt):
2054 import os
2055 from asap._asap import stfiller
2056 first = True
2057 fullnames = []
2058 for name in names:
2059 name = os.path.expandvars(name)
2060 name = os.path.expanduser(name)
2061 if not os.path.exists(name):
2062 msg = "File '%s' does not exists" % (name)
2063 if rcParams['verbose']:
2064 asaplog.push(msg)
2065 #print asaplog.pop().strip()
2066 casalog.post( asaplog.pop().strip(), 'WARN' )
2067 return
2068 raise IOError(msg)
2069 fullnames.append(name)
2070 if average:
2071 asaplog.push('Auto averaging integrations')
2072 stype = int(rcParams['scantable.storage'].lower() == 'disk')
2073 for name in fullnames:
2074 tbl = Scantable(stype)
2075 r = stfiller(tbl)
2076 rx = rcParams['scantable.reference']
2077 r._setreferenceexpr(rx)
2078 msg = "Importing %s..." % (name)
2079 asaplog.push(msg, False)
2080 print_log()
2081 r._open(name, -1, -1, getpt)
2082 r._read()
2083 if average:
2084 tbl = self._math._average((tbl, ), (), 'NONE', 'SCAN')
2085 if not first:
2086 tbl = self._math._merge([self, tbl])
2087 Scantable.__init__(self, tbl)
2088 r._close()
2089 del r, tbl
2090 first = False
2091 if unit is not None:
2092 self.set_fluxunit(unit)
2093 #self.set_freqframe(rcParams['scantable.freqframe'])
2094
2095 def __getitem__(self, key):
2096 if key < 0:
2097 key += self.nrow()
2098 if key >= self.nrow():
2099 raise IndexError("Row index out of range.")
2100 return self._getspectrum(key)
2101
2102 def __setitem__(self, key, value):
2103 if key < 0:
2104 key += self.nrow()
2105 if key >= self.nrow():
2106 raise IndexError("Row index out of range.")
2107 if not hasattr(value, "__len__") or \
2108 len(value) > self.nchan(self.getif(key)):
2109 raise ValueError("Spectrum length doesn't match.")
2110 return self._setspectrum(value, key)
2111
2112 def __len__(self):
2113 return self.nrow()
2114
2115 def __iter__(self):
2116 for i in range(len(self)):
2117 yield self[i]
Note: See TracBrowser for help on using the repository browser.