source: branches/hpc34/python/scantable.py@ 3030

Last change on this file since 3030 was 2640, checked in by WataruKawasaki, 12 years ago

added a new parameter 'csvformat' to sd.scantable.*baseline() and the relevant functions in the C++ side. (2012/08/10 WK)

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