Context Navigation

source: trunk/python/scantable.py@ 799

Visit:

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

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

Download in other formats: