Context Navigation

asapmath.py@ 2141

Visit:

Last change on this file since 2141 was 2139, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: no

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: sdscale unit test

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Bug fix on array operation on scantable.

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

File size: 36.7 KB

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

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

Context Navigation

source: trunk/python/asapmath.py@ 2141

Download in other formats: