Context Navigation

source: trunk/python/scantable.py@ 1858

Visit:

Last change on this file since 1858 was 1857, checked in by Malte Marquarding, 15 years ago
auto_quotient is the first method to use @preserve_selection
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 90.9 KB

Line
1	"""This module defines the scantable class."""
2
3	import os
4	try:
5	from functools import wraps as wraps_dec
6	except ImportError:
7	from asap.compatibility import wraps as wraps_dec
8
9	from asap.env import is_casapy
10	from asap._asap import Scantable
11	from asap._asap import filler
12	from asap.parameters import rcParams
13	from asap.logging import asaplog, print_log, print_log_dec
14	from asap.selector import selector
15	from asap.linecatalog import linecatalog
16	from asap.coordinate import coordinate
17	from asap.utils import _n_bools, mask_not, mask_and, mask_or
18
19
20	def preserve_selection(func):
21	@wraps_dec(func)
22	def wrap(obj, args, *kw):
23	basesel = obj.get_selection()
24	try:
25	val = func(obj, args, *kw)
26	finally:
27	obj.set_selection(basesel)
28	return val
29	return wrap
30
31	def is_scantable(filename):
32	"""Is the given file a scantable?
33
34	Parameters:
35
36	filename: the name of the file/directory to test
37
38	"""
39	return (os.path.isdir(filename)
40	and not os.path.exists(filename+'/table.f1')
41	and os.path.exists(filename+'/table.info'))
42
43
44	class scantable(Scantable):
45	"""\
46	The ASAP container for scans (single-dish data).
47	"""
48
49	@print_log_dec
50	def __init__(self, filename, average=None, unit=None, getpt=None,
51	antenna=None, parallactify=None):
52	"""\
53	Create a scantable from a saved one or make a reference
54
55	Parameters:
56
57	filename: the name of an asap table on disk
58	or
59	the name of a rpfits/sdfits/ms file
60	(integrations within scans are auto averaged
61	and the whole file is read) or
62	[advanced] a reference to an existing scantable
63
64	average: average all integrations withinb a scan on read.
65	The default (True) is taken from .asaprc.
66
67	unit: brightness unit; must be consistent with K or Jy.
68	Over-rides the default selected by the filler
69	(input rpfits/sdfits/ms) or replaces the value
70	in existing scantables
71
72	getpt: for MeasurementSet input data only:
73	If True, all pointing data are filled.
74	The deafult is False, which makes time to load
75	the MS data faster in some cases.
76
77	antenna: Antenna selection. integer (id) or string (name or id).
78
79	parallactify: Indicate that the data had been parallatified. Default
80	is taken from rc file.
81
82	"""
83	if average is None:
84	average = rcParams['scantable.autoaverage']
85	if getpt is None:
86	getpt = True
87	if antenna is not None:
88	asaplog.push("Antenna selection currently unsupported."
89	"Using '0'")
90	print_log('WARN')
91	if antenna is None:
92	antenna = ''
93	elif type(antenna) == int:
94	antenna = '%s' % antenna
95	elif type(antenna) == list:
96	tmpstr = ''
97	for i in range( len(antenna) ):
98	if type(antenna[i]) == int:
99	tmpstr = tmpstr + ('%s,'%(antenna[i]))
100	elif type(antenna[i]) == str:
101	tmpstr=tmpstr+antenna[i]+','
102	else:
103	asaplog.push('Bad antenna selection.')
104	print_log('ERROR')
105	return
106	antenna = tmpstr.rstrip(',')
107	parallactify = parallactify or rcParams['scantable.parallactify']
108	varlist = vars()
109	from asap._asap import stmath
110	self._math = stmath( rcParams['insitu'] )
111	if isinstance(filename, Scantable):
112	Scantable.__init__(self, filename)
113	else:
114	if isinstance(filename, str):
115	filename = os.path.expandvars(filename)
116	filename = os.path.expanduser(filename)
117	if not os.path.exists(filename):
118	s = "File '%s' not found." % (filename)
119	if rcParams['verbose']:
120	asaplog.push(s)
121	print_log('ERROR')
122	return
123	raise IOError(s)
124	if is_scantable(filename):
125	ondisk = rcParams['scantable.storage'] == 'disk'
126	Scantable.__init__(self, filename, ondisk)
127	if unit is not None:
128	self.set_fluxunit(unit)
129	# do not reset to the default freqframe
130	#self.set_freqframe(rcParams['scantable.freqframe'])
131	elif os.path.isdir(filename) \
132	and not os.path.exists(filename+'/table.f1'):
133	msg = "The given file '%s'is not a valid " \
134	"asap table." % (filename)
135	if rcParams['verbose']:
136	#print msg
137	asaplog.push( msg )
138	print_log( 'ERROR' )
139	return
140	else:
141	raise IOError(msg)
142	else:
143	self._fill([filename], unit, average, getpt, antenna)
144	elif (isinstance(filename, list) or isinstance(filename, tuple)) \
145	and isinstance(filename[-1], str):
146	self._fill(filename, unit, average, getpt, antenna)
147	self.parallactify(parallactify)
148	self._add_history("scantable", varlist)
149
150	@print_log_dec
151	def save(self, name=None, format=None, overwrite=False):
152	"""\
153	Store the scantable on disk. This can be an asap (aips++) Table,
154	SDFITS or MS2 format.
155
156	Parameters:
157
158	name: the name of the outputfile. For format "ASCII"
159	this is the root file name (data in 'name'.txt
160	and header in 'name'_header.txt)
161
162	format: an optional file format. Default is ASAP.
163	Allowed are:
164
165	* 'ASAP' (save as ASAP [aips++] Table),
166	* 'SDFITS' (save as SDFITS file)
167	* 'ASCII' (saves as ascii text file)
168	* 'MS2' (saves as an casacore MeasurementSet V2)
169	* 'FITS' (save as image FITS - not readable by class)
170	* 'CLASS' (save as FITS readable by CLASS)
171
172	overwrite: If the file should be overwritten if it exists.
173	The default False is to return with warning
174	without writing the output. USE WITH CARE.
175
176	Example::
177
178	scan.save('myscan.asap')
179	scan.save('myscan.sdfits', 'SDFITS')
180
181	"""
182	from os import path
183	format = format or rcParams['scantable.save']
184	suffix = '.'+format.lower()
185	if name is None or name == "":
186	name = 'scantable'+suffix
187	msg = "No filename given. Using default name %s..." % name
188	asaplog.push(msg)
189	name = path.expandvars(name)
190	if path.isfile(name) or path.isdir(name):
191	if not overwrite:
192	msg = "File %s exists." % name
193	if rcParams['verbose']:
194	#print msg
195	asaplog.push( msg )
196	print_log( 'ERROR' )
197	return
198	else:
199	raise IOError(msg)
200	format2 = format.upper()
201	if format2 == 'ASAP':
202	self._save(name)
203	else:
204	from asap._asap import stwriter as stw
205	writer = stw(format2)
206	writer.write(self, name)
207	return
208
209	def copy(self):
210	"""Return a copy of this scantable.
211
212	Note:
213
214	This makes a full (deep) copy. scan2 = scan1 makes a reference.
215
216	Example::
217
218	copiedscan = scan.copy()
219
220	"""
221	sd = scantable(Scantable._copy(self))
222	return sd
223
224	def drop_scan(self, scanid=None):
225	"""\
226	Return a new scantable where the specified scan number(s) has(have)
227	been dropped.
228
229	Parameters:
230
231	scanid: a (list of) scan number(s)
232
233	"""
234	from asap import _is_sequence_or_number as _is_valid
235	from asap import _to_list
236	from asap import unique
237	if not _is_valid(scanid):
238	if rcParams['verbose']:
239	asaplog.push( 'Please specify a scanno to drop from the scantable' )
240	print_log( 'ERROR' )
241	return
242	else:
243	raise RuntimeError("No scan given")
244	try:
245	scanid = _to_list(scanid)
246	allscans = unique([ self.getscan(i) for i in range(self.nrow())])
247	for sid in scanid: allscans.remove(sid)
248	if len(allscans) == 0:
249	raise ValueError("Can't remove all scans")
250	except ValueError:
251	if rcParams['verbose']:
252	print_log()
253	asaplog.push( "Couldn't find any match." )
254	print_log( 'ERROR' )
255	return
256	else: raise
257	try:
258	sel = selector(scans=allscans)
259	return self._select_copy(sel)
260	except RuntimeError:
261	if rcParams['verbose']:
262	print_log()
263	asaplog.push( "Couldn't find any match." )
264	print_log( 'ERROR' )
265	else:
266	raise
267
268	def _select_copy(self, selection):
269	orig = self.get_selection()
270	self.set_selection(orig+selection)
271	cp = self.copy()
272	self.set_selection(orig)
273	return cp
274
275	def get_scan(self, scanid=None):
276	"""\
277	Return a specific scan (by scanno) or collection of scans (by
278	source name) in a new scantable.
279
280	Note:
281
282	See scantable.drop_scan() for the inverse operation.
283
284	Parameters:
285
286	scanid: a (list of) scanno or a source name, unix-style
287	patterns are accepted for source name matching, e.g.
288	'*_R' gets all 'ref scans
289
290	Example::
291
292	# get all scans containing the source '323p459'
293	newscan = scan.get_scan('323p459')
294	# get all 'off' scans
295	refscans = scan.get_scan('*_R')
296	# get a susbset of scans by scanno (as listed in scan.summary())
297	newscan = scan.get_scan([0, 2, 7, 10])
298
299	"""
300	if scanid is None:
301	if rcParams['verbose']:
302	#print "Please specify a scan no or name to " \
303	# "retrieve from the scantable"
304	asaplog.push( 'Please specify a scan no or name to retrieve'
305	' from the scantable' )
306	print_log( 'ERROR' )
307	return
308	else:
309	raise RuntimeError("No scan given")
310
311	try:
312	bsel = self.get_selection()
313	sel = selector()
314	if type(scanid) is str:
315	sel.set_name(scanid)
316	return self._select_copy(sel)
317	elif type(scanid) is int:
318	sel.set_scans([scanid])
319	return self._select_copy(sel)
320	elif type(scanid) is list:
321	sel.set_scans(scanid)
322	return self._select_copy(sel)
323	else:
324	msg = "Illegal scanid type, use 'int' or 'list' if ints."
325	if rcParams['verbose']:
326	#print msg
327	asaplog.push( msg )
328	print_log( 'ERROR' )
329	else:
330	raise TypeError(msg)
331	except RuntimeError:
332	if rcParams['verbose']:
333	#print "Couldn't find any match."
334	print_log()
335	asaplog.push( "Couldn't find any match." )
336	print_log( 'ERROR' )
337	else: raise
338
339	def __str__(self):
340	return Scantable._summary(self, True)
341
342	def summary(self, filename=None):
343	"""\
344	Print a summary of the contents of this scantable.
345
346	Parameters:
347
348	filename: the name of a file to write the putput to
349	Default - no file output
350
351	"""
352	info = Scantable._summary(self, True)
353	if filename is not None:
354	if filename is "":
355	filename = 'scantable_summary.txt'
356	from os.path import expandvars, isdir
357	filename = expandvars(filename)
358	if not isdir(filename):
359	data = open(filename, 'w')
360	data.write(info)
361	data.close()
362	else:
363	msg = "Illegal file name '%s'." % (filename)
364	if rcParams['verbose']:
365	#print msg
366	asaplog.push( msg )
367	print_log( 'ERROR' )
368	else:
369	raise IOError(msg)
370	if rcParams['verbose']:
371	try:
372	from IPython.genutils import page as pager
373	except ImportError:
374	from pydoc import pager
375	pager(info)
376	else:
377	return info
378
379	def get_spectrum(self, rowno):
380	"""Return the spectrum for the current row in the scantable as a list.
381
382	Parameters:
383
384	rowno: the row number to retrieve the spectrum from
385
386	"""
387	return self._getspectrum(rowno)
388
389	def get_mask(self, rowno):
390	"""Return the mask for the current row in the scantable as a list.
391
392	Parameters:
393
394	rowno: the row number to retrieve the mask from
395
396	"""
397	return self._getmask(rowno)
398
399	def set_spectrum(self, spec, rowno):
400	"""Return the spectrum for the current row in the scantable as a list.
401
402	Parameters:
403
404	spec: the new spectrum
405
406	rowno: the row number to set the spectrum for
407
408	"""
409	assert(len(spec) == self.nchan())
410	return self._setspectrum(spec, rowno)
411
412	def get_coordinate(self, rowno):
413	"""Return the (spectral) coordinate for a a given 'rowno'.
414
415	Note:
416
417	* This coordinate is only valid until a scantable method modifies
418	the frequency axis.
419	* This coordinate does contain the original frequency set-up
420	NOT the new frame. The conversions however are done using the user
421	specified frame (e.g. LSRK/TOPO). To get the 'real' coordinate,
422	use scantable.freq_align first. Without it there is no closure,
423	i.e.::
424
425	c = myscan.get_coordinate(0)
426	c.to_frequency(c.get_reference_pixel()) != c.get_reference_value()
427
428	Parameters:
429
430	rowno: the row number for the spectral coordinate
431
432	"""
433	return coordinate(Scantable.get_coordinate(self, rowno))
434
435	def get_selection(self):
436	"""\
437	Get the selection object currently set on this scantable.
438
439	Example::
440
441	sel = scan.get_selection()
442	sel.set_ifs(0) # select IF 0
443	scan.set_selection(sel) # apply modified selection
444
445	"""
446	return selector(self._getselection())
447
448	def set_selection(self, selection=None, **kw):
449	"""\
450	Select a subset of the data. All following operations on this scantable
451	are only applied to thi selection.
452
453	Parameters:
454
455	selection: a selector object (default unset the selection), or
456	any combination of "pols", "ifs", "beams", "scans",
457	"cycles", "name", "query"
458
459	Examples::
460
461	sel = selector() # create a selection object
462	self.set_scans([0, 3]) # select SCANNO 0 and 3
463	scan.set_selection(sel) # set the selection
464	scan.summary() # will only print summary of scanno 0 an 3
465	scan.set_selection() # unset the selection
466	# or the equivalent
467	scan.set_selection(scans=[0,3])
468	scan.summary() # will only print summary of scanno 0 an 3
469	scan.set_selection() # unset the selection
470
471	"""
472	if selection is None:
473	# reset
474	if len(kw) == 0:
475	selection = selector()
476	else:
477	# try keywords
478	for k in kw:
479	if k not in selector.fields:
480	raise KeyError("Invalid selection key '%s', valid keys are %s" % (k, selector.fields))
481	selection = selector(**kw)
482	self._setselection(selection)
483
484	def get_row(self, row=0, insitu=None):
485	"""\
486	Select a row in the scantable.
487	Return a scantable with single row.
488
489	Parameters:
490
491	row: row no of integration, default is 0.
492	insitu: if False a new scantable is returned. Otherwise, the
493	scaling is done in-situ. The default is taken from .asaprc
494	(False)
495
496	"""
497	if insitu is None: insitu = rcParams['insitu']
498	if not insitu:
499	workscan = self.copy()
500	else:
501	workscan = self
502	# Select a row
503	sel=selector()
504	sel.set_scans([workscan.getscan(row)])
505	sel.set_cycles([workscan.getcycle(row)])
506	sel.set_beams([workscan.getbeam(row)])
507	sel.set_ifs([workscan.getif(row)])
508	sel.set_polarisations([workscan.getpol(row)])
509	sel.set_name(workscan._getsourcename(row))
510	workscan.set_selection(sel)
511	if not workscan.nrow() == 1:
512	msg = "Cloud not identify single row. %d rows selected."%(workscan.nrow())
513	raise RuntimeError(msg)
514	del sel
515	if insitu:
516	self._assign(workscan)
517	else:
518	return workscan
519
520	#def stats(self, stat='stddev', mask=None):
521	def stats(self, stat='stddev', mask=None, form='3.3f'):
522	"""\
523	Determine the specified statistic of the current beam/if/pol
524	Takes a 'mask' as an optional parameter to specify which
525	channels should be excluded.
526
527	Parameters:
528
529	stat: 'min', 'max', 'min_abc', 'max_abc', 'sumsq', 'sum',
530	'mean', 'var', 'stddev', 'avdev', 'rms', 'median'
531
532	mask: an optional mask specifying where the statistic
533	should be determined.
534
535	form: format string to print statistic values
536
537	Example::
538
539	scan.set_unit('channel')
540	msk = scan.create_mask([100, 200], [500, 600])
541	scan.stats(stat='mean', mask=m)
542
543	"""
544	mask = mask or []
545	if not self._check_ifs():
546	raise ValueError("Cannot apply mask as the IFs have different "
547	"number of channels. Please use setselection() "
548	"to select individual IFs")
549	rtnabc = False
550	if stat.lower().endswith('_abc'): rtnabc = True
551	getchan = False
552	if stat.lower().startswith('min') or stat.lower().startswith('max'):
553	chan = self._math._minmaxchan(self, mask, stat)
554	getchan = True
555	statvals = []
556	if not rtnabc: statvals = self._math._stats(self, mask, stat)
557
558	#def cb(i):
559	# return statvals[i]
560
561	#return self._row_callback(cb, stat)
562
563	label=stat
564	#callback=cb
565	out = ""
566	#outvec = []
567	sep = '-'*50
568	for i in range(self.nrow()):
569	refstr = ''
570	statunit= ''
571	if getchan:
572	qx, qy = self.chan2data(rowno=i, chan=chan[i])
573	if rtnabc:
574	statvals.append(qx['value'])
575	refstr = ('(value: %'+form) % (qy['value'])+' ['+qy['unit']+'])'
576	statunit= '['+qx['unit']+']'
577	else:
578	refstr = ('(@ %'+form) % (qx['value'])+' ['+qx['unit']+'])'
579
580	tm = self._gettime(i)
581	src = self._getsourcename(i)
582	out += 'Scan[%d] (%s) ' % (self.getscan(i), src)
583	out += 'Time[%s]:\n' % (tm)
584	if self.nbeam(-1) > 1:
585	out += ' Beam[%d] ' % (self.getbeam(i))
586	if self.nif(-1) > 1: out += ' IF[%d] ' % (self.getif(i))
587	if self.npol(-1) > 1: out += ' Pol[%d] ' % (self.getpol(i))
588	#outvec.append(callback(i))
589	#out += ('= %'+form) % (outvec[i]) +' '+refstr+'\n'
590	out += ('= %'+form) % (statvals[i]) +' '+refstr+'\n'
591	out += sep+"\n"
592
593	if rcParams['verbose']:
594	import os
595	if os.environ.has_key( 'USER' ):
596	usr=os.environ['USER']
597	else:
598	import commands
599	usr=commands.getoutput( 'whoami' )
600	tmpfile='/tmp/tmp_'+usr+'_casapy_asap_scantable_stats'
601	f=open(tmpfile,'w')
602	print >> f, sep
603	print >> f, ' %s %s' % (label, statunit)
604	print >> f, sep
605	print >> f, out
606	f.close()
607	f=open(tmpfile,'r')
608	x=f.readlines()
609	f.close()
610	blanc=''
611	asaplog.push(blanc.join(x), False)
612	#for xx in x:
613	# asaplog.push( xx, False )
614	print_log()
615	return statvals
616
617	def chan2data(self, rowno=0, chan=0):
618	"""\
619	Returns channel/frequency/velocity and spectral value
620	at an arbitrary row and channel in the scantable.
621
622	Parameters:
623
624	rowno: a row number in the scantable. Default is the
625	first row, i.e. rowno=0
626
627	chan: a channel in the scantable. Default is the first
628	channel, i.e. pos=0
629
630	"""
631	if isinstance(rowno, int) and isinstance(chan, int):
632	qx = {'unit': self.get_unit(),
633	'value': self._getabcissa(rowno)[chan]}
634	qy = {'unit': self.get_fluxunit(),
635	'value': self._getspectrum(rowno)[chan]}
636	return qx, qy
637
638	def stddev(self, mask=None):
639	"""\
640	Determine the standard deviation of the current beam/if/pol
641	Takes a 'mask' as an optional parameter to specify which
642	channels should be excluded.
643
644	Parameters:
645
646	mask: an optional mask specifying where the standard
647	deviation should be determined.
648
649	Example::
650
651	scan.set_unit('channel')
652	msk = scan.create_mask([100, 200], [500, 600])
653	scan.stddev(mask=m)
654
655	"""
656	return self.stats(stat='stddev', mask=mask);
657
658
659	def get_column_names(self):
660	"""\
661	Return a list of column names, which can be used for selection.
662	"""
663	return list(Scantable.get_column_names(self))
664
665	def get_tsys(self, row=-1):
666	"""\
667	Return the System temperatures.
668
669	Parameters:
670
671	row: the rowno to get the information for. (default all rows)
672
673	Returns:
674
675	a list of Tsys values for the current selection
676
677	"""
678	if row > -1:
679	return self._get_column(self._gettsys, row)
680	return self._row_callback(self._gettsys, "Tsys")
681
682
683	def get_weather(self, row=-1):
684	"""\
685	Return the weather informations.
686
687	Parameters:
688
689	row: the rowno to get the information for. (default all rows)
690
691	Returns:
692
693	a dict or list of of dicts of values for the current selection
694
695	"""
696
697	values = self._get_column(self._get_weather, row)
698	if row > -1:
699	return {'temperature': values[0],
700	'pressure': values[1], 'humidity' : values[2],
701	'windspeed' : values[3], 'windaz' : values[4]
702	}
703	else:
704	out = []
705	for r in values:
706
707	out.append({'temperature': r[0],
708	'pressure': r[1], 'humidity' : r[2],
709	'windspeed' : r[3], 'windaz' : r[4]
710	})
711	return out
712
713	def _row_callback(self, callback, label):
714	out = ""
715	outvec = []
716	sep = '-'*50
717	for i in range(self.nrow()):
718	tm = self._gettime(i)
719	src = self._getsourcename(i)
720	out += 'Scan[%d] (%s) ' % (self.getscan(i), src)
721	out += 'Time[%s]:\n' % (tm)
722	if self.nbeam(-1) > 1:
723	out += ' Beam[%d] ' % (self.getbeam(i))
724	if self.nif(-1) > 1: out += ' IF[%d] ' % (self.getif(i))
725	if self.npol(-1) > 1: out += ' Pol[%d] ' % (self.getpol(i))
726	outvec.append(callback(i))
727	out += '= %3.3f\n' % (outvec[i])
728	out += sep+'\n'
729	if rcParams['verbose']:
730	asaplog.push(sep)
731	asaplog.push(" %s" % (label))
732	asaplog.push(sep)
733	asaplog.push(out)
734	print_log()
735	return outvec
736
737	def _get_column(self, callback, row=-1):
738	"""
739	"""
740	if row == -1:
741	return [callback(i) for i in range(self.nrow())]
742	else:
743	if 0 <= row < self.nrow():
744	return callback(row)
745
746
747	def get_time(self, row=-1, asdatetime=False):
748	"""\
749	Get a list of time stamps for the observations.
750	Return a datetime object for each integration time stamp in the scantable.
751
752	Parameters:
753
754	row: row no of integration. Default -1 return all rows
755
756	asdatetime: return values as datetime objects rather than strings
757
758	"""
759	from time import strptime
760	from datetime import datetime
761	times = self._get_column(self._gettime, row)
762	if not asdatetime:
763	return times
764	format = "%Y/%m/%d/%H:%M:%S"
765	if isinstance(times, list):
766	return [datetime(*strptime(i, format)[:6]) for i in times]
767	else:
768	return datetime(*strptime(times, format)[:6])
769
770
771	def get_inttime(self, row=-1):
772	"""\
773	Get a list of integration times for the observations.
774	Return a time in seconds for each integration in the scantable.
775
776	Parameters:
777
778	row: row no of integration. Default -1 return all rows.
779
780	"""
781	return self._get_column(self._getinttime, row)
782
783
784	def get_sourcename(self, row=-1):
785	"""\
786	Get a list source names for the observations.
787	Return a string for each integration in the scantable.
788	Parameters:
789
790	row: row no of integration. Default -1 return all rows.
791
792	"""
793	return self._get_column(self._getsourcename, row)
794
795	def get_elevation(self, row=-1):
796	"""\
797	Get a list of elevations for the observations.
798	Return a float for each integration in the scantable.
799
800	Parameters:
801
802	row: row no of integration. Default -1 return all rows.
803
804	"""
805	return self._get_column(self._getelevation, row)
806
807	def get_azimuth(self, row=-1):
808	"""\
809	Get a list of azimuths for the observations.
810	Return a float for each integration in the scantable.
811
812	Parameters:
813	row: row no of integration. Default -1 return all rows.
814
815	"""
816	return self._get_column(self._getazimuth, row)
817
818	def get_parangle(self, row=-1):
819	"""\
820	Get a list of parallactic angles for the observations.
821	Return a float for each integration in the scantable.
822
823	Parameters:
824
825	row: row no of integration. Default -1 return all rows.
826
827	"""
828	return self._get_column(self._getparangle, row)
829
830	def get_direction(self, row=-1):
831	"""
832	Get a list of Positions on the sky (direction) for the observations.
833	Return a string for each integration in the scantable.
834
835	Parameters:
836
837	row: row no of integration. Default -1 return all rows
838
839	"""
840	return self._get_column(self._getdirection, row)
841
842	def get_directionval(self, row=-1):
843	"""\
844	Get a list of Positions on the sky (direction) for the observations.
845	Return a float for each integration in the scantable.
846
847	Parameters:
848
849	row: row no of integration. Default -1 return all rows
850
851	"""
852	return self._get_column(self._getdirectionvec, row)
853
854	@print_log_dec
855	def set_unit(self, unit='channel'):
856	"""\
857	Set the unit for all following operations on this scantable
858
859	Parameters:
860
861	unit: optional unit, default is 'channel'. Use one of '*Hz',
862	'km/s', 'channel' or equivalent ''
863
864	"""
865	varlist = vars()
866	if unit in ['', 'pixel', 'channel']:
867	unit = ''
868	inf = list(self._getcoordinfo())
869	inf[0] = unit
870	self._setcoordinfo(inf)
871	self._add_history("set_unit", varlist)
872
873	@print_log_dec
874	def set_instrument(self, instr):
875	"""\
876	Set the instrument for subsequent processing.
877
878	Parameters:
879
880	instr: Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
881	'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
882
883	"""
884	self._setInstrument(instr)
885	self._add_history("set_instument", vars())
886
887	@print_log_dec
888	def set_feedtype(self, feedtype):
889	"""\
890	Overwrite the feed type, which might not be set correctly.
891
892	Parameters:
893
894	feedtype: 'linear' or 'circular'
895
896	"""
897	self._setfeedtype(feedtype)
898	self._add_history("set_feedtype", vars())
899
900	@print_log_dec
901	def set_doppler(self, doppler='RADIO'):
902	"""\
903	Set the doppler for all following operations on this scantable.
904
905	Parameters:
906
907	doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
908
909	"""
910	varlist = vars()
911	inf = list(self._getcoordinfo())
912	inf[2] = doppler
913	self._setcoordinfo(inf)
914	self._add_history("set_doppler", vars())
915
916	@print_log_dec
917	def set_freqframe(self, frame=None):
918	"""\
919	Set the frame type of the Spectral Axis.
920
921	Parameters:
922
923	frame: an optional frame type, default 'LSRK'. Valid frames are:
924	'TOPO', 'LSRD', 'LSRK', 'BARY',
925	'GEO', 'GALACTO', 'LGROUP', 'CMB'
926
927	Example::
928
929	scan.set_freqframe('BARY')
930
931	"""
932	frame = frame or rcParams['scantable.freqframe']
933	varlist = vars()
934	# "REST" is not implemented in casacore
935	#valid = ['REST', 'TOPO', 'LSRD', 'LSRK', 'BARY', \
936	# 'GEO', 'GALACTO', 'LGROUP', 'CMB']
937	valid = ['TOPO', 'LSRD', 'LSRK', 'BARY', \
938	'GEO', 'GALACTO', 'LGROUP', 'CMB']
939
940	if frame in valid:
941	inf = list(self._getcoordinfo())
942	inf[1] = frame
943	self._setcoordinfo(inf)
944	self._add_history("set_freqframe", varlist)
945	else:
946	msg = "Please specify a valid freq type. Valid types are:\n", valid
947	if rcParams['verbose']:
948	#print msg
949	asaplog.push( msg )
950	print_log( 'ERROR' )
951	else:
952	raise TypeError(msg)
953
954	def set_dirframe(self, frame=""):
955	"""\
956	Set the frame type of the Direction on the sky.
957
958	Parameters:
959
960	frame: an optional frame type, default ''. Valid frames are:
961	'J2000', 'B1950', 'GALACTIC'
962
963	Example:
964
965	scan.set_dirframe('GALACTIC')
966
967	"""
968	varlist = vars()
969	try:
970	Scantable.set_dirframe(self, frame)
971	except RuntimeError, msg:
972	if rcParams['verbose']:
973	#print msg
974	print_log()
975	asaplog.push( str(msg) )
976	print_log( 'ERROR' )
977	else:
978	raise
979	self._add_history("set_dirframe", varlist)
980
981	def get_unit(self):
982	"""\
983	Get the default unit set in this scantable
984
985	Returns:
986
987	A unit string
988
989	"""
990	inf = self._getcoordinfo()
991	unit = inf[0]
992	if unit == '': unit = 'channel'
993	return unit
994
995	@print_log_dec
996	def get_abcissa(self, rowno=0):
997	"""\
998	Get the abcissa in the current coordinate setup for the currently
999	selected Beam/IF/Pol
1000
1001	Parameters:
1002
1003	rowno: an optional row number in the scantable. Default is the
1004	first row, i.e. rowno=0
1005
1006	Returns:
1007
1008	The abcissa values and the format string (as a dictionary)
1009
1010	"""
1011	abc = self._getabcissa(rowno)
1012	lbl = self._getabcissalabel(rowno)
1013	return abc, lbl
1014
1015	def flag(self, mask=None, unflag=False):
1016	"""\
1017	Flag the selected data using an optional channel mask.
1018
1019	Parameters:
1020
1021	mask: an optional channel mask, created with create_mask. Default
1022	(no mask) is all channels.
1023
1024	unflag: if True, unflag the data
1025
1026	"""
1027	varlist = vars()
1028	mask = mask or []
1029	try:
1030	self._flag(mask, unflag)
1031	except RuntimeError, msg:
1032	if rcParams['verbose']:
1033	#print msg
1034	print_log()
1035	asaplog.push( str(msg) )
1036	print_log( 'ERROR' )
1037	return
1038	else: raise
1039	self._add_history("flag", varlist)
1040
1041	def flag_row(self, rows=[], unflag=False):
1042	"""\
1043	Flag the selected data in row-based manner.
1044
1045	Parameters:
1046
1047	rows: list of row numbers to be flagged. Default is no row
1048	(must be explicitly specified to execute row-based flagging).
1049
1050	unflag: if True, unflag the data.
1051
1052	"""
1053	varlist = vars()
1054	try:
1055	self._flag_row(rows, unflag)
1056	except RuntimeError, msg:
1057	if rcParams['verbose']:
1058	print_log()
1059	asaplog.push( str(msg) )
1060	print_log('ERROR')
1061	return
1062	else: raise
1063	self._add_history("flag_row", varlist)
1064
1065	def clip(self, uthres=None, dthres=None, clipoutside=True, unflag=False):
1066	"""\
1067	Flag the selected data outside a specified range (in channel-base)
1068
1069	Parameters:
1070
1071	uthres: upper threshold.
1072
1073	dthres: lower threshold
1074
1075	clipoutside: True for flagging data outside the range [dthres:uthres].
1076	False for glagging data inside the range.
1077
1078	unflag: if True, unflag the data.
1079
1080	"""
1081	varlist = vars()
1082	try:
1083	self._clip(uthres, dthres, clipoutside, unflag)
1084	except RuntimeError, msg:
1085	if rcParams['verbose']:
1086	print_log()
1087	asaplog.push(str(msg))
1088	print_log('ERROR')
1089	return
1090	else: raise
1091	self._add_history("clip", varlist)
1092
1093	@print_log_dec
1094	def lag_flag(self, start, end, unit="MHz", insitu=None):
1095	"""\
1096	Flag the data in 'lag' space by providing a frequency to remove.
1097	Flagged data in the scantable gets interpolated over the region.
1098	No taper is applied.
1099
1100	Parameters:
1101
1102	start: the start frequency (really a period within the
1103	bandwidth) or period to remove
1104
1105	end: the end frequency or period to remove
1106
1107	unit: the frequency unit (default "MHz") or "" for
1108	explicit lag channels
1109
1110	Notes:
1111
1112	It is recommended to flag edges of the band or strong
1113	signals beforehand.
1114
1115	"""
1116	if insitu is None: insitu = rcParams['insitu']
1117	self._math._setinsitu(insitu)
1118	varlist = vars()
1119	base = { "GHz": 1000000000., "MHz": 1000000., "kHz": 1000., "Hz": 1.}
1120	if not (unit == "" or base.has_key(unit)):
1121	raise ValueError("%s is not a valid unit." % unit)
1122	try:
1123	if unit == "":
1124	s = scantable(self._math._lag_flag(self, start, end, "lags"))
1125	else:
1126	s = scantable(self._math._lag_flag(self, start*base[unit],
1127	end*base[unit], "frequency"))
1128	except RuntimeError, msg:
1129	if rcParams['verbose']:
1130	#print msg
1131	print_log()
1132	asaplog.push( str(msg) )
1133	print_log( 'ERROR' )
1134	return
1135	else: raise
1136	s._add_history("lag_flag", varlist)
1137	if insitu:
1138	self._assign(s)
1139	else:
1140	return s
1141
1142	@print_log_dec
1143	def create_mask(self, args, *kwargs):
1144	"""\
1145	Compute and return a mask based on [min, max] windows.
1146	The specified windows are to be INCLUDED, when the mask is
1147	applied.
1148
1149	Parameters:
1150
1151	[min, max], [min2, max2], ...
1152	Pairs of start/end points (inclusive)specifying the regions
1153	to be masked
1154
1155	invert: optional argument. If specified as True,
1156	return an inverted mask, i.e. the regions
1157	specified are EXCLUDED
1158
1159	row: create the mask using the specified row for
1160	unit conversions, default is row=0
1161	only necessary if frequency varies over rows.
1162
1163	Examples::
1164
1165	scan.set_unit('channel')
1166	# a)
1167	msk = scan.create_mask([400, 500], [800, 900])
1168	# masks everything outside 400 and 500
1169	# and 800 and 900 in the unit 'channel'
1170
1171	# b)
1172	msk = scan.create_mask([400, 500], [800, 900], invert=True)
1173	# masks the regions between 400 and 500
1174	# and 800 and 900 in the unit 'channel'
1175
1176	# c)
1177	#mask only channel 400
1178	msk = scan.create_mask([400])
1179
1180	"""
1181	row = kwargs.get("row", 0)
1182	data = self._getabcissa(row)
1183	u = self._getcoordinfo()[0]
1184	if rcParams['verbose']:
1185	if u == "": u = "channel"
1186	msg = "The current mask window unit is %s" % u
1187	i = self._check_ifs()
1188	if not i:
1189	msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
1190	asaplog.push(msg)
1191	n = self.nchan()
1192	msk = _n_bools(n, False)
1193	# test if args is a 'list' or a 'normal *args - UGLY!!!
1194
1195	ws = (isinstance(args[-1][-1], int) or isinstance(args[-1][-1], float)) \
1196	and args or args[0]
1197	for window in ws:
1198	if len(window) == 1:
1199	window = [window[0], window[0]]
1200	if len(window) == 0 or len(window) > 2:
1201	raise ValueError("A window needs to be defined as [start(, end)]")
1202	if window[0] > window[1]:
1203	tmp = window[0]
1204	window[0] = window[1]
1205	window[1] = tmp
1206	for i in range(n):
1207	if data[i] >= window[0] and data[i] <= window[1]:
1208	msk[i] = True
1209	if kwargs.has_key('invert'):
1210	if kwargs.get('invert'):
1211	msk = mask_not(msk)
1212	return msk
1213
1214	def get_masklist(self, mask=None, row=0):
1215	"""\
1216	Compute and return a list of mask windows, [min, max].
1217
1218	Parameters:
1219
1220	mask: channel mask, created with create_mask.
1221
1222	row: calcutate the masklist using the specified row
1223	for unit conversions, default is row=0
1224	only necessary if frequency varies over rows.
1225
1226	Returns:
1227
1228	[min, max], [min2, max2], ...
1229	Pairs of start/end points (inclusive)specifying
1230	the masked regions
1231
1232	"""
1233	if not (isinstance(mask,list) or isinstance(mask, tuple)):
1234	raise TypeError("The mask should be list or tuple.")
1235	if len(mask) < 2:
1236	raise TypeError("The mask elements should be > 1")
1237	if self.nchan() != len(mask):
1238	msg = "Number of channels in scantable != number of mask elements"
1239	raise TypeError(msg)
1240	data = self._getabcissa(row)
1241	u = self._getcoordinfo()[0]
1242	if rcParams['verbose']:
1243	if u == "": u = "channel"
1244	msg = "The current mask window unit is %s" % u
1245	i = self._check_ifs()
1246	if not i:
1247	msg += "\nThis mask is only valid for IF=%d" % (self.getif(i))
1248	asaplog.push(msg)
1249	masklist=[]
1250	ist, ien = None, None
1251	ist, ien=self.get_mask_indices(mask)
1252	if ist is not None and ien is not None:
1253	for i in xrange(len(ist)):
1254	range=[data[ist[i]],data[ien[i]]]
1255	range.sort()
1256	masklist.append([range[0],range[1]])
1257	return masklist
1258
1259	def get_mask_indices(self, mask=None):
1260	"""\
1261	Compute and Return lists of mask start indices and mask end indices.
1262
1263	Parameters:
1264
1265	mask: channel mask, created with create_mask.
1266
1267	Returns:
1268
1269	List of mask start indices and that of mask end indices,
1270	i.e., [istart1,istart2,....], [iend1,iend2,....].
1271
1272	"""
1273	if not (isinstance(mask,list) or isinstance(mask, tuple)):
1274	raise TypeError("The mask should be list or tuple.")
1275	if len(mask) < 2:
1276	raise TypeError("The mask elements should be > 1")
1277	istart=[]
1278	iend=[]
1279	if mask[0]: istart.append(0)
1280	for i in range(len(mask)-1):
1281	if not mask[i] and mask[i+1]:
1282	istart.append(i+1)
1283	elif mask[i] and not mask[i+1]:
1284	iend.append(i)
1285	if mask[len(mask)-1]: iend.append(len(mask)-1)
1286	if len(istart) != len(iend):
1287	raise RuntimeError("Numbers of mask start != mask end.")
1288	for i in range(len(istart)):
1289	if istart[i] > iend[i]:
1290	raise RuntimeError("Mask start index > mask end index")
1291	break
1292	return istart,iend
1293
1294	# def get_restfreqs(self):
1295	# """
1296	# Get the restfrequency(s) stored in this scantable.
1297	# The return value(s) are always of unit 'Hz'
1298	# Parameters:
1299	# none
1300	# Returns:
1301	# a list of doubles
1302	# """
1303	# return list(self._getrestfreqs())
1304
1305	def get_restfreqs(self, ids=None):
1306	"""\
1307	Get the restfrequency(s) stored in this scantable.
1308	The return value(s) are always of unit 'Hz'
1309
1310	Parameters:
1311
1312	ids: (optional) a list of MOLECULE_ID for that restfrequency(s) to
1313	be retrieved
1314
1315	Returns:
1316
1317	dictionary containing ids and a list of doubles for each id
1318
1319	"""
1320	if ids is None:
1321	rfreqs={}
1322	idlist = self.getmolnos()
1323	for i in idlist:
1324	rfreqs[i]=list(self._getrestfreqs(i))
1325	return rfreqs
1326	else:
1327	if type(ids)==list or type(ids)==tuple:
1328	rfreqs={}
1329	for i in ids:
1330	rfreqs[i]=list(self._getrestfreqs(i))
1331	return rfreqs
1332	else:
1333	return list(self._getrestfreqs(ids))
1334	#return list(self._getrestfreqs(ids))
1335
1336	def set_restfreqs(self, freqs=None, unit='Hz'):
1337	"""\
1338	Set or replace the restfrequency specified and
1339	If the 'freqs' argument holds a scalar,
1340	then that rest frequency will be applied to all the selected
1341	data. If the 'freqs' argument holds
1342	a vector, then it MUST be of equal or smaller length than
1343	the number of IFs (and the available restfrequencies will be
1344	replaced by this vector). In this case, all data have
1345	the restfrequency set per IF according
1346	to the corresponding value you give in the 'freqs' vector.
1347	E.g. 'freqs=[1e9, 2e9]' would mean IF 0 gets restfreq 1e9 and
1348	IF 1 gets restfreq 2e9.
1349
1350	You can also specify the frequencies via a linecatalog.
1351
1352	Parameters:
1353
1354	freqs: list of rest frequency values or string idenitfiers
1355
1356	unit: unit for rest frequency (default 'Hz')
1357
1358
1359	Example::
1360
1361	# set the given restfrequency for the all currently selected IFs
1362	scan.set_restfreqs(freqs=1.4e9)
1363	# set restfrequencies for the n IFs (n > 1) in the order of the
1364	# list, i.e
1365	# IF0 -> 1.4e9, IF1 -> 1.41e9, IF3 -> 1.42e9
1366	# len(list_of_restfreqs) == nIF
1367	# for nIF == 1 the following will set multiple restfrequency for
1368	# that IF
1369	scan.set_restfreqs(freqs=[1.4e9, 1.41e9, 1.42e9])
1370	# set multiple restfrequencies per IF. as a list of lists where
1371	# the outer list has nIF elements, the inner s arbitrary
1372	scan.set_restfreqs(freqs=[[1.4e9, 1.41e9], [1.67e9]])
1373
1374	Note:
1375
1376	To do more sophisticate Restfrequency setting, e.g. on a
1377	source and IF basis, use scantable.set_selection() before using
1378	this function::
1379
1380	# provided your scantable is called scan
1381	selection = selector()
1382	selection.set_name("ORION*")
1383	selection.set_ifs([1])
1384	scan.set_selection(selection)
1385	scan.set_restfreqs(freqs=86.6e9)
1386
1387	"""
1388	varlist = vars()
1389	from asap import linecatalog
1390	# simple value
1391	if isinstance(freqs, int) or isinstance(freqs, float):
1392	self._setrestfreqs([freqs], [""], unit)
1393	# list of values
1394	elif isinstance(freqs, list) or isinstance(freqs, tuple):
1395	# list values are scalars
1396	if isinstance(freqs[-1], int) or isinstance(freqs[-1], float):
1397	if len(freqs) == 1:
1398	self._setrestfreqs(freqs, [""], unit)
1399	else:
1400	# allow the 'old' mode of setting mulitple IFs
1401	sel = selector()
1402	savesel = self._getselection()
1403	iflist = self.getifnos()
1404	if len(freqs)>len(iflist):
1405	raise ValueError("number of elements in list of list "
1406	"exeeds the current IF selections")
1407	iflist = self.getifnos()
1408	for i, fval in enumerate(freqs):
1409	sel.set_ifs(iflist[i])
1410	self._setselection(sel)
1411	self._setrestfreqs([fval], [""], unit)
1412	self._setselection(savesel)
1413
1414	# list values are dict, {'value'=, 'name'=)
1415	elif isinstance(freqs[-1], dict):
1416	values = []
1417	names = []
1418	for d in freqs:
1419	values.append(d["value"])
1420	names.append(d["name"])
1421	self._setrestfreqs(values, names, unit)
1422	elif isinstance(freqs[-1], list) or isinstance(freqs[-1], tuple):
1423	sel = selector()
1424	savesel = self._getselection()
1425	iflist = self.getifnos()
1426	if len(freqs)>len(iflist):
1427	raise ValueError("number of elements in list of list exeeds"
1428	" the current IF selections")
1429	for i, fval in enumerate(freqs):
1430	sel.set_ifs(iflist[i])
1431	self._setselection(sel)
1432	self._setrestfreqs(fval, [""], unit)
1433	self._setselection(savesel)
1434	# freqs are to be taken from a linecatalog
1435	elif isinstance(freqs, linecatalog):
1436	sel = selector()
1437	savesel = self._getselection()
1438	for i in xrange(freqs.nrow()):
1439	sel.set_ifs(iflist[i])
1440	self._setselection(sel)
1441	self._setrestfreqs([freqs.get_frequency(i)],
1442	[freqs.get_name(i)], "MHz")
1443	# ensure that we are not iterating past nIF
1444	if i == self.nif()-1: break
1445	self._setselection(savesel)
1446	else:
1447	return
1448	self._add_history("set_restfreqs", varlist)
1449
1450	def shift_refpix(self, delta):
1451	"""\
1452	Shift the reference pixel of the Spectra Coordinate by an
1453	integer amount.
1454
1455	Parameters:
1456
1457	delta: the amount to shift by
1458
1459	Note:
1460
1461	Be careful using this with broadband data.
1462
1463	"""
1464	Scantable.shift_refpix(self, delta)
1465
1466	def history(self, filename=None):
1467	"""\
1468	Print the history. Optionally to a file.
1469
1470	Parameters:
1471
1472	filename: The name of the file to save the history to.
1473
1474	"""
1475	hist = list(self._gethistory())
1476	out = "-"*80
1477	for h in hist:
1478	if h.startswith("---"):
1479	out = "\n".join([out, h])
1480	else:
1481	items = h.split("##")
1482	date = items[0]
1483	func = items[1]
1484	items = items[2:]
1485	out += "\n"+date+"\n"
1486	out += "Function: %s\n Parameters:" % (func)
1487	for i in items:
1488	s = i.split("=")
1489	out += "\n %s = %s" % (s[0], s[1])
1490	out = "\n".join([out, "-"*80])
1491	if filename is not None:
1492	if filename is "":
1493	filename = 'scantable_history.txt'
1494	import os
1495	filename = os.path.expandvars(os.path.expanduser(filename))
1496	if not os.path.isdir(filename):
1497	data = open(filename, 'w')
1498	data.write(out)
1499	data.close()
1500	else:
1501	msg = "Illegal file name '%s'." % (filename)
1502	if rcParams['verbose']:
1503	#print msg
1504	asaplog.push( msg )
1505	print_log( 'ERROR' )
1506	else:
1507	raise IOError(msg)
1508	if rcParams['verbose']:
1509	try:
1510	from IPython.genutils import page as pager
1511	except ImportError:
1512	from pydoc import pager
1513	pager(out)
1514	else:
1515	return out
1516	return
1517	#
1518	# Maths business
1519	#
1520	@print_log_dec
1521	def average_time(self, mask=None, scanav=False, weight='tint', align=False):
1522	"""\
1523	Return the (time) weighted average of a scan.
1524
1525	Note:
1526
1527	in channels only - align if necessary
1528
1529	Parameters:
1530
1531	mask: an optional mask (only used for 'var' and 'tsys'
1532	weighting)
1533
1534	scanav: True averages each scan separately
1535	False (default) averages all scans together,
1536
1537	weight: Weighting scheme.
1538	'none' (mean no weight)
1539	'var' (1/var(spec) weighted)
1540	'tsys' (1/Tsys**2 weighted)
1541	'tint' (integration time weighted)
1542	'tintsys' (Tint/Tsys**2)
1543	'median' ( median averaging)
1544	The default is 'tint'
1545
1546	align: align the spectra in velocity before averaging. It takes
1547	the time of the first spectrum as reference time.
1548
1549	Example::
1550
1551	# time average the scantable without using a mask
1552	newscan = scan.average_time()
1553
1554	"""
1555	varlist = vars()
1556	weight = weight or 'TINT'
1557	mask = mask or ()
1558	scanav = (scanav and 'SCAN') or 'NONE'
1559	scan = (self, )
1560	try:
1561	if align:
1562	scan = (self.freq_align(insitu=False), )
1563	s = None
1564	if weight.upper() == 'MEDIAN':
1565	s = scantable(self._math._averagechannel(scan[0], 'MEDIAN',
1566	scanav))
1567	else:
1568	s = scantable(self._math._average(scan, mask, weight.upper(),
1569	scanav))
1570	except RuntimeError, msg:
1571	if rcParams['verbose']:
1572	#print msg
1573	print_log()
1574	asaplog.push( str(msg) )
1575	print_log( 'ERROR' )
1576	return
1577	else: raise
1578	s._add_history("average_time", varlist)
1579	return s
1580
1581	@print_log_dec
1582	def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None):
1583	"""\
1584	Return a scan where all spectra are converted to either
1585	Jansky or Kelvin depending upon the flux units of the scan table.
1586	By default the function tries to look the values up internally.
1587	If it can't find them (or if you want to over-ride), you must
1588	specify EITHER jyperk OR eta (and D which it will try to look up
1589	also if you don't set it). jyperk takes precedence if you set both.
1590
1591	Parameters:
1592
1593	jyperk: the Jy / K conversion factor
1594
1595	eta: the aperture efficiency
1596
1597	d: the geomtric diameter (metres)
1598
1599	insitu: if False a new scantable is returned.
1600	Otherwise, the scaling is done in-situ
1601	The default is taken from .asaprc (False)
1602
1603	"""
1604	if insitu is None: insitu = rcParams['insitu']
1605	self._math._setinsitu(insitu)
1606	varlist = vars()
1607	jyperk = jyperk or -1.0
1608	d = d or -1.0
1609	eta = eta or -1.0
1610	s = scantable(self._math._convertflux(self, d, eta, jyperk))
1611	s._add_history("convert_flux", varlist)
1612	if insitu: self._assign(s)
1613	else: return s
1614
1615	@print_log_dec
1616	def gain_el(self, poly=None, filename="", method="linear", insitu=None):
1617	"""\
1618	Return a scan after applying a gain-elevation correction.
1619	The correction can be made via either a polynomial or a
1620	table-based interpolation (and extrapolation if necessary).
1621	You specify polynomial coefficients, an ascii table or neither.
1622	If you specify neither, then a polynomial correction will be made
1623	with built in coefficients known for certain telescopes (an error
1624	will occur if the instrument is not known).
1625	The data and Tsys are divided by the scaling factors.
1626
1627	Parameters:
1628
1629	poly: Polynomial coefficients (default None) to compute a
1630	gain-elevation correction as a function of
1631	elevation (in degrees).
1632
1633	filename: The name of an ascii file holding correction factors.
1634	The first row of the ascii file must give the column
1635	names and these MUST include columns
1636	"ELEVATION" (degrees) and "FACTOR" (multiply data
1637	by this) somewhere.
1638	The second row must give the data type of the
1639	column. Use 'R' for Real and 'I' for Integer.
1640	An example file would be
1641	(actual factors are arbitrary) :
1642
1643	TIME ELEVATION FACTOR
1644	R R R
1645	0.1 0 0.8
1646	0.2 20 0.85
1647	0.3 40 0.9
1648	0.4 60 0.85
1649	0.5 80 0.8
1650	0.6 90 0.75
1651
1652	method: Interpolation method when correcting from a table.
1653	Values are "nearest", "linear" (default), "cubic"
1654	and "spline"
1655
1656	insitu: if False a new scantable is returned.
1657	Otherwise, the scaling is done in-situ
1658	The default is taken from .asaprc (False)
1659
1660	"""
1661
1662	if insitu is None: insitu = rcParams['insitu']
1663	self._math._setinsitu(insitu)
1664	varlist = vars()
1665	poly = poly or ()
1666	from os.path import expandvars
1667	filename = expandvars(filename)
1668	s = scantable(self._math._gainel(self, poly, filename, method))
1669	s._add_history("gain_el", varlist)
1670	if insitu:
1671	self._assign(s)
1672	else:
1673	return s
1674
1675	@print_log_dec
1676	def freq_align(self, reftime=None, method='cubic', insitu=None):
1677	"""\
1678	Return a scan where all rows have been aligned in frequency/velocity.
1679	The alignment frequency frame (e.g. LSRK) is that set by function
1680	set_freqframe.
1681
1682	Parameters:
1683
1684	reftime: reference time to align at. By default, the time of
1685	the first row of data is used.
1686
1687	method: Interpolation method for regridding the spectra.
1688	Choose from "nearest", "linear", "cubic" (default)
1689	and "spline"
1690
1691	insitu: if False a new scantable is returned.
1692	Otherwise, the scaling is done in-situ
1693	The default is taken from .asaprc (False)
1694
1695	"""
1696	if insitu is None: insitu = rcParams["insitu"]
1697	self._math._setinsitu(insitu)
1698	varlist = vars()
1699	reftime = reftime or ""
1700	s = scantable(self._math._freq_align(self, reftime, method))
1701	s._add_history("freq_align", varlist)
1702	if insitu: self._assign(s)
1703	else: return s
1704
1705	@print_log_dec
1706	def opacity(self, tau=None, insitu=None):
1707	"""\
1708	Apply an opacity correction. The data
1709	and Tsys are multiplied by the correction factor.
1710
1711	Parameters:
1712
1713	tau: (list of) opacity from which the correction factor is
1714	exp(tau*ZD)
1715	where ZD is the zenith-distance.
1716	If a list is provided, it has to be of length nIF,
1717	nIF*nPol or 1 and in order of IF/POL, e.g.
1718	[opif0pol0, opif0pol1, opif1pol0 ...]
1719	if tau is `None` the opacities are determined from a
1720	model.
1721
1722	insitu: if False a new scantable is returned.
1723	Otherwise, the scaling is done in-situ
1724	The default is taken from .asaprc (False)
1725
1726	"""
1727	if insitu is None: insitu = rcParams['insitu']
1728	self._math._setinsitu(insitu)
1729	varlist = vars()
1730	if not hasattr(tau, "__len__"):
1731	tau = [tau]
1732	s = scantable(self._math._opacity(self, tau))
1733	s._add_history("opacity", varlist)
1734	if insitu: self._assign(s)
1735	else: return s
1736
1737	@print_log_dec
1738	def bin(self, width=5, insitu=None):
1739	"""\
1740	Return a scan where all spectra have been binned up.
1741
1742	Parameters:
1743
1744	width: The bin width (default=5) in pixels
1745
1746	insitu: if False a new scantable is returned.
1747	Otherwise, the scaling is done in-situ
1748	The default is taken from .asaprc (False)
1749
1750	"""
1751	if insitu is None: insitu = rcParams['insitu']
1752	self._math._setinsitu(insitu)
1753	varlist = vars()
1754	s = scantable(self._math._bin(self, width))
1755	s._add_history("bin", varlist)
1756	if insitu:
1757	self._assign(s)
1758	else:
1759	return s
1760
1761	@print_log_dec
1762	def resample(self, width=5, method='cubic', insitu=None):
1763	"""\
1764	Return a scan where all spectra have been binned up.
1765
1766	Parameters:
1767
1768	width: The bin width (default=5) in pixels
1769
1770	method: Interpolation method when correcting from a table.
1771	Values are "nearest", "linear", "cubic" (default)
1772	and "spline"
1773
1774	insitu: if False a new scantable is returned.
1775	Otherwise, the scaling is done in-situ
1776	The default is taken from .asaprc (False)
1777
1778	"""
1779	if insitu is None: insitu = rcParams['insitu']
1780	self._math._setinsitu(insitu)
1781	varlist = vars()
1782	s = scantable(self._math._resample(self, method, width))
1783	s._add_history("resample", varlist)
1784	if insitu: self._assign(s)
1785	else: return s
1786
1787	@print_log_dec
1788	def average_pol(self, mask=None, weight='none'):
1789	"""\
1790	Average the Polarisations together.
1791
1792	Parameters:
1793
1794	mask: An optional mask defining the region, where the
1795	averaging will be applied. The output will have all
1796	specified points masked.
1797
1798	weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
1799	weighted), or 'tsys' (1/Tsys**2 weighted)
1800
1801	"""
1802	varlist = vars()
1803	mask = mask or ()
1804	s = scantable(self._math._averagepol(self, mask, weight.upper()))
1805	s._add_history("average_pol", varlist)
1806	return s
1807
1808	@print_log_dec
1809	def average_beam(self, mask=None, weight='none'):
1810	"""\
1811	Average the Beams together.
1812
1813	Parameters:
1814	mask: An optional mask defining the region, where the
1815	averaging will be applied. The output will have all
1816	specified points masked.
1817
1818	weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
1819	weighted), or 'tsys' (1/Tsys**2 weighted)
1820
1821	"""
1822	varlist = vars()
1823	mask = mask or ()
1824	s = scantable(self._math._averagebeams(self, mask, weight.upper()))
1825	s._add_history("average_beam", varlist)
1826	return s
1827
1828	def parallactify(self, pflag):
1829	"""\
1830	Set a flag to indicate whether this data should be treated as having
1831	been 'parallactified' (total phase == 0.0)
1832
1833	Parameters:
1834
1835	pflag: Bool indicating whether to turn this on (True) or
1836	off (False)
1837
1838	"""
1839	varlist = vars()
1840	self._parallactify(pflag)
1841	self._add_history("parallactify", varlist)
1842
1843	@print_log_dec
1844	def convert_pol(self, poltype=None):
1845	"""\
1846	Convert the data to a different polarisation type.
1847	Note that you will need cross-polarisation terms for most conversions.
1848
1849	Parameters:
1850
1851	poltype: The new polarisation type. Valid types are:
1852	"linear", "circular", "stokes" and "linpol"
1853
1854	"""
1855	varlist = vars()
1856	try:
1857	s = scantable(self._math._convertpol(self, poltype))
1858	except RuntimeError, msg:
1859	if rcParams['verbose']:
1860	#print msg
1861	print_log()
1862	asaplog.push( str(msg) )
1863	print_log( 'ERROR' )
1864	return
1865	else:
1866	raise
1867	s._add_history("convert_pol", varlist)
1868	return s
1869
1870	@print_log_dec
1871	def smooth(self, kernel="hanning", width=5.0, order=2, plot=False, insitu=None):
1872	"""\
1873	Smooth the spectrum by the specified kernel (conserving flux).
1874
1875	Parameters:
1876
1877	kernel: The type of smoothing kernel. Select from
1878	'hanning' (default), 'gaussian', 'boxcar', 'rmedian'
1879	or 'poly'
1880
1881	width: The width of the kernel in pixels. For hanning this is
1882	ignored otherwise it defauls to 5 pixels.
1883	For 'gaussian' it is the Full Width Half
1884	Maximum. For 'boxcar' it is the full width.
1885	For 'rmedian' and 'poly' it is the half width.
1886
1887	order: Optional parameter for 'poly' kernel (default is 2), to
1888	specify the order of the polnomial. Ignored by all other
1889	kernels.
1890
1891	plot: plot the original and the smoothed spectra.
1892	In this each indivual fit has to be approved, by
1893	typing 'y' or 'n'
1894
1895	insitu: if False a new scantable is returned.
1896	Otherwise, the scaling is done in-situ
1897	The default is taken from .asaprc (False)
1898
1899	"""
1900	if insitu is None: insitu = rcParams['insitu']
1901	self._math._setinsitu(insitu)
1902	varlist = vars()
1903
1904	if plot: orgscan = self.copy()
1905
1906	s = scantable(self._math._smooth(self, kernel.lower(), width, order))
1907	s._add_history("smooth", varlist)
1908
1909	if plot:
1910	if rcParams['plotter.gui']:
1911	from asap.asaplotgui import asaplotgui as asaplot
1912	else:
1913	from asap.asaplot import asaplot
1914	self._p=asaplot()
1915	self._p.set_panels()
1916	ylab=s._get_ordinate_label()
1917	#self._p.palette(0,["#777777","red"])
1918	for r in xrange(s.nrow()):
1919	xsm=s._getabcissa(r)
1920	ysm=s._getspectrum(r)
1921	xorg=orgscan._getabcissa(r)
1922	yorg=orgscan._getspectrum(r)
1923	self._p.clear()
1924	self._p.hold()
1925	self._p.set_axes('ylabel',ylab)
1926	self._p.set_axes('xlabel',s._getabcissalabel(r))
1927	self._p.set_axes('title',s._getsourcename(r))
1928	self._p.set_line(label='Original',color="#777777")
1929	self._p.plot(xorg,yorg)
1930	self._p.set_line(label='Smoothed',color="red")
1931	self._p.plot(xsm,ysm)
1932	### Ugly part for legend
1933	for i in [0,1]:
1934	self._p.subplots[0]['lines'].append([self._p.subplots[0]['axes'].lines[i]])
1935	self._p.release()
1936	### Ugly part for legend
1937	self._p.subplots[0]['lines']=[]
1938	res = raw_input("Accept smoothing ([y]/n): ")
1939	if res.upper() == 'N':
1940	s._setspectrum(yorg, r)
1941	self._p.unmap()
1942	self._p = None
1943	del orgscan
1944
1945	if insitu: self._assign(s)
1946	else: return s
1947
1948	@print_log_dec
1949	def poly_baseline(self, mask=None, order=0, plot=False, uselin=False,
1950	insitu=None):
1951	"""\
1952	Return a scan which has been baselined (all rows) by a polynomial.
1953
1954	Parameters:
1955
1956	mask: an optional mask
1957
1958	order: the order of the polynomial (default is 0)
1959
1960	plot: plot the fit and the residual. In this each
1961	indivual fit has to be approved, by typing 'y'
1962	or 'n'
1963
1964	uselin: use linear polynomial fit
1965
1966	insitu: if False a new scantable is returned.
1967	Otherwise, the scaling is done in-situ
1968	The default is taken from .asaprc (False)
1969
1970	Example::
1971
1972	# return a scan baselined by a third order polynomial,
1973	# not using a mask
1974	bscan = scan.poly_baseline(order=3)
1975
1976	"""
1977	if insitu is None: insitu = rcParams['insitu']
1978	if not insitu:
1979	workscan = self.copy()
1980	else:
1981	workscan = self
1982	varlist = vars()
1983	if mask is None:
1984	mask = [True for i in xrange(self.nchan(-1))]
1985
1986	from asap.asapfitter import fitter
1987	try:
1988	f = fitter()
1989	if uselin:
1990	f.set_function(lpoly=order)
1991	else:
1992	f.set_function(poly=order)
1993
1994	rows = range(workscan.nrow())
1995	if len(rows) > 0:
1996	self.blpars = []
1997
1998	for r in rows:
1999	# take into account flagtra info (CAS-1434)
2000	flagtra = workscan._getmask(r)
2001	actualmask = mask[:]
2002	if len(actualmask) == 0:
2003	actualmask = list(flagtra[:])
2004	else:
2005	if len(actualmask) != len(flagtra):
2006	raise RuntimeError, "Mask and flagtra have different length"
2007	else:
2008	for i in range(0, len(actualmask)):
2009	actualmask[i] = actualmask[i] and flagtra[i]
2010	f.set_scan(workscan, actualmask)
2011	f.x = workscan._getabcissa(r)
2012	f.y = workscan._getspectrum(r)
2013	f.data = None
2014	f.fit()
2015	if plot:
2016	f.plot(residual=True)
2017	x = raw_input("Accept fit ( [y]/n ): ")
2018	if x.upper() == 'N':
2019	self.blpars.append(None)
2020	continue
2021	workscan._setspectrum(f.fitter.getresidual(), r)
2022	self.blpars.append(f.get_parameters())
2023
2024	if plot:
2025	f._p.unmap()
2026	f._p = None
2027	workscan._add_history("poly_baseline", varlist)
2028	if insitu:
2029	self._assign(workscan)
2030	else:
2031	return workscan
2032	except RuntimeError:
2033	msg = "The fit failed, possibly because it didn't converge."
2034	if rcParams['verbose']:
2035	#print msg
2036	print_log()
2037	asaplog.push( str(msg) )
2038	print_log( 'ERROR' )
2039	return
2040	else:
2041	raise RuntimeError(msg)
2042
2043
2044	def auto_poly_baseline(self, mask=[], edge=(0, 0), order=0,
2045	threshold=3, chan_avg_limit=1, plot=False,
2046	insitu=None):
2047	"""\
2048	Return a scan which has been baselined (all rows) by a polynomial.
2049	Spectral lines are detected first using linefinder and masked out
2050	to avoid them affecting the baseline solution.
2051
2052	Parameters:
2053
2054	mask: an optional mask retreived from scantable
2055
2056	edge: an optional number of channel to drop at the edge of
2057	spectrum. If only one value is
2058	specified, the same number will be dropped from
2059	both sides of the spectrum. Default is to keep
2060	all channels. Nested tuples represent individual
2061	edge selection for different IFs (a number of spectral
2062	channels can be different)
2063
2064	order: the order of the polynomial (default is 0)
2065
2066	threshold: the threshold used by line finder. It is better to
2067	keep it large as only strong lines affect the
2068	baseline solution.
2069
2070	chan_avg_limit:
2071	a maximum number of consequtive spectral channels to
2072	average during the search of weak and broad lines.
2073	The default is no averaging (and no search for weak
2074	lines). If such lines can affect the fitted baseline
2075	(e.g. a high order polynomial is fitted), increase this
2076	parameter (usually values up to 8 are reasonable). Most
2077	users of this method should find the default value
2078	sufficient.
2079
2080	plot: plot the fit and the residual. In this each
2081	indivual fit has to be approved, by typing 'y'
2082	or 'n'
2083
2084	insitu: if False a new scantable is returned.
2085	Otherwise, the scaling is done in-situ
2086	The default is taken from .asaprc (False)
2087
2088
2089	Example::
2090
2091	scan2 = scan.auto_poly_baseline(order=7, insitu=False)
2092
2093	"""
2094	if insitu is None: insitu = rcParams['insitu']
2095	varlist = vars()
2096	from asap.asapfitter import fitter
2097	from asap.asaplinefind import linefinder
2098	from asap import _is_sequence_or_number as _is_valid
2099
2100	# check whether edge is set up for each IF individually
2101	individualedge = False;
2102	if len(edge) > 1:
2103	if isinstance(edge[0], list) or isinstance(edge[0], tuple):
2104	individualedge = True;
2105
2106	if not _is_valid(edge, int) and not individualedge:
2107	raise ValueError, "Parameter 'edge' has to be an integer or a \
2108	pair of integers specified as a tuple. Nested tuples are allowed \
2109	to make individual selection for different IFs."
2110
2111	curedge = (0, 0)
2112	if individualedge:
2113	for edgepar in edge:
2114	if not _is_valid(edgepar, int):
2115	raise ValueError, "Each element of the 'edge' tuple has \
2116	to be a pair of integers or an integer."
2117	else:
2118	curedge = edge;
2119
2120	# setup fitter
2121	f = fitter()
2122	f.set_function(poly=order)
2123
2124	# setup line finder
2125	fl = linefinder()
2126	fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
2127
2128	if not insitu:
2129	workscan = self.copy()
2130	else:
2131	workscan = self
2132
2133	fl.set_scan(workscan)
2134
2135	rows = range(workscan.nrow())
2136	# Save parameters of baseline fits & masklists as a class attribute.
2137	# NOTICE: It does not reflect changes in scantable!
2138	if len(rows) > 0:
2139	self.blpars=[]
2140	self.masklists=[]
2141	asaplog.push("Processing:")
2142	for r in rows:
2143	msg = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % \
2144	(workscan.getscan(r), workscan.getbeam(r), workscan.getif(r), \
2145	workscan.getpol(r), workscan.getcycle(r))
2146	asaplog.push(msg, False)
2147
2148	# figure out edge parameter
2149	if individualedge:
2150	if len(edge) >= workscan.getif(r):
2151	raise RuntimeError, "Number of edge elements appear to " \
2152	"be less than the number of IFs"
2153	curedge = edge[workscan.getif(r)]
2154
2155	# take into account flagtra info (CAS-1434)
2156	flagtra = workscan._getmask(r)
2157	actualmask = mask[:]
2158	if len(actualmask) == 0:
2159	actualmask = list(flagtra[:])
2160	else:
2161	if len(actualmask) != len(flagtra):
2162	raise RuntimeError, "Mask and flagtra have different length"
2163	else:
2164	for i in range(0, len(actualmask)):
2165	actualmask[i] = actualmask[i] and flagtra[i]
2166
2167	# setup line finder
2168	fl.find_lines(r, actualmask, curedge)
2169	outmask=fl.get_mask()
2170	f.set_scan(workscan, fl.get_mask())
2171	f.x = workscan._getabcissa(r)
2172	f.y = workscan._getspectrum(r)
2173	f.data = None
2174	f.fit()
2175
2176	# Show mask list
2177	masklist=workscan.get_masklist(fl.get_mask(),row=r)
2178	msg = "mask range: "+str(masklist)
2179	asaplog.push(msg, False)
2180
2181	if plot:
2182	f.plot(residual=True)
2183	x = raw_input("Accept fit ( [y]/n ): ")
2184	if x.upper() == 'N':
2185	self.blpars.append(None)
2186	self.masklists.append(None)
2187	continue
2188
2189	workscan._setspectrum(f.fitter.getresidual(), r)
2190	self.blpars.append(f.get_parameters())
2191	self.masklists.append(masklist)
2192	if plot:
2193	f._p.unmap()
2194	f._p = None
2195	workscan._add_history("auto_poly_baseline", varlist)
2196	if insitu:
2197	self._assign(workscan)
2198	else:
2199	return workscan
2200
2201	@print_log_dec
2202	def rotate_linpolphase(self, angle):
2203	"""\
2204	Rotate the phase of the complex polarization O=Q+iU correlation.
2205	This is always done in situ in the raw data. So if you call this
2206	function more than once then each call rotates the phase further.
2207
2208	Parameters:
2209
2210	angle: The angle (degrees) to rotate (add) by.
2211
2212	Example::
2213
2214	scan.rotate_linpolphase(2.3)
2215
2216	"""
2217	varlist = vars()
2218	self._math._rotate_linpolphase(self, angle)
2219	self._add_history("rotate_linpolphase", varlist)
2220	return
2221
2222	@print_log_dec
2223	def rotate_xyphase(self, angle):
2224	"""\
2225	Rotate the phase of the XY correlation. This is always done in situ
2226	in the data. So if you call this function more than once
2227	then each call rotates the phase further.
2228
2229	Parameters:
2230
2231	angle: The angle (degrees) to rotate (add) by.
2232
2233	Example::
2234
2235	scan.rotate_xyphase(2.3)
2236
2237	"""
2238	varlist = vars()
2239	self._math._rotate_xyphase(self, angle)
2240	self._add_history("rotate_xyphase", varlist)
2241	return
2242
2243	@print_log_dec
2244	def swap_linears(self):
2245	"""\
2246	Swap the linear polarisations XX and YY, or better the first two
2247	polarisations as this also works for ciculars.
2248	"""
2249	varlist = vars()
2250	self._math._swap_linears(self)
2251	self._add_history("swap_linears", varlist)
2252	return
2253
2254	@print_log_dec
2255	def invert_phase(self):
2256	"""\
2257	Invert the phase of the complex polarisation
2258	"""
2259	varlist = vars()
2260	self._math._invert_phase(self)
2261	self._add_history("invert_phase", varlist)
2262	return
2263
2264	@print_log_dec
2265	def add(self, offset, insitu=None):
2266	"""\
2267	Return a scan where all spectra have the offset added
2268
2269	Parameters:
2270
2271	offset: the offset
2272
2273	insitu: if False a new scantable is returned.
2274	Otherwise, the scaling is done in-situ
2275	The default is taken from .asaprc (False)
2276
2277	"""
2278	if insitu is None: insitu = rcParams['insitu']
2279	self._math._setinsitu(insitu)
2280	varlist = vars()
2281	s = scantable(self._math._unaryop(self, offset, "ADD", False))
2282	s._add_history("add", varlist)
2283	if insitu:
2284	self._assign(s)
2285	else:
2286	return s
2287
2288	@print_log_dec
2289	def scale(self, factor, tsys=True, insitu=None):
2290	"""\
2291
2292	Return a scan where all spectra are scaled by the give 'factor'
2293
2294	Parameters:
2295
2296	factor: the scaling factor (float or 1D float list)
2297
2298	insitu: if False a new scantable is returned.
2299	Otherwise, the scaling is done in-situ
2300	The default is taken from .asaprc (False)
2301
2302	tsys: if True (default) then apply the operation to Tsys
2303	as well as the data
2304
2305	"""
2306	if insitu is None: insitu = rcParams['insitu']
2307	self._math._setinsitu(insitu)
2308	varlist = vars()
2309	s = None
2310	import numpy
2311	if isinstance(factor, list) or isinstance(factor, numpy.ndarray):
2312	if isinstance(factor[0], list) or isinstance(factor[0], numpy.ndarray):
2313	from asapmath import _array2dOp
2314	s = _array2dOp( self.copy(), factor, "MUL", tsys )
2315	else:
2316	s = scantable( self._math._arrayop( self.copy(), factor, "MUL", tsys ) )
2317	else:
2318	s = scantable(self._math._unaryop(self.copy(), factor, "MUL", tsys))
2319	s._add_history("scale", varlist)
2320	if insitu:
2321	self._assign(s)
2322	else:
2323	return s
2324
2325	def set_sourcetype(self, match, matchtype="pattern",
2326	sourcetype="reference"):
2327	"""\
2328	Set the type of the source to be an source or reference scan
2329	using the provided pattern.
2330
2331	Parameters:
2332
2333	match: a Unix style pattern, regular expression or selector
2334
2335	matchtype: 'pattern' (default) UNIX style pattern or
2336	'regex' regular expression
2337
2338	sourcetype: the type of the source to use (source/reference)
2339
2340	"""
2341	varlist = vars()
2342	basesel = self.get_selection()
2343	stype = -1
2344	if sourcetype.lower().startswith("r"):
2345	stype = 1
2346	elif sourcetype.lower().startswith("s"):
2347	stype = 0
2348	else:
2349	raise ValueError("Illegal sourcetype use s(ource) or r(eference)")
2350	if matchtype.lower().startswith("p"):
2351	matchtype = "pattern"
2352	elif matchtype.lower().startswith("r"):
2353	matchtype = "regex"
2354	else:
2355	raise ValueError("Illegal matchtype, use p(attern) or r(egex)")
2356	sel = selector()
2357	if isinstance(match, selector):
2358	sel = match
2359	else:
2360	sel.set_query("SRCNAME == %s('%s')" % (matchtype, match))
2361	self.set_selection(basesel+sel)
2362	self._setsourcetype(stype)
2363	self.set_selection(basesel)
2364	self._add_history("set_sourcetype", varlist)
2365
2366	@print_log_dec
2367	@preserve_selection
2368	def auto_quotient(self, preserve=True, mode='paired', verify=False):
2369	"""\
2370	This function allows to build quotients automatically.
2371	It assumes the observation to have the same number of
2372	"ons" and "offs"
2373
2374	Parameters:
2375
2376	preserve: you can preserve (default) the continuum or
2377	remove it. The equations used are
2378
2379	preserve: Output = Toff * (on/off) - Toff
2380
2381	remove: Output = Toff * (on/off) - Ton
2382
2383	mode: the on/off detection mode
2384	'paired' (default)
2385	identifies 'off' scans by the
2386	trailing '_R' (Mopra/Parkes) or
2387	'_e'/'_w' (Tid) and matches
2388	on/off pairs from the observing pattern
2389	'time'
2390	finds the closest off in time
2391
2392	.. todo:: verify argument is not implemented
2393
2394	"""
2395	varlist = vars()
2396	modes = ["time", "paired"]
2397	if not mode in modes:
2398	msg = "please provide valid mode. Valid modes are %s" % (modes)
2399	raise ValueError(msg)
2400	s = None
2401	if mode.lower() == "paired":
2402	sel = self.get_selection()
2403	sel.set_query("SRCTYPE==1")
2404	self.set_selection(sel)
2405	offs = self.copy()
2406	sel.set_query("SRCTYPE==0")
2407	self.set_selection(sel)
2408	ons = self.copy()
2409	s = scantable(self._math._quotient(ons, offs, preserve))
2410	elif mode.lower() == "time":
2411	s = scantable(self._math._auto_quotient(self, mode, preserve))
2412	s._add_history("auto_quotient", varlist)
2413	return s
2414
2415	@print_log_dec
2416	def mx_quotient(self, mask = None, weight='median', preserve=True):
2417	"""\
2418	Form a quotient using "off" beams when observing in "MX" mode.
2419
2420	Parameters:
2421
2422	mask: an optional mask to be used when weight == 'stddev'
2423
2424	weight: How to average the off beams. Default is 'median'.
2425
2426	preserve: you can preserve (default) the continuum or
2427	remove it. The equations used are:
2428
2429	preserve: Output = Toff * (on/off) - Toff
2430
2431	remove: Output = Toff * (on/off) - Ton
2432
2433	"""
2434	mask = mask or ()
2435	varlist = vars()
2436	on = scantable(self._math._mx_extract(self, 'on'))
2437	preoff = scantable(self._math._mx_extract(self, 'off'))
2438	off = preoff.average_time(mask=mask, weight=weight, scanav=False)
2439	from asapmath import quotient
2440	q = quotient(on, off, preserve)
2441	q._add_history("mx_quotient", varlist)
2442	return q
2443
2444	@print_log_dec
2445	def freq_switch(self, insitu=None):
2446	"""\
2447	Apply frequency switching to the data.
2448
2449	Parameters:
2450
2451	insitu: if False a new scantable is returned.
2452	Otherwise, the swictching is done in-situ
2453	The default is taken from .asaprc (False)
2454
2455	"""
2456	if insitu is None: insitu = rcParams['insitu']
2457	self._math._setinsitu(insitu)
2458	varlist = vars()
2459	s = scantable(self._math._freqswitch(self))
2460	s._add_history("freq_switch", varlist)
2461	if insitu:
2462	self._assign(s)
2463	else:
2464	return s
2465
2466	@print_log_dec
2467	def recalc_azel(self):
2468	"""Recalculate the azimuth and elevation for each position."""
2469	varlist = vars()
2470	self._recalcazel()
2471	self._add_history("recalc_azel", varlist)
2472	return
2473
2474	@print_log_dec
2475	def __add__(self, other):
2476	varlist = vars()
2477	s = None
2478	if isinstance(other, scantable):
2479	s = scantable(self._math._binaryop(self, other, "ADD"))
2480	elif isinstance(other, float):
2481	s = scantable(self._math._unaryop(self, other, "ADD", False))
2482	else:
2483	raise TypeError("Other input is not a scantable or float value")
2484	s._add_history("operator +", varlist)
2485	return s
2486
2487	@print_log_dec
2488	def __sub__(self, other):
2489	"""
2490	implicit on all axes and on Tsys
2491	"""
2492	varlist = vars()
2493	s = None
2494	if isinstance(other, scantable):
2495	s = scantable(self._math._binaryop(self, other, "SUB"))
2496	elif isinstance(other, float):
2497	s = scantable(self._math._unaryop(self, other, "SUB", False))
2498	else:
2499	raise TypeError("Other input is not a scantable or float value")
2500	s._add_history("operator -", varlist)
2501	return s
2502
2503	@print_log_dec
2504	def __mul__(self, other):
2505	"""
2506	implicit on all axes and on Tsys
2507	"""
2508	varlist = vars()
2509	s = None
2510	if isinstance(other, scantable):
2511	s = scantable(self._math._binaryop(self, other, "MUL"))
2512	elif isinstance(other, float):
2513	s = scantable(self._math._unaryop(self, other, "MUL", False))
2514	else:
2515	raise TypeError("Other input is not a scantable or float value")
2516	s._add_history("operator *", varlist)
2517	return s
2518
2519
2520	@print_log_dec
2521	def __div__(self, other):
2522	"""
2523	implicit on all axes and on Tsys
2524	"""
2525	varlist = vars()
2526	s = None
2527	if isinstance(other, scantable):
2528	s = scantable(self._math._binaryop(self, other, "DIV"))
2529	elif isinstance(other, float):
2530	if other == 0.0:
2531	raise ZeroDivisionError("Dividing by zero is not recommended")
2532	s = scantable(self._math._unaryop(self, other, "DIV", False))
2533	else:
2534	raise TypeError("Other input is not a scantable or float value")
2535	s._add_history("operator /", varlist)
2536	return s
2537
2538	@print_log_dec
2539	def get_fit(self, row=0):
2540	"""\
2541	Print or return the stored fits for a row in the scantable
2542
2543	Parameters:
2544
2545	row: the row which the fit has been applied to.
2546
2547	"""
2548	if row > self.nrow():
2549	return
2550	from asap.asapfit import asapfit
2551	fit = asapfit(self._getfit(row))
2552	if rcParams['verbose']:
2553	#print fit
2554	asaplog.push( '%s' %(fit) )
2555	return
2556	else:
2557	return fit.as_dict()
2558
2559	def flag_nans(self):
2560	"""\
2561	Utility function to flag NaN values in the scantable.
2562	"""
2563	import numpy
2564	basesel = self.get_selection()
2565	for i in range(self.nrow()):
2566	sel = self.get_row_selector(i)
2567	self.set_selection(basesel+sel)
2568	nans = numpy.isnan(self._getspectrum(0))
2569	if numpy.any(nans):
2570	bnans = [ bool(v) for v in nans]
2571	self.flag(bnans)
2572	self.set_selection(basesel)
2573
2574	def get_row_selector(self, rowno):
2575	return selector(beams=self.getbeam(rowno),
2576	ifs=self.getif(rowno),
2577	pols=self.getpol(rowno),
2578	scans=self.getscan(rowno),
2579	cycles=self.getcycle(rowno))
2580
2581	def _add_history(self, funcname, parameters):
2582	if not rcParams['scantable.history']:
2583	return
2584	# create date
2585	sep = "##"
2586	from datetime import datetime
2587	dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
2588	hist = dstr+sep
2589	hist += funcname+sep#cdate+sep
2590	if parameters.has_key('self'): del parameters['self']
2591	for k, v in parameters.iteritems():
2592	if type(v) is dict:
2593	for k2, v2 in v.iteritems():
2594	hist += k2
2595	hist += "="
2596	if isinstance(v2, scantable):
2597	hist += 'scantable'
2598	elif k2 == 'mask':
2599	if isinstance(v2, list) or isinstance(v2, tuple):
2600	hist += str(self._zip_mask(v2))
2601	else:
2602	hist += str(v2)
2603	else:
2604	hist += str(v2)
2605	else:
2606	hist += k
2607	hist += "="
2608	if isinstance(v, scantable):
2609	hist += 'scantable'
2610	elif k == 'mask':
2611	if isinstance(v, list) or isinstance(v, tuple):
2612	hist += str(self._zip_mask(v))
2613	else:
2614	hist += str(v)
2615	else:
2616	hist += str(v)
2617	hist += sep
2618	hist = hist[:-2] # remove trailing '##'
2619	self._addhistory(hist)
2620
2621
2622	def _zip_mask(self, mask):
2623	mask = list(mask)
2624	i = 0
2625	segments = []
2626	while mask[i:].count(1):
2627	i += mask[i:].index(1)
2628	if mask[i:].count(0):
2629	j = i + mask[i:].index(0)
2630	else:
2631	j = len(mask)
2632	segments.append([i, j])
2633	i = j
2634	return segments
2635
2636	def _get_ordinate_label(self):
2637	fu = "("+self.get_fluxunit()+")"
2638	import re
2639	lbl = "Intensity"
2640	if re.match(".K.", fu):
2641	lbl = "Brightness Temperature "+ fu
2642	elif re.match(".Jy.", fu):
2643	lbl = "Flux density "+ fu
2644	return lbl
2645
2646	def _check_ifs(self):
2647	nchans = [self.nchan(i) for i in range(self.nif(-1))]
2648	nchans = filter(lambda t: t > 0, nchans)
2649	return (sum(nchans)/len(nchans) == nchans[0])
2650
2651	def _fill(self, names, unit, average, getpt, antenna):
2652	first = True
2653	fullnames = []
2654	for name in names:
2655	name = os.path.expandvars(name)
2656	name = os.path.expanduser(name)
2657	if not os.path.exists(name):
2658	msg = "File '%s' does not exists" % (name)
2659	if rcParams['verbose']:
2660	asaplog.push(msg)
2661	print_log( 'ERROR' )
2662	return
2663	raise IOError(msg)
2664	fullnames.append(name)
2665	if average:
2666	asaplog.push('Auto averaging integrations')
2667	stype = int(rcParams['scantable.storage'].lower() == 'disk')
2668	for name in fullnames:
2669	tbl = Scantable(stype)
2670	r = filler(tbl)
2671	rx = rcParams['scantable.reference']
2672	r.setreferenceexpr(rx)
2673	msg = "Importing %s..." % (name)
2674	asaplog.push(msg, False)
2675	print_log()
2676	r.open(name)# antenna, -1, -1, getpt)
2677	r.fill()
2678	if average:
2679	tbl = self._math._average((tbl, ), (), 'NONE', 'SCAN')
2680	if not first:
2681	tbl = self._math._merge([self, tbl])
2682	Scantable.__init__(self, tbl)
2683	r.close()
2684	del r, tbl
2685	first = False
2686	if unit is not None:
2687	self.set_fluxunit(unit)
2688	if not is_casapy():
2689	self.set_freqframe(rcParams['scantable.freqframe'])
2690
2691	def __getitem__(self, key):
2692	if key < 0:
2693	key += self.nrow()
2694	if key >= self.nrow():
2695	raise IndexError("Row index out of range.")
2696	return self._getspectrum(key)
2697
2698	def __setitem__(self, key, value):
2699	if key < 0:
2700	key += self.nrow()
2701	if key >= self.nrow():
2702	raise IndexError("Row index out of range.")
2703	if not hasattr(value, "__len__") or \
2704	len(value) > self.nchan(self.getif(key)):
2705	raise ValueError("Spectrum length doesn't match.")
2706	return self._setspectrum(value, key)
2707
2708	def __len__(self):
2709	return self.nrow()
2710
2711	def __iter__(self):
2712	for i in range(len(self)):
2713	yield self[i]

Note: See TracBrowser for help on using the repository browser.

Download in other formats: