Context Navigation

source: trunk/python/scantable.py @ 451

Visit:

Last change on this file since 451 was 451, checked in by kil064, 19 years ago
make default for stokes=False in function save make format cas insensitive in function save
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 26.3 KB

Line
1	from asap._asap import sdtable
2	from asap import rcParams
3	from numarray import ones,zeros
4	import sys
5
6	class scantable(sdtable):
7	"""
8	The ASAP container for scans
9	"""
10
11	def __init__(self, filename, unit=None):
12	"""
13	Create a scantable from a saved one or make a reference
14	Parameters:
15	filename: the name of an asap table on disk
16	or
17	the name of a rpfits/sdfits/ms file
18	(integrations within scans are auto averaged
19	and the whole file is read)
20	or
21	[advanced] a reference to an existing
22	scantable
23	unit: brightness unit; must be consistent with K or Jy.
24	Over-rides the default selected by the reader
25	(input rpfits/sdfits/ms) or replaces the value
26	in existing scantables
27	"""
28	self._vb = rcParams['verbose']
29	self._p = None
30	from os import stat as st
31	import stat
32	if isinstance(filename,sdtable):
33	sdtable.__init__(self, filename)
34	if unit is not None:
35	self.set_fluxunit(unit)
36	else:
37	import os.path
38	if not os.path.exists(filename):
39	print "File '%s' not found." % (filename)
40	return
41	filename = os.path.expandvars(filename)
42	if os.path.isdir(filename):
43	# crude check if asap table
44	if os.path.exists(filename+'/table.info'):
45	sdtable.__init__(self, filename)
46	if unit is not None:
47	self.set_fluxunit(unit)
48	else:
49	print "The given file '%s'is not a valid asap table." % (filename)
50	return
51	else:
52	autoav = rcParams['scantable.autoaverage']
53
54	from asap._asap import sdreader
55	ifSel = -1
56	beamSel = -1
57	r = sdreader(filename,ifSel,beamSel)
58	print 'Importing data...'
59	r._read([-1])
60	tbl = r._getdata()
61	if unit is not None:
62	tbl.set_fluxunit(unit)
63	if autoav:
64	from asap._asap import average
65	tmp = tuple([tbl])
66	print 'Auto averaging integrations...'
67	tbl2 = average(tmp,(),True,'none')
68	sdtable.__init__(self,tbl2)
69	del r,tbl
70	else:
71	sdtable.__init__(self,tbl)
72
73	def save(self, name=None, format=None, stokes=False, overwrite=False):
74	"""
75	Store the scantable on disk. This can be an asap (aips++) Table, SDFITS,
76	Image FITS or MS2 format.
77	Parameters:
78	name: the name of the outputfile. For format="FITS" this
79	is the directory file name into which all the files will
80	be written (default is 'asap_FITS')
81	format: an optional file format. Default is ASAP.
82	Allowed are - 'ASAP' (save as ASAP [aips++] Table),
83	'SDFITS' (save as SDFITS file)
84	'FITS' (saves each row as a FITS Image)
85	'ASCII' (saves as ascii text file)
86	'MS2' (saves as an aips++
87	MeasurementSet V2)
88	stokes: Convert to Stokes parameters (only available
89	currently with FITS and ASCII formats.
90	Default is False.
91	overwrite: If the file should be overwritten if it exists.
92	The default False is to return with warning
93	without writing the output. USE WITH CARE.
94	Example:
95	scan.save('myscan.asap')
96	scan.save('myscan.sdfits','SDFITS')
97	"""
98	from os import path
99	if format is None: format = rcParams['scantable.save']
100	suffix = '.'+format.lower()
101	if name is None or name =="":
102	name = 'scantable'+suffix
103	print "No filename given. Using default name %s..." % name
104	name = path.expandvars(name)
105	if path.isfile(name) or path.isdir(name):
106	if not overwrite:
107	print "File %s already exists." % name
108	return
109	format2 = format.upper()
110	if format2 == 'ASAP':
111	self._save(name)
112	else:
113	from asap._asap import sdwriter as _sw
114	w = _sw(format2)
115	w.write(self, name, stokes)
116	return
117
118	def copy(self):
119	"""
120	Return a copy of this scantable.
121	Parameters:
122	none
123	Example:
124	copiedscan = scan.copy()
125	"""
126	sd = scantable(sdtable._copy(self))
127	return sd
128
129	def get_scan(self, scanid=None):
130	"""
131	Return a specific scan (by scanno) or collection of scans (by
132	source name) in a new scantable.
133	Parameters:
134	scanid: a scanno or a source name
135	Example:
136	scan.get_scan('323p459')
137	# gets all scans containing the source '323p459'
138	"""
139	if scanid is None:
140	print "Please specify a scan no or name to retrieve from the scantable"
141	try:
142	if type(scanid) is str:
143	s = sdtable._getsource(self,scanid)
144	return scantable(s)
145	elif type(scanid) is int:
146	s = sdtable._getscan(self,[scanid])
147	return scantable(s)
148	elif type(scanid) is list:
149	s = sdtable._getscan(self,scanid)
150	return scantable(s)
151	else:
152	print "Illegal scanid type, use 'int' or 'list' if ints."
153	except RuntimeError:
154	print "Couldn't find any match."
155
156	def __str__(self):
157	return sdtable._summary(self,True)
158
159	def summary(self,filename=None, verbose=None):
160	"""
161	Print a summary of the contents of this scantable.
162	Parameters:
163	filename: the name of a file to write the putput to
164	Default - no file output
165	verbose: print extra info such as the frequency table
166	The default (False) is taken from .asaprc
167	"""
168	info = sdtable._summary(self, verbose)
169	if verbose is None: verbose = rcParams['scantable.verbosesummary']
170	if filename is not None:
171	if filename is "":
172	filename = 'scantable_summary.txt'
173	from os.path import expandvars, isdir
174	filename = expandvars(filename)
175	if not isdir(filename):
176	data = open(filename, 'w')
177	data.write(info)
178	data.close()
179	else:
180	print "Illegal file name '%s'." % (filename)
181	print info
182
183	def set_cursor(self, thebeam=0,theif=0,thepol=0):
184	"""
185	Set the spectrum for individual operations.
186	Parameters:
187	thebeam,theif,thepol: a number
188	Example:
189	scan.set_cursor(0,0,1)
190	pol1sig = scan.stats(all=False) # returns std dev for beam=0
191	# if=0, pol=1
192	"""
193	self.setbeam(thebeam)
194	self.setpol(thepol)
195	self.setif(theif)
196	return
197
198	def get_cursor(self):
199	"""
200	Return/print a the current 'cursor' into the Beam/IF/Pol cube.
201	Parameters:
202	none
203	Returns:
204	a list of values (currentBeam,currentIF,currentPol)
205	Example:
206	none
207	"""
208	i = self.getbeam()
209	j = self.getif()
210	k = self.getpol()
211	if self._vb:
212	print "--------------------------------------------------"
213	print " Cursor position"
214	print "--------------------------------------------------"
215	out = 'Beam=%d IF=%d Pol=%d ' % (i,j,k)
216	print out
217	return i,j,k
218
219	def stats(self, stat='stddev', mask=None, allaxes=None):
220	"""
221	Determine the specified statistic of the current beam/if/pol
222	Takes a 'mask' as an optional parameter to specify which
223	channels should be excluded.
224	Parameters:
225	stat: 'min', 'max', 'sumsq', 'sum', 'mean'
226	'var', 'stddev', 'avdev', 'rms', 'median'
227	mask: an optional mask specifying where the statistic
228	should be determined.
229	allaxes: if True apply to all spectra. Otherwise
230	apply only to the selected (beam/pol/if)spectra only.
231	The default is taken from .asaprc (True if none)
232	Example:
233	scan.set_unit('channel')
234	msk = scan.create_mask([100,200],[500,600])
235	scan.stats(stat='mean', mask=m)
236	"""
237	if allaxes is None: allaxes = rcParams['scantable.allaxes']
238	from asap._asap import stats as _stats
239	from numarray import array,zeros,Float
240	if mask == None:
241	mask = ones(self.nchan())
242	axes = ['Beam','IF','Pol','Time']
243
244	beamSel,IFSel,polSel = (self.getbeam(),self.getif(),self.getpol())
245	if allaxes:
246	n = self.nbeam()self.nif()self.npol()*self.nrow()
247	shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
248	arr = array(zeros(n),shape=shp,type=Float)
249
250	for i in range(self.nbeam()):
251	self.setbeam(i)
252	for j in range(self.nif()):
253	self.setif(j)
254	for k in range(self.npol()):
255	self.setpol(k)
256	arr[i,j,k,:] = _stats(self,mask,stat,-1)
257	retval = {'axes': axes, 'data': arr, 'cursor':None}
258	tm = [self._gettime(val) for val in range(self.nrow())]
259	if self._vb:
260	self._print_values(retval,stat,tm)
261	self.setbeam(beamSel)
262	self.setif(IFSel)
263	self.setpol(polSel)
264	return retval
265
266	else:
267	statval = _stats(self,mask,stat,-1)
268	out = ''
269	for l in range(self.nrow()):
270	tm = self._gettime(l)
271	out += 'Time[%s]:\n' % (tm)
272	if self.nbeam() > 1: out += ' Beam[%d] ' % (beamSel)
273	if self.nif() > 1: out += ' IF[%d] ' % (IFSel)
274	if self.npol() > 1: out += ' Pol[%d] ' % (polSel)
275	out += '= %3.3f\n' % (statval[l])
276	out += "--------------------------------------------------\n"
277
278	if self._vb:
279	print "--------------------------------------------------"
280	print " ",stat
281	print "--------------------------------------------------"
282	print out
283	retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
284	return retval
285
286	def stddev(self,mask=None, allaxes=None):
287	"""
288	Determine the standard deviation of the current beam/if/pol
289	Takes a 'mask' as an optional parameter to specify which
290	channels should be excluded.
291	Parameters:
292	mask: an optional mask specifying where the standard
293	deviation should be determined.
294	allaxes: optional flag to show all or a cursor selected
295	spectrum of Beam/IF/Pol. Default is all or taken
296	from .asaprc
297
298	Example:
299	scan.set_unit('channel')
300	msk = scan.create_mask([100,200],[500,600])
301	scan.stddev(mask=m)
302	"""
303	if allaxes is None: allaxes = rcParams['scantable.allaxes']
304	return self.stats(stat='stddev',mask=mask, allaxes=allaxes);
305
306	def get_tsys(self, allaxes=None):
307	"""
308	Return the System temperatures.
309	Parameters:
310	allaxes: if True apply to all spectra. Otherwise
311	apply only to the selected (beam/pol/if)spectra only.
312	The default is taken from .asaprc (True if none)
313	Returns:
314	a list of Tsys values.
315	"""
316	if allaxes is None: allaxes = rcParams['scantable.allaxes']
317	from numarray import array,zeros,Float
318	axes = ['Beam','IF','Pol','Time']
319
320	if allaxes:
321	n = self.nbeam()self.nif()self.npol()*self.nrow()
322	shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
323	arr = array(zeros(n),shape=shp,type=Float)
324
325	for i in range(self.nbeam()):
326	self.setbeam(i)
327	for j in range(self.nif()):
328	self.setif(j)
329	for k in range(self.npol()):
330	self.setpol(k)
331	arr[i,j,k,:] = self._gettsys()
332	retval = {'axes': axes, 'data': arr, 'cursor':None}
333	tm = [self._gettime(val) for val in range(self.nrow())]
334	if self._vb:
335	self._print_values(retval,'Tsys',tm)
336	return retval
337
338	else:
339	i,j,k = (self.getbeam(),self.getif(),self.getpol())
340	statval = self._gettsys()
341	out = ''
342	for l in range(self.nrow()):
343	tm = self._gettime(l)
344	out += 'Time[%s]:\n' % (tm)
345	if self.nbeam() > 1: out += ' Beam[%d] ' % (i)
346	if self.nif() > 1: out += ' IF[%d] ' % (j)
347	if self.npol() > 1: out += ' Pol[%d] ' % (k)
348	out += '= %3.3f\n' % (statval[l])
349	out += "--------------------------------------------------\n"
350
351	if self._vb:
352	print "--------------------------------------------------"
353	print " TSys"
354	print "--------------------------------------------------"
355	print out
356	retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
357	return retval
358
359	def get_time(self, row=-1):
360	"""
361	Get a list of time stamps for the observations.
362	Return a string for each integration in the scantable.
363	Parameters:
364	row: row no of integration. Default -1 return all rows
365	Example:
366	none
367	"""
368	out = []
369	if row == -1:
370	for i in range(self.nrow()):
371	out.append(self._gettime(i))
372	return out
373	else:
374	if row < self.nrow():
375	return self._gettime(row)
376
377	def set_unit(self, unit='channel'):
378	"""
379	Set the unit for all following operations on this scantable
380	Parameters:
381	unit: optional unit, default is 'channel'
382	one of '*Hz','km/s','channel', ''
383	"""
384
385	if unit in ['','pixel', 'channel']:
386	unit = ''
387	inf = list(self._getcoordinfo())
388	inf[0] = unit
389	self._setcoordinfo(inf)
390	if self._p: self.plot()
391
392	def set_instrument (self, instr):
393	"""
394	Set the instrument for subsequent processing
395	Parameters:
396	instr: Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
397	'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
398	"""
399	self._setInstrument(instr)
400
401	def set_doppler(self, doppler='RADIO'):
402	"""
403	Set the doppler for all following operations on this scantable.
404	Parameters:
405	doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
406	"""
407
408	inf = list(self._getcoordinfo())
409	inf[2] = doppler
410	self._setcoordinfo(inf)
411	if self._p: self.plot()
412
413	def set_freqframe(self, frame=None):
414	"""
415	Set the frame type of the Spectral Axis.
416	Parameters:
417	frame: an optional frame type, default 'LSRK'.
418	Examples:
419	scan.set_freqframe('BARY')
420	"""
421	if not frame: frame = rcParams['scantable.freqframe']
422	valid = ['REST','TOPO','LSRD','LSRK','BARY', \
423	'GEO','GALACTO','LGROUP','CMB']
424	if frame in valid:
425	inf = list(self._getcoordinfo())
426	inf[1] = frame
427	self._setcoordinfo(inf)
428	else:
429	print "Please specify a valid freq type. Valid types are:\n",valid
430
431	def get_unit(self):
432	"""
433	Get the default unit set in this scantable
434	Parameters:
435	Returns:
436	A unit string
437	"""
438	inf = self._getcoordinfo()
439	unit = inf[0]
440	if unit == '': unit = 'channel'
441	return unit
442
443	def get_abcissa(self, rowno=0):
444	"""
445	Get the abcissa in the current coordinate setup for the currently
446	selected Beam/IF/Pol
447	Parameters:
448	rowno: an optional row number in the scantable. Default is the
449	first row, i.e. rowno=0
450	Returns:
451	The abcissa values and it's format string (as a dictionary)
452	"""
453	abc = self._getabcissa(rowno)
454	lbl = self._getabcissalabel(rowno)
455	return abc, lbl
456	#return {'abcissa':abc,'label':lbl}
457
458	def create_mask(self, args, *kwargs):
459	"""
460	Compute and return a mask based on [min,max] windows.
461	The specified windows are to be INCLUDED, when the mask is
462	applied.
463	Parameters:
464	[min,max],[min2,max2],...
465	Pairs of start/end points specifying the regions
466	to be masked
467	invert: optional argument. If specified as True,
468	return an inverted mask, i.e. the regions
469	specified are EXCLUDED
470	Example:
471	scan.set_unit('channel')
472
473	a)
474	msk = scan.set_mask([400,500],[800,900])
475	# masks everything outside 400 and 500
476	# and 800 and 900 in the unit 'channel'
477
478	b)
479	msk = scan.set_mask([400,500],[800,900], invert=True)
480	# masks the regions between 400 and 500
481	# and 800 and 900 in the unit 'channel'
482
483	"""
484	u = self._getcoordinfo()[0]
485	if self._vb:
486	if u == "": u = "channel"
487	print "The current mask window unit is", u
488	n = self.nchan()
489	data = self._getabcissa()
490	msk = zeros(n)
491	for window in args:
492	if (len(window) != 2 or window[0] > window[1] ):
493	print "A window needs to be defined as [min,max]"
494	return
495	for i in range(n):
496	if data[i] >= window[0] and data[i] < window[1]:
497	msk[i] = 1
498	if kwargs.has_key('invert'):
499	if kwargs.get('invert'):
500	from numarray import logical_not
501	msk = logical_not(msk)
502	return msk
503
504	def get_restfreqs(self):
505	"""
506	Get the restfrequency(s) stored in this scantable.
507	The return value(s) are always of unit 'Hz'
508	Parameters:
509	none
510	Returns:
511	a list of doubles
512	"""
513	return list(self._getrestfreqs())
514
515	def lines(self):
516	"""
517	Print the list of known spectral lines
518	"""
519	sdtable._lines(self)
520
521	def set_restfreqs(self, freqs=None, unit='Hz', lines=None, source=None, theif=None):
522	"""
523	Select the restfrequency for the specified source and IF OR
524	replace for all IFs. If the 'freqs' argument holds a scalar,
525	then that rest frequency will be applied to the selected
526	data (and added to the list of available rest frequencies).
527	In this way, you can set a rest frequency for each
528	source and IF combination. If the 'freqs' argument holds
529	a vector, then it MUST be of length the number of IFs
530	(and the available restfrequencies will be replaced by
531	this vector). In this case, all data ('source' and
532	'theif' are ignored) have the restfrequency set per IF according
533	to the corresponding value you give in the 'freqs' vector.
534	E.g. 'freqs=[1e9,2e9]' would mean IF 0 gets restfreq 1e9 and
535	IF 1 gets restfreq 2e9.
536
537	You can also specify the frequencies via known line names
538	in the argument 'lines'. Use 'freqs' or 'lines'. 'freqs'
539	takes precedence. See the list of known names via function
540	scantable.lines()
541	Parameters:
542	freqs: list of rest frequencies
543	unit: unit for rest frequency (default 'Hz')
544	lines: list of known spectral lines names (alternative to freqs).
545	See possible list via scantable.lines()
546	source: Source name (blank means all)
547	theif: IF (-1 means all)
548	Example:
549	scan.set_restfreqs(freqs=1.4e9, source='NGC253', theif=2)
550	scan.set_restfreqs(freqs=[1.4e9,1.67e9])
551	"""
552	if source is None:
553	source = ""
554	if theif is None:
555	theif = -1
556	t = type(freqs)
557	if t is int or t is float:
558	freqs = [freqs]
559	if freqs is None:
560	freqs = []
561	t = type(lines)
562	if t is str:
563	lines = [lines]
564	if lines is None:
565	lines = []
566	sdtable._setrestfreqs(self, freqs, unit, lines, source, theif)
567	return
568
569
570	def flag_spectrum(self, thebeam, theif, thepol):
571	"""
572	This flags a selected spectrum in the scan 'for good'.
573	USE WITH CARE - not reversible.
574	Use masks for non-permanent exclusion of channels.
575	Parameters:
576	thebeam,theif,thepol: all have to be explicitly
577	specified
578	Example:
579	scan.flag_spectrum(0,0,1)
580	flags the spectrum for Beam=0, IF=0, Pol=1
581	"""
582	if (thebeam < self.nbeam() and
583	theif < self.nif() and
584	thepol < self.npol()):
585	sdtable.setbeam(self, thebeam)
586	sdtable.setif(self, theif)
587	sdtable.setpol(self, thepol)
588	sdtable._flag(self)
589	else:
590	print "Please specify a valid (Beam/IF/Pol)"
591	return
592
593	def rotate_xyphase (self, angle, allaxes=None):
594	"""
595	Rotate the phase of the XY correlation. This is done in situ
596	in the data. So if you call this function more than once
597	then each call rotates the phase further.
598	Parameters:
599	angle: The angle (degrees) to rotate (add) by.
600	allaxes: If True apply to all spectra. Otherwise
601	apply only to the selected (beam/pol/if)spectra only.
602	The default is taken from .asaprc (True if none)
603	Examples:
604	scan.rotate_xyphase(2.3)
605	"""
606	if allaxes is None: allaxes = rcParams['scantable.allaxes']
607	sdtable._rotate_xyphase(self, angle, allaxes)
608
609	def plot(self, what='spectrum',col='Pol', panel=None):
610	"""
611	Plot the spectra contained in the scan. Alternatively you can also
612	Plot Tsys vs Time
613	Parameters:
614	what: a choice of 'spectrum' (default) or 'tsys'
615	col: which out of Beams/IFs/Pols should be colour stacked
616	panel: set up multiple panels, currently not working.
617	"""
618	print "Warning! Not fully functional. Use plotter.plot() instead"
619
620	validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
621
622	validyax = ['spectrum','tsys']
623	from asap.asaplot import ASAPlot
624	if not self._p:
625	self._p = ASAPlot()
626	#print "Plotting not enabled"
627	#return
628	if self._p.is_dead:
629	del self._p
630	self._p = ASAPlot()
631	npan = 1
632	x = None
633	if what == 'tsys':
634	n = self.nrow()
635	if n < 2:
636	print "Only one integration. Can't plot."
637	return
638	self._p.hold()
639	self._p.clear()
640	if panel == 'Time':
641	npan = self.nrow()
642	self._p.set_panels(rows=npan)
643	xlab,ylab,tlab = None,None,None
644	self._vb = False
645	sel = self.get_cursor()
646	for i in range(npan):
647	if npan > 1:
648	self._p.subplot(i)
649	for j in range(validcol[col]):
650	x = None
651	y = None
652	m = None
653	tlab = self._getsourcename(i)
654	import re
655	tlab = re.sub('_S','',tlab)
656	if col == 'Beam':
657	self.setbeam(j)
658	elif col == 'IF':
659	self.setif(j)
660	elif col == 'Pol':
661	self.setpol(j)
662	if what == 'tsys':
663	x = range(self.nrow())
664	xlab = 'Time [pixel]'
665	m = list(ones(len(x)))
666	y = []
667	ylab = r'$T_{sys}$'
668	for k in range(len(x)):
669	y.append(self._gettsys(k))
670	else:
671	x,xlab = self.get_abcissa(i)
672	y = self._getspectrum(i)
673	ylab = r'Flux'
674	m = self._getmask(i)
675	llab = col+' '+str(j)
676	self._p.set_line(label=llab)
677	self._p.plot(x,y,m)
678	self._p.set_axes('xlabel',xlab)
679	self._p.set_axes('ylabel',ylab)
680	self._p.set_axes('title',tlab)
681	self._p.release()
682	self.set_cursor(sel[0],sel[1],sel[2])
683	self._vb = rcParams['verbose']
684	return
685
686	print out
687
688	def _print_values(self, dat, label='', timestamps=[]):
689	d = dat['data']
690	a = dat['axes']
691	shp = d.getshape()
692	out = ''
693	for i in range(shp[3]):
694	out += '%s [%s]:\n' % (a[3],timestamps[i])
695	t = d[:,:,:,i]
696	for j in range(shp[0]):
697	if shp[0] > 1: out += ' %s[%d] ' % (a[0],j)
698	for k in range(shp[1]):
699	if shp[1] > 1: out += ' %s[%d] ' % (a[1],k)
700	for l in range(shp[2]):
701	if shp[2] > 1: out += ' %s[%d] ' % (a[2],l)
702	out += '= %3.3f\n' % (t[j,k,l])
703	out += "--------------------------------------------------\n"
704	print "--------------------------------------------------"
705	print " ", label
706	print "--------------------------------------------------"
707	print out

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

Download in other formats: