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

Last change on this file since 1655 was 1655, checked in by WataruKawasaki, 15 years ago

New Development: Yes

JIRA Issue: Yes (CAS-1433)

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: A new functionality flag_row() added

Test Programs:

Put in Release Notes: No

Module(s): sdflag

Description: Added a new functionality flag_row() to

enable row_based flagging/unflagging.


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