source: trunk/python/scantable.py@ 3014

Last change on this file since 3014 was 3010, checked in by WataruKawasaki, 10 years ago

New Development: No

JIRA Issue: Yes CAS-6599

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): sd

Description: modified so that sd.scantable.stats() returns None for flagged rows and rows with all channels flagged.


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