source: trunk/python/scantable.py@ 2558

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

New Development: No

JIRA Issue: Yes (related to CAS-3749)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs: Interactive tests

Test I. run sdbaseline with blfunc='poly', verify=True and masklist=[]

on CASA,
or run scantable.auto_poly_baseline and scantable.poly_baseline
with plot=True and mask=None.
--> Verification plots should properly be loaded.

Test II.
(1) run asapfitter.plot, by for example

scan = asap.scantable(infile,False)
f=asap.fitter()
f.set_function(lpoly=2)
f.x = scan._getabcissa(0)
f.y = scan._getspectrum(0)
f.mask=(scan._getmask(0))
f.data=None
f.fit()
f.plot(residual=True)

(2) plot something with sdplot (on CASA) or asapplotter.plot
(3) run asapfitter.plot again

f.plot(residual=True)
--> (3) should work without an error

Put in Release Notes: No

Module(s):

Description:

[I] Fixed bugs which caused scantable.poly_baseline and
scantable.auto_poly_baseline crash when mask=None (default)
and plot=True.
[II] Fixed a bug which caused an error in asapfitter.plot when the
other plotting operation is invoked between multiple run of asapfitter.plot.


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