Context Navigation

scantable.py @ 1613

Visit:

Last change on this file since 1613 was 1613, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No?

Module(s): Module Names change impacts.

Description: Describe your changes here...

Temporary file name has been changed.

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

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

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

Context Navigation

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

Download in other formats: