source: trunk/python/scantable.py@ 2540

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

New Development: No

JIRA Issue: Yes (CAS-3749/TRAC-266)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

interactive tests are required

  1. on CASA

(1) plot something by CASA sdplot task
(2) run CASA sdfit with plotlevel!=0
(3) close plotter window by pressing X button on the window
(4) plot something by sdplot task again.

  1. on stand-alone ASAP

(1) plot something by asapplotter
(2) run asapfitter.auto_fit with plot=T
(3) plot something by asapplotter again

Put in Release Notes: No

Module(s): sdplot, asapplotter

Description:

Fixed a bug which caused sdplot and asapplotter crash after closing
plotter loaded by the other modules (tasks) with X button. It was
a bug caused by the fact ID numbers of figure managers are
identical in both sdplot (asapplotter) and the other plotter.
It is also possible, that user loads plotter with the ID identical
to asapplotter and overrides it. This is because, user can select
arbitrary ID when generating a new plot.

asap.plotter._assert_plotter does more tests to make sure the
previous plotter GUI is alive, and properly reloads plotter if not.

Also some minor fixes to reset linewidth for plot invoked by
modules other than asapplot/sdplot, i.e., scantable, asapfitter,
and asapmath.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 142.0 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 f.mask = mask_and(mask, workscan._getmask(r)) # (CAS-1434)
2967 f.data = None
2968 f.fit()
2969
2970 f.plot(residual=True)
2971 accept_fit = raw_input("Accept fit ( [y]/n ): ")
2972 if accept_fit.upper() == "N":
2973 #workscan._append_blinfo(None, None, None)
2974 continue
2975
2976 blpars = f.get_parameters()
2977 masklist = workscan.get_masklist(f.mask, row=r, silent=True)
2978 #workscan._append_blinfo(blpars, masklist, f.mask)
2979 workscan._setspectrum((f.fitter.getresidual()
2980 if getresidual else f.fitter.getfit()), r)
2981
2982 if outblfile:
2983 rms = workscan.get_rms(f.mask, r)
2984 dataout = \
2985 workscan.format_blparams_row(blpars["params"],
2986 blpars["fixed"],
2987 rms, str(masklist),
2988 r, True)
2989 blf.write(dataout)
2990
2991 f._p.unmap()
2992 f._p = None
2993
2994 if outblfile:
2995 blf.close()
2996 else:
2997 workscan._poly_baseline(mask, order, getresidual,
2998 pack_progress_params(showprogress,
2999 minnrow),
3000 outlog, blfile)
3001
3002 workscan._add_history("poly_baseline", varlist)
3003
3004 if insitu:
3005 self._assign(workscan)
3006 else:
3007 return workscan
3008
3009 except RuntimeError, e:
3010 raise_fitting_failure_exception(e)
3011
3012 @asaplog_post_dec
3013 def auto_poly_baseline(self, mask=None, order=None, edge=None,
3014 threshold=None, chan_avg_limit=None,
3015 plot=None, insitu=None,
3016 getresidual=None, showprogress=None,
3017 minnrow=None, outlog=None, blfile=None):
3018 """\
3019 Return a scan which has been baselined (all rows) by a polynomial.
3020 Spectral lines are detected first using linefinder and masked out
3021 to avoid them affecting the baseline solution.
3022
3023 Parameters:
3024 insitu: if False a new scantable is returned.
3025 Otherwise, the scaling is done in-situ
3026 The default is taken from .asaprc (False)
3027 mask: an optional mask retreived from scantable
3028 order: the order of the polynomial (default is 0)
3029 edge: an optional number of channel to drop at
3030 the edge of spectrum. If only one value is
3031 specified, the same number will be dropped
3032 from both sides of the spectrum. Default
3033 is to keep all channels. Nested tuples
3034 represent individual edge selection for
3035 different IFs (a number of spectral channels
3036 can be different)
3037 threshold: the threshold used by line finder. It is
3038 better to keep it large as only strong lines
3039 affect the baseline solution.
3040 chan_avg_limit: a maximum number of consequtive spectral
3041 channels to average during the search of
3042 weak and broad lines. The default is no
3043 averaging (and no search for weak lines).
3044 If such lines can affect the fitted baseline
3045 (e.g. a high order polynomial is fitted),
3046 increase this parameter (usually values up
3047 to 8 are reasonable). Most users of this
3048 method should find the default value sufficient.
3049 plot: plot the fit and the residual. In this each
3050 indivual fit has to be approved, by typing 'y'
3051 or 'n'
3052 getresidual: if False, returns best-fit values instead of
3053 residual. (default is True)
3054 showprogress: show progress status for large data.
3055 default is True.
3056 minnrow: minimum number of input spectra to show.
3057 default is 1000.
3058 outlog: Output the coefficients of the best-fit
3059 function to logger (default is False)
3060 blfile: Name of a text file in which the best-fit
3061 parameter values to be written
3062 (default is "": no file/logger output)
3063
3064 Example:
3065 bscan = scan.auto_poly_baseline(order=7, insitu=False)
3066 """
3067
3068 try:
3069 varlist = vars()
3070
3071 if insitu is None:
3072 insitu = rcParams['insitu']
3073 if insitu:
3074 workscan = self
3075 else:
3076 workscan = self.copy()
3077
3078 #if mask is None: mask = [True for i in xrange(workscan.nchan())]
3079 if mask is None: mask = []
3080 if order is None: order = 0
3081 if edge is None: edge = (0, 0)
3082 if threshold is None: threshold = 3
3083 if chan_avg_limit is None: chan_avg_limit = 1
3084 if plot is None: plot = False
3085 if getresidual is None: getresidual = True
3086 if showprogress is None: showprogress = True
3087 if minnrow is None: minnrow = 1000
3088 if outlog is None: outlog = False
3089 if blfile is None: blfile = ''
3090
3091 edge = normalise_edge_param(edge)
3092
3093 if plot:
3094 outblfile = (blfile != "") and \
3095 os.path.exists(os.path.expanduser(os.path.expandvars(blfile)))
3096 if outblfile: blf = open(blfile, "a")
3097
3098 from asap.asaplinefind import linefinder
3099 fl = linefinder()
3100 fl.set_options(threshold=threshold, avg_limit=chan_avg_limit)
3101 fl.set_scan(workscan)
3102
3103 f = fitter()
3104 f.set_function(lpoly=order)
3105
3106 rows = xrange(workscan.nrow())
3107 #if len(rows) > 0: workscan._init_blinfo()
3108
3109 for r in rows:
3110 idx = 2*workscan.getif(r)
3111 fl.find_lines(r, mask_and(mask, workscan._getmask(r)),
3112 edge[idx:idx+2]) # (CAS-1434)
3113
3114 f.x = workscan._getabcissa(r)
3115 f.y = workscan._getspectrum(r)
3116 f.mask = fl.get_mask()
3117 f.data = None
3118 f.fit()
3119
3120 f.plot(residual=True)
3121 accept_fit = raw_input("Accept fit ( [y]/n ): ")
3122 if accept_fit.upper() == "N":
3123 #workscan._append_blinfo(None, None, None)
3124 continue
3125
3126 blpars = f.get_parameters()
3127 masklist = workscan.get_masklist(f.mask, row=r, silent=True)
3128 #workscan._append_blinfo(blpars, masklist, f.mask)
3129 workscan._setspectrum(
3130 (f.fitter.getresidual() if getresidual
3131 else f.fitter.getfit()), r
3132 )
3133
3134 if outblfile:
3135 rms = workscan.get_rms(f.mask, r)
3136 dataout = \
3137 workscan.format_blparams_row(blpars["params"],
3138 blpars["fixed"],
3139 rms, str(masklist),
3140 r, True)
3141 blf.write(dataout)
3142
3143 f._p.unmap()
3144 f._p = None
3145
3146 if outblfile: blf.close()
3147 else:
3148 workscan._auto_poly_baseline(mask, order, edge, threshold,
3149 chan_avg_limit, getresidual,
3150 pack_progress_params(showprogress,
3151 minnrow),
3152 outlog, blfile)
3153
3154 workscan._add_history("auto_poly_baseline", varlist)
3155
3156 if insitu:
3157 self._assign(workscan)
3158 else:
3159 return workscan
3160
3161 except RuntimeError, e:
3162 raise_fitting_failure_exception(e)
3163
3164 def _init_blinfo(self):
3165 """\
3166 Initialise the following three auxiliary members:
3167 blpars : parameters of the best-fit baseline,
3168 masklists : mask data (edge positions of masked channels) and
3169 actualmask : mask data (in boolean list),
3170 to keep for use later (including output to logger/text files).
3171 Used by poly_baseline() and auto_poly_baseline() in case of
3172 'plot=True'.
3173 """
3174 self.blpars = []
3175 self.masklists = []
3176 self.actualmask = []
3177 return
3178
3179 def _append_blinfo(self, data_blpars, data_masklists, data_actualmask):
3180 """\
3181 Append baseline-fitting related info to blpars, masklist and
3182 actualmask.
3183 """
3184 self.blpars.append(data_blpars)
3185 self.masklists.append(data_masklists)
3186 self.actualmask.append(data_actualmask)
3187 return
3188
3189 @asaplog_post_dec
3190 def rotate_linpolphase(self, angle):
3191 """\
3192 Rotate the phase of the complex polarization O=Q+iU correlation.
3193 This is always done in situ in the raw data. So if you call this
3194 function more than once then each call rotates the phase further.
3195
3196 Parameters:
3197
3198 angle: The angle (degrees) to rotate (add) by.
3199
3200 Example::
3201
3202 scan.rotate_linpolphase(2.3)
3203
3204 """
3205 varlist = vars()
3206 self._math._rotate_linpolphase(self, angle)
3207 self._add_history("rotate_linpolphase", varlist)
3208 return
3209
3210 @asaplog_post_dec
3211 def rotate_xyphase(self, angle):
3212 """\
3213 Rotate the phase of the XY correlation. This is always done in situ
3214 in the data. So if you call this function more than once
3215 then each call rotates the phase further.
3216
3217 Parameters:
3218
3219 angle: The angle (degrees) to rotate (add) by.
3220
3221 Example::
3222
3223 scan.rotate_xyphase(2.3)
3224
3225 """
3226 varlist = vars()
3227 self._math._rotate_xyphase(self, angle)
3228 self._add_history("rotate_xyphase", varlist)
3229 return
3230
3231 @asaplog_post_dec
3232 def swap_linears(self):
3233 """\
3234 Swap the linear polarisations XX and YY, or better the first two
3235 polarisations as this also works for ciculars.
3236 """
3237 varlist = vars()
3238 self._math._swap_linears(self)
3239 self._add_history("swap_linears", varlist)
3240 return
3241
3242 @asaplog_post_dec
3243 def invert_phase(self):
3244 """\
3245 Invert the phase of the complex polarisation
3246 """
3247 varlist = vars()
3248 self._math._invert_phase(self)
3249 self._add_history("invert_phase", varlist)
3250 return
3251
3252 @asaplog_post_dec
3253 def add(self, offset, insitu=None):
3254 """\
3255 Return a scan where all spectra have the offset added
3256
3257 Parameters:
3258
3259 offset: the offset
3260
3261 insitu: if False a new scantable is returned.
3262 Otherwise, the scaling is done in-situ
3263 The default is taken from .asaprc (False)
3264
3265 """
3266 if insitu is None: insitu = rcParams['insitu']
3267 self._math._setinsitu(insitu)
3268 varlist = vars()
3269 s = scantable(self._math._unaryop(self, offset, "ADD", False))
3270 s._add_history("add", varlist)
3271 if insitu:
3272 self._assign(s)
3273 else:
3274 return s
3275
3276 @asaplog_post_dec
3277 def scale(self, factor, tsys=True, insitu=None):
3278 """\
3279
3280 Return a scan where all spectra are scaled by the given 'factor'
3281
3282 Parameters:
3283
3284 factor: the scaling factor (float or 1D float list)
3285
3286 insitu: if False a new scantable is returned.
3287 Otherwise, the scaling is done in-situ
3288 The default is taken from .asaprc (False)
3289
3290 tsys: if True (default) then apply the operation to Tsys
3291 as well as the data
3292
3293 """
3294 if insitu is None: insitu = rcParams['insitu']
3295 self._math._setinsitu(insitu)
3296 varlist = vars()
3297 s = None
3298 import numpy
3299 if isinstance(factor, list) or isinstance(factor, numpy.ndarray):
3300 if isinstance(factor[0], list) or isinstance(factor[0],
3301 numpy.ndarray):
3302 from asapmath import _array2dOp
3303 s = _array2dOp( self, factor, "MUL", tsys, insitu )
3304 else:
3305 s = scantable( self._math._arrayop( self, factor,
3306 "MUL", tsys ) )
3307 else:
3308 s = scantable(self._math._unaryop(self, factor, "MUL", tsys))
3309 s._add_history("scale", varlist)
3310 if insitu:
3311 self._assign(s)
3312 else:
3313 return s
3314
3315 @preserve_selection
3316 def set_sourcetype(self, match, matchtype="pattern",
3317 sourcetype="reference"):
3318 """\
3319 Set the type of the source to be an source or reference scan
3320 using the provided pattern.
3321
3322 Parameters:
3323
3324 match: a Unix style pattern, regular expression or selector
3325
3326 matchtype: 'pattern' (default) UNIX style pattern or
3327 'regex' regular expression
3328
3329 sourcetype: the type of the source to use (source/reference)
3330
3331 """
3332 varlist = vars()
3333 stype = -1
3334 if sourcetype.lower().startswith("r") or sourcetype.lower() == "off":
3335 stype = 1
3336 elif sourcetype.lower().startswith("s") or sourcetype.lower() == "on":
3337 stype = 0
3338 else:
3339 raise ValueError("Illegal sourcetype use s(ource)/on or r(eference)/off")
3340 if matchtype.lower().startswith("p"):
3341 matchtype = "pattern"
3342 elif matchtype.lower().startswith("r"):
3343 matchtype = "regex"
3344 else:
3345 raise ValueError("Illegal matchtype, use p(attern) or r(egex)")
3346 sel = selector()
3347 if isinstance(match, selector):
3348 sel = match
3349 else:
3350 sel.set_query("SRCNAME=%s('%s')" % (matchtype, match))
3351 self.set_selection(sel)
3352 self._setsourcetype(stype)
3353 self._add_history("set_sourcetype", varlist)
3354
3355 @asaplog_post_dec
3356 @preserve_selection
3357 def auto_quotient(self, preserve=True, mode='paired', verify=False):
3358 """\
3359 This function allows to build quotients automatically.
3360 It assumes the observation to have the same number of
3361 "ons" and "offs"
3362
3363 Parameters:
3364
3365 preserve: you can preserve (default) the continuum or
3366 remove it. The equations used are
3367
3368 preserve: Output = Toff * (on/off) - Toff
3369
3370 remove: Output = Toff * (on/off) - Ton
3371
3372 mode: the on/off detection mode
3373 'paired' (default)
3374 identifies 'off' scans by the
3375 trailing '_R' (Mopra/Parkes) or
3376 '_e'/'_w' (Tid) and matches
3377 on/off pairs from the observing pattern
3378 'time'
3379 finds the closest off in time
3380
3381 .. todo:: verify argument is not implemented
3382
3383 """
3384 varlist = vars()
3385 modes = ["time", "paired"]
3386 if not mode in modes:
3387 msg = "please provide valid mode. Valid modes are %s" % (modes)
3388 raise ValueError(msg)
3389 s = None
3390 if mode.lower() == "paired":
3391 sel = self.get_selection()
3392 sel.set_query("SRCTYPE==psoff")
3393 self.set_selection(sel)
3394 offs = self.copy()
3395 sel.set_query("SRCTYPE==pson")
3396 self.set_selection(sel)
3397 ons = self.copy()
3398 s = scantable(self._math._quotient(ons, offs, preserve))
3399 elif mode.lower() == "time":
3400 s = scantable(self._math._auto_quotient(self, mode, preserve))
3401 s._add_history("auto_quotient", varlist)
3402 return s
3403
3404 @asaplog_post_dec
3405 def mx_quotient(self, mask = None, weight='median', preserve=True):
3406 """\
3407 Form a quotient using "off" beams when observing in "MX" mode.
3408
3409 Parameters:
3410
3411 mask: an optional mask to be used when weight == 'stddev'
3412
3413 weight: How to average the off beams. Default is 'median'.
3414
3415 preserve: you can preserve (default) the continuum or
3416 remove it. The equations used are:
3417
3418 preserve: Output = Toff * (on/off) - Toff
3419
3420 remove: Output = Toff * (on/off) - Ton
3421
3422 """
3423 mask = mask or ()
3424 varlist = vars()
3425 on = scantable(self._math._mx_extract(self, 'on'))
3426 preoff = scantable(self._math._mx_extract(self, 'off'))
3427 off = preoff.average_time(mask=mask, weight=weight, scanav=False)
3428 from asapmath import quotient
3429 q = quotient(on, off, preserve)
3430 q._add_history("mx_quotient", varlist)
3431 return q
3432
3433 @asaplog_post_dec
3434 def freq_switch(self, insitu=None):
3435 """\
3436 Apply frequency switching to the data.
3437
3438 Parameters:
3439
3440 insitu: if False a new scantable is returned.
3441 Otherwise, the swictching is done in-situ
3442 The default is taken from .asaprc (False)
3443
3444 """
3445 if insitu is None: insitu = rcParams['insitu']
3446 self._math._setinsitu(insitu)
3447 varlist = vars()
3448 s = scantable(self._math._freqswitch(self))
3449 s._add_history("freq_switch", varlist)
3450 if insitu:
3451 self._assign(s)
3452 else:
3453 return s
3454
3455 @asaplog_post_dec
3456 def recalc_azel(self):
3457 """Recalculate the azimuth and elevation for each position."""
3458 varlist = vars()
3459 self._recalcazel()
3460 self._add_history("recalc_azel", varlist)
3461 return
3462
3463 @asaplog_post_dec
3464 def __add__(self, other):
3465 varlist = vars()
3466 s = None
3467 if isinstance(other, scantable):
3468 s = scantable(self._math._binaryop(self, other, "ADD"))
3469 elif isinstance(other, float):
3470 s = scantable(self._math._unaryop(self, other, "ADD", False))
3471 elif isinstance(other, list) or isinstance(other, numpy.ndarray):
3472 if isinstance(other[0], list) \
3473 or isinstance(other[0], numpy.ndarray):
3474 from asapmath import _array2dOp
3475 s = _array2dOp( self.copy(), other, "ADD", False )
3476 else:
3477 s = scantable( self._math._arrayop( self.copy(), other,
3478 "ADD", False ) )
3479 else:
3480 raise TypeError("Other input is not a scantable or float value")
3481 s._add_history("operator +", varlist)
3482 return s
3483
3484 @asaplog_post_dec
3485 def __sub__(self, other):
3486 """
3487 implicit on all axes and on Tsys
3488 """
3489 varlist = vars()
3490 s = None
3491 if isinstance(other, scantable):
3492 s = scantable(self._math._binaryop(self, other, "SUB"))
3493 elif isinstance(other, float):
3494 s = scantable(self._math._unaryop(self, other, "SUB", False))
3495 elif isinstance(other, list) or isinstance(other, numpy.ndarray):
3496 if isinstance(other[0], list) \
3497 or isinstance(other[0], numpy.ndarray):
3498 from asapmath import _array2dOp
3499 s = _array2dOp( self.copy(), other, "SUB", False )
3500 else:
3501 s = scantable( self._math._arrayop( self.copy(), other,
3502 "SUB", False ) )
3503 else:
3504 raise TypeError("Other input is not a scantable or float value")
3505 s._add_history("operator -", varlist)
3506 return s
3507
3508 @asaplog_post_dec
3509 def __mul__(self, other):
3510 """
3511 implicit on all axes and on Tsys
3512 """
3513 varlist = vars()
3514 s = None
3515 if isinstance(other, scantable):
3516 s = scantable(self._math._binaryop(self, other, "MUL"))
3517 elif isinstance(other, float):
3518 s = scantable(self._math._unaryop(self, other, "MUL", False))
3519 elif isinstance(other, list) or isinstance(other, numpy.ndarray):
3520 if isinstance(other[0], list) \
3521 or isinstance(other[0], numpy.ndarray):
3522 from asapmath import _array2dOp
3523 s = _array2dOp( self.copy(), other, "MUL", False )
3524 else:
3525 s = scantable( self._math._arrayop( self.copy(), other,
3526 "MUL", False ) )
3527 else:
3528 raise TypeError("Other input is not a scantable or float value")
3529 s._add_history("operator *", varlist)
3530 return s
3531
3532
3533 @asaplog_post_dec
3534 def __div__(self, other):
3535 """
3536 implicit on all axes and on Tsys
3537 """
3538 varlist = vars()
3539 s = None
3540 if isinstance(other, scantable):
3541 s = scantable(self._math._binaryop(self, other, "DIV"))
3542 elif isinstance(other, float):
3543 if other == 0.0:
3544 raise ZeroDivisionError("Dividing by zero is not recommended")
3545 s = scantable(self._math._unaryop(self, other, "DIV", False))
3546 elif isinstance(other, list) or isinstance(other, numpy.ndarray):
3547 if isinstance(other[0], list) \
3548 or isinstance(other[0], numpy.ndarray):
3549 from asapmath import _array2dOp
3550 s = _array2dOp( self.copy(), other, "DIV", False )
3551 else:
3552 s = scantable( self._math._arrayop( self.copy(), other,
3553 "DIV", False ) )
3554 else:
3555 raise TypeError("Other input is not a scantable or float value")
3556 s._add_history("operator /", varlist)
3557 return s
3558
3559 @asaplog_post_dec
3560 def get_fit(self, row=0):
3561 """\
3562 Print or return the stored fits for a row in the scantable
3563
3564 Parameters:
3565
3566 row: the row which the fit has been applied to.
3567
3568 """
3569 if row > self.nrow():
3570 return
3571 from asap.asapfit import asapfit
3572 fit = asapfit(self._getfit(row))
3573 asaplog.push( '%s' %(fit) )
3574 return fit.as_dict()
3575
3576 @preserve_selection
3577 def flag_nans(self):
3578 """\
3579 Utility function to flag NaN values in the scantable.
3580 """
3581 import numpy
3582 basesel = self.get_selection()
3583 for i in range(self.nrow()):
3584 sel = self.get_row_selector(i)
3585 self.set_selection(basesel+sel)
3586 nans = numpy.isnan(self._getspectrum(0))
3587 if numpy.any(nans):
3588 bnans = [ bool(v) for v in nans]
3589 self.flag(bnans)
3590
3591 def get_row_selector(self, rowno):
3592 return selector(rows=[rowno])
3593
3594 def _add_history(self, funcname, parameters):
3595 if not rcParams['scantable.history']:
3596 return
3597 # create date
3598 sep = "##"
3599 from datetime import datetime
3600 dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
3601 hist = dstr+sep
3602 hist += funcname+sep#cdate+sep
3603 if parameters.has_key('self'):
3604 del parameters['self']
3605 for k, v in parameters.iteritems():
3606 if type(v) is dict:
3607 for k2, v2 in v.iteritems():
3608 hist += k2
3609 hist += "="
3610 if isinstance(v2, scantable):
3611 hist += 'scantable'
3612 elif k2 == 'mask':
3613 if isinstance(v2, list) or isinstance(v2, tuple):
3614 hist += str(self._zip_mask(v2))
3615 else:
3616 hist += str(v2)
3617 else:
3618 hist += str(v2)
3619 else:
3620 hist += k
3621 hist += "="
3622 if isinstance(v, scantable):
3623 hist += 'scantable'
3624 elif k == 'mask':
3625 if isinstance(v, list) or isinstance(v, tuple):
3626 hist += str(self._zip_mask(v))
3627 else:
3628 hist += str(v)
3629 else:
3630 hist += str(v)
3631 hist += sep
3632 hist = hist[:-2] # remove trailing '##'
3633 self._addhistory(hist)
3634
3635
3636 def _zip_mask(self, mask):
3637 mask = list(mask)
3638 i = 0
3639 segments = []
3640 while mask[i:].count(1):
3641 i += mask[i:].index(1)
3642 if mask[i:].count(0):
3643 j = i + mask[i:].index(0)
3644 else:
3645 j = len(mask)
3646 segments.append([i, j])
3647 i = j
3648 return segments
3649
3650 def _get_ordinate_label(self):
3651 fu = "("+self.get_fluxunit()+")"
3652 import re
3653 lbl = "Intensity"
3654 if re.match(".K.", fu):
3655 lbl = "Brightness Temperature "+ fu
3656 elif re.match(".Jy.", fu):
3657 lbl = "Flux density "+ fu
3658 return lbl
3659
3660 def _check_ifs(self):
3661# return len(set([self.nchan(i) for i in self.getifnos()])) == 1
3662 nchans = [self.nchan(i) for i in self.getifnos()]
3663 nchans = filter(lambda t: t > 0, nchans)
3664 return (sum(nchans)/len(nchans) == nchans[0])
3665
3666 @asaplog_post_dec
3667 def _fill(self, names, unit, average, opts={}):
3668 first = True
3669 fullnames = []
3670 for name in names:
3671 name = os.path.expandvars(name)
3672 name = os.path.expanduser(name)
3673 if not os.path.exists(name):
3674 msg = "File '%s' does not exists" % (name)
3675 raise IOError(msg)
3676 fullnames.append(name)
3677 if average:
3678 asaplog.push('Auto averaging integrations')
3679 stype = int(rcParams['scantable.storage'].lower() == 'disk')
3680 for name in fullnames:
3681 tbl = Scantable(stype)
3682 if is_ms( name ):
3683 r = msfiller( tbl )
3684 else:
3685 r = filler( tbl )
3686 msg = "Importing %s..." % (name)
3687 asaplog.push(msg, False)
3688 r.open(name, opts)
3689 rx = rcParams['scantable.reference']
3690 r.setreferenceexpr(rx)
3691 r.fill()
3692 if average:
3693 tbl = self._math._average((tbl, ), (), 'NONE', 'SCAN')
3694 if not first:
3695 tbl = self._math._merge([self, tbl])
3696 Scantable.__init__(self, tbl)
3697 r.close()
3698 del r, tbl
3699 first = False
3700 #flush log
3701 asaplog.post()
3702 if unit is not None:
3703 self.set_fluxunit(unit)
3704 if not is_casapy():
3705 self.set_freqframe(rcParams['scantable.freqframe'])
3706
3707
3708 def __getitem__(self, key):
3709 if key < 0:
3710 key += self.nrow()
3711 if key >= self.nrow():
3712 raise IndexError("Row index out of range.")
3713 return self._getspectrum(key)
3714
3715 def __setitem__(self, key, value):
3716 if key < 0:
3717 key += self.nrow()
3718 if key >= self.nrow():
3719 raise IndexError("Row index out of range.")
3720 if not hasattr(value, "__len__") or \
3721 len(value) > self.nchan(self.getif(key)):
3722 raise ValueError("Spectrum length doesn't match.")
3723 return self._setspectrum(value, key)
3724
3725 def __len__(self):
3726 return self.nrow()
3727
3728 def __iter__(self):
3729 for i in range(len(self)):
3730 yield self[i]
Note: See TracBrowser for help on using the repository browser.