source: trunk/python/scantable.py@ 2669

Last change on this file since 2669 was 2645, checked in by WataruKawasaki, 12 years ago

New Development: Yes

JIRA Issue: Yes CAS-4145

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: Yes

Module(s): scantable

Description: added scantable methods [auto_]chebyshev_baseline() to subtract baseline using Chebyshev polynomials.


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