Context Navigation

source: trunk/python/scantable.py@ 729

Visit:

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

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

Download in other formats: