source: trunk/python/scantable.py@ 2630

Last change on this file since 2630 was 2611, checked in by Kana Sugimoto, 12 years ago

New Development: Yes

JIRA Issue: Yes (CAS-2820/Trac-277)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: a new method, scantable.parse_idx_selection(mode, expr)

Test Programs: comming soon

Put in Release Notes: No

Module(s): scantable and CASA SD tasks

Description:

Added a new method, scantable.parse_idx_selection(mode, expr), in scantable class
to enable CASA type selection syntax of SCANNO, IFNO, POLNO, BEAMNO and row IDs
in scantable.
This method parses a comma separated selection expression (expr) in a string to
a list of IDs to set to selector. Index ranges in scantable are taken into account
when parsing.
Available modes are, 'scan' (parses SCANNO selection),
'if' (prases IFNO selection), 'pol' (parses POLNO selection), 'beam' (parses BEAMNO
selection), and 'row' (parses row id selction).

Usage:
stab = asap.scantable("SCANTABLE_NAME")
stab.scannos()
# assume this returns successive numbers from 2 to 20.
scanlist = stab.parse_idx_selection("scan","<5,10~13,>18")
# returns [2,3,4, 10,11,12,13, 19,20]
stab.set_selection(scans=scanlist)

Also made "<=" and ">=" available in scantable._parse_selection().

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 145.6 KB
Line 
1"""This module defines the scantable class."""
2
3import os
4import tempfile
5import numpy
6try:
7 from functools import wraps as wraps_dec
8except ImportError:
9 from asap.compatibility import wraps as wraps_dec
10
11from asap.env import is_casapy
12from asap._asap import Scantable
13from asap._asap import filler, msfiller
14from asap.parameters import rcParams
15from asap.logging import asaplog, asaplog_post_dec
16from asap.selector import selector
17from asap.linecatalog import linecatalog
18from asap.coordinate import coordinate
19from asap.utils import _n_bools, mask_not, mask_and, mask_or, page
20from asap.asapfitter import fitter
21
22def preserve_selection(func):
23 @wraps_dec(func)
24 def wrap(obj, *args, **kw):
25 basesel = obj.get_selection()
26 try:
27 val = func(obj, *args, **kw)
28 finally:
29 obj.set_selection(basesel)
30 return val
31 return wrap
32
33def is_scantable(filename):
34 """Is the given file a scantable?
35
36 Parameters:
37
38 filename: the name of the file/directory to test
39
40 """
41 if ( os.path.isdir(filename)
42 and os.path.exists(filename+'/table.info')
43 and os.path.exists(filename+'/table.dat') ):
44 f=open(filename+'/table.info')
45 l=f.readline()
46 f.close()
47 #if ( l.find('Scantable') != -1 ):
48 if ( l.find('Measurement Set') == -1 ):
49 return True
50 else:
51 return False
52 else:
53 return False
54## return (os.path.isdir(filename)
55## and not os.path.exists(filename+'/table.f1')
56## and os.path.exists(filename+'/table.info'))
57
58def is_ms(filename):
59 """Is the given file a MeasurementSet?
60
61 Parameters:
62
63 filename: the name of the file/directory to test
64
65 """
66 if ( os.path.isdir(filename)
67 and os.path.exists(filename+'/table.info')
68 and os.path.exists(filename+'/table.dat') ):
69 f=open(filename+'/table.info')
70 l=f.readline()
71 f.close()
72 if ( l.find('Measurement Set') != -1 ):
73 return True
74 else:
75 return False
76 else:
77 return False
78
79def normalise_edge_param(edge):
80 """\
81 Convert a given edge value to a one-dimensional array that can be
82 given to baseline-fitting/subtraction functions.
83 The length of the output value will be an even because values for
84 the both sides of spectra are to be contained for each IF. When
85 the length is 2, the values will be applied to all IFs. If the length
86 is larger than 2, it will be 2*ifnos().
87 Accepted format of edge include:
88 * an integer - will be used for both sides of spectra of all IFs.
89 e.g. 10 is converted to [10,10]
90 * an empty list/tuple [] - converted to [0, 0] and used for all IFs.
91 * a list/tuple containing an integer - same as the above case.
92 e.g. [10] is converted to [10,10]
93 * a list/tuple containing two integers - will be used for all IFs.
94 e.g. [5,10] is output as it is. no need to convert.
95 * a list/tuple of lists/tuples containing TWO integers -
96 each element of edge will be used for each IF.
97 e.g. [[5,10],[15,20]] - [5,10] for IF[0] and [15,20] for IF[1].
98
99 If an element contains the same integer values, the input 'edge'
100 parameter can be given in a simpler shape in the following cases:
101 ** when len(edge)!=2
102 any elements containing the same values can be replaced
103 to single integers.
104 e.g. [[15,15]] can be simplified to [15] (or [15,15] or 15 also).
105 e.g. [[1,1],[2,2],[3,3]] can be simplified to [1,2,3].
106 ** when len(edge)=2
107 care is needed for this case: ONLY ONE of the
108 elements can be a single integer,
109 e.g. [[5,5],[10,10]] can be simplified to [5,[10,10]]
110 or [[5,5],10], but can NOT be simplified to [5,10].
111 when [5,10] given, it is interpreted as
112 [[5,10],[5,10],[5,10],...] instead, as shown before.
113 """
114 from asap import _is_sequence_or_number as _is_valid
115 if isinstance(edge, list) or isinstance(edge, tuple):
116 for edgepar in edge:
117 if not _is_valid(edgepar, int):
118 raise ValueError, "Each element of the 'edge' tuple has \
119 to be a pair of integers or an integer."
120 if isinstance(edgepar, list) or isinstance(edgepar, tuple):
121 if len(edgepar) != 2:
122 raise ValueError, "Each element of the 'edge' tuple has \
123 to be a pair of integers or an integer."
124 else:
125 if not _is_valid(edge, int):
126 raise ValueError, "Parameter 'edge' has to be an integer or a \
127 pair of integers specified as a tuple. \
128 Nested tuples are allowed \
129 to make individual selection for different IFs."
130
131
132 if isinstance(edge, int):
133 edge = [ edge, edge ] # e.g. 3 => [3,3]
134 elif isinstance(edge, list) or isinstance(edge, tuple):
135 if len(edge) == 0:
136 edge = [0, 0] # e.g. [] => [0,0]
137 elif len(edge) == 1:
138 if isinstance(edge[0], int):
139 edge = [ edge[0], edge[0] ] # e.g. [1] => [1,1]
140
141 commonedge = True
142 if len(edge) > 2: commonedge = False
143 else:
144 for edgepar in edge:
145 if isinstance(edgepar, list) or isinstance(edgepar, tuple):
146 commonedge = False
147 break
148
149 if commonedge:
150 if len(edge) > 1:
151 norm_edge = edge
152 else:
153 norm_edge = edge + edge
154 else:
155 norm_edge = []
156 for edgepar in edge:
157 if isinstance(edgepar, int):
158 norm_edge += [edgepar, edgepar]
159 else:
160 norm_edge += edgepar
161
162 return norm_edge
163
164def raise_fitting_failure_exception(e):
165 msg = "The fit failed, possibly because it didn't converge."
166 if rcParams["verbose"]:
167 asaplog.push(str(e))
168 asaplog.push(str(msg))
169 else:
170 raise RuntimeError(str(e)+'\n'+msg)
171
172def pack_progress_params(showprogress, minnrow):
173 return str(showprogress).lower() + ',' + str(minnrow)
174
175class scantable(Scantable):
176 """\
177 The ASAP container for scans (single-dish data).
178 """
179
180 @asaplog_post_dec
181 def __init__(self, filename, average=None, unit=None, parallactify=None,
182 **args):
183 """\
184 Create a scantable from a saved one or make a reference
185
186 Parameters:
187
188 filename: the name of an asap table on disk
189 or
190 the name of a rpfits/sdfits/ms file
191 (integrations within scans are auto averaged
192 and the whole file is read) or
193 [advanced] a reference to an existing scantable
194
195 average: average all integrations withinb a scan on read.
196 The default (True) is taken from .asaprc.
197
198 unit: brightness unit; must be consistent with K or Jy.
199 Over-rides the default selected by the filler
200 (input rpfits/sdfits/ms) or replaces the value
201 in existing scantables
202
203 antenna: for MeasurementSet input data only:
204 Antenna selection. integer (id) or string
205 (name or id).
206
207 parallactify: Indicate that the data had been parallactified.
208 Default (false) is taken from rc file.
209
210 """
211 if average is None:
212 average = rcParams['scantable.autoaverage']
213 parallactify = parallactify or rcParams['scantable.parallactify']
214 varlist = vars()
215 from asap._asap import stmath
216 self._math = stmath( rcParams['insitu'] )
217 if isinstance(filename, Scantable):
218 Scantable.__init__(self, filename)
219 else:
220 if isinstance(filename, str):
221 filename = os.path.expandvars(filename)
222 filename = os.path.expanduser(filename)
223 if not os.path.exists(filename):
224 s = "File '%s' not found." % (filename)
225 raise IOError(s)
226 if is_scantable(filename):
227 ondisk = rcParams['scantable.storage'] == 'disk'
228 Scantable.__init__(self, filename, ondisk)
229 if unit is not None:
230 self.set_fluxunit(unit)
231 if average:
232 self._assign( self.average_time( scanav=True ) )
233 # do not reset to the default freqframe
234 #self.set_freqframe(rcParams['scantable.freqframe'])
235 elif is_ms(filename):
236 # Measurement Set
237 opts={'ms': {}}
238 mskeys=['getpt','antenna']
239 for key in mskeys:
240 if key in args.keys():
241 opts['ms'][key] = args[key]
242 self._fill([filename], unit, average, opts)
243 elif os.path.isfile(filename):
244 self._fill([filename], unit, average)
245 # only apply to new data not "copy constructor"
246 self.parallactify(parallactify)
247 else:
248 msg = "The given file '%s'is not a valid " \
249 "asap table." % (filename)
250 raise IOError(msg)
251 elif (isinstance(filename, list) or isinstance(filename, tuple)) \
252 and isinstance(filename[-1], str):
253 self._fill(filename, unit, average)
254 self.parallactify(parallactify)
255 self._add_history("scantable", varlist)
256
257 @asaplog_post_dec
258 def save(self, name=None, format=None, overwrite=False):
259 """\
260 Store the scantable on disk. This can be an asap (aips++) Table,
261 SDFITS or MS2 format.
262
263 Parameters:
264
265 name: the name of the outputfile. For format 'ASCII'
266 this is the root file name (data in 'name'.txt
267 and header in 'name'_header.txt)
268
269 format: an optional file format. Default is ASAP.
270 Allowed are:
271
272 * 'ASAP' (save as ASAP [aips++] Table),
273 * 'SDFITS' (save as SDFITS file)
274 * 'ASCII' (saves as ascii text file)
275 * 'MS2' (saves as an casacore MeasurementSet V2)
276 * 'FITS' (save as image FITS - not readable by
277 class)
278 * 'CLASS' (save as FITS readable by CLASS)
279
280 overwrite: If the file should be overwritten if it exists.
281 The default False is to return with warning
282 without writing the output. USE WITH CARE.
283
284 Example::
285
286 scan.save('myscan.asap')
287 scan.save('myscan.sdfits', 'SDFITS')
288
289 """
290 from os import path
291 format = format or rcParams['scantable.save']
292 suffix = '.'+format.lower()
293 if name is None or name == "":
294 name = 'scantable'+suffix
295 msg = "No filename given. Using default name %s..." % name
296 asaplog.push(msg)
297 name = path.expandvars(name)
298 if path.isfile(name) or path.isdir(name):
299 if not overwrite:
300 msg = "File %s exists." % name
301 raise IOError(msg)
302 format2 = format.upper()
303 if format2 == 'ASAP':
304 self._save(name)
305 elif format2 == 'MS2':
306 msopt = {'ms': {'overwrite': overwrite } }
307 from asap._asap import mswriter
308 writer = mswriter( self )
309 writer.write( name, msopt )
310 else:
311 from asap._asap import stwriter as stw
312 writer = stw(format2)
313 writer.write(self, name)
314 return
315
316 def copy(self):
317 """Return a copy of this scantable.
318
319 *Note*:
320
321 This makes a full (deep) copy. scan2 = scan1 makes a reference.
322
323 Example::
324
325 copiedscan = scan.copy()
326
327 """
328 sd = scantable(Scantable._copy(self))
329 return sd
330
331 def drop_scan(self, scanid=None):
332 """\
333 Return a new scantable where the specified scan number(s) has(have)
334 been dropped.
335
336 Parameters:
337
338 scanid: a (list of) scan number(s)
339
340 """
341 from asap import _is_sequence_or_number as _is_valid
342 from asap import _to_list
343 from asap import unique
344 if not _is_valid(scanid):
345 raise RuntimeError( 'Please specify a scanno to drop from the'
346 ' scantable' )
347 scanid = _to_list(scanid)
348 allscans = unique([ self.getscan(i) for i in range(self.nrow())])
349 for sid in scanid: allscans.remove(sid)
350 if len(allscans) == 0:
351 raise ValueError("Can't remove all scans")
352 sel = selector(scans=allscans)
353 return self._select_copy(sel)
354
355 def _select_copy(self, selection):
356 orig = self.get_selection()
357 self.set_selection(orig+selection)
358 cp = self.copy()
359 self.set_selection(orig)
360 return cp
361
362 def get_scan(self, scanid=None):
363 """\
364 Return a specific scan (by scanno) or collection of scans (by
365 source name) in a new scantable.
366
367 *Note*:
368
369 See scantable.drop_scan() for the inverse operation.
370
371 Parameters:
372
373 scanid: a (list of) scanno or a source name, unix-style
374 patterns are accepted for source name matching, e.g.
375 '*_R' gets all 'ref scans
376
377 Example::
378
379 # get all scans containing the source '323p459'
380 newscan = scan.get_scan('323p459')
381 # get all 'off' scans
382 refscans = scan.get_scan('*_R')
383 # get a susbset of scans by scanno (as listed in scan.summary())
384 newscan = scan.get_scan([0, 2, 7, 10])
385
386 """
387 if scanid is None:
388 raise RuntimeError( 'Please specify a scan no or name to '
389 'retrieve from the scantable' )
390 try:
391 bsel = self.get_selection()
392 sel = selector()
393 if type(scanid) is str:
394 sel.set_name(scanid)
395 return self._select_copy(sel)
396 elif type(scanid) is int:
397 sel.set_scans([scanid])
398 return self._select_copy(sel)
399 elif type(scanid) is list:
400 sel.set_scans(scanid)
401 return self._select_copy(sel)
402 else:
403 msg = "Illegal scanid type, use 'int' or 'list' if ints."
404 raise TypeError(msg)
405 except RuntimeError:
406 raise
407
408 def __str__(self):
409 tempFile = tempfile.NamedTemporaryFile()
410 Scantable._summary(self, tempFile.name)
411 tempFile.seek(0)
412 asaplog.clear()
413 return tempFile.file.read()
414
415 @asaplog_post_dec
416 def summary(self, filename=None):
417 """\
418 Print a summary of the contents of this scantable.
419
420 Parameters:
421
422 filename: the name of a file to write the putput to
423 Default - no file output
424
425 """
426 if filename is not None:
427 if filename is "":
428 filename = 'scantable_summary.txt'
429 from os.path import expandvars, isdir
430 filename = expandvars(filename)
431 if isdir(filename):
432 msg = "Illegal file name '%s'." % (filename)
433 raise IOError(msg)
434 else:
435 filename = ""
436 Scantable._summary(self, filename)
437
438 def get_spectrum(self, rowno):
439 """Return the spectrum for the current row in the scantable as a list.
440
441 Parameters:
442
443 rowno: the row number to retrieve the spectrum from
444
445 """
446 return self._getspectrum(rowno)
447
448 def get_mask(self, rowno):
449 """Return the mask for the current row in the scantable as a list.
450
451 Parameters:
452
453 rowno: the row number to retrieve the mask from
454
455 """
456 return self._getmask(rowno)
457
458 def set_spectrum(self, spec, rowno):
459 """Set the spectrum for the current row in the scantable.
460
461 Parameters:
462
463 spec: the new spectrum
464
465 rowno: the row number to set the spectrum for
466
467 """
468 assert(len(spec) == self.nchan(self.getif(rowno)))
469 return self._setspectrum(spec, rowno)
470
471 def get_coordinate(self, rowno):
472 """Return the (spectral) coordinate for a a given 'rowno'.
473
474 *Note*:
475
476 * This coordinate is only valid until a scantable method modifies
477 the frequency axis.
478 * This coordinate does contain the original frequency set-up
479 NOT the new frame. The conversions however are done using the user
480 specified frame (e.g. LSRK/TOPO). To get the 'real' coordinate,
481 use scantable.freq_align first. Without it there is no closure,
482 i.e.::
483
484 c = myscan.get_coordinate(0)
485 c.to_frequency(c.get_reference_pixel()) != c.get_reference_value()
486
487 Parameters:
488
489 rowno: the row number for the spectral coordinate
490
491 """
492 return coordinate(Scantable.get_coordinate(self, rowno))
493
494 def get_selection(self):
495 """\
496 Get the selection object currently set on this scantable.
497
498 Example::
499
500 sel = scan.get_selection()
501 sel.set_ifs(0) # select IF 0
502 scan.set_selection(sel) # apply modified selection
503
504 """
505 return selector(self._getselection())
506
507 def set_selection(self, selection=None, **kw):
508 """\
509 Select a subset of the data. All following operations on this scantable
510 are only applied to thi selection.
511
512 Parameters:
513
514 selection: a selector object (default unset the selection), or
515 any combination of 'pols', 'ifs', 'beams', 'scans',
516 'cycles', 'name', 'query'
517
518 Examples::
519
520 sel = selector() # create a selection object
521 self.set_scans([0, 3]) # select SCANNO 0 and 3
522 scan.set_selection(sel) # set the selection
523 scan.summary() # will only print summary of scanno 0 an 3
524 scan.set_selection() # unset the selection
525 # or the equivalent
526 scan.set_selection(scans=[0,3])
527 scan.summary() # will only print summary of scanno 0 an 3
528 scan.set_selection() # unset the selection
529
530 """
531 if selection is None:
532 # reset
533 if len(kw) == 0:
534 selection = selector()
535 else:
536 # try keywords
537 for k in kw:
538 if k not in selector.fields:
539 raise KeyError("Invalid selection key '%s', "
540 "valid keys are %s" % (k,
541 selector.fields))
542 selection = selector(**kw)
543 self._setselection(selection)
544
545 def get_row(self, row=0, insitu=None):
546 """\
547 Select a row in the scantable.
548 Return a scantable with single row.
549
550 Parameters:
551
552 row: row no of integration, default is 0.
553 insitu: if False a new scantable is returned. Otherwise, the
554 scaling is done in-situ. The default is taken from .asaprc
555 (False)
556
557 """
558 if insitu is None:
559 insitu = rcParams['insitu']
560 if not insitu:
561 workscan = self.copy()
562 else:
563 workscan = self
564 # Select a row
565 sel = selector()
566 sel.set_rows([row])
567 workscan.set_selection(sel)
568 if not workscan.nrow() == 1:
569 msg = "Could not identify single row. %d rows selected." \
570 % (workscan.nrow())
571 raise RuntimeError(msg)
572 if insitu:
573 self._assign(workscan)
574 else:
575 return workscan
576
577 @asaplog_post_dec
578 def stats(self, stat='stddev', mask=None, form='3.3f', row=None):
579 """\
580 Determine the specified statistic of the current beam/if/pol
581 Takes a 'mask' as an optional parameter to specify which
582 channels should be excluded.
583
584 Parameters:
585
586 stat: 'min', 'max', 'min_abc', 'max_abc', 'sumsq', 'sum',
587 'mean', 'var', 'stddev', 'avdev', 'rms', 'median'
588
589 mask: an optional mask specifying where the statistic
590 should be determined.
591
592 form: format string to print statistic values
593
594 row: row number of spectrum to process.
595 (default is None: for all rows)
596
597 Example:
598 scan.set_unit('channel')
599 msk = scan.create_mask([100, 200], [500, 600])
600 scan.stats(stat='mean', mask=m)
601
602 """
603 mask = mask or []
604 if not self._check_ifs():
605 raise ValueError("Cannot apply mask as the IFs have different "
606 "number of channels. Please use setselection() "
607 "to select individual IFs")
608 rtnabc = False
609 if stat.lower().endswith('_abc'): rtnabc = True
610 getchan = False
611 if stat.lower().startswith('min') or stat.lower().startswith('max'):
612 chan = self._math._minmaxchan(self, mask, stat)
613 getchan = True
614 statvals = []
615 if not rtnabc:
616 if row == None:
617 statvals = self._math._stats(self, mask, stat)
618 else:
619 statvals = self._math._statsrow(self, mask, stat, int(row))
620
621 #def cb(i):
622 # return statvals[i]
623
624 #return self._row_callback(cb, stat)
625
626 label=stat
627 #callback=cb
628 out = ""
629 #outvec = []
630 sep = '-'*50
631
632 if row == None:
633 rows = xrange(self.nrow())
634 elif isinstance(row, int):
635 rows = [ row ]
636
637 for i in rows:
638 refstr = ''
639 statunit= ''
640 if getchan:
641 qx, qy = self.chan2data(rowno=i, chan=chan[i])
642 if rtnabc:
643 statvals.append(qx['value'])
644 refstr = ('(value: %'+form) % (qy['value'])+' ['+qy['unit']+'])'
645 statunit= '['+qx['unit']+']'
646 else:
647 refstr = ('(@ %'+form) % (qx['value'])+' ['+qx['unit']+'])'
648
649 tm = self._gettime(i)
650 src = self._getsourcename(i)
651 out += 'Scan[%d] (%s) ' % (self.getscan(i), src)
652 out += 'Time[%s]:\n' % (tm)
653 if self.nbeam(-1) > 1: out += ' Beam[%d] ' % (self.getbeam(i))
654 if self.nif(-1) > 1: out += ' IF[%d] ' % (self.getif(i))
655 if self.npol(-1) > 1: out += ' Pol[%d] ' % (self.getpol(i))
656 #outvec.append(callback(i))
657 if len(rows) > 1:
658 # out += ('= %'+form) % (outvec[i]) +' '+refstr+'\n'
659 out += ('= %'+form) % (statvals[i]) +' '+refstr+'\n'
660 else:
661 # out += ('= %'+form) % (outvec[0]) +' '+refstr+'\n'
662 out += ('= %'+form) % (statvals[0]) +' '+refstr+'\n'
663 out += sep+"\n"
664
665 import os
666 if os.environ.has_key( 'USER' ):
667 usr = os.environ['USER']
668 else:
669 import commands
670 usr = commands.getoutput( 'whoami' )
671 tmpfile = '/tmp/tmp_'+usr+'_casapy_asap_scantable_stats'
672 f = open(tmpfile,'w')
673 print >> f, sep
674 print >> f, ' %s %s' % (label, statunit)
675 print >> f, sep
676 print >> f, out
677 f.close()
678 f = open(tmpfile,'r')
679 x = f.readlines()
680 f.close()
681 asaplog.push(''.join(x), False)
682
683 return statvals
684
685 def chan2data(self, rowno=0, chan=0):
686 """\
687 Returns channel/frequency/velocity and spectral value
688 at an arbitrary row and channel in the scantable.
689
690 Parameters:
691
692 rowno: a row number in the scantable. Default is the
693 first row, i.e. rowno=0
694
695 chan: a channel in the scantable. Default is the first
696 channel, i.e. pos=0
697
698 """
699 if isinstance(rowno, int) and isinstance(chan, int):
700 qx = {'unit': self.get_unit(),
701 'value': self._getabcissa(rowno)[chan]}
702 qy = {'unit': self.get_fluxunit(),
703 'value': self._getspectrum(rowno)[chan]}
704 return qx, qy
705
706 def stddev(self, mask=None):
707 """\
708 Determine the standard deviation of the current beam/if/pol
709 Takes a 'mask' as an optional parameter to specify which
710 channels should be excluded.
711
712 Parameters:
713
714 mask: an optional mask specifying where the standard
715 deviation should be determined.
716
717 Example::
718
719 scan.set_unit('channel')
720 msk = scan.create_mask([100, 200], [500, 600])
721 scan.stddev(mask=m)
722
723 """
724 return self.stats(stat='stddev', mask=mask);
725
726
727 def get_column_names(self):
728 """\
729 Return a list of column names, which can be used for selection.
730 """
731 return list(Scantable.get_column_names(self))
732
733 def get_tsys(self, row=-1):
734 """\
735 Return the System temperatures.
736
737 Parameters:
738
739 row: the rowno to get the information for. (default all rows)
740
741 Returns:
742
743 a list of Tsys values for the current selection
744
745 """
746 if row > -1:
747 return self._get_column(self._gettsys, row)
748 return self._row_callback(self._gettsys, "Tsys")
749
750 def get_tsysspectrum(self, row=-1):
751 """\
752 Return the channel dependent system temperatures.
753
754 Parameters:
755
756 row: the rowno to get the information for. (default all rows)
757
758 Returns:
759
760 a list of Tsys values for the current selection
761
762 """
763 return self._get_column( self._gettsysspectrum, row )
764
765 def get_weather(self, row=-1):
766 """\
767 Return the weather informations.
768
769 Parameters:
770
771 row: the rowno to get the information for. (default all rows)
772
773 Returns:
774
775 a dict or list of of dicts of values for the current selection
776
777 """
778
779 values = self._get_column(self._get_weather, row)
780 if row > -1:
781 return {'temperature': values[0],
782 'pressure': values[1], 'humidity' : values[2],
783 'windspeed' : values[3], 'windaz' : values[4]
784 }
785 else:
786 out = []
787 for r in values:
788
789 out.append({'temperature': r[0],
790 'pressure': r[1], 'humidity' : r[2],
791 'windspeed' : r[3], 'windaz' : r[4]
792 })
793 return out
794
795 def _row_callback(self, callback, label):
796 out = ""
797 outvec = []
798 sep = '-'*50
799 for i in range(self.nrow()):
800 tm = self._gettime(i)
801 src = self._getsourcename(i)
802 out += 'Scan[%d] (%s) ' % (self.getscan(i), src)
803 out += 'Time[%s]:\n' % (tm)
804 if self.nbeam(-1) > 1:
805 out += ' Beam[%d] ' % (self.getbeam(i))
806 if self.nif(-1) > 1: out += ' IF[%d] ' % (self.getif(i))
807 if self.npol(-1) > 1: out += ' Pol[%d] ' % (self.getpol(i))
808 outvec.append(callback(i))
809 out += '= %3.3f\n' % (outvec[i])
810 out += sep+'\n'
811
812 asaplog.push(sep)
813 asaplog.push(" %s" % (label))
814 asaplog.push(sep)
815 asaplog.push(out)
816 asaplog.post()
817 return outvec
818
819 def _get_column(self, callback, row=-1, *args):
820 """
821 """
822 if row == -1:
823 return [callback(i, *args) for i in range(self.nrow())]
824 else:
825 if 0 <= row < self.nrow():
826 return callback(row, *args)
827
828
829 def get_time(self, row=-1, asdatetime=False, prec=-1):
830 """\
831 Get a list of time stamps for the observations.
832 Return a datetime object or a string (default) for each
833 integration time stamp in the scantable.
834
835 Parameters:
836
837 row: row no of integration. Default -1 return all rows
838
839 asdatetime: return values as datetime objects rather than strings
840
841 prec: number of digits shown. Default -1 to automatic
842 calculation.
843 Note this number is equals to the digits of MVTime,
844 i.e., 0<prec<3: dates with hh:: only,
845 <5: with hh:mm:, <7 or 0: with hh:mm:ss,
846 and 6> : with hh:mm:ss.tt... (prec-6 t's added)
847
848 """
849 from datetime import datetime
850 if prec < 0:
851 # automagically set necessary precision +1
852 prec = 7 - \
853 numpy.floor(numpy.log10(numpy.min(self.get_inttime(row))))
854 prec = max(6, int(prec))
855 else:
856 prec = max(0, prec)
857 if asdatetime:
858 #precision can be 1 millisecond at max
859 prec = min(12, prec)
860
861 times = self._get_column(self._gettime, row, prec)
862 if not asdatetime:
863 return times
864 format = "%Y/%m/%d/%H:%M:%S.%f"
865 if prec < 7:
866 nsub = 1 + (((6-prec)/2) % 3)
867 substr = [".%f","%S","%M"]
868 for i in range(nsub):
869 format = format.replace(substr[i],"")
870 if isinstance(times, list):
871 return [datetime.strptime(i, format) for i in times]
872 else:
873 return datetime.strptime(times, format)
874
875
876 def get_inttime(self, row=-1):
877 """\
878 Get a list of integration times for the observations.
879 Return a time in seconds for each integration in the scantable.
880
881 Parameters:
882
883 row: row no of integration. Default -1 return all rows.
884
885 """
886 return self._get_column(self._getinttime, row)
887
888
889 def get_sourcename(self, row=-1):
890 """\
891 Get a list source names for the observations.
892 Return a string for each integration in the scantable.
893 Parameters:
894
895 row: row no of integration. Default -1 return all rows.
896
897 """
898 return self._get_column(self._getsourcename, row)
899
900 def get_elevation(self, row=-1):
901 """\
902 Get a list of elevations for the observations.
903 Return a float for each integration in the scantable.
904
905 Parameters:
906
907 row: row no of integration. Default -1 return all rows.
908
909 """
910 return self._get_column(self._getelevation, row)
911
912 def get_azimuth(self, row=-1):
913 """\
914 Get a list of azimuths for the observations.
915 Return a float for each integration in the scantable.
916
917 Parameters:
918 row: row no of integration. Default -1 return all rows.
919
920 """
921 return self._get_column(self._getazimuth, row)
922
923 def get_parangle(self, row=-1):
924 """\
925 Get a list of parallactic angles for the observations.
926 Return a float for each integration in the scantable.
927
928 Parameters:
929
930 row: row no of integration. Default -1 return all rows.
931
932 """
933 return self._get_column(self._getparangle, row)
934
935 def get_direction(self, row=-1):
936 """
937 Get a list of Positions on the sky (direction) for the observations.
938 Return a string for each integration in the scantable.
939
940 Parameters:
941
942 row: row no of integration. Default -1 return all rows
943
944 """
945 return self._get_column(self._getdirection, row)
946
947 def get_directionval(self, row=-1):
948 """\
949 Get a list of Positions on the sky (direction) for the observations.
950 Return a float for each integration in the scantable.
951
952 Parameters:
953
954 row: row no of integration. Default -1 return all rows
955
956 """
957 return self._get_column(self._getdirectionvec, row)
958
959 @asaplog_post_dec
960 def set_unit(self, unit='channel'):
961 """\
962 Set the unit for all following operations on this scantable
963
964 Parameters:
965
966 unit: optional unit, default is 'channel'. Use one of '*Hz',
967 'km/s', 'channel' or equivalent ''
968
969 """
970 varlist = vars()
971 if unit in ['', 'pixel', 'channel']:
972 unit = ''
973 inf = list(self._getcoordinfo())
974 inf[0] = unit
975 self._setcoordinfo(inf)
976 self._add_history("set_unit", varlist)
977
978 @asaplog_post_dec
979 def set_instrument(self, instr):
980 """\
981 Set the instrument for subsequent processing.
982
983 Parameters:
984
985 instr: Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
986 'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
987
988 """
989 self._setInstrument(instr)
990 self._add_history("set_instument", vars())
991
992 @asaplog_post_dec
993 def set_feedtype(self, feedtype):
994 """\
995 Overwrite the feed type, which might not be set correctly.
996
997 Parameters:
998
999 feedtype: 'linear' or 'circular'
1000
1001 """
1002 self._setfeedtype(feedtype)
1003 self._add_history("set_feedtype", vars())
1004
1005 @asaplog_post_dec
1006 def set_doppler(self, doppler='RADIO'):
1007 """\
1008 Set the doppler for all following operations on this scantable.
1009
1010 Parameters:
1011
1012 doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
1013
1014 """
1015 varlist = vars()
1016 inf = list(self._getcoordinfo())
1017 inf[2] = doppler
1018 self._setcoordinfo(inf)
1019 self._add_history("set_doppler", vars())
1020
1021 @asaplog_post_dec
1022 def set_freqframe(self, frame=None):
1023 """\
1024 Set the frame type of the Spectral Axis.
1025
1026 Parameters:
1027
1028 frame: an optional frame type, default 'LSRK'. Valid frames are:
1029 'TOPO', 'LSRD', 'LSRK', 'BARY',
1030 'GEO', 'GALACTO', 'LGROUP', 'CMB'
1031
1032 Example::
1033
1034 scan.set_freqframe('BARY')
1035
1036 """
1037 frame = frame or rcParams['scantable.freqframe']
1038 varlist = vars()
1039 # "REST" is not implemented in casacore
1040 #valid = ['REST', 'TOPO', 'LSRD', 'LSRK', 'BARY', \
1041 # 'GEO', 'GALACTO', 'LGROUP', 'CMB']
1042 valid = ['TOPO', 'LSRD', 'LSRK', 'BARY', \
1043 'GEO', 'GALACTO', 'LGROUP', 'CMB']
1044
1045 if frame in valid:
1046 inf = list(self._getcoordinfo())
1047 inf[1] = frame
1048 self._setcoordinfo(inf)
1049 self._add_history("set_freqframe", varlist)
1050 else:
1051 msg = "Please specify a valid freq type. Valid types are:\n", valid
1052 raise TypeError(msg)
1053
1054 @asaplog_post_dec
1055 def set_dirframe(self, frame=""):
1056 """\
1057 Set the frame type of the Direction on the sky.
1058
1059 Parameters:
1060
1061 frame: an optional frame type, default ''. Valid frames are:
1062 'J2000', 'B1950', 'GALACTIC'
1063
1064 Example:
1065
1066 scan.set_dirframe('GALACTIC')
1067
1068 """
1069 varlist = vars()
1070 Scantable.set_dirframe(self, frame)
1071 self._add_history("set_dirframe", varlist)
1072
1073 def get_unit(self):
1074 """\
1075 Get the default unit set in this scantable
1076
1077 Returns:
1078
1079 A unit string
1080
1081 """
1082 inf = self._getcoordinfo()
1083 unit = inf[0]
1084 if unit == '': unit = 'channel'
1085 return unit
1086
1087 @asaplog_post_dec
1088 def get_abcissa(self, rowno=0):
1089 """\
1090 Get the abcissa in the current coordinate setup for the currently
1091 selected Beam/IF/Pol
1092
1093 Parameters:
1094
1095 rowno: an optional row number in the scantable. Default is the
1096 first row, i.e. rowno=0
1097
1098 Returns:
1099
1100 The abcissa values and the format string (as a dictionary)
1101
1102 """
1103 abc = self._getabcissa(rowno)
1104 lbl = self._getabcissalabel(rowno)
1105 return abc, lbl
1106
1107 @asaplog_post_dec
1108 def flag(self, mask=None, unflag=False, row=-1):
1109 """\
1110 Flag the selected data using an optional channel mask.
1111
1112 Parameters:
1113
1114 mask: an optional channel mask, created with create_mask. Default
1115 (no mask) is all channels.
1116
1117 unflag: if True, unflag the data
1118
1119 row: an optional row number in the scantable.
1120 Default -1 flags all rows
1121
1122 """
1123 varlist = vars()
1124 mask = mask or []
1125 self._flag(row, mask, unflag)
1126 self._add_history("flag", varlist)
1127
1128 @asaplog_post_dec
1129 def flag_row(self, rows=None, unflag=False):
1130 """\
1131 Flag the selected data in row-based manner.
1132
1133 Parameters:
1134
1135 rows: list of row numbers to be flagged. Default is no row
1136 (must be explicitly specified to execute row-based
1137 flagging).
1138
1139 unflag: if True, unflag the data.
1140
1141 """
1142 varlist = vars()
1143 if rows is None:
1144 rows = []
1145 self._flag_row(rows, unflag)
1146 self._add_history("flag_row", varlist)
1147
1148 @asaplog_post_dec
1149 def clip(self, uthres=None, dthres=None, clipoutside=True, unflag=False):
1150 """\
1151 Flag the selected data outside a specified range (in channel-base)
1152
1153 Parameters:
1154
1155 uthres: upper threshold.
1156
1157 dthres: lower threshold
1158
1159 clipoutside: True for flagging data outside the range
1160 [dthres:uthres].
1161 False for flagging data inside the range.
1162
1163 unflag: if True, unflag the data.
1164
1165 """
1166 varlist = vars()
1167 self._clip(uthres, dthres, clipoutside, unflag)
1168 self._add_history("clip", varlist)
1169
1170 @asaplog_post_dec
1171 def lag_flag(self, start, end, unit="MHz", insitu=None):
1172 """\
1173 Flag the data in 'lag' space by providing a frequency to remove.
1174 Flagged data in the scantable get interpolated over the region.
1175 No taper is applied.
1176
1177 Parameters:
1178
1179 start: the start frequency (really a period within the
1180 bandwidth) or period to remove
1181
1182 end: the end frequency or period to remove
1183
1184 unit: the frequency unit (default 'MHz') or '' for
1185 explicit lag channels
1186
1187 *Notes*:
1188
1189 It is recommended to flag edges of the band or strong
1190 signals beforehand.
1191
1192 """
1193 if insitu is None: insitu = rcParams['insitu']
1194 self._math._setinsitu(insitu)
1195 varlist = vars()
1196 base = { "GHz": 1000000000., "MHz": 1000000., "kHz": 1000., "Hz": 1.}
1197 if not (unit == "" or base.has_key(unit)):
1198 raise ValueError("%s is not a valid unit." % unit)
1199 if unit == "":
1200 s = scantable(self._math._lag_flag(self, start, end, "lags"))
1201 else:
1202 s = scantable(self._math._lag_flag(self, start*base[unit],
1203 end*base[unit], "frequency"))
1204 s._add_history("lag_flag", varlist)
1205 if insitu:
1206 self._assign(s)
1207 else:
1208 return s
1209
1210 @asaplog_post_dec
1211 def fft(self, rowno=None, mask=None, getrealimag=False):
1212 """\
1213 Apply FFT to the spectra.
1214 Flagged data in the scantable get interpolated over the region.
1215
1216 Parameters:
1217
1218 rowno: The row number(s) to be processed. int, list
1219 and tuple are accepted. By default (None), FFT
1220 is applied to the whole data.
1221
1222 mask: Auxiliary channel mask(s). Given as a boolean
1223 list, it is applied to all specified rows.
1224 A list of boolean lists can also be used to
1225 apply different masks. In the latter case, the
1226 length of 'mask' must be the same as that of
1227 'rowno'. The default is None.
1228
1229 getrealimag: If True, returns the real and imaginary part
1230 values of the complex results.
1231 If False (the default), returns the amplitude
1232 (absolute value) normalised with Ndata/2 and
1233 phase (argument, in unit of radian).
1234
1235 Returns:
1236
1237 A list of dictionaries containing the results for each spectrum.
1238 Each dictionary contains two values, the real and the imaginary
1239 parts when getrealimag = True, or the amplitude(absolute value)
1240 and the phase(argument) when getrealimag = False. The key for
1241 these values are 'real' and 'imag', or 'ampl' and 'phase',
1242 respectively.
1243 """
1244 if rowno is None:
1245 rowno = []
1246 if isinstance(rowno, int):
1247 rowno = [rowno]
1248 elif not (isinstance(rowno, list) or isinstance(rowno, tuple)):
1249 raise TypeError("The row number(s) must be int, list or tuple.")
1250 if len(rowno) == 0: rowno = [i for i in xrange(self.nrow())]
1251
1252 usecommonmask = True
1253
1254 if mask is None:
1255 mask = []
1256 if isinstance(mask, list) or isinstance(mask, tuple):
1257 if len(mask) == 0:
1258 mask = [[]]
1259 else:
1260 if isinstance(mask[0], bool):
1261 if len(mask) != self.nchan(self.getif(rowno[0])):
1262 raise ValueError("The spectra and the mask have "
1263 "different length.")
1264 mask = [mask]
1265 elif isinstance(mask[0], list) or isinstance(mask[0], tuple):
1266 usecommonmask = False
1267 if len(mask) != len(rowno):
1268 raise ValueError("When specifying masks for each "
1269 "spectrum, the numbers of them "
1270 "must be identical.")
1271 for i in xrange(mask):
1272 if len(mask[i]) != self.nchan(self.getif(rowno[i])):
1273 raise ValueError("The spectra and the mask have "
1274 "different length.")
1275 else:
1276 raise TypeError("The mask must be a boolean list or "
1277 "a list of boolean list.")
1278 else:
1279 raise TypeError("The mask must be a boolean list or a list of "
1280 "boolean list.")
1281
1282 res = []
1283
1284 imask = 0
1285 for whichrow in rowno:
1286 fspec = self._fft(whichrow, mask[imask], getrealimag)
1287 nspec = len(fspec)
1288
1289 i = 0
1290 v1 = []
1291 v2 = []
1292 reselem = {"real":[],"imag":[]} if getrealimag \
1293 else {"ampl":[],"phase":[]}
1294
1295 while (i < nspec):
1296 v1.append(fspec[i])
1297 v2.append(fspec[i+1])
1298 i += 2
1299
1300 if getrealimag:
1301 reselem["real"] += v1
1302 reselem["imag"] += v2
1303 else:
1304 reselem["ampl"] += v1
1305 reselem["phase"] += v2
1306
1307 res.append(reselem)
1308
1309 if not usecommonmask:
1310 imask += 1
1311
1312 return res
1313
1314 @asaplog_post_dec
1315 def create_mask(self, *args, **kwargs):
1316 """\
1317 Compute and return a mask based on [min, max] windows.
1318 The specified windows are to be INCLUDED, when the mask is
1319 applied.
1320
1321 Parameters:
1322
1323 [min, max], [min2, max2], ...
1324 Pairs of start/end points (inclusive)specifying the regions
1325 to be masked
1326
1327 invert: optional argument. If specified as True,
1328 return an inverted mask, i.e. the regions
1329 specified are EXCLUDED
1330
1331 row: create the mask using the specified row for
1332 unit conversions, default is row=0
1333 only necessary if frequency varies over rows.
1334
1335 Examples::
1336
1337 scan.set_unit('channel')
1338 # a)
1339 msk = scan.create_mask([400, 500], [800, 900])
1340 # masks everything outside 400 and 500
1341 # and 800 and 900 in the unit 'channel'
1342
1343 # b)
1344 msk = scan.create_mask([400, 500], [800, 900], invert=True)
1345 # masks the regions between 400 and 500
1346 # and 800 and 900 in the unit 'channel'
1347
1348 # c)
1349 #mask only channel 400
1350 msk = scan.create_mask([400])
1351
1352 """
1353 row = kwargs.get("row", 0)
1354 data = self._getabcissa(row)
1355 u = self._getcoordinfo()[0]
1356 if u == "":
1357 u = "channel"
1358 msg = "The current mask window unit is %s" % u
1359 i = self._check_ifs()
1360 if not i:
1361 msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
1362 asaplog.push(msg)
1363 n = len(data)
1364 msk = _n_bools(n, False)
1365 # test if args is a 'list' or a 'normal *args - UGLY!!!
1366
1367 ws = (isinstance(args[-1][-1], int)
1368 or isinstance(args[-1][-1], float)) and args or args[0]
1369 for window in ws:
1370 if len(window) == 1:
1371 window = [window[0], window[0]]
1372 if len(window) == 0 or len(window) > 2:
1373 raise ValueError("A window needs to be defined as "
1374 "[start(, end)]")
1375 if window[0] > window[1]:
1376 tmp = window[0]
1377 window[0] = window[1]
1378 window[1] = tmp
1379 for i in range(n):
1380 if data[i] >= window[0] and data[i] <= window[1]:
1381 msk[i] = True
1382 if kwargs.has_key('invert'):
1383 if kwargs.get('invert'):
1384 msk = mask_not(msk)
1385 return msk
1386
1387 def get_masklist(self, mask=None, row=0, silent=False):
1388 """\
1389 Compute and return a list of mask windows, [min, max].
1390
1391 Parameters:
1392
1393 mask: channel mask, created with create_mask.
1394
1395 row: calcutate the masklist using the specified row
1396 for unit conversions, default is row=0
1397 only necessary if frequency varies over rows.
1398
1399 Returns:
1400
1401 [min, max], [min2, max2], ...
1402 Pairs of start/end points (inclusive)specifying
1403 the masked regions
1404
1405 """
1406 if not (isinstance(mask,list) or isinstance(mask, tuple)):
1407 raise TypeError("The mask should be list or tuple.")
1408 if len(mask) <= 0:
1409 raise TypeError("The mask elements should be > 0")
1410 data = self._getabcissa(row)
1411 if len(data) != len(mask):
1412 msg = "Number of channels in scantable != number of mask elements"
1413 raise TypeError(msg)
1414 u = self._getcoordinfo()[0]
1415 if u == "":
1416 u = "channel"
1417 msg = "The current mask window unit is %s" % u
1418 i = self._check_ifs()
1419 if not i:
1420 msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
1421 if not silent:
1422 asaplog.push(msg)
1423 masklist = []
1424 ist, ien = None, None
1425 ist, ien=self.get_mask_indices(mask)
1426 if ist is not None and ien is not None:
1427 for i in xrange(len(ist)):
1428 range=[data[ist[i]],data[ien[i]]]
1429 range.sort()
1430 masklist.append([range[0],range[1]])
1431 return masklist
1432
1433 def get_mask_indices(self, mask=None):
1434 """\
1435 Compute and Return lists of mask start indices and mask end indices.
1436
1437 Parameters:
1438
1439 mask: channel mask, created with create_mask.
1440
1441 Returns:
1442
1443 List of mask start indices and that of mask end indices,
1444 i.e., [istart1,istart2,....], [iend1,iend2,....].
1445
1446 """
1447 if not (isinstance(mask,list) or isinstance(mask, tuple)):
1448 raise TypeError("The mask should be list or tuple.")
1449 if len(mask) <= 0:
1450 raise TypeError("The mask elements should be > 0")
1451 istart = []
1452 iend = []
1453 if mask[0]:
1454 istart.append(0)
1455 for i in range(len(mask)-1):
1456 if not mask[i] and mask[i+1]:
1457 istart.append(i+1)
1458 elif mask[i] and not mask[i+1]:
1459 iend.append(i)
1460 if mask[len(mask)-1]:
1461 iend.append(len(mask)-1)
1462 if len(istart) != len(iend):
1463 raise RuntimeError("Numbers of mask start != mask end.")
1464 for i in range(len(istart)):
1465 if istart[i] > iend[i]:
1466 raise RuntimeError("Mask start index > mask end index")
1467 break
1468 return istart,iend
1469
1470 @asaplog_post_dec
1471 def parse_maskexpr(self, maskstring):
1472 """
1473 Parse CASA type mask selection syntax (IF dependent).
1474
1475 Parameters:
1476 maskstring : A string mask selection expression.
1477 A comma separated selections mean different IF -
1478 channel combinations. IFs and channel selections
1479 are partitioned by a colon, ':'.
1480 examples:
1481 '' = all IFs (all channels)
1482 '<2,4~6,9' = IFs 0,1,4,5,6,9 (all channels)
1483 '3:3~45;60' = channels 3 to 45 and 60 in IF 3
1484 '0~1:2~6,8' = channels 2 to 6 in IFs 0,1, and
1485 all channels in IF8
1486 Returns:
1487 A dictionary of selected (valid) IF and masklist pairs,
1488 e.g. {'0': [[50,250],[350,462]], '2': [[100,400],[550,974]]}
1489 """
1490 if not isinstance(maskstring,str):
1491 asaplog.post()
1492 asaplog.push("Mask expression should be a string.")
1493 asaplog.post("ERROR")
1494
1495 valid_ifs = self.getifnos()
1496 frequnit = self.get_unit()
1497 seldict = {}
1498 if maskstring == "":
1499 maskstring = str(valid_ifs)[1:-1]
1500 ## split each selection "IF range[:CHAN range]"
1501 sellist = maskstring.split(',')
1502 for currselstr in sellist:
1503 selset = currselstr.split(':')
1504 # spw and mask string (may include ~, < or >)
1505 spwmasklist = self._parse_selection(selset[0], typestr='integer',
1506 minval=min(valid_ifs),
1507 maxval=max(valid_ifs))
1508 for spwlist in spwmasklist:
1509 selspws = []
1510 for ispw in range(spwlist[0],spwlist[1]+1):
1511 # Put into the list only if ispw exists
1512 if valid_ifs.count(ispw):
1513 selspws.append(ispw)
1514 del spwmasklist, spwlist
1515
1516 # parse frequency mask list
1517 if len(selset) > 1:
1518 freqmasklist = self._parse_selection(selset[1], typestr='float',
1519 offset=0.)
1520 else:
1521 # want to select the whole spectrum
1522 freqmasklist = [None]
1523
1524 ## define a dictionary of spw - masklist combination
1525 for ispw in selspws:
1526 #print "working on", ispw
1527 spwstr = str(ispw)
1528 if len(selspws) == 0:
1529 # empty spw
1530 continue
1531 else:
1532 ## want to get min and max of the spw and
1533 ## offset to set for '<' and '>'
1534 if frequnit == 'channel':
1535 minfreq = 0
1536 maxfreq = self.nchan(ifno=ispw)
1537 offset = 0.5
1538 else:
1539 ## This is ugly part. need improvement
1540 for ifrow in xrange(self.nrow()):
1541 if self.getif(ifrow) == ispw:
1542 #print "IF",ispw,"found in row =",ifrow
1543 break
1544 freqcoord = self.get_coordinate(ifrow)
1545 freqs = self._getabcissa(ifrow)
1546 minfreq = min(freqs)
1547 maxfreq = max(freqs)
1548 if len(freqs) == 1:
1549 offset = 0.5
1550 elif frequnit.find('Hz') > 0:
1551 offset = abs(freqcoord.to_frequency(1,
1552 unit=frequnit)
1553 -freqcoord.to_frequency(0,
1554 unit=frequnit)
1555 )*0.5
1556 elif frequnit.find('m/s') > 0:
1557 offset = abs(freqcoord.to_velocity(1,
1558 unit=frequnit)
1559 -freqcoord.to_velocity(0,
1560 unit=frequnit)
1561 )*0.5
1562 else:
1563 asaplog.post()
1564 asaplog.push("Invalid frequency unit")
1565 asaplog.post("ERROR")
1566 del freqs, freqcoord, ifrow
1567 for freq in freqmasklist:
1568 selmask = freq or [minfreq, maxfreq]
1569 if selmask[0] == None:
1570 ## selection was "<freq[1]".
1571 if selmask[1] < minfreq:
1572 ## avoid adding region selection
1573 selmask = None
1574 else:
1575 selmask = [minfreq,selmask[1]-offset]
1576 elif selmask[1] == None:
1577 ## selection was ">freq[0]"
1578 if selmask[0] > maxfreq:
1579 ## avoid adding region selection
1580 selmask = None
1581 else:
1582 selmask = [selmask[0]+offset,maxfreq]
1583 if selmask:
1584 if not seldict.has_key(spwstr):
1585 # new spw selection
1586 seldict[spwstr] = []
1587 seldict[spwstr] += [selmask]
1588 del minfreq,maxfreq,offset,freq,selmask
1589 del spwstr
1590 del freqmasklist
1591 del valid_ifs
1592 if len(seldict) == 0:
1593 asaplog.post()
1594 asaplog.push("No valid selection in the mask expression: "
1595 +maskstring)
1596 asaplog.post("WARN")
1597 return None
1598 msg = "Selected masklist:\n"
1599 for sif, lmask in seldict.iteritems():
1600 msg += " IF"+sif+" - "+str(lmask)+"\n"
1601 asaplog.push(msg)
1602 return seldict
1603
1604 @asaplog_post_dec
1605 def parse_idx_selection(self, mode, selexpr):
1606 """
1607 Parse CASA type mask selection syntax of SCANNO, IFNO, POLNO,
1608 BEAMNO, and row number
1609
1610 Parameters:
1611 mode : which column to select.
1612 ['scan',|'if'|'pol'|'beam'|'row']
1613 selexpr : A comma separated selection expression.
1614 examples:
1615 '' = all (returns [])
1616 '<2,4~6,9' = indices less than 2, 4 to 6 and 9
1617 (returns [0,1,4,5,6,9])
1618 Returns:
1619 A List of selected indices
1620 """
1621 if selexpr == "":
1622 return []
1623 valid_modes = {'s': 'scan', 'i': 'if', 'p': 'pol',
1624 'b': 'beam', 'r': 'row'}
1625 smode = mode.lower()[0]
1626 if not (smode in valid_modes.keys()):
1627 msg = "Invalid mode '%s'. Valid modes are %s" %\
1628 (mode, str(valid_modes.values()))
1629 asaplog.post()
1630 asaplog.push(msg)
1631 asaplog.post("ERROR")
1632 mode = valid_modes[smode]
1633 minidx = None
1634 maxidx = None
1635 if smode == 'r':
1636 minidx = 0
1637 maxidx = self.nrow()
1638 else:
1639 idx = getattr(self,"get"+mode+"nos")()
1640 minidx = min(idx)
1641 maxidx = max(idx)
1642 del idx
1643 sellist = selexpr.split(',')
1644 idxlist = []
1645 for currselstr in sellist:
1646 # single range (may include ~, < or >)
1647 currlist = self._parse_selection(currselstr, typestr='integer',
1648 minval=minidx,maxval=maxidx)
1649 for thelist in currlist:
1650 idxlist += range(thelist[0],thelist[1]+1)
1651 msg = "Selected %s: %s" % (mode.upper()+"NO", str(idxlist))
1652 asaplog.push(msg)
1653 return idxlist
1654
1655 def _parse_selection(self, selstr, typestr='float', offset=0.,
1656 minval=None, maxval=None):
1657 """
1658 Parameters:
1659 selstr : The Selection string, e.g., '<3;5~7;100~103;9'
1660 typestr : The type of the values in returned list
1661 ('integer' or 'float')
1662 offset : The offset value to subtract from or add to
1663 the boundary value if the selection string
1664 includes '<' or '>' [Valid only for typestr='float']
1665 minval, maxval : The minimum/maximum values to set if the
1666 selection string includes '<' or '>'.
1667 The list element is filled with None by default.
1668 Returns:
1669 A list of min/max pair of selections.
1670 Example:
1671 _parse_selection('<3;5~7;9',typestr='int',minval=0)
1672 --> returns [[0,2],[5,7],[9,9]]
1673 _parse_selection('<3;5~7;9',typestr='float',offset=0.5,minval=0)
1674 --> returns [[0.,2.5],[5.0,7.0],[9.,9.]]
1675 """
1676 selgroups = selstr.split(';')
1677 sellists = []
1678 if typestr.lower().startswith('int'):
1679 formatfunc = int
1680 offset = 1
1681 else:
1682 formatfunc = float
1683
1684 for currsel in selgroups:
1685 if currsel.find('~') > 0:
1686 # val0 <= x <= val1
1687 minsel = formatfunc(currsel.split('~')[0].strip())
1688 maxsel = formatfunc(currsel.split('~')[1].strip())
1689 elif currsel.strip().find('<=') > -1:
1690 bound = currsel.split('<=')
1691 try: # try "x <= val"
1692 minsel = minval
1693 maxsel = formatfunc(bound[1].strip())
1694 except ValueError: # now "val <= x"
1695 minsel = formatfunc(bound[0].strip())
1696 maxsel = maxval
1697 elif currsel.strip().find('>=') > -1:
1698 bound = currsel.split('>=')
1699 try: # try "x >= val"
1700 minsel = formatfunc(bound[1].strip())
1701 maxsel = maxval
1702 except ValueError: # now "val >= x"
1703 minsel = minval
1704 maxsel = formatfunc(bound[0].strip())
1705 elif currsel.strip().find('<') > -1:
1706 bound = currsel.split('<')
1707 try: # try "x < val"
1708 minsel = minval
1709 maxsel = formatfunc(bound[1].strip()) \
1710 - formatfunc(offset)
1711 except ValueError: # now "val < x"
1712 minsel = formatfunc(bound[0].strip()) \
1713 + formatfunc(offset)
1714 maxsel = maxval
1715 elif currsel.strip().find('>') > -1:
1716 bound = currsel.split('>')
1717 try: # try "x > val"
1718 minsel = formatfunc(bound[1].strip()) \
1719 + formatfunc(offset)
1720 maxsel = maxval
1721 except ValueError: # now "val > x"
1722 minsel = minval
1723 maxsel = formatfunc(bound[0].strip()) \
1724 - formatfunc(offset)
1725 else:
1726 minsel = formatfunc(currsel)
1727 maxsel = formatfunc(currsel)
1728 sellists.append([minsel,maxsel])
1729 return sellists
1730
1731# def get_restfreqs(self):
1732# """
1733# Get the restfrequency(s) stored in this scantable.
1734# The return value(s) are always of unit 'Hz'
1735# Parameters:
1736# none
1737# Returns:
1738# a list of doubles
1739# """
1740# return list(self._getrestfreqs())
1741
1742 def get_restfreqs(self, ids=None):
1743 """\
1744 Get the restfrequency(s) stored in this scantable.
1745 The return value(s) are always of unit 'Hz'
1746
1747 Parameters:
1748
1749 ids: (optional) a list of MOLECULE_ID for that restfrequency(s) to
1750 be retrieved
1751
1752 Returns:
1753
1754 dictionary containing ids and a list of doubles for each id
1755
1756 """
1757 if ids is None:
1758 rfreqs = {}
1759 idlist = self.getmolnos()
1760 for i in idlist:
1761 rfreqs[i] = list(self._getrestfreqs(i))
1762 return rfreqs
1763 else:
1764 if type(ids) == list or type(ids) == tuple:
1765 rfreqs = {}
1766 for i in ids:
1767 rfreqs[i] = list(self._getrestfreqs(i))
1768 return rfreqs
1769 else:
1770 return list(self._getrestfreqs(ids))
1771
1772 @asaplog_post_dec
1773 def set_restfreqs(self, freqs=None, unit='Hz'):
1774 """\
1775 Set or replace the restfrequency specified and
1776 if the 'freqs' argument holds a scalar,
1777 then that rest frequency will be applied to all the selected
1778 data. If the 'freqs' argument holds
1779 a vector, then it MUST be of equal or smaller length than
1780 the number of IFs (and the available restfrequencies will be
1781 replaced by this vector). In this case, *all* data have
1782 the restfrequency set per IF according
1783 to the corresponding value you give in the 'freqs' vector.
1784 E.g. 'freqs=[1e9, 2e9]' would mean IF 0 gets restfreq 1e9 and
1785 IF 1 gets restfreq 2e9.
1786
1787 You can also specify the frequencies via a linecatalog.
1788
1789 Parameters:
1790
1791 freqs: list of rest frequency values or string idenitfiers
1792
1793 unit: unit for rest frequency (default 'Hz')
1794
1795
1796 Example::
1797
1798 # set the given restfrequency for the all currently selected IFs
1799 scan.set_restfreqs(freqs=1.4e9)
1800 # set restfrequencies for the n IFs (n > 1) in the order of the
1801 # list, i.e
1802 # IF0 -> 1.4e9, IF1 -> 1.41e9, IF3 -> 1.42e9
1803 # len(list_of_restfreqs) == nIF
1804 # for nIF == 1 the following will set multiple restfrequency for
1805 # that IF
1806 scan.set_restfreqs(freqs=[1.4e9, 1.41e9, 1.42e9])
1807 # set multiple restfrequencies per IF. as a list of lists where
1808 # the outer list has nIF elements, the inner s arbitrary
1809 scan.set_restfreqs(freqs=[[1.4e9, 1.41e9], [1.67e9]])
1810
1811 *Note*:
1812
1813 To do more sophisticate Restfrequency setting, e.g. on a
1814 source and IF basis, use scantable.set_selection() before using
1815 this function::
1816
1817 # provided your scantable is called scan
1818 selection = selector()
1819 selection.set_name('ORION*')
1820 selection.set_ifs([1])
1821 scan.set_selection(selection)
1822 scan.set_restfreqs(freqs=86.6e9)
1823
1824 """
1825 varlist = vars()
1826 from asap import linecatalog
1827 # simple value
1828 if isinstance(freqs, int) or isinstance(freqs, float):
1829 self._setrestfreqs([freqs], [""], unit)
1830 # list of values
1831 elif isinstance(freqs, list) or isinstance(freqs, tuple):
1832 # list values are scalars
1833 if isinstance(freqs[-1], int) or isinstance(freqs[-1], float):
1834 if len(freqs) == 1:
1835 self._setrestfreqs(freqs, [""], unit)
1836 else:
1837 # allow the 'old' mode of setting mulitple IFs
1838 savesel = self._getselection()
1839 sel = self.get_selection()
1840 iflist = self.getifnos()
1841 if len(freqs)>len(iflist):
1842 raise ValueError("number of elements in list of list "
1843 "exeeds the current IF selections")
1844 iflist = self.getifnos()
1845 for i, fval in enumerate(freqs):
1846 sel.set_ifs(iflist[i])
1847 self._setselection(sel)
1848 self._setrestfreqs([fval], [""], unit)
1849 self._setselection(savesel)
1850
1851 # list values are dict, {'value'=, 'name'=)
1852 elif isinstance(freqs[-1], dict):
1853 values = []
1854 names = []
1855 for d in freqs:
1856 values.append(d["value"])
1857 names.append(d["name"])
1858 self._setrestfreqs(values, names, unit)
1859 elif isinstance(freqs[-1], list) or isinstance(freqs[-1], tuple):
1860 savesel = self._getselection()
1861 sel = self.get_selection()
1862 iflist = self.getifnos()
1863 if len(freqs)>len(iflist):
1864 raise ValueError("number of elements in list of list exeeds"
1865 " the current IF selections")
1866 for i, fval in enumerate(freqs):
1867 sel.set_ifs(iflist[i])
1868 self._setselection(sel)
1869 self._setrestfreqs(fval, [""], unit)
1870 self._setselection(savesel)
1871 # freqs are to be taken from a linecatalog
1872 elif isinstance(freqs, linecatalog):
1873 savesel = self._getselection()
1874 sel = self.get_selection()
1875 for i in xrange(freqs.nrow()):
1876 sel.set_ifs(iflist[i])
1877 self._setselection(sel)
1878 self._setrestfreqs([freqs.get_frequency(i)],
1879 [freqs.get_name(i)], "MHz")
1880 # ensure that we are not iterating past nIF
1881 if i == self.nif()-1: break
1882 self._setselection(savesel)
1883 else:
1884 return
1885 self._add_history("set_restfreqs", varlist)
1886
1887 @asaplog_post_dec
1888 def shift_refpix(self, delta):
1889 """\
1890 Shift the reference pixel of the Spectra Coordinate by an
1891 integer amount.
1892
1893 Parameters:
1894
1895 delta: the amount to shift by
1896
1897 *Note*:
1898
1899 Be careful using this with broadband data.
1900
1901 """
1902 varlist = vars()
1903 Scantable.shift_refpix(self, delta)
1904 s._add_history("shift_refpix", varlist)
1905
1906 @asaplog_post_dec
1907 def history(self, filename=None):
1908 """\
1909 Print the history. Optionally to a file.
1910
1911 Parameters:
1912
1913 filename: The name of the file to save the history to.
1914
1915 """
1916 hist = list(self._gethistory())
1917 out = "-"*80
1918 for h in hist:
1919 if h.startswith("---"):
1920 out = "\n".join([out, h])
1921 else:
1922 items = h.split("##")
1923 date = items[0]
1924 func = items[1]
1925 items = items[2:]
1926 out += "\n"+date+"\n"
1927 out += "Function: %s\n Parameters:" % (func)
1928 for i in items:
1929 if i == '':
1930 continue
1931 s = i.split("=")
1932 out += "\n %s = %s" % (s[0], s[1])
1933 out = "\n".join([out, "-"*80])
1934 if filename is not None:
1935 if filename is "":
1936 filename = 'scantable_history.txt'
1937 import os
1938 filename = os.path.expandvars(os.path.expanduser(filename))
1939 if not os.path.isdir(filename):
1940 data = open(filename, 'w')
1941 data.write(out)
1942 data.close()
1943 else:
1944 msg = "Illegal file name '%s'." % (filename)
1945 raise IOError(msg)
1946 return page(out)
1947
1948 #
1949 # Maths business
1950 #
1951 @asaplog_post_dec
1952 def average_time(self, mask=None, scanav=False, weight='tint', align=False):
1953 """\
1954 Return the (time) weighted average of a scan. Scans will be averaged
1955 only if the source direction (RA/DEC) is within 1' otherwise
1956
1957 *Note*:
1958
1959 in channels only - align if necessary
1960
1961 Parameters:
1962
1963 mask: an optional mask (only used for 'var' and 'tsys'
1964 weighting)
1965
1966 scanav: True averages each scan separately
1967 False (default) averages all scans together,
1968
1969 weight: Weighting scheme.
1970 'none' (mean no weight)
1971 'var' (1/var(spec) weighted)
1972 'tsys' (1/Tsys**2 weighted)
1973 'tint' (integration time weighted)
1974 'tintsys' (Tint/Tsys**2)
1975 'median' ( median averaging)
1976 The default is 'tint'
1977
1978 align: align the spectra in velocity before averaging. It takes
1979 the time of the first spectrum as reference time.
1980
1981 Example::
1982
1983 # time average the scantable without using a mask
1984 newscan = scan.average_time()
1985
1986 """
1987 varlist = vars()
1988 weight = weight or 'TINT'
1989 mask = mask or ()
1990 scanav = (scanav and 'SCAN') or 'NONE'
1991 scan = (self, )
1992
1993 if align:
1994 scan = (self.freq_align(insitu=False), )
1995 s = None
1996 if weight.upper() == 'MEDIAN':
1997 s = scantable(self._math._averagechannel(scan[0], 'MEDIAN',
1998 scanav))
1999 else:
2000 s = scantable(self._math._average(scan, mask, weight.upper(),
2001 scanav))
2002 s._add_history("average_time", varlist)
2003 return s
2004
2005 @asaplog_post_dec
2006 def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None):
2007 """\
2008 Return a scan where all spectra are converted to either
2009 Jansky or Kelvin depending upon the flux units of the scan table.
2010 By default the function tries to look the values up internally.
2011 If it can't find them (or if you want to over-ride), you must
2012 specify EITHER jyperk OR eta (and D which it will try to look up
2013 also if you don't set it). jyperk takes precedence if you set both.
2014
2015 Parameters:
2016
2017 jyperk: the Jy / K conversion factor
2018
2019 eta: the aperture efficiency
2020
2021 d: the geometric diameter (metres)
2022
2023 insitu: if False a new scantable is returned.
2024 Otherwise, the scaling is done in-situ
2025 The default is taken from .asaprc (False)
2026
2027 """
2028 if insitu is None: insitu = rcParams['insitu']
2029 self._math._setinsitu(insitu)
2030 varlist = vars()
2031 jyperk = jyperk or -1.0
2032 d = d or -1.0
2033 eta = eta or -1.0
2034 s = scantable(self._math._convertflux(self, d, eta, jyperk))
2035 s._add_history("convert_flux", varlist)
2036 if insitu: self._assign(s)
2037 else: return s
2038
2039 @asaplog_post_dec
2040 def gain_el(self, poly=None, filename="", method="linear", insitu=None):
2041 """\
2042 Return a scan after applying a gain-elevation correction.
2043 The correction can be made via either a polynomial or a
2044 table-based interpolation (and extrapolation if necessary).
2045 You specify polynomial coefficients, an ascii table or neither.
2046 If you specify neither, then a polynomial correction will be made
2047 with built in coefficients known for certain telescopes (an error
2048 will occur if the instrument is not known).
2049 The data and Tsys are *divided* by the scaling factors.
2050
2051 Parameters:
2052
2053 poly: Polynomial coefficients (default None) to compute a
2054 gain-elevation correction as a function of
2055 elevation (in degrees).
2056
2057 filename: The name of an ascii file holding correction factors.
2058 The first row of the ascii file must give the column
2059 names and these MUST include columns
2060 'ELEVATION' (degrees) and 'FACTOR' (multiply data
2061 by this) somewhere.
2062 The second row must give the data type of the
2063 column. Use 'R' for Real and 'I' for Integer.
2064 An example file would be
2065 (actual factors are arbitrary) :
2066
2067 TIME ELEVATION FACTOR
2068 R R R
2069 0.1 0 0.8
2070 0.2 20 0.85
2071 0.3 40 0.9
2072 0.4 60 0.85
2073 0.5 80 0.8
2074 0.6 90 0.75
2075
2076 method: Interpolation method when correcting from a table.
2077 Values are 'nearest', 'linear' (default), 'cubic'
2078 and 'spline'
2079
2080 insitu: if False a new scantable is returned.
2081 Otherwise, the scaling is done in-situ
2082 The default is taken from .asaprc (False)
2083
2084 """
2085
2086 if insitu is None: insitu = rcParams['insitu']
2087 self._math._setinsitu(insitu)
2088 varlist = vars()
2089 poly = poly or ()
2090 from os.path import expandvars
2091 filename = expandvars(filename)
2092 s = scantable(self._math._gainel(self, poly, filename, method))
2093 s._add_history("gain_el", varlist)
2094 if insitu:
2095 self._assign(s)
2096 else:
2097 return s
2098
2099 @asaplog_post_dec
2100 def freq_align(self, reftime=None, method='cubic', insitu=None):
2101 """\
2102 Return a scan where all rows have been aligned in frequency/velocity.
2103 The alignment frequency frame (e.g. LSRK) is that set by function
2104 set_freqframe.
2105
2106 Parameters:
2107
2108 reftime: reference time to align at. By default, the time of
2109 the first row of data is used.
2110
2111 method: Interpolation method for regridding the spectra.
2112 Choose from 'nearest', 'linear', 'cubic' (default)
2113 and 'spline'
2114
2115 insitu: if False a new scantable is returned.
2116 Otherwise, the scaling is done in-situ
2117 The default is taken from .asaprc (False)
2118
2119 """
2120 if insitu is None: insitu = rcParams["insitu"]
2121 oldInsitu = self._math._insitu()
2122 self._math._setinsitu(insitu)
2123 varlist = vars()
2124 reftime = reftime or ""
2125 s = scantable(self._math._freq_align(self, reftime, method))
2126 s._add_history("freq_align", varlist)
2127 self._math._setinsitu(oldInsitu)
2128 if insitu:
2129 self._assign(s)
2130 else:
2131 return s
2132
2133 @asaplog_post_dec
2134 def opacity(self, tau=None, insitu=None):
2135 """\
2136 Apply an opacity correction. The data
2137 and Tsys are multiplied by the correction factor.
2138
2139 Parameters:
2140
2141 tau: (list of) opacity from which the correction factor is
2142 exp(tau*ZD)
2143 where ZD is the zenith-distance.
2144 If a list is provided, it has to be of length nIF,
2145 nIF*nPol or 1 and in order of IF/POL, e.g.
2146 [opif0pol0, opif0pol1, opif1pol0 ...]
2147 if tau is `None` the opacities are determined from a
2148 model.
2149
2150 insitu: if False a new scantable is returned.
2151 Otherwise, the scaling is done in-situ
2152 The default is taken from .asaprc (False)
2153
2154 """
2155 if insitu is None:
2156 insitu = rcParams['insitu']
2157 self._math._setinsitu(insitu)
2158 varlist = vars()
2159 if not hasattr(tau, "__len__"):
2160 tau = [tau]
2161 s = scantable(self._math._opacity(self, tau))
2162 s._add_history("opacity", varlist)
2163 if insitu:
2164 self._assign(s)
2165 else:
2166 return s
2167
2168 @asaplog_post_dec
2169 def bin(self, width=5, insitu=None):
2170 """\
2171 Return a scan where all spectra have been binned up.
2172
2173 Parameters:
2174
2175 width: The bin width (default=5) in pixels
2176
2177 insitu: if False a new scantable is returned.
2178 Otherwise, the scaling is done in-situ
2179 The default is taken from .asaprc (False)
2180
2181 """
2182 if insitu is None:
2183 insitu = rcParams['insitu']
2184 self._math._setinsitu(insitu)
2185 varlist = vars()
2186 s = scantable(self._math._bin(self, width))
2187 s._add_history("bin", varlist)
2188 if insitu:
2189 self._assign(s)
2190 else:
2191 return s
2192
2193 @asaplog_post_dec
2194 def resample(self, width=5, method='cubic', insitu=None):
2195 """\
2196 Return a scan where all spectra have been binned up.
2197
2198 Parameters:
2199
2200 width: The bin width (default=5) in pixels
2201
2202 method: Interpolation method when correcting from a table.
2203 Values are 'nearest', 'linear', 'cubic' (default)
2204 and 'spline'
2205
2206 insitu: if False a new scantable is returned.
2207 Otherwise, the scaling is done in-situ
2208 The default is taken from .asaprc (False)
2209
2210 """
2211 if insitu is None:
2212 insitu = rcParams['insitu']
2213 self._math._setinsitu(insitu)
2214 varlist = vars()
2215 s = scantable(self._math._resample(self, method, width))
2216 s._add_history("resample", varlist)
2217 if insitu:
2218 self._assign(s)
2219 else:
2220 return s
2221
2222 @asaplog_post_dec
2223 def average_pol(self, mask=None, weight='none'):
2224 """\
2225 Average the Polarisations together.
2226
2227 Parameters:
2228
2229 mask: An optional mask defining the region, where the
2230 averaging will be applied. The output will have all
2231 specified points masked.
2232
2233 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
2234 weighted), or 'tsys' (1/Tsys**2 weighted)
2235
2236 """
2237 varlist = vars()
2238 mask = mask or ()
2239 s = scantable(self._math._averagepol(self, mask, weight.upper()))
2240 s._add_history("average_pol", varlist)
2241 return s
2242
2243 @asaplog_post_dec
2244 def average_beam(self, mask=None, weight='none'):
2245 """\
2246 Average the Beams together.
2247
2248 Parameters:
2249 mask: An optional mask defining the region, where the
2250 averaging will be applied. The output will have all
2251 specified points masked.
2252
2253 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
2254 weighted), or 'tsys' (1/Tsys**2 weighted)
2255
2256 """
2257 varlist = vars()
2258 mask = mask or ()
2259 s = scantable(self._math._averagebeams(self, mask, weight.upper()))
2260 s._add_history("average_beam", varlist)
2261 return s
2262
2263 def parallactify(self, pflag):
2264 """\
2265 Set a flag to indicate whether this data should be treated as having
2266 been 'parallactified' (total phase == 0.0)
2267
2268 Parameters:
2269
2270 pflag: Bool indicating whether to turn this on (True) or
2271 off (False)
2272
2273 """
2274 varlist = vars()
2275 self._parallactify(pflag)
2276 self._add_history("parallactify", varlist)
2277
2278 @asaplog_post_dec
2279 def convert_pol(self, poltype=None):
2280 """\
2281 Convert the data to a different polarisation type.
2282 Note that you will need cross-polarisation terms for most conversions.
2283
2284 Parameters:
2285
2286 poltype: The new polarisation type. Valid types are:
2287 'linear', 'circular', 'stokes' and 'linpol'
2288
2289 """
2290 varlist = vars()
2291 s = scantable(self._math._convertpol(self, poltype))
2292 s._add_history("convert_pol", varlist)
2293 return s
2294
2295 @asaplog_post_dec
2296 def smooth(self, kernel="hanning", width=5.0, order=2, plot=False,
2297 insitu=None):
2298 """\
2299 Smooth the spectrum by the specified kernel (conserving flux).
2300
2301 Parameters:
2302
2303 kernel: The type of smoothing kernel. Select from
2304 'hanning' (default), 'gaussian', 'boxcar', 'rmedian'
2305 or 'poly'
2306
2307 width: The width of the kernel in pixels. For hanning this is
2308 ignored otherwise it defauls to 5 pixels.
2309 For 'gaussian' it is the Full Width Half
2310 Maximum. For 'boxcar' it is the full width.
2311 For 'rmedian' and 'poly' it is the half width.
2312
2313 order: Optional parameter for 'poly' kernel (default is 2), to
2314 specify the order of the polnomial. Ignored by all other
2315 kernels.
2316
2317 plot: plot the original and the smoothed spectra.
2318 In this each indivual fit has to be approved, by
2319 typing 'y' or 'n'
2320
2321 insitu: if False a new scantable is returned.
2322 Otherwise, the scaling is done in-situ
2323 The default is taken from .asaprc (False)
2324
2325 """
2326 if insitu is None: insitu = rcParams['insitu']
2327 self._math._setinsitu(insitu)
2328 varlist = vars()
2329
2330 if plot: orgscan = self.copy()
2331
2332 s = scantable(self._math._smooth(self, kernel.lower(), width, order))
2333 s._add_history("smooth", varlist)
2334
2335 action = 'H'
2336 if plot:
2337 from asap.asapplotter import new_asaplot
2338 theplot = new_asaplot(rcParams['plotter.gui'])
2339 from matplotlib import rc as rcp
2340 rcp('lines', linewidth=1)
2341 theplot.set_panels()
2342 ylab=s._get_ordinate_label()
2343 #theplot.palette(0,["#777777","red"])
2344 for r in xrange(s.nrow()):
2345 xsm=s._getabcissa(r)
2346 ysm=s._getspectrum(r)
2347 xorg=orgscan._getabcissa(r)
2348 yorg=orgscan._getspectrum(r)
2349 if action != "N": #skip plotting if rejecting all
2350 theplot.clear()
2351 theplot.hold()
2352 theplot.set_axes('ylabel',ylab)
2353 theplot.set_axes('xlabel',s._getabcissalabel(r))
2354 theplot.set_axes('title',s._getsourcename(r))
2355 theplot.set_line(label='Original',color="#777777")
2356 theplot.plot(xorg,yorg)
2357 theplot.set_line(label='Smoothed',color="red")
2358 theplot.plot(xsm,ysm)
2359 ### Ugly part for legend
2360 for i in [0,1]:
2361 theplot.subplots[0]['lines'].append(
2362 [theplot.subplots[0]['axes'].lines[i]]
2363 )
2364 theplot.release()
2365 ### Ugly part for legend
2366 theplot.subplots[0]['lines']=[]
2367 res = self._get_verify_action("Accept smoothing?",action)
2368 #print "IF%d, POL%d: got result = %s" %(s.getif(r),s.getpol(r),res)
2369 if r == 0: action = None
2370 #res = raw_input("Accept smoothing ([y]/n): ")
2371 if res.upper() == 'N':
2372 # reject for the current rows
2373 s._setspectrum(yorg, r)
2374 elif res.upper() == 'R':
2375 # reject all the following rows
2376 action = "N"
2377 s._setspectrum(yorg, r)
2378 elif res.upper() == 'A':
2379 # accept all the following rows
2380 break
2381 theplot.quit()
2382 del theplot
2383 del orgscan
2384
2385 if insitu: self._assign(s)
2386 else: return s
2387
2388 @asaplog_post_dec
2389 def regrid_channel(self, width=5, plot=False, insitu=None):
2390 """\
2391 Regrid the spectra by the specified channel width
2392
2393 Parameters:
2394
2395 width: The channel width (float) of regridded spectra
2396 in the current spectral unit.
2397
2398 plot: [NOT IMPLEMENTED YET]
2399 plot the original and the regridded spectra.
2400 In this each indivual fit has to be approved, by
2401 typing 'y' or 'n'
2402
2403 insitu: if False a new scantable is returned.
2404 Otherwise, the scaling is done in-situ
2405 The default is taken from .asaprc (False)
2406
2407 """
2408 if insitu is None: insitu = rcParams['insitu']
2409 varlist = vars()
2410
2411 if plot:
2412 asaplog.post()
2413 asaplog.push("Verification plot is not implemtnetd yet.")
2414 asaplog.post("WARN")
2415
2416 s = self.copy()
2417 s._regrid_specchan(width)
2418
2419 s._add_history("regrid_channel", varlist)
2420
2421# if plot:
2422# from asap.asapplotter import new_asaplot
2423# theplot = new_asaplot(rcParams['plotter.gui'])
2424# from matplotlib import rc as rcp
2425# rcp('lines', linewidth=1)
2426# theplot.set_panels()
2427# ylab=s._get_ordinate_label()
2428# #theplot.palette(0,["#777777","red"])
2429# for r in xrange(s.nrow()):
2430# xsm=s._getabcissa(r)
2431# ysm=s._getspectrum(r)
2432# xorg=orgscan._getabcissa(r)
2433# yorg=orgscan._getspectrum(r)
2434# theplot.clear()
2435# theplot.hold()
2436# theplot.set_axes('ylabel',ylab)
2437# theplot.set_axes('xlabel',s._getabcissalabel(r))
2438# theplot.set_axes('title',s._getsourcename(r))
2439# theplot.set_line(label='Original',color="#777777")
2440# theplot.plot(xorg,yorg)
2441# theplot.set_line(label='Smoothed',color="red")
2442# theplot.plot(xsm,ysm)
2443# ### Ugly part for legend
2444# for i in [0,1]:
2445# theplot.subplots[0]['lines'].append(
2446# [theplot.subplots[0]['axes'].lines[i]]
2447# )
2448# theplot.release()
2449# ### Ugly part for legend
2450# theplot.subplots[0]['lines']=[]
2451# res = raw_input("Accept smoothing ([y]/n): ")
2452# if res.upper() == 'N':
2453# s._setspectrum(yorg, r)
2454# theplot.quit()
2455# del theplot
2456# del orgscan
2457
2458 if insitu: self._assign(s)
2459 else: return s
2460
2461 @asaplog_post_dec
2462 def _parse_wn(self, wn):
2463 if isinstance(wn, list) or isinstance(wn, tuple):
2464 return wn
2465 elif isinstance(wn, int):
2466 return [ wn ]
2467 elif isinstance(wn, str):
2468 if '-' in wn: # case 'a-b' : return [a,a+1,...,b-1,b]
2469 val = wn.split('-')
2470 val = [int(val[0]), int(val[1])]
2471 val.sort()
2472 res = [i for i in xrange(val[0], val[1]+1)]
2473 elif wn[:2] == '<=' or wn[:2] == '=<': # cases '<=a','=<a' : return [0,1,...,a-1,a]
2474 val = int(wn[2:])+1
2475 res = [i for i in xrange(val)]
2476 elif wn[-2:] == '>=' or wn[-2:] == '=>': # cases 'a>=','a=>' : return [0,1,...,a-1,a]
2477 val = int(wn[:-2])+1
2478 res = [i for i in xrange(val)]
2479 elif wn[0] == '<': # case '<a' : return [0,1,...,a-2,a-1]
2480 val = int(wn[1:])
2481 res = [i for i in xrange(val)]
2482 elif wn[-1] == '>': # case 'a>' : return [0,1,...,a-2,a-1]
2483 val = int(wn[:-1])
2484 res = [i for i in xrange(val)]
2485 elif wn[:2] == '>=' or wn[:2] == '=>': # cases '>=a','=>a' : return [a,-999], which is
2486 # then interpreted in C++
2487 # side as [a,a+1,...,a_nyq]
2488 # (CAS-3759)
2489 val = int(wn[2:])
2490 res = [val, -999]
2491 #res = [i for i in xrange(val, self.nchan()/2+1)]
2492 elif wn[-2:] == '<=' or wn[-2:] == '=<': # cases 'a<=','a=<' : return [a,-999], which is
2493 # then interpreted in C++
2494 # side as [a,a+1,...,a_nyq]
2495 # (CAS-3759)
2496 val = int(wn[:-2])
2497 res = [val, -999]
2498 #res = [i for i in xrange(val, self.nchan()/2+1)]
2499 elif wn[0] == '>': # case '>a' : return [a+1,-999], which is
2500 # then interpreted in C++
2501 # side as [a+1,a+2,...,a_nyq]
2502 # (CAS-3759)
2503 val = int(wn[1:])+1
2504 res = [val, -999]
2505 #res = [i for i in xrange(val, self.nchan()/2+1)]
2506 elif wn[-1] == '<': # case 'a<' : return [a+1,-999], which is
2507 # then interpreted in C++
2508 # side as [a+1,a+2,...,a_nyq]
2509 # (CAS-3759)
2510 val = int(wn[:-1])+1
2511 res = [val, -999]
2512 #res = [i for i in xrange(val, self.nchan()/2+1)]
2513
2514 return res
2515 else:
2516 msg = 'wrong value given for addwn/rejwn'
2517 raise RuntimeError(msg)
2518
2519
2520 @asaplog_post_dec
2521 def sinusoid_baseline(self, insitu=None, mask=None, applyfft=None,
2522 fftmethod=None, fftthresh=None,
2523 addwn=None, rejwn=None, clipthresh=None,
2524 clipniter=None, plot=None,
2525 getresidual=None, showprogress=None,
2526 minnrow=None, outlog=None, blfile=None):
2527 """\
2528 Return a scan which has been baselined (all rows) with sinusoidal
2529 functions.
2530
2531 Parameters:
2532 insitu: if False a new scantable is returned.
2533 Otherwise, the scaling is done in-situ
2534 The default is taken from .asaprc (False)
2535 mask: an optional mask
2536 applyfft: if True use some method, such as FFT, to find
2537 strongest sinusoidal components in the wavenumber
2538 domain to be used for baseline fitting.
2539 default is True.
2540 fftmethod: method to find the strong sinusoidal components.
2541 now only 'fft' is available and it is the default.
2542 fftthresh: the threshold to select wave numbers to be used for
2543 fitting from the distribution of amplitudes in the
2544 wavenumber domain.
2545 both float and string values accepted.
2546 given a float value, the unit is set to sigma.
2547 for string values, allowed formats include:
2548 'xsigma' or 'x' (= x-sigma level. e.g.,
2549 '3sigma'), or
2550 'topx' (= the x strongest ones, e.g. 'top5').
2551 default is 3.0 (unit: sigma).
2552 addwn: the additional wave numbers to be used for fitting.
2553 list or integer value is accepted to specify every
2554 wave numbers. also string value can be used in case
2555 you need to specify wave numbers in a certain range,
2556 e.g., 'a-b' (= a, a+1, a+2, ..., b-1, b),
2557 '<a' (= 0,1,...,a-2,a-1),
2558 '>=a' (= a, a+1, ... up to the maximum wave
2559 number corresponding to the Nyquist
2560 frequency for the case of FFT).
2561 default is [0].
2562 rejwn: the wave numbers NOT to be used for fitting.
2563 can be set just as addwn but has higher priority:
2564 wave numbers which are specified both in addwn
2565 and rejwn will NOT be used. default is [].
2566 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2567 clipniter: maximum number of iteration of 'clipthresh'-sigma
2568 clipping (default is 0)
2569 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2570 plot the fit and the residual. In this each
2571 indivual fit has to be approved, by typing 'y'
2572 or 'n'
2573 getresidual: if False, returns best-fit values instead of
2574 residual. (default is True)
2575 showprogress: show progress status for large data.
2576 default is True.
2577 minnrow: minimum number of input spectra to show.
2578 default is 1000.
2579 outlog: Output the coefficients of the best-fit
2580 function to logger (default is False)
2581 blfile: Name of a text file in which the best-fit
2582 parameter values to be written
2583 (default is '': no file/logger output)
2584
2585 Example:
2586 # return a scan baselined by a combination of sinusoidal curves
2587 # having wave numbers in spectral window up to 10,
2588 # also with 3-sigma clipping, iteration up to 4 times
2589 bscan = scan.sinusoid_baseline(addwn='<=10',clipthresh=3.0,clipniter=4)
2590
2591 Note:
2592 The best-fit parameter values output in logger and/or blfile are now
2593 based on specunit of 'channel'.
2594 """
2595
2596 try:
2597 varlist = vars()
2598
2599 if insitu is None: insitu = rcParams['insitu']
2600 if insitu:
2601 workscan = self
2602 else:
2603 workscan = self.copy()
2604
2605 #if mask is None: mask = [True for i in xrange(workscan.nchan())]
2606 if mask is None: mask = []
2607 if applyfft is None: applyfft = True
2608 if fftmethod is None: fftmethod = 'fft'
2609 if fftthresh is None: fftthresh = 3.0
2610 if addwn is None: addwn = [0]
2611 if rejwn is None: rejwn = []
2612 if clipthresh is None: clipthresh = 3.0
2613 if clipniter is None: clipniter = 0
2614 if plot is None: plot = False
2615 if getresidual is None: getresidual = True
2616 if showprogress is None: showprogress = True
2617 if minnrow is None: minnrow = 1000
2618 if outlog is None: outlog = False
2619 if blfile is None: blfile = ''
2620
2621 #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
2622 workscan._sinusoid_baseline(mask, applyfft, fftmethod.lower(),
2623 str(fftthresh).lower(),
2624 workscan._parse_wn(addwn),
2625 workscan._parse_wn(rejwn), clipthresh,
2626 clipniter, getresidual,
2627 pack_progress_params(showprogress,
2628 minnrow), outlog,
2629 blfile)
2630 workscan._add_history('sinusoid_baseline', varlist)
2631
2632 if insitu:
2633 self._assign(workscan)
2634 else:
2635 return workscan
2636
2637 except RuntimeError, e:
2638 raise_fitting_failure_exception(e)
2639
2640
2641 @asaplog_post_dec
2642 def auto_sinusoid_baseline(self, insitu=None, mask=None, applyfft=None,
2643 fftmethod=None, fftthresh=None,
2644 addwn=None, rejwn=None, clipthresh=None,
2645 clipniter=None, edge=None, threshold=None,
2646 chan_avg_limit=None, plot=None,
2647 getresidual=None, showprogress=None,
2648 minnrow=None, outlog=None, blfile=None):
2649 """\
2650 Return a scan which has been baselined (all rows) with sinusoidal
2651 functions.
2652 Spectral lines are detected first using linefinder and masked out
2653 to avoid them affecting the baseline solution.
2654
2655 Parameters:
2656 insitu: if False a new scantable is returned.
2657 Otherwise, the scaling is done in-situ
2658 The default is taken from .asaprc (False)
2659 mask: an optional mask retreived from scantable
2660 applyfft: if True use some method, such as FFT, to find
2661 strongest sinusoidal components in the wavenumber
2662 domain to be used for baseline fitting.
2663 default is True.
2664 fftmethod: method to find the strong sinusoidal components.
2665 now only 'fft' is available and it is the default.
2666 fftthresh: the threshold to select wave numbers to be used for
2667 fitting from the distribution of amplitudes in the
2668 wavenumber domain.
2669 both float and string values accepted.
2670 given a float value, the unit is set to sigma.
2671 for string values, allowed formats include:
2672 'xsigma' or 'x' (= x-sigma level. e.g.,
2673 '3sigma'), or
2674 'topx' (= the x strongest ones, e.g. 'top5').
2675 default is 3.0 (unit: sigma).
2676 addwn: the additional wave numbers to be used for fitting.
2677 list or integer value is accepted to specify every
2678 wave numbers. also string value can be used in case
2679 you need to specify wave numbers in a certain range,
2680 e.g., 'a-b' (= a, a+1, a+2, ..., b-1, b),
2681 '<a' (= 0,1,...,a-2,a-1),
2682 '>=a' (= a, a+1, ... up to the maximum wave
2683 number corresponding to the Nyquist
2684 frequency for the case of FFT).
2685 default is [0].
2686 rejwn: the wave numbers NOT to be used for fitting.
2687 can be set just as addwn but has higher priority:
2688 wave numbers which are specified both in addwn
2689 and rejwn will NOT be used. default is [].
2690 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2691 clipniter: maximum number of iteration of 'clipthresh'-sigma
2692 clipping (default is 0)
2693 edge: an optional number of channel to drop at
2694 the edge of spectrum. If only one value is
2695 specified, the same number will be dropped
2696 from both sides of the spectrum. Default
2697 is to keep all channels. Nested tuples
2698 represent individual edge selection for
2699 different IFs (a number of spectral channels
2700 can be different)
2701 threshold: the threshold used by line finder. It is
2702 better to keep it large as only strong lines
2703 affect the baseline solution.
2704 chan_avg_limit: a maximum number of consequtive spectral
2705 channels to average during the search of
2706 weak and broad lines. The default is no
2707 averaging (and no search for weak lines).
2708 If such lines can affect the fitted baseline
2709 (e.g. a high order polynomial is fitted),
2710 increase this parameter (usually values up
2711 to 8 are reasonable). Most users of this
2712 method should find the default value sufficient.
2713 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2714 plot the fit and the residual. In this each
2715 indivual fit has to be approved, by typing 'y'
2716 or 'n'
2717 getresidual: if False, returns best-fit values instead of
2718 residual. (default is True)
2719 showprogress: show progress status for large data.
2720 default is True.
2721 minnrow: minimum number of input spectra to show.
2722 default is 1000.
2723 outlog: Output the coefficients of the best-fit
2724 function to logger (default is False)
2725 blfile: Name of a text file in which the best-fit
2726 parameter values to be written
2727 (default is "": no file/logger output)
2728
2729 Example:
2730 bscan = scan.auto_sinusoid_baseline(addwn='<=10', insitu=False)
2731
2732 Note:
2733 The best-fit parameter values output in logger and/or blfile are now
2734 based on specunit of 'channel'.
2735 """
2736
2737 try:
2738 varlist = vars()
2739
2740 if insitu is None: insitu = rcParams['insitu']
2741 if insitu:
2742 workscan = self
2743 else:
2744 workscan = self.copy()
2745
2746 #if mask is None: mask = [True for i in xrange(workscan.nchan())]
2747 if mask is None: mask = []
2748 if applyfft is None: applyfft = True
2749 if fftmethod is None: fftmethod = 'fft'
2750 if fftthresh is None: fftthresh = 3.0
2751 if addwn is None: addwn = [0]
2752 if rejwn is None: rejwn = []
2753 if clipthresh is None: clipthresh = 3.0
2754 if clipniter is None: clipniter = 0
2755 if edge is None: edge = (0,0)
2756 if threshold is None: threshold = 3
2757 if chan_avg_limit is None: chan_avg_limit = 1
2758 if plot is None: plot = False
2759 if getresidual is None: getresidual = True
2760 if showprogress is None: showprogress = True
2761 if minnrow is None: minnrow = 1000
2762 if outlog is None: outlog = False
2763 if blfile is None: blfile = ''
2764
2765 #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
2766 workscan._auto_sinusoid_baseline(mask, applyfft,
2767 fftmethod.lower(),
2768 str(fftthresh).lower(),
2769 workscan._parse_wn(addwn),
2770 workscan._parse_wn(rejwn),
2771 clipthresh, clipniter,
2772 normalise_edge_param(edge),
2773 threshold, chan_avg_limit,
2774 getresidual,
2775 pack_progress_params(showprogress,
2776 minnrow),
2777 outlog, blfile)
2778 workscan._add_history("auto_sinusoid_baseline", varlist)
2779
2780 if insitu:
2781 self._assign(workscan)
2782 else:
2783 return workscan
2784
2785 except RuntimeError, e:
2786 raise_fitting_failure_exception(e)
2787
2788 @asaplog_post_dec
2789 def cspline_baseline(self, insitu=None, mask=None, npiece=None,
2790 clipthresh=None, clipniter=None, plot=None,
2791 getresidual=None, showprogress=None, minnrow=None,
2792 outlog=None, blfile=None):
2793 """\
2794 Return a scan which has been baselined (all rows) by cubic spline
2795 function (piecewise cubic polynomial).
2796
2797 Parameters:
2798 insitu: If False a new scantable is returned.
2799 Otherwise, the scaling is done in-situ
2800 The default is taken from .asaprc (False)
2801 mask: An optional mask
2802 npiece: Number of pieces. (default is 2)
2803 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2804 clipniter: maximum number of iteration of 'clipthresh'-sigma
2805 clipping (default is 0)
2806 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2807 plot the fit and the residual. In this each
2808 indivual fit has to be approved, by typing 'y'
2809 or 'n'
2810 getresidual: if False, returns best-fit values instead of
2811 residual. (default is True)
2812 showprogress: show progress status for large data.
2813 default is True.
2814 minnrow: minimum number of input spectra to show.
2815 default is 1000.
2816 outlog: Output the coefficients of the best-fit
2817 function to logger (default is False)
2818 blfile: Name of a text file in which the best-fit
2819 parameter values to be written
2820 (default is "": no file/logger output)
2821
2822 Example:
2823 # return a scan baselined by a cubic spline consisting of 2 pieces
2824 # (i.e., 1 internal knot),
2825 # also with 3-sigma clipping, iteration up to 4 times
2826 bscan = scan.cspline_baseline(npiece=2,clipthresh=3.0,clipniter=4)
2827
2828 Note:
2829 The best-fit parameter values output in logger and/or blfile are now
2830 based on specunit of 'channel'.
2831 """
2832
2833 try:
2834 varlist = vars()
2835
2836 if insitu is None: insitu = rcParams['insitu']
2837 if insitu:
2838 workscan = self
2839 else:
2840 workscan = self.copy()
2841
2842 #if mask is None: mask = [True for i in xrange(workscan.nchan())]
2843 if mask is None: mask = []
2844 if npiece is None: npiece = 2
2845 if clipthresh is None: clipthresh = 3.0
2846 if clipniter is None: clipniter = 0
2847 if plot is None: plot = False
2848 if getresidual is None: getresidual = True
2849 if showprogress is None: showprogress = True
2850 if minnrow is None: minnrow = 1000
2851 if outlog is None: outlog = False
2852 if blfile is None: blfile = ''
2853
2854 #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
2855 workscan._cspline_baseline(mask, npiece, clipthresh, clipniter,
2856 getresidual,
2857 pack_progress_params(showprogress,
2858 minnrow), outlog,
2859 blfile)
2860 workscan._add_history("cspline_baseline", varlist)
2861
2862 if insitu:
2863 self._assign(workscan)
2864 else:
2865 return workscan
2866
2867 except RuntimeError, e:
2868 raise_fitting_failure_exception(e)
2869
2870 @asaplog_post_dec
2871 def auto_cspline_baseline(self, insitu=None, mask=None, npiece=None,
2872 clipthresh=None, clipniter=None,
2873 edge=None, threshold=None, chan_avg_limit=None,
2874 getresidual=None, plot=None,
2875 showprogress=None, minnrow=None, outlog=None,
2876 blfile=None):
2877 """\
2878 Return a scan which has been baselined (all rows) by cubic spline
2879 function (piecewise cubic polynomial).
2880 Spectral lines are detected first using linefinder and masked out
2881 to avoid them affecting the baseline solution.
2882
2883 Parameters:
2884 insitu: if False a new scantable is returned.
2885 Otherwise, the scaling is done in-situ
2886 The default is taken from .asaprc (False)
2887 mask: an optional mask retreived from scantable
2888 npiece: Number of pieces. (default is 2)
2889 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2890 clipniter: maximum number of iteration of 'clipthresh'-sigma
2891 clipping (default is 0)
2892 edge: an optional number of channel to drop at
2893 the edge of spectrum. If only one value is
2894 specified, the same number will be dropped
2895 from both sides of the spectrum. Default
2896 is to keep all channels. Nested tuples
2897 represent individual edge selection for
2898 different IFs (a number of spectral channels
2899 can be different)
2900 threshold: the threshold used by line finder. It is
2901 better to keep it large as only strong lines
2902 affect the baseline solution.
2903 chan_avg_limit: a maximum number of consequtive spectral
2904 channels to average during the search of
2905 weak and broad lines. The default is no
2906 averaging (and no search for weak lines).
2907 If such lines can affect the fitted baseline
2908 (e.g. a high order polynomial is fitted),
2909 increase this parameter (usually values up
2910 to 8 are reasonable). Most users of this
2911 method should find the default value sufficient.
2912 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2913 plot the fit and the residual. In this each
2914 indivual fit has to be approved, by typing 'y'
2915 or 'n'
2916 getresidual: if False, returns best-fit values instead of
2917 residual. (default is True)
2918 showprogress: show progress status for large data.
2919 default is True.
2920 minnrow: minimum number of input spectra to show.
2921 default is 1000.
2922 outlog: Output the coefficients of the best-fit
2923 function to logger (default is False)
2924 blfile: Name of a text file in which the best-fit
2925 parameter values to be written
2926 (default is "": no file/logger output)
2927
2928 Example:
2929 bscan = scan.auto_cspline_baseline(npiece=3, insitu=False)
2930
2931 Note:
2932 The best-fit parameter values output in logger and/or blfile are now
2933 based on specunit of 'channel'.
2934 """
2935
2936 try:
2937 varlist = vars()
2938
2939 if insitu is None: insitu = rcParams['insitu']
2940 if insitu:
2941 workscan = self
2942 else:
2943 workscan = self.copy()
2944
2945 #if mask is None: mask = [True for i in xrange(workscan.nchan())]
2946 if mask is None: mask = []
2947 if npiece is None: npiece = 2
2948 if clipthresh is None: clipthresh = 3.0
2949 if clipniter is None: clipniter = 0
2950 if edge is None: edge = (0, 0)
2951 if threshold is None: threshold = 3
2952 if chan_avg_limit is None: chan_avg_limit = 1
2953 if plot is None: plot = False
2954 if getresidual is None: getresidual = True
2955 if showprogress is None: showprogress = True
2956 if minnrow is None: minnrow = 1000
2957 if outlog is None: outlog = False
2958 if blfile is None: blfile = ''
2959
2960 #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
2961 workscan._auto_cspline_baseline(mask, npiece, clipthresh,
2962 clipniter,
2963 normalise_edge_param(edge),
2964 threshold,
2965 chan_avg_limit, getresidual,
2966 pack_progress_params(showprogress,
2967 minnrow),
2968 outlog, blfile)
2969 workscan._add_history("auto_cspline_baseline", varlist)
2970
2971 if insitu:
2972 self._assign(workscan)
2973 else:
2974 return workscan
2975
2976 except RuntimeError, e:
2977 raise_fitting_failure_exception(e)
2978
2979 @asaplog_post_dec
2980 def poly_baseline(self, mask=None, order=None, insitu=None, plot=None,
2981 getresidual=None, showprogress=None, minnrow=None,
2982 outlog=None, blfile=None):
2983 """\
2984 Return a scan which has been baselined (all rows) by a polynomial.
2985 Parameters:
2986 insitu: if False a new scantable is returned.
2987 Otherwise, the scaling is done in-situ
2988 The default is taken from .asaprc (False)
2989 mask: an optional mask
2990 order: the order of the polynomial (default is 0)
2991 plot: plot the fit and the residual. In this each
2992 indivual fit has to be approved, by typing 'y'
2993 or 'n'
2994 getresidual: if False, returns best-fit values instead of
2995 residual. (default is True)
2996 showprogress: show progress status for large data.
2997 default is True.
2998 minnrow: minimum number of input spectra to show.
2999 default is 1000.
3000 outlog: Output the coefficients of the best-fit
3001 function to logger (default is False)
3002 blfile: Name of a text file in which the best-fit
3003 parameter values to be written
3004 (default is "": no file/logger output)
3005
3006 Example:
3007 # return a scan baselined by a third order polynomial,
3008 # not using a mask
3009 bscan = scan.poly_baseline(order=3)
3010 """
3011
3012 try:
3013 varlist = vars()
3014
3015 if insitu is None:
3016 insitu = rcParams["insitu"]
3017 if insitu:
3018 workscan = self
3019 else:
3020 workscan = self.copy()
3021
3022 #if mask is None: mask = [True for i in \
3023 # xrange(workscan.nchan())]
3024 if mask is None: mask = []
3025 if order is None: order = 0
3026 if plot is None: plot = False
3027 if getresidual is None: getresidual = True
3028 if showprogress is None: showprogress = True
3029 if minnrow is None: minnrow = 1000
3030 if outlog is None: outlog = False
3031 if blfile is None: blfile = ""
3032
3033 if plot:
3034 outblfile = (blfile != "") and \
3035 os.path.exists(os.path.expanduser(
3036 os.path.expandvars(blfile))
3037 )
3038 if outblfile:
3039 blf = open(blfile, "a")
3040
3041 f = fitter()
3042 f.set_function(lpoly=order)
3043
3044 rows = xrange(workscan.nrow())
3045 #if len(rows) > 0: workscan._init_blinfo()
3046
3047 action = "H"
3048 for r in rows:
3049 f.x = workscan._getabcissa(r)
3050 f.y = workscan._getspectrum(r)
3051 if mask:
3052 f.mask = mask_and(mask, workscan._getmask(r)) # (CAS-1434)
3053 else: # mask=None
3054 f.mask = workscan._getmask(r)
3055
3056 f.data = None
3057 f.fit()
3058
3059 if action != "Y": # skip plotting when accepting all
3060 f.plot(residual=True)
3061 #accept_fit = raw_input("Accept fit ( [y]/n ): ")
3062 #if accept_fit.upper() == "N":
3063 # #workscan._append_blinfo(None, None, None)
3064 # continue
3065 accept_fit = self._get_verify_action("Accept fit?",action)
3066 if r == 0: action = None
3067 if accept_fit.upper() == "N":
3068 continue
3069 elif accept_fit.upper() == "R":
3070 break
3071 elif accept_fit.upper() == "A":
3072 action = "Y"
3073
3074 blpars = f.get_parameters()
3075 masklist = workscan.get_masklist(f.mask, row=r, silent=True)
3076 #workscan._append_blinfo(blpars, masklist, f.mask)
3077 workscan._setspectrum((f.fitter.getresidual()
3078 if getresidual else f.fitter.getfit()), r)
3079
3080 if outblfile:
3081 rms = workscan.get_rms(f.mask, r)
3082 dataout = \
3083 workscan.format_blparams_row(blpars["params"],
3084 blpars["fixed"],
3085 rms, str(masklist),
3086 r, True)
3087 blf.write(dataout)
3088
3089 f._p.unmap()
3090 f._p = None
3091
3092 if outblfile:
3093 blf.close()
3094 else:
3095 workscan._poly_baseline(mask, order, getresidual,
3096 pack_progress_params(showprogress,
3097 minnrow),
3098 outlog, blfile)
3099
3100 workscan._add_history("poly_baseline", varlist)
3101
3102 if insitu:
3103 self._assign(workscan)
3104 else:
3105 return workscan
3106
3107 except RuntimeError, e:
3108 raise_fitting_failure_exception(e)
3109
3110 @asaplog_post_dec
3111 def auto_poly_baseline(self, mask=None, order=None, edge=None,
3112 threshold=None, chan_avg_limit=None,
3113 plot=None, insitu=None,
3114 getresidual=None, showprogress=None,
3115 minnrow=None, outlog=None, blfile=None):
3116 """\
3117 Return a scan which has been baselined (all rows) by a polynomial.
3118 Spectral lines are detected first using linefinder and masked out
3119 to avoid them affecting the baseline solution.
3120
3121 Parameters:
3122 insitu: if False a new scantable is returned.
3123 Otherwise, the scaling is done in-situ
3124 The default is taken from .asaprc (False)
3125 mask: an optional mask retreived from scantable
3126 order: the order of the polynomial (default is 0)
3127 edge: an optional number of channel to drop at
3128 the edge of spectrum. If only one value is
3129 specified, the same number will be dropped
3130 from both sides of the spectrum. Default
3131 is to keep all channels. Nested tuples
3132 represent individual edge selection for
3133 different IFs (a number of spectral channels
3134 can be different)
3135 threshold: the threshold used by line finder. It is
3136 better to keep it large as only strong lines
3137 affect the baseline solution.
3138 chan_avg_limit: a maximum number of consequtive spectral
3139 channels to average during the search of
3140 weak and broad lines. The default is no
3141 averaging (and no search for weak lines).
3142 If such lines can affect the fitted baseline
3143 (e.g. a high order polynomial is fitted),
3144 increase this parameter (usually values up
3145 to 8 are reasonable). Most users of this
3146 method should find the default value sufficient.
3147 plot: plot the fit and the residual. In this each
3148 indivual fit has to be approved, by typing 'y'
3149 or 'n'
3150 getresidual: if False, returns best-fit values instead of
3151 residual. (default is True)
3152 showprogress: show progress status for large data.
3153 default is True.
3154 minnrow: minimum number of input spectra to show.
3155 default is 1000.
3156 outlog: Output the coefficients of the best-fit
3157 function to logger (default is False)
3158 blfile: Name of a text file in which the best-fit
3159 parameter values to be written
3160 (default is "": no file/logger output)
3161
3162 Example:
3163 bscan = scan.auto_poly_baseline(order=7, insitu=False)
3164 """
3165
3166 try:
3167 varlist = vars()
3168
3169 if insitu is None:
3170 insitu = rcParams['insitu']
3171 if insitu:
3172 workscan = self
3173 else:
3174 workscan = self.copy()
3175
3176 #if mask is None: mask = [True for i in xrange(workscan.nchan())]
3177 if mask is None: mask = []
3178 if order is None: order = 0
3179 if edge is None: edge = (0, 0)
3180 if threshold is None: threshold = 3
3181 if chan_avg_limit is None: chan_avg_limit = 1
3182 if plot is None: plot = False
3183 if getresidual is None: getresidual = True
3184 if showprogress is None: showprogress = True
3185 if minnrow is None: minnrow = 1000
3186 if outlog is None: outlog = False
3187 if blfile is None: blfile = ''
3188
3189 edge = normalise_edge_param(edge)
3190
3191 if plot:
3192 outblfile = (blfile != "") and \
3193 os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
3194 if outblfile: blf = open(blfile, "a")
3195
3196 from asap.asaplinefind import linefinder
3197 fl = linefinder()
3198 fl.set_options(threshold=threshold, avg_limit=chan_avg_limit)
3199 fl.set_scan(workscan)
3200
3201 f = fitter()
3202 f.set_function(lpoly=order)
3203
3204 rows = xrange(workscan.nrow())
3205 #if len(rows) > 0: workscan._init_blinfo()
3206
3207 action = "H"
3208 for r in rows:
3209 idx = 2*workscan.getif(r)
3210 if mask:
3211 msk = mask_and(mask, workscan._getmask(r)) # (CAS-1434)
3212 else: # mask=None
3213 msk = workscan._getmask(r)
3214 fl.find_lines(r, msk, edge[idx:idx+2])
3215
3216 f.x = workscan._getabcissa(r)
3217 f.y = workscan._getspectrum(r)
3218 f.mask = fl.get_mask()
3219 f.data = None
3220 f.fit()
3221
3222 if action != "Y": # skip plotting when accepting all
3223 f.plot(residual=True)
3224 #accept_fit = raw_input("Accept fit ( [y]/n ): ")
3225 accept_fit = self._get_verify_action("Accept fit?",action)
3226 if r == 0: action = None
3227 if accept_fit.upper() == "N":
3228 #workscan._append_blinfo(None, None, None)
3229 continue
3230 elif accept_fit.upper() == "R":
3231 break
3232 elif accept_fit.upper() == "A":
3233 action = "Y"
3234
3235 blpars = f.get_parameters()
3236 masklist = workscan.get_masklist(f.mask, row=r, silent=True)
3237 #workscan._append_blinfo(blpars, masklist, f.mask)
3238 workscan._setspectrum(
3239 (f.fitter.getresidual() if getresidual
3240 else f.fitter.getfit()), r
3241 )
3242
3243 if outblfile:
3244 rms = workscan.get_rms(f.mask, r)
3245 dataout = \
3246 workscan.format_blparams_row(blpars["params"],
3247 blpars["fixed"],
3248 rms, str(masklist),
3249 r, True)
3250 blf.write(dataout)
3251
3252 f._p.unmap()
3253 f._p = None
3254
3255 if outblfile: blf.close()
3256 else:
3257 workscan._auto_poly_baseline(mask, order, edge, threshold,
3258 chan_avg_limit, getresidual,
3259 pack_progress_params(showprogress,
3260 minnrow),
3261 outlog, blfile)
3262
3263 workscan._add_history("auto_poly_baseline", varlist)
3264
3265 if insitu:
3266 self._assign(workscan)
3267 else:
3268 return workscan
3269
3270 except RuntimeError, e:
3271 raise_fitting_failure_exception(e)
3272
3273 def _init_blinfo(self):
3274 """\
3275 Initialise the following three auxiliary members:
3276 blpars : parameters of the best-fit baseline,
3277 masklists : mask data (edge positions of masked channels) and
3278 actualmask : mask data (in boolean list),
3279 to keep for use later (including output to logger/text files).
3280 Used by poly_baseline() and auto_poly_baseline() in case of
3281 'plot=True'.
3282 """
3283 self.blpars = []
3284 self.masklists = []
3285 self.actualmask = []
3286 return
3287
3288 def _append_blinfo(self, data_blpars, data_masklists, data_actualmask):
3289 """\
3290 Append baseline-fitting related info to blpars, masklist and
3291 actualmask.
3292 """
3293 self.blpars.append(data_blpars)
3294 self.masklists.append(data_masklists)
3295 self.actualmask.append(data_actualmask)
3296 return
3297
3298 @asaplog_post_dec
3299 def rotate_linpolphase(self, angle):
3300 """\
3301 Rotate the phase of the complex polarization O=Q+iU correlation.
3302 This is always done in situ in the raw data. So if you call this
3303 function more than once then each call rotates the phase further.
3304
3305 Parameters:
3306
3307 angle: The angle (degrees) to rotate (add) by.
3308
3309 Example::
3310
3311 scan.rotate_linpolphase(2.3)
3312
3313 """
3314 varlist = vars()
3315 self._math._rotate_linpolphase(self, angle)
3316 self._add_history("rotate_linpolphase", varlist)
3317 return
3318
3319 @asaplog_post_dec
3320 def rotate_xyphase(self, angle):
3321 """\
3322 Rotate the phase of the XY correlation. This is always done in situ
3323 in the data. So if you call this function more than once
3324 then each call rotates the phase further.
3325
3326 Parameters:
3327
3328 angle: The angle (degrees) to rotate (add) by.
3329
3330 Example::
3331
3332 scan.rotate_xyphase(2.3)
3333
3334 """
3335 varlist = vars()
3336 self._math._rotate_xyphase(self, angle)
3337 self._add_history("rotate_xyphase", varlist)
3338 return
3339
3340 @asaplog_post_dec
3341 def swap_linears(self):
3342 """\
3343 Swap the linear polarisations XX and YY, or better the first two
3344 polarisations as this also works for ciculars.
3345 """
3346 varlist = vars()
3347 self._math._swap_linears(self)
3348 self._add_history("swap_linears", varlist)
3349 return
3350
3351 @asaplog_post_dec
3352 def invert_phase(self):
3353 """\
3354 Invert the phase of the complex polarisation
3355 """
3356 varlist = vars()
3357 self._math._invert_phase(self)
3358 self._add_history("invert_phase", varlist)
3359 return
3360
3361 @asaplog_post_dec
3362 def add(self, offset, insitu=None):
3363 """\
3364 Return a scan where all spectra have the offset added
3365
3366 Parameters:
3367
3368 offset: the offset
3369
3370 insitu: if False a new scantable is returned.
3371 Otherwise, the scaling is done in-situ
3372 The default is taken from .asaprc (False)
3373
3374 """
3375 if insitu is None: insitu = rcParams['insitu']
3376 self._math._setinsitu(insitu)
3377 varlist = vars()
3378 s = scantable(self._math._unaryop(self, offset, "ADD", False))
3379 s._add_history("add", varlist)
3380 if insitu:
3381 self._assign(s)
3382 else:
3383 return s
3384
3385 @asaplog_post_dec
3386 def scale(self, factor, tsys=True, insitu=None):
3387 """\
3388
3389 Return a scan where all spectra are scaled by the given 'factor'
3390
3391 Parameters:
3392
3393 factor: the scaling factor (float or 1D float list)
3394
3395 insitu: if False a new scantable is returned.
3396 Otherwise, the scaling is done in-situ
3397 The default is taken from .asaprc (False)
3398
3399 tsys: if True (default) then apply the operation to Tsys
3400 as well as the data
3401
3402 """
3403 if insitu is None: insitu = rcParams['insitu']
3404 self._math._setinsitu(insitu)
3405 varlist = vars()
3406 s = None
3407 import numpy
3408 if isinstance(factor, list) or isinstance(factor, numpy.ndarray):
3409 if isinstance(factor[0], list) or isinstance(factor[0],
3410 numpy.ndarray):
3411 from asapmath import _array2dOp
3412 s = _array2dOp( self, factor, "MUL", tsys, insitu )
3413 else:
3414 s = scantable( self._math._arrayop( self, factor,
3415 "MUL", tsys ) )
3416 else:
3417 s = scantable(self._math._unaryop(self, factor, "MUL", tsys))
3418 s._add_history("scale", varlist)
3419 if insitu:
3420 self._assign(s)
3421 else:
3422 return s
3423
3424 @preserve_selection
3425 def set_sourcetype(self, match, matchtype="pattern",
3426 sourcetype="reference"):
3427 """\
3428 Set the type of the source to be an source or reference scan
3429 using the provided pattern.
3430
3431 Parameters:
3432
3433 match: a Unix style pattern, regular expression or selector
3434
3435 matchtype: 'pattern' (default) UNIX style pattern or
3436 'regex' regular expression
3437
3438 sourcetype: the type of the source to use (source/reference)
3439
3440 """
3441 varlist = vars()
3442 stype = -1
3443 if sourcetype.lower().startswith("r") or sourcetype.lower() == "off":
3444 stype = 1
3445 elif sourcetype.lower().startswith("s") or sourcetype.lower() == "on":
3446 stype = 0
3447 else:
3448 raise ValueError("Illegal sourcetype use s(ource)/on or r(eference)/off")
3449 if matchtype.lower().startswith("p"):
3450 matchtype = "pattern"
3451 elif matchtype.lower().startswith("r"):
3452 matchtype = "regex"
3453 else:
3454 raise ValueError("Illegal matchtype, use p(attern) or r(egex)")
3455 sel = selector()
3456 if isinstance(match, selector):
3457 sel = match
3458 else:
3459 sel.set_query("SRCNAME=%s('%s')" % (matchtype, match))
3460 self.set_selection(sel)
3461 self._setsourcetype(stype)
3462 self._add_history("set_sourcetype", varlist)
3463
3464 @asaplog_post_dec
3465 @preserve_selection
3466 def auto_quotient(self, preserve=True, mode='paired', verify=False):
3467 """\
3468 This function allows to build quotients automatically.
3469 It assumes the observation to have the same number of
3470 "ons" and "offs"
3471
3472 Parameters:
3473
3474 preserve: you can preserve (default) the continuum or
3475 remove it. The equations used are
3476
3477 preserve: Output = Toff * (on/off) - Toff
3478
3479 remove: Output = Toff * (on/off) - Ton
3480
3481 mode: the on/off detection mode
3482 'paired' (default)
3483 identifies 'off' scans by the
3484 trailing '_R' (Mopra/Parkes) or
3485 '_e'/'_w' (Tid) and matches
3486 on/off pairs from the observing pattern
3487 'time'
3488 finds the closest off in time
3489
3490 .. todo:: verify argument is not implemented
3491
3492 """
3493 varlist = vars()
3494 modes = ["time", "paired"]
3495 if not mode in modes:
3496 msg = "please provide valid mode. Valid modes are %s" % (modes)
3497 raise ValueError(msg)
3498 s = None
3499 if mode.lower() == "paired":
3500 sel = self.get_selection()
3501 sel.set_query("SRCTYPE==psoff")
3502 self.set_selection(sel)
3503 offs = self.copy()
3504 sel.set_query("SRCTYPE==pson")
3505 self.set_selection(sel)
3506 ons = self.copy()
3507 s = scantable(self._math._quotient(ons, offs, preserve))
3508 elif mode.lower() == "time":
3509 s = scantable(self._math._auto_quotient(self, mode, preserve))
3510 s._add_history("auto_quotient", varlist)
3511 return s
3512
3513 @asaplog_post_dec
3514 def mx_quotient(self, mask = None, weight='median', preserve=True):
3515 """\
3516 Form a quotient using "off" beams when observing in "MX" mode.
3517
3518 Parameters:
3519
3520 mask: an optional mask to be used when weight == 'stddev'
3521
3522 weight: How to average the off beams. Default is 'median'.
3523
3524 preserve: you can preserve (default) the continuum or
3525 remove it. The equations used are:
3526
3527 preserve: Output = Toff * (on/off) - Toff
3528
3529 remove: Output = Toff * (on/off) - Ton
3530
3531 """
3532 mask = mask or ()
3533 varlist = vars()
3534 on = scantable(self._math._mx_extract(self, 'on'))
3535 preoff = scantable(self._math._mx_extract(self, 'off'))
3536 off = preoff.average_time(mask=mask, weight=weight, scanav=False)
3537 from asapmath import quotient
3538 q = quotient(on, off, preserve)
3539 q._add_history("mx_quotient", varlist)
3540 return q
3541
3542 @asaplog_post_dec
3543 def freq_switch(self, insitu=None):
3544 """\
3545 Apply frequency switching to the data.
3546
3547 Parameters:
3548
3549 insitu: if False a new scantable is returned.
3550 Otherwise, the swictching is done in-situ
3551 The default is taken from .asaprc (False)
3552
3553 """
3554 if insitu is None: insitu = rcParams['insitu']
3555 self._math._setinsitu(insitu)
3556 varlist = vars()
3557 s = scantable(self._math._freqswitch(self))
3558 s._add_history("freq_switch", varlist)
3559 if insitu:
3560 self._assign(s)
3561 else:
3562 return s
3563
3564 @asaplog_post_dec
3565 def recalc_azel(self):
3566 """Recalculate the azimuth and elevation for each position."""
3567 varlist = vars()
3568 self._recalcazel()
3569 self._add_history("recalc_azel", varlist)
3570 return
3571
3572 @asaplog_post_dec
3573 def __add__(self, other):
3574 """
3575 implicit on all axes and on Tsys
3576 """
3577 varlist = vars()
3578 s = self.__op( other, "ADD" )
3579 s._add_history("operator +", varlist)
3580 return s
3581
3582 @asaplog_post_dec
3583 def __sub__(self, other):
3584 """
3585 implicit on all axes and on Tsys
3586 """
3587 varlist = vars()
3588 s = self.__op( other, "SUB" )
3589 s._add_history("operator -", varlist)
3590 return s
3591
3592 @asaplog_post_dec
3593 def __mul__(self, other):
3594 """
3595 implicit on all axes and on Tsys
3596 """
3597 varlist = vars()
3598 s = self.__op( other, "MUL" ) ;
3599 s._add_history("operator *", varlist)
3600 return s
3601
3602
3603 @asaplog_post_dec
3604 def __div__(self, other):
3605 """
3606 implicit on all axes and on Tsys
3607 """
3608 varlist = vars()
3609 s = self.__op( other, "DIV" )
3610 s._add_history("operator /", varlist)
3611 return s
3612
3613 @asaplog_post_dec
3614 def __op( self, other, mode ):
3615 s = None
3616 if isinstance(other, scantable):
3617 s = scantable(self._math._binaryop(self, other, mode))
3618 elif isinstance(other, float):
3619 if other == 0.0:
3620 raise ZeroDivisionError("Dividing by zero is not recommended")
3621 s = scantable(self._math._unaryop(self, other, mode, False))
3622 elif isinstance(other, list) or isinstance(other, numpy.ndarray):
3623 if isinstance(other[0], list) \
3624 or isinstance(other[0], numpy.ndarray):
3625 from asapmath import _array2dOp
3626 s = _array2dOp( self, other, mode, False )
3627 else:
3628 s = scantable( self._math._arrayop( self, other,
3629 mode, False ) )
3630 else:
3631 raise TypeError("Other input is not a scantable or float value")
3632 return s
3633
3634 @asaplog_post_dec
3635 def get_fit(self, row=0):
3636 """\
3637 Print or return the stored fits for a row in the scantable
3638
3639 Parameters:
3640
3641 row: the row which the fit has been applied to.
3642
3643 """
3644 if row > self.nrow():
3645 return
3646 from asap.asapfit import asapfit
3647 fit = asapfit(self._getfit(row))
3648 asaplog.push( '%s' %(fit) )
3649 return fit.as_dict()
3650
3651 @preserve_selection
3652 def flag_nans(self):
3653 """\
3654 Utility function to flag NaN values in the scantable.
3655 """
3656 import numpy
3657 basesel = self.get_selection()
3658 for i in range(self.nrow()):
3659 sel = self.get_row_selector(i)
3660 self.set_selection(basesel+sel)
3661 nans = numpy.isnan(self._getspectrum(0))
3662 if numpy.any(nans):
3663 bnans = [ bool(v) for v in nans]
3664 self.flag(bnans)
3665
3666 def get_row_selector(self, rowno):
3667 return selector(rows=[rowno])
3668
3669 def _add_history(self, funcname, parameters):
3670 if not rcParams['scantable.history']:
3671 return
3672 # create date
3673 sep = "##"
3674 from datetime import datetime
3675 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
3676 hist = dstr+sep
3677 hist += funcname+sep#cdate+sep
3678 if parameters.has_key('self'):
3679 del parameters['self']
3680 for k, v in parameters.iteritems():
3681 if type(v) is dict:
3682 for k2, v2 in v.iteritems():
3683 hist += k2
3684 hist += "="
3685 if isinstance(v2, scantable):
3686 hist += 'scantable'
3687 elif k2 == 'mask':
3688 if isinstance(v2, list) or isinstance(v2, tuple):
3689 hist += str(self._zip_mask(v2))
3690 else:
3691 hist += str(v2)
3692 else:
3693 hist += str(v2)
3694 else:
3695 hist += k
3696 hist += "="
3697 if isinstance(v, scantable):
3698 hist += 'scantable'
3699 elif k == 'mask':
3700 if isinstance(v, list) or isinstance(v, tuple):
3701 hist += str(self._zip_mask(v))
3702 else:
3703 hist += str(v)
3704 else:
3705 hist += str(v)
3706 hist += sep
3707 hist = hist[:-2] # remove trailing '##'
3708 self._addhistory(hist)
3709
3710
3711 def _zip_mask(self, mask):
3712 mask = list(mask)
3713 i = 0
3714 segments = []
3715 while mask[i:].count(1):
3716 i += mask[i:].index(1)
3717 if mask[i:].count(0):
3718 j = i + mask[i:].index(0)
3719 else:
3720 j = len(mask)
3721 segments.append([i, j])
3722 i = j
3723 return segments
3724
3725 def _get_ordinate_label(self):
3726 fu = "("+self.get_fluxunit()+")"
3727 import re
3728 lbl = "Intensity"
3729 if re.match(".K.", fu):
3730 lbl = "Brightness Temperature "+ fu
3731 elif re.match(".Jy.", fu):
3732 lbl = "Flux density "+ fu
3733 return lbl
3734
3735 def _check_ifs(self):
3736# return len(set([self.nchan(i) for i in self.getifnos()])) == 1
3737 nchans = [self.nchan(i) for i in self.getifnos()]
3738 nchans = filter(lambda t: t > 0, nchans)
3739 return (sum(nchans)/len(nchans) == nchans[0])
3740
3741 @asaplog_post_dec
3742 def _fill(self, names, unit, average, opts={}):
3743 first = True
3744 fullnames = []
3745 for name in names:
3746 name = os.path.expandvars(name)
3747 name = os.path.expanduser(name)
3748 if not os.path.exists(name):
3749 msg = "File '%s' does not exists" % (name)
3750 raise IOError(msg)
3751 fullnames.append(name)
3752 if average:
3753 asaplog.push('Auto averaging integrations')
3754 stype = int(rcParams['scantable.storage'].lower() == 'disk')
3755 for name in fullnames:
3756 tbl = Scantable(stype)
3757 if is_ms( name ):
3758 r = msfiller( tbl )
3759 else:
3760 r = filler( tbl )
3761 msg = "Importing %s..." % (name)
3762 asaplog.push(msg, False)
3763 r.open(name, opts)
3764 rx = rcParams['scantable.reference']
3765 r.setreferenceexpr(rx)
3766 r.fill()
3767 if average:
3768 tbl = self._math._average((tbl, ), (), 'NONE', 'SCAN')
3769 if not first:
3770 tbl = self._math._merge([self, tbl])
3771 Scantable.__init__(self, tbl)
3772 r.close()
3773 del r, tbl
3774 first = False
3775 #flush log
3776 asaplog.post()
3777 if unit is not None:
3778 self.set_fluxunit(unit)
3779 if not is_casapy():
3780 self.set_freqframe(rcParams['scantable.freqframe'])
3781
3782 def _get_verify_action( self, msg, action=None ):
3783 valid_act = ['Y', 'N', 'A', 'R']
3784 if not action or not isinstance(action, str):
3785 action = raw_input("%s [Y/n/a/r] (h for help): " % msg)
3786 if action == '':
3787 return "Y"
3788 elif (action.upper()[0] in valid_act):
3789 return action.upper()[0]
3790 elif (action.upper()[0] in ['H','?']):
3791 print "Available actions of verification [Y|n|a|r]"
3792 print " Y : Yes for current data (default)"
3793 print " N : No for current data"
3794 print " A : Accept all in the following and exit from verification"
3795 print " R : Reject all in the following and exit from verification"
3796 print " H or ?: help (show this message)"
3797 return self._get_verify_action(msg)
3798 else:
3799 return 'Y'
3800
3801 def __getitem__(self, key):
3802 if key < 0:
3803 key += self.nrow()
3804 if key >= self.nrow():
3805 raise IndexError("Row index out of range.")
3806 return self._getspectrum(key)
3807
3808 def __setitem__(self, key, value):
3809 if key < 0:
3810 key += self.nrow()
3811 if key >= self.nrow():
3812 raise IndexError("Row index out of range.")
3813 if not hasattr(value, "__len__") or \
3814 len(value) > self.nchan(self.getif(key)):
3815 raise ValueError("Spectrum length doesn't match.")
3816 return self._setspectrum(value, key)
3817
3818 def __len__(self):
3819 return self.nrow()
3820
3821 def __iter__(self):
3822 for i in range(len(self)):
3823 yield self[i]
Note: See TracBrowser for help on using the repository browser.