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

Last change on this file since 1688 was 1685, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-1810

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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