Context Navigation

asapmath.py@ 2798

Visit:

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

New Development: No

JIRA Issue: Yes CSV-2532 (may be related to CSV-1908 and CSV-2161)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: In tool level, added parameter 'freq_tolsr' to

scantable constructor and function sd.splitant.

Test Programs: test_sdsave, test_importasdm_sd

Put in Release Notes: Yes

Module(s): Module Names change impacts.

Description: Describe your changes here...

In importing MS to Scantable, frequency frame information is
imported as is by default, i.e., base frame in Scantable is
TOPO for ALMA data, which is forcibly converted to LSRK with
wrong time and direction reference.

Some functions have a boolean parameter 'freq_tolsr' that controls
the above behavior. If freq_tolsr is False (default), frequency
is imported as is, while frequency is converted to LSRK (wrongly)
when it is True.

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

File size: 37.8 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.asapplotter import new_asaplot
[2538]	6	from matplotlib import rc as rcp
[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.
[2472]	28	compel: True forces to average overwrapped IFs.
[113]	29	Example:
	30	# return a time averaged scan from scana and scanb
	31	# without using a mask
[129]	32	scanav = average_time(scana,scanb)
[1589]	33	# or equivalent
	34	# scanav = average_time([scana, scanb])
[113]	35	# return the (time) averaged scan, i.e. the average of
	36	# all correlator cycles
[558]	37	scanav = average_time(scan, scanav=True)
[101]	38	"""
[930]	39	scanav = False
[143]	40	if kwargs.has_key('scanav'):
[930]	41	scanav = kwargs.get('scanav')
[524]	42	weight = 'tint'
[143]	43	if kwargs.has_key('weight'):
	44	weight = kwargs.get('weight')
	45	mask = ()
	46	if kwargs.has_key('mask'):
	47	mask = kwargs.get('mask')
[930]	48	align = False
	49	if kwargs.has_key('align'):
	50	align = kwargs.get('align')
[1819]	51	compel = False
	52	if kwargs.has_key('compel'):
	53	compel = kwargs.get('compel')
[489]	54	varlist = vars()
[665]	55	if isinstance(args[0],list):
[981]	56	lst = args[0]
[665]	57	elif isinstance(args[0],tuple):
[981]	58	lst = list(args[0])
[665]	59	else:
[981]	60	lst = list(args)
[720]	61
[489]	62	del varlist["kwargs"]
	63	varlist["args"] = "%d scantables" % len(lst)
[981]	64	# need special formatting here for history...
[720]	65
[876]	66	from asap._asap import stmath
	67	stm = stmath()
[113]	68	for s in lst:
[101]	69	if not isinstance(s,scantable):
[720]	70	msg = "Please give a list of scantables"
[1859]	71	raise TypeError(msg)
[945]	72	if scanav: scanav = "SCAN"
	73	else: scanav = "NONE"
[981]	74	alignedlst = []
	75	if align:
	76	refepoch = lst[0].get_time(0)
	77	for scan in lst:
	78	alignedlst.append(scan.freq_align(refepoch,insitu=False))
	79	else:
[1059]	80	alignedlst = lst
[1232]	81	if weight.upper() == 'MEDIAN':
	82	# median doesn't support list of scantables - merge first
	83	merged = None
	84	if len(alignedlst) > 1:
	85	merged = merge(alignedlst)
	86	else:
	87	merged = alignedlst[0]
	88	s = scantable(stm._averagechannel(merged, 'MEDIAN', scanav))
	89	del merged
	90	else:
[1819]	91	#s = scantable(stm._average(alignedlst, mask, weight.upper(), scanav))
	92	s = scantable(stm._new_average(alignedlst, compel, mask, weight.upper(), scanav))
[489]	93	s._add_history("average_time",varlist)
[1859]	94
[489]	95	return s
[101]	96
[1862]	97	@asaplog_post_dec
[1074]	98	def quotient(source, reference, preserve=True):
	99	"""
	100	Return the quotient of a 'source' (signal) scan and a 'reference' scan.
	101	The reference can have just one scan, even if the signal has many. Otherwise
	102	they must have the same number of scans.
	103	The cursor of the output scan is set to 0
	104	Parameters:
	105	source: the 'on' scan
	106	reference: the 'off' scan
	107	preserve: you can preserve (default) the continuum or
	108	remove it. The equations used are
	109	preserve: Output = Toff * (on/off) - Toff
	110	remove: Output = Toff * (on/off) - Ton
	111	"""
	112	varlist = vars()
	113	from asap._asap import stmath
	114	stm = stmath()
	115	stm._setinsitu(False)
	116	s = scantable(stm._quotient(source, reference, preserve))
	117	s._add_history("quotient",varlist)
	118	return s
[101]	119
[1862]	120	@asaplog_post_dec
[1391]	121	def dototalpower(calon, caloff, tcalval=0.0):
	122	"""
	123	Do calibration for CAL on,off signals.
	124	Adopted from GBTIDL dototalpower
	125	Parameters:
	126	calon: the 'cal on' subintegration
	127	caloff: the 'cal off' subintegration
	128	tcalval: user supplied Tcal value
	129	"""
	130	varlist = vars()
	131	from asap._asap import stmath
	132	stm = stmath()
	133	stm._setinsitu(False)
	134	s = scantable(stm._dototalpower(calon, caloff, tcalval))
	135	s._add_history("dototalpower",varlist)
	136	return s
	137
[1862]	138	@asaplog_post_dec
[1391]	139	def dosigref(sig, ref, smooth, tsysval=0.0, tauval=0.0):
	140	"""
	141	Calculate a quotient (sig-ref/ref * Tsys)
	142	Adopted from GBTIDL dosigref
	143	Parameters:
	144	sig: on source data
	145	ref: reference data
	146	smooth: width of box car smoothing for reference
	147	tsysval: user specified Tsys (scalar only)
	148	tauval: user specified Tau (required if tsysval is set)
	149	"""
	150	varlist = vars()
	151	from asap._asap import stmath
	152	stm = stmath()
	153	stm._setinsitu(False)
	154	s = scantable(stm._dosigref(sig, ref, smooth, tsysval, tauval))
	155	s._add_history("dosigref",varlist)
	156	return s
	157
[1862]	158	@asaplog_post_dec
[1819]	159	def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
[1391]	160	"""
	161	Calibrate GBT position switched data
	162	Adopted from GBTIDL getps
[1920]	163	Currently calps identify the scans as position switched data if source
	164	type enum is pson or psoff. The data must contains 'CAL' signal
	165	on/off in each integration. To identify 'CAL' on state, the source type
[2472]	166	enum of poncal and poffcal need to be present.
[1391]	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')
[2150]	313	p=new_asaplot()
[2535]	314	rcp('lines', linewidth=1)
[1819]	315	#nr=min(6,len(ifnos)*len(polnos))
	316	nr=len(ifnos)*len(polnos)
	317	titles=[]
	318	btics=[]
	319	if nr<4:
	320	p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False)
	321	for i in range(2*nr):
	322	b=False
	323	if i >= 2*nr-2:
	324	b=True
	325	btics.append(b)
	326	elif nr==4:
	327	p.set_panels(rows=2,cols=4,nplots=8,ganged=False)
	328	for i in range(2*nr):
	329	b=False
	330	if i >= 2*nr-4:
	331	b=True
	332	btics.append(b)
	333	elif nr<7:
	334	p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
	335	for i in range(2*nr):
	336	if i >= 2*nr-4:
	337	b=True
	338	btics.append(b)
	339	else:
[1861]	340	asaplog.post()
[1819]	341	asaplog.push('Only first 6 [if,pol] pairs are plotted.')
[1861]	342	asaplog.post('WARN')
[1819]	343	nr=6
	344	for i in range(2*nr):
	345	b=False
	346	if i >= 2*nr-4:
	347	b=True
	348	btics.append(b)
	349	p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
	350	for i in range(nr):
	351	p.subplot(2*i)
	352	p.color=0
	353	title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
	354	titles.append(title)
	355	#p.set_axes('title',title,fontsize=40)
	356	ymin=1.0e100
	357	ymax=-1.0e100
[2348]	358	nchan=s.nchan(ifnos[int(i/len(polnos))])
[1819]	359	edge=int(nchan*0.01)
	360	for j in range(4):
	361	spmin=min(precal[keys[j]][i][edge:nchan-edge])
	362	spmax=max(precal[keys[j]][i][edge:nchan-edge])
	363	ymin=min(ymin,spmin)
	364	ymax=max(ymax,spmax)
	365	for j in range(4):
	366	if i==0:
	367	p.set_line(label=keys[j])
	368	else:
	369	p.legend()
	370	p.plot(precal[keys[j]][i])
	371	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	372	if not btics[2*i]:
	373	p.axes.set_xticks([])
	374	p.subplot(2*i+1)
	375	p.color=0
	376	title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
	377	titles.append(title)
	378	#p.set_axes('title',title)
	379	p.legend()
	380	ymin=postcal[i][edge:nchan-edge].min()
	381	ymax=postcal[i][edge:nchan-edge].max()
	382	p.plot(postcal[i])
	383	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	384	if not btics[2*i+1]:
	385	p.axes.set_xticks([])
	386	for i in range(2*nr):
	387	p.subplot(i)
	388	p.set_axes('title',titles[i],fontsize='medium')
	389	x=raw_input('Accept calibration ([y]/n): ' )
	390	if x.upper() == 'N':
[2150]	391	p.quit()
[1819]	392	del p
	393	return scabtab
[2150]	394	p.quit()
[1819]	395	del p
	396	###
[1391]	397	ress._add_history("calps", varlist)
	398	return ress
	399
[1862]	400	@asaplog_post_dec
[1819]	401	def calnod(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
[1391]	402	"""
	403	Do full (but a pair of scans at time) processing of GBT Nod data
	404	calibration.
	405	Adopted from GBTIDL's getnod
	406	Parameters:
	407	scantab: scantable
	408	scannos: a pair of scan numbers, or the first scan number of the pair
	409	smooth: box car smoothing order
	410	tsysval: optional user specified Tsys value
	411	tauval: optional user specified tau value (not implemented yet)
	412	tcalval: optional user specified Tcal value
[1920]	413	verify: Verify calibration if true
[1391]	414	"""
	415	varlist = vars()
	416	from asap._asap import stmath
[1819]	417	from asap._asap import srctype
[1391]	418	stm = stmath()
	419	stm._setinsitu(False)
	420
	421	# check for the appropriate data
[1819]	422	## s = scantab.get_scan('_nod')
	423	## if s is None:
	424	## msg = "The input data appear to contain no Nod observing mode data."
[1859]	425	## raise TypeError(msg)
[1819]	426	s = scantab.copy()
	427	sel = selector()
	428	sel.set_types( srctype.nod )
	429	try:
	430	s.set_selection( sel )
	431	except Exception, e:
[1391]	432	msg = "The input data appear to contain no Nod observing mode data."
[1859]	433	raise TypeError(msg)
[1819]	434	sel.reset()
	435	del sel
	436	del s
[1391]	437
	438	# need check correspondance of each beam with sig-ref ...
	439	# check for timestamps, scan numbers, subscan id (not available in
	440	# ASAP data format...). Assume 1st scan of the pair (beam 0 - sig
	441	# and beam 1 - ref...)
	442	# First scan number of paired scans or list of pairs of
	443	# scan numbers (has to have even number of pairs.)
	444
	445	#data splitting
	446	scan1no = scan2no = 0
	447
	448	if len(scannos)==1:
	449	scan1no = scannos[0]
	450	scan2no = scannos[0]+1
	451	pairScans = [scan1no, scan2no]
	452	else:
	453	#if len(scannos)>2:
	454	# msg = "calnod can only process a pair of nod scans at time."
[1859]	455	# raise TypeError(msg)
[1391]	456	#
	457	#if len(scannos)==2:
	458	# scan1no = scannos[0]
	459	# scan2no = scannos[1]
	460	pairScans = list(scannos)
	461
	462	if tsysval>0.0:
	463	if tauval<=0.0:
	464	msg = "Need to supply a valid tau to use the supplied Tsys"
[1859]	465	raise TypeError(msg)
[1391]	466	else:
	467	scantab.recalc_azel()
	468	resspec = scantable(stm._donod(scantab, pairScans, smooth, tsysval,tauval,tcalval))
[1819]	469	###
	470	if verify:
	471	# get data
	472	import numpy
	473	precal={}
	474	postcal=[]
	475	keys=['','_calon']
	476	types=[srctype.nod,srctype.nodcal]
	477	ifnos=list(scantab.getifnos())
	478	polnos=list(scantab.getpolnos())
	479	sel=selector()
	480	ss = scantab.copy()
	481	for i in range(2):
	482	#ss=scantab.get_scan('*'+keys[i])
	483	ll=[]
	484	ll2=[]
	485	for j in range(len(ifnos)):
	486	for k in range(len(polnos)):
	487	sel.set_ifs(ifnos[j])
	488	sel.set_polarizations(polnos[k])
	489	sel.set_scans(pairScans[0])
	490	sel.set_types(types[i])
	491	try:
	492	ss.set_selection(sel)
	493	except:
	494	continue
	495	ll.append(numpy.array(ss._getspectrum(0)))
	496	sel.reset()
	497	ss.set_selection()
	498	sel.set_ifs(ifnos[j])
	499	sel.set_polarizations(polnos[k])
	500	sel.set_scans(pairScans[1])
	501	sel.set_types(types[i])
	502	try:
	503	ss.set_selection(sel)
	504	except:
	505	ll.pop()
	506	continue
	507	ll2.append(numpy.array(ss._getspectrum(0)))
	508	sel.reset()
	509	ss.set_selection()
	510	key='%s%s' %(pairScans[0],keys[i])
	511	precal[key]=ll
	512	key='%s%s' %(pairScans[1],keys[i])
	513	precal[key]=ll2
	514	#del ss
	515	keys=precal.keys()
	516	for j in range(len(ifnos)):
	517	for k in range(len(polnos)):
	518	sel.set_ifs(ifnos[j])
	519	sel.set_polarizations(polnos[k])
	520	sel.set_scans(pairScans[0])
	521	try:
	522	resspec.set_selection(sel)
	523	except:
	524	continue
	525	postcal.append(numpy.array(resspec._getspectrum(0)))
	526	sel.reset()
	527	resspec.set_selection()
	528	del sel
	529	# plot
[1861]	530	asaplog.post()
[1819]	531	asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
[1861]	532	asaplog.post('WARN')
[2150]	533	p=new_asaplot()
[2535]	534	rcp('lines', linewidth=1)
[1819]	535	#nr=min(6,len(ifnos)*len(polnos))
	536	nr=len(ifnos)*len(polnos)
	537	titles=[]
	538	btics=[]
	539	if nr<4:
	540	p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False)
	541	for i in range(2*nr):
	542	b=False
	543	if i >= 2*nr-2:
	544	b=True
	545	btics.append(b)
	546	elif nr==4:
	547	p.set_panels(rows=2,cols=4,nplots=8,ganged=False)
	548	for i in range(2*nr):
	549	b=False
	550	if i >= 2*nr-4:
	551	b=True
	552	btics.append(b)
	553	elif nr<7:
	554	p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
	555	for i in range(2*nr):
	556	if i >= 2*nr-4:
	557	b=True
	558	btics.append(b)
	559	else:
[1861]	560	asaplog.post()
[1819]	561	asaplog.push('Only first 6 [if,pol] pairs are plotted.')
[1861]	562	asaplog.post('WARN')
[1819]	563	nr=6
	564	for i in range(2*nr):
	565	b=False
	566	if i >= 2*nr-4:
	567	b=True
	568	btics.append(b)
	569	p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
	570	for i in range(nr):
	571	p.subplot(2*i)
	572	p.color=0
	573	title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
	574	titles.append(title)
	575	#p.set_axes('title',title,fontsize=40)
	576	ymin=1.0e100
	577	ymax=-1.0e100
[2348]	578	nchan=scantab.nchan(ifnos[int(i/len(polnos))])
[1819]	579	edge=int(nchan*0.01)
	580	for j in range(4):
	581	spmin=min(precal[keys[j]][i][edge:nchan-edge])
	582	spmax=max(precal[keys[j]][i][edge:nchan-edge])
	583	ymin=min(ymin,spmin)
	584	ymax=max(ymax,spmax)
	585	for j in range(4):
	586	if i==0:
	587	p.set_line(label=keys[j])
	588	else:
	589	p.legend()
	590	p.plot(precal[keys[j]][i])
	591	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	592	if not btics[2*i]:
	593	p.axes.set_xticks([])
	594	p.subplot(2*i+1)
	595	p.color=0
	596	title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
	597	titles.append(title)
	598	#p.set_axes('title',title)
	599	p.legend()
	600	ymin=postcal[i][edge:nchan-edge].min()
	601	ymax=postcal[i][edge:nchan-edge].max()
	602	p.plot(postcal[i])
	603	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	604	if not btics[2*i+1]:
	605	p.axes.set_xticks([])
	606	for i in range(2*nr):
	607	p.subplot(i)
	608	p.set_axes('title',titles[i],fontsize='medium')
	609	x=raw_input('Accept calibration ([y]/n): ' )
	610	if x.upper() == 'N':
[2150]	611	p.quit()
[1819]	612	del p
	613	return scabtab
[2150]	614	p.quit()
[1819]	615	del p
	616	###
[1391]	617	resspec._add_history("calnod",varlist)
	618	return resspec
	619
[1862]	620	@asaplog_post_dec
[1819]	621	def calfs(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
[1391]	622	"""
	623	Calibrate GBT frequency switched data.
	624	Adopted from GBTIDL getfs.
[1920]	625	Currently calfs identify the scans as frequency switched data if source
	626	type enum is fson and fsoff. The data must contains 'CAL' signal
	627	on/off in each integration. To identify 'CAL' on state, the source type
[2472]	628	enum of foncal and foffcal need to be present.
[1391]	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')
[2150]	744	p=new_asaplot()
[2535]	745	rcp('lines', linewidth=1)
[1819]	746	#nr=min(6,len(ifnos)*len(polnos))
	747	nr=len(ifnos)/2*len(polnos)
	748	titles=[]
	749	btics=[]
	750	if nr>3:
[1861]	751	asaplog.post()
[1819]	752	asaplog.push('Only first 3 [if,pol] pairs are plotted.')
[1861]	753	asaplog.post('WARN')
[1819]	754	nr=3
	755	p.set_panels(rows=nr,cols=3,nplots=3*nr,ganged=False)
	756	for i in range(3*nr):
	757	b=False
	758	if i >= 3*nr-3:
	759	b=True
	760	btics.append(b)
	761	for i in range(nr):
	762	p.subplot(3*i)
	763	p.color=0
	764	title='raw data IF%s,%s POL%s' % (ifnos[2int(i/len(polnos))],ifnos[2int(i/len(polnos))+1],polnos[i%len(polnos)])
	765	titles.append(title)
	766	#p.set_axes('title',title,fontsize=40)
	767	ymin=1.0e100
	768	ymax=-1.0e100
[2348]	769	nchan=s.nchan(ifnos[2*int(i/len(polnos))])
[1819]	770	edge=int(nchan*0.01)
	771	for j in range(4):
	772	spmin=min(precal[keys[j]][i][edge:nchan-edge])
	773	spmax=max(precal[keys[j]][i][edge:nchan-edge])
	774	ymin=min(ymin,spmin)
	775	ymax=max(ymax,spmax)
	776	for j in range(4):
	777	if i==0:
	778	p.set_line(label=keys[j])
	779	else:
	780	p.legend()
	781	p.plot(precal[keys[j]][i])
	782	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	783	if not btics[3*i]:
	784	p.axes.set_xticks([])
	785	p.subplot(3*i+1)
	786	p.color=0
	787	title='sig data IF%s POL%s' % (ifnos[2*int(i/len(polnos))],polnos[i%len(polnos)])
	788	titles.append(title)
	789	#p.set_axes('title',title)
	790	p.legend()
	791	ymin=postcal[2*i][edge:nchan-edge].min()
	792	ymax=postcal[2*i][edge:nchan-edge].max()
	793	p.plot(postcal[2*i])
	794	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	795	if not btics[3*i+1]:
	796	p.axes.set_xticks([])
	797	p.subplot(3*i+2)
	798	p.color=0
	799	title='ref data IF%s POL%s' % (ifnos[2*int(i/len(polnos))+1],polnos[i%len(polnos)])
	800	titles.append(title)
	801	#p.set_axes('title',title)
	802	p.legend()
	803	ymin=postcal[2*i+1][edge:nchan-edge].min()
	804	ymax=postcal[2*i+1][edge:nchan-edge].max()
	805	p.plot(postcal[2*i+1])
	806	p.axes.set_ylim(ymin-0.1abs(ymin),ymax+0.1abs(ymax))
	807	if not btics[3*i+2]:
	808	p.axes.set_xticks([])
	809	for i in range(3*nr):
	810	p.subplot(i)
	811	p.set_axes('title',titles[i],fontsize='medium')
	812	x=raw_input('Accept calibration ([y]/n): ' )
	813	if x.upper() == 'N':
[2150]	814	p.quit()
[1819]	815	del p
	816	return scabtab
[2150]	817	p.quit()
[1819]	818	del p
	819	###
[1391]	820	resspec._add_history("calfs",varlist)
	821	return resspec
	822
[1862]	823	@asaplog_post_dec
[918]	824	def merge(*args):
[945]	825	"""
[1362]	826	Merge a list of scanatables, or comma-sperated scantables into one
	827	scnatble.
	828	Parameters:
	829	A list [scan1, scan2] or scan1, scan2.
	830	Example:
	831	myscans = [scan1, scan2]
[1589]	832	allscans = merge(myscans)
	833	# or equivalent
	834	sameallscans = merge(scan1, scan2)
[945]	835	"""
[918]	836	varlist = vars()
	837	if isinstance(args[0],list):
	838	lst = tuple(args[0])
	839	elif isinstance(args[0],tuple):
	840	lst = args[0]
	841	else:
	842	lst = tuple(args)
	843	varlist["args"] = "%d scantables" % len(lst)
	844	# need special formatting her for history...
	845	from asap._asap import stmath
	846	stm = stmath()
	847	for s in lst:
	848	if not isinstance(s,scantable):
	849	msg = "Please give a list of scantables"
[1859]	850	raise TypeError(msg)
[918]	851	s = scantable(stm._merge(lst))
	852	s._add_history("merge", varlist)
	853	return s
[1819]	854
[1862]	855	@asaplog_post_dec
[1819]	856	def calibrate( scantab, scannos=[], calmode='none', verify=None ):
	857	"""
	858	Calibrate data.
[1826]	859
[1819]	860	Parameters:
	861	scantab: scantable
	862	scannos: list of scan number
	863	calmode: calibration mode
[1826]	864	verify: verify calibration
[1819]	865	"""
[2102]	866	import re
[1819]	867	antname = scantab.get_antennaname()
	868	if ( calmode == 'nod' ):
	869	asaplog.push( 'Calibrating nod data.' )
	870	scal = calnod( scantab, scannos=scannos, verify=verify )
	871	elif ( calmode == 'quotient' ):
	872	asaplog.push( 'Calibrating using quotient.' )
	873	scal = scantab.auto_quotient( verify=verify )
	874	elif ( calmode == 'ps' ):
	875	asaplog.push( 'Calibrating %s position-switched data.' % antname )
	876	if ( antname.find( 'APEX' ) != -1 ):
	877	scal = apexcal( scantab, scannos, calmode, verify )
[2102]	878	elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1
	879	or re.match('DV[0-9][0-9]$',antname) is not None
	880	or re.match('PM[0-9][0-9]$',antname) is not None
	881	or re.match('CM[0-9][0-9]$',antname) is not None
	882	or re.match('DA[0-9][0-9]$',antname) is not None ):
[1819]	883	scal = almacal( scantab, scannos, calmode, verify )
	884	else:
	885	scal = calps( scantab, scannos=scannos, verify=verify )
	886	elif ( calmode == 'fs' or calmode == 'fsotf' ):
	887	asaplog.push( 'Calibrating %s frequency-switched data.' % antname )
	888	if ( antname.find( 'APEX' ) != -1 ):
	889	scal = apexcal( scantab, scannos, calmode, verify )
	890	elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ):
	891	scal = almacal( scantab, scannos, calmode, verify )
	892	else:
	893	scal = calfs( scantab, scannos=scannos, verify=verify )
	894	elif ( calmode == 'otf' ):
	895	asaplog.push( 'Calibrating %s On-The-Fly data.' % antname )
	896	scal = almacal( scantab, scannos, calmode, verify )
	897	else:
	898	asaplog.push( 'No calibration.' )
	899	scal = scantab.copy()
	900
[1826]	901	return scal
[1819]	902
	903	def apexcal( scantab, scannos=[], calmode='none', verify=False ):
	904	"""
	905	Calibrate APEX data
	906
	907	Parameters:
	908	scantab: scantable
	909	scannos: list of scan number
	910	calmode: calibration mode
	911
[1826]	912	verify: verify calibration
[1819]	913	"""
	914	from asap._asap import stmath
	915	stm = stmath()
	916	antname = scantab.get_antennaname()
[2580]	917	selection=selector()
	918	selection.set_scans(scannos)
	919	orig = scantab.get_selection()
	920	scantab.set_selection(orig+selection)
	921	## ssub = scantab.get_scan( scannos )
	922	## scal = scantable( stm.cwcal( ssub, calmode, antname ) )
	923	scal = scantable( stm.cwcal( scantab, calmode, antname ) )
	924	scantab.set_selection(orig)
[1819]	925	return scal
	926
	927	def almacal( scantab, scannos=[], calmode='none', verify=False ):
	928	"""
	929	Calibrate ALMA data
	930
	931	Parameters:
	932	scantab: scantable
	933	scannos: list of scan number
	934	calmode: calibration mode
	935
[1826]	936	verify: verify calibration
[1819]	937	"""
	938	from asap._asap import stmath
	939	stm = stmath()
[2580]	940	selection=selector()
	941	selection.set_scans(scannos)
	942	orig = scantab.get_selection()
	943	scantab.set_selection(orig+selection)
	944	## ssub = scantab.get_scan( scannos )
	945	## scal = scantable( stm.almacal( ssub, calmode ) )
	946	scal = scantable( stm.almacal( scantab, calmode ) )
	947	scantab.set_selection(orig)
[1819]	948	return scal
	949
[1862]	950	@asaplog_post_dec
[2754]	951	def splitant(filename, outprefix='',overwrite=False, getpt=True, freq_tolsr=False):
[1819]	952	"""
	953	Split Measurement set by antenna name, save data as a scantables,
	954	and return a list of filename.
[1826]	955	Notice this method can only be available from CASA.
[1819]	956	Prameter
[1826]	957	filename: the name of Measurement set to be read.
[1819]	958	outprefix: the prefix of output scantable name.
	959	the names of output scantable will be
	960	outprefix.antenna1, outprefix.antenna2, ....
	961	If not specified, outprefix = filename is assumed.
	962	overwrite If the file should be overwritten if it exists.
	963	The default False is to return with warning
	964	without writing the output. USE WITH CARE.
[2754]	965	getpt Whether to import direction from MS/POINTING
	966	table or not. Default is True (import direction).
	967	freq_tolsr Whether to convert frequency frame information
	968	to LSRK or not. Default is False (import as is).
[1826]	969
[1819]	970	"""
	971	# Import the table toolkit from CASA
[2646]	972	from casac import casac
[1918]	973	from asap.scantable import is_ms
[2646]	974	tb = casac.table()
[1819]	975	# Check the input filename
	976	if isinstance(filename, str):
	977	import os.path
	978	filename = os.path.expandvars(filename)
	979	filename = os.path.expanduser(filename)
	980	if not os.path.exists(filename):
	981	s = "File '%s' not found." % (filename)
	982	raise IOError(s)
	983	# check if input file is MS
[1883]	984	#if not os.path.isdir(filename) \
	985	# or not os.path.exists(filename+'/ANTENNA') \
	986	# or not os.path.exists(filename+'/table.f1'):
	987	if not is_ms(filename):
[1819]	988	s = "File '%s' is not a Measurement set." % (filename)
	989	raise IOError(s)
	990	else:
	991	s = "The filename should be string. "
	992	raise TypeError(s)
	993	# Check out put file name
	994	outname=''
	995	if len(outprefix) > 0: prefix=outprefix+'.'
	996	else:
	997	prefix=filename.rstrip('/')
	998	# Now do the actual splitting.
	999	outfiles=[]
[2034]	1000	tb.open(tablename=filename,nomodify=True)
	1001	ant1=tb.getcol('ANTENNA1',0,-1,1)
[2366]	1002	#anttab=tb.getkeyword('ANTENNA').split()[-1]
	1003	anttab=tb.getkeyword('ANTENNA').lstrip('Table: ')
[2034]	1004	tb.close()
	1005	#tb.open(tablename=filename+'/ANTENNA',nomodify=True)
	1006	tb.open(tablename=anttab,nomodify=True)
[1819]	1007	nant=tb.nrows()
	1008	antnames=tb.getcol('NAME',0,nant,1)
	1009	tb.close()
[1880]	1010	tmpname='asapmath.splitant.tmp'
[1819]	1011	for antid in set(ant1):
[1883]	1012	tb.open(tablename=filename,nomodify=True)
	1013	tbsel=tb.query('ANTENNA1 == %s && ANTENNA2 == %s'%(antid,antid),tmpname)
[2754]	1014	scan=scantable(tmpname,average=False,antenna=int(antid),getpt=getpt,freq_tolsr=freq_tolsr,)
[1918]	1015	outname=prefix+antnames[antid]+'.asap'
	1016	scan.save(outname,format='ASAP',overwrite=overwrite)
[1880]	1017	tbsel.close()
[1883]	1018	tb.close()
[1880]	1019	del tbsel
[1819]	1020	del scan
	1021	outfiles.append(outname)
[1880]	1022	os.system('rm -rf '+tmpname)
	1023	del tb
[1819]	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@ 2798

Download in other formats: