Context Navigation

source: trunk/python/scantable.py @ 402

Visit:

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

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

Download in other formats: