Context Navigation

source: trunk/python/scantable.py @ 1554

Visit:

Last change on this file since 1554 was 1554, checked in by Malte Marquarding, 15 years ago
Support maks windows of style [channel]
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 70.1 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	Parameters:
1217	poltype: The new polarisation type. Valid types are:
1218	"linear", "stokes" and "circular"
1219	"""
1220	varlist = vars()
1221	try:
1222	s = scantable(self._math._convertpol(self, poltype))
1223	except RuntimeError, msg:
1224	if rcParams['verbose']:
1225	print msg
1226	return
1227	else:
1228	raise
1229	s._add_history("convert_pol", varlist)
1230	print_log()
1231	return s
1232
1233	def smooth(self, kernel="hanning", width=5.0, insitu=None):
1234	"""
1235	Smooth the spectrum by the specified kernel (conserving flux).
1236	Parameters:
1237	kernel: The type of smoothing kernel. Select from
1238	'hanning' (default), 'gaussian', 'boxcar' and
1239	'rmedian'
1240	width: The width of the kernel in pixels. For hanning this is
1241	ignored otherwise it defauls to 5 pixels.
1242	For 'gaussian' it is the Full Width Half
1243	Maximum. For 'boxcar' it is the full width.
1244	For 'rmedian' it is the half width.
1245	insitu: if False a new scantable is returned.
1246	Otherwise, the scaling is done in-situ
1247	The default is taken from .asaprc (False)
1248	Example:
1249	none
1250	"""
1251	if insitu is None: insitu = rcParams['insitu']
1252	self._math._setinsitu(insitu)
1253	varlist = vars()
1254	s = scantable(self._math._smooth(self, kernel.lower(), width))
1255	s._add_history("smooth", varlist)
1256	print_log()
1257	if insitu: self._assign(s)
1258	else: return s
1259
1260
1261	def poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None):
1262	"""
1263	Return a scan which has been baselined (all rows) by a polynomial.
1264	Parameters:
1265	mask: an optional mask
1266	order: the order of the polynomial (default is 0)
1267	plot: plot the fit and the residual. In this each
1268	indivual fit has to be approved, by typing 'y'
1269	or 'n'
1270	uselin: use linear polynomial fit
1271	insitu: if False a new scantable is returned.
1272	Otherwise, the scaling is done in-situ
1273	The default is taken from .asaprc (False)
1274	Example:
1275	# return a scan baselined by a third order polynomial,
1276	# not using a mask
1277	bscan = scan.poly_baseline(order=3)
1278	"""
1279	if insitu is None: insitu = rcParams['insitu']
1280	varlist = vars()
1281	if mask is None:
1282	mask = [True for i in xrange(self.nchan(-1))]
1283	from asap.asapfitter import fitter
1284	try:
1285	f = fitter()
1286	f.set_scan(self, mask)
1287	#f.set_function(poly=order)
1288	if uselin:
1289	f.set_function(lpoly=order)
1290	else:
1291	f.set_function(poly=order)
1292	s = f.auto_fit(insitu, plot=plot)
1293	s._add_history("poly_baseline", varlist)
1294	print_log()
1295	if insitu: self._assign(s)
1296	else: return s
1297	except RuntimeError:
1298	msg = "The fit failed, possibly because it didn't converge."
1299	if rcParams['verbose']:
1300	print msg
1301	return
1302	else:
1303	raise RuntimeError(msg)
1304
1305
1306	def auto_poly_baseline(self, mask=[], edge=(0, 0), order=0,
1307	threshold=3, chan_avg_limit=1, plot=False,
1308	insitu=None):
1309	"""
1310	Return a scan which has been baselined (all rows) by a polynomial.
1311	Spectral lines are detected first using linefinder and masked out
1312	to avoid them affecting the baseline solution.
1313
1314	Parameters:
1315	mask: an optional mask retreived from scantable
1316	edge: an optional number of channel to drop at
1317	the edge of spectrum. If only one value is
1318	specified, the same number will be dropped from
1319	both sides of the spectrum. Default is to keep
1320	all channels. Nested tuples represent individual
1321	edge selection for different IFs (a number of spectral
1322	channels can be different)
1323	order: the order of the polynomial (default is 0)
1324	threshold: the threshold used by line finder. It is better to
1325	keep it large as only strong lines affect the
1326	baseline solution.
1327	chan_avg_limit:
1328	a maximum number of consequtive spectral channels to
1329	average during the search of weak and broad lines.
1330	The default is no averaging (and no search for weak
1331	lines). If such lines can affect the fitted baseline
1332	(e.g. a high order polynomial is fitted), increase this
1333	parameter (usually values up to 8 are reasonable). Most
1334	users of this method should find the default value
1335	sufficient.
1336	plot: plot the fit and the residual. In this each
1337	indivual fit has to be approved, by typing 'y'
1338	or 'n'
1339	insitu: if False a new scantable is returned.
1340	Otherwise, the scaling is done in-situ
1341	The default is taken from .asaprc (False)
1342
1343	Example:
1344	scan2=scan.auto_poly_baseline(order=7)
1345	"""
1346	if insitu is None: insitu = rcParams['insitu']
1347	varlist = vars()
1348	from asap.asapfitter import fitter
1349	from asap.asaplinefind import linefinder
1350	from asap import _is_sequence_or_number as _is_valid
1351
1352	# check whether edge is set up for each IF individually
1353	individualedge = False;
1354	if len(edge) > 1:
1355	if isinstance(edge[0], list) or isinstance(edge[0], tuple):
1356	individualedge = True;
1357
1358	if not _is_valid(edge, int) and not individualedge:
1359	raise ValueError, "Parameter 'edge' has to be an integer or a \
1360	pair of integers specified as a tuple. Nested tuples are allowed \
1361	to make individual selection for different IFs."
1362
1363	curedge = (0, 0)
1364	if individualedge:
1365	for edgepar in edge:
1366	if not _is_valid(edgepar, int):
1367	raise ValueError, "Each element of the 'edge' tuple has \
1368	to be a pair of integers or an integer."
1369	else:
1370	curedge = edge;
1371
1372	# setup fitter
1373	f = fitter()
1374	f.set_function(poly=order)
1375
1376	# setup line finder
1377	fl = linefinder()
1378	fl.set_options(threshold=threshold,avg_limit=chan_avg_limit)
1379
1380	if not insitu:
1381	workscan = self.copy()
1382	else:
1383	workscan = self
1384
1385	fl.set_scan(workscan)
1386
1387	rows = range(workscan.nrow())
1388	asaplog.push("Processing:")
1389	for r in rows:
1390	msg = " Scan[%d] Beam[%d] IF[%d] Pol[%d] Cycle[%d]" % \
1391	(workscan.getscan(r), workscan.getbeam(r), workscan.getif(r), \
1392	workscan.getpol(r), workscan.getcycle(r))
1393	asaplog.push(msg, False)
1394
1395	# figure out edge parameter
1396	if individualedge:
1397	if len(edge) >= workscan.getif(r):
1398	raise RuntimeError, "Number of edge elements appear to " \
1399	"be less than the number of IFs"
1400	curedge = edge[workscan.getif(r)]
1401
1402	# setup line finder
1403	fl.find_lines(r, mask, curedge)
1404	f.set_data(workscan._getabcissa(r), workscan._getspectrum(r),
1405	mask_and(workscan._getmask(r), fl.get_mask()))
1406	f.fit()
1407	x = f.get_parameters()
1408	if plot:
1409	f.plot(residual=True)
1410	x = raw_input("Accept fit ( [y]/n ): ")
1411	if x.upper() == 'N':
1412	continue
1413	workscan._setspectrum(f.fitter.getresidual(), r)
1414	if plot:
1415	f._p.unmap()
1416	f._p = None
1417	workscan._add_history("auto_poly_baseline", varlist)
1418	if insitu:
1419	self._assign(workscan)
1420	else:
1421	return workscan
1422
1423	def rotate_linpolphase(self, angle):
1424	"""
1425	Rotate the phase of the complex polarization O=Q+iU correlation.
1426	This is always done in situ in the raw data. So if you call this
1427	function more than once then each call rotates the phase further.
1428	Parameters:
1429	angle: The angle (degrees) to rotate (add) by.
1430	Examples:
1431	scan.rotate_linpolphase(2.3)
1432	"""
1433	varlist = vars()
1434	self._math._rotate_linpolphase(self, angle)
1435	self._add_history("rotate_linpolphase", varlist)
1436	print_log()
1437	return
1438
1439
1440	def rotate_xyphase(self, angle):
1441	"""
1442	Rotate the phase of the XY correlation. This is always done in situ
1443	in the data. So if you call this function more than once
1444	then each call rotates the phase further.
1445	Parameters:
1446	angle: The angle (degrees) to rotate (add) by.
1447	Examples:
1448	scan.rotate_xyphase(2.3)
1449	"""
1450	varlist = vars()
1451	self._math._rotate_xyphase(self, angle)
1452	self._add_history("rotate_xyphase", varlist)
1453	print_log()
1454	return
1455
1456	def swap_linears(self):
1457	"""
1458	Swap the linear polarisations XX and YY, or better the first two
1459	polarisations as this also works for ciculars.
1460	"""
1461	varlist = vars()
1462	self._math._swap_linears(self)
1463	self._add_history("swap_linears", varlist)
1464	print_log()
1465	return
1466
1467	def invert_phase(self):
1468	"""
1469	Invert the phase of the complex polarisation
1470	"""
1471	varlist = vars()
1472	self._math._invert_phase(self)
1473	self._add_history("invert_phase", varlist)
1474	print_log()
1475	return
1476
1477	def add(self, offset, insitu=None):
1478	"""
1479	Return a scan where all spectra have the offset added
1480	Parameters:
1481	offset: the offset
1482	insitu: if False a new scantable is returned.
1483	Otherwise, the scaling is done in-situ
1484	The default is taken from .asaprc (False)
1485	"""
1486	if insitu is None: insitu = rcParams['insitu']
1487	self._math._setinsitu(insitu)
1488	varlist = vars()
1489	s = scantable(self._math._unaryop(self, offset, "ADD", False))
1490	s._add_history("add", varlist)
1491	print_log()
1492	if insitu:
1493	self._assign(s)
1494	else:
1495	return s
1496
1497	def scale(self, factor, tsys=True, insitu=None):
1498	"""
1499	Return a scan where all spectra are scaled by the give 'factor'
1500	Parameters:
1501	factor: the scaling factor
1502	insitu: if False a new scantable is returned.
1503	Otherwise, the scaling is done in-situ
1504	The default is taken from .asaprc (False)
1505	tsys: if True (default) then apply the operation to Tsys
1506	as well as the data
1507	"""
1508	if insitu is None: insitu = rcParams['insitu']
1509	self._math._setinsitu(insitu)
1510	varlist = vars()
1511	s = scantable(self._math._unaryop(self, factor, "MUL", tsys))
1512	s._add_history("scale", varlist)
1513	print_log()
1514	if insitu:
1515	self._assign(s)
1516	else:
1517	return s
1518
1519	def set_sourcetype(self, match, matchtype="pattern",
1520	sourcetype="reference"):
1521	"""
1522	Set the type of the source to be an source or reference scan
1523	using the provided pattern:
1524	Parameters:
1525	match: a Unix style pattern, regular expression or selector
1526	matchtype: 'pattern' (default) UNIX style pattern or
1527	'regex' regular expression
1528	sourcetype: the type of the source to use (source/reference)
1529	"""
1530	varlist = vars()
1531	basesel = self.get_selection()
1532	stype = -1
1533	if sourcetype.lower().startswith("r"):
1534	stype = 1
1535	elif sourcetype.lower().startswith("s"):
1536	stype = 0
1537	else:
1538	raise ValueError("Illegal sourcetype use s(ource) or r(eference)")
1539	if matchtype.lower().startswith("p"):
1540	matchtype = "pattern"
1541	elif matchtype.lower().startswith("r"):
1542	matchtype = "regex"
1543	else:
1544	raise ValueError("Illegal matchtype, use p(attern) or r(egex)")
1545	sel = selector()
1546	if isinstance(match, selector):
1547	sel = match
1548	else:
1549	sel.set_query("SRCNAME == %s('%s')" % (matchtype, match))
1550	self.set_selection(basesel+sel)
1551	self._setsourcetype(stype)
1552	self.set_selection(basesel)
1553	s._add_history("set_sourcetype", varlist)
1554
1555	def auto_quotient(self, preserve=True, mode='paired'):
1556	"""
1557	This function allows to build quotients automatically.
1558	It assumes the observation to have the same numer of
1559	"ons" and "offs"
1560	Parameters:
1561	preserve: you can preserve (default) the continuum or
1562	remove it. The equations used are
1563	preserve: Output = Toff * (on/off) - Toff
1564	remove: Output = Toff * (on/off) - Ton
1565	mode: the on/off detection mode
1566	'paired' (default)
1567	identifies 'off' scans by the
1568	trailing '_R' (Mopra/Parkes) or
1569	'_e'/'_w' (Tid) and matches
1570	on/off pairs from the observing pattern
1571	'time'
1572	finds the closest off in time
1573
1574	"""
1575	modes = ["time", "paired"]
1576	if not mode in modes:
1577	msg = "please provide valid mode. Valid modes are %s" % (modes)
1578	raise ValueError(msg)
1579	varlist = vars()
1580	s = None
1581	if mode.lower() == "paired":
1582	basesel = self.get_selection()
1583	sel = selector()+basesel
1584	sel.set_query("SRCTYPE==1")
1585	self.set_selection(sel)
1586	offs = self.copy()
1587	sel.set_query("SRCTYPE==0")
1588	self.set_selection(sel)
1589	ons = self.copy()
1590	s = scantable(self._math._quotient(ons, offs, preserve))
1591	self.set_selection(basesel)
1592	elif mode.lower() == "time":
1593	s = scantable(self._math._auto_quotient(self, mode, preserve))
1594	s._add_history("auto_quotient", varlist)
1595	print_log()
1596	return s
1597
1598	def mx_quotient(self, mask = None, weight='median', preserve=True):
1599	"""
1600	Form a quotient using "off" beams when observing in "MX" mode.
1601	Parameters:
1602	mask: an optional mask to be used when weight == 'stddev'
1603	weight: How to average the off beams. Default is 'median'.
1604	preserve: you can preserve (default) the continuum or
1605	remove it. The equations used are
1606	preserve: Output = Toff * (on/off) - Toff
1607	remove: Output = Toff * (on/off) - Ton
1608	"""
1609	if mask is None: mask = ()
1610	varlist = vars()
1611	on = scantable(self._math._mx_extract(self, 'on'))
1612	preoff = scantable(self._math._mx_extract(self, 'off'))
1613	off = preoff.average_time(mask=mask, weight=weight, scanav=False)
1614	from asapmath import quotient
1615	q = quotient(on, off, preserve)
1616	q._add_history("mx_quotient", varlist)
1617	print_log()
1618	return q
1619
1620	def freq_switch(self, insitu=None):
1621	"""
1622	Apply frequency switching to the data.
1623	Parameters:
1624	insitu: if False a new scantable is returned.
1625	Otherwise, the swictching is done in-situ
1626	The default is taken from .asaprc (False)
1627	Example:
1628	none
1629	"""
1630	if insitu is None: insitu = rcParams['insitu']
1631	self._math._setinsitu(insitu)
1632	varlist = vars()
1633	s = scantable(self._math._freqswitch(self))
1634	s._add_history("freq_switch", varlist)
1635	print_log()
1636	if insitu: self._assign(s)
1637	else: return s
1638
1639	def recalc_azel(self):
1640	"""
1641	Recalculate the azimuth and elevation for each position.
1642	Parameters:
1643	none
1644	Example:
1645	"""
1646	varlist = vars()
1647	self._recalcazel()
1648	self._add_history("recalc_azel", varlist)
1649	print_log()
1650	return
1651
1652	def __add__(self, other):
1653	varlist = vars()
1654	s = None
1655	if isinstance(other, scantable):
1656	s = scantable(self._math._binaryop(self, other, "ADD"))
1657	elif isinstance(other, float):
1658	s = scantable(self._math._unaryop(self, other, "ADD", False))
1659	else:
1660	raise TypeError("Other input is not a scantable or float value")
1661	s._add_history("operator +", varlist)
1662	print_log()
1663	return s
1664
1665	def __sub__(self, other):
1666	"""
1667	implicit on all axes and on Tsys
1668	"""
1669	varlist = vars()
1670	s = None
1671	if isinstance(other, scantable):
1672	s = scantable(self._math._binaryop(self, other, "SUB"))
1673	elif isinstance(other, float):
1674	s = scantable(self._math._unaryop(self, other, "SUB", False))
1675	else:
1676	raise TypeError("Other input is not a scantable or float value")
1677	s._add_history("operator -", varlist)
1678	print_log()
1679	return s
1680
1681	def __mul__(self, other):
1682	"""
1683	implicit on all axes and on Tsys
1684	"""
1685	varlist = vars()
1686	s = None
1687	if isinstance(other, scantable):
1688	s = scantable(self._math._binaryop(self, other, "MUL"))
1689	elif isinstance(other, float):
1690	s = scantable(self._math._unaryop(self, other, "MUL", False))
1691	else:
1692	raise TypeError("Other input is not a scantable or float value")
1693	s._add_history("operator *", varlist)
1694	print_log()
1695	return s
1696
1697
1698	def __div__(self, other):
1699	"""
1700	implicit on all axes and on Tsys
1701	"""
1702	varlist = vars()
1703	s = None
1704	if isinstance(other, scantable):
1705	s = scantable(self._math._binaryop(self, other, "DIV"))
1706	elif isinstance(other, float):
1707	if other == 0.0:
1708	raise ZeroDivisionError("Dividing by zero is not recommended")
1709	s = scantable(self._math._unaryop(self, other, "DIV", False))
1710	else:
1711	raise TypeError("Other input is not a scantable or float value")
1712	s._add_history("operator /", varlist)
1713	print_log()
1714	return s
1715
1716	def get_fit(self, row=0):
1717	"""
1718	Print or return the stored fits for a row in the scantable
1719	Parameters:
1720	row: the row which the fit has been applied to.
1721	"""
1722	if row > self.nrow():
1723	return
1724	from asap.asapfit import asapfit
1725	fit = asapfit(self._getfit(row))
1726	if rcParams['verbose']:
1727	print fit
1728	return
1729	else:
1730	return fit.as_dict()
1731
1732	def flag_nans(self):
1733	"""
1734	Utility function to flag NaN values in the scantable.
1735	"""
1736	import numpy
1737	basesel = self.get_selection()
1738	for i in range(self.nrow()):
1739	sel = selector()+basesel
1740	sel.set_scans(self.getscan(i))
1741	sel.set_beams(self.getbeam(i))
1742	sel.set_ifs(self.getif(i))
1743	sel.set_polarisations(self.getpol(i))
1744	self.set_selection(sel)
1745	nans = numpy.isnan(self._getspectrum(0))
1746	if numpy.any(nans):
1747	bnans = [ bool(v) for v in nans]
1748	self.flag(bnans)
1749	self.set_selection(basesel)
1750
1751
1752	def _add_history(self, funcname, parameters):
1753	if not rcParams['scantable.history']:
1754	return
1755	# create date
1756	sep = "##"
1757	from datetime import datetime
1758	dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
1759	hist = dstr+sep
1760	hist += funcname+sep#cdate+sep
1761	if parameters.has_key('self'): del parameters['self']
1762	for k, v in parameters.iteritems():
1763	if type(v) is dict:
1764	for k2, v2 in v.iteritems():
1765	hist += k2
1766	hist += "="
1767	if isinstance(v2, scantable):
1768	hist += 'scantable'
1769	elif k2 == 'mask':
1770	if isinstance(v2, list) or isinstance(v2, tuple):
1771	hist += str(self._zip_mask(v2))
1772	else:
1773	hist += str(v2)
1774	else:
1775	hist += str(v2)
1776	else:
1777	hist += k
1778	hist += "="
1779	if isinstance(v, scantable):
1780	hist += 'scantable'
1781	elif k == 'mask':
1782	if isinstance(v, list) or isinstance(v, tuple):
1783	hist += str(self._zip_mask(v))
1784	else:
1785	hist += str(v)
1786	else:
1787	hist += str(v)
1788	hist += sep
1789	hist = hist[:-2] # remove trailing '##'
1790	self._addhistory(hist)
1791
1792
1793	def _zip_mask(self, mask):
1794	mask = list(mask)
1795	i = 0
1796	segments = []
1797	while mask[i:].count(1):
1798	i += mask[i:].index(1)
1799	if mask[i:].count(0):
1800	j = i + mask[i:].index(0)
1801	else:
1802	j = len(mask)
1803	segments.append([i, j])
1804	i = j
1805	return segments
1806
1807	def _get_ordinate_label(self):
1808	fu = "("+self.get_fluxunit()+")"
1809	import re
1810	lbl = "Intensity"
1811	if re.match(".K.", fu):
1812	lbl = "Brightness Temperature "+ fu
1813	elif re.match(".Jy.", fu):
1814	lbl = "Flux density "+ fu
1815	return lbl
1816
1817	def _check_ifs(self):
1818	nchans = [self.nchan(i) for i in range(self.nif(-1))]
1819	nchans = filter(lambda t: t > 0, nchans)
1820	return (sum(nchans)/len(nchans) == nchans[0])
1821
1822	def _fill(self, names, unit, average):
1823	import os
1824	from asap._asap import stfiller
1825	first = True
1826	fullnames = []
1827	for name in names:
1828	name = os.path.expandvars(name)
1829	name = os.path.expanduser(name)
1830	if not os.path.exists(name):
1831	msg = "File '%s' does not exists" % (name)
1832	if rcParams['verbose']:
1833	asaplog.push(msg)
1834	print asaplog.pop().strip()
1835	return
1836	raise IOError(msg)
1837	fullnames.append(name)
1838	if average:
1839	asaplog.push('Auto averaging integrations')
1840	stype = int(rcParams['scantable.storage'].lower() == 'disk')
1841	for name in fullnames:
1842	tbl = Scantable(stype)
1843	r = stfiller(tbl)
1844	rx = rcParams['scantable.reference']
1845	r._setreferenceexpr(rx)
1846	msg = "Importing %s..." % (name)
1847	asaplog.push(msg, False)
1848	print_log()
1849	r._open(name, -1, -1)
1850	r._read()
1851	if average:
1852	tbl = self._math._average((tbl, ), (), 'NONE', 'SCAN')
1853	if not first:
1854	tbl = self._math._merge([self, tbl])
1855	Scantable.__init__(self, tbl)
1856	r._close()
1857	del r, tbl
1858	first = False
1859	if unit is not None:
1860	self.set_fluxunit(unit)
1861	self.set_freqframe(rcParams['scantable.freqframe'])
1862
1863	def __getitem__(self, key):
1864	if key < 0:
1865	key += self.nrow()
1866	if key >= self.nrow():
1867	raise IndexError("Row index out of range.")
1868	return self._getspectrum(key)
1869
1870	def __setitem__(self, key, value):
1871	if key < 0:
1872	key += self.nrow()
1873	if key >= self.nrow():
1874	raise IndexError("Row index out of range.")
1875	if not hasattr(value, "__len__") or \
1876	len(value) > self.nchan(self.getif(key)):
1877	raise ValueError("Spectrum length doesn't match.")
1878	return self._setspectrum(value, key)
1879
1880	def __len__(self):
1881	return self.nrow()
1882
1883	def __iter__(self):
1884	for i in range(len(self)):
1885	yield self[i]

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

Download in other formats: