Context Navigation

scantable.py @ 1612

Visit:

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

New Development: No

JIRA Issue: Yes CAS-729, CAS-1147

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes

Module(s): Module Names change impacts.

Description: Describe your changes here...

I have changed that almost all log messages are output to casapy.log,
not to the terminal window. After this change, asap becomes to depend on casapy
and is not running in standalone, because asap have to import taskinit module
to access casalogger.

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

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

Context Navigation

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

Download in other formats: