Context Navigation

asapmath.py@ 2924

Visit:

Last change on this file since 2924 was 2904, checked in by Takeshi Nakazato, 11 years ago

New Development: No

JIRA Issue: Yes CAS-5875

Ready for Test: Yes

Interface Changes: /No

What Interface Changed: Please list interface changes

Test Programs: test_sdcoadd

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Support for numeric value for freq_tol in sd.merge.

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

File size: 37.7 KB

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

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

Context Navigation

source: trunk/python/asapmath.py@ 2924

Download in other formats: