Context Navigation

scantable.py@ 1751

Visit:

Last change on this file since 1751 was 669, checked in by mar637, 19 years ago
fixed typo in import rotate_xyphase
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 57.2 KB

Rev	Line
[102]	1	from asap._asap import sdtable
[226]	2	from asap import rcParams
[189]	3	from numarray import ones,zeros
[102]	4	import sys
	5
	6	class scantable(sdtable):
	7	"""
	8	The ASAP container for scans
	9	"""
[113]	10
[484]	11	def __init__(self, filename, average=None, unit=None):
[102]	12	"""
	13	Create a scantable from a saved one or make a reference
	14	Parameters:
[181]	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
[102]	22	scantable
[484]	23	average: average all integrations withinb a scan on read.
	24	The default (True) is taken from .asaprc.
	25	unit: brightness unit; must be consistent with K or Jy.
[340]	26	Over-rides the default selected by the reader
	27	(input rpfits/sdfits/ms) or replaces the value
	28	in existing scantables
[489]	29	"""
	30	if average is None or type(average) is not bool:
	31	average = rcParams['scantable.autoaverage']
	32
[484]	33	varlist = vars()
[226]	34	self._vb = rcParams['verbose']
[113]	35	self._p = None
[484]	36
[181]	37	if isinstance(filename,sdtable):
	38	sdtable.__init__(self, filename)
[340]	39	if unit is not None:
[346]	40	self.set_fluxunit(unit)
[181]	41	else:
[411]	42	import os.path
	43	if not os.path.exists(filename):
	44	print "File '%s' not found." % (filename)
[226]	45	return
[411]	46	filename = os.path.expandvars(filename)
	47	if os.path.isdir(filename):
[181]	48	# crude check if asap table
[411]	49	if os.path.exists(filename+'/table.info'):
[181]	50	sdtable.__init__(self, filename)
[340]	51	if unit is not None:
	52	self.set_fluxunit(unit)
[181]	53	else:
[411]	54	print "The given file '%s'is not a valid asap table." % (filename)
[226]	55	return
[484]	56	else:
[181]	57	from asap._asap import sdreader
[340]	58	ifSel = -1
	59	beamSel = -1
[345]	60	r = sdreader(filename,ifSel,beamSel)
[181]	61	print 'Importing data...'
[411]	62	r._read([-1])
	63	tbl = r._getdata()
[346]	64	if unit is not None:
	65	tbl.set_fluxunit(unit)
[489]	66	if average:
	67	from asap._asap import average as _av
[226]	68	print 'Auto averaging integrations...'
[513]	69	tbl2 = _av((tbl,),(),True,'none')
[226]	70	sdtable.__init__(self,tbl2)
[513]	71	del tbl2
[226]	72	else:
	73	sdtable.__init__(self,tbl)
[513]	74	del r,tbl
[484]	75	self._add_history("scantable", varlist)
[102]	76
[451]	77	def save(self, name=None, format=None, stokes=False, overwrite=False):
[116]	78	"""
[280]	79	Store the scantable on disk. This can be an asap (aips++) Table, SDFITS,
	80	Image FITS or MS2 format.
[116]	81	Parameters:
[196]	82	name: the name of the outputfile. For format="FITS" this
[489]	83	is the directory file name into which all the files
[497]	84	will be written (default is 'asap_FITS'). For format
	85	"ASCII" this is the root file name (data in 'name'.txt
	86	and header in 'name'_header.txt)
[116]	87	format: an optional file format. Default is ASAP.
[280]	88	Allowed are - 'ASAP' (save as ASAP [aips++] Table),
[194]	89	'SDFITS' (save as SDFITS file)
	90	'FITS' (saves each row as a FITS Image)
[200]	91	'ASCII' (saves as ascii text file)
[226]	92	'MS2' (saves as an aips++
	93	MeasurementSet V2)
[489]	94	stokes: Convert to Stokes parameters (only available
	95	currently with FITS and ASCII formats.
	96	Default is False.
[411]	97	overwrite: If the file should be overwritten if it exists.
[256]	98	The default False is to return with warning
[411]	99	without writing the output. USE WITH CARE.
[116]	100	Example:
	101	scan.save('myscan.asap')
	102	scan.save('myscan.sdfits','SDFITS')
	103	"""
[411]	104	from os import path
[226]	105	if format is None: format = rcParams['scantable.save']
[256]	106	suffix = '.'+format.lower()
	107	if name is None or name =="":
	108	name = 'scantable'+suffix
[283]	109	print "No filename given. Using default name %s..." % name
[411]	110	name = path.expandvars(name)
[256]	111	if path.isfile(name) or path.isdir(name):
	112	if not overwrite:
	113	print "File %s already exists." % name
	114	return
[451]	115	format2 = format.upper()
	116	if format2 == 'ASAP':
[116]	117	self._save(name)
	118	else:
	119	from asap._asap import sdwriter as _sw
[451]	120	w = _sw(format2)
[445]	121	w.write(self, name, stokes)
[116]	122	return
	123
[102]	124	def copy(self):
	125	"""
	126	Return a copy of this scantable.
	127	Parameters:
[113]	128	none
[102]	129	Example:
	130	copiedscan = scan.copy()
	131	"""
[113]	132	sd = scantable(sdtable._copy(self))
	133	return sd
	134
[102]	135	def get_scan(self, scanid=None):
	136	"""
	137	Return a specific scan (by scanno) or collection of scans (by
	138	source name) in a new scantable.
	139	Parameters:
[513]	140	scanid: a (list of) scanno or a source name, unix-style
	141	patterns are accepted for source name matching, e.g.
	142	'*_R' gets all 'ref scans
[102]	143	Example:
[513]	144	# get all scans containing the source '323p459'
	145	newscan = scan.get_scan('323p459')
	146	# get all 'off' scans
	147	refscans = scan.get_scan('*_R')
	148	# get a susbset of scans by scanno (as listed in scan.summary())
	149	newscan = scan.get_scan([0,2,7,10])
[102]	150	"""
	151	if scanid is None:
	152	print "Please specify a scan no or name to retrieve from the scantable"
	153	try:
	154	if type(scanid) is str:
[113]	155	s = sdtable._getsource(self,scanid)
[102]	156	return scantable(s)
	157	elif type(scanid) is int:
[381]	158	s = sdtable._getscan(self,[scanid])
	159	return scantable(s)
	160	elif type(scanid) is list:
[113]	161	s = sdtable._getscan(self,scanid)
[102]	162	return scantable(s)
[381]	163	else:
	164	print "Illegal scanid type, use 'int' or 'list' if ints."
[102]	165	except RuntimeError:
	166	print "Couldn't find any match."
	167
	168	def __str__(self):
[381]	169	return sdtable._summary(self,True)
[102]	170
[381]	171	def summary(self,filename=None, verbose=None):
[102]	172	"""
	173	Print a summary of the contents of this scantable.
	174	Parameters:
	175	filename: the name of a file to write the putput to
	176	Default - no file output
[381]	177	verbose: print extra info such as the frequency table
	178	The default (False) is taken from .asaprc
[102]	179	"""
[381]	180	info = sdtable._summary(self, verbose)
	181	if verbose is None: verbose = rcParams['scantable.verbosesummary']
[102]	182	if filename is not None:
[256]	183	if filename is "":
	184	filename = 'scantable_summary.txt'
[415]	185	from os.path import expandvars, isdir
[411]	186	filename = expandvars(filename)
[415]	187	if not isdir(filename):
[413]	188	data = open(filename, 'w')
	189	data.write(info)
	190	data.close()
	191	else:
	192	print "Illegal file name '%s'." % (filename)
[102]	193	print info
[484]	194
[553]	195	def set_cursor(self, beam=0, IF=0, pol=0):
[102]	196	"""
	197	Set the spectrum for individual operations.
	198	Parameters:
[553]	199	beam, IF, pol: a number
[102]	200	Example:
[256]	201	scan.set_cursor(0,0,1)
[135]	202	pol1sig = scan.stats(all=False) # returns std dev for beam=0
[181]	203	# if=0, pol=1
[102]	204	"""
[484]	205	varlist = vars()
[553]	206	self.setbeam(beam)
	207	self.setpol(pol)
	208	self.setif(IF)
[484]	209	self._add_history("set_cursor",varlist)
[102]	210	return
	211
[256]	212	def get_cursor(self):
[113]	213	"""
	214	Return/print a the current 'cursor' into the Beam/IF/Pol cube.
	215	Parameters:
	216	none
	217	Returns:
	218	a list of values (currentBeam,currentIF,currentPol)
	219	Example:
	220	none
	221	"""
[102]	222	i = self.getbeam()
	223	j = self.getif()
	224	k = self.getpol()
[113]	225	if self._vb:
[181]	226	print "--------------------------------------------------"
[256]	227	print " Cursor position"
[181]	228	print "--------------------------------------------------"
[226]	229	out = 'Beam=%d IF=%d Pol=%d ' % (i,j,k)
[113]	230	print out
[102]	231	return i,j,k
	232
[436]	233	def stats(self, stat='stddev', mask=None, allaxes=None):
[102]	234	"""
[135]	235	Determine the specified statistic of the current beam/if/pol
[102]	236	Takes a 'mask' as an optional parameter to specify which
	237	channels should be excluded.
	238	Parameters:
[135]	239	stat: 'min', 'max', 'sumsq', 'sum', 'mean'
	240	'var', 'stddev', 'avdev', 'rms', 'median'
	241	mask: an optional mask specifying where the statistic
[102]	242	should be determined.
[436]	243	allaxes: if True apply to all spectra. Otherwise
	244	apply only to the selected (beam/pol/if)spectra only.
	245	The default is taken from .asaprc (True if none)
[102]	246	Example:
[113]	247	scan.set_unit('channel')
[102]	248	msk = scan.create_mask([100,200],[500,600])
[135]	249	scan.stats(stat='mean', mask=m)
[102]	250	"""
[436]	251	if allaxes is None: allaxes = rcParams['scantable.allaxes']
[135]	252	from asap._asap import stats as _stats
[256]	253	from numarray import array,zeros,Float
[102]	254	if mask == None:
	255	mask = ones(self.nchan())
[256]	256	axes = ['Beam','IF','Pol','Time']
	257
[436]	258	beamSel,IFSel,polSel = (self.getbeam(),self.getif(),self.getpol())
	259	if allaxes:
[256]	260	n = self.nbeam()self.nif()self.npol()*self.nrow()
	261	shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
	262	arr = array(zeros(n),shape=shp,type=Float)
	263
	264	for i in range(self.nbeam()):
	265	self.setbeam(i)
	266	for j in range(self.nif()):
	267	self.setif(j)
	268	for k in range(self.npol()):
	269	self.setpol(k)
	270	arr[i,j,k,:] = _stats(self,mask,stat,-1)
	271	retval = {'axes': axes, 'data': arr, 'cursor':None}
	272	tm = [self._gettime(val) for val in range(self.nrow())]
[181]	273	if self._vb:
[256]	274	self._print_values(retval,stat,tm)
[436]	275	self.setbeam(beamSel)
	276	self.setif(IFSel)
	277	self.setpol(polSel)
[256]	278	return retval
[102]	279
	280	else:
[241]	281	statval = _stats(self,mask,stat,-1)
[181]	282	out = ''
	283	for l in range(self.nrow()):
	284	tm = self._gettime(l)
	285	out += 'Time[%s]:\n' % (tm)
[436]	286	if self.nbeam() > 1: out += ' Beam[%d] ' % (beamSel)
	287	if self.nif() > 1: out += ' IF[%d] ' % (IFSel)
	288	if self.npol() > 1: out += ' Pol[%d] ' % (polSel)
[181]	289	out += '= %3.3f\n' % (statval[l])
	290	out += "--------------------------------------------------\n"
[226]	291
[181]	292	if self._vb:
	293	print "--------------------------------------------------"
	294	print " ",stat
	295	print "--------------------------------------------------"
	296	print out
[256]	297	retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
	298	return retval
[102]	299
[436]	300	def stddev(self,mask=None, allaxes=None):
[135]	301	"""
	302	Determine the standard deviation of the current beam/if/pol
	303	Takes a 'mask' as an optional parameter to specify which
	304	channels should be excluded.
	305	Parameters:
	306	mask: an optional mask specifying where the standard
	307	deviation should be determined.
[436]	308	allaxes: optional flag to show all or a cursor selected
[226]	309	spectrum of Beam/IF/Pol. Default is all or taken
	310	from .asaprc
[135]	311
	312	Example:
	313	scan.set_unit('channel')
	314	msk = scan.create_mask([100,200],[500,600])
	315	scan.stddev(mask=m)
	316	"""
[436]	317	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	318	return self.stats(stat='stddev',mask=mask, allaxes=allaxes);
[135]	319
[436]	320	def get_tsys(self, allaxes=None):
[113]	321	"""
	322	Return the System temperatures.
	323	Parameters:
[436]	324	allaxes: if True apply to all spectra. Otherwise
	325	apply only to the selected (beam/pol/if)spectra only.
	326	The default is taken from .asaprc (True if none)
[113]	327	Returns:
	328	a list of Tsys values.
	329	"""
[436]	330	if allaxes is None: allaxes = rcParams['scantable.allaxes']
[256]	331	from numarray import array,zeros,Float
	332	axes = ['Beam','IF','Pol','Time']
	333
[436]	334	if allaxes:
[256]	335	n = self.nbeam()self.nif()self.npol()*self.nrow()
	336	shp = [self.nbeam(),self.nif(),self.npol(),self.nrow()]
	337	arr = array(zeros(n),shape=shp,type=Float)
	338
	339	for i in range(self.nbeam()):
	340	self.setbeam(i)
	341	for j in range(self.nif()):
	342	self.setif(j)
	343	for k in range(self.npol()):
	344	self.setpol(k)
	345	arr[i,j,k,:] = self._gettsys()
	346	retval = {'axes': axes, 'data': arr, 'cursor':None}
	347	tm = [self._gettime(val) for val in range(self.nrow())]
[181]	348	if self._vb:
[256]	349	self._print_values(retval,'Tsys',tm)
	350	return retval
	351
[102]	352	else:
[256]	353	i,j,k = (self.getbeam(),self.getif(),self.getpol())
	354	statval = self._gettsys()
