Context Navigation

source: branches/Release12/python/scantable.py@ 787

Visit:

Last change on this file since 787 was 787, checked in by phi196, 19 years ago
Added get_elevation, azimuth and parangle
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 60.9 KB

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

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

Download in other formats: