Context Navigation

asapmath.py@ 1738

Visit:

Last change on this file since 1738 was 1693, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-1908

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes

Module(s): Module Names change impacts.

Description: Describe your changes here...

Changed a tagging as the source type is tagged by using SRCTYPE, not
an extra string in SRCNAME. To do this, I have defined a selection method
by SRCTYPE in STSelector class. I have newly added python_SrcType.cpp
that defines a Python interface of SrcType enums which is defined
in atnf/PKSIO/SrcType.h.

Since I have added new file in the src directory, I have modified src/Makefile
to compile new file.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 39.6 KB

Rev	Line
[1085]	1	from asap.scantable import scantable
[258]	2	from asap import rcParams
[720]	3	from asap import print_log
[1389]	4	from asap import selector
[1614]	5	from asap import asaplog
[1631]	6	from asap import asaplotgui
[101]	7
[143]	8	def average_time(args, *kwargs):
[101]	9	"""
[113]	10	Return the (time) average of a scan or list of scans. [in channels only]
[305]	11	The cursor of the output scan is set to 0
[113]	12	Parameters:
[1361]	13	one scan or comma separated scans or a list of scans
[143]	14	mask: an optional mask (only used for 'var' and 'tsys' weighting)
[558]	15	scanav: True averages each scan separately.
	16	False (default) averages all scans together,
[1232]	17	weight: Weighting scheme.
	18	'none' (mean no weight)
	19	'var' (1/var(spec) weighted)
	20	'tsys' (1/Tsys**2 weighted)
	21	'tint' (integration time weighted)
	22	'tintsys' (Tint/Tsys**2)
	23	'median' ( median averaging)
[930]	24	align: align the spectra in velocity before averaging. It takes
	25	the time of the first spectrum in the first scantable
	26	as reference time.
[113]	27	Example:
	28	# return a time averaged scan from scana and scanb
	29	# without using a mask
[129]	30	scanav = average_time(scana,scanb)
[1361]	31	# or equivalent
	32	# scanav = average_time([scana, scanb])
[113]	33	# return the (time) averaged scan, i.e. the average of
	34	# all correlator cycles
[558]	35	scanav = average_time(scan, scanav=True)
[101]	36	"""
[930]	37	scanav = False
[143]	38	if kwargs.has_key('scanav'):
[930]	39	scanav = kwargs.get('scanav')
[524]	40	weight = 'tint'
[143]	41	if kwargs.has_key('weight'):
	42	weight = kwargs.get('weight')
	43	mask = ()
	44	if kwargs.has_key('mask'):
	45	mask = kwargs.get('mask')
[930]	46	align = False
	47	if kwargs.has_key('align'):
	48	align = kwargs.get('align')
[1446]	49	compel = False
	50	if kwargs.has_key('compel'):
	51	compel = kwargs.get('compel')
[489]	52	varlist = vars()
[665]	53	if isinstance(args[0],list):
[981]	54	lst = args[0]
[665]	55	elif isinstance(args[0],tuple):
[981]	56	lst = list(args[0])
[665]	57	else:
[981]	58	lst = list(args)
[720]	59
[489]	60	del varlist["kwargs"]
	61	varlist["args"] = "%d scantables" % len(lst)
[981]	62	# need special formatting here for history...
[720]	63
[876]	64	from asap._asap import stmath
	65	stm = stmath()
[113]	66	for s in lst:
[101]	67	if not isinstance(s,scantable):
[720]	68	msg = "Please give a list of scantables"
	69	if rcParams['verbose']:
[1612]	70	#print msg
[1614]	71	asaplog.push(msg)
	72	print_log('ERROR')
[720]	73	return
	74	else:
	75	raise TypeError(msg)
[945]	76	if scanav: scanav = "SCAN"
	77	else: scanav = "NONE"
[981]	78	alignedlst = []
	79	if align:
	80	refepoch = lst[0].get_time(0)
	81	for scan in lst:
	82	alignedlst.append(scan.freq_align(refepoch,insitu=False))
	83	else:
[1059]	84	alignedlst = lst
[1232]	85	if weight.upper() == 'MEDIAN':
	86	# median doesn't support list of scantables - merge first
	87	merged = None
	88	if len(alignedlst) > 1:
	89	merged = merge(alignedlst)
	90	else:
	91	merged = alignedlst[0]
	92	s = scantable(stm._averagechannel(merged, 'MEDIAN', scanav))
	93	del merged
	94	else:
[1446]	95	#s = scantable(stm._average(alignedlst, mask, weight.upper(), scanav))
	96	s = scantable(stm._new_average(alignedlst, compel, mask, weight.upper(), scanav))
[489]	97	s._add_history("average_time",varlist)
[720]	98	print_log()
[489]	99	return s
[101]	100
[1074]	101	def quotient(source, reference, preserve=True):
	102	"""
	103	Return the quotient of a 'source' (signal) scan and a 'reference' scan.
	104	The reference can have just one scan, even if the signal has many. Otherwise
	105	they must have the same number of scans.
	106	The cursor of the output scan is set to 0
	107	Parameters:
	108	source: the 'on' scan
	109	reference: the 'off' scan
	110	preserve: you can preserve (default) the continuum or
	111	remove it. The equations used are
	112	preserve: Output = Toff * (on/off) - Toff
	113	remove: Output = Toff * (on/off) - Ton
	114	"""
	115	varlist = vars()
	116	from asap._asap import stmath
	117	stm = stmath()
	118	stm._setinsitu(False)
	119	s = scantable(stm._quotient(source, reference, preserve))
	120	s._add_history("quotient",varlist)
	121	print_log()
	122	return s
[101]	123
[1389]	124	def dototalpower(calon, caloff, tcalval=0.0):
	125	"""
	126	Do calibration for CAL on,off signals.
	127	Adopted from GBTIDL dototalpower
	128	Parameters:
	129	calon: the 'cal on' subintegration
	130	caloff: the 'cal off' subintegration
	131	tcalval: user supplied Tcal value
	132	"""
	133	varlist = vars()
	134	from asap._asap import stmath
	135	stm = stmath()
	136	stm._setinsitu(False)
	137	s = scantable(stm._dototalpower(calon, caloff, tcalval))
	138	s._add_history("dototalpower",varlist)
	139	print_log()
	140	return s
	141
	142	def dosigref(sig, ref, smooth, tsysval=0.0, tauval=0.0):
	143	"""
	144	Calculate a quotient (sig-ref/ref * Tsys)
	145	Adopted from GBTIDL dosigref
	146	Parameters:
	147	sig: on source data
	148	ref: reference data
	149	smooth: width of box car smoothing for reference
	150	tsysval: user specified Tsys (scalar only)
	151	tauval: user specified Tau (required if tsysval is set)
	152	"""
	153	varlist = vars()
	154	from asap._asap import stmath
	155	stm = stmath()
	156	stm._setinsitu(False)
	157	s = scantable(stm._dosigref(sig, ref, smooth, tsysval, tauval))
	158	s._add_history("dosigref",varlist)
	159	print_log()
	160	return s
	161
[1631]	162	def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
[1389]	163	"""
	164	Calibrate GBT position switched data
	165	Adopted from GBTIDL getps
	166	Currently calps identify the scans as position switched data if they
	167	contain '_ps' in the source name. The data must contains 'CAL' signal
	168	on/off in each integration. To identify 'CAL' on state, the word, 'calon'
	169	need to be present in the source name field.
	170	(GBT MS data reading process to scantable automatically append these
	171	id names to the source names)
	172
	173	Parameters:
	174	scantab: scantable
	175	scannos: list of scan numbers
	176	smooth: optional box smoothing order for the reference
	177	(default is 1 = no smoothing)
	178	tsysval: optional user specified Tsys (default is 0.0,
	179	use Tsys in the data)
	180	tauval: optional user specified Tau
	181	tcalval: optional user specified Tcal (default is 0.0,
	182	use Tcal value in the data)
	183	"""
	184	varlist = vars()
	185	# check for the appropriate data
[1693]	186	## s = scantab.get_scan('_ps')
	187	## if s is None:
	188	## msg = "The input data appear to contain no position-switch mode data."
	189	## if rcParams['verbose']:
	190	## #print msg
	191	## asaplog.push(msg)
	192	## print_log('ERROR')
	193	## return
	194	## else:
	195	## raise TypeError(msg)
	196	s = scantab.copy()
	197	from asap._asap import srctype
	198	sel = selector()
	199	sel.set_types( srctype.pson )
	200	try:
	201	scantab.set_selection( sel )
	202	except Exception, e:
[1389]	203	msg = "The input data appear to contain no position-switch mode data."
	204	if rcParams['verbose']:
[1612]	205	#print msg
[1614]	206	asaplog.push(msg)
	207	print_log('ERROR')
[1389]	208	return
	209	else:
	210	raise TypeError(msg)
[1693]	211	s.set_selection()
	212	sel.reset()
[1389]	213	ssub = s.get_scan(scannos)
	214	if ssub is None:
	215	msg = "No data was found with given scan numbers!"
	216	if rcParams['verbose']:
[1612]	217	#print msg
[1614]	218	asaplog.push(msg)
	219	print_log('ERROR')
[1389]	220	return
	221	else:
	222	raise TypeError(msg)
[1693]	223	#ssubon = ssub.get_scan('*calon')
	224	#ssuboff = ssub.get_scan('*[^calon]')
	225	sel.set_types( [srctype.poncal,srctype.poffcal] )
	226	ssub.set_selection( sel )
	227	ssubon = ssub.copy()
	228	ssub.set_selection()
	229	sel.reset()
	230	sel.set_types( [srctype.pson,srctype.psoff] )
	231	ssub.set_selection( sel )
	232	ssuboff = ssub.copy()
	233	ssub.set_selection()
	234	sel.reset()
[1389]	235	if ssubon.nrow() != ssuboff.nrow():
	236	msg = "mismatch in numbers of CAL on/off scans. Cannot calibrate. Check the scan numbers."
	237	if rcParams['verbose']:
[1612]	238	#print msg
[1614]	239	asaplog.push(msg)
	240	print_log('ERROR')
[1389]	241	return
	242	else:
	243	raise TypeError(msg)
	244	cals = dototalpower(ssubon, ssuboff, tcalval)
[1693]	245	#sig = cals.get_scan('*ps')
	246	#ref = cals.get_scan('*psr')
	247	sel.set_types( srctype.pson )
	248	cals.set_selection( sel )
	249	sig = cals.copy()
	250	cals.set_selection()
	251	sel.reset()
	252	sel.set_types( srctype.psoff )
	253	cals.set_selection( sel )
	254	ref = cals.copy()
	255	cals.set_selection()
	256	sel.reset()
[1389]	257	if sig.nscan() != ref.nscan():
	258	msg = "mismatch in numbers of on/off scans. Cannot calibrate. Check the scan numbers."
	259	if rcParams['verbose']:
[1612]	260	#print msg
[1614]	261	asaplog.push(msg)
	262	print_log('ERROR')
[1389]	263	return
	264	else:
	265	raise TypeError(msg)
	266
	267	#for user supplied Tsys
	268	if tsysval>0.0:
	269	if tauval<=0.0:
	270	msg = "Need to supply a valid tau to use the supplied Tsys"
	271	if rcParams['verbose']:
[1612]	272	#print msg
[1614]	273	asaplog.push(msg)
	274	print_log('ERROR')
[1389]	275	return
	276	else:
	277	raise TypeError(msg)
	278	else:
	279	sig.recalc_azel()
	280	ref.recalc_azel()
	281	#msg = "Use of user specified Tsys is not fully implemented yet."
	282	#if rcParams['verbose']:
	283	# print msg
	284	# return
	285	#else:
	286	# raise TypeError(msg)
	287	# use get_elevation to get elevation and
	288	# calculate a scaling factor using the formula
	289	# -> tsys use to dosigref
	290
	291	#ress = dosigref(sig, ref, smooth, tsysval)
	292	ress = dosigref(sig, ref, smooth, tsysval, tauval)
[1631]	293	###
	294	if verify:
	295	# get data
	296	import numpy
	297	precal={}
	298	postcal=[]
	299	keys=['ps','ps_calon','psr','psr_calon']
[1693]	300	types=[srctype.pson,srctype.poncal,srctype.psoff,srctype.poffcal]
[1631]	301	ifnos=list(ssub.getifnos())
	302	polnos=list(ssub.getpolnos())
	303	sel=selector()
	304	for i in range(2):
[1693]	305	#ss=ssuboff.get_scan(''+keys[2i])
[1631]	306	ll=[]
	307	for j in range(len(ifnos)):
	308	for k in range(len(polnos)):
	309	sel.set_ifs(ifnos[j])
	310	sel.set_polarizations(polnos[k])
[1693]	311	sel.set_types(types[2*i])
[1631]	312	try:
[1693]	313	#ss.set_selection(sel)
	314	ssuboff.set_selection(sel)
[1631]	315	except:
	316	continue
[1693]	317	#ll.append(numpy.array(ss._getspectrum(0)))
	318	ll.append(numpy.array(ssuboff._getspectrum(0)))
[1631]	319	sel.reset()
[1693]	320	ssuboff.set_selection()
[1631]	321	precal[keys[2*i]]=ll
[1693]	322	#del ss
	323	#ss=ssubon.get_scan(''+keys[2i+1])
[1631]	324	ll=[]
	325	for j in range(len(ifnos)):
	326	for k in range(len(polnos)):
	327	sel.set_ifs(ifnos[j])
	328	sel.set_polarizations(polnos[k])
[1693]	329	sel.set_types(types[2*i+1])
[1631]	330	try:
[1693]	331	#ss.set_selection(sel)
	332	ssubon.set_selection(sel)
[1631]	333	except:
	334	continue
[1693]	335	#ll.append(numpy.array(ss._getspectrum(0)))
	336	ll.append(numpy.array(ssubon._getspectrum(0)))
[1631]	337	sel.reset()
[1693]	338	ssubon.set_selection()
[1631]	339	precal[keys[2*i+1]]=ll
[1693]	340	#del ss
[1631]	341	for j in range(len(ifnos)):
	342	for k in range(len(polnos)):
	343	sel.set_ifs(ifnos[j])
	344	sel.set_polarizations(polnos[k])
	345	try:
	346	ress.set_selection(sel)
	347	except:
	348	continue
	349	postcal.append(numpy.array(ress._getspectrum(0)))
	350	sel.reset()
	351	ress.set_selection()
	352	del sel
	353	# plot
	354	print_log()
	355	asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
	356	print_log('WARN')
	357	p=asaplotgui.asaplotgui()
	358	#nr=min(6,len(ifnos)*len(polnos))
	359	nr=len(ifnos)*len(polnos)
	360	titles=[]
	361	btics=[]
	362	if nr<4:
	363	p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False)
	364	for i in range(2*nr):
	365	b=False
	366	if i >= 2*nr-2:
	367	b=True
	368	btics.append(b)
	369	elif nr==4:
	370	p.set_panels(rows=2,cols=4,nplots=8,ganged=False)
	371	for i in range(2*nr):
	372	b=False
	373	if i >= 2*nr-4:
	374	b=True
	375	btics.append(b)
	376	elif nr<7:
	377	p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
	378	for i in range(2*nr):
	379	if i >= 2*nr-4:
	380	b=True
	381	btics.append(b)
	382	else:
	383	print_log()
	384	asaplog.push('Only first 6 [if,pol] pairs are plotted.')
	385	print_log('WARN')
	386	nr=6
	387	for i in range(2*nr):
	388	b=False
	389	if i >= 2*nr-4:
	390	b=True
	391	btics.append(b)
	392	p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
	393	for i in range(nr):
	394	p.subplot(2*i)
	395	p.color=0
	396	title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
	397	titles.append(title)
	398	#p.set_axes('title',title,fontsize=40)
	399	ymin=1.0e100
	400	ymax=-1.0e100
	401	nchan=s.nchan()
	402	edge=int(nchan*0.01)
	403	for j in range(4):
	404	spmin=min(precal[keys[j]][i][edge:nchan-edge])
	405	spmax=max(precal[keys[j]][i][edge:nchan-edge])
	406	ymin=min(ymin,spmin)
	407	ymax=max(ymax,spmax)
	408	for j in range(4):
	409	if i==0:
	410	p.set_line(label=keys[j])
	411	else:
	412	p.legend()
	413	p.plot(precal[keys[j]][i])
	414	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	415	if not btics[2*i]:
	416	p.axes.set_xticks([])
	417	p.subplot(2*i+1)
	418	p.color=0
	419	title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
	420	titles.append(title)
	421	#p.set_axes('title',title)
	422	p.legend()
	423	ymin=postcal[i][edge:nchan-edge].min()
	424	ymax=postcal[i][edge:nchan-edge].max()
	425	p.plot(postcal[i])
	426	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	427	if not btics[2*i+1]:
	428	p.axes.set_xticks([])
	429	for i in range(2*nr):
	430	p.subplot(i)
	431	p.set_axes('title',titles[i],fontsize='medium')
	432	x=raw_input('Accept calibration ([y]/n): ' )
	433	if x.upper() == 'N':
	434	p.unmap()
	435	del p
	436	return scabtab
	437	p.unmap()
	438	del p
	439	###
[1389]	440	ress._add_history("calps", varlist)
	441	print_log()
	442	return ress
	443
[1631]	444	def calnod(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
[1389]	445	"""
	446	Do full (but a pair of scans at time) processing of GBT Nod data
	447	calibration.
	448	Adopted from GBTIDL's getnod
	449	Parameters:
	450	scantab: scantable
	451	scannos: a pair of scan numbers, or the first scan number of the pair
	452	smooth: box car smoothing order
	453	tsysval: optional user specified Tsys value
	454	tauval: optional user specified tau value (not implemented yet)
	455	tcalval: optional user specified Tcal value
	456	"""
	457	varlist = vars()
	458	from asap._asap import stmath
[1693]	459	from asap._asap import srctype
[1389]	460	stm = stmath()
	461	stm._setinsitu(False)
	462
	463	# check for the appropriate data
[1693]	464	## s = scantab.get_scan('_nod')
	465	## if s is None:
	466	## msg = "The input data appear to contain no Nod observing mode data."
	467	## if rcParams['verbose']:
	468	## #print msg
	469	## asaplog.push(msg)
	470	## print_log('ERROR')
	471	## return
	472	## else:
	473	## raise TypeError(msg)
	474	s = scantab.copy()
	475	sel = selector()
	476	sel.set_types( srctype.nod )
	477	try:
	478	s.set_selection( sel )
	479	except Exception, e:
[1389]	480	msg = "The input data appear to contain no Nod observing mode data."
	481	if rcParams['verbose']:
[1612]	482	#print msg
[1614]	483	asaplog.push(msg)
	484	print_log('ERROR')
[1389]	485	return
	486	else:
	487	raise TypeError(msg)
[1693]	488	sel.reset()
	489	del sel
	490	del s
[1389]	491
	492	# need check correspondance of each beam with sig-ref ...
	493	# check for timestamps, scan numbers, subscan id (not available in
	494	# ASAP data format...). Assume 1st scan of the pair (beam 0 - sig
	495	# and beam 1 - ref...)
	496	# First scan number of paired scans or list of pairs of
	497	# scan numbers (has to have even number of pairs.)
	498
	499	#data splitting
	500	scan1no = scan2no = 0
	501
	502	if len(scannos)==1:
	503	scan1no = scannos[0]
	504	scan2no = scannos[0]+1
	505	pairScans = [scan1no, scan2no]
	506	else:
	507	#if len(scannos)>2:
	508	# msg = "calnod can only process a pair of nod scans at time."
	509	# if rcParams['verbose']:
	510	# print msg
	511	# return
	512	# else:
	513	# raise TypeError(msg)
	514	#
	515	#if len(scannos)==2:
	516	# scan1no = scannos[0]
	517	# scan2no = scannos[1]
	518	pairScans = list(scannos)
	519
	520	if tsysval>0.0:
	521	if tauval<=0.0:
	522	msg = "Need to supply a valid tau to use the supplied Tsys"
	523	if rcParams['verbose']:
[1612]	524	#print msg
[1614]	525	asaplog.push(msg)
	526	print_log('ERROR')
[1389]	527	return
	528	else:
	529	raise TypeError(msg)
	530	else:
	531	scantab.recalc_azel()
	532	resspec = scantable(stm._donod(scantab, pairScans, smooth, tsysval,tauval,tcalval))
[1631]	533	###
	534	if verify:
	535	# get data
	536	import numpy
	537	precal={}
	538	postcal=[]
[1693]	539	keys=['','_calon']
	540	types=[srctype.nod,srctype.nodcal]
[1631]	541	ifnos=list(scantab.getifnos())
	542	polnos=list(scantab.getpolnos())
	543	sel=selector()
[1693]	544	ss = scantab.copy()
[1631]	545	for i in range(2):
[1693]	546	#ss=scantab.get_scan('*'+keys[i])
[1631]	547	ll=[]
	548	ll2=[]
	549	for j in range(len(ifnos)):
	550	for k in range(len(polnos)):
	551	sel.set_ifs(ifnos[j])
	552	sel.set_polarizations(polnos[k])
	553	sel.set_scans(pairScans[0])
[1693]	554	sel.set_types(types[i])
[1631]	555	try:
	556	ss.set_selection(sel)
	557	except:
	558	continue
	559	ll.append(numpy.array(ss._getspectrum(0)))
	560	sel.reset()
	561	ss.set_selection()
	562	sel.set_ifs(ifnos[j])
	563	sel.set_polarizations(polnos[k])
	564	sel.set_scans(pairScans[1])
[1693]	565	sel.set_types(types[i])
[1631]	566	try:
	567	ss.set_selection(sel)
	568	except:
	569	ll.pop()
	570	continue
	571	ll2.append(numpy.array(ss._getspectrum(0)))
	572	sel.reset()
	573	ss.set_selection()
[1693]	574	key='%s%s' %(pairScans[0],keys[i])
[1631]	575	precal[key]=ll
[1693]	576	key='%s%s' %(pairScans[1],keys[i])
[1631]	577	precal[key]=ll2
[1693]	578	#del ss
[1631]	579	keys=precal.keys()
	580	for j in range(len(ifnos)):
	581	for k in range(len(polnos)):
	582	sel.set_ifs(ifnos[j])
	583	sel.set_polarizations(polnos[k])
	584	sel.set_scans(pairScans[0])
	585	try:
	586	resspec.set_selection(sel)
	587	except:
	588	continue
	589	postcal.append(numpy.array(resspec._getspectrum(0)))
	590	sel.reset()
	591	resspec.set_selection()
	592	del sel
	593	# plot
	594	print_log()
	595	asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
	596	print_log('WARN')
	597	p=asaplotgui.asaplotgui()
	598	#nr=min(6,len(ifnos)*len(polnos))
	599	nr=len(ifnos)*len(polnos)
	600	titles=[]
	601	btics=[]
	602	if nr<4:
	603	p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False)
	604	for i in range(2*nr):
	605	b=False
	606	if i >= 2*nr-2:
	607	b=True
	608	btics.append(b)
	609	elif nr==4:
	610	p.set_panels(rows=2,cols=4,nplots=8,ganged=False)
	611	for i in range(2*nr):
	612	b=False
	613	if i >= 2*nr-4:
	614	b=True
	615	btics.append(b)
	616	elif nr<7:
	617	p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
	618	for i in range(2*nr):
	619	if i >= 2*nr-4:
	620	b=True
	621	btics.append(b)
	622	else:
	623	print_log()
	624	asaplog.push('Only first 6 [if,pol] pairs are plotted.')
	625	print_log('WARN')
	626	nr=6
	627	for i in range(2*nr):
	628	b=False
	629	if i >= 2*nr-4:
	630	b=True
	631	btics.append(b)
	632	p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
	633	for i in range(nr):
	634	p.subplot(2*i)
	635	p.color=0
	636	title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
	637	titles.append(title)
	638	#p.set_axes('title',title,fontsize=40)
	639	ymin=1.0e100
	640	ymax=-1.0e100
	641	nchan=scantab.nchan()
	642	edge=int(nchan*0.01)
	643	for j in range(4):
	644	spmin=min(precal[keys[j]][i][edge:nchan-edge])
	645	spmax=max(precal[keys[j]][i][edge:nchan-edge])
	646	ymin=min(ymin,spmin)
	647	ymax=max(ymax,spmax)
	648	for j in range(4):
	649	if i==0:
	650	p.set_line(label=keys[j])
	651	else:
	652	p.legend()
	653	p.plot(precal[keys[j]][i])
	654	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	655	if not btics[2*i]:
	656	p.axes.set_xticks([])
	657	p.subplot(2*i+1)
	658	p.color=0
	659	title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
	660	titles.append(title)
	661	#p.set_axes('title',title)
	662	p.legend()
	663	ymin=postcal[i][edge:nchan-edge].min()
	664	ymax=postcal[i][edge:nchan-edge].max()
	665	p.plot(postcal[i])
	666	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	667	if not btics[2*i+1]:
	668	p.axes.set_xticks([])
	669	for i in range(2*nr):
	670	p.subplot(i)
	671	p.set_axes('title',titles[i],fontsize='medium')
	672	x=raw_input('Accept calibration ([y]/n): ' )
	673	if x.upper() == 'N':
	674	p.unmap()
	675	del p
	676	return scabtab
	677	p.unmap()
	678	del p
	679	###
[1389]	680	resspec._add_history("calnod",varlist)
	681	print_log()
	682	return resspec
	683
[1631]	684	def calfs(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
[1389]	685	"""
	686	Calibrate GBT frequency switched data.
	687	Adopted from GBTIDL getfs.
	688	Currently calfs identify the scans as frequency switched data if they
	689	contain '_fs' in the source name. The data must contains 'CAL' signal
	690	on/off in each integration. To identify 'CAL' on state, the word, 'calon'
	691	need to be present in the source name field.
	692	(GBT MS data reading via scantable automatically append these
	693	id names to the source names)
	694
	695	Parameters:
	696	scantab: scantable
	697	scannos: list of scan numbers
	698	smooth: optional box smoothing order for the reference
	699	(default is 1 = no smoothing)
	700	tsysval: optional user specified Tsys (default is 0.0,
	701	use Tsys in the data)
	702	tauval: optional user specified Tau
	703	"""
	704	varlist = vars()
	705	from asap._asap import stmath
[1693]	706	from asap._asap import srctype
[1389]	707	stm = stmath()
	708	stm._setinsitu(False)
	709
	710	# check = scantab.get_scan('_fs')
	711	# if check is None:
	712	# msg = "The input data appear to contain no Nod observing mode data."
	713	# if rcParams['verbose']:
	714	# print msg
	715	# return
	716	# else:
	717	# raise TypeError(msg)
	718	s = scantab.get_scan(scannos)
	719	del scantab
	720
	721	resspec = scantable(stm._dofs(s, scannos, smooth, tsysval,tauval,tcalval))
[1631]	722	###
	723	if verify:
	724	# get data
	725	ssub = s.get_scan(scannos)
[1693]	726	#ssubon = ssub.get_scan('*calon')
	727	#ssuboff = ssub.get_scan('*[^calon]')
	728	sel = selector()
	729	sel.set_types( [srctype.foncal,srctype.foffcal] )
	730	ssub.set_selection( sel )
	731	ssubon = ssub.copy()
	732	ssub.set_selection()
	733	sel.reset()
	734	sel.set_types( [srctype.fson,srctype.fsoff] )
	735	ssub.set_selection( sel )
	736	ssuboff = ssub.copy()
	737	ssub.set_selection()
	738	sel.reset()
[1631]	739	import numpy
	740	precal={}
	741	postcal=[]
	742	keys=['fs','fs_calon','fsr','fsr_calon']
[1693]	743	types=[srctype.fson,srctype.foncal,srctype.fsoff,srctype.foffcal]
[1631]	744	ifnos=list(ssub.getifnos())
	745	polnos=list(ssub.getpolnos())
	746	for i in range(2):
[1693]	747	#ss=ssuboff.get_scan(''+keys[2i])
[1631]	748	ll=[]
	749	for j in range(len(ifnos)):
	750	for k in range(len(polnos)):
	751	sel.set_ifs(ifnos[j])
	752	sel.set_polarizations(polnos[k])
[1693]	753	sel.set_types(types[2*i])
[1631]	754	try:
[1693]	755	#ss.set_selection(sel)
	756	ssuboff.set_selection(sel)
[1631]	757	except:
	758	continue
	759	ll.append(numpy.array(ss._getspectrum(0)))
	760	sel.reset()
[1693]	761	#ss.set_selection()
	762	ssuboff.set_selection()
[1631]	763	precal[keys[2*i]]=ll
[1693]	764	#del ss
	765	#ss=ssubon.get_scan(''+keys[2i+1])
[1631]	766	ll=[]
	767	for j in range(len(ifnos)):
	768	for k in range(len(polnos)):
	769	sel.set_ifs(ifnos[j])
	770	sel.set_polarizations(polnos[k])
[1693]	771	sel.set_types(types[2*i+1])
[1631]	772	try:
[1693]	773	#ss.set_selection(sel)
	774	ssubon.set_selection(sel)
[1631]	775	except:
	776	continue
	777	ll.append(numpy.array(ss._getspectrum(0)))
	778	sel.reset()
[1693]	779	#ss.set_selection()
	780	ssubon.set_selection()
[1631]	781	precal[keys[2*i+1]]=ll
[1693]	782	#del ss
	783	#sig=resspec.get_scan('*_fs')
	784	#ref=resspec.get_scan('*_fsr')
	785	sel.set_types( srctype.fson )
	786	resspec.set_selection( sel )
	787	sig=resspec.copy()
	788	resspec.set_selection()
	789	sel.reset()
	790	sel.set_type( srctype.fsoff )
	791	resspec.set_selection( sel )
	792	ref=resspec.copy()
	793	resspec.set_selection()
	794	sel.reset()
[1631]	795	for k in range(len(polnos)):
	796	for j in range(len(ifnos)):
	797	sel.set_ifs(ifnos[j])
	798	sel.set_polarizations(polnos[k])
	799	try:
	800	sig.set_selection(sel)
	801	postcal.append(numpy.array(sig._getspectrum(0)))
	802	except:
	803	ref.set_selection(sel)
	804	postcal.append(numpy.array(ref._getspectrum(0)))
	805	sel.reset()
	806	resspec.set_selection()
	807	del sel
	808	# plot
	809	print_log()
	810	asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
	811	print_log('WARN')
	812	p=asaplotgui.asaplotgui()
	813	#nr=min(6,len(ifnos)*len(polnos))
	814	nr=len(ifnos)/2*len(polnos)
	815	titles=[]
	816	btics=[]
	817	if nr>3:
	818	print_log()
	819	asaplog.push('Only first 3 [if,pol] pairs are plotted.')
	820	print_log('WARN')
	821	nr=3
	822	p.set_panels(rows=nr,cols=3,nplots=3*nr,ganged=False)
	823	for i in range(3*nr):
	824	b=False
	825	if i >= 3*nr-3:
	826	b=True
	827	btics.append(b)
	828	for i in range(nr):
	829	p.subplot(3*i)
	830	p.color=0
	831	title='raw data IF%s,%s POL%s' % (ifnos[2int(i/len(polnos))],ifnos[2int(i/len(polnos))+1],polnos[i%len(polnos)])
	832	titles.append(title)
	833	#p.set_axes('title',title,fontsize=40)
	834	ymin=1.0e100
	835	ymax=-1.0e100
	836	nchan=s.nchan()
	837	edge=int(nchan*0.01)
	838	for j in range(4):
	839	spmin=min(precal[keys[j]][i][edge:nchan-edge])
	840	spmax=max(precal[keys[j]][i][edge:nchan-edge])
	841	ymin=min(ymin,spmin)
	842	ymax=max(ymax,spmax)
	843	for j in range(4):
	844	if i==0:
	845	p.set_line(label=keys[j])
	846	else:
	847	p.legend()
	848	p.plot(precal[keys[j]][i])
	849	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	850	if not btics[3*i]:
	851	p.axes.set_xticks([])
	852	p.subplot(3*i+1)
	853	p.color=0
	854	title='sig data IF%s POL%s' % (ifnos[2*int(i/len(polnos))],polnos[i%len(polnos)])
	855	titles.append(title)
	856	#p.set_axes('title',title)
	857	p.legend()
	858	ymin=postcal[2*i][edge:nchan-edge].min()
	859	ymax=postcal[2*i][edge:nchan-edge].max()
	860	p.plot(postcal[2*i])
	861	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	862	if not btics[3*i+1]:
	863	p.axes.set_xticks([])
	864	p.subplot(3*i+2)
	865	p.color=0
	866	title='ref data IF%s POL%s' % (ifnos[2*int(i/len(polnos))+1],polnos[i%len(polnos)])
	867	titles.append(title)
	868	#p.set_axes('title',title)
	869	p.legend()
	870	ymin=postcal[2*i+1][edge:nchan-edge].min()
	871	ymax=postcal[2*i+1][edge:nchan-edge].max()
	872	p.plot(postcal[2*i+1])
	873	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	874	if not btics[3*i+2]:
	875	p.axes.set_xticks([])
	876	for i in range(3*nr):
	877	p.subplot(i)
	878	p.set_axes('title',titles[i],fontsize='medium')
	879	x=raw_input('Accept calibration ([y]/n): ' )
	880	if x.upper() == 'N':
	881	p.unmap()
	882	del p
	883	return scabtab
	884	p.unmap()
	885	del p
	886	###
[1389]	887	resspec._add_history("calfs",varlist)
	888	print_log()
	889	return resspec
	890
[296]	891	def simple_math(left, right, op='add', tsys=True):
[242]	892	"""
[720]	893	Apply simple mathematical binary operations to two
[242]	894	scan tables, returning the result in a new scan table.
	895	The operation is applied to both the correlations and the TSys data
[305]	896	The cursor of the output scan is set to 0
[242]	897	Parameters:
	898	left: the 'left' scan
	899	right: the 'right' scan
	900	op: the operation: 'add' (default), 'sub', 'mul', 'div'
[296]	901	tsys: if True (default) then apply the operation to Tsys
	902	as well as the data
[242]	903	"""
[1612]	904	#print "simple_math is deprecated use +=/* instead."
[1614]	905	asaplog.push( "simple_math is deprecated use +=/* instead." )
	906	print_log('WARN')
[918]	907
	908	def merge(*args):
[945]	909	"""
[1362]	910	Merge a list of scanatables, or comma-sperated scantables into one
	911	scnatble.
	912	Parameters:
	913	A list [scan1, scan2] or scan1, scan2.
	914	Example:
	915	myscans = [scan1, scan2]
	916	allscans = merge(myscans)
	917	# or equivalent
	918	sameallscans = merge(scan1, scan2)
[945]	919	"""
[918]	920	varlist = vars()
	921	if isinstance(args[0],list):
	922	lst = tuple(args[0])
	923	elif isinstance(args[0],tuple):
	924	lst = args[0]
	925	else:
	926	lst = tuple(args)
	927	varlist["args"] = "%d scantables" % len(lst)
	928	# need special formatting her for history...
	929	from asap._asap import stmath
	930	stm = stmath()
	931	for s in lst:
	932	if not isinstance(s,scantable):
	933	msg = "Please give a list of scantables"
	934	if rcParams['verbose']:
[1612]	935	#print msg
[1614]	936	asaplog.push(msg)
	937	print_log('ERROR')
[918]	938	return
	939	else:
	940	raise TypeError(msg)
	941	s = scantable(stm._merge(lst))
	942	s._add_history("merge", varlist)
	943	print_log()
	944	return s
[1074]	945
[1633]	946	def calibrate( scantab, scannos=[], calmode='none', verify=None ):
	947	"""
	948	Calibrate data.
[1631]	949
[1633]	950	Parameters:
	951	scantab: scantable
	952	scannos: list of scan number
	953	calmode: calibration mode
	954	verify: verify calibration
	955	"""
	956	antname = scantab.get_antennaname()
	957	if ( calmode == 'nod' ):
	958	asaplog.push( 'Calibrating nod data.' )
	959	print_log()
	960	scal = calnod( scantab, scannos=scannos, verify=verify )
	961	elif ( calmode == 'quotient' ):
	962	asaplog.push( 'Calibrating using quotient.' )
	963	print_log()
	964	scal = scantab.auto_quotient( verify=verify )
	965	elif ( calmode == 'ps' ):
	966	asaplog.push( 'Calibrating %s position-switched data.' % antname )
	967	print_log()
[1673]	968	if ( antname.find( 'APEX' ) != -1 ):
	969	scal = apexcal( scantab, scannos, calmode, verify )
	970	elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ):
[1633]	971	scal = almacal( scantab, scannos, calmode, verify )
	972	else:
	973	scal = calps( scantab, scannos=scannos, verify=verify )
	974	elif ( calmode == 'fs' or calmode == 'fsotf' ):
	975	asaplog.push( 'Calibrating %s frequency-switched data.' % antname )
	976	print_log()
[1673]	977	if ( antname.find( 'APEX' ) != -1 ):
	978	scal = apexcal( scantab, scannos, calmode, verify )
	979	elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ):
[1633]	980	scal = almacal( scantab, scannos, calmode, verify )
	981	else:
	982	scal = calfs( scantab, scannos=scannos, verify=verify )
	983	elif ( calmode == 'otf' ):
	984	asaplog.push( 'Calibrating %s On-The-Fly data.' % antname )
	985	print_log()
	986	scal = almacal( scantab, scannos, calmode, verify )
	987	else:
	988	asaplog.push( 'No calibration.' )
	989	scal = scantab.copy()
	990
	991	return scal
	992
[1673]	993	def apexcal( scantab, scannos=[], calmode='none', verify=False ):
[1633]	994	"""
[1673]	995	Calibrate APEX data
[1633]	996
	997	Parameters:
	998	scantab: scantable
	999	scannos: list of scan number
	1000	calmode: calibration mode
[1650]	1001
[1633]	1002	verify: verify calibration
	1003	"""
	1004	from asap._asap import stmath
	1005	stm = stmath()
	1006	antname = scantab.get_antennaname()
	1007	ssub = scantab.get_scan( scannos )
	1008	scal = scantable( stm.cwcal( ssub, calmode, antname ) )
	1009	return scal
[1650]	1010
[1673]	1011	def almacal( scantab, scannos=[], calmode='none', verify=False ):
	1012	"""
	1013	Calibrate ALMA data
	1014
	1015	Parameters:
	1016	scantab: scantable
	1017	scannos: list of scan number
	1018	calmode: calibration mode
	1019
	1020	verify: verify calibration
	1021	"""
	1022	from asap._asap import stmath
	1023	stm = stmath()
	1024	ssub = scantab.get_scan( scannos )
	1025	scal = scantable( stm.almacal( ssub, calmode ) )
	1026	return scal
	1027
[1650]	1028	def splitant(filename, outprefix='',overwrite=False):
	1029	"""
	1030	Split Measurement set by antenna name, save data as a scantables,
	1031	and return a list of filename.
	1032	Notice this method can only be available from CASA.
	1033	Prameter
	1034	filename: the name of Measurement set to be read.
	1035	outprefix: the prefix of output scantable name.
	1036	the names of output scantable will be
	1037	outprefix.antenna1, outprefix.antenna2, ....
[1651]	1038	If not specified, outprefix = filename is assumed.
[1650]	1039	overwrite If the file should be overwritten if it exists.
	1040	The default False is to return with warning
	1041	without writing the output. USE WITH CARE.
	1042
	1043	"""
	1044	# Import the table toolkit from CASA
	1045	try:
[1686]	1046	import casac
[1650]	1047	except ImportError:
[1686]	1048	if rcParams['verbose']:
	1049	#print "failed to load casa"
	1050	print_log()
	1051	asaplog.push("failed to load casa")
	1052	print_log('ERROR')
	1053	else: raise
	1054	return False
[1650]	1055	try:
[1686]	1056	tbtool = casac.homefinder.find_home_by_name('tableHome')
	1057	tb = tbtool.create()
	1058	tb2 = tbtool.create()
[1650]	1059	except:
[1686]	1060	if rcParams['verbose']:
	1061	#print "failed to load a table tool:\n", e
	1062	print_log()
	1063	asaplog.push("failed to load table tool")
	1064	print_log('ERROR')
	1065	else: raise
	1066	return False
[1650]	1067	# Check the input filename
	1068	if isinstance(filename, str):
	1069	import os.path
	1070	filename = os.path.expandvars(filename)
	1071	filename = os.path.expanduser(filename)
	1072	if not os.path.exists(filename):
	1073	s = "File '%s' not found." % (filename)
	1074	if rcParams['verbose']:
	1075	print_log()
	1076	asaplog.push(s)
	1077	print_log('ERROR')
	1078	return
	1079	raise IOError(s)
	1080	# check if input file is MS
	1081	if not os.path.isdir(filename) \
	1082	or not os.path.exists(filename+'/ANTENNA') \
	1083	or not os.path.exists(filename+'/table.f1'):
	1084	s = "File '%s' is not a Measurement set." % (filename)
	1085	if rcParams['verbose']:
	1086	print_log()
	1087	asaplog.push(s)
	1088	print_log('ERROR')
	1089	return
	1090	raise IOError(s)
	1091	else:
	1092	s = "The filename should be string. "
	1093	if rcParams['verbose']:
	1094	print_log()
	1095	asaplog.push(s)
	1096	print_log('ERROR')
	1097	return
	1098	raise TypeError(s)
	1099	# Check out put file name
	1100	outname=''
	1101	if len(outprefix) > 0: prefix=outprefix+'.'
[1651]	1102	else:
[1686]	1103	prefix=filename.rstrip('/')
[1650]	1104	# Now do the actual splitting.
	1105	outfiles=[]
	1106	tb.open(tablename=filename+'/ANTENNA',nomodify=True)
	1107	nant=tb.nrows()
	1108	antnames=tb.getcol('NAME',0,nant,1)
	1109	antpos=tb.getcol('POSITION',0,nant,1).transpose()
	1110	tb.close()
	1111	tb.open(tablename=filename,nomodify=True)
	1112	ant1=tb.getcol('ANTENNA1',0,-1,1)
[1686]	1113	tb.close()
[1650]	1114	for antid in set(ant1):
[1686]	1115	scan=scantable(filename,average=False,getpt=True,antenna=int(antid))
[1650]	1116	outname=prefix+antnames[antid]+'.asap'
	1117	scan.save(outname,format='ASAP',overwrite=overwrite)
[1686]	1118	del scan
[1650]	1119	outfiles.append(outname)
	1120	del tb, tb2
	1121	return outfiles
[1679]	1122
	1123	def _array2dOp( scan, value, mode="ADD", tsys=False ):
	1124	"""
	1125	This function is workaround on the basic operation of scantable
	1126	with 2 dimensional float list.
	1127
	1128	scan: scantable operand
	1129	value: float list operand
	1130	mode: operation mode (ADD, SUB, MUL, DIV)
	1131	tsys: if True, operate tsys as well
	1132	"""
	1133	nrow = scan.nrow()
	1134	s = None
[1682]	1135	if len( value ) == 1:
	1136	from asap._asap import stmath
	1137	stm = stmath()
	1138	s = scantable( stm._arrayop( scan.copy(), value[0], mode, tsys ) )
	1139	del stm
	1140	elif len( value ) != nrow:
	1141	asaplog.push( 'len(value) must be 1 or conform to scan.nrow()' )
[1679]	1142	print_log( 'ERROR' )
	1143	else:
	1144	from asap._asap import stmath
	1145	stm = stmath()
	1146	# insitu must be True
	1147	stm._setinsitu( True )
	1148	s = scan.copy()
	1149	sel = selector()
	1150	for irow in range( nrow ):
	1151	sel.set_rows( irow )
	1152	s.set_selection( sel )
	1153	if len( value[irow] ) == 1:
	1154	stm._unaryop( s, value[irow][0], mode, tsys )
	1155	else:
	1156	stm._arrayop( s, value[irow], mode, tsys, 'channel' )
	1157	s.set_selection()
	1158	sel.reset()
	1159	del sel
	1160	del stm
	1161	return s
	1162
	1163
	1164

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

Context Navigation

source: branches/alma/python/asapmath.py@ 1738

Download in other formats: