Context Navigation

source: trunk/python/scantable.py @ 1574

Visit:

Last change on this file since 1574 was 1574, checked in by Malte Marquarding, 15 years ago
Ticket #167: python part of running polynomial smoothing
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 70.4 KB

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

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

Download in other formats: