Context Navigation

source: trunk/python/scantable.py @ 445

Visit:

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

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

Download in other formats: