Context Navigation

source: branches/alma/python/scantable.py @ 1701

Visit:

Last change on this file since 1701 was 1701, checked in by WataruKawasaki, 14 years ago

New Development: Yes

JIRA Issue: Yes (CAS-1800 + CAS-1807)

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: added new methods to scantable and fitter.

Test Programs:

Put in Release Notes: No

Module(s): sdfit, sdflag

Description: added new methods 'scantable.clip' and 'fitter.set_lorentz_parameters'.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 87.0 KB

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

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

Download in other formats: