source: trunk/python/scantable.py@ 2601

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

New Development: No

JIRA Issue: No (a bug fix)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): scantable, and sdtasks with restfreq parameter

Description:

Fixed a bug in scantable.set_restfreq(), i.e., the method ignored user selection
of scantable and caused unexpected modification of restvfrequencies.
It now modifies rest frequencies of selected data only.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 139.8 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("Invalid mask expression")
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
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 offset=1, 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 def _parse_selection(self, selstr, typestr='float', offset=0.,
1605 minval=None, maxval=None):
1606 """
1607 Parameters:
1608 selstr : The Selection string, e.g., '<3;5~7;100~103;9'
1609 typestr : The type of the values in returned list
1610 ('integer' or 'float')
1611 offset : The offset value to subtract from or add to
1612 the boundary value if the selection string
1613 includes '<' or '>'
1614 minval, maxval : The minimum/maximum values to set if the
1615 selection string includes '<' or '>'.
1616 The list element is filled with None by default.
1617 Returns:
1618 A list of min/max pair of selections.
1619 Example:
1620 _parseSelection('<3;5~7;9',typestr='int',offset=1,minval=0)
1621 returns [[0,2],[5,7],[9,9]]
1622 """
1623 selgroups = selstr.split(';')
1624 sellists = []
1625 if typestr.lower().startswith('int'):
1626 formatfunc = int
1627 else:
1628 formatfunc = float
1629
1630 for currsel in selgroups:
1631 if currsel.find('~') > 0:
1632 minsel = formatfunc(currsel.split('~')[0].strip())
1633 maxsel = formatfunc(currsel.split('~')[1].strip())
1634 elif currsel.strip().startswith('<'):
1635 minsel = minval
1636 maxsel = formatfunc(currsel.split('<')[1].strip()) \
1637 - formatfunc(offset)
1638 elif currsel.strip().startswith('>'):
1639 minsel = formatfunc(currsel.split('>')[1].strip()) \
1640 + formatfunc(offset)
1641 maxsel = maxval
1642 else:
1643 minsel = formatfunc(currsel)
1644 maxsel = formatfunc(currsel)
1645 sellists.append([minsel,maxsel])
1646 return sellists
1647
1648# def get_restfreqs(self):
1649# """
1650# Get the restfrequency(s) stored in this scantable.
1651# The return value(s) are always of unit 'Hz'
1652# Parameters:
1653# none
1654# Returns:
1655# a list of doubles
1656# """
1657# return list(self._getrestfreqs())
1658
1659 def get_restfreqs(self, ids=None):
1660 """\
1661 Get the restfrequency(s) stored in this scantable.
1662 The return value(s) are always of unit 'Hz'
1663
1664 Parameters:
1665
1666 ids: (optional) a list of MOLECULE_ID for that restfrequency(s) to
1667 be retrieved
1668
1669 Returns:
1670
1671 dictionary containing ids and a list of doubles for each id
1672
1673 """
1674 if ids is None:
1675 rfreqs = {}
1676 idlist = self.getmolnos()
1677 for i in idlist:
1678 rfreqs[i] = list(self._getrestfreqs(i))
1679 return rfreqs
1680 else:
1681 if type(ids) == list or type(ids) == tuple:
1682 rfreqs = {}
1683 for i in ids:
1684 rfreqs[i] = list(self._getrestfreqs(i))
1685 return rfreqs
1686 else:
1687 return list(self._getrestfreqs(ids))
1688
1689 @asaplog_post_dec
1690 def set_restfreqs(self, freqs=None, unit='Hz'):
1691 """\
1692 Set or replace the restfrequency specified and
1693 if the 'freqs' argument holds a scalar,
1694 then that rest frequency will be applied to all the selected
1695 data. If the 'freqs' argument holds
1696 a vector, then it MUST be of equal or smaller length than
1697 the number of IFs (and the available restfrequencies will be
1698 replaced by this vector). In this case, *all* data have
1699 the restfrequency set per IF according
1700 to the corresponding value you give in the 'freqs' vector.
1701 E.g. 'freqs=[1e9, 2e9]' would mean IF 0 gets restfreq 1e9 and
1702 IF 1 gets restfreq 2e9.
1703
1704 You can also specify the frequencies via a linecatalog.
1705
1706 Parameters:
1707
1708 freqs: list of rest frequency values or string idenitfiers
1709
1710 unit: unit for rest frequency (default 'Hz')
1711
1712
1713 Example::
1714
1715 # set the given restfrequency for the all currently selected IFs
1716 scan.set_restfreqs(freqs=1.4e9)
1717 # set restfrequencies for the n IFs (n > 1) in the order of the
1718 # list, i.e
1719 # IF0 -> 1.4e9, IF1 -> 1.41e9, IF3 -> 1.42e9
1720 # len(list_of_restfreqs) == nIF
1721 # for nIF == 1 the following will set multiple restfrequency for
1722 # that IF
1723 scan.set_restfreqs(freqs=[1.4e9, 1.41e9, 1.42e9])
1724 # set multiple restfrequencies per IF. as a list of lists where
1725 # the outer list has nIF elements, the inner s arbitrary
1726 scan.set_restfreqs(freqs=[[1.4e9, 1.41e9], [1.67e9]])
1727
1728 *Note*:
1729
1730 To do more sophisticate Restfrequency setting, e.g. on a
1731 source and IF basis, use scantable.set_selection() before using
1732 this function::
1733
1734 # provided your scantable is called scan
1735 selection = selector()
1736 selection.set_name('ORION*')
1737 selection.set_ifs([1])
1738 scan.set_selection(selection)
1739 scan.set_restfreqs(freqs=86.6e9)
1740
1741 """
1742 varlist = vars()
1743 from asap import linecatalog
1744 # simple value
1745 if isinstance(freqs, int) or isinstance(freqs, float):
1746 self._setrestfreqs([freqs], [""], unit)
1747 # list of values
1748 elif isinstance(freqs, list) or isinstance(freqs, tuple):
1749 # list values are scalars
1750 if isinstance(freqs[-1], int) or isinstance(freqs[-1], float):
1751 if len(freqs) == 1:
1752 self._setrestfreqs(freqs, [""], unit)
1753 else:
1754 # allow the 'old' mode of setting mulitple IFs
1755 savesel = self._getselection()
1756 sel = self.get_selection()
1757 iflist = self.getifnos()
1758 if len(freqs)>len(iflist):
1759 raise ValueError("number of elements in list of list "
1760 "exeeds the current IF selections")
1761 iflist = self.getifnos()
1762 for i, fval in enumerate(freqs):
1763 sel.set_ifs(iflist[i])
1764 self._setselection(sel)
1765 self._setrestfreqs([fval], [""], unit)
1766 self._setselection(savesel)
1767
1768 # list values are dict, {'value'=, 'name'=)
1769 elif isinstance(freqs[-1], dict):
1770 values = []
1771 names = []
1772 for d in freqs:
1773 values.append(d["value"])
1774 names.append(d["name"])
1775 self._setrestfreqs(values, names, unit)
1776 elif isinstance(freqs[-1], list) or isinstance(freqs[-1], tuple):
1777 savesel = self._getselection()
1778 sel = self.get_selection()
1779 iflist = self.getifnos()
1780 if len(freqs)>len(iflist):
1781 raise ValueError("number of elements in list of list exeeds"
1782 " the current IF selections")
1783 for i, fval in enumerate(freqs):
1784 sel.set_ifs(iflist[i])
1785 self._setselection(sel)
1786 self._setrestfreqs(fval, [""], unit)
1787 self._setselection(savesel)
1788 # freqs are to be taken from a linecatalog
1789 elif isinstance(freqs, linecatalog):
1790 savesel = self._getselection()
1791 sel = self.get_selection()
1792 for i in xrange(freqs.nrow()):
1793 sel.set_ifs(iflist[i])
1794 self._setselection(sel)
1795 self._setrestfreqs([freqs.get_frequency(i)],
1796 [freqs.get_name(i)], "MHz")
1797 # ensure that we are not iterating past nIF
1798 if i == self.nif()-1: break
1799 self._setselection(savesel)
1800 else:
1801 return
1802 self._add_history("set_restfreqs", varlist)
1803
1804 @asaplog_post_dec
1805 def shift_refpix(self, delta):
1806 """\
1807 Shift the reference pixel of the Spectra Coordinate by an
1808 integer amount.
1809
1810 Parameters:
1811
1812 delta: the amount to shift by
1813
1814 *Note*:
1815
1816 Be careful using this with broadband data.
1817
1818 """
1819 varlist = vars()
1820 Scantable.shift_refpix(self, delta)
1821 s._add_history("shift_refpix", varlist)
1822
1823 @asaplog_post_dec
1824 def history(self, filename=None):
1825 """\
1826 Print the history. Optionally to a file.
1827
1828 Parameters:
1829
1830 filename: The name of the file to save the history to.
1831
1832 """
1833 hist = list(self._gethistory())
1834 out = "-"*80
1835 for h in hist:
1836 if h.startswith("---"):
1837 out = "\n".join([out, h])
1838 else:
1839 items = h.split("##")
1840 date = items[0]
1841 func = items[1]
1842 items = items[2:]
1843 out += "\n"+date+"\n"
1844 out += "Function: %s\n Parameters:" % (func)
1845 for i in items:
1846 if i == '':
1847 continue
1848 s = i.split("=")
1849 out += "\n %s = %s" % (s[0], s[1])
1850 out = "\n".join([out, "-"*80])
1851 if filename is not None:
1852 if filename is "":
1853 filename = 'scantable_history.txt'
1854 import os
1855 filename = os.path.expandvars(os.path.expanduser(filename))
1856 if not os.path.isdir(filename):
1857 data = open(filename, 'w')
1858 data.write(out)
1859 data.close()
1860 else:
1861 msg = "Illegal file name '%s'." % (filename)
1862 raise IOError(msg)
1863 return page(out)
1864
1865 #
1866 # Maths business
1867 #
1868 @asaplog_post_dec
1869 def average_time(self, mask=None, scanav=False, weight='tint', align=False):
1870 """\
1871 Return the (time) weighted average of a scan. Scans will be averaged
1872 only if the source direction (RA/DEC) is within 1' otherwise
1873
1874 *Note*:
1875
1876 in channels only - align if necessary
1877
1878 Parameters:
1879
1880 mask: an optional mask (only used for 'var' and 'tsys'
1881 weighting)
1882
1883 scanav: True averages each scan separately
1884 False (default) averages all scans together,
1885
1886 weight: Weighting scheme.
1887 'none' (mean no weight)
1888 'var' (1/var(spec) weighted)
1889 'tsys' (1/Tsys**2 weighted)
1890 'tint' (integration time weighted)
1891 'tintsys' (Tint/Tsys**2)
1892 'median' ( median averaging)
1893 The default is 'tint'
1894
1895 align: align the spectra in velocity before averaging. It takes
1896 the time of the first spectrum as reference time.
1897
1898 Example::
1899
1900 # time average the scantable without using a mask
1901 newscan = scan.average_time()
1902
1903 """
1904 varlist = vars()
1905 weight = weight or 'TINT'
1906 mask = mask or ()
1907 scanav = (scanav and 'SCAN') or 'NONE'
1908 scan = (self, )
1909
1910 if align:
1911 scan = (self.freq_align(insitu=False), )
1912 s = None
1913 if weight.upper() == 'MEDIAN':
1914 s = scantable(self._math._averagechannel(scan[0], 'MEDIAN',
1915 scanav))
1916 else:
1917 s = scantable(self._math._average(scan, mask, weight.upper(),
1918 scanav))
1919 s._add_history("average_time", varlist)
1920 return s
1921
1922 @asaplog_post_dec
1923 def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None):
1924 """\
1925 Return a scan where all spectra are converted to either
1926 Jansky or Kelvin depending upon the flux units of the scan table.
1927 By default the function tries to look the values up internally.
1928 If it can't find them (or if you want to over-ride), you must
1929 specify EITHER jyperk OR eta (and D which it will try to look up
1930 also if you don't set it). jyperk takes precedence if you set both.
1931
1932 Parameters:
1933
1934 jyperk: the Jy / K conversion factor
1935
1936 eta: the aperture efficiency
1937
1938 d: the geometric diameter (metres)
1939
1940 insitu: if False a new scantable is returned.
1941 Otherwise, the scaling is done in-situ
1942 The default is taken from .asaprc (False)
1943
1944 """
1945 if insitu is None: insitu = rcParams['insitu']
1946 self._math._setinsitu(insitu)
1947 varlist = vars()
1948 jyperk = jyperk or -1.0
1949 d = d or -1.0
1950 eta = eta or -1.0
1951 s = scantable(self._math._convertflux(self, d, eta, jyperk))
1952 s._add_history("convert_flux", varlist)
1953 if insitu: self._assign(s)
1954 else: return s
1955
1956 @asaplog_post_dec
1957 def gain_el(self, poly=None, filename="", method="linear", insitu=None):
1958 """\
1959 Return a scan after applying a gain-elevation correction.
1960 The correction can be made via either a polynomial or a
1961 table-based interpolation (and extrapolation if necessary).
1962 You specify polynomial coefficients, an ascii table or neither.
1963 If you specify neither, then a polynomial correction will be made
1964 with built in coefficients known for certain telescopes (an error
1965 will occur if the instrument is not known).
1966 The data and Tsys are *divided* by the scaling factors.
1967
1968 Parameters:
1969
1970 poly: Polynomial coefficients (default None) to compute a
1971 gain-elevation correction as a function of
1972 elevation (in degrees).
1973
1974 filename: The name of an ascii file holding correction factors.
1975 The first row of the ascii file must give the column
1976 names and these MUST include columns
1977 'ELEVATION' (degrees) and 'FACTOR' (multiply data
1978 by this) somewhere.
1979 The second row must give the data type of the
1980 column. Use 'R' for Real and 'I' for Integer.
1981 An example file would be
1982 (actual factors are arbitrary) :
1983
1984 TIME ELEVATION FACTOR
1985 R R R
1986 0.1 0 0.8
1987 0.2 20 0.85
1988 0.3 40 0.9
1989 0.4 60 0.85
1990 0.5 80 0.8
1991 0.6 90 0.75
1992
1993 method: Interpolation method when correcting from a table.
1994 Values are 'nearest', 'linear' (default), 'cubic'
1995 and 'spline'
1996
1997 insitu: if False a new scantable is returned.
1998 Otherwise, the scaling is done in-situ
1999 The default is taken from .asaprc (False)
2000
2001 """
2002
2003 if insitu is None: insitu = rcParams['insitu']
2004 self._math._setinsitu(insitu)
2005 varlist = vars()
2006 poly = poly or ()
2007 from os.path import expandvars
2008 filename = expandvars(filename)
2009 s = scantable(self._math._gainel(self, poly, filename, method))
2010 s._add_history("gain_el", varlist)
2011 if insitu:
2012 self._assign(s)
2013 else:
2014 return s
2015
2016 @asaplog_post_dec
2017 def freq_align(self, reftime=None, method='cubic', insitu=None):
2018 """\
2019 Return a scan where all rows have been aligned in frequency/velocity.
2020 The alignment frequency frame (e.g. LSRK) is that set by function
2021 set_freqframe.
2022
2023 Parameters:
2024
2025 reftime: reference time to align at. By default, the time of
2026 the first row of data is used.
2027
2028 method: Interpolation method for regridding the spectra.
2029 Choose from 'nearest', 'linear', 'cubic' (default)
2030 and 'spline'
2031
2032 insitu: if False a new scantable is returned.
2033 Otherwise, the scaling is done in-situ
2034 The default is taken from .asaprc (False)
2035
2036 """
2037 if insitu is None: insitu = rcParams["insitu"]
2038 oldInsitu = self._math._insitu()
2039 self._math._setinsitu(insitu)
2040 varlist = vars()
2041 reftime = reftime or ""
2042 s = scantable(self._math._freq_align(self, reftime, method))
2043 s._add_history("freq_align", varlist)
2044 self._math._setinsitu(oldInsitu)
2045 if insitu:
2046 self._assign(s)
2047 else:
2048 return s
2049
2050 @asaplog_post_dec
2051 def opacity(self, tau=None, insitu=None):
2052 """\
2053 Apply an opacity correction. The data
2054 and Tsys are multiplied by the correction factor.
2055
2056 Parameters:
2057
2058 tau: (list of) opacity from which the correction factor is
2059 exp(tau*ZD)
2060 where ZD is the zenith-distance.
2061 If a list is provided, it has to be of length nIF,
2062 nIF*nPol or 1 and in order of IF/POL, e.g.
2063 [opif0pol0, opif0pol1, opif1pol0 ...]
2064 if tau is `None` the opacities are determined from a
2065 model.
2066
2067 insitu: if False a new scantable is returned.
2068 Otherwise, the scaling is done in-situ
2069 The default is taken from .asaprc (False)
2070
2071 """
2072 if insitu is None:
2073 insitu = rcParams['insitu']
2074 self._math._setinsitu(insitu)
2075 varlist = vars()
2076 if not hasattr(tau, "__len__"):
2077 tau = [tau]
2078 s = scantable(self._math._opacity(self, tau))
2079 s._add_history("opacity", varlist)
2080 if insitu:
2081 self._assign(s)
2082 else:
2083 return s
2084
2085 @asaplog_post_dec
2086 def bin(self, width=5, insitu=None):
2087 """\
2088 Return a scan where all spectra have been binned up.
2089
2090 Parameters:
2091
2092 width: The bin width (default=5) in pixels
2093
2094 insitu: if False a new scantable is returned.
2095 Otherwise, the scaling is done in-situ
2096 The default is taken from .asaprc (False)
2097
2098 """
2099 if insitu is None:
2100 insitu = rcParams['insitu']
2101 self._math._setinsitu(insitu)
2102 varlist = vars()
2103 s = scantable(self._math._bin(self, width))
2104 s._add_history("bin", varlist)
2105 if insitu:
2106 self._assign(s)
2107 else:
2108 return s
2109
2110 @asaplog_post_dec
2111 def resample(self, width=5, method='cubic', insitu=None):
2112 """\
2113 Return a scan where all spectra have been binned up.
2114
2115 Parameters:
2116
2117 width: The bin width (default=5) in pixels
2118
2119 method: Interpolation method when correcting from a table.
2120 Values are 'nearest', 'linear', 'cubic' (default)
2121 and 'spline'
2122
2123 insitu: if False a new scantable is returned.
2124 Otherwise, the scaling is done in-situ
2125 The default is taken from .asaprc (False)
2126
2127 """
2128 if insitu is None:
2129 insitu = rcParams['insitu']
2130 self._math._setinsitu(insitu)
2131 varlist = vars()
2132 s = scantable(self._math._resample(self, method, width))
2133 s._add_history("resample", varlist)
2134 if insitu:
2135 self._assign(s)
2136 else:
2137 return s
2138
2139 @asaplog_post_dec
2140 def average_pol(self, mask=None, weight='none'):
2141 """\
2142 Average the Polarisations together.
2143
2144 Parameters:
2145
2146 mask: An optional mask defining the region, where the
2147 averaging will be applied. The output will have all
2148 specified points masked.
2149
2150 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
2151 weighted), or 'tsys' (1/Tsys**2 weighted)
2152
2153 """
2154 varlist = vars()
2155 mask = mask or ()
2156 s = scantable(self._math._averagepol(self, mask, weight.upper()))
2157 s._add_history("average_pol", varlist)
2158 return s
2159
2160 @asaplog_post_dec
2161 def average_beam(self, mask=None, weight='none'):
2162 """\
2163 Average the Beams together.
2164
2165 Parameters:
2166 mask: An optional mask defining the region, where the
2167 averaging will be applied. The output will have all
2168 specified points masked.
2169
2170 weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
2171 weighted), or 'tsys' (1/Tsys**2 weighted)
2172
2173 """
2174 varlist = vars()
2175 mask = mask or ()
2176 s = scantable(self._math._averagebeams(self, mask, weight.upper()))
2177 s._add_history("average_beam", varlist)
2178 return s
2179
2180 def parallactify(self, pflag):
2181 """\
2182 Set a flag to indicate whether this data should be treated as having
2183 been 'parallactified' (total phase == 0.0)
2184
2185 Parameters:
2186
2187 pflag: Bool indicating whether to turn this on (True) or
2188 off (False)
2189
2190 """
2191 varlist = vars()
2192 self._parallactify(pflag)
2193 self._add_history("parallactify", varlist)
2194
2195 @asaplog_post_dec
2196 def convert_pol(self, poltype=None):
2197 """\
2198 Convert the data to a different polarisation type.
2199 Note that you will need cross-polarisation terms for most conversions.
2200
2201 Parameters:
2202
2203 poltype: The new polarisation type. Valid types are:
2204 'linear', 'circular', 'stokes' and 'linpol'
2205
2206 """
2207 varlist = vars()
2208 s = scantable(self._math._convertpol(self, poltype))
2209 s._add_history("convert_pol", varlist)
2210 return s
2211
2212 @asaplog_post_dec
2213 def smooth(self, kernel="hanning", width=5.0, order=2, plot=False,
2214 insitu=None):
2215 """\
2216 Smooth the spectrum by the specified kernel (conserving flux).
2217
2218 Parameters:
2219
2220 kernel: The type of smoothing kernel. Select from
2221 'hanning' (default), 'gaussian', 'boxcar', 'rmedian'
2222 or 'poly'
2223
2224 width: The width of the kernel in pixels. For hanning this is
2225 ignored otherwise it defauls to 5 pixels.
2226 For 'gaussian' it is the Full Width Half
2227 Maximum. For 'boxcar' it is the full width.
2228 For 'rmedian' and 'poly' it is the half width.
2229
2230 order: Optional parameter for 'poly' kernel (default is 2), to
2231 specify the order of the polnomial. Ignored by all other
2232 kernels.
2233
2234 plot: plot the original and the smoothed spectra.
2235 In this each indivual fit has to be approved, by
2236 typing 'y' or 'n'
2237
2238 insitu: if False a new scantable is returned.
2239 Otherwise, the scaling is done in-situ
2240 The default is taken from .asaprc (False)
2241
2242 """
2243 if insitu is None: insitu = rcParams['insitu']
2244 self._math._setinsitu(insitu)
2245 varlist = vars()
2246
2247 if plot: orgscan = self.copy()
2248
2249 s = scantable(self._math._smooth(self, kernel.lower(), width, order))
2250 s._add_history("smooth", varlist)
2251
2252 if plot:
2253 from asap.asapplotter import new_asaplot
2254 theplot = new_asaplot(rcParams['plotter.gui'])
2255 from matplotlib import rc as rcp
2256 rcp('lines', linewidth=1)
2257 theplot.set_panels()
2258 ylab=s._get_ordinate_label()
2259 #theplot.palette(0,["#777777","red"])
2260 for r in xrange(s.nrow()):
2261 xsm=s._getabcissa(r)
2262 ysm=s._getspectrum(r)
2263 xorg=orgscan._getabcissa(r)
2264 yorg=orgscan._getspectrum(r)
2265 theplot.clear()
2266 theplot.hold()
2267 theplot.set_axes('ylabel',ylab)
2268 theplot.set_axes('xlabel',s._getabcissalabel(r))
2269 theplot.set_axes('title',s._getsourcename(r))
2270 theplot.set_line(label='Original',color="#777777")
2271 theplot.plot(xorg,yorg)
2272 theplot.set_line(label='Smoothed',color="red")
2273 theplot.plot(xsm,ysm)
2274 ### Ugly part for legend
2275 for i in [0,1]:
2276 theplot.subplots[0]['lines'].append(
2277 [theplot.subplots[0]['axes'].lines[i]]
2278 )
2279 theplot.release()
2280 ### Ugly part for legend
2281 theplot.subplots[0]['lines']=[]
2282 res = raw_input("Accept smoothing ([y]/n): ")
2283 if res.upper() == 'N':
2284 s._setspectrum(yorg, r)
2285 theplot.quit()
2286 del theplot
2287 del orgscan
2288
2289 if insitu: self._assign(s)
2290 else: return s
2291
2292 @asaplog_post_dec
2293 def regrid_channel(self, width=5, plot=False, insitu=None):
2294 """\
2295 Regrid the spectra by the specified channel width
2296
2297 Parameters:
2298
2299 width: The channel width (float) of regridded spectra
2300 in the current spectral unit.
2301
2302 plot: [NOT IMPLEMENTED YET]
2303 plot the original and the regridded spectra.
2304 In this each indivual fit has to be approved, by
2305 typing 'y' or 'n'
2306
2307 insitu: if False a new scantable is returned.
2308 Otherwise, the scaling is done in-situ
2309 The default is taken from .asaprc (False)
2310
2311 """
2312 if insitu is None: insitu = rcParams['insitu']
2313 varlist = vars()
2314
2315 if plot:
2316 asaplog.post()
2317 asaplog.push("Verification plot is not implemtnetd yet.")
2318 asaplog.post("WARN")
2319
2320 s = self.copy()
2321 s._regrid_specchan(width)
2322
2323 s._add_history("regrid_channel", varlist)
2324
2325# if plot:
2326# from asap.asapplotter import new_asaplot
2327# theplot = new_asaplot(rcParams['plotter.gui'])
2328# from matplotlib import rc as rcp
2329# rcp('lines', linewidth=1)
2330# theplot.set_panels()
2331# ylab=s._get_ordinate_label()
2332# #theplot.palette(0,["#777777","red"])
2333# for r in xrange(s.nrow()):
2334# xsm=s._getabcissa(r)
2335# ysm=s._getspectrum(r)
2336# xorg=orgscan._getabcissa(r)
2337# yorg=orgscan._getspectrum(r)
2338# theplot.clear()
2339# theplot.hold()
2340# theplot.set_axes('ylabel',ylab)
2341# theplot.set_axes('xlabel',s._getabcissalabel(r))
2342# theplot.set_axes('title',s._getsourcename(r))
2343# theplot.set_line(label='Original',color="#777777")
2344# theplot.plot(xorg,yorg)
2345# theplot.set_line(label='Smoothed',color="red")
2346# theplot.plot(xsm,ysm)
2347# ### Ugly part for legend
2348# for i in [0,1]:
2349# theplot.subplots[0]['lines'].append(
2350# [theplot.subplots[0]['axes'].lines[i]]
2351# )
2352# theplot.release()
2353# ### Ugly part for legend
2354# theplot.subplots[0]['lines']=[]
2355# res = raw_input("Accept smoothing ([y]/n): ")
2356# if res.upper() == 'N':
2357# s._setspectrum(yorg, r)
2358# theplot.quit()
2359# del theplot
2360# del orgscan
2361
2362 if insitu: self._assign(s)
2363 else: return s
2364
2365 @asaplog_post_dec
2366 def _parse_wn(self, wn):
2367 if isinstance(wn, list) or isinstance(wn, tuple):
2368 return wn
2369 elif isinstance(wn, int):
2370 return [ wn ]
2371 elif isinstance(wn, str):
2372 if '-' in wn: # case 'a-b' : return [a,a+1,...,b-1,b]
2373 val = wn.split('-')
2374 val = [int(val[0]), int(val[1])]
2375 val.sort()
2376 res = [i for i in xrange(val[0], val[1]+1)]
2377 elif wn[:2] == '<=' or wn[:2] == '=<': # cases '<=a','=<a' : return [0,1,...,a-1,a]
2378 val = int(wn[2:])+1
2379 res = [i for i in xrange(val)]
2380 elif wn[-2:] == '>=' or wn[-2:] == '=>': # cases 'a>=','a=>' : return [0,1,...,a-1,a]
2381 val = int(wn[:-2])+1
2382 res = [i for i in xrange(val)]
2383 elif wn[0] == '<': # case '<a' : return [0,1,...,a-2,a-1]
2384 val = int(wn[1:])
2385 res = [i for i in xrange(val)]
2386 elif wn[-1] == '>': # case 'a>' : return [0,1,...,a-2,a-1]
2387 val = int(wn[:-1])
2388 res = [i for i in xrange(val)]
2389 elif wn[:2] == '>=' or wn[:2] == '=>': # cases '>=a','=>a' : return [a,-999], which is
2390 # then interpreted in C++
2391 # side as [a,a+1,...,a_nyq]
2392 # (CAS-3759)
2393 val = int(wn[2:])
2394 res = [val, -999]
2395 #res = [i for i in xrange(val, self.nchan()/2+1)]
2396 elif wn[-2:] == '<=' or wn[-2:] == '=<': # cases 'a<=','a=<' : return [a,-999], which is
2397 # then interpreted in C++
2398 # side as [a,a+1,...,a_nyq]
2399 # (CAS-3759)
2400 val = int(wn[:-2])
2401 res = [val, -999]
2402 #res = [i for i in xrange(val, self.nchan()/2+1)]
2403 elif wn[0] == '>': # case '>a' : return [a+1,-999], which is
2404 # then interpreted in C++
2405 # side as [a+1,a+2,...,a_nyq]
2406 # (CAS-3759)
2407 val = int(wn[1:])+1
2408 res = [val, -999]
2409 #res = [i for i in xrange(val, self.nchan()/2+1)]
2410 elif wn[-1] == '<': # case 'a<' : return [a+1,-999], which is
2411 # then interpreted in C++
2412 # side as [a+1,a+2,...,a_nyq]
2413 # (CAS-3759)
2414 val = int(wn[:-1])+1
2415 res = [val, -999]
2416 #res = [i for i in xrange(val, self.nchan()/2+1)]
2417
2418 return res
2419 else:
2420 msg = 'wrong value given for addwn/rejwn'
2421 raise RuntimeError(msg)
2422
2423
2424 @asaplog_post_dec
2425 def sinusoid_baseline(self, insitu=None, mask=None, applyfft=None,
2426 fftmethod=None, fftthresh=None,
2427 addwn=None, rejwn=None, clipthresh=None,
2428 clipniter=None, plot=None,
2429 getresidual=None, showprogress=None,
2430 minnrow=None, outlog=None, blfile=None):
2431 """\
2432 Return a scan which has been baselined (all rows) with sinusoidal
2433 functions.
2434
2435 Parameters:
2436 insitu: if False a new scantable is returned.
2437 Otherwise, the scaling is done in-situ
2438 The default is taken from .asaprc (False)
2439 mask: an optional mask
2440 applyfft: if True use some method, such as FFT, to find
2441 strongest sinusoidal components in the wavenumber
2442 domain to be used for baseline fitting.
2443 default is True.
2444 fftmethod: method to find the strong sinusoidal components.
2445 now only 'fft' is available and it is the default.
2446 fftthresh: the threshold to select wave numbers to be used for
2447 fitting from the distribution of amplitudes in the
2448 wavenumber domain.
2449 both float and string values accepted.
2450 given a float value, the unit is set to sigma.
2451 for string values, allowed formats include:
2452 'xsigma' or 'x' (= x-sigma level. e.g.,
2453 '3sigma'), or
2454 'topx' (= the x strongest ones, e.g. 'top5').
2455 default is 3.0 (unit: sigma).
2456 addwn: the additional wave numbers to be used for fitting.
2457 list or integer value is accepted to specify every
2458 wave numbers. also string value can be used in case
2459 you need to specify wave numbers in a certain range,
2460 e.g., 'a-b' (= a, a+1, a+2, ..., b-1, b),
2461 '<a' (= 0,1,...,a-2,a-1),
2462 '>=a' (= a, a+1, ... up to the maximum wave
2463 number corresponding to the Nyquist
2464 frequency for the case of FFT).
2465 default is [0].
2466 rejwn: the wave numbers NOT to be used for fitting.
2467 can be set just as addwn but has higher priority:
2468 wave numbers which are specified both in addwn
2469 and rejwn will NOT be used. default is [].
2470 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2471 clipniter: maximum number of iteration of 'clipthresh'-sigma
2472 clipping (default is 0)
2473 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2474 plot the fit and the residual. In this each
2475 indivual fit has to be approved, by typing 'y'
2476 or 'n'
2477 getresidual: if False, returns best-fit values instead of
2478 residual. (default is True)
2479 showprogress: show progress status for large data.
2480 default is True.
2481 minnrow: minimum number of input spectra to show.
2482 default is 1000.
2483 outlog: Output the coefficients of the best-fit
2484 function to logger (default is False)
2485 blfile: Name of a text file in which the best-fit
2486 parameter values to be written
2487 (default is '': no file/logger output)
2488
2489 Example:
2490 # return a scan baselined by a combination of sinusoidal curves
2491 # having wave numbers in spectral window up to 10,
2492 # also with 3-sigma clipping, iteration up to 4 times
2493 bscan = scan.sinusoid_baseline(addwn='<=10',clipthresh=3.0,clipniter=4)
2494
2495 Note:
2496 The best-fit parameter values output in logger and/or blfile are now
2497 based on specunit of 'channel'.
2498 """
2499
2500 try:
2501 varlist = vars()
2502
2503 if insitu is None: insitu = rcParams['insitu']
2504 if insitu:
2505 workscan = self
2506 else:
2507 workscan = self.copy()
2508
2509 #if mask is None: mask = [True for i in xrange(workscan.nchan())]
2510 if mask is None: mask = []
2511 if applyfft is None: applyfft = True
2512 if fftmethod is None: fftmethod = 'fft'
2513 if fftthresh is None: fftthresh = 3.0
2514 if addwn is None: addwn = [0]
2515 if rejwn is None: rejwn = []
2516 if clipthresh is None: clipthresh = 3.0
2517 if clipniter is None: clipniter = 0
2518 if plot is None: plot = False
2519 if getresidual is None: getresidual = True
2520 if showprogress is None: showprogress = True
2521 if minnrow is None: minnrow = 1000
2522 if outlog is None: outlog = False
2523 if blfile is None: blfile = ''
2524
2525 #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
2526 workscan._sinusoid_baseline(mask, applyfft, fftmethod.lower(),
2527 str(fftthresh).lower(),
2528 workscan._parse_wn(addwn),
2529 workscan._parse_wn(rejwn), clipthresh,
2530 clipniter, getresidual,
2531 pack_progress_params(showprogress,
2532 minnrow), outlog,
2533 blfile)
2534 workscan._add_history('sinusoid_baseline', varlist)
2535
2536 if insitu:
2537 self._assign(workscan)
2538 else:
2539 return workscan
2540
2541 except RuntimeError, e:
2542 raise_fitting_failure_exception(e)
2543
2544
2545 @asaplog_post_dec
2546 def auto_sinusoid_baseline(self, insitu=None, mask=None, applyfft=None,
2547 fftmethod=None, fftthresh=None,
2548 addwn=None, rejwn=None, clipthresh=None,
2549 clipniter=None, edge=None, threshold=None,
2550 chan_avg_limit=None, plot=None,
2551 getresidual=None, showprogress=None,
2552 minnrow=None, outlog=None, blfile=None):
2553 """\
2554 Return a scan which has been baselined (all rows) with sinusoidal
2555 functions.
2556 Spectral lines are detected first using linefinder and masked out
2557 to avoid them affecting the baseline solution.
2558
2559 Parameters:
2560 insitu: if False a new scantable is returned.
2561 Otherwise, the scaling is done in-situ
2562 The default is taken from .asaprc (False)
2563 mask: an optional mask retreived from scantable
2564 applyfft: if True use some method, such as FFT, to find
2565 strongest sinusoidal components in the wavenumber
2566 domain to be used for baseline fitting.
2567 default is True.
2568 fftmethod: method to find the strong sinusoidal components.
2569 now only 'fft' is available and it is the default.
2570 fftthresh: the threshold to select wave numbers to be used for
2571 fitting from the distribution of amplitudes in the
2572 wavenumber domain.
2573 both float and string values accepted.
2574 given a float value, the unit is set to sigma.
2575 for string values, allowed formats include:
2576 'xsigma' or 'x' (= x-sigma level. e.g.,
2577 '3sigma'), or
2578 'topx' (= the x strongest ones, e.g. 'top5').
2579 default is 3.0 (unit: sigma).
2580 addwn: the additional wave numbers to be used for fitting.
2581 list or integer value is accepted to specify every
2582 wave numbers. also string value can be used in case
2583 you need to specify wave numbers in a certain range,
2584 e.g., 'a-b' (= a, a+1, a+2, ..., b-1, b),
2585 '<a' (= 0,1,...,a-2,a-1),
2586 '>=a' (= a, a+1, ... up to the maximum wave
2587 number corresponding to the Nyquist
2588 frequency for the case of FFT).
2589 default is [0].
2590 rejwn: the wave numbers NOT to be used for fitting.
2591 can be set just as addwn but has higher priority:
2592 wave numbers which are specified both in addwn
2593 and rejwn will NOT be used. default is [].
2594 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2595 clipniter: maximum number of iteration of 'clipthresh'-sigma
2596 clipping (default is 0)
2597 edge: an optional number of channel to drop at
2598 the edge of spectrum. If only one value is
2599 specified, the same number will be dropped
2600 from both sides of the spectrum. Default
2601 is to keep all channels. Nested tuples
2602 represent individual edge selection for
2603 different IFs (a number of spectral channels
2604 can be different)
2605 threshold: the threshold used by line finder. It is
2606 better to keep it large as only strong lines
2607 affect the baseline solution.
2608 chan_avg_limit: a maximum number of consequtive spectral
2609 channels to average during the search of
2610 weak and broad lines. The default is no
2611 averaging (and no search for weak lines).
2612 If such lines can affect the fitted baseline
2613 (e.g. a high order polynomial is fitted),
2614 increase this parameter (usually values up
2615 to 8 are reasonable). Most users of this
2616 method should find the default value sufficient.
2617 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2618 plot the fit and the residual. In this each
2619 indivual fit has to be approved, by typing 'y'
2620 or 'n'
2621 getresidual: if False, returns best-fit values instead of
2622 residual. (default is True)
2623 showprogress: show progress status for large data.
2624 default is True.
2625 minnrow: minimum number of input spectra to show.
2626 default is 1000.
2627 outlog: Output the coefficients of the best-fit
2628 function to logger (default is False)
2629 blfile: Name of a text file in which the best-fit
2630 parameter values to be written
2631 (default is "": no file/logger output)
2632
2633 Example:
2634 bscan = scan.auto_sinusoid_baseline(addwn='<=10', insitu=False)
2635
2636 Note:
2637 The best-fit parameter values output in logger and/or blfile are now
2638 based on specunit of 'channel'.
2639 """
2640
2641 try:
2642 varlist = vars()
2643
2644 if insitu is None: insitu = rcParams['insitu']
2645 if insitu:
2646 workscan = self
2647 else:
2648 workscan = self.copy()
2649
2650 #if mask is None: mask = [True for i in xrange(workscan.nchan())]
2651 if mask is None: mask = []
2652 if applyfft is None: applyfft = True
2653 if fftmethod is None: fftmethod = 'fft'
2654 if fftthresh is None: fftthresh = 3.0
2655 if addwn is None: addwn = [0]
2656 if rejwn is None: rejwn = []
2657 if clipthresh is None: clipthresh = 3.0
2658 if clipniter is None: clipniter = 0
2659 if edge is None: edge = (0,0)
2660 if threshold is None: threshold = 3
2661 if chan_avg_limit is None: chan_avg_limit = 1
2662 if plot is None: plot = False
2663 if getresidual is None: getresidual = True
2664 if showprogress is None: showprogress = True
2665 if minnrow is None: minnrow = 1000
2666 if outlog is None: outlog = False
2667 if blfile is None: blfile = ''
2668
2669 #CURRENTLY, PLOT=true is UNAVAILABLE UNTIL sinusoidal fitting is implemented as a fitter method.
2670 workscan._auto_sinusoid_baseline(mask, applyfft,
2671 fftmethod.lower(),
2672 str(fftthresh).lower(),
2673 workscan._parse_wn(addwn),
2674 workscan._parse_wn(rejwn),
2675 clipthresh, clipniter,
2676 normalise_edge_param(edge),
2677 threshold, chan_avg_limit,
2678 getresidual,
2679 pack_progress_params(showprogress,
2680 minnrow),
2681 outlog, blfile)
2682 workscan._add_history("auto_sinusoid_baseline", varlist)
2683
2684 if insitu:
2685 self._assign(workscan)
2686 else:
2687 return workscan
2688
2689 except RuntimeError, e:
2690 raise_fitting_failure_exception(e)
2691
2692 @asaplog_post_dec
2693 def cspline_baseline(self, insitu=None, mask=None, npiece=None,
2694 clipthresh=None, clipniter=None, plot=None,
2695 getresidual=None, showprogress=None, minnrow=None,
2696 outlog=None, blfile=None):
2697 """\
2698 Return a scan which has been baselined (all rows) by cubic spline
2699 function (piecewise cubic polynomial).
2700
2701 Parameters:
2702 insitu: If False a new scantable is returned.
2703 Otherwise, the scaling is done in-situ
2704 The default is taken from .asaprc (False)
2705 mask: An optional mask
2706 npiece: Number of pieces. (default is 2)
2707 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2708 clipniter: maximum number of iteration of 'clipthresh'-sigma
2709 clipping (default is 0)
2710 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2711 plot the fit and the residual. In this each
2712 indivual fit has to be approved, by typing 'y'
2713 or 'n'
2714 getresidual: if False, returns best-fit values instead of
2715 residual. (default is True)
2716 showprogress: show progress status for large data.
2717 default is True.
2718 minnrow: minimum number of input spectra to show.
2719 default is 1000.
2720 outlog: Output the coefficients of the best-fit
2721 function to logger (default is False)
2722 blfile: Name of a text file in which the best-fit
2723 parameter values to be written
2724 (default is "": no file/logger output)
2725
2726 Example:
2727 # return a scan baselined by a cubic spline consisting of 2 pieces
2728 # (i.e., 1 internal knot),
2729 # also with 3-sigma clipping, iteration up to 4 times
2730 bscan = scan.cspline_baseline(npiece=2,clipthresh=3.0,clipniter=4)
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 npiece is None: npiece = 2
2749 if clipthresh is None: clipthresh = 3.0
2750 if clipniter is None: clipniter = 0
2751 if plot is None: plot = False
2752 if getresidual is None: getresidual = True
2753 if showprogress is None: showprogress = True
2754 if minnrow is None: minnrow = 1000
2755 if outlog is None: outlog = False
2756 if blfile is None: blfile = ''
2757
2758 #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
2759 workscan._cspline_baseline(mask, npiece, clipthresh, clipniter,
2760 getresidual,
2761 pack_progress_params(showprogress,
2762 minnrow), outlog,
2763 blfile)
2764 workscan._add_history("cspline_baseline", varlist)
2765
2766 if insitu:
2767 self._assign(workscan)
2768 else:
2769 return workscan
2770
2771 except RuntimeError, e:
2772 raise_fitting_failure_exception(e)
2773
2774 @asaplog_post_dec
2775 def auto_cspline_baseline(self, insitu=None, mask=None, npiece=None,
2776 clipthresh=None, clipniter=None,
2777 edge=None, threshold=None, chan_avg_limit=None,
2778 getresidual=None, plot=None,
2779 showprogress=None, minnrow=None, outlog=None,
2780 blfile=None):
2781 """\
2782 Return a scan which has been baselined (all rows) by cubic spline
2783 function (piecewise cubic polynomial).
2784 Spectral lines are detected first using linefinder and masked out
2785 to avoid them affecting the baseline solution.
2786
2787 Parameters:
2788 insitu: if False a new scantable is returned.
2789 Otherwise, the scaling is done in-situ
2790 The default is taken from .asaprc (False)
2791 mask: an optional mask retreived from scantable
2792 npiece: Number of pieces. (default is 2)
2793 clipthresh: Clipping threshold. (default is 3.0, unit: sigma)
2794 clipniter: maximum number of iteration of 'clipthresh'-sigma
2795 clipping (default is 0)
2796 edge: an optional number of channel to drop at
2797 the edge of spectrum. If only one value is
2798 specified, the same number will be dropped
2799 from both sides of the spectrum. Default
2800 is to keep all channels. Nested tuples
2801 represent individual edge selection for
2802 different IFs (a number of spectral channels
2803 can be different)
2804 threshold: the threshold used by line finder. It is
2805 better to keep it large as only strong lines
2806 affect the baseline solution.
2807 chan_avg_limit: a maximum number of consequtive spectral
2808 channels to average during the search of
2809 weak and broad lines. The default is no
2810 averaging (and no search for weak lines).
2811 If such lines can affect the fitted baseline
2812 (e.g. a high order polynomial is fitted),
2813 increase this parameter (usually values up
2814 to 8 are reasonable). Most users of this
2815 method should find the default value sufficient.
2816 plot: *** CURRENTLY UNAVAILABLE, ALWAYS FALSE ***
2817 plot the fit and the residual. In this each
2818 indivual fit has to be approved, by typing 'y'
2819 or 'n'
2820 getresidual: if False, returns best-fit values instead of
2821 residual. (default is True)
2822 showprogress: show progress status for large data.
2823 default is True.
2824 minnrow: minimum number of input spectra to show.
2825 default is 1000.
2826 outlog: Output the coefficients of the best-fit
2827 function to logger (default is False)
2828 blfile: Name of a text file in which the best-fit
2829 parameter values to be written
2830 (default is "": no file/logger output)
2831
2832 Example:
2833 bscan = scan.auto_cspline_baseline(npiece=3, insitu=False)
2834
2835 Note:
2836 The best-fit parameter values output in logger and/or blfile are now
2837 based on specunit of 'channel'.
2838 """
2839
2840 try:
2841 varlist = vars()
2842
2843 if insitu is None: insitu = rcParams['insitu']
2844 if insitu:
2845 workscan = self
2846 else:
2847 workscan = self.copy()
2848
2849 #if mask is None: mask = [True for i in xrange(workscan.nchan())]
2850 if mask is None: mask = []
2851 if npiece is None: npiece = 2
2852 if clipthresh is None: clipthresh = 3.0
2853 if clipniter is None: clipniter = 0
2854 if edge is None: edge = (0, 0)
2855 if threshold is None: threshold = 3
2856 if chan_avg_limit is None: chan_avg_limit = 1
2857 if plot is None: plot = False
2858 if getresidual is None: getresidual = True
2859 if showprogress is None: showprogress = True
2860 if minnrow is None: minnrow = 1000
2861 if outlog is None: outlog = False
2862 if blfile is None: blfile = ''
2863
2864 #CURRENTLY, PLOT=true UNAVAILABLE UNTIL cubic spline fitting is implemented as a fitter method.
2865 workscan._auto_cspline_baseline(mask, npiece, clipthresh,
2866 clipniter,
2867 normalise_edge_param(edge),
2868 threshold,
2869 chan_avg_limit, getresidual,
2870 pack_progress_params(showprogress,
2871 minnrow),
2872 outlog, blfile)
2873 workscan._add_history("auto_cspline_baseline", varlist)
2874
2875 if insitu:
2876 self._assign(workscan)
2877 else:
2878 return workscan
2879
2880 except RuntimeError, e:
2881 raise_fitting_failure_exception(e)
2882
2883 @asaplog_post_dec
2884 def poly_baseline(self, mask=None, order=None, insitu=None, plot=None,
2885 getresidual=None, showprogress=None, minnrow=None,
2886 outlog=None, blfile=None):
2887 """\
2888 Return a scan which has been baselined (all rows) by a polynomial.
2889 Parameters:
2890 insitu: if False a new scantable is returned.
2891 Otherwise, the scaling is done in-situ
2892 The default is taken from .asaprc (False)
2893 mask: an optional mask
2894 order: the order of the polynomial (default is 0)
2895 plot: plot the fit and the residual. In this each
2896 indivual fit has to be approved, by typing 'y'
2897 or 'n'
2898 getresidual: if False, returns best-fit values instead of
2899 residual. (default is True)
2900 showprogress: show progress status for large data.
2901 default is True.
2902 minnrow: minimum number of input spectra to show.
2903 default is 1000.
2904 outlog: Output the coefficients of the best-fit
2905 function to logger (default is False)
2906 blfile: Name of a text file in which the best-fit
2907 parameter values to be written
2908 (default is "": no file/logger output)
2909
2910 Example:
2911 # return a scan baselined by a third order polynomial,
2912 # not using a mask
2913 bscan = scan.poly_baseline(order=3)
2914 """
2915
2916 try:
2917 varlist = vars()
2918
2919 if insitu is None:
2920 insitu = rcParams["insitu"]
2921 if insitu:
2922 workscan = self
2923 else:
2924 workscan = self.copy()
2925
2926 #if mask is None: mask = [True for i in \
2927 # xrange(workscan.nchan())]
2928 if mask is None: mask = []
2929 if order is None: order = 0
2930 if plot is None: plot = False
2931 if getresidual is None: getresidual = True
2932 if showprogress is None: showprogress = True
2933 if minnrow is None: minnrow = 1000
2934 if outlog is None: outlog = False
2935 if blfile is None: blfile = ""
2936
2937 if plot:
2938 outblfile = (blfile != "") and \
2939 os.path.exists(os.path.expanduser(
2940 os.path.expandvars(blfile))
2941 )
2942 if outblfile:
2943 blf = open(blfile, "a")
2944
2945 f = fitter()
2946 f.set_function(lpoly=order)
2947
2948 rows = xrange(workscan.nrow())
2949 #if len(rows) > 0: workscan._init_blinfo()
2950
2951 for r in rows:
2952 f.x = workscan._getabcissa(r)
2953 f.y = workscan._getspectrum(r)
2954 if mask:
2955 f.mask = mask_and(mask, workscan._getmask(r)) # (CAS-1434)
2956 else: # mask=None
2957 f.mask = workscan._getmask(r)
2958
2959 f.data = None
2960 f.fit()
2961
2962 f.plot(residual=True)
2963 accept_fit = raw_input("Accept fit ( [y]/n ): ")
2964 if accept_fit.upper() == "N":
2965 #workscan._append_blinfo(None, None, None)
2966 continue
2967
2968 blpars = f.get_parameters()
2969 masklist = workscan.get_masklist(f.mask, row=r, silent=True)
2970 #workscan._append_blinfo(blpars, masklist, f.mask)
2971 workscan._setspectrum((f.fitter.getresidual()
2972 if getresidual else f.fitter.getfit()), r)
2973
2974 if outblfile:
2975 rms = workscan.get_rms(f.mask, r)
2976 dataout = \
2977 workscan.format_blparams_row(blpars["params"],
2978 blpars["fixed"],
2979 rms, str(masklist),
2980 r, True)
2981 blf.write(dataout)
2982
2983 f._p.unmap()
2984 f._p = None
2985
2986 if outblfile:
2987 blf.close()
2988 else:
2989 workscan._poly_baseline(mask, order, getresidual,
2990 pack_progress_params(showprogress,
2991 minnrow),
2992 outlog, blfile)
2993
2994 workscan._add_history("poly_baseline", varlist)
2995
2996 if insitu:
2997 self._assign(workscan)
2998 else:
2999 return workscan
3000
3001 except RuntimeError, e:
3002 raise_fitting_failure_exception(e)
3003
3004 @asaplog_post_dec
3005 def auto_poly_baseline(self, mask=None, order=None, edge=None,
3006 threshold=None, chan_avg_limit=None,
3007 plot=None, insitu=None,
3008 getresidual=None, showprogress=None,
3009 minnrow=None, outlog=None, blfile=None):
3010 """\
3011 Return a scan which has been baselined (all rows) by a polynomial.
3012 Spectral lines are detected first using linefinder and masked out
3013 to avoid them affecting the baseline solution.
3014
3015 Parameters:
3016 insitu: if False a new scantable is returned.
3017 Otherwise, the scaling is done in-situ
3018 The default is taken from .asaprc (False)
3019 mask: an optional mask retreived from scantable
3020 order: the order of the polynomial (default is 0)
3021 edge: an optional number of channel to drop at
3022 the edge of spectrum. If only one value is
3023 specified, the same number will be dropped
3024 from both sides of the spectrum. Default
3025 is to keep all channels. Nested tuples
3026 represent individual edge selection for
3027 different IFs (a number of spectral channels
3028 can be different)
3029 threshold: the threshold used by line finder. It is
3030 better to keep it large as only strong lines
3031 affect the baseline solution.
3032 chan_avg_limit: a maximum number of consequtive spectral
3033 channels to average during the search of
3034 weak and broad lines. The default is no
3035 averaging (and no search for weak lines).
3036 If such lines can affect the fitted baseline
3037 (e.g. a high order polynomial is fitted),
3038 increase this parameter (usually values up
3039 to 8 are reasonable). Most users of this
3040 method should find the default value sufficient.
3041 plot: plot the fit and the residual. In this each
3042 indivual fit has to be approved, by typing 'y'
3043 or 'n'
3044 getresidual: if False, returns best-fit values instead of
3045 residual. (default is True)
3046 showprogress: show progress status for large data.
3047 default is True.
3048 minnrow: minimum number of input spectra to show.
3049 default is 1000.
3050 outlog: Output the coefficients of the best-fit
3051 function to logger (default is False)
3052 blfile: Name of a text file in which the best-fit
3053 parameter values to be written
3054 (default is "": no file/logger output)
3055
3056 Example:
3057 bscan = scan.auto_poly_baseline(order=7, insitu=False)
3058 """
3059
3060 try:
3061 varlist = vars()
3062
3063 if insitu is None:
3064 insitu = rcParams['insitu']
3065 if insitu:
3066 workscan = self
3067 else:
3068 workscan = self.copy()
3069
3070 #if mask is None: mask = [True for i in xrange(workscan.nchan())]
3071 if mask is None: mask = []
3072 if order is None: order = 0
3073 if edge is None: edge = (0, 0)
3074 if threshold is None: threshold = 3
3075 if chan_avg_limit is None: chan_avg_limit = 1
3076 if plot is None: plot = False
3077 if getresidual is None: getresidual = True
3078 if showprogress is None: showprogress = True
3079 if minnrow is None: minnrow = 1000
3080 if outlog is None: outlog = False
3081 if blfile is None: blfile = ''
3082
3083 edge = normalise_edge_param(edge)
3084
3085 if plot:
3086 outblfile = (blfile != "") and \
3087 os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
3088 if outblfile: blf = open(blfile, "a")
3089
3090 from asap.asaplinefind import linefinder
3091 fl = linefinder()
3092 fl.set_options(threshold=threshold, avg_limit=chan_avg_limit)
3093 fl.set_scan(workscan)
3094
3095 f = fitter()
3096 f.set_function(lpoly=order)
3097
3098 rows = xrange(workscan.nrow())
3099 #if len(rows) > 0: workscan._init_blinfo()
3100
3101 for r in rows:
3102 idx = 2*workscan.getif(r)
3103 if mask:
3104 msk = mask_and(mask, workscan._getmask(r)) # (CAS-1434)
3105 else: # mask=None
3106 msk = workscan._getmask(r)
3107 fl.find_lines(r, msk, edge[idx:idx+2])
3108
3109 f.x = workscan._getabcissa(r)
3110 f.y = workscan._getspectrum(r)
3111 f.mask = fl.get_mask()
3112 f.data = None
3113 f.fit()
3114
3115 f.plot(residual=True)
3116 accept_fit = raw_input("Accept fit ( [y]/n ): ")
3117 if accept_fit.upper() == "N":
3118 #workscan._append_blinfo(None, None, None)
3119 continue
3120
3121 blpars = f.get_parameters()
3122 masklist = workscan.get_masklist(f.mask, row=r, silent=True)
3123 #workscan._append_blinfo(blpars, masklist, f.mask)
3124 workscan._setspectrum(
3125 (f.fitter.getresidual() if getresidual
3126 else f.fitter.getfit()), r
3127 )
3128
3129 if outblfile:
3130 rms = workscan.get_rms(f.mask, r)
3131 dataout = \
3132 workscan.format_blparams_row(blpars["params"],
3133 blpars["fixed"],
3134 rms, str(masklist),
3135 r, True)
3136 blf.write(dataout)
3137
3138 f._p.unmap()
3139 f._p = None
3140
3141 if outblfile: blf.close()
3142 else:
3143 workscan._auto_poly_baseline(mask, order, edge, threshold,
3144 chan_avg_limit, getresidual,
3145 pack_progress_params(showprogress,
3146 minnrow),
3147 outlog, blfile)
3148
3149 workscan._add_history("auto_poly_baseline", varlist)
3150
3151 if insitu:
3152 self._assign(workscan)
3153 else:
3154 return workscan
3155
3156 except RuntimeError, e:
3157 raise_fitting_failure_exception(e)
3158
3159 def _init_blinfo(self):
3160 """\
3161 Initialise the following three auxiliary members:
3162 blpars : parameters of the best-fit baseline,
3163 masklists : mask data (edge positions of masked channels) and
3164 actualmask : mask data (in boolean list),
3165 to keep for use later (including output to logger/text files).
3166 Used by poly_baseline() and auto_poly_baseline() in case of
3167 'plot=True'.
3168 """
3169 self.blpars = []
3170 self.masklists = []
3171 self.actualmask = []
3172 return
3173
3174 def _append_blinfo(self, data_blpars, data_masklists, data_actualmask):
3175 """\
3176 Append baseline-fitting related info to blpars, masklist and
3177 actualmask.
3178 """
3179 self.blpars.append(data_blpars)
3180 self.masklists.append(data_masklists)
3181 self.actualmask.append(data_actualmask)
3182 return
3183
3184 @asaplog_post_dec
3185 def rotate_linpolphase(self, angle):
3186 """\
3187 Rotate the phase of the complex polarization O=Q+iU correlation.
3188 This is always done in situ in the raw data. So if you call this
3189 function more than once then each call rotates the phase further.
3190
3191 Parameters:
3192
3193 angle: The angle (degrees) to rotate (add) by.
3194
3195 Example::
3196
3197 scan.rotate_linpolphase(2.3)
3198
3199 """
3200 varlist = vars()
3201 self._math._rotate_linpolphase(self, angle)
3202 self._add_history("rotate_linpolphase", varlist)
3203 return
3204
3205 @asaplog_post_dec
3206 def rotate_xyphase(self, angle):
3207 """\
3208 Rotate the phase of the XY correlation. This is always done in situ
3209 in the data. So if you call this function more than once
3210 then each call rotates the phase further.
3211
3212 Parameters:
3213
3214 angle: The angle (degrees) to rotate (add) by.
3215
3216 Example::
3217
3218 scan.rotate_xyphase(2.3)
3219
3220 """
3221 varlist = vars()
3222 self._math._rotate_xyphase(self, angle)
3223 self._add_history("rotate_xyphase", varlist)
3224 return
3225
3226 @asaplog_post_dec
3227 def swap_linears(self):
3228 """\
3229 Swap the linear polarisations XX and YY, or better the first two
3230 polarisations as this also works for ciculars.
3231 """
3232 varlist = vars()
3233 self._math._swap_linears(self)
3234 self._add_history("swap_linears", varlist)
3235 return
3236
3237 @asaplog_post_dec
3238 def invert_phase(self):
3239 """\
3240 Invert the phase of the complex polarisation
3241 """
3242 varlist = vars()
3243 self._math._invert_phase(self)
3244 self._add_history("invert_phase", varlist)
3245 return
3246
3247 @asaplog_post_dec
3248 def add(self, offset, insitu=None):
3249 """\
3250 Return a scan where all spectra have the offset added
3251
3252 Parameters:
3253
3254 offset: the offset
3255
3256 insitu: if False a new scantable is returned.
3257 Otherwise, the scaling is done in-situ
3258 The default is taken from .asaprc (False)
3259
3260 """
3261 if insitu is None: insitu = rcParams['insitu']
3262 self._math._setinsitu(insitu)
3263 varlist = vars()
3264 s = scantable(self._math._unaryop(self, offset, "ADD", False))
3265 s._add_history("add", varlist)
3266 if insitu:
3267 self._assign(s)
3268 else:
3269 return s
3270
3271 @asaplog_post_dec
3272 def scale(self, factor, tsys=True, insitu=None):
3273 """\
3274
3275 Return a scan where all spectra are scaled by the given 'factor'
3276
3277 Parameters:
3278
3279 factor: the scaling factor (float or 1D float list)
3280
3281 insitu: if False a new scantable is returned.
3282 Otherwise, the scaling is done in-situ
3283 The default is taken from .asaprc (False)
3284
3285 tsys: if True (default) then apply the operation to Tsys
3286 as well as the data
3287
3288 """
3289 if insitu is None: insitu = rcParams['insitu']
3290 self._math._setinsitu(insitu)
3291 varlist = vars()
3292 s = None
3293 import numpy
3294 if isinstance(factor, list) or isinstance(factor, numpy.ndarray):
3295 if isinstance(factor[0], list) or isinstance(factor[0],
3296 numpy.ndarray):
3297 from asapmath import _array2dOp
3298 s = _array2dOp( self, factor, "MUL", tsys, insitu )
3299 else:
3300 s = scantable( self._math._arrayop( self, factor,
3301 "MUL", tsys ) )
3302 else:
3303 s = scantable(self._math._unaryop(self, factor, "MUL", tsys))
3304 s._add_history("scale", varlist)
3305 if insitu:
3306 self._assign(s)
3307 else:
3308 return s
3309
3310 @preserve_selection
3311 def set_sourcetype(self, match, matchtype="pattern",
3312 sourcetype="reference"):
3313 """\
3314 Set the type of the source to be an source or reference scan
3315 using the provided pattern.
3316
3317 Parameters:
3318
3319 match: a Unix style pattern, regular expression or selector
3320
3321 matchtype: 'pattern' (default) UNIX style pattern or
3322 'regex' regular expression
3323
3324 sourcetype: the type of the source to use (source/reference)
3325
3326 """
3327 varlist = vars()
3328 stype = -1
3329 if sourcetype.lower().startswith("r") or sourcetype.lower() == "off":
3330 stype = 1
3331 elif sourcetype.lower().startswith("s") or sourcetype.lower() == "on":
3332 stype = 0
3333 else:
3334 raise ValueError("Illegal sourcetype use s(ource)/on or r(eference)/off")
3335 if matchtype.lower().startswith("p"):
3336 matchtype = "pattern"
3337 elif matchtype.lower().startswith("r"):
3338 matchtype = "regex"
3339 else:
3340 raise ValueError("Illegal matchtype, use p(attern) or r(egex)")
3341 sel = selector()
3342 if isinstance(match, selector):
3343 sel = match
3344 else:
3345 sel.set_query("SRCNAME=%s('%s')" % (matchtype, match))
3346 self.set_selection(sel)
3347 self._setsourcetype(stype)
3348 self._add_history("set_sourcetype", varlist)
3349
3350 @asaplog_post_dec
3351 @preserve_selection
3352 def auto_quotient(self, preserve=True, mode='paired', verify=False):
3353 """\
3354 This function allows to build quotients automatically.
3355 It assumes the observation to have the same number of
3356 "ons" and "offs"
3357
3358 Parameters:
3359
3360 preserve: you can preserve (default) the continuum or
3361 remove it. The equations used are
3362
3363 preserve: Output = Toff * (on/off) - Toff
3364
3365 remove: Output = Toff * (on/off) - Ton
3366
3367 mode: the on/off detection mode
3368 'paired' (default)
3369 identifies 'off' scans by the
3370 trailing '_R' (Mopra/Parkes) or
3371 '_e'/'_w' (Tid) and matches
3372 on/off pairs from the observing pattern
3373 'time'
3374 finds the closest off in time
3375
3376 .. todo:: verify argument is not implemented
3377
3378 """
3379 varlist = vars()
3380 modes = ["time", "paired"]
3381 if not mode in modes:
3382 msg = "please provide valid mode. Valid modes are %s" % (modes)
3383 raise ValueError(msg)
3384 s = None
3385 if mode.lower() == "paired":
3386 sel = self.get_selection()
3387 sel.set_query("SRCTYPE==psoff")
3388 self.set_selection(sel)
3389 offs = self.copy()
3390 sel.set_query("SRCTYPE==pson")
3391 self.set_selection(sel)
3392 ons = self.copy()
3393 s = scantable(self._math._quotient(ons, offs, preserve))
3394 elif mode.lower() == "time":
3395 s = scantable(self._math._auto_quotient(self, mode, preserve))
3396 s._add_history("auto_quotient", varlist)
3397 return s
3398
3399 @asaplog_post_dec
3400 def mx_quotient(self, mask = None, weight='median', preserve=True):
3401 """\
3402 Form a quotient using "off" beams when observing in "MX" mode.
3403
3404 Parameters:
3405
3406 mask: an optional mask to be used when weight == 'stddev'
3407
3408 weight: How to average the off beams. Default is 'median'.
3409
3410 preserve: you can preserve (default) the continuum or
3411 remove it. The equations used are:
3412
3413 preserve: Output = Toff * (on/off) - Toff
3414
3415 remove: Output = Toff * (on/off) - Ton
3416
3417 """
3418 mask = mask or ()
3419 varlist = vars()
3420 on = scantable(self._math._mx_extract(self, 'on'))
3421 preoff = scantable(self._math._mx_extract(self, 'off'))
3422 off = preoff.average_time(mask=mask, weight=weight, scanav=False)
3423 from asapmath import quotient
3424 q = quotient(on, off, preserve)
3425 q._add_history("mx_quotient", varlist)
3426 return q
3427
3428 @asaplog_post_dec
3429 def freq_switch(self, insitu=None):
3430 """\
3431 Apply frequency switching to the data.
3432
3433 Parameters:
3434
3435 insitu: if False a new scantable is returned.
3436 Otherwise, the swictching is done in-situ
3437 The default is taken from .asaprc (False)
3438
3439 """
3440 if insitu is None: insitu = rcParams['insitu']
3441 self._math._setinsitu(insitu)
3442 varlist = vars()
3443 s = scantable(self._math._freqswitch(self))
3444 s._add_history("freq_switch", varlist)
3445 if insitu:
3446 self._assign(s)
3447 else:
3448 return s
3449
3450 @asaplog_post_dec
3451 def recalc_azel(self):
3452 """Recalculate the azimuth and elevation for each position."""
3453 varlist = vars()
3454 self._recalcazel()
3455 self._add_history("recalc_azel", varlist)
3456 return
3457
3458 @asaplog_post_dec
3459 def __add__(self, other):
3460 """
3461 implicit on all axes and on Tsys
3462 """
3463 varlist = vars()
3464 s = self.__op( other, "ADD" )
3465 s._add_history("operator +", varlist)
3466 return s
3467
3468 @asaplog_post_dec
3469 def __sub__(self, other):
3470 """
3471 implicit on all axes and on Tsys
3472 """
3473 varlist = vars()
3474 s = self.__op( other, "SUB" )
3475 s._add_history("operator -", varlist)
3476 return s
3477
3478 @asaplog_post_dec
3479 def __mul__(self, other):
3480 """
3481 implicit on all axes and on Tsys
3482 """
3483 varlist = vars()
3484 s = self.__op( other, "MUL" ) ;
3485 s._add_history("operator *", varlist)
3486 return s
3487
3488
3489 @asaplog_post_dec
3490 def __div__(self, other):
3491 """
3492 implicit on all axes and on Tsys
3493 """
3494 varlist = vars()
3495 s = self.__op( other, "DIV" )
3496 s._add_history("operator /", varlist)
3497 return s
3498
3499 @asaplog_post_dec
3500 def __op( self, other, mode ):
3501 s = None
3502 if isinstance(other, scantable):
3503 s = scantable(self._math._binaryop(self, other, mode))
3504 elif isinstance(other, float):
3505 if other == 0.0:
3506 raise ZeroDivisionError("Dividing by zero is not recommended")
3507 s = scantable(self._math._unaryop(self, other, mode, False))
3508 elif isinstance(other, list) or isinstance(other, numpy.ndarray):
3509 if isinstance(other[0], list) \
3510 or isinstance(other[0], numpy.ndarray):
3511 from asapmath import _array2dOp
3512 s = _array2dOp( self, other, mode, False )
3513 else:
3514 s = scantable( self._math._arrayop( self, other,
3515 mode, False ) )
3516 else:
3517 raise TypeError("Other input is not a scantable or float value")
3518 return s
3519
3520 @asaplog_post_dec
3521 def get_fit(self, row=0):
3522 """\
3523 Print or return the stored fits for a row in the scantable
3524
3525 Parameters:
3526
3527 row: the row which the fit has been applied to.
3528
3529 """
3530 if row > self.nrow():
3531 return
3532 from asap.asapfit import asapfit
3533 fit = asapfit(self._getfit(row))
3534 asaplog.push( '%s' %(fit) )
3535 return fit.as_dict()
3536
3537 @preserve_selection
3538 def flag_nans(self):
3539 """\
3540 Utility function to flag NaN values in the scantable.
3541 """
3542 import numpy
3543 basesel = self.get_selection()
3544 for i in range(self.nrow()):
3545 sel = self.get_row_selector(i)
3546 self.set_selection(basesel+sel)
3547 nans = numpy.isnan(self._getspectrum(0))
3548 if numpy.any(nans):
3549 bnans = [ bool(v) for v in nans]
3550 self.flag(bnans)
3551
3552 def get_row_selector(self, rowno):
3553 return selector(rows=[rowno])
3554
3555 def _add_history(self, funcname, parameters):
3556 if not rcParams['scantable.history']:
3557 return
3558 # create date
3559 sep = "##"
3560 from datetime import datetime
3561 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
3562 hist = dstr+sep
3563 hist += funcname+sep#cdate+sep
3564 if parameters.has_key('self'):
3565 del parameters['self']
3566 for k, v in parameters.iteritems():
3567 if type(v) is dict:
3568 for k2, v2 in v.iteritems():
3569 hist += k2
3570 hist += "="
3571 if isinstance(v2, scantable):
3572 hist += 'scantable'
3573 elif k2 == 'mask':
3574 if isinstance(v2, list) or isinstance(v2, tuple):
3575 hist += str(self._zip_mask(v2))
3576 else:
3577 hist += str(v2)
3578 else:
3579 hist += str(v2)
3580 else:
3581 hist += k
3582 hist += "="
3583 if isinstance(v, scantable):
3584 hist += 'scantable'
3585 elif k == 'mask':
3586 if isinstance(v, list) or isinstance(v, tuple):
3587 hist += str(self._zip_mask(v))
3588 else:
3589 hist += str(v)
3590 else:
3591 hist += str(v)
3592 hist += sep
3593 hist = hist[:-2] # remove trailing '##'
3594 self._addhistory(hist)
3595
3596
3597 def _zip_mask(self, mask):
3598 mask = list(mask)
3599 i = 0
3600 segments = []
3601 while mask[i:].count(1):
3602 i += mask[i:].index(1)
3603 if mask[i:].count(0):
3604 j = i + mask[i:].index(0)
3605 else:
3606 j = len(mask)
3607 segments.append([i, j])
3608 i = j
3609 return segments
3610
3611 def _get_ordinate_label(self):
3612 fu = "("+self.get_fluxunit()+")"
3613 import re
3614 lbl = "Intensity"
3615 if re.match(".K.", fu):
3616 lbl = "Brightness Temperature "+ fu
3617 elif re.match(".Jy.", fu):
3618 lbl = "Flux density "+ fu
3619 return lbl
3620
3621 def _check_ifs(self):
3622# return len(set([self.nchan(i) for i in self.getifnos()])) == 1
3623 nchans = [self.nchan(i) for i in self.getifnos()]
3624 nchans = filter(lambda t: t > 0, nchans)
3625 return (sum(nchans)/len(nchans) == nchans[0])
3626
3627 @asaplog_post_dec
3628 def _fill(self, names, unit, average, opts={}):
3629 first = True
3630 fullnames = []
3631 for name in names:
3632 name = os.path.expandvars(name)
3633 name = os.path.expanduser(name)
3634 if not os.path.exists(name):
3635 msg = "File '%s' does not exists" % (name)
3636 raise IOError(msg)
3637 fullnames.append(name)
3638 if average:
3639 asaplog.push('Auto averaging integrations')
3640 stype = int(rcParams['scantable.storage'].lower() == 'disk')
3641 for name in fullnames:
3642 tbl = Scantable(stype)
3643 if is_ms( name ):
3644 r = msfiller( tbl )
3645 else:
3646 r = filler( tbl )
3647 msg = "Importing %s..." % (name)
3648 asaplog.push(msg, False)
3649 r.open(name, opts)
3650 rx = rcParams['scantable.reference']
3651 r.setreferenceexpr(rx)
3652 r.fill()
3653 if average:
3654 tbl = self._math._average((tbl, ), (), 'NONE', 'SCAN')
3655 if not first:
3656 tbl = self._math._merge([self, tbl])
3657 Scantable.__init__(self, tbl)
3658 r.close()
3659 del r, tbl
3660 first = False
3661 #flush log
3662 asaplog.post()
3663 if unit is not None:
3664 self.set_fluxunit(unit)
3665 if not is_casapy():
3666 self.set_freqframe(rcParams['scantable.freqframe'])
3667
3668
3669 def __getitem__(self, key):
3670 if key < 0:
3671 key += self.nrow()
3672 if key >= self.nrow():
3673 raise IndexError("Row index out of range.")
3674 return self._getspectrum(key)
3675
3676 def __setitem__(self, key, value):
3677 if key < 0:
3678 key += self.nrow()
3679 if key >= self.nrow():
3680 raise IndexError("Row index out of range.")
3681 if not hasattr(value, "__len__") or \
3682 len(value) > self.nchan(self.getif(key)):
3683 raise ValueError("Spectrum length doesn't match.")
3684 return self._setspectrum(value, key)
3685
3686 def __len__(self):
3687 return self.nrow()
3688
3689 def __iter__(self):
3690 for i in range(len(self)):
3691 yield self[i]
Note: See TracBrowser for help on using the repository browser.