source: trunk/python/scantable.py@ 2765

Last change on this file since 2765 was 2761, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Added new option, freqref, to NRO filler. Posible values are:
1) 'rest' to import frequency in REST frame, which results in an exactly
same frequency label as NEWSTAR, and 2) 'vref' to import frequency
in the frame that source velocity refers, which results in the same
velocity label as NEWSTAR. The option must be given to scantable
constructor.


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