Context Navigation

asapmath.py@ 2413

Visit:

Last change on this file since 2413 was 1782, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: No (merge)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): single dish package

Description:

Merged some unimplemented developments from the trunk.
[Details]

Restored PKSmsg.cc, PKSmsg.h to alma branch from trunk
Updated pks_matchs.cc and src/SConscript to ones in trunk
uncommented "@print_log_dec"
implemented bug fixes to add(), sub(), mul(), and div() in scantable.py
implemented get_region() and get_point() in asaplotbase.py
activated interactive keyword for text(), arrow(), annotate(), axvline(), axhline(), axvspan(), and axhspan() in asapplotter.py

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

File size: 39.1 KB

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

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

Context Navigation

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

Download in other formats: