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

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

New Development: No

JIRA Issue: Yes CAS-729, CAS-1147

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes

Module(s): Module Names change impacts.

Description: Describe your changes here...

I have changed that almost all log messages are output to casapy.log,
not to the terminal window. After this change, asap becomes to depend on casapy
and is not running in standalone, because asap have to import taskinit module
to access casalogger.


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