source: trunk/python/scantable.py@ 2259

Last change on this file since 2259 was 2189, checked in by WataruKawasaki, 13 years ago

New Development: No

JIRA Issue: Yes CAS-3149

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: scantable.*_baseline() parameter

Test Programs:

Put in Release Notes: No

Module(s):

Description: Added two parameters 'showprogress' and 'minnrow' to scantable.*_baseline() to enable to show progress status during time-consuming processes.


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