Context Navigation

source: trunk/python/scantable.py @ 436

Visit:

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

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

Download in other formats: