source: branches/casa-release-4_3-test02/python/scantable.py@ 3125

Last change on this file since 3125 was 2987, checked in by WataruKawasaki, 10 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes:

Module(s): sd

Description: fixed parse_idx_selection for mode='row' so that it returns a list of existing row numbers only.


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