[181]	355	out = ''
	356	for l in range(self.nrow()):
	357	tm = self._gettime(l)
	358	out += 'Time[%s]:\n' % (tm)
	359	if self.nbeam() > 1: out += ' Beam[%d] ' % (i)
	360	if self.nif() > 1: out += ' IF[%d] ' % (j)
	361	if self.npol() > 1: out += ' Pol[%d] ' % (k)
[256]	362	out += '= %3.3f\n' % (statval[l])
	363	out += "--------------------------------------------------\n"
	364
[181]	365	if self._vb:
	366	print "--------------------------------------------------"
[256]	367	print " TSys"
[181]	368	print "--------------------------------------------------"
	369	print out
[256]	370	retval = {'axes': axes, 'data': array(statval), 'cursor':(i,j,k)}
	371	return retval
[113]	372
[407]	373	def get_time(self, row=-1):
[113]	374	"""
	375	Get a list of time stamps for the observations.
[181]	376	Return a string for each integration in the scantable.
[113]	377	Parameters:
[407]	378	row: row no of integration. Default -1 return all rows
[113]	379	Example:
	380	none
	381	"""
	382	out = []
[407]	383	if row == -1:
	384	for i in range(self.nrow()):
	385	out.append(self._gettime(i))
	386	return out
	387	else:
	388	if row < self.nrow():
	389	return self._gettime(row)
[102]	390
	391	def set_unit(self, unit='channel'):
	392	"""
	393	Set the unit for all following operations on this scantable
	394	Parameters:
	395	unit: optional unit, default is 'channel'
[113]	396	one of '*Hz','km/s','channel', ''
[102]	397	"""
[484]	398	varlist = vars()
[113]	399	if unit in ['','pixel', 'channel']:
	400	unit = ''
	401	inf = list(self._getcoordinfo())
	402	inf[0] = unit
	403	self._setcoordinfo(inf)
	404	if self._p: self.plot()
[484]	405	self._add_history("set_unit",varlist)
[113]	406
[484]	407	def set_instrument(self, instr):
[358]	408	"""
	409	Set the instrument for subsequent processing
	410	Parameters:
[407]	411	instr: Select from 'ATPKSMB', 'ATPKSHOH', 'ATMOPRA',
	412	'DSS-43' (Tid), 'CEDUNA', and 'HOBART'
[358]	413	"""
	414	self._setInstrument(instr)
[484]	415	self._add_history("set_instument",vars())
[358]	416
[276]	417	def set_doppler(self, doppler='RADIO'):
	418	"""
	419	Set the doppler for all following operations on this scantable.
	420	Parameters:
	421	doppler: One of 'RADIO', 'OPTICAL', 'Z', 'BETA', 'GAMMA'
	422	"""
[484]	423	varlist = vars()
[276]	424	inf = list(self._getcoordinfo())
	425	inf[2] = doppler
	426	self._setcoordinfo(inf)
	427	if self._p: self.plot()
[484]	428	self._add_history("set_doppler",vars())
[432]	429
[226]	430	def set_freqframe(self, frame=None):
[113]	431	"""
	432	Set the frame type of the Spectral Axis.
	433	Parameters:
[591]	434	frame: an optional frame type, default 'LSRK'. Valid frames are:
	435	'REST','TOPO','LSRD','LSRK','BARY',
	436	'GEO','GALACTO','LGROUP','CMB'
[113]	437	Examples:
	438	scan.set_freqframe('BARY')
	439	"""
[484]	440	if frame is None: frame = rcParams['scantable.freqframe']
	441	varlist = vars()
[113]	442	valid = ['REST','TOPO','LSRD','LSRK','BARY', \
	443	'GEO','GALACTO','LGROUP','CMB']
[591]	444
[113]	445	if frame in valid:
	446	inf = list(self._getcoordinfo())
	447	inf[1] = frame
	448	self._setcoordinfo(inf)
[484]	449	self._add_history("set_freqframe",varlist)
[102]	450	else:
[113]	451	print "Please specify a valid freq type. Valid types are:\n",valid
	452
	453	def get_unit(self):
	454	"""
	455	Get the default unit set in this scantable
	456	Parameters:
	457	Returns:
	458	A unit string
	459	"""
	460	inf = self._getcoordinfo()
	461	unit = inf[0]
	462	if unit == '': unit = 'channel'
	463	return unit
[102]	464
[158]	465	def get_abcissa(self, rowno=0):
[102]	466	"""
[158]	467	Get the abcissa in the current coordinate setup for the currently
[113]	468	selected Beam/IF/Pol
	469	Parameters:
[226]	470	rowno: an optional row number in the scantable. Default is the
	471	first row, i.e. rowno=0
[113]	472	Returns:
[407]	473	The abcissa values and it's format string (as a dictionary)
[113]	474	"""
[256]	475	abc = self._getabcissa(rowno)
[407]	476	lbl = self._getabcissalabel(rowno)
[158]	477	return abc, lbl
[407]	478	#return {'abcissa':abc,'label':lbl}
[113]	479
	480	def create_mask(self, args, *kwargs):
	481	"""
[102]	482	Compute and return a mask based on [min,max] windows.
[189]	483	The specified windows are to be INCLUDED, when the mask is
[113]	484	applied.
[102]	485	Parameters:
	486	[min,max],[min2,max2],...
	487	Pairs of start/end points specifying the regions
	488	to be masked
[189]	489	invert: optional argument. If specified as True,
	490	return an inverted mask, i.e. the regions
	491	specified are EXCLUDED
[513]	492	row: create the mask using the specified row for
	493	unit conversions, default is row=0
	494	only necessary if frequency varies over rows.
[102]	495	Example:
[113]	496	scan.set_unit('channel')
	497
	498	a)
[102]	499	msk = scan.set_mask([400,500],[800,900])
[189]	500	# masks everything outside 400 and 500
[113]	501	# and 800 and 900 in the unit 'channel'
	502
	503	b)
	504	msk = scan.set_mask([400,500],[800,900], invert=True)
[189]	505	# masks the regions between 400 and 500
[113]	506	# and 800 and 900 in the unit 'channel'
	507
[102]	508	"""
[513]	509	row = 0
	510	if kwargs.has_key("row"):
	511	row = kwargs.get("row")
	512	data = self._getabcissa(row)
[113]	513	u = self._getcoordinfo()[0]
	514	if self._vb:
	515	if u == "": u = "channel"
	516	print "The current mask window unit is", u
[102]	517	n = self.nchan()
[189]	518	msk = zeros(n)
[484]	519	for window in args:
[102]	520	if (len(window) != 2 or window[0] > window[1] ):
	521	print "A window needs to be defined as [min,max]"
	522	return
	523	for i in range(n):
[113]	524	if data[i] >= window[0] and data[i] < window[1]:
[189]	525	msk[i] = 1
[113]	526	if kwargs.has_key('invert'):
	527	if kwargs.get('invert'):
	528	from numarray import logical_not
[513]	529	msk = logical_not(msk)
[102]	530	return msk
[113]	531
[256]	532	def get_restfreqs(self):
	533	"""
	534	Get the restfrequency(s) stored in this scantable.
	535	The return value(s) are always of unit 'Hz'
	536	Parameters:
	537	none
	538	Returns:
	539	a list of doubles
	540	"""
	541	return list(self._getrestfreqs())
[102]	542
[402]	543	def lines(self):
[391]	544	"""
[402]	545	Print the list of known spectral lines
	546	"""
	547	sdtable._lines(self)
	548
[484]	549	def set_restfreqs(self, freqs=None, unit='Hz', lines=None, source=None,
	550	theif=None):
[402]	551	"""
[393]	552	Select the restfrequency for the specified source and IF OR
	553	replace for all IFs. If the 'freqs' argument holds a scalar,
	554	then that rest frequency will be applied to the selected
	555	data (and added to the list of available rest frequencies).
	556	In this way, you can set a rest frequency for each
	557	source and IF combination. If the 'freqs' argument holds
	558	a vector, then it MUST be of length the number of IFs
	559	(and the available restfrequencies will be replaced by
	560	this vector). In this case, all data ('source' and
	561	'theif' are ignored) have the restfrequency set per IF according
	562	to the corresponding value you give in the 'freqs' vector.
	563	E.g. 'freqs=[1e9,2e9]' would mean IF 0 gets restfreq 1e9 and
	564	IF 1 gets restfreq 2e9.
[402]	565
	566	You can also specify the frequencies via known line names
	567	in the argument 'lines'. Use 'freqs' or 'lines'. 'freqs'
	568	takes precedence. See the list of known names via function
	569	scantable.lines()
[391]	570	Parameters:
[402]	571	freqs: list of rest frequencies
[393]	572	unit: unit for rest frequency (default 'Hz')
[402]	573	lines: list of known spectral lines names (alternative to freqs).
	574	See possible list via scantable.lines()
[391]	575	source: Source name (blank means all)
	576	theif: IF (-1 means all)
	577	Example:
[393]	578	scan.set_restfreqs(freqs=1.4e9, source='NGC253', theif=2)
	579	scan.set_restfreqs(freqs=[1.4e9,1.67e9])
[391]	580	"""
[484]	581	varlist = vars()
[391]	582	if source is None:
	583	source = ""
	584	if theif is None:
	585	theif = -1
[393]	586	t = type(freqs)
	587	if t is int or t is float:
	588	freqs = [freqs]
[402]	589	if freqs is None:
	590	freqs = []
	591	t = type(lines)
	592	if t is str:
	593	lines = [lines]
	594	if lines is None:
	595	lines = []
	596	sdtable._setrestfreqs(self, freqs, unit, lines, source, theif)
[484]	597	self._add_history("set_restfreqs", varlist)
	598
[391]	599
	600
[102]	601	def flag_spectrum(self, thebeam, theif, thepol):
	602	"""
	603	This flags a selected spectrum in the scan 'for good'.
	604	USE WITH CARE - not reversible.
	605	Use masks for non-permanent exclusion of channels.
	606	Parameters:
	607	thebeam,theif,thepol: all have to be explicitly
	608	specified
	609	Example:
	610	scan.flag_spectrum(0,0,1)
	611	flags the spectrum for Beam=0, IF=0, Pol=1
	612	"""
[407]	613	if (thebeam < self.nbeam() and
	614	theif < self.nif() and
	615	thepol < self.npol()):
	616	sdtable.setbeam(self, thebeam)
	617	sdtable.setif(self, theif)
	618	sdtable.setpol(self, thepol)
	619	sdtable._flag(self)
[484]	620	self._add_history("flag_spectrum", vars())
[102]	621	else:
	622	print "Please specify a valid (Beam/IF/Pol)"
	623	return
[113]	624
	625	def plot(self, what='spectrum',col='Pol', panel=None):
	626	"""
	627	Plot the spectra contained in the scan. Alternatively you can also
	628	Plot Tsys vs Time
	629	Parameters:
	630	what: a choice of 'spectrum' (default) or 'tsys'
[135]	631	col: which out of Beams/IFs/Pols should be colour stacked
[113]	632	panel: set up multiple panels, currently not working.
	633	"""
[283]	634	print "Warning! Not fully functional. Use plotter.plot() instead"
	635
[113]	636	validcol = {'Beam':self.nbeam(),'IF':self.nif(),'Pol':self.npol()}
[124]	637
[113]	638	validyax = ['spectrum','tsys']
[189]	639	from asap.asaplot import ASAPlot
[113]	640	if not self._p:
	641	self._p = ASAPlot()
[189]	642	#print "Plotting not enabled"
	643	#return
	644	if self._p.is_dead:
	645	del self._p
	646	self._p = ASAPlot()
[113]	647	npan = 1
	648	x = None
	649	if what == 'tsys':
	650	n = self.nrow()
	651	if n < 2:
	652	print "Only one integration. Can't plot."
	653	return
[124]	654	self._p.hold()
	655	self._p.clear()
[113]	656	if panel == 'Time':
	657	npan = self.nrow()
[124]	658	self._p.set_panels(rows=npan)
[113]	659	xlab,ylab,tlab = None,None,None
[226]	660	self._vb = False
[256]	661	sel = self.get_cursor()
[113]	662	for i in range(npan):
[226]	663	if npan > 1:
	664	self._p.subplot(i)
[113]	665	for j in range(validcol[col]):
	666	x = None
	667	y = None
	668	m = None
	669	tlab = self._getsourcename(i)
[124]	670	import re
	671	tlab = re.sub('_S','',tlab)
[113]	672	if col == 'Beam':
	673	self.setbeam(j)
	674	elif col == 'IF':
	675	self.setif(j)
	676	elif col == 'Pol':
	677	self.setpol(j)
	678	if what == 'tsys':
	679	x = range(self.nrow())
	680	xlab = 'Time [pixel]'
	681	m = list(ones(len(x)))
	682	y = []
	683	ylab = r'$T_{sys}$'
	684	for k in range(len(x)):
	685	y.append(self._gettsys(k))
	686	else:
[158]	687	x,xlab = self.get_abcissa(i)
[256]	688	y = self._getspectrum(i)
[113]	689	ylab = r'Flux'
[256]	690	m = self._getmask(i)
[113]	691	llab = col+' '+str(j)
	692	self._p.set_line(label=llab)
	693	self._p.plot(x,y,m)
[124]	694	self._p.set_axes('xlabel',xlab)
	695	self._p.set_axes('ylabel',ylab)
	696	self._p.set_axes('title',tlab)
[113]	697	self._p.release()
[256]	698	self.set_cursor(sel[0],sel[1],sel[2])
[226]	699	self._vb = rcParams['verbose']
[113]	700	return
[256]	701
	702	print out
	703
	704	def _print_values(self, dat, label='', timestamps=[]):
	705	d = dat['data']
	706	a = dat['axes']
	707	shp = d.getshape()
	708	out = ''
	709	for i in range(shp[3]):
	710	out += '%s [%s]:\n' % (a[3],timestamps[i])
	711	t = d[:,:,:,i]
	712	for j in range(shp[0]):
	713	if shp[0] > 1: out += ' %s[%d] ' % (a[0],j)
	714	for k in range(shp[1]):
	715	if shp[1] > 1: out += ' %s[%d] ' % (a[1],k)
	716	for l in range(shp[2]):
	717	if shp[2] > 1: out += ' %s[%d] ' % (a[2],l)
	718	out += '= %3.3f\n' % (t[j,k,l])
[484]	719	out += "-"*80
	720	out += "\n"
	721	print "-"*80
[256]	722	print " ", label
[484]	723	print "-"*80
[256]	724	print out
[484]	725
	726	def history(self):
	727	hist = list(self._gethistory())
	728	print "-"*80
	729	for h in hist:
[489]	730	if h.startswith("---"):
	731	print h
	732	else:
	733	items = h.split("##")
	734	date = items[0]
	735	func = items[1]
	736	items = items[2:]
	737	print date
	738	print "Function: %s\n Parameters:" % (func)
	739	for i in items:
	740	s = i.split("=")
	741	print " %s = %s" % (s[0],s[1])
	742	print "-"*80
[484]	743	return
	744
[513]	745	#
	746	# Maths business
	747	#
	748
[558]	749	def average_time(self, mask=None, scanav=False, weight='tint'):
[513]	750	"""
	751	Return the (time) average of a scan, or apply it 'insitu'.
	752	Note:
	753	in channels only
	754	The cursor of the output scan is set to 0.
	755	Parameters:
	756	one scan or comma separated scans
	757	mask: an optional mask (only used for 'var' and 'tsys'
	758	weighting)
[558]	759	scanav: True averages each scan separately
	760	False (default) averages all scans together,
[524]	761	weight: Weighting scheme. 'none', 'var' (1/var(spec)
[535]	762	weighted), 'tsys' (1/Tsys**2 weighted), 'tint'
	763	(integration time weighted) or 'tintsys' (Tint/Tsys**2).
	764	The default is 'tint'
[513]	765	Example:
	766	# time average the scantable without using a mask
	767	newscan = scan.average_time()
	768	"""
	769	varlist = vars()
[524]	770	if weight is None: weight = 'tint'
[513]	771	if mask is None: mask = ()
	772	from asap._asap import average as _av
	773	s = scantable(_av((self,), mask, scanav, weight))
	774	s._add_history("average_time",varlist)
	775	return s
	776
	777	def convert_flux(self, jyperk=None, eta=None, d=None, insitu=None,
	778	allaxes=None):
	779	"""
	780	Return a scan where all spectra are converted to either
	781	Jansky or Kelvin depending upon the flux units of the scan table.
	782	By default the function tries to look the values up internally.
	783	If it can't find them (or if you want to over-ride), you must
	784	specify EITHER jyperk OR eta (and D which it will try to look up
	785	also if you don't set it). jyperk takes precedence if you set both.
	786	Parameters:
	787	jyperk: the Jy / K conversion factor
	788	eta: the aperture efficiency
	789	d: the geomtric diameter (metres)
	790	insitu: if False a new scantable is returned.
	791	Otherwise, the scaling is done in-situ
	792	The default is taken from .asaprc (False)
	793	allaxes: if True apply to all spectra. Otherwise
	794	apply only to the selected (beam/pol/if)spectra only
	795	The default is taken from .asaprc (True if none)
	796	"""
	797	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	798	if insitu is None: insitu = rcParams['insitu']
	799	varlist = vars()
	800	if jyperk is None: jyperk = -1.0
	801	if d is None: d = -1.0
	802	if eta is None: eta = -1.0
	803	if not insitu:
	804	from asap._asap import convertflux as _convert
	805	s = scantable(_convert(self, d, eta, jyperk, allaxes))
	806	s._add_history("convert_flux", varlist)
	807	return s
	808	else:
	809	from asap._asap import convertflux_insitu as _convert
	810	_convert(self, d, eta, jyperk, allaxes)
	811	self._add_history("convert_flux", varlist)
	812	return
	813
	814	def gain_el(self, poly=None, filename="", method="linear",
	815	insitu=None, allaxes=None):
	816	"""
	817	Return a scan after applying a gain-elevation correction.
	818	The correction can be made via either a polynomial or a
	819	table-based interpolation (and extrapolation if necessary).
	820	You specify polynomial coefficients, an ascii table or neither.
	821	If you specify neither, then a polynomial correction will be made
	822	with built in coefficients known for certain telescopes (an error
	823	will occur if the instrument is not known).
	824	The data and Tsys are divided by the scaling factors.
	825	Parameters:
	826	poly: Polynomial coefficients (default None) to compute a
	827	gain-elevation correction as a function of
	828	elevation (in degrees).
	829	filename: The name of an ascii file holding correction factors.
	830	The first row of the ascii file must give the column
	831	names and these MUST include columns
	832	"ELEVATION" (degrees) and "FACTOR" (multiply data
	833	by this) somewhere.
	834	The second row must give the data type of the
	835	column. Use 'R' for Real and 'I' for Integer.
	836	An example file would be
	837	(actual factors are arbitrary) :
	838
	839	TIME ELEVATION FACTOR
	840	R R R
	841	0.1 0 0.8
	842	0.2 20 0.85
	843	0.3 40 0.9
	844	0.4 60 0.85
	845	0.5 80 0.8
	846	0.6 90 0.75
	847	method: Interpolation method when correcting from a table.
	848	Values are "nearest", "linear" (default), "cubic"
	849	and "spline"
	850	insitu: if False a new scantable is returned.
	851	Otherwise, the scaling is done in-situ
	852	The default is taken from .asaprc (False)
	853	allaxes: If True apply to all spectra. Otherwise
	854	apply only to the selected (beam/pol/if) spectra only
	855	The default is taken from .asaprc (True if none)
	856	"""
	857
	858	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	859	if insitu is None: insitu = rcParams['insitu']
	860	varlist = vars()
	861	if poly is None:
	862	poly = ()
	863	from os.path import expandvars
	864	filename = expandvars(filename)
	865	if not insitu:
	866	from asap._asap import gainel as _gainEl
	867	s = scantable(_gainEl(self, poly, filename, method, allaxes))
	868	s._add_history("gain_el", varlist)
	869	return s
	870	else:
	871	from asap._asap import gainel_insitu as _gainEl
	872	_gainEl(self, poly, filename, method, allaxes)
	873	self._add_history("gain_el", varlist)
	874	return
	875
	876	def freq_align(self, reftime=None, method='cubic', perif=False,
	877	insitu=None):
	878	"""
	879	Return a scan where all rows have been aligned in frequency/velocity.
	880	The alignment frequency frame (e.g. LSRK) is that set by function
	881	set_freqframe.
	882	Parameters:
	883	reftime: reference time to align at. By default, the time of
	884	the first row of data is used.
	885	method: Interpolation method for regridding the spectra.
	886	Choose from "nearest", "linear", "cubic" (default)
	887	and "spline"
	888	perif: Generate aligners per freqID (no doppler tracking) or
	889	per IF (scan-based doppler tracking)
	890	insitu: if False a new scantable is returned.
	891	Otherwise, the scaling is done in-situ
	892	The default is taken from .asaprc (False)
	893	"""
	894	if insitu is None: insitu = rcParams['insitu']
	895	varlist = vars()
	896	if reftime is None: reftime = ''
	897	perfreqid = not perif
	898	if not insitu:
	899	from asap._asap import freq_align as _align
	900	s = scantable(_align(self, reftime, method, perfreqid))
	901	s._add_history("freq_align", varlist)
	902	return s
	903	else:
	904	from asap._asap import freq_align_insitu as _align
	905	_align(self, reftime, method, perfreqid)
	906	self._add_history("freq_align", varlist)
	907	return
	908
	909	def opacity(self, tau, insitu=None, allaxes=None):
	910	"""
	911	Apply an opacity correction. The data
	912	and Tsys are multiplied by the correction factor.
	913	Parameters:
	914	tau: Opacity from which the correction factor is
	915	exp(tau*ZD)
	916	where ZD is the zenith-distance
	917	insitu: if False a new scantable is returned.
	918	Otherwise, the scaling is done in-situ
	919	The default is taken from .asaprc (False)
	920	allaxes: if True apply to all spectra. Otherwise
	921	apply only to the selected (beam/pol/if)spectra only
	922	The default is taken from .asaprc (True if none)
	923	"""
	924	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	925	if insitu is None: insitu = rcParams['insitu']
	926	varlist = vars()
	927	if not insitu:
	928	from asap._asap import opacity as _opacity
	929	s = scantable(_opacity(self, tau, allaxes))
	930	s._add_history("opacity", varlist)
	931	return s
	932	else:
	933	from asap._asap import opacity_insitu as _opacity
	934	_opacity(self, tau, allaxes)
	935	self._add_history("opacity", varlist)
	936	return
	937
	938	def bin(self, width=5, insitu=None):
	939	"""
	940	Return a scan where all spectra have been binned up.
	941	width: The bin width (default=5) in pixels
	942	insitu: if False a new scantable is returned.
	943	Otherwise, the scaling is done in-situ
	944	The default is taken from .asaprc (False)
	945	"""
	946	if insitu is None: insitu = rcParams['insitu']
	947	varlist = vars()
	948	if not insitu:
	949	from asap._asap import bin as _bin
	950	s = scantable(_bin(self, width))
	951	s._add_history("bin",varlist)
	952	return s
	953	else:
	954	from asap._asap import bin_insitu as _bin
	955	_bin(self, width)
	956	self._add_history("bin",varlist)
	957	return
	958
	959
	960	def resample(self, width=5, method='cubic', insitu=None):
	961	"""
	962	Return a scan where all spectra have been binned up
	963	width: The bin width (default=5) in pixels
	964	method: Interpolation method when correcting from a table.
	965	Values are "nearest", "linear", "cubic" (default)
	966	and "spline"
	967	insitu: if False a new scantable is returned.
	968	Otherwise, the scaling is done in-situ
	969	The default is taken from .asaprc (False)
	970	"""
	971	if insitu is None: insitu = rcParams['insitu']
	972	varlist = vars()
	973	if not insitu:
	974	from asap._asap import resample as _resample
	975	s = scantable(_resample(self, method, width))
	976	s._add_history("resample",varlist)
	977	return s
	978	else:
	979	from asap._asap import resample_insitu as _resample
	980	_resample(self, method, width)
	981	self._add_history("resample",varlist)
	982	return
	983
	984	def average_pol(self, mask=None, weight='none', insitu=None):
	985	"""
	986	Average the Polarisations together.
	987	The polarisation cursor of the output scan is set to 0
	988	Parameters:
	989	mask: An optional mask defining the region, where the
	990	averaging will be applied. The output will have all
	991	specified points masked.
[532]	992	weight: Weighting scheme. 'none' (default), 'var' (1/var(spec)
	993	weighted), or 'tsys' (1/Tsys**2 weighted)
[513]	994	insitu: if False a new scantable is returned.
	995	Otherwise, the scaling is done in-situ
	996	The default is taken from .asaprc (False)
	997	"""
	998	if insitu is None: insitu = rcParams['insitu']
	999	varlist = vars()
	1000
	1001	if mask is None:
	1002	mask = ()
	1003	if not insitu:
	1004	from asap._asap import averagepol as _avpol
	1005	s = scantable(_avpol(self, mask, weight))
	1006	s._add_history("average_pol",varlist)
	1007	return s
	1008	else:
	1009	from asap._asap import averagepol_insitu as _avpol
	1010	_avpol(self, mask, weight)
	1011	self._add_history("average_pol",varlist)
	1012	return
	1013
	1014	def smooth(self, kernel="hanning", width=5.0, insitu=None, allaxes=None):
	1015	"""
	1016	Smooth the spectrum by the specified kernel (conserving flux).
	1017	Parameters:
	1018	scan: The input scan
	1019	kernel: The type of smoothing kernel. Select from
	1020	'hanning' (default), 'gaussian' and 'boxcar'.
	1021	The first three characters are sufficient.
	1022	width: The width of the kernel in pixels. For hanning this is
	1023	ignored otherwise it defauls to 5 pixels.
	1024	For 'gaussian' it is the Full Width Half
	1025	Maximum. For 'boxcar' it is the full width.
	1026	insitu: if False a new scantable is returned.
	1027	Otherwise, the scaling is done in-situ
	1028	The default is taken from .asaprc (False)
	1029	allaxes: If True (default) apply to all spectra. Otherwise
	1030	apply only to the selected (beam/pol/if)spectra only
	1031	The default is taken from .asaprc (True if none)
	1032	Example:
	1033	none
	1034	"""
	1035	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	1036	if insitu is None: insitu = rcParams['insitu']
	1037	varlist = vars()
	1038	if not insitu:
	1039	from asap._asap import smooth as _smooth
	1040	s = scantable(_smooth(self,kernel,width,allaxes))
	1041	s._add_history("smooth", varlist)
	1042	return s
	1043	else:
	1044	from asap._asap import smooth_insitu as _smooth
	1045	_smooth(self,kernel,width,allaxes)
	1046	self._add_history("smooth", varlist)
	1047	return
	1048
	1049	def poly_baseline(self, mask=None, order=0, insitu=None):
	1050	"""
	1051	Return a scan which has been baselined (all rows) by a polynomial.
	1052	Parameters:
	1053	scan: a scantable
	1054	mask: an optional mask
	1055	order: the order of the polynomial (default is 0)
	1056	insitu: if False a new scantable is returned.
	1057	Otherwise, the scaling is done in-situ
	1058	The default is taken from .asaprc (False)
	1059	Example:
	1060	# return a scan baselined by a third order polynomial,
	1061	# not using a mask
	1062	bscan = scan.poly_baseline(order=3)
[579]	1063	"""
[513]	1064	if insitu is None: insitu = rcParams['insitu']
	1065	varlist = vars()
	1066	if mask is None:
	1067	from numarray import ones
[626]	1068	mask = list(ones(self.nchan()))
[513]	1069	from asap.asapfitter import fitter
	1070	f = fitter()
	1071	f._verbose(True)
	1072	f.set_scan(self, mask)
	1073	f.set_function(poly=order)
	1074	sf = f.auto_fit(insitu)
	1075	if insitu:
	1076	self._add_history("poly_baseline", varlist)
	1077	return
	1078	else:
	1079	sf._add_history("poly_baseline", varlist)
	1080	return sf
	1081
	1082	def auto_poly_baseline(self, mask=None, edge=(0,0), order=0,
	1083	threshold=3,insitu=None):
	1084	"""
	1085	Return a scan which has been baselined (all rows) by a polynomial.
	1086	Spectral lines are detected first using linefinder and masked out
	1087	to avoid them affecting the baseline solution.
	1088
	1089	Parameters:
	1090	scan: a scantable
	1091	mask: an optional mask retreived from scantable
	1092	edge: an optional number of channel to drop at
	1093	the edge of spectrum. If only one value is
	1094	specified, the same number will be dropped from
	1095	both sides of the spectrum. Default is to keep
	1096	all channels
	1097	order: the order of the polynomial (default is 0)
	1098	threshold: the threshold used by line finder. It is better to
	1099	keep it large as only strong lines affect the
	1100	baseline solution.
	1101	insitu: if False a new scantable is returned.
	1102	Otherwise, the scaling is done in-situ
	1103	The default is taken from .asaprc (False)
	1104
	1105	Example:
[516]	1106	scan2=scan.auto_poly_baseline(order=7)
[513]	1107	"""
	1108	if insitu is None: insitu = rcParams['insitu']
	1109	varlist = vars()
	1110	from asap.asapfitter import fitter
	1111	from asap.asaplinefind import linefinder
	1112	from asap import _is_sequence_or_number as _is_valid
	1113
	1114	if not _is_valid(edge, int):
[516]	1115	raise RuntimeError, "Parameter 'edge' has to be an integer or a \
	1116	pair of integers specified as a tuple"
[513]	1117
	1118	# setup fitter
	1119	f = fitter()
	1120	f._verbose(True)
	1121	f.set_function(poly=order)
	1122
	1123	# setup line finder
	1124	fl=linefinder()
	1125	fl.set_options(threshold=threshold)
	1126
	1127	if not insitu:
	1128	workscan=self.copy()
	1129	else:
	1130	workscan=self
	1131
	1132	vb=workscan._vb
	1133	# remember the verbose parameter and selection
	1134	workscan._vb=False
	1135	sel=workscan.get_cursor()
	1136	rows=range(workscan.nrow())
	1137	for i in range(workscan.nbeam()):
	1138	workscan.setbeam(i)
	1139	for j in range(workscan.nif()):
	1140	workscan.setif(j)
	1141	for k in range(workscan.npol()):
	1142	workscan.setpol(k)
	1143	if f._vb:
	1144	print "Processing:"
	1145	print 'Beam[%d], IF[%d], Pol[%d]' % (i,j,k)
	1146	for iRow in rows:
	1147	fl.set_scan(workscan,mask,edge)
	1148	fl.find_lines(iRow)
	1149	f.set_scan(workscan, fl.get_mask())
	1150	f.x=workscan._getabcissa(iRow)
	1151	f.y=workscan._getspectrum(iRow)
	1152	f.data=None
	1153	f.fit()
	1154	x=f.get_parameters()
	1155	workscan._setspectrum(f.fitter.getresidual(),iRow)
	1156	workscan.set_cursor(sel[0],sel[1],sel[2])
	1157	workscan._vb = vb
	1158	if not insitu:
	1159	return workscan
	1160
	1161	def rotate_linpolphase(self, angle, allaxes=None):
	1162	"""
	1163	Rotate the phase of the complex polarization O=Q+iU correlation.
	1164	This is always done in situ in the raw data. So if you call this
	1165	function more than once then each call rotates the phase further.
	1166	Parameters:
	1167	angle: The angle (degrees) to rotate (add) by.
	1168	allaxes: If True apply to all spectra. Otherwise
	1169	apply only to the selected (beam/pol/if)spectra only.
	1170	The default is taken from .asaprc (True if none)
	1171	Examples:
	1172	scan.rotate_linpolphase(2.3)
	1173	"""
	1174	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	1175	varlist = vars()
	1176	from asap._asap import _rotate_linpolphase as _rotate
	1177	_rotate(self, angle, allaxes)
	1178	self._add_history("rotate_linpolphase", varlist)
	1179	return
	1180
	1181
	1182	def rotate_xyphase(self, angle, allaxes=None):
	1183	"""
	1184	Rotate the phase of the XY correlation. This is always done in situ
	1185	in the data. So if you call this function more than once
	1186	then each call rotates the phase further.
	1187	Parameters:
	1188	angle: The angle (degrees) to rotate (add) by.
	1189	allaxes: If True apply to all spectra. Otherwise
	1190	apply only to the selected (beam/pol/if)spectra only.
	1191	The default is taken from .asaprc (True if none)
	1192	Examples:
	1193	scan.rotate_xyphase(2.3)
	1194	"""
	1195	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	1196	varlist = vars()
[669]	1197	from asap._asap import _rotate_xyphase
[513]	1198	_rotate_xyphase(self, angle, allaxes)
	1199	self._add_history("rotate_xyphase", varlist)
	1200	return
	1201
	1202
	1203	def add(self, offset, insitu=None, allaxes=None):
	1204	"""
	1205	Return a scan where all spectra have the offset added
	1206	Parameters:
	1207	offset: the offset
	1208	insitu: if False a new scantable is returned.
	1209	Otherwise, the scaling is done in-situ
	1210	The default is taken from .asaprc (False)
	1211	allaxes: if True apply to all spectra. Otherwise
	1212	apply only to the selected (beam/pol/if)spectra only
	1213	The default is taken from .asaprc (True if none)
	1214	"""
	1215	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	1216	if insitu is None: insitu = rcParams['insitu']
	1217	varlist = vars()
	1218	if not insitu:
	1219	from asap._asap import add as _add
	1220	s = scantable(_add(self, offset, allaxes))
	1221	s._add_history("add",varlist)
	1222	return s
	1223	else:
	1224	from asap._asap import add_insitu as _add
	1225	_add(self, offset, allaxes)
	1226	self._add_history("add",varlist)
	1227	return
	1228
	1229	def scale(self, factor, insitu=None, allaxes=None, tsys=True):
	1230	"""
	1231	Return a scan where all spectra are scaled by the give 'factor'
	1232	Parameters:
	1233	factor: the scaling factor
	1234	insitu: if False a new scantable is returned.
	1235	Otherwise, the scaling is done in-situ
	1236	The default is taken from .asaprc (False)
	1237	allaxes: if True apply to all spectra. Otherwise
	1238	apply only to the selected (beam/pol/if)spectra only.
	1239	The default is taken from .asaprc (True if none)
	1240	tsys: if True (default) then apply the operation to Tsys
	1241	as well as the data
	1242	"""
	1243	if allaxes is None: allaxes = rcParams['scantable.allaxes']
	1244	if insitu is None: insitu = rcParams['insitu']
	1245	varlist = vars()
	1246	if not insitu:
	1247	from asap._asap import scale as _scale
	1248	s = scantable(_scale(self, factor, allaxes, tsys))
	1249	s._add_history("scale",varlist)
	1250	return s
	1251	else:
	1252	from asap._asap import scale_insitu as _scale
	1253	_scale(self, factor, allaxes, tsys)
	1254	self._add_history("scale",varlist)
	1255	return
	1256
[666]	1257	def auto_quotient(self, mode='suffix', preserve=True):
	1258	"""
	1259	This function allows to build quotients automatically.
	1260	It assumes the observation to have the same numer of
	1261	"ons" and "offs"
	1262	It will support "closest off in time" in the future
	1263	Parameters:
	1264	mode: the on/off detection mode; 'suffix' (default)
	1265	'suffix' identifies 'off' scans by the
	1266	trailing '_R' (Mopra/Parkes) or
	1267	'_e'/'_w' (Tid)
	1268	preserve: you can preserve (default) the continuum or
	1269	remove it. The equations used are
	1270	preserve: Output = Toff * (on/off) - Toff
	1271	remove: Output = Tref * (on/off) - Ton
	1272	"""
	1273	modes = ["suffix","time"]
	1274	print mode
	1275	if not mode in modes:
	1276	print "please provide valid mode. Valid modes are %s" % (modes)
	1277	return None
	1278	from asap._asap import quotient as _quot
	1279	if mode == "suffix":
	1280	srcs = self.get_scan("*[^_ewR]")
	1281	refs = self.get_scan("*[_ewR]")
	1282	return scantable(_quot(srcs,refs, preserve))
	1283	else:
	1284	print "not yet implemented"
	1285	return None
	1286
[513]	1287	def quotient(self, other, isreference=True, preserve=True):
	1288	"""
[582]	1289	Return the quotient of a 'source' (on) scan and a 'reference' (off)
	1290	scan.
[513]	1291	The reference can have just one row, even if the signal has many.
	1292	Otherwise they must have the same number of rows.
	1293	The cursor of the output scan is set to 0
	1294	Parameters:
[582]	1295	other: the 'other' scan
[513]	1296	isreference: if the 'other' scan is the reference (default)
	1297	or source
	1298	preserve: you can preserve (default) the continuum or
	1299	remove it. The equations used are
[582]	1300	preserve: Output = Toff * (on/off) - Toff
	1301	remove: Output = Tref * (on/off) - Ton
[513]	1302	Example:
	1303	# src is a scantable for an 'on' scan, ref for an 'off' scantable
	1304	q1 = src.quotient(ref)
	1305	q2 = ref.quotient(src, isreference=False)
	1306	# gives the same result
	1307	"""
	1308	from asap._asap import quotient as _quot
	1309	if isreference:
	1310	return scantable(_quot(self, other, preserve))
	1311	else:
	1312	return scantable(_quot(other, self, preserve))
	1313
	1314	def __add__(self, other):
	1315	varlist = vars()
	1316	s = None
	1317	if isinstance(other, scantable):
	1318	from asap._asap import b_operate as _bop
	1319	s = scantable(_bop(self, other, 'add', True))
	1320	elif isinstance(other, float):
	1321	from asap._asap import add as _add
	1322	s = scantable(_add(self, other, True))
	1323	else:
	1324	print "Other input is not a scantable or float value"
	1325	return
	1326	s._add_history("operator +", varlist)
	1327	return s
	1328
	1329
	1330	def __sub__(self, other):
	1331	"""
	1332	implicit on all axes and on Tsys
	1333	"""
	1334	varlist = vars()
	1335	s = None
	1336	if isinstance(other, scantable):
	1337	from asap._asap import b_operate as _bop
	1338	s = scantable(_bop(self, other, 'sub', True))
	1339	elif isinstance(other, float):
	1340	from asap._asap import add as _add
	1341	s = scantable(_add(self, -other, True))
	1342	else:
	1343	print "Other input is not a scantable or float value"
	1344	return
	1345	s._add_history("operator -", varlist)
	1346	return s
	1347
	1348	def __mul__(self, other):
	1349	"""
	1350	implicit on all axes and on Tsys
	1351	"""
	1352	varlist = vars()
	1353	s = None
	1354	if isinstance(other, scantable):
	1355	from asap._asap import b_operate as _bop
	1356	s = scantable(_bop(self, other, 'mul', True))
	1357	elif isinstance(other, float):
	1358	if other == 0.0:
	1359	print "Multiplying by zero is not recommended."
	1360	return
	1361	from asap._asap import scale as _sca
	1362	s = scantable(_sca(self, other, True))
	1363	else:
	1364	print "Other input is not a scantable or float value"
	1365	return
	1366	s._add_history("operator *", varlist)
	1367	return s
	1368
	1369
	1370	def __div__(self, other):
	1371	"""
	1372	implicit on all axes and on Tsys
	1373	"""
	1374	varlist = vars()
	1375	s = None
	1376	if isinstance(other, scantable):
	1377	from asap._asap import b_operate as _bop
	1378	s = scantable(_bop(self, other, 'div', True))
	1379	elif isinstance(other, float):
	1380	if other == 0.0:
	1381	print "Dividing by zero is not recommended"
	1382	return
	1383	from asap._asap import scale as _sca
	1384	s = scantable(_sca(self, 1.0/other, True))
	1385	else:
	1386	print "Other input is not a scantable or float value"
	1387	return
	1388	s._add_history("operator /", varlist)
	1389	return s
	1390
[530]	1391	def get_fit(self, row=0):
	1392	"""
	1393	Print or return the stored fits for a row in the scantable
	1394	Parameters:
	1395	row: the row which the fit has been applied to.
	1396	"""
	1397	if row > self.nrow():
	1398	return
	1399	from asap import asapfit
	1400	fit = asapfit(self._getfit(row))
	1401	if self._vb:
	1402	print fit
	1403	return
	1404	else:
	1405	return fit.as_dict()
	1406
[484]	1407	def _add_history(self, funcname, parameters):
	1408	# create date
	1409	sep = "##"
	1410	from datetime import datetime
	1411	dstr = datetime.now().strftime('%Y/%m/%d %H:%M:%S')
	1412	hist = dstr+sep
	1413	hist += funcname+sep#cdate+sep
	1414	if parameters.has_key('self'): del parameters['self']
	1415	for k,v in parameters.iteritems():
	1416	if type(v) is dict:
	1417	for k2,v2 in v.iteritems():
	1418	hist += k2
	1419	hist += "="
	1420	if isinstance(v2,scantable):
	1421	hist += 'scantable'
	1422	elif k2 == 'mask':
[513]	1423	if isinstance(v2,list) or isinstance(v2,tuple):
	1424	hist += str(self._zip_mask(v2))
	1425	else:
	1426	hist += str(v2)
[484]	1427	else:
[513]	1428	hist += str(v2)
[484]	1429	else:
	1430	hist += k
	1431	hist += "="
	1432	if isinstance(v,scantable):
	1433	hist += 'scantable'
	1434	elif k == 'mask':
[513]	1435	if isinstance(v,list) or isinstance(v,tuple):
	1436	hist += str(self._zip_mask(v))
	1437	else:
	1438	hist += str(v)
[484]	1439	else:
	1440	hist += str(v)
	1441	hist += sep
	1442	hist = hist[:-2] # remove trailing '##'
	1443	self._addhistory(hist)
	1444
	1445
	1446	def _zip_mask(self, mask):
	1447	mask = list(mask)
	1448	i = 0
	1449	segments = []
	1450	while mask[i:].count(1):
	1451	i += mask[i:].index(1)
	1452	if mask[i:].count(0):
	1453	j = i + mask[i:].index(0)
	1454	else:
	1455	j = len(mask)
	1456	segments.append([i,j])
[513]	1457	i = j
[484]	1458	return segments
[626]	1459	def _get_ordinate_label(self):
	1460	fu = "("+self.get_fluxunit()+")"
	1461	import re
	1462	lbl = "Intensity"
	1463	if re.match(".K.",fu):
	1464	lbl = "Brightness Temperature "+ fu
	1465	elif re.match(".Jy.",fu):
	1466	lbl = "Flux density "+ fu
	1467	return lbl
	1468

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

Context Navigation

source: branches/Release-2-fixes/python/scantable.py@ 1751

Download in other formats: