source: trunk/python/scantable.py@ 2965

Last change on this file since 2965 was 2957, checked in by WataruKawasaki, 10 years ago

New Development: No

JIRA Issue: Yes CAS-6599

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: function paramter

Test Programs:

Put in Release Notes: No

Module(s): sd

Description: add a paramter to sd.scantable.stats() to skip outputting formatted statistical values to temporary text file.


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