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

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

New Development: No

JIRA Issue: Yes CAS-1810

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: Added 'antenna' parameter to scantable constructor

Test Programs: List test programs

Put in Release Notes: Yes

Module(s): atnf

Description: Describe your changes here...

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