Context Navigation

asapmath.py@ 2563

Visit:

Last change on this file since 2563 was 2556, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

More speed-up of almacal: skip unecessary data copy.

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

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

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

Context Navigation

source: branches/hpc33/python/asapmath.py@ 2563

Download in other formats: