schain Commit - r502:8b1242a708e8 · Jicamarca Repository

Processing Modules added:...

Julio Valdez -

r502:8b1242a708e8

parent child

Context file:

r502:8b1242a708e8

Collapse all files

schainpy/model/graphics/jroplot_correlation.py created 644 +196 0

@@ -0,0 +1,196
	1	import os
	2	import datetime
	3	import numpy
	4	import copy
	5
	6	from figure import Figure, isRealtime
	7
	8	class CorrelationPlot(Figure):
	9
	10	isConfig = None
	11	__nsubplots = None
	12
	13	WIDTHPROF = None
	14	HEIGHTPROF = None
	15	PREFIX = 'corr'
	16
	17	def __init__(self):
	18
	19	self.isConfig = False
	20	self.__nsubplots = 1
	21
	22	self.WIDTH = 280
	23	self.HEIGHT = 250
	24	self.WIDTHPROF = 120
	25	self.HEIGHTPROF = 0
	26	self.counter_imagwr = 0
	27
	28	self.PLOT_CODE = 1
	29	self.FTP_WEI = None
	30	self.EXP_CODE = None
	31	self.SUB_EXP_CODE = None
	32	self.PLOT_POS = None
	33
	34	def getSubplots(self):
	35
	36	ncol = int(numpy.sqrt(self.nplots)+0.9)
	37	nrow = int(self.nplots*1./ncol + 0.9)
	38
	39	return nrow, ncol
	40
	41	def setup(self, id, nplots, wintitle, showprofile=False, show=True):
	42
	43	showprofile = False
	44	self.__showprofile = showprofile
	45	self.nplots = nplots
	46
	47	ncolspan = 1
	48	colspan = 1
	49	if showprofile:
	50	ncolspan = 3
	51	colspan = 2
	52	self.__nsubplots = 2
	53
	54	self.createFigure(id = id,
	55	wintitle = wintitle,
	56	widthplot = self.WIDTH + self.WIDTHPROF,
	57	heightplot = self.HEIGHT + self.HEIGHTPROF,
	58	show=show)
	59
	60	nrow, ncol = self.getSubplots()
	61
	62	counter = 0
	63	for y in range(nrow):
	64	for x in range(ncol):
	65
	66	if counter >= self.nplots:
	67	break
	68
	69	self.addAxes(nrow, ncolncolspan, y, xncolspan, colspan, 1)
	70
	71	if showprofile:
	72	self.addAxes(nrow, ncolncolspan, y, xncolspan+colspan, 1, 1)
	73
	74	counter += 1
	75
	76	def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
	77	xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
	78	save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
	79	server=None, folder=None, username=None, password=None,
	80	ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
	81
	82	"""
	83
	84	Input:
	85	dataOut :
	86	id :
	87	wintitle :
	88	channelList :
	89	showProfile :
	90	xmin : None,
	91	xmax : None,
	92	ymin : None,
	93	ymax : None,
	94	zmin : None,
	95	zmax : None
	96	"""
	97
	98	if dataOut.flagNoData:
	99	return None
	100
	101	if realtime:
	102	if not(isRealtime(utcdatatime = dataOut.utctime)):
	103	print 'Skipping this plot function'
	104	return
	105
	106	if channelList == None:
	107	channelIndexList = dataOut.channelIndexList
	108	else:
	109	channelIndexList = []
	110	for channel in channelList:
	111	if channel not in dataOut.channelList:
	112	raise ValueError, "Channel %d is not in dataOut.channelList"
	113	channelIndexList.append(dataOut.channelList.index(channel))
	114
	115	factor = dataOut.normFactor
	116	lenfactor = factor.shape[1]
	117	x = dataOut.getLagTRange(1)
	118	y = dataOut.getHeiRange()
	119
	120	z = copy.copy(dataOut.data_corr[:,:,0,:])
	121	for i in range(dataOut.data_corr.shape[0]):
	122	z[i,:,:] = z[i,:,:]/factor[i,:]
	123	zdB = numpy.abs(z)
	124
	125	avg = numpy.average(z, axis=1)
	126	# avg = numpy.nanmean(z, axis=1)
	127	# noise = dataOut.noise/factor
	128
	129	#thisDatetime = dataOut.datatime
	130	thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
	131	title = wintitle + " Correlation"
	132	xlabel = "Lag T (s)"
	133	ylabel = "Range (Km)"
	134
	135	if not self.isConfig:
	136
	137	nplots = dataOut.data_corr.shape[0]
	138
	139	self.setup(id=id,
	140	nplots=nplots,
	141	wintitle=wintitle,
	142	showprofile=showprofile,
	143	show=show)
	144
	145	if xmin == None: xmin = numpy.nanmin(x)
	146	if xmax == None: xmax = numpy.nanmax(x)
	147	if ymin == None: ymin = numpy.nanmin(y)
	148	if ymax == None: ymax = numpy.nanmax(y)
	149	if zmin == None: zmin = 0
	150	if zmax == None: zmax = 1
	151
	152	self.FTP_WEI = ftp_wei
	153	self.EXP_CODE = exp_code
	154	self.SUB_EXP_CODE = sub_exp_code
	155	self.PLOT_POS = plot_pos
	156
	157	self.isConfig = True
	158
	159	self.setWinTitle(title)
	160
	161	for i in range(self.nplots):
	162	str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
	163	title = "Channel %d and %d: : %s" %(dataOut.pairsList[i][0]+1,dataOut.pairsList[i][1]+1 , str_datetime)
	164	axes = self.axesList[i*self.__nsubplots]
	165	axes.pcolor(x, y, zdB[i,:,:],
	166	xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
	167	xlabel=xlabel, ylabel=ylabel, title=title,
	168	ticksize=9, cblabel='')
	169
	170	# if self.__showprofile:
	171	# axes = self.axesList[i*self.__nsubplots +1]
	172	# axes.pline(avgdB[i], y,
	173	# xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
	174	# xlabel='dB', ylabel='', title='',
	175	# ytick_visible=False,
	176	# grid='x')
	177	#
	178	# noiseline = numpy.repeat(noisedB[i], len(y))
	179	# axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
	180
	181	self.draw()
	182
	183	if figfile == None:
	184	str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
	185	figfile = self.getFilename(name = str_datetime)
	186
	187	if figpath != '':
	188	self.counter_imagwr += 1
	189	if (self.counter_imagwr>=wr_period):
	190	# store png plot to local folder
	191	self.saveFigure(figpath, figfile)
	192	# store png plot to FTP server according to RT-Web format
	193	name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
	194	ftp_filename = os.path.join(figpath, name)
	195	self.saveFigure(figpath, ftp_filename)
	196	self.counter_imagwr = 0

schainpy/model/graphics/jroplot_parameters.py created 644 +579 0

This diff has been collapsed as it changes many lines, (579 lines changed) Show them Hide them
	@@ -0,0 +1,579
		1	import os
		2	import datetime
		3	import numpy
		4
		5	from figure import Figure, isRealtime
		6
		7	class MomentsPlot(Figure):
		8
		9	isConfig = None
		10	__nsubplots = None
		11
		12	WIDTHPROF = None
		13	HEIGHTPROF = None
		14	PREFIX = 'prm'
		15
		16	def __init__(self):
		17
		18	self.isConfig = False
		19	self.__nsubplots = 1
		20
		21	self.WIDTH = 280
		22	self.HEIGHT = 250
		23	self.WIDTHPROF = 120
		24	self.HEIGHTPROF = 0
		25	self.counter_imagwr = 0
		26
		27	self.PLOT_CODE = 1
		28	self.FTP_WEI = None
		29	self.EXP_CODE = None
		30	self.SUB_EXP_CODE = None
		31	self.PLOT_POS = None
		32
		33	def getSubplots(self):
		34
		35	ncol = int(numpy.sqrt(self.nplots)+0.9)
		36	nrow = int(self.nplots*1./ncol + 0.9)
		37
		38	return nrow, ncol
		39
		40	def setup(self, id, nplots, wintitle, showprofile=True, show=True):
		41
		42	self.__showprofile = showprofile
		43	self.nplots = nplots
		44
		45	ncolspan = 1
		46	colspan = 1
		47	if showprofile:
		48	ncolspan = 3
		49	colspan = 2
		50	self.__nsubplots = 2
		51
		52	self.createFigure(id = id,
		53	wintitle = wintitle,
		54	widthplot = self.WIDTH + self.WIDTHPROF,
		55	heightplot = self.HEIGHT + self.HEIGHTPROF,
		56	show=show)
		57
		58	nrow, ncol = self.getSubplots()
		59
		60	counter = 0
		61	for y in range(nrow):
		62	for x in range(ncol):
		63
		64	if counter >= self.nplots:
		65	break
		66
		67	self.addAxes(nrow, ncolncolspan, y, xncolspan, colspan, 1)
		68
		69	if showprofile:
		70	self.addAxes(nrow, ncolncolspan, y, xncolspan+colspan, 1, 1)
		71
		72	counter += 1
		73
		74	def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
		75	xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
		76	save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
		77	server=None, folder=None, username=None, password=None,
		78	ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
		79
		80	"""
		81
		82	Input:
		83	dataOut :
		84	id :
		85	wintitle :
		86	channelList :
		87	showProfile :
		88	xmin : None,
		89	xmax : None,
		90	ymin : None,
		91	ymax : None,
		92	zmin : None,
		93	zmax : None
		94	"""
		95
		96	if dataOut.flagNoData:
		97	return None
		98
		99	if realtime:
		100	if not(isRealtime(utcdatatime = dataOut.utctime)):
		101	print 'Skipping this plot function'
		102	return
		103
		104	if channelList == None:
		105	channelIndexList = dataOut.channelIndexList
		106	else:
		107	channelIndexList = []
		108	for channel in channelList:
		109	if channel not in dataOut.channelList:
		110	raise ValueError, "Channel %d is not in dataOut.channelList"
		111	channelIndexList.append(dataOut.channelList.index(channel))
		112
		113	factor = dataOut.normFactor
		114	x = dataOut.abscissaRange
		115	y = dataOut.heightRange
		116
		117	z = dataOut.data_pre[channelIndexList,:,:]/factor
		118	z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
		119	avg = numpy.average(z, axis=1)
		120	noise = dataOut.noise/factor
		121
		122	zdB = 10*numpy.log10(z)
		123	avgdB = 10*numpy.log10(avg)
		124	noisedB = 10*numpy.log10(noise)
		125
		126	#thisDatetime = dataOut.datatime
		127	thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
		128	title = wintitle + " Parameters"
		129	xlabel = "Velocity (m/s)"
		130	ylabel = "Range (Km)"
		131
		132	if not self.isConfig:
		133
		134	nplots = len(channelIndexList)
		135
		136	self.setup(id=id,
		137	nplots=nplots,
		138	wintitle=wintitle,
		139	showprofile=showprofile,
		140	show=show)
		141
		142	if xmin == None: xmin = numpy.nanmin(x)
		143	if xmax == None: xmax = numpy.nanmax(x)
		144	if ymin == None: ymin = numpy.nanmin(y)
		145	if ymax == None: ymax = numpy.nanmax(y)
		146	if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
		147	if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
		148
		149	self.FTP_WEI = ftp_wei
		150	self.EXP_CODE = exp_code
		151	self.SUB_EXP_CODE = sub_exp_code
		152	self.PLOT_POS = plot_pos
		153
		154	self.isConfig = True
		155
		156	self.setWinTitle(title)
		157
		158	for i in range(self.nplots):
		159	str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
		160	title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i]+1, noisedB[i], str_datetime)
		161	axes = self.axesList[i*self.__nsubplots]
		162	axes.pcolor(x, y, zdB[i,:,:],
		163	xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
		164	xlabel=xlabel, ylabel=ylabel, title=title,
		165	ticksize=9, cblabel='')
		166	#Mean Line
		167	mean = dataOut.data_param[i, 1, :]
		168	axes.addpline(mean, y, idline=0, color="black", linestyle="solid", lw=1)
		169
		170	if self.__showprofile:
		171	axes = self.axesList[i*self.__nsubplots +1]
		172	axes.pline(avgdB[i], y,
		173	xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
		174	xlabel='dB', ylabel='', title='',
		175	ytick_visible=False,
		176	grid='x')
		177
		178	noiseline = numpy.repeat(noisedB[i], len(y))
		179	axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
		180
		181	self.draw()
		182
		183	if figfile == None:
		184	str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
		185	figfile = self.getFilename(name = str_datetime)
		186
		187	if figpath != '':
		188	self.counter_imagwr += 1
		189	if (self.counter_imagwr>=wr_period):
		190	# store png plot to local folder
		191	self.saveFigure(figpath, figfile)
		192	# store png plot to FTP server according to RT-Web format
		193	name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
		194	ftp_filename = os.path.join(figpath, name)
		195	self.saveFigure(figpath, ftp_filename)
		196	self.counter_imagwr = 0
		197
		198	class SkyMapPlot(Figure):
		199
		200	__isConfig = None
		201	__nsubplots = None
		202
		203	WIDTHPROF = None
		204	HEIGHTPROF = None
		205	PREFIX = 'prm'
		206
		207	def __init__(self):
		208
		209	self.__isConfig = False
		210	self.__nsubplots = 1
		211
		212	# self.WIDTH = 280
		213	# self.HEIGHT = 250
		214	self.WIDTH = 600
		215	self.HEIGHT = 600
		216	self.WIDTHPROF = 120
		217	self.HEIGHTPROF = 0
		218	self.counter_imagwr = 0
		219
		220	self.PLOT_CODE = 1
		221	self.FTP_WEI = None
		222	self.EXP_CODE = None
		223	self.SUB_EXP_CODE = None
		224	self.PLOT_POS = None
		225
		226	def getSubplots(self):
		227
		228	ncol = int(numpy.sqrt(self.nplots)+0.9)
		229	nrow = int(self.nplots*1./ncol + 0.9)
		230
		231	return nrow, ncol
		232
		233	def setup(self, id, nplots, wintitle, showprofile=False, show=True):
		234
		235	self.__showprofile = showprofile
		236	self.nplots = nplots
		237
		238	ncolspan = 1
		239	colspan = 1
		240
		241	self.createFigure(id = id,
		242	wintitle = wintitle,
		243	widthplot = self.WIDTH, #+ self.WIDTHPROF,
		244	heightplot = self.HEIGHT,# + self.HEIGHTPROF,
		245	show=show)
		246
		247	nrow, ncol = 1,1
		248	counter = 0
		249	x = 0
		250	y = 0
		251	self.addAxes(1, 1, 0, 0, 1, 1, True)
		252
		253	def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
		254	xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
		255	save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
		256	server=None, folder=None, username=None, password=None,
		257	ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
		258
		259	"""
		260
		261	Input:
		262	dataOut :
		263	id :
		264	wintitle :
		265	channelList :
		266	showProfile :
		267	xmin : None,
		268	xmax : None,
		269	ymin : None,
		270	ymax : None,
		271	zmin : None,
		272	zmax : None
		273	"""
		274
		275	arrayParameters = dataOut.data_param
		276	error = arrayParameters[:,-1]
		277	indValid = numpy.where(error == 0)[0]
		278	finalMeteor = arrayParameters[indValid,:]
		279	finalAzimuth = finalMeteor[:,4]
		280	finalZenith = finalMeteor[:,5]
		281
		282	x = finalAzimuth*numpy.pi/180
		283	y = finalZenith
		284
		285
		286	#thisDatetime = dataOut.datatime
		287	thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
		288	title = wintitle + " Parameters"
		289	xlabel = "Zonal Zenith Angle (deg) "
		290	ylabel = "Meridional Zenith Angle (deg)"
		291
		292	if not self.__isConfig:
		293
		294	nplots = 1
		295
		296	self.setup(id=id,
		297	nplots=nplots,
		298	wintitle=wintitle,
		299	showprofile=showprofile,
		300	show=show)
		301
		302	self.FTP_WEI = ftp_wei
		303	self.EXP_CODE = exp_code
		304	self.SUB_EXP_CODE = sub_exp_code
		305	self.PLOT_POS = plot_pos
		306	self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
		307	self.firstdate = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
		308	self.__isConfig = True
		309
		310	self.setWinTitle(title)
		311
		312	i = 0
		313	str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
		314
		315	axes = self.axesList[i*self.__nsubplots]
		316	nevents = axes.x_buffer.shape[0] + x.shape[0]
		317	title = "Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n" %(self.firstdate,str_datetime,nevents)
		318	axes.polar(x, y,
		319	title=title, xlabel=xlabel, ylabel=ylabel,
		320	ticksize=9, cblabel='')
		321
		322	self.draw()
		323
		324	if save:
		325
		326	self.counter_imagwr += 1
		327	if (self.counter_imagwr==wr_period):
		328
		329	if figfile == None:
		330	figfile = self.getFilename(name = self.name)
		331	self.saveFigure(figpath, figfile)
		332
		333	if ftp:
		334	#provisionalmente envia archivos en el formato de la web en tiempo real
		335	name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
		336	path = '%s%03d' %(self.PREFIX, self.id)
		337	ftp_file = os.path.join(path,'ftp','%s.png'%name)
		338	self.saveFigure(figpath, ftp_file)
		339	ftp_filename = os.path.join(figpath,ftp_file)
		340
		341
		342	try:
		343	self.sendByFTP(ftp_filename, server, folder, username, password)
		344	except:
		345	self.counter_imagwr = 0
		346	raise ValueError, 'Error FTP'
		347
		348	self.counter_imagwr = 0
		349
		350
		351	class WindProfilerPlot(Figure):
		352
		353	__isConfig = None
		354	__nsubplots = None
		355
		356	WIDTHPROF = None
		357	HEIGHTPROF = None
		358	PREFIX = 'wind'
		359
		360	def __init__(self):
		361
		362	self.timerange = 26060
		363	self.isConfig = False
		364	self.__nsubplots = 1
		365
		366	self.WIDTH = 800
		367	self.HEIGHT = 150
		368	self.WIDTHPROF = 120
		369	self.HEIGHTPROF = 0
		370	self.counter_imagwr = 0
		371
		372	self.PLOT_CODE = 0
		373	self.FTP_WEI = None
		374	self.EXP_CODE = None
		375	self.SUB_EXP_CODE = None
		376	self.PLOT_POS = None
		377	self.tmin = None
		378	self.tmax = None
		379
		380	self.xmin = None
		381	self.xmax = None
		382
		383	self.figfile = None
		384
		385	def getSubplots(self):
		386
		387	ncol = 1
		388	nrow = self.nplots
		389
		390	return nrow, ncol
		391
		392	def setup(self, id, nplots, wintitle, showprofile=True, show=True):
		393
		394	self.__showprofile = showprofile
		395	self.nplots = nplots
		396
		397	ncolspan = 1
		398	colspan = 1
		399
		400	self.createFigure(id = id,
		401	wintitle = wintitle,
		402	widthplot = self.WIDTH + self.WIDTHPROF,
		403	heightplot = self.HEIGHT + self.HEIGHTPROF,
		404	show=show)
		405
		406	nrow, ncol = self.getSubplots()
		407
		408	counter = 0
		409	for y in range(nrow):
		410	if counter >= self.nplots:
		411	break
		412
		413	self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
		414	counter += 1
		415
		416	def run(self, dataOut, id, wintitle="", channelList=None,
		417	xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
		418	zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None,
		419	timerange=None, SNRthresh = None,
		420	save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
		421	server=None, folder=None, username=None, password=None,
		422	ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
		423	"""
		424
		425	Input:
		426	dataOut :
		427	id :
		428	wintitle :
		429	channelList :
		430	showProfile :
		431	xmin : None,
		432	xmax : None,
		433	ymin : None,
		434	ymax : None,
		435	zmin : None,
		436	zmax : None
		437	"""
		438
		439	if channelList == None:
		440	channelIndexList = dataOut.channelIndexList
		441	else:
		442	channelIndexList = []
		443	for channel in channelList:
		444	if channel not in dataOut.channelList:
		445	raise ValueError, "Channel %d is not in dataOut.channelList"
		446	channelIndexList.append(dataOut.channelList.index(channel))
		447
		448	if timerange != None:
		449	self.timerange = timerange
		450
		451	tmin = None
		452	tmax = None
		453
		454	x = dataOut.getTimeRange1()
		455	y = dataOut.heightRange
		456
		457	z = dataOut.winds
		458	nplots = z.shape[0] #Number of wind dimensions estimated
		459	nplotsw = nplots
		460
		461	#If there is a SNR function defined
		462	if dataOut.SNR != None:
		463	nplots += 1
		464	SNR = dataOut.SNR
		465	SNRavg = numpy.average(SNR, axis=0)
		466
		467	SNRdB = 10*numpy.log10(SNR)
		468	SNRavgdB = 10*numpy.log10(SNRavg)
		469
		470	if SNRthresh == None: SNRthresh = -5.0
		471	ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
		472
		473	for i in range(nplotsw):
		474	z[i,ind] = numpy.nan
		475
		476
		477	showprofile = False
		478	# thisDatetime = dataOut.datatime
		479	thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
		480	title = wintitle + "Wind"
		481	xlabel = ""
		482	ylabel = "Range (Km)"
		483
		484	if not self.__isConfig:
		485
		486
		487
		488	self.setup(id=id,
		489	nplots=nplots,
		490	wintitle=wintitle,
		491	showprofile=showprofile,
		492	show=show)
		493
		494	self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
		495
		496	if ymin == None: ymin = numpy.nanmin(y)
		497	if ymax == None: ymax = numpy.nanmax(y)
		498
		499	if zmax == None: zmax = numpy.nanmax(abs(z[range(2),:]))
		500	#if numpy.isnan(zmax): zmax = 50
		501	if zmin == None: zmin = -zmax
		502
		503	if nplotsw == 3:
		504	if zmax_ver == None: zmax_ver = numpy.nanmax(abs(z[2,:]))
		505	if zmin_ver == None: zmin_ver = -zmax_ver
		506
		507	if dataOut.SNR != None:
		508	if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
		509	if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
		510
		511	self.FTP_WEI = ftp_wei
		512	self.EXP_CODE = exp_code
		513	self.SUB_EXP_CODE = sub_exp_code
		514	self.PLOT_POS = plot_pos
		515
		516	self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
		517	self.__isConfig = True
		518
		519
		520	self.setWinTitle(title)
		521
		522	if ((self.xmax - x[1]) < (x[1]-x[0])):
		523	x[1] = self.xmax
		524
		525	strWind = ['Zonal', 'Meridional', 'Vertical']
		526	strCb = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
		527	zmaxVector = [zmax, zmax, zmax_ver]
		528	zminVector = [zmin, zmin, zmin_ver]
		529	windFactor = [1,1,100]
		530
		531	for i in range(nplotsw):
		532
		533	title = "%s Wind: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
		534	axes = self.axesList[i*self.__nsubplots]
		535
		536	z1 = z[i,:].reshape((1,-1))*windFactor[i]
		537
		538	axes.pcolorbuffer(x, y, z1,
		539	xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
		540	xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
		541	ticksize=9, cblabel=strCb[i], cbsize="1%", colormap="RdBu_r" )
		542
		543	if dataOut.SNR != None:
		544	i += 1
		545	title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
		546	axes = self.axesList[i*self.__nsubplots]
		547
		548	SNRavgdB = SNRavgdB.reshape((1,-1))
		549
		550	axes.pcolorbuffer(x, y, SNRavgdB,
		551	xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
		552	xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
		553	ticksize=9, cblabel='', cbsize="1%", colormap="jet")
		554
		555	self.draw()
		556
		557	if x[1] >= self.axesList[0].xmax:
		558	self.counter_imagwr = wr_period
		559	self.__isConfig = False
		560
		561
		562	if self.figfile == None:
		563	str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
		564	self.figfile = self.getFilename(name = str_datetime)
		565
		566	if figpath != '':
		567
		568	self.counter_imagwr += 1
		569	if (self.counter_imagwr>=wr_period):
		570	# store png plot to local folder
		571	self.saveFigure(figpath, self.figfile)
		572	# store png plot to FTP server according to RT-Web format
		573	name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
		574	ftp_filename = os.path.join(figpath, name)
		575	self.saveFigure(figpath, ftp_filename)
		576
		577	self.counter_imagwr = 0
		578
		579	No newline at end of file

schainpy/model/proc/jroproc_correlation.py created 644 +246 0

@@ -0,0 +1,246
	1	import numpy
	2
	3	from jroproc_base import ProcessingUnit, Operation
	4	from model.data.jrodata import Correlation
	5
	6	class CorrelationProc(ProcessingUnit):
	7
	8	def __init__(self):
	9
	10	ProcessingUnit.__init__(self)
	11
	12	self.objectDict = {}
	13	self.buffer = None
	14	self.firstdatatime = None
	15	self.profIndex = 0
	16	self.dataOut = Correlation()
	17
	18	def __updateObjFromInput(self):
	19
	20	self.dataOut.timeZone = self.dataIn.timeZone
	21	self.dataOut.dstFlag = self.dataIn.dstFlag
	22	self.dataOut.errorCount = self.dataIn.errorCount
	23	self.dataOut.useLocalTime = self.dataIn.useLocalTime
	24
	25	self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
	26	self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
	27	self.dataOut.channelList = self.dataIn.channelList
	28	self.dataOut.heightList = self.dataIn.heightList
	29	self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
	30	# self.dataOut.nHeights = self.dataIn.nHeights
	31	# self.dataOut.nChannels = self.dataIn.nChannels
	32	self.dataOut.nBaud = self.dataIn.nBaud
	33	self.dataOut.nCode = self.dataIn.nCode
	34	self.dataOut.code = self.dataIn.code
	35	# self.dataOut.nProfiles = self.dataOut.nFFTPoints
	36	self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
	37	self.dataOut.utctime = self.firstdatatime
	38	self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
	39	self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
	40	# self.dataOut.nCohInt = self.dataIn.nCohInt
	41	# self.dataOut.nIncohInt = 1
	42	self.dataOut.ippSeconds = self.dataIn.ippSeconds
	43	# self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
	44
	45	self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nPoints
	46
	47
	48	def removeDC(self, jspectra):
	49
	50	nChannel = jspectra.shape[0]
	51
	52	for i in range(nChannel):
	53	jspectra_tmp = jspectra[i,:,:]
	54	jspectra_DC = numpy.mean(jspectra_tmp,axis = 0)
	55
	56	jspectra_tmp = jspectra_tmp - jspectra_DC
	57	jspectra[i,:,:] = jspectra_tmp
	58
	59	return jspectra
	60
	61
	62	def removeNoise(self, mode = 2):
	63	indR = numpy.where(self.dataOut.lagR == 0)[0][0]
	64	indT = numpy.where(self.dataOut.lagT == 0)[0][0]
	65
	66	jspectra = self.dataOut.data_corr[:,:,indR,:]
	67
	68	num_chan = jspectra.shape[0]
	69	num_hei = jspectra.shape[2]
	70
	71	freq_dc = indT
	72	ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
	73
	74	NPot = self.dataOut.getNoise(mode)
	75	jspectra[:,freq_dc,:] = jspectra[:,freq_dc,:] - NPot
	76	SPot = jspectra[:,freq_dc,:]
	77	pairsAutoCorr = self.dataOut.getPairsAutoCorr()
	78	# self.dataOut.signalPotency = SPot
	79	self.dataOut.noise = NPot
	80	self.dataOut.SNR = (SPot/NPot)[pairsAutoCorr]
	81	self.dataOut.data_corr[:,:,indR,:] = jspectra
	82
	83	return 1
	84
	85
	86	def calculateNormFactor(self):
	87
	88	pairsList = self.dataOut.pairsList
	89	pairsAutoCorr = self.dataOut.pairsAutoCorr
	90	nHeights = self.dataOut.nHeights
	91	nPairs = len(pairsList)
	92	normFactor = numpy.zeros((nPairs,nHeights))
	93
	94	indR = numpy.where(self.dataOut.lagR == 0)[0][0]
	95	indT = numpy.where(self.dataOut.lagT == 0)[0][0]
	96
	97	for l in range(len(pairsList)):
	98	firstChannel = pairsList[l][0]
	99	secondChannel = pairsList[l][1]
	100
	101	AC1 = pairsAutoCorr[firstChannel]
	102	AC2 = pairsAutoCorr[secondChannel]
	103
	104	if (AC1 >= 0 and AC2 >= 0):
	105
	106	data1 = numpy.abs(self.dataOut.data_corr[AC1,:,indR,:])
	107	data2 = numpy.abs(self.dataOut.data_corr[AC2,:,indR,:])
	108	maxim1 = data1.max(axis = 0)
	109	maxim2 = data1.max(axis = 0)
	110	maxim = numpy.sqrt(maxim1*maxim2)
	111	else:
	112	#In case there is no autocorrelation for the pair
	113	data = numpy.abs(self.dataOut.data_corr[l,:,indR,:])
	114	maxim = numpy.max(data, axis = 0)
	115
	116	normFactor[l,:] = maxim
	117
	118	self.dataOut.normFactor = normFactor
	119
	120	return 1
	121
	122	def run(self, lagT=None, lagR=None, pairsList=None,
	123	nPoints=None, nAvg=None, bufferSize=None,
	124	fullT = False, fullR = False, removeDC = False):
	125
	126	self.dataOut.flagNoData = True
	127
	128	if self.dataIn.type == "Correlation":
	129
	130	self.dataOut.copy(self.dataIn)
	131
	132	return
	133
	134	if self.dataIn.type == "Voltage":
	135
	136	if pairsList == None:
	137	pairsList = [numpy.array([0,0])]
	138
	139	if nPoints == None:
	140	nPoints = 128
	141	#------------------------------------------------------------
	142	#Condicionales para calcular Correlaciones en Tiempo y Rango
	143	if fullT:
	144	lagT = numpy.arange(nPoints*2 - 1) - nPoints + 1
	145	elif lagT == None:
	146	lagT = numpy.array([0])
	147	else:
	148	lagT = numpy.array(lagT)
	149
	150	if fullR:
	151	lagR = numpy.arange(self.dataOut.nHeights)
	152	elif lagR == None:
	153	lagR = numpy.array([0])
	154	#-------------------------------------------------------------
	155
	156	if nAvg == None:
	157	nAvg = 1
	158
	159	if bufferSize == None:
	160	bufferSize = 0
	161
	162	deltaH = self.dataIn.heightList[1] - self.dataIn.heightList[0]
	163	self.dataOut.lagR = numpy.round(numpy.array(lagR)/deltaH)
	164	self.dataOut.pairsList = pairsList
	165	self.dataOut.nPoints = nPoints
	166	# channels = numpy.sort(list(set(list(itertools.chain.from_iterable(pairsList)))))
	167
	168	if self.buffer == None:
	169
	170	self.buffer = numpy.zeros((self.dataIn.nChannels,self.dataIn.nProfiles,self.dataIn.nHeights),dtype='complex')
	171
	172
	173	self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()[:,:]
	174
	175	self.profIndex += 1
	176
	177	if self.firstdatatime == None:
	178
	179	self.firstdatatime = self.dataIn.utctime
	180
	181	if self.profIndex == nPoints:
	182
	183	tmp = self.buffer[:,0:nPoints,:]
	184	self.buffer = None
	185	self.buffer = tmp
	186
	187	#--------------- Remover DC ------------
	188	if removeDC:
	189	self.buffer = self.removeDC(self.buffer)
	190	#---------------------------------------------
	191	self.dataOut.data_volts = self.buffer
	192	self.__updateObjFromInput()
	193	self.dataOut.data_corr = numpy.zeros((len(pairsList),
	194	len(lagT),len(lagR),
	195	self.dataIn.nHeights),
	196	dtype='complex')
	197
	198	for l in range(len(pairsList)):
	199
	200	firstChannel = pairsList[l][0]
	201	secondChannel = pairsList[l][1]
	202
	203	tmp = None
	204	tmp = numpy.zeros((len(lagT),len(lagR),self.dataIn.nHeights),dtype='complex')
	205
	206	for t in range(len(lagT)):
	207
	208	for r in range(len(lagR)):
	209
	210	idxT = lagT[t]
	211	idxR = lagR[r]
	212
	213	if idxT >= 0:
	214	vStacked = numpy.vstack((self.buffer[secondChannel,idxT:,:],
	215	numpy.zeros((idxT,self.dataIn.nHeights),dtype='complex')))
	216	else:
	217	vStacked = numpy.vstack((numpy.zeros((-idxT,self.dataIn.nHeights),dtype='complex'),
	218	self.buffer[secondChannel,:(nPoints + idxT),:]))
	219
	220	if idxR >= 0:
	221	hStacked = numpy.hstack((vStacked[:,idxR:],numpy.zeros((nPoints,idxR),dtype='complex')))
	222	else:
	223	hStacked = numpy.hstack((numpy.zeros((nPoints,-idxR),dtype='complex'),vStacked[:,(self.dataOut.nHeights + idxR)]))
	224
	225
	226	tmp[t,r,:] = numpy.sum((numpy.conjugate(self.buffer[firstChannel,:,:])*hStacked),axis=0)
	227
	228
	229	hStacked = None
	230	vStacked = None
	231
	232	self.dataOut.data_corr[l,:,:,:] = tmp[:,:,:]
	233
	234	#Se Calcula los factores de Normalizacion
	235	self.dataOut.pairsAutoCorr = self.dataOut.getPairsAutoCorr()
	236	self.dataOut.lagT = lagTself.dataIn.ippSecondsself.dataIn.nCohInt
	237	self.dataOut.lagR = lagR
	238
	239	self.calculateNormFactor()
	240
	241	self.dataOut.flagNoData = False
	242	self.buffer = None
	243	self.firstdatatime = None
	244	self.profIndex = 0
	245
	246	return No newline at end of file

schainpy/model/proc/jroproc_parameters.py created 644 +1521 0

This diff has been collapsed as it changes many lines, (1521 lines changed) Show them Hide them
	@@ -0,0 +1,1521
		1	import numpy
		2	import math
		3	from scipy import optimize
		4	from scipy import interpolate
		5	from scipy import signal
		6	from scipy import stats
		7	import re
		8	import datetime
		9	import copy
		10
		11
		12	from jroproc_base import ProcessingUnit, Operation
		13	from model.data.jrodata import Parameters
		14
		15
		16	class ParametersProc(ProcessingUnit):
		17
		18	nSeconds = None
		19
		20	def __init__(self):
		21	ProcessingUnit.__init__(self)
		22
		23	self.objectDict = {}
		24	self.buffer = None
		25	self.firstdatatime = None
		26	self.profIndex = 0
		27	self.dataOut = Parameters()
		28
		29	def __updateObjFromInput(self):
		30
		31	self.dataOut.inputUnit = self.dataIn.type
		32
		33	self.dataOut.timeZone = self.dataIn.timeZone
		34	self.dataOut.dstFlag = self.dataIn.dstFlag
		35	self.dataOut.errorCount = self.dataIn.errorCount
		36	self.dataOut.useLocalTime = self.dataIn.useLocalTime
		37
		38	self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
		39	self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
		40	self.dataOut.channelList = self.dataIn.channelList
		41	self.dataOut.heightList = self.dataIn.heightList
		42	self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
		43	# self.dataOut.nHeights = self.dataIn.nHeights
		44	# self.dataOut.nChannels = self.dataIn.nChannels
		45	self.dataOut.nBaud = self.dataIn.nBaud
		46	self.dataOut.nCode = self.dataIn.nCode
		47	self.dataOut.code = self.dataIn.code
		48	# self.dataOut.nProfiles = self.dataOut.nFFTPoints
		49	self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
		50	self.dataOut.utctime = self.firstdatatime
		51	self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
		52	self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
		53	# self.dataOut.nCohInt = self.dataIn.nCohInt
		54	# self.dataOut.nIncohInt = 1
		55	self.dataOut.ippSeconds = self.dataIn.ippSeconds
		56	# self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
		57	self.dataOut.timeInterval = self.dataIn.timeInterval
		58	self.dataOut.heightRange = self.dataIn.getHeiRange()
		59	self.dataOut.frequency = self.dataIn.frequency
		60
		61	def run(self, nSeconds = None, nProfiles = None):
		62
		63	self.dataOut.flagNoData = True
		64
		65	if self.firstdatatime == None:
		66	self.firstdatatime = self.dataIn.utctime
		67
		68	#---------------------- Voltage Data ---------------------------
		69
		70	if self.dataIn.type == "Voltage":
		71	if nSeconds != None:
		72	self.nSeconds = nSeconds
		73	self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))
		74
		75	if self.buffer == None:
		76	self.buffer = numpy.zeros((self.dataIn.nChannels,
		77	self.nProfiles,
		78	self.dataIn.nHeights),
		79	dtype='complex')
		80
		81	self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
		82	self.profIndex += 1
		83
		84	if self.profIndex == self.nProfiles:
		85
		86	self.__updateObjFromInput()
		87	self.dataOut.data_pre = self.buffer.copy()
		88	self.dataOut.paramInterval = nSeconds
		89	self.dataOut.flagNoData = False
		90
		91	self.buffer = None
		92	self.firstdatatime = None
		93	self.profIndex = 0
		94
		95	#---------------------- Spectra Data ---------------------------
		96
		97	if self.dataIn.type == "Spectra":
		98	self.dataOut.data_pre = self.dataIn.data_spc.copy()
		99	self.dataOut.abscissaRange = self.dataIn.getVelRange(1)
		100	self.dataOut.noise = self.dataIn.getNoise()
		101	self.dataOut.normFactor = self.dataIn.normFactor
		102
		103	self.__updateObjFromInput()
		104	self.dataOut.flagNoData = False
		105	self.firstdatatime = None
		106
		107	#---------------------- Correlation Data ---------------------------
		108
		109	if self.dataIn.type == "Correlation":
		110	lagRRange = self.dataIn.lagR
		111	indR = numpy.where(lagRRange == 0)[0][0]
		112
		113	self.dataOut.data_pre = self.dataIn.data_corr.copy()[:,:,indR,:]
		114	self.dataOut.abscissaRange = self.dataIn.getLagTRange(1)
		115	self.dataOut.noise = self.dataIn.noise
		116	self.dataOut.normFactor = self.dataIn.normFactor
		117	self.dataOut.SNR = self.dataIn.SNR
		118	self.dataOut.pairsList = self.dataIn.pairsList
		119
		120	self.__updateObjFromInput()
		121	self.dataOut.flagNoData = False
		122	self.firstdatatime = None
		123
		124	#------------------- Get Moments ----------------------------------
		125	def GetMoments(self, channelList = None):
		126	'''
		127	Function GetMoments()
		128
		129	Input:
		130	channelList : simple channel list to select e.g. [2,3,7]
		131	self.dataOut.data_pre
		132	self.dataOut.abscissaRange
		133	self.dataOut.noise
		134
		135	Affected:
		136	self.dataOut.data_param
		137	self.dataOut.SNR
		138
		139	'''
		140	data = self.dataOut.data_pre
		141	absc = self.dataOut.abscissaRange[:-1]
		142	noise = self.dataOut.noise
		143
		144	data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
		145
		146	if channelList== None: channelList = self.dataOut.channelList
		147
		148	for ind in channelList:
		149	data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
		150
		151	self.dataOut.data_param = data_param[:,1:]
		152	self.dataOut.SNR = data_param[:,0]
		153	return
		154
		155	def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
		156
		157	if (nicoh == None): nicoh = 1
		158	if (graph == None): graph = 0
		159	if (smooth == None): smooth = 0
		160	elif (self.smooth < 3): smooth = 0
		161
		162	if (type1 == None): type1 = 0
		163	if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
		164	if (snrth == None): snrth = -3
		165	if (dc == None): dc = 0
		166	if (aliasing == None): aliasing = 0
		167	if (oldfd == None): oldfd = 0
		168	if (wwauto == None): wwauto = 0
		169
		170	if (n0 < 1.e-20): n0 = 1.e-20
		171
		172	freq = oldfreq
		173	vec_power = numpy.zeros(oldspec.shape[1])
		174	vec_fd = numpy.zeros(oldspec.shape[1])
		175	vec_w = numpy.zeros(oldspec.shape[1])
		176	vec_snr = numpy.zeros(oldspec.shape[1])
		177
		178	for ind in range(oldspec.shape[1]):
		179
		180	spec = oldspec[:,ind]
		181	aux = spec*fwindow
		182	max_spec = aux.max()
		183	m = list(aux).index(max_spec)
		184
		185	#Smooth
		186	if (smooth == 0): spec2 = spec
		187	else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
		188
		189	# Calculo de Momentos
		190	bb = spec2[range(m,spec2.size)]
		191	bb = (bb<n0).nonzero()
		192	bb = bb[0]
		193
		194	ss = spec2[range(0,m + 1)]
		195	ss = (ss<n0).nonzero()
		196	ss = ss[0]
		197
		198	if (bb.size == 0):
		199	bb0 = spec.size - 1 - m
		200	else:
		201	bb0 = bb[0] - 1
		202	if (bb0 < 0):
		203	bb0 = 0
		204
		205	if (ss.size == 0): ss1 = 1
		206	else: ss1 = max(ss) + 1
		207
		208	if (ss1 > m): ss1 = m
		209
		210	valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
		211	power = ((spec2[valid] - n0)*fwindow[valid]).sum()
		212	fd = ((spec2[valid]- n0)freq[valid]fwindow[valid]).sum()/power
		213	w = math.sqrt(((spec2[valid] - n0)fwindow[valid](freq[valid]- fd)**2).sum()/power)
		214	snr = (spec2.mean()-n0)/n0
		215
		216	if (snr < 1.e-20) :
		217	snr = 1.e-20
		218
		219	vec_power[ind] = power
		220	vec_fd[ind] = fd
		221	vec_w[ind] = w
		222	vec_snr[ind] = snr
		223
		224	moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
		225	return moments
		226
		227	#------------------- Get Lags ----------------------------------
		228
		229	def GetLags(self):
		230	'''
		231	Function GetMoments()
		232
		233	Input:
		234	self.dataOut.data_pre
		235	self.dataOut.abscissaRange
		236	self.dataOut.noise
		237	self.dataOut.normFactor
		238	self.dataOut.SNR
		239	self.dataOut.pairsList
		240	self.dataOut.nChannels
		241
		242	Affected:
		243	self.dataOut.data_param
		244
		245	'''
		246	data = self.dataOut.data_pre
		247	normFactor = self.dataOut.normFactor
		248	nHeights = self.dataOut.nHeights
		249	absc = self.dataOut.abscissaRange[:-1]
		250	noise = self.dataOut.noise
		251	SNR = self.dataOut.SNR
		252	pairsList = self.dataOut.pairsList
		253	nChannels = self.dataOut.nChannels
		254	pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
		255	self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
		256
		257	dataNorm = numpy.abs(data)
		258	for l in range(len(pairsList)):
		259	dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]
		260
		261	self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)
		262	self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, pairsAutoCorr, absc)
		263	return
		264
		265	def __getPairsAutoCorr(self, pairsList, nChannels):
		266
		267	pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
		268
		269	for l in range(len(pairsList)):
		270	firstChannel = pairsList[l][0]
		271	secondChannel = pairsList[l][1]
		272
		273	#Obteniendo pares de Autocorrelacion
		274	if firstChannel == secondChannel:
		275	pairsAutoCorr[firstChannel] = int(l)
		276
		277	pairsAutoCorr = pairsAutoCorr.astype(int)
		278
		279	pairsCrossCorr = range(len(pairsList))
		280	pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
		281
		282	return pairsAutoCorr, pairsCrossCorr
		283
		284	def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):
		285
		286	Pt0 = data.shape[1]/2
		287	#Funcion de Autocorrelacion
		288	dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0)
		289
		290	#Obtencion Indice de TauCross
		291	indCross = data[pairsCrossCorr,:,:].argmax(axis = 1)
		292	#Obtencion Indice de TauAuto
		293	indAuto = numpy.zeros(indCross.shape,dtype = 'int')
		294	CCValue = data[pairsCrossCorr,Pt0,:]
		295	for i in range(pairsCrossCorr.size):
		296	indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0)
		297
		298	#Obtencion de TauCross y TauAuto
		299	tauCross = lagTRange[indCross]
		300	tauAuto = lagTRange[indAuto]
		301
		302	Nan1, Nan2 = numpy.where(tauCross == lagTRange[0])
		303
		304	tauCross[Nan1,Nan2] = numpy.nan
		305	tauAuto[Nan1,Nan2] = numpy.nan
		306	tau = numpy.vstack((tauCross,tauAuto))
		307
		308	return tau
		309
		310	def __calculateLag1Phase(self, data, pairs, lagTRange):
		311	data1 = stats.nanmean(data[pairs,:,:], axis = 0)
		312	lag1 = numpy.where(lagTRange == 0)[0][0] + 1
		313
		314	phase = numpy.angle(data1[lag1,:])
		315
		316	return phase
		317	#------------------- Detect Meteors ------------------------------
		318
		319	def DetectMeteors(self, hei_ref = None, tauindex = 0,
		320	predefinedPhaseShifts = None, centerReceiverIndex = 2,
		321	cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
		322	noise_timeStep = 4, noise_multiple = 4,
		323	multDet_timeLimit = 1, multDet_rangeLimit = 3,
		324	phaseThresh = 20, SNRThresh = 8,
		325	hmin = 70, hmax=110, azimuth = 0) :
		326
		327	'''
		328	Function DetectMeteors()
		329	Project developed with paper:
		330	HOLDSWORTH ET AL. 2004
		331
		332	Input:
		333	self.dataOut.data_pre
		334
		335	centerReceiverIndex: From the channels, which is the center receiver
		336
		337	hei_ref: Height reference for the Beacon signal extraction
		338	tauindex:
		339	predefinedPhaseShifts: Predefined phase offset for the voltge signals
		340
		341	cohDetection: Whether to user Coherent detection or not
		342	cohDet_timeStep: Coherent Detection calculation time step
		343	cohDet_thresh: Coherent Detection phase threshold to correct phases
		344
		345	noise_timeStep: Noise calculation time step
		346	noise_multiple: Noise multiple to define signal threshold
		347
		348	multDet_timeLimit: Multiple Detection Removal time limit in seconds
		349	multDet_rangeLimit: Multiple Detection Removal range limit in km
		350
		351	phaseThresh: Maximum phase difference between receiver to be consider a meteor
		352	SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
		353
		354	hmin: Minimum Height of the meteor to use it in the further wind estimations
		355	hmax: Maximum Height of the meteor to use it in the further wind estimations
		356	azimuth: Azimuth angle correction
		357
		358	Affected:
		359	self.dataOut.data_param
		360
		361	Rejection Criteria (Errors):
		362	0: No error; analysis OK
		363	1: SNR < SNR threshold
		364	2: angle of arrival (AOA) ambiguously determined
		365	3: AOA estimate not feasible
		366	4: Large difference in AOAs obtained from different antenna baselines
		367	5: echo at start or end of time series
		368	6: echo less than 5 examples long; too short for analysis
		369	7: echo rise exceeds 0.3s
		370	8: echo decay time less than twice rise time
		371	9: large power level before echo
		372	10: large power level after echo
		373	11: poor fit to amplitude for estimation of decay time
		374	12: poor fit to CCF phase variation for estimation of radial drift velocity
		375	13: height unresolvable echo: not valid height within 70 to 110 km
		376	14: height ambiguous echo: more then one possible height within 70 to 110 km
		377	15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
		378	16: oscilatory echo, indicating event most likely not an underdense echo
		379
		380	17: phase difference in meteor Reestimation
		381
		382	Data Storage:
		383	Meteors for Wind Estimation (8):
		384	Day Hour \| Range Height
		385	Azimuth Zenith errorCosDir
		386	VelRad errorVelRad
		387	TypeError
		388
		389	'''
		390	#Get Beacon signal
		391	newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
		392
		393	if hei_ref != None:
		394	newheis = numpy.where(self.dataOut.heightList>hei_ref)
		395
		396	heiRang = self.dataOut.getHeiRange()
		397	#Pairs List
		398	pairslist = []
		399	nChannel = self.dataOut.nChannels
		400	for i in range(nChannel):
		401	if i != centerReceiverIndex:
		402	pairslist.append((centerReceiverIndex,i))
		403
		404	#**************REMOVING HARDWARE PHASE DIFFERENCES*************
		405	# see if the user put in pre defined phase shifts
		406	voltsPShift = self.dataOut.data_pre.copy()
		407
		408	if predefinedPhaseShifts != None:
		409	hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
		410	else:
		411	#get hardware phase shifts using beacon signal
		412	hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
		413	hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
		414
		415	voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
		416	for i in range(self.dataOut.data_pre.shape[0]):
		417	voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
		418	#****************END OF REMOVING HARDWARE PHASE DIFFERENCES*******
		419
		420	#Remove DC
		421	voltsDC = numpy.mean(voltsPShift,1)
		422	voltsDC = numpy.mean(voltsDC,1)
		423	for i in range(voltsDC.shape[0]):
		424	voltsPShift[i] = voltsPShift[i] - voltsDC[i]
		425
		426	#Don't considerate last heights, theyre used to calculate Hardware Phase Shift
		427	voltsPShift = voltsPShift[:,:,:newheis[0][0]]
		428
		429	#********** FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) ********
		430	#Coherent Detection
		431	if cohDetection:
		432	#use coherent detection to get the net power
		433	cohDet_thresh = cohDet_thresh*numpy.pi/180
		434	voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist, cohDet_thresh)
		435
		436	#Non-coherent detection!
		437	powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
		438	#******** END OF COH/NON-COH POWER CALCULATION********************
		439
		440	#******** FIND THE NOISE LEVEL AND POSSIBLE METEORS **************
		441	#Get noise
		442	noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)
		443	# noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
		444	#Get signal threshold
		445	signalThresh = noise_multiple*noise
		446	#Meteor echoes detection
		447	listMeteors = self.__findMeteors(powerNet, signalThresh)
		448	#***** END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION ********
		449
		450	#************ REMOVE MULTIPLE DETECTIONS (3.5) *************************
		451	#Parameters
		452	heiRange = self.dataOut.getHeiRange()
		453	rangeInterval = heiRange[1] - heiRange[0]
		454	rangeLimit = multDet_rangeLimit/rangeInterval
		455	timeLimit = multDet_timeLimit/self.dataOut.timeInterval
		456	#Multiple detection removals
		457	listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
		458	#********** END OF REMOVE MULTIPLE DETECTIONS ********************
		459
		460	#******************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ******************
		461	#Parameters
		462	phaseThresh = phaseThresh*numpy.pi/180
		463	thresh = [phaseThresh, noise_multiple, SNRThresh]
		464	#Meteor reestimation (Errors N 1, 6, 12, 17)
		465	listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
		466	# listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
		467	#Estimation of decay times (Errors N 7, 8, 11)
		468	listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
		469	#***************** END OF METEOR REESTIMATION *****************
		470
		471	#******************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) ************************
		472	#Calculating Radial Velocity (Error N 15)
		473	radialStdThresh = 10
		474	listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
		475
		476	if len(listMeteors4) > 0:
		477	#Setting New Array
		478	date = repr(self.dataOut.datatime)
		479	arrayMeteors4, arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
		480
		481	#Calculate AOA (Error N 3, 4)
		482	#JONES ET AL. 1998
		483	AOAthresh = numpy.pi/8
		484	error = arrayParameters[:,-1]
		485	phases = -arrayMeteors4[:,9:13]
		486	pairsList = []
		487	pairsList.append((0,3))
		488	pairsList.append((1,2))
		489	arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
		490
		491	#Calculate Heights (Error N 13 and 14)
		492	error = arrayParameters[:,-1]
		493	Ranges = arrayParameters[:,2]
		494	zenith = arrayParameters[:,5]
		495	arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
		496	#******************* END OF PARAMETERS CALCULATION ************************
		497
		498	#*************************+ SAVE DATA IN HDF5 FORMAT ********************
		499	self.dataOut.data_param = arrayParameters
		500
		501	return
		502
		503	def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
		504
		505	minIndex = min(newheis[0])
		506	maxIndex = max(newheis[0])
		507
		508	voltage = voltage0[:,:,minIndex:maxIndex+1]
		509	nLength = voltage.shape[1]/n
		510	nMin = 0
		511	nMax = 0
		512	phaseOffset = numpy.zeros((len(pairslist),n))
		513
		514	for i in range(n):
		515	nMax += nLength
		516	phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
		517	phaseCCF = numpy.mean(phaseCCF, axis = 2)
		518	phaseOffset[:,i] = phaseCCF.transpose()
		519	nMin = nMax
		520	# phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
		521
		522	#Remove Outliers
		523	factor = 2
		524	wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
		525	dw = numpy.std(wt,axis = 1)
		526	dw = dw.reshape((dw.size,1))
		527	ind = numpy.where(numpy.logical_or(wt>dwfactor,wt<-dwfactor))
		528	phaseOffset[ind] = numpy.nan
		529	phaseOffset = stats.nanmean(phaseOffset, axis=1)
		530
		531	return phaseOffset
		532
		533	def __shiftPhase(self, data, phaseShift):
		534	#this will shift the phase of a complex number
		535	dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
		536	return dataShifted
		537
		538	def __estimatePhaseDifference(self, array, pairslist):
		539	nChannel = array.shape[0]
		540	nHeights = array.shape[2]
		541	numPairs = len(pairslist)
		542	# phaseCCF = numpy.zeros((nChannel, 5, nHeights))
		543	phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
		544
		545	#Correct phases
		546	derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
		547	indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
		548
		549	if indDer[0].shape[0] > 0:
		550	for i in range(indDer[0].shape[0]):
		551	signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
		552	phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo2numpy.pi
		553
		554	# for j in range(numSides):
		555	# phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
		556	# phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
		557	#
		558	#Linear
		559	phaseInt = numpy.zeros((numPairs,1))
		560	angAllCCF = phaseCCF[:,[0,1,3,4],0]
		561	for j in range(numPairs):
		562	fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
		563	phaseInt[j] = fit[1]
		564	#Phase Differences
		565	phaseDiff = phaseInt - phaseCCF[:,2,:]
		566	phaseArrival = phaseInt.reshape(phaseInt.size)
		567
		568	#Dealias
		569	indAlias = numpy.where(phaseArrival > numpy.pi)
		570	phaseArrival[indAlias] -= 2*numpy.pi
		571	indAlias = numpy.where(phaseArrival < -numpy.pi)
		572	phaseArrival[indAlias] += 2*numpy.pi
		573
		574	return phaseDiff, phaseArrival
		575
		576	def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
		577	#this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
		578	#find the phase shifts of each channel over 1 second intervals
		579	#only look at ranges below the beacon signal
		580	numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
		581	numBlocks = int(volts.shape[1]/numProfPerBlock)
		582	numHeights = volts.shape[2]
		583	nChannel = volts.shape[0]
		584	voltsCohDet = volts.copy()
		585
		586	pairsarray = numpy.array(pairslist)
		587	indSides = pairsarray[:,1]
		588	# indSides = numpy.array(range(nChannel))
		589	# indSides = numpy.delete(indSides, indCenter)
		590	#
		591	# listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
		592	listBlocks = numpy.array_split(volts, numBlocks, 1)
		593
		594	startInd = 0
		595	endInd = 0
		596
		597	for i in range(numBlocks):
		598	startInd = endInd
		599	endInd = endInd + listBlocks[i].shape[1]
		600
		601	arrayBlock = listBlocks[i]
		602	# arrayBlockCenter = listCenter[i]
		603
		604	#Estimate the Phase Difference
		605	phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
		606	#Phase Difference RMS
		607	arrayPhaseRMS = numpy.abs(phaseDiff)
		608	phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
		609	indPhase = numpy.where(phaseRMSaux==4)
		610	#Shifting
		611	if indPhase[0].shape[0] > 0:
		612	for j in range(indSides.size):
		613	arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
		614	voltsCohDet[:,startInd:endInd,:] = arrayBlock
		615
		616	return voltsCohDet
		617
		618	def __calculateCCF(self, volts, pairslist ,laglist):
		619
		620	nHeights = volts.shape[2]
		621	nPoints = volts.shape[1]
		622	voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
		623
		624	for i in range(len(pairslist)):
		625	volts1 = volts[pairslist[i][0]]
		626	volts2 = volts[pairslist[i][1]]
		627
		628	for t in range(len(laglist)):
		629	idxT = laglist[t]
		630	if idxT >= 0:
		631	vStacked = numpy.vstack((volts2[idxT:,:],
		632	numpy.zeros((idxT, nHeights),dtype='complex')))
		633	else:
		634	vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
		635	volts2[:(nPoints + idxT),:]))
		636	voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
		637
		638	vStacked = None
		639	return voltsCCF
		640
		641	def __getNoise(self, power, timeSegment, timeInterval):
		642	numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
		643	numBlocks = int(power.shape[0]/numProfPerBlock)
		644	numHeights = power.shape[1]
		645
		646	listPower = numpy.array_split(power, numBlocks, 0)
		647	noise = numpy.zeros((power.shape[0], power.shape[1]))
		648	noise1 = numpy.zeros((power.shape[0], power.shape[1]))
		649
		650	startInd = 0
		651	endInd = 0
		652
		653	for i in range(numBlocks): #split por canal
		654	startInd = endInd
		655	endInd = endInd + listPower[i].shape[0]
		656
		657	arrayBlock = listPower[i]
		658	noiseAux = numpy.mean(arrayBlock, 0)
		659	# noiseAux = numpy.median(noiseAux)
		660	# noiseAux = numpy.mean(arrayBlock)
		661	noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
		662
		663	noiseAux1 = numpy.mean(arrayBlock)
		664	noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
		665
		666	return noise, noise1
		667
		668	def __findMeteors(self, power, thresh):
		669	nProf = power.shape[0]
		670	nHeights = power.shape[1]
		671	listMeteors = []
		672
		673	for i in range(nHeights):
		674	powerAux = power[:,i]
		675	threshAux = thresh[:,i]
		676
		677	indUPthresh = numpy.where(powerAux > threshAux)[0]
		678	indDNthresh = numpy.where(powerAux <= threshAux)[0]
		679
		680	j = 0
		681
		682	while (j < indUPthresh.size - 2):
		683	if (indUPthresh[j + 2] == indUPthresh[j] + 2):
		684	indDNAux = numpy.where(indDNthresh > indUPthresh[j])
		685	indDNthresh = indDNthresh[indDNAux]
		686
		687	if (indDNthresh.size > 0):
		688	indEnd = indDNthresh[0] - 1
		689	indInit = indUPthresh[j]
		690
		691	meteor = powerAux[indInit:indEnd + 1]
		692	indPeak = meteor.argmax() + indInit
		693	FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
		694
		695	listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
		696	j = numpy.where(indUPthresh == indEnd)[0] + 1
		697	else: j+=1
		698	else: j+=1
		699
		700	return listMeteors
		701
		702	def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
		703
		704	arrayMeteors = numpy.asarray(listMeteors)
		705	listMeteors1 = []
		706
		707	while arrayMeteors.shape[0] > 0:
		708	FLAs = arrayMeteors[:,4]
		709	maxFLA = FLAs.argmax()
		710	listMeteors1.append(arrayMeteors[maxFLA,:])
		711
		712	MeteorInitTime = arrayMeteors[maxFLA,1]
		713	MeteorEndTime = arrayMeteors[maxFLA,3]
		714	MeteorHeight = arrayMeteors[maxFLA,0]
		715
		716	#Check neighborhood
		717	maxHeightIndex = MeteorHeight + rangeLimit
		718	minHeightIndex = MeteorHeight - rangeLimit
		719	minTimeIndex = MeteorInitTime - timeLimit
		720	maxTimeIndex = MeteorEndTime + timeLimit
		721
		722	#Check Heights
		723	indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
		724	indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
		725	indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
		726
		727	arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
		728
		729	return listMeteors1
		730
		731	def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
		732	numHeights = volts.shape[2]
		733	nChannel = volts.shape[0]
		734
		735	thresholdPhase = thresh[0]
		736	thresholdNoise = thresh[1]
		737	thresholdDB = float(thresh[2])
		738
		739	thresholdDB1 = 10**(thresholdDB/10)
		740	pairsarray = numpy.array(pairslist)
		741	indSides = pairsarray[:,1]
		742
		743	pairslist1 = list(pairslist)
		744	pairslist1.append((0,1))
		745	pairslist1.append((3,4))
		746
		747	listMeteors1 = []
		748	listPowerSeries = []
		749	listVoltageSeries = []
		750	#volts has the war data
		751
		752	if frequency == 30e6:
		753	timeLag = 4510*-3
		754	else:
		755	timeLag = 1510*-3
		756	lag = numpy.ceil(timeLag/timeInterval)
		757
		758	for i in range(len(listMeteors)):
		759
		760	###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
		761	meteorAux = numpy.zeros(16)
		762
		763	#Loading meteor Data (mHeight, mStart, mPeak, mEnd)
		764	mHeight = listMeteors[i][0]
		765	mStart = listMeteors[i][1]
		766	mPeak = listMeteors[i][2]
		767	mEnd = listMeteors[i][3]
		768
		769	#get the volt data between the start and end times of the meteor
		770	meteorVolts = volts[:,mStart:mEnd+1,mHeight]
		771	meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
		772
		773	#3.6. Phase Difference estimation
		774	phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
		775
		776	#3.7. Phase difference removal & meteor start, peak and end times reestimated
		777	#meteorVolts0.- all Channels, all Profiles
		778	meteorVolts0 = volts[:,:,mHeight]
		779	meteorThresh = noise[:,mHeight]*thresholdNoise
		780	meteorNoise = noise[:,mHeight]
		781	meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
		782	powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
		783
		784	#Times reestimation
		785	mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
		786	if mStart1.size > 0:
		787	mStart1 = mStart1[-1] + 1
		788
		789	else:
		790	mStart1 = mPeak
		791
		792	mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
		793	mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
		794	if mEndDecayTime1.size == 0:
		795	mEndDecayTime1 = powerNet0.size
		796	else:
		797	mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
		798	# mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
		799
		800	#meteorVolts1.- all Channels, from start to end
		801	meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
		802	meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
		803	if meteorVolts2.shape[1] == 0:
		804	meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
		805	meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
		806	meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
		807	##################### END PARAMETERS REESTIMATION #########################
		808
		809	##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
		810	# if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
		811	if meteorVolts2.shape[1] > 0:
		812	#Phase Difference re-estimation
		813	phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
		814	# phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
		815	meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
		816	phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
		817	meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
		818
		819	#Phase Difference RMS
		820	phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
		821	powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
		822	#Data from Meteor
		823	mPeak1 = powerNet1.argmax() + mStart1
		824	mPeakPower1 = powerNet1.max()
		825	noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
		826	mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
		827	Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
		828	Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
		829	PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
		830	#Vectorize
		831	meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
		832	meteorAux[7:11] = phaseDiffint[0:4]
		833
		834	#Rejection Criterions
		835	if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
		836	meteorAux[-1] = 17
		837	elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
		838	meteorAux[-1] = 1
		839
		840
		841	else:
		842	meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
		843	meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
		844	PowerSeries = 0
		845
		846	listMeteors1.append(meteorAux)
		847	listPowerSeries.append(PowerSeries)
		848	listVoltageSeries.append(meteorVolts1)
		849
		850	return listMeteors1, listPowerSeries, listVoltageSeries
		851
		852	def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
		853
		854	threshError = 10
		855	#Depending if it is 30 or 50 MHz
		856	if frequency == 30e6:
		857	timeLag = 4510*-3
		858	else:
		859	timeLag = 1510*-3
		860	lag = numpy.ceil(timeLag/timeInterval)
		861
		862	listMeteors1 = []
		863
		864	for i in range(len(listMeteors)):
		865	meteorPower = listPower[i]
		866	meteorAux = listMeteors[i]
		867
		868	if meteorAux[-1] == 0:
		869
		870	try:
		871	indmax = meteorPower.argmax()
		872	indlag = indmax + lag
		873
		874	y = meteorPower[indlag:]
		875	x = numpy.arange(0, y.size)*timeLag
		876
		877	#first guess
		878	a = y[0]
		879	tau = timeLag
		880	#exponential fit
		881	popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
		882	y1 = self.__exponential_function(x, *popt)
		883	#error estimation
		884	error = sum((y - y1)*2)/(numpy.var(y)(y.size - popt.size))
		885
		886	decayTime = popt[1]
		887	riseTime = indmax*timeInterval
		888	meteorAux[11:13] = [decayTime, error]
		889
		890	#Table items 7, 8 and 11
		891	if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
		892	meteorAux[-1] = 7
		893	elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
		894	meteorAux[-1] = 8
		895	if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
		896	meteorAux[-1] = 11
		897
		898
		899	except:
		900	meteorAux[-1] = 11
		901
		902
		903	listMeteors1.append(meteorAux)
		904
		905	return listMeteors1
		906
		907	#Exponential Function
		908
		909	def __exponential_function(self, x, a, tau):
		910	y = a*numpy.exp(-x/tau)
		911	return y
		912
		913	def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
		914
		915	pairslist1 = list(pairslist)
		916	pairslist1.append((0,1))
		917	pairslist1.append((3,4))
		918	numPairs = len(pairslist1)
		919	#Time Lag
		920	timeLag = 4510*-3
		921	c = 3e8
		922	lag = numpy.ceil(timeLag/timeInterval)
		923	freq = 30e6
		924
		925	listMeteors1 = []
		926
		927	for i in range(len(listMeteors)):
		928	meteor = listMeteors[i]
		929	meteorAux = numpy.hstack((meteor[:-1], 0, 0, meteor[-1]))
		930	if meteor[-1] == 0:
		931	mStart = listMeteors[i][1]
		932	mPeak = listMeteors[i][2]
		933	mLag = mPeak - mStart + lag
		934
		935	#get the volt data between the start and end times of the meteor
		936	meteorVolts = listVolts[i]
		937	meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
		938
		939	#Get CCF
		940	allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
		941
		942	#Method 2
		943	slopes = numpy.zeros(numPairs)
		944	time = numpy.array([-2,-1,1,2])*timeInterval
		945	angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
		946
		947	#Correct phases
		948	derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
		949	indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
		950
		951	if indDer[0].shape[0] > 0:
		952	for i in range(indDer[0].shape[0]):
		953	signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
		954	angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo2numpy.pi
		955
		956	# fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
		957	for j in range(numPairs):
		958	fit = stats.linregress(time, angAllCCF[j,:])
		959	slopes[j] = fit[0]
		960
		961	#Remove Outlier
		962	# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
		963	# slopes = numpy.delete(slopes,indOut)
		964	# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
		965	# slopes = numpy.delete(slopes,indOut)
		966
		967	radialVelocity = -numpy.mean(slopes)(0.25/numpy.pi)(c/freq)
		968	radialError = numpy.std(slopes)(0.25/numpy.pi)(c/freq)
		969	meteorAux[-2] = radialError
		970	meteorAux[-3] = radialVelocity
		971
		972	#Setting Error
		973	#Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
		974	if numpy.abs(radialVelocity) > 200:
		975	meteorAux[-1] = 15
		976	#Number 12: Poor fit to CCF variation for estimation of radial drift velocity
		977	elif radialError > radialStdThresh:
		978	meteorAux[-1] = 12
		979
		980	listMeteors1.append(meteorAux)
		981	return listMeteors1
		982
		983	def __setNewArrays(self, listMeteors, date, heiRang):
		984
		985	#New arrays
		986	arrayMeteors = numpy.array(listMeteors)
		987	arrayParameters = numpy.zeros((len(listMeteors),10))
		988
		989	#Date inclusion
		990	date = re.findall(r'\((.*?)\)', date)
		991	date = date[0].split(',')
		992	date = map(int, date)
		993	date = [date[0]10000 + date[1]100 + date[2], date[3]10000 + date[4]100 + date[5]]
		994	arrayDate = numpy.tile(date, (len(listMeteors), 1))
		995
		996	#Meteor array
		997	arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
		998	arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
		999
		1000	#Parameters Array
		1001	arrayParameters[:,0:3] = arrayMeteors[:,0:3]
		1002	arrayParameters[:,-3:] = arrayMeteors[:,-3:]
		1003
		1004	return arrayMeteors, arrayParameters
		1005
		1006	def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
		1007
		1008	arrayAOA = numpy.zeros((phases.shape[0],3))
		1009	cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
		1010
		1011	arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
		1012	cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
		1013	arrayAOA[:,2] = cosDirError
		1014
		1015	azimuthAngle = arrayAOA[:,0]
		1016	zenithAngle = arrayAOA[:,1]
		1017
		1018	#Setting Error
		1019	#Number 3: AOA not fesible
		1020	indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
		1021	error[indInvalid] = 3
		1022	#Number 4: Large difference in AOAs obtained from different antenna baselines
		1023	indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
		1024	error[indInvalid] = 4
		1025	return arrayAOA, error
		1026
		1027	def __getDirectionCosines(self, arrayPhase, pairsList):
		1028
		1029	#Initializing some variables
		1030	ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])2numpy.pi
		1031	ang_aux = ang_aux.reshape(1,ang_aux.size)
		1032
		1033	cosdir = numpy.zeros((arrayPhase.shape[0],2))
		1034	cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
		1035
		1036
		1037	for i in range(2):
		1038	#First Estimation
		1039	phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
		1040	#Dealias
		1041	indcsi = numpy.where(phi0_aux > numpy.pi)
		1042	phi0_aux[indcsi] -= 2*numpy.pi
		1043	indcsi = numpy.where(phi0_aux < -numpy.pi)
		1044	phi0_aux[indcsi] += 2*numpy.pi
		1045	#Direction Cosine 0
		1046	cosdir0[:,i] = -(phi0_aux)/(2numpy.pi0.5)
		1047
		1048	#Most-Accurate Second Estimation
		1049	phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
		1050	phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
		1051	#Direction Cosine 1
		1052	cosdir1 = -(phi1_aux + ang_aux)/(2numpy.pi4.5)
		1053
		1054	#Searching the correct Direction Cosine
		1055	cosdir0_aux = cosdir0[:,i]
		1056	cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
		1057	#Minimum Distance
		1058	cosDiff = (cosdir1 - cosdir0_aux)**2
		1059	indcos = cosDiff.argmin(axis = 1)
		1060	#Saving Value obtained
		1061	cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
		1062
		1063	return cosdir0, cosdir
		1064
		1065	def __calculateAOA(self, cosdir, azimuth):
		1066	cosdirX = cosdir[:,0]
		1067	cosdirY = cosdir[:,1]
		1068
		1069	zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX2 - cosdirY2))*180/numpy.pi
		1070	azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
		1071	angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
		1072
		1073	return angles
		1074
		1075	def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
		1076
		1077	Ramb = 375 #Ramb = c/(2*PRF)
		1078	Re = 6371 #Earth Radius
		1079	heights = numpy.zeros(Ranges.shape)
		1080
		1081	R_aux = numpy.array([0,1,2])*Ramb
		1082	R_aux = R_aux.reshape(1,R_aux.size)
		1083
		1084	Ranges = Ranges.reshape(Ranges.size,1)
		1085
		1086	Ri = Ranges + R_aux
		1087	hi = numpy.sqrt(Re2 + Ri2 + (2Renumpy.cos(zenithnumpy.pi/180)Ri.transpose()).transpose()) - Re
		1088
		1089	#Check if there is a height between 70 and 110 km
		1090	h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
		1091	ind_h = numpy.where(h_bool == 1)[0]
		1092
		1093	hCorr = hi[ind_h, :]
		1094	ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
		1095
		1096	hCorr = hi[ind_hCorr]
		1097	heights[ind_h] = hCorr
		1098
		1099	#Setting Error
		1100	#Number 13: Height unresolvable echo: not valid height within 70 to 110 km
		1101	#Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
		1102
		1103	indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
		1104	error[indInvalid2] = 14
		1105	indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
		1106	error[indInvalid1] = 13
		1107
		1108	return heights, error
		1109
		1110
		1111	class WindProfiler(Operation):
		1112
		1113	__isConfig = False
		1114
		1115	__initime = None
		1116	__lastdatatime = None
		1117	__integrationtime = None
		1118
		1119	__buffer = None
		1120
		1121	__dataReady = False
		1122
		1123	__firstdata = None
		1124
		1125	n = None
		1126
		1127	def __init__(self):
		1128	Operation.__init__(self)
		1129
		1130	def __calculateAngles(self, theta_x, theta_y, azimuth):
		1131
		1132	dir_cosw = numpy.sqrt(1-theta_x2-theta_y2)
		1133	zenith_arr = numpy.arccos(dir_cosw)
		1134	azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
		1135
		1136	dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
		1137	dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
		1138
		1139	return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
		1140
		1141	def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
		1142
		1143	#
		1144	if horOnly:
		1145	A = numpy.c_[dir_cosu,dir_cosv]
		1146	else:
		1147	A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
		1148	A = numpy.asmatrix(A)
		1149	A1 = numpy.linalg.inv(A.transpose()A)A.transpose()
		1150
		1151	return A1
		1152
		1153	def __correctValues(self, heiRang, phi, velRadial, SNR):
		1154	listPhi = phi.tolist()
		1155	maxid = listPhi.index(max(listPhi))
		1156	minid = listPhi.index(min(listPhi))
		1157
		1158	rango = range(len(phi))
		1159	# rango = numpy.delete(rango,maxid)
		1160
		1161	heiRang1 = heiRang*math.cos(phi[maxid])
		1162	heiRangAux = heiRang*math.cos(phi[minid])
		1163	indOut = (heiRang1 < heiRangAux[0]).nonzero()
		1164	heiRang1 = numpy.delete(heiRang1,indOut)
		1165
		1166	velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
		1167	SNR1 = numpy.zeros([len(phi),len(heiRang1)])
		1168
		1169	for i in rango:
		1170	x = heiRang*math.cos(phi[i])
		1171	y1 = velRadial[i,:]
		1172	f1 = interpolate.interp1d(x,y1,kind = 'cubic')
		1173
		1174	x1 = heiRang1
		1175	y11 = f1(x1)
		1176
		1177	y2 = SNR[i,:]
		1178	f2 = interpolate.interp1d(x,y2,kind = 'cubic')
		1179	y21 = f2(x1)
		1180
		1181	velRadial1[i,:] = y11
		1182	SNR1[i,:] = y21
		1183
		1184	return heiRang1, velRadial1, SNR1
		1185
		1186	def __calculateVelUVW(self, A, velRadial):
		1187
		1188	#Operacion Matricial
		1189	# velUVW = numpy.zeros((velRadial.shape[1],3))
		1190	# for ind in range(velRadial.shape[1]):
		1191	# velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
		1192	# velUVW = velUVW.transpose()
		1193	velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
		1194	velUVW[:,:] = numpy.dot(A,velRadial)
		1195
		1196
		1197	return velUVW
		1198
		1199	def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
		1200	"""
		1201	Function that implements Doppler Beam Swinging (DBS) technique.
		1202
		1203	Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
		1204	Direction correction (if necessary), Ranges and SNR
		1205
		1206	Output: Winds estimation (Zonal, Meridional and Vertical)
		1207
		1208	Parameters affected: Winds, height range, SNR
		1209	"""
		1210	azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth)
		1211	heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0)
		1212	A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
		1213
		1214	#Calculo de Componentes de la velocidad con DBS
		1215	winds = self.__calculateVelUVW(A,velRadial1)
		1216
		1217	return winds, heiRang1, SNR1
		1218
		1219	def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):
		1220
		1221	posx = numpy.asarray(posx)
		1222	posy = numpy.asarray(posy)
		1223
		1224	#Rotacion Inversa para alinear con el azimuth
		1225	if azimuth!= None:
		1226	azimuth = azimuth*math.pi/180
		1227	posx1 = posxmath.cos(azimuth) + posymath.sin(azimuth)
		1228	posy1 = -posxmath.sin(azimuth) + posymath.cos(azimuth)
		1229	else:
		1230	posx1 = posx
		1231	posy1 = posy
		1232
		1233	#Calculo de Distancias
		1234	distx = numpy.zeros(pairsCrossCorr.size)
		1235	disty = numpy.zeros(pairsCrossCorr.size)
		1236	dist = numpy.zeros(pairsCrossCorr.size)
		1237	ang = numpy.zeros(pairsCrossCorr.size)
		1238
		1239	for i in range(pairsCrossCorr.size):
		1240	distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]
		1241	disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]
		1242	dist[i] = numpy.sqrt(distx[i]2 + disty[i]2)
		1243	ang[i] = numpy.arctan2(disty[i],distx[i])
		1244	#Calculo de Matrices
		1245	nPairs = len(pairs)
		1246	ang1 = numpy.zeros((nPairs, 2, 1))
		1247	dist1 = numpy.zeros((nPairs, 2, 1))
		1248
		1249	for j in range(nPairs):
		1250	dist1[j,0,0] = dist[pairs[j][0]]
		1251	dist1[j,1,0] = dist[pairs[j][1]]
		1252	ang1[j,0,0] = ang[pairs[j][0]]
		1253	ang1[j,1,0] = ang[pairs[j][1]]
		1254
		1255	return distx,disty, dist1,ang1
		1256
		1257	def __calculateVelVer(self, phase, lagTRange, _lambda):
		1258
		1259	Ts = lagTRange[1] - lagTRange[0]
		1260	velW = -_lambdaphase/(4math.pi*Ts)
		1261
		1262	return velW
		1263
		1264	def __calculateVelHorDir(self, dist, tau1, tau2, ang):
		1265	nPairs = tau1.shape[0]
		1266	vel = numpy.zeros((nPairs,3,tau1.shape[2]))
		1267
		1268	angCos = numpy.cos(ang)
		1269	angSin = numpy.sin(ang)
		1270
		1271	vel0 = disttau1/(2tau2**2)
		1272	vel[:,0,:] = (vel0*angCos).sum(axis = 1)
		1273	vel[:,1,:] = (vel0*angSin).sum(axis = 1)
		1274
		1275	ind = numpy.where(numpy.isinf(vel))
		1276	vel[ind] = numpy.nan
		1277
		1278	return vel
		1279
		1280	def __getPairsAutoCorr(self, pairsList, nChannels):
		1281
		1282	pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
		1283
		1284	for l in range(len(pairsList)):
		1285	firstChannel = pairsList[l][0]
		1286	secondChannel = pairsList[l][1]
		1287
		1288	#Obteniendo pares de Autocorrelacion
		1289	if firstChannel == secondChannel:
		1290	pairsAutoCorr[firstChannel] = int(l)
		1291
		1292	pairsAutoCorr = pairsAutoCorr.astype(int)
		1293
		1294	pairsCrossCorr = range(len(pairsList))
		1295	pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
		1296
		1297	return pairsAutoCorr, pairsCrossCorr
		1298
		1299	def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
		1300	"""
		1301	Function that implements Spaced Antenna (SA) technique.
		1302
		1303	Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
		1304	Direction correction (if necessary), Ranges and SNR
		1305
		1306	Output: Winds estimation (Zonal, Meridional and Vertical)
		1307
		1308	Parameters affected: Winds
		1309	"""
		1310	#Cross Correlation pairs obtained
		1311	pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
		1312	pairsArray = numpy.array(pairsList)[pairsCrossCorr]
		1313	pairsSelArray = numpy.array(pairsSelected)
		1314	pairs = []
		1315
		1316	#Wind estimation pairs obtained
		1317	for i in range(pairsSelArray.shape[0]/2):
		1318	ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
		1319	ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
		1320	pairs.append((ind1,ind2))
		1321
		1322	indtau = tau.shape[0]/2
		1323	tau1 = tau[:indtau,:]
		1324	tau2 = tau[indtau:-1,:]
		1325	tau1 = tau1[pairs,:]
		1326	tau2 = tau2[pairs,:]
		1327	phase1 = tau[-1,:]
		1328
		1329	#---------------------------------------------------------------------
		1330	#Metodo Directo
		1331	distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)
		1332	winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
		1333	winds = stats.nanmean(winds, axis=0)
		1334	#---------------------------------------------------------------------
		1335	#Metodo General
		1336	# distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
		1337	# #Calculo Coeficientes de Funcion de Correlacion
		1338	# F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
		1339	# #Calculo de Velocidades
		1340	# winds = self.calculateVelUV(F,G,A,B,H)
		1341
		1342	#---------------------------------------------------------------------
		1343	winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
		1344	winds = correctFactor*winds
		1345	return winds
		1346
		1347	def __checkTime(self, currentTime, paramInterval, windsInterval):
		1348
		1349	dataTime = currentTime + paramInterval
		1350	deltaTime = dataTime - self.__initime
		1351
		1352	if deltaTime >= windsInterval or deltaTime < 0:
		1353	self.__dataReady = True
		1354	return
		1355
		1356	def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
		1357	'''
		1358	Function that implements winds estimation technique with detected meteors.
		1359
		1360	Input: Detected meteors, Minimum meteor quantity to wind estimation
		1361
		1362	Output: Winds estimation (Zonal and Meridional)
		1363
		1364	Parameters affected: Winds
		1365	'''
		1366	#Settings
		1367	nInt = (heightMax - heightMin)/2
		1368	winds = numpy.zeros((2,nInt))*numpy.nan
		1369
		1370	#Filter errors
		1371	error = numpy.where(arrayMeteor[:,-1] == 0)[0]
		1372	finalMeteor = arrayMeteor[error,:]
		1373
		1374	#Meteor Histogram
		1375	finalHeights = finalMeteor[:,3]
		1376	hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
		1377	nMeteorsPerI = hist[0]
		1378	heightPerI = hist[1]
		1379
		1380	#Sort of meteors
		1381	indSort = finalHeights.argsort()
		1382	finalMeteor2 = finalMeteor[indSort,:]
		1383
		1384	# Calculating winds
		1385	ind1 = 0
		1386	ind2 = 0
		1387
		1388	for i in range(nInt):
		1389	nMet = nMeteorsPerI[i]
		1390	ind1 = ind2
		1391	ind2 = ind1 + nMet
		1392
		1393	meteorAux = finalMeteor2[ind1:ind2,:]
		1394
		1395	if meteorAux.shape[0] >= meteorThresh:
		1396	vel = meteorAux[:, 7]
		1397	zen = meteorAux[:, 5]*numpy.pi/180
		1398	azim = meteorAux[:, 4]*numpy.pi/180
		1399
		1400	n = numpy.cos(zen)
		1401	# m = (1 - n2)/(1 - numpy.tan(azim)2)
		1402	# l = m*numpy.tan(azim)
		1403	l = numpy.sin(zen)*numpy.sin(azim)
		1404	m = numpy.sin(zen)*numpy.cos(azim)
		1405
		1406	A = numpy.vstack((l, m)).transpose()
		1407	A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
		1408	windsAux = numpy.dot(A1, vel)
		1409
		1410	winds[0,i] = windsAux[0]
		1411	winds[1,i] = windsAux[1]
		1412
		1413	return winds, heightPerI[:-1]
		1414
		1415	def run(self, dataOut, technique, **kwargs):
		1416
		1417	param = dataOut.data_param
		1418	if dataOut.abscissaRange != None:
		1419	absc = dataOut.abscissaRange[:-1]
		1420	noise = dataOut.noise
		1421	heightRange = dataOut.getHeiRange()
		1422	SNR = dataOut.SNR
		1423
		1424	if technique == 'DBS':
		1425
		1426	theta_x = numpy.array(kwargs['dirCosx'])
		1427	theta_y = numpy.array(kwargs['dirCosy'])
		1428	azimuth = kwargs['azimuth']
		1429	if kwargs.has_key('horizontalOnly'):
		1430	horizontalOnly = kwargs['horizontalOnly']
		1431	else: horizontalOnly = False
		1432	if kwargs.has_key('correctFactor'):
		1433	correctFactor = kwargs['correctFactor']
		1434	else: correctFactor = 1
		1435	if kwargs.has_key('channelList'):
		1436	channelList = kwargs['channelList']
		1437	if len(channelList) == 2:
		1438	horizontalOnly = True
		1439	arrayChannel = numpy.array(channelList)
		1440	param = param[arrayChannel,:,:]
		1441	theta_x = theta_x[arrayChannel]
		1442	theta_y = theta_y[arrayChannel]
		1443
		1444	velRadial0 = param[:,1,:] #Radial velocity
		1445	dataOut.winds, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function
		1446	dataOut.initUtcTime = dataOut.ltctime
		1447	dataOut.windsInterval = dataOut.timeInterval
		1448
		1449	elif technique == 'SA':
		1450
		1451	#Parameters
		1452	position_x = kwargs['positionX']
		1453	position_y = kwargs['positionY']
		1454	azimuth = kwargs['azimuth']
		1455
		1456	if kwargs.has_key('crosspairsList'):
		1457	pairs = kwargs['crosspairsList']
		1458	else:
		1459	pairs = None
		1460
		1461	if kwargs.has_key('correctFactor'):
		1462	correctFactor = kwargs['correctFactor']
		1463	else:
		1464	correctFactor = 1
		1465
		1466	tau = dataOut.data_param
		1467	_lambda = dataOut.C/dataOut.frequency
		1468	pairsList = dataOut.pairsList
		1469	nChannels = dataOut.nChannels
		1470
		1471	dataOut.winds = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
		1472	dataOut.initUtcTime = dataOut.ltctime
		1473	dataOut.windsInterval = dataOut.timeInterval
		1474
		1475	elif technique == 'Meteors':
		1476	dataOut.flagNoData = True
		1477	self.__dataReady = False
		1478
		1479	if kwargs.has_key('nHours'):
		1480	nHours = kwargs['nHours']
		1481	else:
		1482	nHours = 1
		1483
		1484	if kwargs.has_key('meteorsPerBin'):
		1485	meteorThresh = kwargs['meteorsPerBin']
		1486	else:
		1487	meteorThresh = 6
		1488
		1489	if kwargs.has_key('hmin'):
		1490	hmin = kwargs['hmin']
		1491	else: hmin = 70
		1492	if kwargs.has_key('hmax'):
		1493	hmax = kwargs['hmax']
		1494	else: hmax = 110
		1495
		1496	dataOut.windsInterval = nHours*3600
		1497
		1498	if self.__isConfig == False:
		1499	# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
		1500	#Get Initial LTC time
		1501	self.__initime = (dataOut.datatime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
		1502	self.__isConfig = True
		1503
		1504	if self.__buffer == None:
		1505	self.__buffer = dataOut.data_param
		1506	self.__firstdata = copy.copy(dataOut)
		1507
		1508	else:
		1509	self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
		1510
		1511	self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.windsInterval) #Check if the buffer is ready
		1512
		1513	if self.__dataReady:
		1514	dataOut.initUtcTime = self.__initime
		1515	self.__initime = self.__initime + dataOut.windsInterval #to erase time offset
		1516
		1517	dataOut.winds, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
		1518	dataOut.flagNoData = False
		1519	self.__buffer = None
		1520
		1521	return No newline at end of file

schainpy/test/Meteor_JASMET30MHz_Winds.py created 644 +103 0

@@ -0,0 +1,103
	1	# DIAS 19 Y 20 FEB 2014
	2	# Comprobacion de Resultados DBS con SA
	3
	4	import os, sys
	5
	6	path = os.path.split(os.getcwd())[0]
	7	sys.path.append(path)
	8
	9	from controller import *
	10
	11	desc = "JASMET Experiment Test"
	12	filename = "JASMETtest.xml"
	13
	14	controllerObj = Project()
	15
	16	controllerObj.setup(id = '191', name='test01', description=desc)
	17
	18	#Experimentos
	19
	20	#2014051 20 Feb 2014
	21	path = '/home/soporte/Data/JASMET/JASMET_30/2014106'
	22	pathFigure = '/home/soporte/workspace/Graficos/JASMET/prueba1'
	23
	24	startTime = '00:00:00'
	25	endTime = '23:59:59'
	26	xmin ='19.0'
	27	xmax = '34.0'
	28
	29	#------------------------------------------------------------------------------------------------
	30	readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
	31	path=path,
	32	startDate='2014/04/15',
	33	endDate='2014/04/16',
	34	startTime=startTime,
	35	endTime=endTime,
	36	online=0,
	37	delay=5,
	38	walk=0)
	39
	40	opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
	41
	42
	43	#--------------------------------------------------------------------------------------------------
	44
	45	procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
	46
	47	opObj00 = procUnitConfObj0.addOperation(name='selectChannels')
	48	opObj00.addParameter(name='channelList', value='0, 1, 2, 3, 4', format='intlist')
	49
	50	opObj01 = procUnitConfObj0.addOperation(name='setRadarFrequency')
	51	opObj01.addParameter(name='frequency', value='30.e6', format='float')
	52
	53	#opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
	54	#--------------------------------------------------------------------------------------------------
	55
	56	procUnitConfObj1 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj0.getId())
	57	procUnitConfObj1.addParameter(name='nSeconds', value='100', format='int')
	58
	59	opObj10 = procUnitConfObj1.addOperation(name='DetectMeteors')
	60	opObj10.addParameter(name='predefinedPhaseShifts', value='-89.5, 41.5, 0.0, -138.0, -85.5', format='floatlist')
	61	opObj10.addParameter(name='cohDetection', value='0', format='bool')
	62	opObj10.addParameter(name='noise_multiple', value='4', format='int')
	63	opObj10.addParameter(name='SNRThresh', value='5', format='float')
	64	opObj10.addParameter(name='phaseThresh', value='20', format='float')
	65	opObj10.addParameter(name='azimuth', value='45', format='float')
	66	opObj10.addParameter(name='hmin', value='68', format='float')
	67	opObj10.addParameter(name='hmax', value='112', format='float')
	68
	69	opObj13 = procUnitConfObj1.addOperation(name='SkyMapPlot', optype='other')
	70	opObj13.addParameter(name='id', value='1', format='int')
	71	opObj13.addParameter(name='wintitle', value='Sky Map', format='str')
	72	opObj13.addParameter(name='save', value='1', format='bool')
	73	opObj13.addParameter(name='figpath', value=pathFigure, format='str')
	74
	75	opObj11 = procUnitConfObj1.addOperation(name='WindProfiler', optype='other')
	76	opObj11.addParameter(name='technique', value='Meteors', format='str')
	77	opObj11.addParameter(name='nHours', value='1.0', format='float')
	78	opObj11.addParameter(name='SNRThresh', value='12.0', format='float')
	79
	80	opObj12 = procUnitConfObj1.addOperation(name='WindProfilerPlot', optype='other')
	81	opObj12.addParameter(name='id', value='2', format='int')
	82	opObj12.addParameter(name='wintitle', value='Wind Profiler', format='str')
	83	opObj12.addParameter(name='save', value='1', format='bool')
	84	opObj12.addParameter(name='figpath', value = pathFigure, format='str')
	85	opObj12.addParameter(name='zmin', value='-120', format='int')
	86	opObj12.addParameter(name='zmax', value='120', format='int')
	87	# opObj12.addParameter(name='zmin_ver', value='-0.8', format='float')
	88	# opObj12.addParameter(name='zmax_ver', value='0.8', format='float')
	89	# opObj23.addParameter(name='SNRmin', value='-10', format='int')
	90	# opObj23.addParameter(name='SNRmax', value='60', format='int')
	91	# opObj23.addParameter(name='SNRthresh', value='0', format='float')
	92	opObj12.addParameter(name='xmin', value=xmin, format='float')
	93	opObj12.addParameter(name='xmax', value=xmax, format='float')
	94
	95	#--------------------------------------------------------------------------------------------------
	96	print "Escribiendo el archivo XML"
	97	controllerObj.writeXml(filename)
	98	print "Leyendo el archivo XML"
	99	controllerObj.readXml(filename)
	100
	101	controllerObj.createObjects()
	102	controllerObj.connectObjects()
	103	controllerObj.run() No newline at end of file

schainpy/test/Meteor_JASMET50MHz_Winds.py created 644 +105 0

@@ -0,0 +1,105
	1	# DIAS 19 Y 20 FEB 2014
	2	# Comprobacion de Resultados DBS con SA
	3
	4	import os, sys
	5
	6	path = os.path.split(os.getcwd())[0]
	7	sys.path.append(path)
	8
	9	from controller import *
	10
	11	desc = "JASMET Experiment Test"
	12	filename = "JASMETtest.xml"
	13
	14	controllerObj = Project()
	15
	16	controllerObj.setup(id = '191', name='test01', description=desc)
	17
	18	#Experimentos
	19
	20	#2014051 20 Feb 2014
	21	path = '/home/soporte/Data/JASMET/JASMET_50/2014106'
	22	pathFigure = '/home/soporte/workspace/Graficos/JASMET'
	23
	24	startTime = '00:00:00'
	25	endTime = '23:59:59'
	26	xmin ='19.0'
	27	xmax = '33.0'
	28
	29
	30	#------------------------------------------------------------------------------------------------
	31	readUnitConfObj = controllerObj.addReadUnit(datatype='Voltage',
	32	path=path,
	33	startDate='2014/04/15',
	34	endDate='2014/04/16',
	35	startTime=startTime,
	36	endTime=endTime,
	37	online=0,
	38	delay=5,
	39	walk=0)
	40
	41	opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
	42
	43
	44	#--------------------------------------------------------------------------------------------------
	45
	46	procUnitConfObj0 = controllerObj.addProcUnit(datatype='Voltage', inputId=readUnitConfObj.getId())
	47
	48	opObj00 = procUnitConfObj0.addOperation(name='selectChannels')
	49	opObj00.addParameter(name='channelList', value='0, 1, 2, 3, 4', format='intlist')
	50
	51	opObj01 = procUnitConfObj0.addOperation(name='setRadarFrequency')
	52	opObj01.addParameter(name='frequency', value='50.e6', format='float')
	53
	54	#opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
	55	#--------------------------------------------------------------------------------------------------
	56
	57	procUnitConfObj1 = controllerObj.addProcUnit(datatype='Parameters', inputId=procUnitConfObj0.getId())
	58	procUnitConfObj1.addParameter(name='nSeconds', value='100', format='int')
	59
	60	opObj10 = procUnitConfObj1.addOperation(name='DetectMeteors')
	61	opObj10.addParameter(name='predefinedPhaseShifts', value='-17.0, -139.0, 0.0, 48.0, -130.0', format='floatlist')
	62	opObj10.addParameter(name='cohDetection', value='0', format='bool')
	63	opObj10.addParameter(name='noise_multiple', value='4', format='int')
	64	opObj10.addParameter(name='SNRThresh', value='5', format='float')
	65	opObj10.addParameter(name='phaseThresh', value='20', format='float')
	66	opObj10.addParameter(name='azimuth', value='45', format='float')
	67	opObj10.addParameter(name='hmin', value='68', format='float')
	68	opObj10.addParameter(name='hmax', value='112', format='float')
	69	opObj10.addParameter(name='savefile', value='1', format='bool')
	70	opObj10.addParameter(name='filename', value=filehdf5, format='str')
	71
	72	opObj13 = procUnitConfObj1.addOperation(name='SkyMapPlot', optype='other')
	73	opObj13.addParameter(name='id', value='1', format='int')
	74	opObj13.addParameter(name='wintitle', value='Sky Map', format='str')
	75	opObj13.addParameter(name='save', value='1', format='bool')
	76	opObj13.addParameter(name='figpath', value=pathFigure, format='str')
	77
	78	opObj11 = procUnitConfObj1.addOperation(name='WindProfiler', optype='other')
	79	opObj11.addParameter(name='technique', value='Meteors', format='str')
	80	opObj11.addParameter(name='nHours', value='1.0', format='float')
	81
	82	opObj12 = procUnitConfObj1.addOperation(name='WindProfilerPlot', optype='other')
	83	opObj12.addParameter(name='id', value='2', format='int')
	84	opObj12.addParameter(name='wintitle', value='Wind Profiler', format='str')
	85	opObj12.addParameter(name='save', value='1', format='bool')
	86	opObj12.addParameter(name='figpath', value = pathFigure, format='str')
	87	opObj12.addParameter(name='zmin', value='-120', format='int')
	88	opObj12.addParameter(name='zmax', value='120', format='int')
	89	# opObj12.addParameter(name='zmin_ver', value='-0.8', format='float')
	90	# opObj12.addParameter(name='zmax_ver', value='0.8', format='float')
	91	# opObj23.addParameter(name='SNRmin', value='-10', format='int')
	92	# opObj23.addParameter(name='SNRmax', value='60', format='int')
	93	# opObj23.addParameter(name='SNRthresh', value='0', format='float')
	94	opObj12.addParameter(name='xmin', value=xmin, format='float')
	95	opObj12.addParameter(name='xmax', value=xmax, format='float')
	96
	97	#--------------------------------------------------------------------------------------------------
	98	print "Escribiendo el archivo XML"
	99	controllerObj.writeXml(filename)
	100	print "Leyendo el archivo XML"
	101	controllerObj.readXml(filename)
	102
	103	controllerObj.createObjects()
	104	controllerObj.connectObjects()
	105	controllerObj.run() No newline at end of file

schainpy/test/WindProfiler_DBS01.py created 644 +144 0

@@ -0,0 +1,144
	1	# DIAS 19 Y 20 FEB 2014
	2	# Comprobacion de Resultados DBS con SA
	3
	4	import os, sys
	5
	6	path = os.path.split(os.getcwd())[0]
	7	sys.path.append(path)
	8
	9	from controller import *
	10
	11	desc = "DBS Experiment Test"
	12	filename = "DBStest.xml"
	13
	14	controllerObj = Project()
	15
	16	controllerObj.setup(id = '191', name='test01', description=desc)
	17
	18	#Experimentos
	19
	20	#2014050 19 Feb 2014
	21	# path = '/home/soporte/Documents/MST_Data/DBS/d2014050'
	22	# pathFigure = '/home/soporte/workspace/Graficos/DBS/d2014050p/'
	23	# xmin = '15.5'
	24	# xmax = '23.99999999'
	25	# startTime = '17:25:00'
	26	# filehdf5 = "DBS_2014050.hdf5"
	27
	28	#2014051 20 Feb 2014
	29	path = '/home/soporte/Data/MST/DBS/d2014051'
	30	pathFigure = '/home/soporte/workspace/Graficos/DBS/prueba1/'
	31	xmin = '0.0'
	32	xmax = '8.0'
	33	startTime = '00:00:00'
	34	filehdf5 = "DBS_2014051.hdf5"
	35
	36
	37
	38	#------------------------------------------------------------------------------------------------
	39	readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
	40	path=path,
	41	startDate='2014/01/31',
	42	endDate='2014/03/31',
	43	startTime=startTime,
	44	endTime='23:59:59',
	45	online=0,
	46	delay=5,
	47	walk=0)
	48
	49	opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
	50
	51
	52	#--------------------------------------------------------------------------------------------------
	53
	54	procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
	55
	56	opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
	57
	58	opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other')
	59	opObj11.addParameter(name='n', value='256', format='int')
	60	# opObj11.addParameter(name='n', value='16', format='int')
	61
	62	opObj11 = procUnitConfObj0.addOperation(name='selectHeightsByIndex')
	63	opObj11.addParameter(name='minIndex', value='10', format='float')
	64	opObj11.addParameter(name='maxIndex', value='60', format='float')
	65
	66	#---------------------------------------------------------------------------------------------------
	67
	68	procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId())
	69	procUnitConfObj1.addParameter(name='nFFTPoints', value='64', format='int')
	70	procUnitConfObj1.addParameter(name='nProfiles', value='64', format='int')
	71	# procUnitConfObj1.addParameter(name='ippFactor', value='2', format='int')
	72	procUnitConfObj1.addParameter(name='pairsList', value='(0,0),(0,1),(2,1)', format='pairsList')
	73
	74	opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
	75	opObj11.addParameter(name='n', value='5', format='int')
	76
	77	opObj14 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
	78	opObj14.addParameter(name='id', value='1', format='int')
	79	opObj14.addParameter(name='wintitle', value='Con interf', format='str')
	80	opObj14.addParameter(name='save', value='1', format='bool')
	81	opObj14.addParameter(name='figpath', value=pathFigure, format='str')
	82	opObj14.addParameter(name='zmin', value='5', format='int')
	83	opObj14.addParameter(name='zmax', value='90', format='int')
	84
	85	opObj12 = procUnitConfObj1.addOperation(name='removeInterference')
	86	opObj13 = procUnitConfObj1.addOperation(name='removeDC')
	87	opObj13.addParameter(name='mode', value='1', format='int')
	88
	89	opObj12 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
	90	opObj12.addParameter(name='id', value='2', format='int')
	91	opObj12.addParameter(name='wintitle', value='RTI Plot', format='str')
	92	opObj12.addParameter(name='save', value='1', format='bool')
	93	opObj12.addParameter(name='figpath', value = pathFigure, format='str')
	94	opObj12.addParameter(name='xmin', value=xmin, format='float')
	95	opObj12.addParameter(name='xmax', value=xmax, format='float')
	96	opObj12.addParameter(name='zmin', value='5', format='int')
	97	opObj12.addParameter(name='zmax', value='90', format='int')
	98
	99	#--------------------------------------------------------------------------------------------------
	100
	101	procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
	102	opObj20 = procUnitConfObj2.addOperation(name='GetMoments')
	103
	104	opObj21 = procUnitConfObj2.addOperation(name='MomentsPlot', optype='other')
	105	opObj21.addParameter(name='id', value='3', format='int')
	106	opObj21.addParameter(name='wintitle', value='Moments Plot', format='str')
	107	opObj21.addParameter(name='save', value='1', format='bool')
	108	opObj21.addParameter(name='figpath', value=pathFigure, format='str')
	109	opObj21.addParameter(name='zmin', value='5', format='int')
	110	opObj21.addParameter(name='zmax', value='90', format='int')
	111
	112	opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other')
	113	opObj22.addParameter(name='technique', value='DBS', format='str')
	114	opObj22.addParameter(name='azimuth', value='51.06', format='float')
	115	opObj22.addParameter(name='correctFactor', value='-1', format='float')
	116	opObj22.addParameter(name='dirCosx', value='0.041016, 0, -0.054688', format='floatlist')
	117	opObj22.addParameter(name='dirCosy', value='-0.041016, 0.025391, -0.023438', format='floatlist')
	118	# opObj22.addParameter(name='horizontalOnly', value='1', format='bool')
	119	# opObj22.addParameter(name='channelList', value='1,2', format='intlist')
	120
	121	opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
	122	opObj23.addParameter(name='id', value='4', format='int')
	123	opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str')
	124	opObj23.addParameter(name='save', value='1', format='bool')
	125	opObj23.addParameter(name='figpath', value = pathFigure, format='str')
	126	opObj23.addParameter(name='zmin', value='-10', format='int')
	127	opObj23.addParameter(name='zmax', value='10', format='int')
	128	opObj23.addParameter(name='zmin_ver', value='-80', format='float')
	129	opObj23.addParameter(name='zmax_ver', value='80', format='float')
	130	opObj23.addParameter(name='SNRmin', value='-10', format='int')
	131	opObj23.addParameter(name='SNRmax', value='60', format='int')
	132	opObj23.addParameter(name='SNRthresh', value='0', format='float')
	133	opObj23.addParameter(name='xmin', value=xmin, format='float')
	134	opObj23.addParameter(name='xmax', value=xmax, format='float')
	135
	136	#--------------------------------------------------------------------------------------------------
	137	print "Escribiendo el archivo XML"
	138	controllerObj.writeXml(filename)
	139	print "Leyendo el archivo XML"
	140	controllerObj.readXml(filename)
	141
	142	controllerObj.createObjects()
	143	controllerObj.connectObjects()
	144	controllerObj.run() No newline at end of file

schainpy/test/WindProfiler_SA01.py created 644 +139 0

@@ -0,0 +1,139
	1	# DIAS 19 Y 20 FEB 2014
	2	# Comprobacion de Resultados DBS con SA
	3
	4	import os, sys
	5
	6	path = os.path.split(os.getcwd())[0]
	7	sys.path.append(path)
	8
	9	from controller import *
	10
	11	desc = "SA Experiment Test"
	12	filename = "SA2014050.xml"
	13
	14	controllerObj = Project()
	15
	16	controllerObj.setup(id = '191', name='test01', description=desc)
	17
	18
	19	#Experimentos
	20
	21	#2014050 19 Feb 2014
	22	# path = '/home/soporte/Documents/MST_Data/SA/d2014050'
	23	# pathFigure = '/home/soporte/workspace/Graficos/SA/d2014050_prueba/'
	24	# xmin = '15.5'
	25	# xmax = '23.99999999'
	26	# startTime = '15:30:00'
	27	# filehdf5 = "SA_2014050.hdf5"
	28
	29	#2014051 20 Feb 2014
	30	path = '/home/soporte/Data/MST/SA/d2014051'
	31	pathFigure = '/home/soporte/workspace/Graficos/SA/prueba1/'
	32	xmin = '0.0'
	33	xmax = '8.0'
	34	startTime = '06:00:00'
	35	filehdf5 = "SA_2014051.hdf5"
	36
	37	readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
	38	path=path,
	39	startDate='2014/01/01',
	40	endDate='2014/03/31',
	41	startTime=startTime,
	42	endTime='23:59:59',
	43	online=0,
	44	delay=5,
	45	walk=0)
	46
	47	opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
	48
	49
	50	#--------------------------------------------------------------------------------------------------
	51
	52	procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
	53
	54	opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
	55
	56	opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other')
	57	opObj11.addParameter(name='n', value='600', format='int')
	58	# opObj11.addParameter(name='n', value='10', format='int')
	59
	60	opObj11 = procUnitConfObj0.addOperation(name='selectHeightsByIndex')
	61	opObj11.addParameter(name='minIndex', value='10', format='float')
	62	opObj11.addParameter(name='maxIndex', value='60', format='float')
	63	#---------------------------------------------------------------------------------------------------
	64	procUnitConfObj1 = controllerObj.addProcUnit(datatype='CorrelationProc', inputId=procUnitConfObj0.getId())
	65	# procUnitConfObj1.addParameter(name='pairsList', value='(0,0),(1,1),(2,2),(3,3),(1,0),(2,3)', format='pairsList')
	66	procUnitConfObj1.addParameter(name='pairsList', value='(0,0),(1,1),(2,2),(3,3),(0,3),(0,2),(1,3),(1,2)', format='pairsList')
	67	procUnitConfObj1.addParameter(name='fullT', value='1', format='bool')
	68	procUnitConfObj1.addParameter(name='removeDC', value='1', format='bool')
	69	#procUnitConfObj1.addParameter(name='lagT', value='0,1,2,3', format='intlist')
	70
	71	opObj12 = procUnitConfObj1.addOperation(name='CorrelationPlot', optype='other')
	72	opObj12.addParameter(name='id', value='1', format='int')
	73	opObj12.addParameter(name='wintitle', value='CrossCorrelation Plot', format='str')
	74	opObj12.addParameter(name='save', value='1', format='bool')
	75	opObj12.addParameter(name='zmin', value='0', format='int')
	76	opObj12.addParameter(name='zmax', value='1', format='int')
	77	opObj12.addParameter(name='figpath', value = pathFigure, format='str')
	78
	79	opObj12 = procUnitConfObj1.addOperation(name='removeNoise')
	80	opObj12.addParameter(name='mode', value='2', format='int')
	81	opObj12 = procUnitConfObj1.addOperation(name='calculateNormFactor')
	82
	83	opObj12 = procUnitConfObj1.addOperation(name='CorrelationPlot', optype='other')
	84	opObj12.addParameter(name='id', value='2', format='int')
	85	opObj12.addParameter(name='wintitle', value='CrossCorrelation Plot', format='str')
	86	opObj12.addParameter(name='save', value='1', format='bool')
	87	opObj12.addParameter(name='zmin', value='0', format='int')
	88	opObj12.addParameter(name='zmax', value='1', format='int')
	89	opObj12.addParameter(name='figpath', value = pathFigure, format='str')
	90
	91	#---------------------------------------------------------------------------------------------------
	92	procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
	93	opObj20 = procUnitConfObj2.addOperation(name='GetLags')
	94
	95	opObj21 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other')
	96	opObj21.addParameter(name='technique', value='SA', format='str')
	97	# opObj21.addParameter(name='correctFactor', value='-1', format='float')
	98	opObj21.addParameter(name='positionX', value='36,0,36,0', format='floatlist')
	99	opObj21.addParameter(name='positionY', value='36,0,0,36', format='floatlist')
	100	opObj21.addParameter(name='azimuth', value='51.06', format='float')
	101	opObj21.addParameter(name='crosspairsList', value='(0,3),(0,2),(1,3),(1,2)', format='pairsList')#COrregir
	102	#
	103	opObj22 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
	104	opObj22.addParameter(name='id', value='4', format='int')
	105	opObj22.addParameter(name='wintitle', value='Wind Profiler', format='str')
	106	opObj22.addParameter(name='save', value='1', format='bool')
	107	opObj22.addParameter(name='figpath', value = pathFigure, format='str')
	108	opObj22.addParameter(name='zmin', value='-15', format='int')
	109	opObj22.addParameter(name='zmax', value='15', format='int')
	110	opObj22.addParameter(name='zmin_ver', value='-80', format='float')
	111	opObj22.addParameter(name='zmax_ver', value='80', format='float')
	112	opObj22.addParameter(name='SNRmin', value='-20', format='int')
	113	opObj22.addParameter(name='SNRmax', value='40', format='int')
	114	opObj22.addParameter(name='SNRthresh', value='-3.5', format='float')
	115	opObj22.addParameter(name='xmin', value=xmin, format='float')
	116	opObj22.addParameter(name='xmax', value=xmax, format='float')
	117	# #-----------------------------------------------------------------------------------
	118	#
	119	# procUnitConfObj2 = controllerObj.addProcUnit(datatype='Spectra', inputId=procUnitConfObj0.getId())
	120	# procUnitConfObj2.addParameter(name='nFFTPoints', value='128', format='int')
	121	# procUnitConfObj2.addParameter(name='nProfiles', value='128', format='int')
	122	# procUnitConfObj2.addParameter(name='pairsList', value='(0,0),(0,1),(2,1)', format='pairsList')
	123	#
	124	# opObj22 = procUnitConfObj2.addOperation(name='SpectraPlot', optype='other')
	125	# opObj22.addParameter(name='id', value='5', format='int')
	126	# opObj22.addParameter(name='wintitle', value='Spectra Plot', format='str')
	127	# opObj22.addParameter(name='save', value='1', format='bool')
	128	# opObj22.addParameter(name='figpath', value = pathFigure, format='str')
	129
	130	#-----------------------------------------------------------------------------------
	131
	132	print "Escribiendo el archivo XML"
	133	controllerObj.writeXml(filename)
	134	print "Leyendo el archivo XML"
	135	controllerObj.readXml(filename)
	136
	137	controllerObj.createObjects()
	138	controllerObj.connectObjects()
	139	controllerObj.run() No newline at end of file

schainpy/test/amisr_windEstimation.py created 644 +131 0

@@ -0,0 +1,131
	1	import os, sys
	2
	3	path = os.path.split(os.getcwd())[0]
	4	sys.path.append(path)
	5
	6	from controller import *
	7
	8	desc = "AMISR Experiment"
	9
	10	filename = "amisr_reader.xml"
	11
	12	controllerObj = Project()
	13
	14	controllerObj.setup(id = '191', name='test01', description=desc)
	15
	16
	17	path = os.path.join(os.environ['HOME'],'Development/amisr/data')
	18	path = '/home/soporte/Data/AMISR'
	19	figpath = os.path.join(os.environ['HOME'],'Pictures/amisr')
	20
	21	pathFigure = '/home/soporte/workspace/Graficos/DBS/amisr/'
	22	xmin = '12.0'
	23	xmax = '14.0'
	24
	25	readUnitConfObj = controllerObj.addReadUnit(datatype='AMISRReader',
	26	path=path,
	27	startDate='2014/10/21',
	28	endDate='2014/10/21',
	29	startTime='00:00:00',
	30	endTime='23:59:59',
	31	walk=1,
	32	timezone='lt')
	33
	34	#AMISR Processing Unit
	35	procUnitAMISRBeam0 = controllerObj.addProcUnit(datatype='AMISRProc', inputId=readUnitConfObj.getId())
	36
	37	opObj11 = procUnitAMISRBeam0.addOperation(name='PrintInfo', optype='other')
	38
	39	#Reshaper
	40	opObj11 = procUnitAMISRBeam0.addOperation(name='ProfileToChannels', optype='other')
	41
	42
	43	#Beam Selector
	44	#opObj11 = procUnitAMISRBeam0.addOperation(name='BeamSelector', optype='other')
	45	#opObj11.addParameter(name='beam', value='0', format='int')
	46
	47	#Voltage Processing Unit
	48	procUnitConfObjBeam0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=procUnitAMISRBeam0.getId())
	49	#Coherent Integration
	50	opObj11 = procUnitConfObjBeam0.addOperation(name='CohInt', optype='other')
	51	opObj11.addParameter(name='n', value='8', format='int')
	52	#Spectra Unit Processing, getting spectras with nProfiles and nFFTPoints
	53	procUnitConfObjSpectraBeam0 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjBeam0.getId())
	54	procUnitConfObjSpectraBeam0.addParameter(name='nFFTPoints', value=32, format='int')
	55	procUnitConfObjSpectraBeam0.addParameter(name='nProfiles', value=32, format='int')
	56	#RemoveDc
	57	opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='removeDC')
	58
	59	#Noise Estimation
	60	opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='getNoise')
	61	opObj11.addParameter(name='minHei', value='5', format='float')
	62	opObj11.addParameter(name='maxHei', value='20', format='float')
	63
	64	#SpectraPlot
	65	opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='SpectraPlot', optype='other')
	66	opObj11.addParameter(name='id', value='100', format='int')
	67	opObj11.addParameter(name='wintitle', value='AMISR Beam 0', format='str')
	68	opObj11.addParameter(name='zmin', value='30', format='int')
	69	opObj11.addParameter(name='zmax', value='80', format='int')
	70
	71	#RTIPlot
	72	#title0 = 'RTI AMISR Beam 0'
	73	#opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='RTIPlot', optype='other')
	74	#opObj11.addParameter(name='id', value='200', format='int')
	75	#opObj11.addParameter(name='wintitle', value=title0, format='str')
	76	#opObj11.addParameter(name='showprofile', value='0', format='int')
	77	##Setting RTI time using xmin,xmax
	78	#opObj11.addParameter(name='xmin', value='15', format='int')
	79	#opObj11.addParameter(name='xmax', value='23', format='int')
	80	#Setting dB range with zmin, zmax
	81	#opObj11.addParameter(name='zmin', value='45', format='int')
	82	#opObj11.addParameter(name='zmax', value='70', format='int')
	83	#Save RTI
	84	#figfile0 = 'amisr_rti_beam0.png'
	85	#opObj11.addParameter(name='figpath', value=figpath, format='str')
	86	#opObj11.addParameter(name='figfile', value=figfile0, format='str')
	87
	88	#-----------------------------------------------------------------------------------------------
	89	procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObjSpectraBeam0 .getId())
	90	opObj20 = procUnitConfObj2.addOperation(name='GetMoments')
	91
	92	# opObj21 = procUnitConfObj2.addOperation(name='MomentsPlot', optype='other')
	93	# opObj21.addParameter(name='id', value='3', format='int')
	94	# opObj21.addParameter(name='wintitle', value='Moments Plot', format='str')
	95	# opObj21.addParameter(name='save', value='1', format='bool')
	96	# opObj21.addParameter(name='figpath', value=pathFigure, format='str')
	97	# opObj21.addParameter(name='zmin', value='5', format='int')
	98	# opObj21.addParameter(name='zmax', value='90', format='int')
	99
	100	opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other')
	101	opObj22.addParameter(name='technique', value='DBS', format='str')
	102	opObj22.addParameter(name='azimuth', value='51.06', format='float')
	103	opObj22.addParameter(name='correctFactor', value='-1', format='float')
	104	opObj22.addParameter(name='dirCosx', value='9.63770247e-01, 5.93066547e-17, 1, 5.93066547e-17,-9.68583161e-01', format='floatlist')
	105	opObj22.addParameter(name='dirCosy', value=' 0,-0.96858316,0,0.96858316,0', format='floatlist')
	106	# opObj22.addParameter(name='horizontalOnly', value='1', format='bool')
	107	# opObj22.addParameter(name='channelList', value='1,2', format='intlist')
	108
	109	opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
	110	opObj23.addParameter(name='id', value='4', format='int')
	111	opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str')
	112	opObj23.addParameter(name='save', value='1', format='bool')
	113	opObj23.addParameter(name='figpath', value = pathFigure, format='str')
	114	opObj23.addParameter(name='zmin', value='-10', format='int')
	115	opObj23.addParameter(name='zmax', value='10', format='int')
	116	opObj23.addParameter(name='zmin_ver', value='-80', format='float')
	117	opObj23.addParameter(name='zmax_ver', value='80', format='float')
	118	opObj23.addParameter(name='SNRmin', value='-10', format='int')
	119	opObj23.addParameter(name='SNRmax', value='60', format='int')
	120	opObj23.addParameter(name='SNRthresh', value='0', format='float')
	121	opObj23.addParameter(name='xmin', value=xmin, format='float')
	122	opObj23.addParameter(name='xmax', value=xmax, format='float')
	123
	124	print "Escribiendo el archivo XML"
	125	controllerObj.writeXml(filename)
	126	print "Leyendo el archivo XML"
	127	controllerObj.readXml(filename)
	128
	129	controllerObj.createObjects()
	130	controllerObj.connectObjects()
	131	controllerObj.run()

schainpy/model/data/jrodata.py +221 0

             '''
             $Author: murco $
             $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
             '''
             import copy
             import numpy
             import datetime
             from jroheaderIO import SystemHeader, RadarControllerHeader
             def getNumpyDtype(dataTypeCode):
                 if dataTypeCode == 0:
                     numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
                 elif dataTypeCode == 1:
                     numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
                 elif dataTypeCode == 2:
                     numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
                 elif dataTypeCode == 3:
                     numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
                 elif dataTypeCode == 4:
                     numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
                 elif dataTypeCode == 5:
                     numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
                 else:
                     raise ValueError, 'dataTypeCode was not defined'
                 return numpyDtype
             def getDataTypeCode(numpyDtype):
                 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
                     datatype = 0
                 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
                     datatype = 1
                 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
                     datatype = 2
                 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
                     datatype = 3
                 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
                     datatype = 4
                 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
                     datatype = 5
                 else:
                     datatype = None
                 return datatype
             def hildebrand_sekhon(data, navg):
                 data = data.copy()
                 sortdata = numpy.sort(data,axis=None)
                 lenOfData = len(sortdata)
                 nums_min = lenOfData/10
                 if (lenOfData/10) > 2:
                     nums_min = lenOfData/10
                 else:
                     nums_min = 2
                 sump = 0.
                 sumq = 0.
                 j = 0
                 cont = 1
                 while((cont==1)and(j<lenOfData)):
                     sump += sortdata[j]
                     sumq += sortdata[j]**2
                     j += 1
                     if j > nums_min:
                         rtest = float(j)/(j-1) + 1.0/navg
                         if ((sumq*j) > (rtest*sump**2)):
                             j = j - 1
                             sump  = sump - sortdata[j]
                             sumq =  sumq - sortdata[j]**2
                             cont = 0
                 lnoise = sump /j
                 stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
                 return lnoise
             class Beam:
                 def __init__(self):
                     self.codeList = []
                     self.azimuthList = []
                     self.zenithList = []
             class GenericData(object):
                 flagNoData = True
                 def __init__(self):
                     raise ValueError, "This class has not been implemented"
                 def copy(self, inputObj=None):
                     if inputObj == None:
                         return copy.deepcopy(self)
                     for key in inputObj.__dict__.keys():
                         self.__dict__[key] = inputObj.__dict__[key]
                 def deepcopy(self):
                     return copy.deepcopy(self)
                 def isEmpty(self):
                     return self.flagNoData
             class JROData(GenericData):
             #    m_BasicHeader = BasicHeader()
             #    m_ProcessingHeader = ProcessingHeader()
                 systemHeaderObj = SystemHeader()
                 radarControllerHeaderObj = RadarControllerHeader()
             #    data = None
                 type = None
                 datatype = None #dtype but in string
             #     dtype = None
             #    nChannels = None
             #    nHeights = None
                 nProfiles = None
                 heightList = None
                 channelList = None
                 flagTimeBlock = False
                 useLocalTime = False
                 utctime = None
                 timeZone = None
                 dstFlag = None
                 errorCount = None
                 blocksize = None
                 nCode = None
                 nBaud = None
                 code = None
                 flagDecodeData = False #asumo q la data no esta decodificada
                 flagDeflipData = False #asumo q la data no esta sin flip
                 flagShiftFFT = False
             #     ippSeconds = None
                 timeInterval = None
                 nCohInt = None
                 noise = None
                 windowOfFilter = 1
                 #Speed of ligth
                 C = 3e8
                 frequency = 49.92e6
                 realtime = False
                 beacon_heiIndexList = None
                 last_block = None
                 blocknow = None
                 azimuth = None
                 zenith = None
                 beam = Beam()
                 def __init__(self):
                     raise ValueError, "This class has not been implemented"
                 def getNoise(self):
                     raise ValueError, "Not implemented"
                 def getNChannels(self):
                     return len(self.channelList)
                 def getChannelIndexList(self):
                     return range(self.nChannels)
                 def getNHeights(self):
                     return len(self.heightList)
                 def getHeiRange(self, extrapoints=0):
                     heis = self.heightList
             #        deltah = self.heightList[1] - self.heightList[0]
             #
             #        heis.append(self.heightList[-1])
                     return heis
                 def getltctime(self):
                     if self.useLocalTime:
                         return self.utctime - self.timeZone*60
                     return self.utctime
                 def getDatatime(self):
                     datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
                     return datatimeValue
                 def getTimeRange(self):
                     datatime = []
                     datatime.append(self.ltctime)
                     datatime.append(self.ltctime + self.timeInterval)
                     datatime = numpy.array(datatime)
                     return datatime
                 def getFmax(self):
                     PRF = 1./(self.ippSeconds * self.nCohInt)
                     fmax = PRF/2.
                     return fmax
                 def getVmax(self):
                     _lambda = self.C/self.frequency
                     vmax = self.getFmax() * _lambda
                     return vmax
                 def get_ippSeconds(self):
                     '''
                     '''
                     return self.radarControllerHeaderObj.ippSeconds
                 def set_ippSeconds(self, ippSeconds):
                     '''
                     '''
                     self.radarControllerHeaderObj.ippSeconds = ippSeconds
                     return
                 def get_dtype(self):
                     '''
                     '''
                     return getNumpyDtype(self.datatype)
                 def set_dtype(self, numpyDtype):
                     '''
                     '''
                     self.datatype = getDataTypeCode(numpyDtype)
                 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
                 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
                 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
                 #noise = property(getNoise, "I'm the 'nHeights' property.")
                 datatime = property(getDatatime, "I'm the 'datatime' property")
                 ltctime = property(getltctime, "I'm the 'ltctime' property")
                 ippSeconds = property(get_ippSeconds, set_ippSeconds)
                 dtype = property(get_dtype, set_dtype)
             class Voltage(JROData):
                 #data es un numpy array de 2 dmensiones (canales, alturas)
                 data = None
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Voltage"
                     self.data = None
             #         self.dtype = None
             #        self.nChannels = 0
             #        self.nHeights = 0
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #        self.channelIndexList = None
                     self.flagNoData = True
                     self.flagTimeBlock = False
                     self.utctime = None
                     self.timeZone = None
                     self.dstFlag = None
                     self.errorCount = None
                     self.nCohInt = None
                     self.blocksize = None
                     self.flagDecodeData = False #asumo q la data no esta decodificada
                     self.flagDeflipData = False #asumo q la data no esta sin flip
                     self.flagShiftFFT = False
                 def getNoisebyHildebrand(self):
                     """
                     Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
                     Return:
                         noiselevel
                     """
                     for channel in range(self.nChannels):
                         daux = self.data_spc[channel,:,:]
                         self.noise[channel] = hildebrand_sekhon(daux, self.nCohInt)
                     return self.noise
                 def getNoise(self, type = 1):
                     self.noise = numpy.zeros(self.nChannels)
                     if type == 1:
                         noise = self.getNoisebyHildebrand()
                     return 10*numpy.log10(noise)
                 noise = property(getNoise, "I'm the 'nHeights' property.")
             class Spectra(JROData):
                 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
                 data_spc = None
                 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
                 data_cspc = None
                 #data es un numpy array de 2 dmensiones (canales, alturas)
                 data_dc = None
                 nFFTPoints = None
             #     nPairs = None
                 pairsList = None
                 nIncohInt = None
                 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
                 nCohInt = None #se requiere para determinar el valor de timeInterval
                 ippFactor = None
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Spectra"
             #        self.data = None
             #         self.dtype = None
             #        self.nChannels = 0
             #        self.nHeights = 0
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #        self.channelIndexList = None
                     self.pairsList = None
                     self.flagNoData = True
                     self.flagTimeBlock = False
                     self.utctime = None
                     self.nCohInt = None
                     self.nIncohInt = None
                     self.blocksize = None
                     self.nFFTPoints = None
                     self.wavelength = None
                     self.flagDecodeData = False #asumo q la data no esta decodificada
                     self.flagDeflipData = False #asumo q la data no esta sin flip
                     self.flagShiftFFT = False
                     self.ippFactor = 1
                     #self.noise = None
                     self.beacon_heiIndexList = []
                     self.noise_estimation = None
                 def getNoisebyHildebrand(self):
                     """
                     Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
                     Return:
                         noiselevel
                     """
                     noise = numpy.zeros(self.nChannels)
                     for channel in range(self.nChannels):
                         daux = self.data_spc[channel,:,:]
                         noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
                     return noise
                 def getNoise(self):
                     if self.noise_estimation != None:
                         return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
                     else:
                         noise = self.getNoisebyHildebrand()
                         return noise
                 def getFreqRange(self, extrapoints=0):
                     deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
                     freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
                     return freqrange
                 def getVelRange(self, extrapoints=0):
                     deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
                     velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
                     return velrange
                 def getNPairs(self):
                     return len(self.pairsList)
                 def getPairsIndexList(self):
                     return range(self.nPairs)
                 def getNormFactor(self):
                     pwcode = 1
                     if self.flagDecodeData:
                         pwcode = numpy.sum(self.code[0]**2)
                     #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
                     normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
                     return normFactor
                 def getFlagCspc(self):
                     if self.data_cspc == None:
                         return True
                     return False
                 def getFlagDc(self):
                     if self.data_dc == None:
                         return True
                     return False
                 nPairs = property(getNPairs, "I'm the 'nPairs' property.")
                 pairsIndexList = property(getPairsIndexList, "I'm the 'pairsIndexList' property.")
                 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
                 flag_cspc = property(getFlagCspc)
                 flag_dc = property(getFlagDc)
                 noise = property(getNoise, "I'm the 'nHeights' property.")
             class SpectraHeis(Spectra):
                 data_spc = None
                 data_cspc = None
                 data_dc = None
                 nFFTPoints = None
             #     nPairs = None
                 pairsList = None
                 nIncohInt = None
                 def __init__(self):
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "SpectraHeis"
             #         self.dtype = None
             #        self.nChannels = 0
             #        self.nHeights = 0
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #        self.channelIndexList = None
                     self.flagNoData = True
                     self.flagTimeBlock = False
             #         self.nPairs = 0
                     self.utctime = None
                     self.blocksize = None
                 def getNormFactor(self):
                     pwcode = 1
                     if self.flagDecodeData:
                         pwcode = numpy.sum(self.code[0]**2)
                     normFactor = self.nIncohInt*self.nCohInt*pwcode
                     return normFactor
                 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
             class Fits:
                 heightList = None
                 channelList = None
                 flagNoData = True
                 flagTimeBlock = False
                 useLocalTime = False
                 utctime = None
                 timeZone = None
             #     ippSeconds = None
                 timeInterval = None
                 nCohInt = None
                 nIncohInt = None
                 noise = None
                 windowOfFilter = 1
                 #Speed of ligth
                 C = 3e8
                 frequency = 49.92e6
                 realtime = False
                 def __init__(self):
                     self.type = "Fits"
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #         self.channelIndexList = None
                     self.flagNoData = True
                     self.utctime = None
                     self.nCohInt = None
                     self.nIncohInt = None
                     self.useLocalTime = True
             #         self.utctime = None
             #         self.timeZone = None
             #         self.ltctime = None
             #         self.timeInterval = None
             #         self.header = None
             #         self.data_header = None
             #         self.data = None
             #         self.datatime = None
             #         self.flagNoData = False
             #         self.expName = ''
             #         self.nChannels = None
             #         self.nSamples = None
             #         self.dataBlocksPerFile = None
             #         self.comments = ''
             #
                 def getltctime(self):
                     if self.useLocalTime:
                         return self.utctime - self.timeZone*60
                     return self.utctime
                 def getDatatime(self):
                     datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
                     return datatime
                 def getTimeRange(self):
                     datatime = []
                     datatime.append(self.ltctime)
                     datatime.append(self.ltctime + self.timeInterval)
                     datatime = numpy.array(datatime)
                     return datatime
                 def getHeiRange(self):
                     heis = self.heightList
                     return heis
                 def isEmpty(self):
                     return self.flagNoData
                 def getNHeights(self):
                     return len(self.heightList)
                 def getNChannels(self):
                     return len(self.channelList)
                 def getChannelIndexList(self):
                     return range(self.nChannels)
                 def getNoise(self, type = 1):
                     self.noise = numpy.zeros(self.nChannels)
                     if type == 1:
                         noise = self.getNoisebyHildebrand()
                     if type == 2:
                         noise = self.getNoisebySort()
                     if type == 3:
                         noise = self.getNoisebyWindow()
                     return noise
                 datatime = property(getDatatime, "I'm the 'datatime' property")
                 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
                 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
                 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
                 noise = property(getNoise, "I'm the 'nHeights' property.")
                 datatime = property(getDatatime, "I'm the 'datatime' property")
                 ltctime = property(getltctime, "I'm the 'ltctime' property")
+            class Correlation(JROData):
+                noise = None
+                SNR = None
+                pairsAutoCorr = None    #Pairs of Autocorrelation
+                #--------------------------------------------------
+                data_corr = None
+                data_volt = None
+                lagT = None # each element value is a profileIndex
+                lagR = None # each element value is in km
+                pairsList = None
+                calculateVelocity = None
+                nPoints = None
+                nAvg = None
+                bufferSize = None
+                def __init__(self):
+                    '''
+                    Constructor
+                    '''
+                    self.radarControllerHeaderObj = RadarControllerHeader()
+                    self.systemHeaderObj = SystemHeader()
+                    self.type = "Correlation"
+                    self.data = None
+                    self.dtype = None
+                    self.nProfiles = None
+                    self.heightList = None
+                    self.channelList = None
+                    self.flagNoData = True
+                    self.flagTimeBlock = False
+                    self.utctime = None
+                    self.timeZone = None
+                    self.dstFlag = None
+                    self.errorCount = None
+                    self.blocksize = None
+                    self.flagDecodeData = False #asumo q la data no esta decodificada
+                    self.flagDeflipData = False #asumo q la data no esta sin flip
+                    self.pairsList = None
+                    self.nPoints = None
+                def getLagTRange(self, extrapoints=0):
+                    lagTRange = self.lagT
+                    diff = lagTRange[1] - lagTRange[0]
+                    extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1]
+                    lagTRange = numpy.hstack((lagTRange, extra))
+                    return lagTRange
+                def getLagRRange(self, extrapoints=0):
+                    return self.lagR
+                def getPairsList(self):
+                    return self.pairsList
+                def getCalculateVelocity(self):
+                    return self.calculateVelocity
+                def getNPoints(self):
+                    return self.nPoints
+                def getNAvg(self):
+                    return self.nAvg
+                def getBufferSize(self):
+                    return self.bufferSize
+                def getPairsAutoCorr(self):
+                    pairsList = self.pairsList
+                    pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan
+                    for l in range(len(pairsList)):
+                        firstChannel = pairsList[l][0]
+                        secondChannel = pairsList[l][1]
+                        #Obteniendo pares de Autocorrelacion
+                        if firstChannel == secondChannel:
+                            pairsAutoCorr[firstChannel] = int(l)
+                    pairsAutoCorr = pairsAutoCorr.astype(int)
+                    return pairsAutoCorr
+                def getNoise(self, mode = 2):
+                    indR = numpy.where(self.lagR == 0)[0][0]
+                    indT = numpy.where(self.lagT == 0)[0][0]
+                    jspectra0 = self.data_corr[:,:,indR,:]
+                    jspectra = copy.copy(jspectra0)
+                    num_chan = jspectra.shape[0]
+                    num_hei = jspectra.shape[2]
+                    freq_dc = jspectra.shape[1]/2
+                    ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
+                    if ind_vel[0]<0:
+                        ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
+                    if mode == 1:
+                        jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
+                    if mode == 2:
+                        vel = numpy.array([-2,-1,1,2])
+                        xx = numpy.zeros([4,4])
+                        for fil in range(4):
+                            xx[fil,:] = vel[fil]**numpy.asarray(range(4))
+                        xx_inv = numpy.linalg.inv(xx)
+                        xx_aux = xx_inv[0,:]
+                        for ich in range(num_chan):
+                            yy = jspectra[ich,ind_vel,:]
+                            jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
+                            junkid = jspectra[ich,freq_dc,:]<=0
+                            cjunkid = sum(junkid)
+                            if cjunkid.any():
+                                jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
+                    noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
+                    return noise
+            #     pairsList = property(getPairsList, "I'm the 'pairsList' property.")
+            #     nPoints = property(getNPoints, "I'm the 'nPoints' property.")
+                calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.")
+                nAvg = property(getNAvg, "I'm the 'nAvg' property.")
+                bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.")
+            class Parameters(JROData):
+                inputUnit = None        #Type of data to be processed
+                operation = None        #Type of operation to parametrize
+                data_param = None       #Parameters obtained
+                data_pre = None         #Data Pre Parametrization
+                heightRange = None      #Heights
+                abscissaRange = None    #Abscissa, can be velocities, lags or time
+                noise = None            #Noise Potency
+                SNR = None              #Signal to Noise Ratio
+                pairsList = None        #List of Pairs for Cross correlations or Cross spectrum
+                initUtcTime = None      #Initial UTC time
+                paramInterval = None    #Time interval to calculate Parameters in seconds
+                windsInterval = None     #Time interval to calculate Winds in seconds
+                normFactor = None       #Normalization Factor
+                winds = None            #Wind estimations
+                def __init__(self):
+                    '''
+                    Constructor
+                    '''
+                    self.radarControllerHeaderObj = RadarControllerHeader()
+                    self.systemHeaderObj = SystemHeader()
+                    self.type = "Parameters"
+                def getTimeRange1(self):
+                    datatime = []
+                    datatime.append(self.initUtcTime)
+                    datatime.append(self.initUtcTime + self.windsInterval - 1)
+                    datatime = numpy.array(datatime)
+                    return datatime

schainpy/model/graphics/figure.py +17 0

             import os
             import numpy
             import time, datetime
             import mpldriver
             import Queue
             import threading
             def isRealtime(utcdatatime):
                 utcnow = time.mktime(time.localtime())
                 delta = abs(utcnow - utcdatatime) # abs
                 if delta >= 30.:
                     return False
                 return True
             class Figure:
                 __driver = mpldriver
                 __isConfigThread = False
                 fig = None
                 id = None
                 wintitle = None
                 width = None
                 height = None
                 nplots = None
                 timerange = None
                 axesObjList = []
                 WIDTH = None
                 HEIGHT = None
                 PREFIX = 'fig'
                 xmin = None
                 xmax = None
                 def __init__(self):
                     raise ValueError, "This method is not implemented"
                 def __del__(self):
                     self.__driver.closeFigure()
                 def getFilename(self, name, ext='.png'):
                     path = '%s%03d' %(self.PREFIX, self.id)
                     filename = '%s_%s%s' %(self.PREFIX, name, ext)
                     return os.path.join(path, filename)
                 def getAxesObjList(self):
                     return self.axesObjList
                 def getSubplots(self):
                     raise ValueError, "Abstract method: This method should be defined"
                 def getScreenDim(self, widthplot, heightplot):
                     nrow, ncol = self.getSubplots()
                     widthscreen = widthplot*ncol
                     heightscreen = heightplot*nrow
                     return widthscreen, heightscreen
                 def getTimeLim(self, x, xmin=None, xmax=None, timerange=None):
                     if self.xmin != None and self.xmax != None:
                         if timerange == None:
                             timerange = self.xmax - self.xmin
                         xmin = self.xmin + timerange
                         xmax = self.xmax + timerange
                         return xmin, xmax
                     if timerange != None and self.xmin == None and self.xmax == None:
                         txmin = x[0] - x[0]%timerange
                         thisdatetime = datetime.datetime.utcfromtimestamp(txmin)
                         thisdate = datetime.datetime.combine(thisdatetime.date(), datetime.time(0,0,0))
                         xmin = (thisdatetime - thisdate).seconds/(60*60.)
                         xmax = xmin + timerange/(60*60.)
                     if timerange == None:
                         txmin = numpy.min(x)
                         thisdatetime = datetime.datetime.utcfromtimestamp(txmin)
                         thisdate = datetime.datetime.combine(thisdatetime.date(), datetime.time(0,0,0))
                     mindt = thisdate + datetime.timedelta(hours=xmin) - datetime.timedelta(seconds=time.timezone)
                     xmin_sec = time.mktime(mindt.timetuple())
                     maxdt = thisdate + datetime.timedelta(hours=xmax) - datetime.timedelta(seconds=time.timezone)
                     xmax_sec = time.mktime(maxdt.timetuple())
                     return xmin_sec, xmax_sec
             #         if timerange != None:
             #             txmin = x[0] - x[0]%timerange
             #         else:
             #             txmin = numpy.min(x)
             #
             #         thisdatetime = datetime.datetime.utcfromtimestamp(txmin)
             #         thisdate = datetime.datetime.combine(thisdatetime.date(), datetime.time(0,0,0))
             #
             #         ####################################################
             #         #If the x is out of xrange
             #         if xmax != None:
             #             if xmax < (thisdatetime - thisdate).seconds/(60*60.):
             #                 xmin = None
             #                 xmax = None
             #
             #         if xmin == None:
             #             td = thisdatetime - thisdate
             #             xmin = td.seconds/(60*60.)
             #
             #         if xmax == None:
             #             xmax = xmin + self.timerange/(60*60.)
             #
             #         mindt = thisdate + datetime.timedelta(hours=xmin) - datetime.timedelta(seconds=time.timezone)
             #         tmin = time.mktime(mindt.timetuple())
             #
             #         maxdt = thisdate + datetime.timedelta(hours=xmax) - datetime.timedelta(seconds=time.timezone)
             #         tmax = time.mktime(maxdt.timetuple())
             #
             #         #self.timerange = tmax - tmin
             #
             #         return tmin, tmax
                 def init(self, id, nplots, wintitle):
                     raise ValueError, "This method has been replaced with createFigure"
                 def createFigure(self, id, wintitle, widthplot=None, heightplot=None, show=True):
                     """
                     Crea la figura de acuerdo al driver y parametros seleccionados seleccionados.
                     Las dimensiones de la pantalla es calculada a partir de los atributos self.WIDTH
                     y self.HEIGHT y el numero de subplots (nrow, ncol)
                     Input:
                         id    :    Los parametros necesarios son
                         wintitle    :
                     """
                     if widthplot == None:
                         widthplot = self.WIDTH
                     if heightplot == None:
                         heightplot = self.HEIGHT
                     self.id = id
                     self.wintitle = wintitle
                     self.widthscreen, self.heightscreen = self.getScreenDim(widthplot, heightplot)
                     self.fig = self.__driver.createFigure(id=self.id,
                                                           wintitle=self.wintitle,
                                                           width=self.widthscreen,
                                                           height=self.heightscreen,
                                                           show=show)
                     self.axesObjList = []
                 def setDriver(self, driver=mpldriver):
                     self.__driver = driver
                 def setTitle(self, title):
                     self.__driver.setTitle(self.fig, title)
                 def setWinTitle(self, title):
                     self.__driver.setWinTitle(self.fig, title=title)
                 def setTextFromAxes(self, text):
                     raise ValueError, "Este metodo ha sido reemplazaado con el metodo setText de la clase Axes"
                 def makeAxes(self, nrow, ncol, xpos, ypos, colspan, rowspan):
                     raise ValueError, "Este metodo ha sido reemplazaado con el metodo addAxes"
                 def addAxes(self, *args):
                     """
                     Input:
                         *args    :    Los parametros necesarios son
                                       nrow, ncol, xpos, ypos, colspan, rowspan
                     """
                     axesObj = Axes(self.fig, *args)
                     self.axesObjList.append(axesObj)
                 def saveFigure(self, figpath, figfile, *args):
                     filename = os.path.join(figpath, figfile)
                     fullpath = os.path.split(filename)[0]
                     if not os.path.exists(fullpath):
                         subpath = os.path.split(fullpath)[0]
                         if not os.path.exists(subpath):
                             os.mkdir(subpath)
                         os.mkdir(fullpath)
                     self.__driver.saveFigure(self.fig, filename, *args)
                 def getNameToFtp(self, thisDatetime, FTP_WEI, EXP_CODE, SUB_EXP_CODE, PLOT_CODE, PLOT_POS):
                     YEAR_STR = '%4.4d'%thisDatetime.timetuple().tm_year
                     DOY_STR = '%3.3d'%thisDatetime.timetuple().tm_yday
                     FTP_WEI = '%2.2d'%FTP_WEI
                     EXP_CODE = '%3.3d'%EXP_CODE
                     SUB_EXP_CODE = '%2.2d'%SUB_EXP_CODE
                     PLOT_CODE = '%2.2d'%PLOT_CODE
                     PLOT_POS = '%2.2d'%PLOT_POS
                     name = YEAR_STR + DOY_STR + FTP_WEI + EXP_CODE + SUB_EXP_CODE + PLOT_CODE + PLOT_POS
                     return name
                 def draw(self):
                     self.__driver.draw(self.fig)
                 def run(self):
                     raise ValueError, "This method is not implemented"
                 def close(self):
                     self.__driver.show(True)
                 axesList = property(getAxesObjList)
             class Axes:
                 __driver = mpldriver
                 fig = None
                 ax = None
                 plot = None
                 __missing = 1E30
                 __firsttime = None
                 __showprofile = False
                 xmin = None
                 xmax = None
                 ymin = None
                 ymax = None
                 zmin = None
                 zmax = None
                 x_buffer = None
                 z_buffer = None
                 decimationx = None
                 decimationy = None
                 __MAXNUMX = 300
                 __MAXNUMY = 150
                 def __init__(self, *args):
                     """
                     Input:
                         *args    :    Los parametros necesarios son
                                       fig, nrow, ncol, xpos, ypos, colspan, rowspan
                     """
                     ax = self.__driver.createAxes(*args)
                     self.fig = args[0]
                     self.ax = ax
                     self.plot = None
                     self.__firsttime = True
                     self.idlineList = []
                     self.x_buffer = numpy.array([])
                     self.z_buffer = numpy.array([])
                 def setText(self, text):
                     self.__driver.setAxesText(self.ax, text)
                 def setXAxisAsTime(self):
                     pass
                 def pline(self, x, y,
                            xmin=None, xmax=None,
                            ymin=None, ymax=None,
                            xlabel='', ylabel='',
                            title='',
                            **kwargs):
                     """
                     Input:
                         x        :
                         y        :
                         xmin    :
                         xmax    :
                         ymin    :
                         ymax    :
                         xlabel    :
                         ylabel    :
                         title    :
                         **kwargs :    Los parametros aceptados son
                                       ticksize
                                       ytick_visible
                     """
                     if self.__firsttime:
                         if xmin == None: xmin = numpy.nanmin(x)
                         if xmax == None: xmax = numpy.nanmax(x)
                         if ymin == None: ymin = numpy.nanmin(y)
                         if ymax == None: ymax = numpy.nanmax(y)
                         self.plot = self.__driver.createPline(self.ax, x, y,
                                                               xmin, xmax,
                                                               ymin, ymax,
                                                               xlabel=xlabel,
                                                                 ylabel=ylabel,
                                                                 title=title,
                                                               **kwargs)
                         self.idlineList.append(0)
                         self.__firsttime = False
                         return
                     self.__driver.pline(self.plot, x, y, xlabel=xlabel,
                                                                 ylabel=ylabel,
                                                                 title=title)
                 def addpline(self, x, y, idline, **kwargs):
                     lines = self.ax.lines
                     if idline in self.idlineList:
                         self.__driver.set_linedata(self.ax, x, y, idline)
                     if  idline not in(self.idlineList):
                         self.__driver.addpline(self.ax, x, y, **kwargs)
                         self.idlineList.append(idline)
                     return
                 def pmultiline(self, x, y,
                                xmin=None, xmax=None,
                                ymin=None, ymax=None,
                                xlabel='', ylabel='',
                                title='',
                                **kwargs):
                     if self.__firsttime:
                         if xmin == None: xmin = numpy.nanmin(x)
                         if xmax == None: xmax = numpy.nanmax(x)
                         if ymin == None: ymin = numpy.nanmin(y)
                         if ymax == None: ymax = numpy.nanmax(y)
                         self.plot = self.__driver.createPmultiline(self.ax, x, y,
                                                               xmin, xmax,
                                                               ymin, ymax,
                                                               xlabel=xlabel,
                                                               ylabel=ylabel,
                                                               title=title,
                                                               **kwargs)
                         self.__firsttime = False
                         return
                     self.__driver.pmultiline(self.plot, x, y, xlabel=xlabel,
                                                                 ylabel=ylabel,
                                                                 title=title)
                 def pmultilineyaxis(self, x, y,
                                xmin=None, xmax=None,
                                ymin=None, ymax=None,
                                xlabel='', ylabel='',
                                title='',
                                **kwargs):
                     if self.__firsttime:
                         if xmin == None: xmin = numpy.nanmin(x)
                         if xmax == None: xmax = numpy.nanmax(x)
                         if ymin == None: ymin = numpy.nanmin(y)
                         if ymax == None: ymax = numpy.nanmax(y)
                         self.plot = self.__driver.createPmultilineYAxis(self.ax, x, y,
                                                               xmin, xmax,
                                                               ymin, ymax,
                                                               xlabel=xlabel,
                                                               ylabel=ylabel,
                                                               title=title,
                                                               **kwargs)
                         if self.xmin == None: self.xmin = xmin
                         if self.xmax == None: self.xmax = xmax
                         if self.ymin == None: self.ymin = ymin
                         if self.ymax == None: self.ymax = ymax
                         self.__firsttime = False
                         return
                     self.__driver.pmultilineyaxis(self.plot, x, y, xlabel=xlabel,
                                                                 ylabel=ylabel,
                                                                 title=title)
                 def pcolor(self, x, y, z,
                            xmin=None, xmax=None,
                            ymin=None, ymax=None,
                            zmin=None, zmax=None,
                            xlabel='', ylabel='',
                            title='', rti = False, colormap='jet',
                            **kwargs):
                     """
                     Input:
                         x        :
                         y        :
                         x        :
                         xmin    :
                         xmax    :
                         ymin    :
                         ymax    :
                         zmin    :
                         zmax    :
                         xlabel    :
                         ylabel    :
                         title    :
                         **kwargs :    Los parametros aceptados son
                                       ticksize=9,
                                       cblabel=''
                                       rti = True or False
                     """
                     if self.__firsttime:
                         if xmin == None: xmin = numpy.nanmin(x)
                         if xmax == None: xmax = numpy.nanmax(x)
                         if ymin == None: ymin = numpy.nanmin(y)
                         if ymax == None: ymax = numpy.nanmax(y)
                         if zmin == None: zmin = numpy.nanmin(z)
                         if zmax == None: zmax = numpy.nanmax(z)
                         self.plot = self.__driver.createPcolor(self.ax, x, y, z,
                                                                xmin, xmax,
                                                                ymin, ymax,
                                                                zmin, zmax,
                                                                xlabel=xlabel,
                                                                 ylabel=ylabel,
                                                                 title=title,
                                                                 colormap=colormap,
                                                                 **kwargs)
                         if self.xmin == None: self.xmin = xmin
                         if self.xmax == None: self.xmax = xmax
                         if self.ymin == None: self.ymin = ymin
                         if self.ymax == None: self.ymax = ymax
                         if self.zmin == None: self.zmin = zmin
                         if self.zmax == None: self.zmax = zmax
                         self.__firsttime = False
                         return
                     if rti:
                         self.__driver.addpcolor(self.ax, x, y, z, self.zmin, self.zmax,
                                                 xlabel=xlabel,
                                                 ylabel=ylabel,
                                                 title=title,
                                                 colormap=colormap)
                         return
                     self.__driver.pcolor(self.plot, z,
                                          xlabel=xlabel,
                                          ylabel=ylabel,
                                          title=title)
                 def pcolorbuffer(self, x, y, z,
                            xmin=None, xmax=None,
                            ymin=None, ymax=None,
                            zmin=None, zmax=None,
                            xlabel='', ylabel='',
                            title='', rti = True, colormap='jet',
                            maxNumX = None, maxNumY = None,
                            **kwargs):
                     if maxNumX == None:
                         maxNumX = self.__MAXNUMX
                     if maxNumY == None:
                         maxNumY = self.__MAXNUMY
                     if self.__firsttime:
                         self.z_buffer = z
                         self.x_buffer = numpy.hstack((self.x_buffer, x))
                         if xmin == None: xmin = numpy.nanmin(x)
                         if xmax == None: xmax = numpy.nanmax(x)
                         if ymin == None: ymin = numpy.nanmin(y)
                         if ymax == None: ymax = numpy.nanmax(y)
                         if zmin == None: zmin = numpy.nanmin(z)
                         if zmax == None: zmax = numpy.nanmax(z)
                         self.plot = self.__driver.createPcolor(self.ax, self.x_buffer, y, z,
                                                                xmin, xmax,
                                                                ymin, ymax,
                                                                zmin, zmax,
                                                                xlabel=xlabel,
                                                                 ylabel=ylabel,
                                                                 title=title,
                                                                 colormap=colormap,
                                                                 **kwargs)
                         if self.xmin == None: self.xmin = xmin
                         if self.xmax == None: self.xmax = xmax
                         if self.ymin == None: self.ymin = ymin
                         if self.ymax == None: self.ymax = ymax
                         if self.zmin == None: self.zmin = zmin
                         if self.zmax == None: self.zmax = zmax
                         self.__firsttime = False
                         return
                     self.x_buffer = numpy.hstack((self.x_buffer, x[-1]))
                     self.z_buffer = numpy.hstack((self.z_buffer, z))
                     if self.decimationx == None:
                         deltax = float(self.xmax - self.xmin)/maxNumX
                         deltay = float(self.ymax - self.ymin)/maxNumY
                         resolutionx = self.x_buffer[2]-self.x_buffer[0]
                         resolutiony = y[1]-y[0]
                         self.decimationx = numpy.ceil(deltax / resolutionx)
                         self.decimationy = numpy.ceil(deltay / resolutiony)
                     z_buffer = self.z_buffer.reshape(-1,len(y))
                     x_buffer = self.x_buffer[::self.decimationx]
                     y_buffer = y[::self.decimationy]
                     z_buffer = z_buffer[::self.decimationx, ::self.decimationy]
                     #===================================================
                     x_buffer, y_buffer, z_buffer = self.__fillGaps(x_buffer, y_buffer, z_buffer)
                     self.__driver.addpcolorbuffer(self.ax, x_buffer, y_buffer, z_buffer, self.zmin, self.zmax,
                                             xlabel=xlabel,
                                             ylabel=ylabel,
                                             title=title,
                                             colormap=colormap)
+                def polar(self, x, y,
+                           title='', xlabel='',ylabel='',**kwargs):
+                    if self.__firsttime:
+                        self.plot = self.__driver.createPolar(self.ax, x, y, title = title, xlabel = xlabel, ylabel = ylabel)
+                        self.__firsttime = False
+                        self.x_buffer = x
+                        self.y_buffer = y
+                        return
+                    self.x_buffer = numpy.hstack((self.x_buffer,x))
+                    self.y_buffer = numpy.hstack((self.y_buffer,y))
+                    self.__driver.polar(self.plot, self.x_buffer, self.y_buffer, xlabel=xlabel,
+                                                        ylabel=ylabel,
+                                                        title=title)
                 def __fillGaps(self, x_buffer, y_buffer, z_buffer):
                     deltas = x_buffer[1:] - x_buffer[0:-1]
                     x_median = numpy.median(deltas)
                     index = numpy.where(deltas >= 2*x_median)
                     if len(index[0]) != 0:
                         z_buffer[index[0],::] = self.__missing
                         z_buffer = numpy.ma.masked_inside(z_buffer,0.99*self.__missing,1.01*self.__missing)
                     return x_buffer, y_buffer, z_buffer
   No newline at end of file

schainpy/model/graphics/jroplot.py +2 0

             from jroplot_voltage import *
             from jroplot_spectra import *
             from jroplot_heispectra import *
+            from jroplot_correlation import *
+            from jroplot_parameters import *
  No newline at end of file

schainpy/model/graphics/jroplot_spectra.py +1 -1

             '''
             @author: Daniel Suarez
             '''
             import os
             import datetime
             import numpy
             from figure import Figure, isRealtime
             class SpectraPlot(Figure):
                 isConfig = None
                 __nsubplots = None
                 WIDTHPROF = None
                 HEIGHTPROF = None
                 PREFIX = 'spc'
                 def __init__(self):
                     self.isConfig = False
                     self.__nsubplots = 1
                     self.WIDTH = 280
                     self.HEIGHT = 250
                     self.WIDTHPROF = 120
                     self.HEIGHTPROF = 0
                     self.counter_imagwr = 0
                     self.PLOT_CODE = 1
                     self.FTP_WEI = None
                     self.EXP_CODE = None
                     self.SUB_EXP_CODE = None
                     self.PLOT_POS = None
                 def getSubplots(self):
                     ncol = int(numpy.sqrt(self.nplots)+0.9)
                     nrow = int(self.nplots*1./ncol + 0.9)
                     return nrow, ncol
                 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
                     self.__showprofile = showprofile
                     self.nplots = nplots
                     ncolspan = 1
                     colspan = 1
                     if showprofile:
                         ncolspan = 3
                         colspan = 2
                         self.__nsubplots = 2
                     self.createFigure(id = id,
                                       wintitle = wintitle,
                                       widthplot = self.WIDTH + self.WIDTHPROF,
                                       heightplot = self.HEIGHT + self.HEIGHTPROF,
                                       show=show)
                     nrow, ncol = self.getSubplots()
                     counter = 0
                     for y in range(nrow):
                         for x in range(ncol):
                             if counter >= self.nplots:
                                 break
                             self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
                             if showprofile:
                                 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
                             counter += 1
                 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
                         xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
                         save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
                         server=None, folder=None, username=None, password=None,
                         ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
                     """
                     Input:
                         dataOut         :
                         id        :
                         wintitle        :
                         channelList     :
                         showProfile     :
                         xmin            :    None,
                         xmax            :    None,
                         ymin            :    None,
                         ymax            :    None,
                         zmin            :    None,
                         zmax            :    None
                     """
                     if dataOut.flagNoData:
                         return None
                     if realtime:
                         if not(isRealtime(utcdatatime = dataOut.utctime)):
                             print 'Skipping this plot function'
                             return
                     if channelList == None:
                         channelIndexList = dataOut.channelIndexList
                     else:
                         channelIndexList = []
                         for channel in channelList:
                             if channel not in dataOut.channelList:
                                 raise ValueError, "Channel %d is not in dataOut.channelList"
                             channelIndexList.append(dataOut.channelList.index(channel))
                     factor = dataOut.normFactor
                     x = dataOut.getVelRange(1)
                     y = dataOut.getHeiRange()
                     z = dataOut.data_spc[channelIndexList,:,:]/factor
                     z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
                     avg = numpy.average(z, axis=1)
-                    avg = numpy.nanmean(z, axis=1)
+                    #avg = numpy.nanmean(z, axis=1)
                     noise = dataOut.noise/factor
                     zdB = 10*numpy.log10(z)
                     avgdB = 10*numpy.log10(avg)
                     noisedB = 10*numpy.log10(noise)
                     #thisDatetime = dataOut.datatime
                     thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
                     title = wintitle + " Spectra"
                     if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
                         title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
                     xlabel = "Velocity (m/s)"
                     ylabel = "Range (Km)"
                     if not self.isConfig:
                         nplots = len(channelIndexList)
                         self.setup(id=id,
                                    nplots=nplots,
                                    wintitle=wintitle,
                                    showprofile=showprofile,
                                    show=show)
                         if xmin == None: xmin = numpy.nanmin(x)
                         if xmax == None: xmax = numpy.nanmax(x)
                         if ymin == None: ymin = numpy.nanmin(y)
                         if ymax == None: ymax = numpy.nanmax(y)
                         if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
                         if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
                         self.FTP_WEI = ftp_wei
                         self.EXP_CODE = exp_code
                         self.SUB_EXP_CODE = sub_exp_code
                         self.PLOT_POS = plot_pos
                         self.isConfig = True
                     self.setWinTitle(title)
                     for i in range(self.nplots):
                         str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
                         title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i]+1, noisedB[i], str_datetime)
                         if len(dataOut.beam.codeList) != 0:
                             title = "Ch%d:%4.2fdB,%2.2f,%2.2f:%s" %(dataOut.channelList[i]+1, noisedB[i], dataOut.beam.azimuthList[i], dataOut.beam.zenithList[i], str_datetime)
                         axes = self.axesList[i*self.__nsubplots]
                         axes.pcolor(x, y, zdB[i,:,:],
                                     xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
                                     xlabel=xlabel, ylabel=ylabel, title=title,
                                     ticksize=9, cblabel='')
                         if self.__showprofile:
                             axes = self.axesList[i*self.__nsubplots +1]
                             axes.pline(avgdB[i], y,
                                     xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
                                     xlabel='dB', ylabel='', title='',
                                     ytick_visible=False,
                                     grid='x')
                             noiseline = numpy.repeat(noisedB[i], len(y))
                             axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
                     self.draw()
                     if figfile == None:
                         str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         figfile = self.getFilename(name = str_datetime)
                         name = str_datetime
                         if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
                             name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith)
                         figfile = self.getFilename(name)
                     if figpath != '':
                         self.counter_imagwr += 1
                         if (self.counter_imagwr>=wr_period):
                             # store png plot to local folder
                             self.saveFigure(figpath, figfile)
                             # store png plot to FTP server according to RT-Web format
                             name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
                             ftp_filename = os.path.join(figpath, name)
                             self.saveFigure(figpath, ftp_filename)
                             self.counter_imagwr = 0
             class CrossSpectraPlot(Figure):
                 isConfig = None
                 __nsubplots = None
                 WIDTH = None
                 HEIGHT = None
                 WIDTHPROF = None
                 HEIGHTPROF = None
                 PREFIX = 'cspc'
                 def __init__(self):
                     self.isConfig = False
                     self.__nsubplots = 4
                     self.counter_imagwr = 0
                     self.WIDTH = 250
                     self.HEIGHT = 250
                     self.WIDTHPROF = 0
                     self.HEIGHTPROF = 0
                     self.PLOT_CODE = 1
                     self.FTP_WEI = None
                     self.EXP_CODE = None
                     self.SUB_EXP_CODE = None
                     self.PLOT_POS = None
                 def getSubplots(self):
                     ncol = 4
                     nrow = self.nplots
                     return nrow, ncol
                 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
                     self.__showprofile = showprofile
                     self.nplots = nplots
                     ncolspan = 1
                     colspan = 1
                     self.createFigure(id = id,
                                       wintitle = wintitle,
                                       widthplot = self.WIDTH + self.WIDTHPROF,
                                       heightplot = self.HEIGHT + self.HEIGHTPROF,
                                       show=True)
                     nrow, ncol = self.getSubplots()
                     counter = 0
                     for y in range(nrow):
                         for x in range(ncol):
                             self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
                             counter += 1
                 def run(self, dataOut, id, wintitle="", pairsList=None,
                         xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
                         save=False, figpath='', figfile=None, ftp=False, wr_period=1,
                         power_cmap='jet', coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
                         server=None, folder=None, username=None, password=None,
                         ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
                     """
                     Input:
                         dataOut         :
                         id        :
                         wintitle        :
                         channelList     :
                         showProfile     :
                         xmin            :    None,
                         xmax            :    None,
                         ymin            :    None,
                         ymax            :    None,
                         zmin            :    None,
                         zmax            :    None
                     """
                     if pairsList == None:
                         pairsIndexList = dataOut.pairsIndexList
                     else:
                         pairsIndexList = []
                         for pair in pairsList:
                             if pair not in dataOut.pairsList:
                                 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
                             pairsIndexList.append(dataOut.pairsList.index(pair))
                     if pairsIndexList == []:
                         return
                     if len(pairsIndexList) > 4:
                         pairsIndexList = pairsIndexList[0:4]
                     factor = dataOut.normFactor
                     x = dataOut.getVelRange(1)
                     y = dataOut.getHeiRange()
                     z = dataOut.data_spc[:,:,:]/factor
             #        z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
                     avg = numpy.abs(numpy.average(z, axis=1))
                     noise = dataOut.noise()/factor
                     zdB = 10*numpy.log10(z)
                     avgdB = 10*numpy.log10(avg)
                     noisedB = 10*numpy.log10(noise)
                     #thisDatetime = dataOut.datatime
                     thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
                     title = wintitle + " Cross-Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
                     xlabel = "Velocity (m/s)"
                     ylabel = "Range (Km)"
                     if not self.isConfig:
                         nplots = len(pairsIndexList)
                         self.setup(id=id,
                                    nplots=nplots,
                                    wintitle=wintitle,
                                    showprofile=False,
                                    show=show)
                         if xmin == None: xmin = numpy.nanmin(x)
                         if xmax == None: xmax = numpy.nanmax(x)
                         if ymin == None: ymin = numpy.nanmin(y)
                         if ymax == None: ymax = numpy.nanmax(y)
                         if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
                         if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
                         self.FTP_WEI = ftp_wei
                         self.EXP_CODE = exp_code
                         self.SUB_EXP_CODE = sub_exp_code
                         self.PLOT_POS = plot_pos
                         self.isConfig = True
                     self.setWinTitle(title)
                     for i in range(self.nplots):
                         pair = dataOut.pairsList[pairsIndexList[i]]
                         str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
                         title = "Ch%d: %4.2fdB: %s" %(pair[0], noisedB[pair[0]], str_datetime)
                         zdB = 10.*numpy.log10(dataOut.data_spc[pair[0],:,:]/factor)
                         axes0 = self.axesList[i*self.__nsubplots]
                         axes0.pcolor(x, y, zdB,
                                     xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
                                     xlabel=xlabel, ylabel=ylabel, title=title,
                                     ticksize=9, colormap=power_cmap, cblabel='')
                         title = "Ch%d: %4.2fdB: %s" %(pair[1], noisedB[pair[1]], str_datetime)
                         zdB = 10.*numpy.log10(dataOut.data_spc[pair[1],:,:]/factor)
                         axes0 = self.axesList[i*self.__nsubplots+1]
                         axes0.pcolor(x, y, zdB,
                                     xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
                                     xlabel=xlabel, ylabel=ylabel, title=title,
                                     ticksize=9, colormap=power_cmap, cblabel='')
                         coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[pair[0],:,:]*dataOut.data_spc[pair[1],:,:])
                         coherence = numpy.abs(coherenceComplex)
             #            phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
                         phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
                         title = "Coherence %d%d" %(pair[0], pair[1])
                         axes0 = self.axesList[i*self.__nsubplots+2]
                         axes0.pcolor(x, y, coherence,
                                     xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=0, zmax=1,
                                     xlabel=xlabel, ylabel=ylabel, title=title,
                                     ticksize=9, colormap=coherence_cmap, cblabel='')
                         title = "Phase %d%d" %(pair[0], pair[1])
                         axes0 = self.axesList[i*self.__nsubplots+3]
                         axes0.pcolor(x, y, phase,
                                     xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=-180, zmax=180,
                                     xlabel=xlabel, ylabel=ylabel, title=title,
                                     ticksize=9, colormap=phase_cmap, cblabel='')
                     self.draw()
                     if figfile == None:
                         str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         figfile = self.getFilename(name = str_datetime)
                     if figpath != '':
                         self.counter_imagwr += 1
                         if (self.counter_imagwr>=wr_period):
                             # store png plot to local folder
                             self.saveFigure(figpath, figfile)
                             # store png plot to FTP server according to RT-Web format
                             name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
                             ftp_filename = os.path.join(figpath, name)
                             self.saveFigure(figpath, ftp_filename)
                             self.counter_imagwr = 0
             class RTIPlot(Figure):
                 isConfig = None
                 __nsubplots = None
                 WIDTHPROF = None
                 HEIGHTPROF = None
                 PREFIX = 'rti'
                 def __init__(self):
                     self.timerange = 2*60*60
                     self.isConfig = False
                     self.__nsubplots = 1
                     self.WIDTH = 800
                     self.HEIGHT = 150
                     self.WIDTHPROF = 120
                     self.HEIGHTPROF = 0
                     self.counter_imagwr = 0
                     self.PLOT_CODE = 0
                     self.FTP_WEI = None
                     self.EXP_CODE = None
                     self.SUB_EXP_CODE = None
                     self.PLOT_POS = None
                     self.tmin = None
                     self.tmax = None
                     self.xmin = None
                     self.xmax = None
                     self.figfile = None
                 def getSubplots(self):
                     ncol = 1
                     nrow = self.nplots
                     return nrow, ncol
                 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
                     self.__showprofile = showprofile
                     self.nplots = nplots
                     ncolspan = 1
                     colspan = 1
                     if showprofile:
                         ncolspan = 7
                         colspan = 6
                         self.__nsubplots = 2
                     self.createFigure(id = id,
                                       wintitle = wintitle,
                                       widthplot = self.WIDTH + self.WIDTHPROF,
                                       heightplot = self.HEIGHT + self.HEIGHTPROF,
                                       show=show)
                     nrow, ncol = self.getSubplots()
                     counter = 0
                     for y in range(nrow):
                         for x in range(ncol):
                             if counter >= self.nplots:
                                 break
                             self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
                             if showprofile:
                                 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
                             counter += 1
                 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
                         xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
                         timerange=None,
                         save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
                         server=None, folder=None, username=None, password=None,
                         ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
                     """
                     Input:
                         dataOut         :
                         id        :
                         wintitle        :
                         channelList     :
                         showProfile     :
                         xmin            :    None,
                         xmax            :    None,
                         ymin            :    None,
                         ymax            :    None,
                         zmin            :    None,
                         zmax            :    None
                     """
                     if channelList == None:
                         channelIndexList = dataOut.channelIndexList
                     else:
                         channelIndexList = []
                         for channel in channelList:
                             if channel not in dataOut.channelList:
                                 raise ValueError, "Channel %d is not in dataOut.channelList"
                             channelIndexList.append(dataOut.channelList.index(channel))
                     if timerange != None:
                         self.timerange = timerange
                     #tmin = None
                     #tmax = None
                     factor = dataOut.normFactor
                     x = dataOut.getTimeRange()
                     y = dataOut.getHeiRange()
                     z = dataOut.data_spc[channelIndexList,:,:]/factor
                     z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
                     avg = numpy.average(z, axis=1)
                     avgdB = 10.*numpy.log10(avg)
             #        thisDatetime = dataOut.datatime
                     thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
                     title = wintitle + " RTI" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
                     xlabel = ""
                     ylabel = "Range (Km)"
                     if not self.isConfig:
                         nplots = len(channelIndexList)
                         self.setup(id=id,
                                    nplots=nplots,
                                    wintitle=wintitle,
                                    showprofile=showprofile,
                                    show=show)
                         self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
             #             if timerange != None:
             #                 self.timerange = timerange
             #                 self.xmin, self.tmax = self.getTimeLim(x, xmin, xmax, timerange)
                         if ymin == None: ymin = numpy.nanmin(y)
                         if ymax == None: ymax = numpy.nanmax(y)
                         if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
                         if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
                         self.FTP_WEI = ftp_wei
                         self.EXP_CODE = exp_code
                         self.SUB_EXP_CODE = sub_exp_code
                         self.PLOT_POS = plot_pos
                         self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         self.isConfig = True
                         self.figfile = figfile
                     self.setWinTitle(title)
                     if ((self.xmax - x[1]) < (x[1]-x[0])):
                         x[1] = self.xmax
                     for i in range(self.nplots):
                         title = "Channel %d: %s" %(dataOut.channelList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
                         if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
                             title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
                         axes = self.axesList[i*self.__nsubplots]
                         zdB = avgdB[i].reshape((1,-1))
                         axes.pcolorbuffer(x, y, zdB,
                                     xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
                                     xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
                                     ticksize=9, cblabel='', cbsize="1%")
                         if self.__showprofile:
                             axes = self.axesList[i*self.__nsubplots +1]
                             axes.pline(avgdB[i], y,
                                     xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
                                     xlabel='dB', ylabel='', title='',
                                     ytick_visible=False,
                                     grid='x')
                     self.draw()
                     if x[1] >= self.axesList[0].xmax:
                         self.counter_imagwr = wr_period
                         self.__isConfig = False
                     if self.figfile == None:
                         str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         self.figfile = self.getFilename(name = str_datetime)
                     if figpath != '':
                         self.counter_imagwr += 1
                         if (self.counter_imagwr>=wr_period):
                             # store png plot to local folder
                             self.saveFigure(figpath, self.figfile)
                             # store png plot to FTP server according to RT-Web format
                             name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
                             ftp_filename = os.path.join(figpath, name)
                             self.saveFigure(figpath, ftp_filename)
                             self.counter_imagwr = 0
             class CoherenceMap(Figure):
                 isConfig = None
                 __nsubplots = None
                 WIDTHPROF = None
                 HEIGHTPROF = None
                 PREFIX = 'cmap'
                 def __init__(self):
                     self.timerange = 2*60*60
                     self.isConfig = False
                     self.__nsubplots = 1
                     self.WIDTH = 800
                     self.HEIGHT = 150
                     self.WIDTHPROF = 120
                     self.HEIGHTPROF = 0
                     self.counter_imagwr = 0
                     self.PLOT_CODE = 3
                     self.FTP_WEI = None
                     self.EXP_CODE = None
                     self.SUB_EXP_CODE = None
                     self.PLOT_POS = None
                     self.counter_imagwr = 0
                     self.xmin = None
                     self.xmax = None
                 def getSubplots(self):
                     ncol = 1
                     nrow = self.nplots*2
                     return nrow, ncol
                 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
                     self.__showprofile = showprofile
                     self.nplots = nplots
                     ncolspan = 1
                     colspan = 1
                     if showprofile:
                         ncolspan = 7
                         colspan = 6
                         self.__nsubplots = 2
                     self.createFigure(id = id,
                                       wintitle = wintitle,
                                       widthplot = self.WIDTH + self.WIDTHPROF,
                                       heightplot = self.HEIGHT + self.HEIGHTPROF,
                                       show=True)
                     nrow, ncol = self.getSubplots()
                     for y in range(nrow):
                         for x in range(ncol):
                             self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
                             if showprofile:
                                 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
                 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
                         xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
                         timerange=None,
                         save=False, figpath='', figfile=None, ftp=False, wr_period=1,
                         coherence_cmap='jet', phase_cmap='RdBu_r', show=True,
                         server=None, folder=None, username=None, password=None,
                         ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
                     if pairsList == None:
                         pairsIndexList = dataOut.pairsIndexList
                     else:
                         pairsIndexList = []
                         for pair in pairsList:
                             if pair not in dataOut.pairsList:
                                 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
                             pairsIndexList.append(dataOut.pairsList.index(pair))
                     if timerange != None:
                         self.timerange = timerange
                     if pairsIndexList == []:
                         return
                     if len(pairsIndexList) > 4:
                         pairsIndexList = pairsIndexList[0:4]
             #         tmin = None
             #         tmax = None
                     x = dataOut.getTimeRange()
                     y = dataOut.getHeiRange()
                     #thisDatetime = dataOut.datatime
                     thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
                     title = wintitle + " CoherenceMap" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
                     xlabel = ""
                     ylabel = "Range (Km)"
                     if not self.isConfig:
                         nplots = len(pairsIndexList)
                         self.setup(id=id,
                                    nplots=nplots,
                                    wintitle=wintitle,
                                    showprofile=showprofile,
                                    show=show)
                         #tmin, tmax = self.getTimeLim(x, xmin, xmax)
                         self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
                         if ymin == None: ymin = numpy.nanmin(y)
                         if ymax == None: ymax = numpy.nanmax(y)
                         if zmin == None: zmin = 0.
                         if zmax == None: zmax = 1.
                         self.FTP_WEI = ftp_wei
                         self.EXP_CODE = exp_code
                         self.SUB_EXP_CODE = sub_exp_code
                         self.PLOT_POS = plot_pos
                         self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         self.isConfig = True
                     self.setWinTitle(title)
                     if ((self.xmax - x[1]) < (x[1]-x[0])):
                         x[1] = self.xmax
                     for i in range(self.nplots):
                         pair = dataOut.pairsList[pairsIndexList[i]]
                         ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i],:,:],axis=0)
                         powa = numpy.average(dataOut.data_spc[pair[0],:,:],axis=0)
                         powb = numpy.average(dataOut.data_spc[pair[1],:,:],axis=0)
                         avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
                         coherence = numpy.abs(avgcoherenceComplex)
                         z = coherence.reshape((1,-1))
                         counter = 0
                         title = "Coherence %d%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
                         axes = self.axesList[i*self.__nsubplots*2]
                         axes.pcolorbuffer(x, y, z,
                                     xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
                                     xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
                                     ticksize=9, cblabel='', colormap=coherence_cmap, cbsize="1%")
                         if self.__showprofile:
                             counter += 1
                             axes = self.axesList[i*self.__nsubplots*2 + counter]
                             axes.pline(coherence, y,
                                     xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
                                     xlabel='', ylabel='', title='', ticksize=7,
                                     ytick_visible=False, nxticks=5,
                                     grid='x')
                         counter += 1
                         phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
                         z = phase.reshape((1,-1))
                         title = "Phase %d%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
                         axes = self.axesList[i*self.__nsubplots*2 + counter]
                         axes.pcolorbuffer(x, y, z,
                                     xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=-180, zmax=180,
                                     xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
                                     ticksize=9, cblabel='', colormap=phase_cmap, cbsize="1%")
                         if self.__showprofile:
                             counter += 1
                             axes = self.axesList[i*self.__nsubplots*2 + counter]
                             axes.pline(phase, y,
                                     xmin=-180, xmax=180, ymin=ymin, ymax=ymax,
                                     xlabel='', ylabel='', title='', ticksize=7,
                                     ytick_visible=False, nxticks=4,
                                     grid='x')
                     self.draw()
                     if x[1] >= self.axesList[0].xmax:
                         self.counter_imagwr = wr_period
                         self.__isConfig = False
                     if figfile == None:
                         str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         figfile = self.getFilename(name = str_datetime)
                     if figpath != '':
                         self.counter_imagwr += 1
                         if (self.counter_imagwr>=wr_period):
                             # store png plot to local folder
                             self.saveFigure(figpath, figfile)
                             # store png plot to FTP server according to RT-Web format
                             name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
                             ftp_filename = os.path.join(figpath, name)
                             self.saveFigure(figpath, ftp_filename)
                             self.counter_imagwr = 0
             class PowerProfile(Figure):
                 isConfig = None
                 __nsubplots = None
                 WIDTHPROF = None
                 HEIGHTPROF = None
                 PREFIX = 'spcprofile'
                 def __init__(self):
                     self.isConfig = False
                     self.__nsubplots = 1
                     self.WIDTH = 300
                     self.HEIGHT = 500
                     self.counter_imagwr = 0
                 def getSubplots(self):
                     ncol = 1
                     nrow = 1
                     return nrow, ncol
                 def setup(self, id, nplots, wintitle, show):
                     self.nplots = nplots
                     ncolspan = 1
                     colspan = 1
                     self.createFigure(id = id,
                                       wintitle = wintitle,
                                       widthplot = self.WIDTH,
                                       heightplot = self.HEIGHT,
                                       show=show)
                     nrow, ncol = self.getSubplots()
                     counter = 0
                     for y in range(nrow):
                         for x in range(ncol):
                             self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
                 def run(self, dataOut, id, wintitle="", channelList=None,
                         xmin=None, xmax=None, ymin=None, ymax=None,
                         save=False, figpath='', figfile=None, show=True, wr_period=1,
                         server=None, folder=None, username=None, password=None,):
                     if dataOut.flagNoData:
                         return None
                     if channelList == None:
                         channelIndexList = dataOut.channelIndexList
                         channelList = dataOut.channelList
                     else:
                         channelIndexList = []
                         for channel in channelList:
                             if channel not in dataOut.channelList:
                                 raise ValueError, "Channel %d is not in dataOut.channelList"
                             channelIndexList.append(dataOut.channelList.index(channel))
                     try:
                         factor = dataOut.normFactor
                     except:
                         factor = 1
                     y = dataOut.getHeiRange()
                     #for voltage
                     if dataOut.type == 'Voltage':
                         x = dataOut.data[channelIndexList,:] * numpy.conjugate(dataOut.data[channelIndexList,:])
                         x = x.real
                         x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
                     #for spectra
                     if dataOut.type == 'Spectra':
                         x = dataOut.data_spc[channelIndexList,:,:]/factor
                         x = numpy.where(numpy.isfinite(x), x, numpy.NAN)
                         x = numpy.average(x, axis=1)
                     xdB = 10*numpy.log10(x)
                     thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
                     title = wintitle + " Power Profile %s" %(thisDatetime.strftime("%d-%b-%Y"))
                     xlabel = "dB"
                     ylabel = "Range (Km)"
                     if not self.isConfig:
                         nplots = 1
                         self.setup(id=id,
                                    nplots=nplots,
                                    wintitle=wintitle,
                                    show=show)
                         if ymin == None: ymin = numpy.nanmin(y)
                         if ymax == None: ymax = numpy.nanmax(y)
                         if xmin == None: xmin = numpy.nanmin(xdB)*0.9
                         if xmax == None: xmax = numpy.nanmax(xdB)*0.9
                         self.__isConfig = True
                     self.setWinTitle(title)
                     title = "Power Profile: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
                     axes = self.axesList[0]
                     legendlabels = ["channel %d"%x for x in channelList]
                     axes.pmultiline(xdB, y,
                             xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
                             xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels,
                             ytick_visible=True, nxticks=5,
                             grid='x')
                     self.draw()
                     if figfile == None:
                         str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         figfile = self.getFilename(name = str_datetime)
                     if figpath != '':
                         self.counter_imagwr += 1
                         if (self.counter_imagwr>=wr_period):
                             # store png plot to local folder
                             self.saveFigure(figpath, figfile)
                             # store png plot to FTP server according to RT-Web format
                             #name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
                             #ftp_filename = os.path.join(figpath, name)
                             #self.saveFigure(figpath, ftp_filename)
                             self.counter_imagwr = 0
             class Noise(Figure):
                 isConfig = None
                 __nsubplots = None
                 PREFIX = 'noise'
                 def __init__(self):
                     self.timerange = 24*60*60
                     self.isConfig = False
                     self.__nsubplots = 1
                     self.counter_imagwr = 0
                     self.WIDTH = 600
                     self.HEIGHT = 300
                     self.WIDTHPROF = 120
                     self.HEIGHTPROF = 0
                     self.xdata = None
                     self.ydata = None
                     self.PLOT_CODE = 17
                     self.FTP_WEI = None
                     self.EXP_CODE = None
                     self.SUB_EXP_CODE = None
                     self.PLOT_POS = None
                     self.figfile = None
                 def getSubplots(self):
                     ncol = 1
                     nrow = 1
                     return nrow, ncol
                 def openfile(self, filename):
                     f = open(filename,'w+')
                     f.write('\n\n')
                     f.write('JICAMARCA RADIO OBSERVATORY - Noise \n')
                     f.write('DD MM YYYY  HH MM SS   Channel0    Channel1    Channel2    Channel3\n\n' )
                     f.close()
                 def save_data(self, filename_phase, data, data_datetime):
                     f=open(filename_phase,'a')
                     timetuple_data = data_datetime.timetuple()
                     day = str(timetuple_data.tm_mday)
                     month = str(timetuple_data.tm_mon)
                     year = str(timetuple_data.tm_year)
                     hour = str(timetuple_data.tm_hour)
                     minute = str(timetuple_data.tm_min)
                     second = str(timetuple_data.tm_sec)
                     f.write(day+' '+month+' '+year+'  '+hour+' '+minute+' '+second+'   '+str(data[0])+'   '+str(data[1])+'   '+str(data[2])+'   '+str(data[3])+'\n')
                     f.close()
                 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
                     self.__showprofile = showprofile
                     self.nplots = nplots
                     ncolspan = 7
                     colspan = 6
                     self.__nsubplots = 2
                     self.createFigure(id = id,
                                       wintitle = wintitle,
                                       widthplot = self.WIDTH+self.WIDTHPROF,
                                       heightplot = self.HEIGHT+self.HEIGHTPROF,
                                       show=show)
                     nrow, ncol = self.getSubplots()
                     self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
                 def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True',
                         xmin=None, xmax=None, ymin=None, ymax=None,
                         timerange=None,
                         save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
                         server=None, folder=None, username=None, password=None,
                         ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
                     if channelList == None:
                         channelIndexList = dataOut.channelIndexList
                         channelList = dataOut.channelList
                     else:
                         channelIndexList = []
                         for channel in channelList:
                             if channel not in dataOut.channelList:
                                 raise ValueError, "Channel %d is not in dataOut.channelList"
                             channelIndexList.append(dataOut.channelList.index(channel))
                     if timerange != None:
                         self.timerange = timerange
                     tmin = None
                     tmax = None
                     x = dataOut.getTimeRange()
                     y = dataOut.getHeiRange()
                     factor = dataOut.normFactor
                     noise = dataOut.noise()/factor
                     noisedB = 10*numpy.log10(noise)
                     #thisDatetime = dataOut.datatime
                     thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
                     title = wintitle + " Noise" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
                     xlabel = ""
                     ylabel = "Intensity (dB)"
                     if not self.isConfig:
                         nplots = 1
                         self.setup(id=id,
                                    nplots=nplots,
                                    wintitle=wintitle,
                                    showprofile=showprofile,
                                    show=show)
                         tmin, tmax = self.getTimeLim(x, xmin, xmax)
                         if ymin == None: ymin = numpy.nanmin(noisedB) - 10.0
                         if ymax == None: ymax = numpy.nanmax(noisedB) + 10.0
                         self.FTP_WEI = ftp_wei
                         self.EXP_CODE = exp_code
                         self.SUB_EXP_CODE = sub_exp_code
                         self.PLOT_POS = plot_pos
                         self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         self.isConfig = True
                         self.figfile = figfile
                         self.xdata = numpy.array([])
                         self.ydata = numpy.array([])
                         #open file beacon phase
                         path = '%s%03d' %(self.PREFIX, self.id)
                         noise_file = os.path.join(path,'%s.txt'%self.name)
                         self.filename_noise = os.path.join(figpath,noise_file)
                         self.openfile(self.filename_noise)
                     #store data beacon phase
                     self.save_data(self.filename_noise, noisedB, thisDatetime)
                     self.setWinTitle(title)
                     title = "Noise %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
                     legendlabels = ["channel %d"%(idchannel+1) for idchannel in channelList]
                     axes = self.axesList[0]
                     self.xdata = numpy.hstack((self.xdata, x[0:1]))
                     if len(self.ydata)==0:
                         self.ydata = noisedB[channelIndexList].reshape(-1,1)
                     else:
                         self.ydata = numpy.hstack((self.ydata, noisedB[channelIndexList].reshape(-1,1)))
                     axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
                                 xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax,
                                 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
                                 XAxisAsTime=True, grid='both'
                                 )
                     self.draw()
                     if x[1] >= self.axesList[0].xmax:
                         self.counter_imagwr = wr_period
                         del self.xdata
                         del self.ydata
                         self.__isConfig = False
                     if self.figfile == None:
                         str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         self.figfile = self.getFilename(name = str_datetime)
                     if figpath != '':
                         self.counter_imagwr += 1
                         if (self.counter_imagwr>=wr_period):
                             # store png plot to local folder
                             self.saveFigure(figpath, self.figfile)
                             # store png plot to FTP server according to RT-Web format
                             name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
                             ftp_filename = os.path.join(figpath, name)
                             self.saveFigure(figpath, ftp_filename)
                             self.counter_imagwr = 0
             class BeaconPhase(Figure):
                 __isConfig = None
                 __nsubplots = None
                 PREFIX = 'beacon_phase'
                 def __init__(self):
                     self.timerange = 24*60*60
                     self.__isConfig = False
                     self.__nsubplots = 1
                     self.counter_imagwr = 0
                     self.WIDTH = 600
                     self.HEIGHT = 300
                     self.WIDTHPROF = 120
                     self.HEIGHTPROF = 0
                     self.xdata = None
                     self.ydata = None
                     self.PLOT_CODE = 18
                     self.FTP_WEI = None
                     self.EXP_CODE = None
                     self.SUB_EXP_CODE = None
                     self.PLOT_POS = None
                     self.filename_phase = None
                     self.figfile = None
                 def getSubplots(self):
                     ncol = 1
                     nrow = 1
                     return nrow, ncol
                 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
                     self.__showprofile = showprofile
                     self.nplots = nplots
                     ncolspan = 7
                     colspan = 6
                     self.__nsubplots = 2
                     self.createFigure(id = id,
                                       wintitle = wintitle,
                                       widthplot = self.WIDTH+self.WIDTHPROF,
                                       heightplot = self.HEIGHT+self.HEIGHTPROF,
                                       show=show)
                     nrow, ncol = self.getSubplots()
                     self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
                 def save_phase(self, filename_phase):
                     f = open(filename_phase,'w+')
                     f.write('\n\n')
                     f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
                     f.write('DD MM YYYY  HH MM SS   pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
                     f.close()
                 def save_data(self, filename_phase, data, data_datetime):
                     f=open(filename_phase,'a')
                     timetuple_data = data_datetime.timetuple()
                     day = str(timetuple_data.tm_mday)
                     month = str(timetuple_data.tm_mon)
                     year = str(timetuple_data.tm_year)
                     hour = str(timetuple_data.tm_hour)
                     minute = str(timetuple_data.tm_min)
                     second = str(timetuple_data.tm_sec)
                     f.write(day+' '+month+' '+year+'  '+hour+' '+minute+' '+second+'   '+str(data[0])+'   '+str(data[1])+'   '+str(data[2])+'   '+str(data[3])+'\n')
                     f.close()
                 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
                         xmin=None, xmax=None, ymin=None, ymax=None,
                         timerange=None,
                         save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1,
                         server=None, folder=None, username=None, password=None,
                         ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
                     if pairsList == None:
                         pairsIndexList = dataOut.pairsIndexList
                     else:
                         pairsIndexList = []
                         for pair in pairsList:
                             if pair not in dataOut.pairsList:
                                 raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair)
                             pairsIndexList.append(dataOut.pairsList.index(pair))
                     if pairsIndexList == []:
                         return
             #         if len(pairsIndexList) > 4:
             #             pairsIndexList = pairsIndexList[0:4]
                     if timerange != None:
                         self.timerange = timerange
                     tmin = None
                     tmax = None
                     x = dataOut.getTimeRange()
                     y = dataOut.getHeiRange()
                     #thisDatetime = dataOut.datatime
                     thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
                     title = wintitle + " Phase of Beacon Signal" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
                     xlabel = "Local Time"
                     ylabel = "Phase"
                     nplots = len(pairsIndexList)
                     #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
                     phase_beacon = numpy.zeros(len(pairsIndexList))
                     for i in range(nplots):
                         pair = dataOut.pairsList[pairsIndexList[i]]
                         ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i],:,:],axis=0)
                         powa = numpy.average(dataOut.data_spc[pair[0],:,:],axis=0)
                         powb = numpy.average(dataOut.data_spc[pair[1],:,:],axis=0)
                         avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
                         phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
                         #print "Phase %d%d" %(pair[0], pair[1])
                         #print phase[dataOut.beacon_heiIndexList]
                         phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
                     if not self.__isConfig:
                         nplots = len(pairsIndexList)
                         self.setup(id=id,
                                    nplots=nplots,
                                    wintitle=wintitle,
                                    showprofile=showprofile,
                                    show=show)
                         tmin, tmax = self.getTimeLim(x, xmin, xmax)
                         if ymin == None: ymin = numpy.nanmin(phase_beacon) - 10.0
                         if ymax == None: ymax = numpy.nanmax(phase_beacon) + 10.0
                         self.FTP_WEI = ftp_wei
                         self.EXP_CODE = exp_code
                         self.SUB_EXP_CODE = sub_exp_code
                         self.PLOT_POS = plot_pos
                         self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         self.__isConfig = True
                         self.figfile = figfile
                         self.xdata = numpy.array([])
                         self.ydata = numpy.array([])
                         #open file beacon phase
                         path = '%s%03d' %(self.PREFIX, self.id)
                         beacon_file = os.path.join(path,'%s.txt'%self.name)
                         self.filename_phase = os.path.join(figpath,beacon_file)
                         #self.save_phase(self.filename_phase)
                     #store data beacon phase
                     #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
                     self.setWinTitle(title)
                     title = "Beacon Signal %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
                     legendlabels = ["pairs %d%d"%(pair[0], pair[1]) for pair in dataOut.pairsList]
                     axes = self.axesList[0]
                     self.xdata = numpy.hstack((self.xdata, x[0:1]))
                     if len(self.ydata)==0:
                         self.ydata = phase_beacon.reshape(-1,1)
                     else:
                         self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
                     axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
                                 xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax,
                                 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
                                 XAxisAsTime=True, grid='both'
                                 )
                     self.draw()
                     if x[1] >= self.axesList[0].xmax:
                         self.counter_imagwr = wr_period
                         del self.xdata
                         del self.ydata
                         self.__isConfig = False
                     if self.figfile == None:
                         str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         self.figfile = self.getFilename(name = str_datetime)
                     if figpath != '':
                         self.counter_imagwr += 1
                         if (self.counter_imagwr>=wr_period):
                             # store png plot to local folder
                             self.saveFigure(figpath, self.figfile)
                             # store png plot to FTP server according to RT-Web format
                             name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
                             ftp_filename = os.path.join(figpath, name)
                             self.saveFigure(figpath, ftp_filename)
                             self.counter_imagwr = 0

schainpy/model/graphics/mpldriver.py +46 -2

             import numpy
             import datetime
             import sys
             import matplotlib
             if 'linux' in sys.platform:
                 matplotlib.use("TKAgg")
             if 'darwin' in sys.platform:
                 matplotlib.use("TKAgg")
             import matplotlib.pyplot
             from mpl_toolkits.axes_grid1 import make_axes_locatable
             from matplotlib.ticker import *
             ###########################################
             #Actualizacion de las funciones del driver
             ###########################################
             def createFigure(id, wintitle, width, height, facecolor="w", show=True):
                 matplotlib.pyplot.ioff()
                 fig = matplotlib.pyplot.figure(num=id, facecolor=facecolor)
                 fig.canvas.manager.set_window_title(wintitle)
                 fig.canvas.manager.resize(width, height)
                 matplotlib.pyplot.ion()
                 if show:
                     matplotlib.pyplot.show()
                 return fig
             def closeFigure(show=True):
                 matplotlib.pyplot.ioff()
                 if show:
                     matplotlib.pyplot.show()
                 return
             def saveFigure(fig, filename):
                 matplotlib.pyplot.ioff()
                 fig.savefig(filename)
                 matplotlib.pyplot.ion()
             def setWinTitle(fig, title):
                 fig.canvas.manager.set_window_title(title)
             def setTitle(fig, title):
                 fig.suptitle(title)
-            def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan):
+            def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan, polar=False):
                 matplotlib.pyplot.ioff()
                 matplotlib.pyplot.figure(fig.number)
                 axes = matplotlib.pyplot.subplot2grid((nrow, ncol),
                                                     (xpos, ypos),
                                                     colspan=colspan,
-                                                    rowspan=rowspan)
+                                                    rowspan=rowspan,
+                                                    polar=polar)
                 matplotlib.pyplot.ion()
                 return axes
             def setAxesText(ax, text):
                 ax.annotate(text,
                             xy = (.1, .99),
                             xycoords = 'figure fraction',
                             horizontalalignment = 'left',
                             verticalalignment = 'top',
                             fontsize = 10)
             def printLabels(ax, xlabel, ylabel, title):
                 ax.set_xlabel(xlabel, size=11)
                 ax.set_ylabel(ylabel, size=11)
                 ax.set_title(title, size=8)
             def createPline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='',
                             ticksize=9, xtick_visible=True, ytick_visible=True,
                             nxticks=4, nyticks=10,
                             grid=None,color='blue'):
                 """
                 Input:
                     grid    :    None, 'both', 'x', 'y'
                 """
                 matplotlib.pyplot.ioff()
                 ax.set_xlim([xmin,xmax])
                 ax.set_ylim([ymin,ymax])
                 printLabels(ax, xlabel, ylabel, title)
                 ######################################################
                 if (xmax-xmin)<=1:
                     xtickspos = numpy.linspace(xmin,xmax,nxticks)
                     xtickspos = numpy.array([float("%.1f"%i) for i in xtickspos])
                     ax.set_xticks(xtickspos)
                 else:
                     xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
             #         xtickspos = numpy.arange(nxticks)*float(xmax-xmin)/float(nxticks) + int(xmin)
                     ax.set_xticks(xtickspos)
                 for tick in ax.get_xticklabels():
                     tick.set_visible(xtick_visible)
                 for tick in ax.xaxis.get_major_ticks():
                     tick.label.set_fontsize(ticksize)
                 ######################################################
                 for tick in ax.get_yticklabels():
                     tick.set_visible(ytick_visible)
                 for tick in ax.yaxis.get_major_ticks():
                     tick.label.set_fontsize(ticksize)
                 ax.plot(x, y, color=color)
                 iplot = ax.lines[-1]
                 ######################################################
                 if '0.' in matplotlib.__version__[0:2]:
                     print "The matplotlib version has to be updated to 1.1 or newer"
                     return iplot
                 if '1.0.' in matplotlib.__version__[0:4]:
                     print "The matplotlib version has to be updated to 1.1 or newer"
                     return iplot
                 if grid != None:
                     ax.grid(b=True, which='major', axis=grid)
                 matplotlib.pyplot.tight_layout()
                 matplotlib.pyplot.ion()
                 return iplot
             def set_linedata(ax, x, y, idline):
                 ax.lines[idline].set_data(x,y)
             def pline(iplot, x, y, xlabel='', ylabel='', title=''):
                 ax = iplot.get_axes()
                 printLabels(ax, xlabel, ylabel, title)
                 set_linedata(ax, x, y, idline=0)
             def addpline(ax, x, y, color, linestyle, lw):
                 ax.plot(x,y,color=color,linestyle=linestyle,lw=lw)
             def createPcolor(ax, x, y, z, xmin, xmax, ymin, ymax, zmin, zmax,
                              xlabel='', ylabel='', title='', ticksize = 9,
                              colormap='jet',cblabel='', cbsize="5%",
                              XAxisAsTime=False):
                 matplotlib.pyplot.ioff()
                 divider = make_axes_locatable(ax)
                 ax_cb = divider.new_horizontal(size=cbsize, pad=0.05)
                 fig = ax.get_figure()
                 fig.add_axes(ax_cb)
                 ax.set_xlim([xmin,xmax])
                 ax.set_ylim([ymin,ymax])
                 printLabels(ax, xlabel, ylabel, title)
                 imesh = ax.pcolormesh(x,y,z.T, vmin=zmin, vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
                 cb =  matplotlib.pyplot.colorbar(imesh, cax=ax_cb)
                 cb.set_label(cblabel)
             #    for tl in ax_cb.get_yticklabels():
             #        tl.set_visible(True)
                 for tick in ax.yaxis.get_major_ticks():
                     tick.label.set_fontsize(ticksize)
                 for tick in ax.xaxis.get_major_ticks():
                     tick.label.set_fontsize(ticksize)
                 for tick in cb.ax.get_yticklabels():
                     tick.set_fontsize(ticksize)
                 ax_cb.yaxis.tick_right()
                 if '0.' in matplotlib.__version__[0:2]:
                     print "The matplotlib version has to be updated to 1.1 or newer"
                     return imesh
                 if '1.0.' in matplotlib.__version__[0:4]:
                     print "The matplotlib version has to be updated to 1.1 or newer"
                     return imesh
                 matplotlib.pyplot.tight_layout()
                 if XAxisAsTime:
                     func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
                     ax.xaxis.set_major_formatter(FuncFormatter(func))
                     ax.xaxis.set_major_locator(LinearLocator(7))
                 matplotlib.pyplot.ion()
                 return imesh
             def pcolor(imesh, z, xlabel='', ylabel='', title=''):
                 z = z.T
                 ax = imesh.get_axes()
                 printLabels(ax, xlabel, ylabel, title)
                 imesh.set_array(z.ravel())
             def addpcolor(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
                 printLabels(ax, xlabel, ylabel, title)
                 ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
             def addpcolorbuffer(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
                 printLabels(ax, xlabel, ylabel, title)
                 ax.collections.remove(ax.collections[0])
                 ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
             def createPmultiline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
                             ticksize=9, xtick_visible=True, ytick_visible=True,
                             nxticks=4, nyticks=10,
                             grid=None):
                 """
                 Input:
                     grid    :    None, 'both', 'x', 'y'
                 """
                 matplotlib.pyplot.ioff()
                 lines = ax.plot(x.T, y)
                 leg = ax.legend(lines, legendlabels, loc='upper right')
                 leg.get_frame().set_alpha(0.5)
                 ax.set_xlim([xmin,xmax])
                 ax.set_ylim([ymin,ymax])
                 printLabels(ax, xlabel, ylabel, title)
                 xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
                 ax.set_xticks(xtickspos)
                 for tick in ax.get_xticklabels():
                     tick.set_visible(xtick_visible)
                 for tick in ax.xaxis.get_major_ticks():
                     tick.label.set_fontsize(ticksize)
                 for tick in ax.get_yticklabels():
                     tick.set_visible(ytick_visible)
                 for tick in ax.yaxis.get_major_ticks():
                     tick.label.set_fontsize(ticksize)
                 iplot = ax.lines[-1]
                 if '0.' in matplotlib.__version__[0:2]:
                     print "The matplotlib version has to be updated to 1.1 or newer"
                     return iplot
                 if '1.0.' in matplotlib.__version__[0:4]:
                     print "The matplotlib version has to be updated to 1.1 or newer"
                     return iplot
                 if grid != None:
                     ax.grid(b=True, which='major', axis=grid)
                 matplotlib.pyplot.tight_layout()
                 matplotlib.pyplot.ion()
                 return iplot
             def pmultiline(iplot, x, y, xlabel='', ylabel='', title=''):
                 ax = iplot.get_axes()
                 printLabels(ax, xlabel, ylabel, title)
                 for i in range(len(ax.lines)):
                     line = ax.lines[i]
                     line.set_data(x[i,:],y)
             def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
                             ticksize=9, xtick_visible=True, ytick_visible=True,
                             nxticks=4, nyticks=10, marker='.', markersize=10, linestyle="None",
                             grid=None, XAxisAsTime=False):
                 """
                 Input:
                     grid    :    None, 'both', 'x', 'y'
                 """
                 matplotlib.pyplot.ioff()
             #    lines = ax.plot(x, y.T, marker=marker,markersize=markersize,linestyle=linestyle)
                 lines = ax.plot(x, y.T, linestyle='None', marker='.', markersize=markersize)
                 leg = ax.legend(lines, legendlabels, loc='upper left', bbox_to_anchor=(1.01, 1.00), numpoints=1, handlelength=1.5, \
                                 handletextpad=0.5, borderpad=0.5, labelspacing=0.5, borderaxespad=0.)
                 for label in leg.get_texts(): label.set_fontsize(9)
                 ax.set_xlim([xmin,xmax])
                 ax.set_ylim([ymin,ymax])
                 printLabels(ax, xlabel, ylabel, title)
             #    xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
             #    ax.set_xticks(xtickspos)
                 for tick in ax.get_xticklabels():
                     tick.set_visible(xtick_visible)
                 for tick in ax.xaxis.get_major_ticks():
                     tick.label.set_fontsize(ticksize)
                 for tick in ax.get_yticklabels():
                     tick.set_visible(ytick_visible)
                 for tick in ax.yaxis.get_major_ticks():
                     tick.label.set_fontsize(ticksize)
                 iplot = ax.lines[-1]
                 if '0.' in matplotlib.__version__[0:2]:
                     print "The matplotlib version has to be updated to 1.1 or newer"
                     return iplot
                 if '1.0.' in matplotlib.__version__[0:4]:
                     print "The matplotlib version has to be updated to 1.1 or newer"
                     return iplot
                 if grid != None:
                     ax.grid(b=True, which='major', axis=grid)
                 matplotlib.pyplot.tight_layout()
                 if XAxisAsTime:
                     func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
                     ax.xaxis.set_major_formatter(FuncFormatter(func))
                     ax.xaxis.set_major_locator(LinearLocator(7))
                 matplotlib.pyplot.ion()
                 return iplot
             def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''):
                 ax = iplot.get_axes()
                 printLabels(ax, xlabel, ylabel, title)
                 for i in range(len(ax.lines)):
                     line = ax.lines[i]
                     line.set_data(x,y[i,:])
+            def createPolar(ax, x, y,
+                            xlabel='', ylabel='', title='', ticksize = 9,
+                             colormap='jet',cblabel='', cbsize="5%",
+                             XAxisAsTime=False):
+                matplotlib.pyplot.ioff()
+                ax.plot(x,y,'bo', markersize=5)
+            #     ax.set_rmax(90)
+                ax.set_ylim(0,90)
+                ax.set_yticks(numpy.arange(0,90,20))
+                ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11')
+            #     ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical')
+                printLabels(ax, xlabel, '', title)
+                iplot = ax.lines[-1]
+                if '0.' in matplotlib.__version__[0:2]:
+                    print "The matplotlib version has to be updated to 1.1 or newer"
+                    return iplot
+                if '1.0.' in matplotlib.__version__[0:4]:
+                    print "The matplotlib version has to be updated to 1.1 or newer"
+                    return iplot
+            #     if grid != None:
+            #         ax.grid(b=True, which='major', axis=grid)
+                matplotlib.pyplot.tight_layout()
+                matplotlib.pyplot.ion()
+                return iplot
+            def polar(iplot, x, y, xlabel='', ylabel='', title=''):
+                ax = iplot.get_axes()
+            #     ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11')
+                printLabels(ax, xlabel, '', title)
+                set_linedata(ax, x, y, idline=0)
             def draw(fig):
                 if type(fig) == 'int':
                     raise ValueError, "This parameter should be of tpye matplotlib figure"
                 fig.canvas.draw()

schainpy/model/proc/jroproc_spectra.py +1 0

             import numpy
+            import math
             from jroproc_base import ProcessingUnit, Operation
             from model.data.jrodata import Spectra
             from model.data.jrodata import hildebrand_sekhon
             class SpectraProc(ProcessingUnit):
                 def __init__(self):
                     ProcessingUnit.__init__(self)
                     self.buffer = None
                     self.firstdatatime = None
                     self.profIndex = 0
                     self.dataOut = Spectra()
                     self.id_min = None
                     self.id_max = None
                 def __updateObjFromInput(self):
                     self.dataOut.timeZone = self.dataIn.timeZone
                     self.dataOut.dstFlag = self.dataIn.dstFlag
                     self.dataOut.errorCount = self.dataIn.errorCount
                     self.dataOut.useLocalTime = self.dataIn.useLocalTime
                     self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
                     self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
                     self.dataOut.channelList = self.dataIn.channelList
                     self.dataOut.heightList = self.dataIn.heightList
                     self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
             #        self.dataOut.nHeights = self.dataIn.nHeights
             #        self.dataOut.nChannels = self.dataIn.nChannels
                     self.dataOut.nBaud = self.dataIn.nBaud
                     self.dataOut.nCode = self.dataIn.nCode
                     self.dataOut.code = self.dataIn.code
                     self.dataOut.nProfiles = self.dataOut.nFFTPoints
             #        self.dataOut.channelIndexList = self.dataIn.channelIndexList
                     self.dataOut.flagTimeBlock = self.dataIn.flagTimeBlock
                     self.dataOut.utctime = self.firstdatatime
                     self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
                     self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
             #        self.dataOut.flagShiftFFT = self.dataIn.flagShiftFFT
                     self.dataOut.nCohInt = self.dataIn.nCohInt
                     self.dataOut.nIncohInt = 1
             #         self.dataOut.ippSeconds = self.dataIn.ippSeconds
                     self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
                     self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nFFTPoints*self.dataOut.nIncohInt
                     self.dataOut.frequency = self.dataIn.frequency
                     self.dataOut.realtime = self.dataIn.realtime
                     self.dataOut.azimuth = self.dataIn.azimuth
                     self.dataOut.zenith = self.dataIn.zenith
                     self.dataOut.beam.codeList = self.dataIn.beam.codeList
                     self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
                     self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
                 def __getFft(self):
                     """
                     Convierte valores de Voltaje a Spectra
                     Affected:
                         self.dataOut.data_spc
                         self.dataOut.data_cspc
                         self.dataOut.data_dc
                         self.dataOut.heightList
                         self.profIndex
                         self.buffer
                         self.dataOut.flagNoData
                     """
                     fft_volt = numpy.fft.fft(self.buffer,n=self.dataOut.nFFTPoints,axis=1)
                     fft_volt = fft_volt.astype(numpy.dtype('complex'))
                     dc = fft_volt[:,0,:]
                     #calculo de self-spectra
                     fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
                     spc = fft_volt * numpy.conjugate(fft_volt)
                     spc = spc.real
                     blocksize = 0
                     blocksize += dc.size
                     blocksize += spc.size
                     cspc = None
                     pairIndex = 0
                     if self.dataOut.pairsList != None:
                         #calculo de cross-spectra
                         cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
                         for pair in self.dataOut.pairsList:
                             cspc[pairIndex,:,:] = fft_volt[pair[0],:,:] * numpy.conjugate(fft_volt[pair[1],:,:])
                             pairIndex += 1
                         blocksize += cspc.size
                     self.dataOut.data_spc = spc
                     self.dataOut.data_cspc = cspc
                     self.dataOut.data_dc = dc
                     self.dataOut.blockSize = blocksize
                     self.dataOut.flagShiftFFT = False
                 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], ippFactor=None):
                     self.dataOut.flagNoData = True
                     if self.dataIn.type == "Spectra":
                         self.dataOut.copy(self.dataIn)
                         return True
                     if self.dataIn.type == "Voltage":
                         if nFFTPoints == None:
                             raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
                         if nProfiles == None:
                             raise ValueError, "This SpectraProc.run() need nProfiles input variable"
                         if ippFactor == None:
                             ippFactor = 1
                         self.dataOut.ippFactor = ippFactor
                         self.dataOut.nFFTPoints = nFFTPoints
                         self.dataOut.pairsList = pairsList
                         if self.buffer == None:
                             self.buffer = numpy.zeros((self.dataIn.nChannels,
                                                        nProfiles,
                                                        self.dataIn.nHeights),
                                                     dtype='complex')
                             self.id_min = 0
                             self.id_max = self.dataIn.data.shape[1]
                         if len(self.dataIn.data.shape) == 2:
                             self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
                             self.profIndex += 1
                         else:
                             if self.dataIn.data.shape[1] == nProfiles:
                                 self.buffer = self.dataIn.data.copy()
                                 self.profIndex = nProfiles
                             elif self.dataIn.data.shape[1] < nProfiles:
                                 self.buffer[:,self.id_min:self.id_max,:] = self.dataIn.data
                                 self.profIndex += self.dataIn.data.shape[1]
                                 self.id_min += self.dataIn.data.shape[1]
                                 self.id_max += self.dataIn.data.shape[1]
                             else:
                                 raise ValueError, "The type object %s has %d profiles, it should be equal to %d profiles"%(self.dataIn.type,self.dataIn.data.shape[1],nProfiles)
                                 self.dataOut.flagNoData = True
                                 return 0
                         if self.firstdatatime == None:
                             self.firstdatatime = self.dataIn.utctime
                         if self.profIndex == nProfiles:
                             self.__updateObjFromInput()
                             self.__getFft()
                             self.dataOut.flagNoData = False
                             self.buffer = None
                             self.firstdatatime = None
                             self.profIndex = 0
                         return True
                     raise ValueError, "The type object %s is not valid"%(self.dataIn.type)
                 def selectChannels(self, channelList):
                     channelIndexList = []
                     for channel in channelList:
                         index = self.dataOut.channelList.index(channel)
                         channelIndexList.append(index)
                     self.selectChannelsByIndex(channelIndexList)
                 def selectChannelsByIndex(self, channelIndexList):
                     """
                     Selecciona un bloque de datos en base a canales segun el channelIndexList
                     Input:
                         channelIndexList    :    lista sencilla de canales a seleccionar por ej. [2,3,7]
                     Affected:
                         self.dataOut.data_spc
                         self.dataOut.channelIndexList
                         self.dataOut.nChannels
                     Return:
                         None
                     """
                     for channelIndex in channelIndexList:
                         if channelIndex not in self.dataOut.channelIndexList:
                             print channelIndexList
                             raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
             #         nChannels = len(channelIndexList)
                     data_spc = self.dataOut.data_spc[channelIndexList,:]
                     self.dataOut.data_spc = data_spc
                     self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
             #        self.dataOut.nChannels = nChannels
                     return 1
                 def selectHeights(self, minHei, maxHei):
                     """
                     Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
                     minHei <= height <= maxHei
                     Input:
                         minHei    :    valor minimo de altura a considerar
                         maxHei    :    valor maximo de altura a considerar
                     Affected:
                         Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
                     Return:
 si el metodo se ejecuto con exito caso contrario devuelve 0
                     """
                     if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
                         raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
                     if (maxHei > self.dataOut.heightList[-1]):
                         maxHei = self.dataOut.heightList[-1]
             #            raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
                     minIndex = 0
                     maxIndex = 0
                     heights = self.dataOut.heightList
                     inda = numpy.where(heights >= minHei)
                     indb = numpy.where(heights <= maxHei)
                     try:
                         minIndex = inda[0][0]
                     except:
                         minIndex = 0
                     try:
                         maxIndex = indb[0][-1]
                     except:
                         maxIndex = len(heights)
                     self.selectHeightsByIndex(minIndex, maxIndex)
                     return 1
                 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
                     newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
                     if hei_ref != None:
                         newheis = numpy.where(self.dataOut.heightList>hei_ref)
                     minIndex = min(newheis[0])
                     maxIndex = max(newheis[0])
                     data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
                     heightList = self.dataOut.heightList[minIndex:maxIndex+1]
                     # determina indices
                     nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
                     avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
                     beacon_dB = numpy.sort(avg_dB)[-nheis:]
                     beacon_heiIndexList = []
                     for val in avg_dB.tolist():
                         if val >= beacon_dB[0]:
                             beacon_heiIndexList.append(avg_dB.tolist().index(val))
                     #data_spc = data_spc[:,:,beacon_heiIndexList]
                     data_cspc = None
                     if self.dataOut.data_cspc != None:
                         data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
                         #data_cspc = data_cspc[:,:,beacon_heiIndexList]
                     data_dc = None
                     if self.dataOut.data_dc != None:
                         data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
                         #data_dc = data_dc[:,beacon_heiIndexList]
                     self.dataOut.data_spc = data_spc
                     self.dataOut.data_cspc = data_cspc
                     self.dataOut.data_dc = data_dc
                     self.dataOut.heightList = heightList
                     self.dataOut.beacon_heiIndexList = beacon_heiIndexList
                     return 1
                 def selectHeightsByIndex(self, minIndex, maxIndex):
                     """
                     Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
                     minIndex <= index <= maxIndex
                     Input:
                         minIndex    :    valor de indice minimo de altura a considerar
                         maxIndex    :    valor de indice maximo de altura a considerar
                     Affected:
                         self.dataOut.data_spc
                         self.dataOut.data_cspc
                         self.dataOut.data_dc
                         self.dataOut.heightList
                     Return:
 si el metodo se ejecuto con exito caso contrario devuelve 0
                     """
                     if (minIndex < 0) or (minIndex > maxIndex):
                         raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
                     if (maxIndex >= self.dataOut.nHeights):
                         maxIndex = self.dataOut.nHeights-1
             #            raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
             #         nHeights = maxIndex - minIndex + 1
                     #Spectra
                     data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
                     data_cspc = None
                     if self.dataOut.data_cspc != None:
                         data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
                     data_dc = None
                     if self.dataOut.data_dc != None:
                         data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
                     self.dataOut.data_spc = data_spc
                     self.dataOut.data_cspc = data_cspc
                     self.dataOut.data_dc = data_dc
                     self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
                     return 1
                 def removeDC(self, mode = 2):
                     jspectra = self.dataOut.data_spc
                     jcspectra = self.dataOut.data_cspc
                     num_chan = jspectra.shape[0]
                     num_hei = jspectra.shape[2]
                     if jcspectra != None:
                         jcspectraExist = True
                         num_pairs = jcspectra.shape[0]
                     else:   jcspectraExist = False
                     freq_dc = jspectra.shape[1]/2
                     ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
                     if ind_vel[0]<0:
                         ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
                     if mode == 1:
                         jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
                         if jcspectraExist:
                             jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
                     if mode == 2:
                         vel = numpy.array([-2,-1,1,2])
                         xx = numpy.zeros([4,4])
                         for fil in range(4):
                             xx[fil,:] = vel[fil]**numpy.asarray(range(4))
                         xx_inv = numpy.linalg.inv(xx)
                         xx_aux = xx_inv[0,:]
                         for ich in range(num_chan):
                             yy = jspectra[ich,ind_vel,:]
                             jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
                             junkid = jspectra[ich,freq_dc,:]<=0
                             cjunkid = sum(junkid)
                             if cjunkid.any():
                                 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
                         if jcspectraExist:
                             for ip in range(num_pairs):
                                 yy = jcspectra[ip,ind_vel,:]
                                 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
                     self.dataOut.data_spc = jspectra
                     self.dataOut.data_cspc = jcspectra
                     return 1
                 def removeInterference(self,  interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
                     jspectra = self.dataOut.data_spc
                     jcspectra = self.dataOut.data_cspc
                     jnoise = self.dataOut.getNoise()
                     num_incoh = self.dataOut.nIncohInt
                     num_channel  = jspectra.shape[0]
                     num_prof  = jspectra.shape[1]
                     num_hei   = jspectra.shape[2]
                     #hei_interf
                     if hei_interf == None:
                         count_hei = num_hei/2   #Como es entero no importa
                         hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
                         hei_interf = numpy.asarray(hei_interf)[0]
                     #nhei_interf
                     if (nhei_interf == None):
                         nhei_interf = 5
                     if (nhei_interf < 1):
                         nhei_interf = 1
                     if (nhei_interf > count_hei):
                         nhei_interf = count_hei
                     if (offhei_interf == None):
                         offhei_interf = 0
                     ind_hei = range(num_hei)
             #         mask_prof = numpy.asarray(range(num_prof - 2)) + 1
             #         mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
                     mask_prof = numpy.asarray(range(num_prof))
                     num_mask_prof = mask_prof.size
                     comp_mask_prof = [0, num_prof/2]
                     #noise_exist:    Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
                     if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
                         jnoise = numpy.nan
                     noise_exist = jnoise[0] < numpy.Inf
                     #Subrutina de Remocion de la Interferencia
                     for ich in range(num_channel):
                         #Se ordena los espectros segun su potencia (menor a mayor)
                         power = jspectra[ich,mask_prof,:]
                         power = power[:,hei_interf]
                         power = power.sum(axis = 0)
                         psort = power.ravel().argsort()
                         #Se estima la interferencia promedio en los Espectros de Potencia empleando
                         junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
                         if noise_exist:
                         #    tmp_noise = jnoise[ich] / num_prof
                             tmp_noise = jnoise[ich]
                         junkspc_interf = junkspc_interf - tmp_noise
                         #junkspc_interf[:,comp_mask_prof] = 0
                         jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
                         jspc_interf = jspc_interf.transpose()
                         #Calculando el espectro de interferencia promedio
                         noiseid =  numpy.where(jspc_interf <= tmp_noise/ math.sqrt(num_incoh))
                         noiseid = noiseid[0]
                         cnoiseid = noiseid.size
                         interfid = numpy.where(jspc_interf > tmp_noise/ math.sqrt(num_incoh))
                         interfid = interfid[0]
                         cinterfid = interfid.size
                         if (cnoiseid > 0):   jspc_interf[noiseid] = 0
                         #Expandiendo los perfiles a limpiar
                         if (cinterfid > 0):
                             new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
                             new_interfid = numpy.asarray(new_interfid)
                             new_interfid = {x for x in new_interfid}
                             new_interfid = numpy.array(list(new_interfid))
                             new_cinterfid = new_interfid.size
                         else: new_cinterfid = 0
                         for ip in range(new_cinterfid):
                             ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
                             jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
                         jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
                         #Removiendo la interferencia del punto de mayor interferencia
                         ListAux = jspc_interf[mask_prof].tolist()
                         maxid = ListAux.index(max(ListAux))
                         if cinterfid > 0:
                             for ip in range(cinterfid*(interf == 2) - 1):
                                 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/math.sqrt(num_incoh))).nonzero()
                                 cind = len(ind)
                                 if (cind > 0):
                                     jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/math.sqrt(num_incoh))
                             ind = numpy.array([-2,-1,1,2])
                             xx = numpy.zeros([4,4])
                             for id1 in range(4):
                                 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
                             xx_inv = numpy.linalg.inv(xx)
                             xx = xx_inv[:,0]
                             ind = (ind + maxid + num_mask_prof)%num_mask_prof
                             yy = jspectra[ich,mask_prof[ind],:]
                             jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
                         indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/math.sqrt(num_incoh))).nonzero()
                         jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/math.sqrt(num_incoh))
                     #Remocion de Interferencia en el Cross Spectra
                     if jcspectra == None: return jspectra, jcspectra
                     num_pairs = jcspectra.size/(num_prof*num_hei)
                     jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
                     for ip in range(num_pairs):
                         #-------------------------------------------
                         cspower = numpy.abs(jcspectra[ip,mask_prof,:])
                         cspower = cspower[:,hei_interf]
                         cspower = cspower.sum(axis = 0)
                         cspsort = cspower.ravel().argsort()
                         junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
                         junkcspc_interf = junkcspc_interf.transpose()
                         jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
                         ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
                         median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
                         median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
                         junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
                         for iprof in range(num_prof):
                             ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
                             jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
                         #Removiendo la Interferencia
                         jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
                         ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
                         maxid = ListAux.index(max(ListAux))
                         ind = numpy.array([-2,-1,1,2])
                         xx = numpy.zeros([4,4])
                         for id1 in range(4):
                             xx[:,id1] = ind[id1]**numpy.asarray(range(4))
                         xx_inv = numpy.linalg.inv(xx)
                         xx = xx_inv[:,0]
                         ind = (ind + maxid + num_mask_prof)%num_mask_prof
                         yy = jcspectra[ip,mask_prof[ind],:]
                         jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
                     #Guardar Resultados
                     self.dataOut.data_spc = jspectra
                     self.dataOut.data_cspc = jcspectra
                     return 1
                 def setRadarFrequency(self, frequency=None):
                     if frequency != None:
                         self.dataOut.frequency = frequency
                     return 1
                 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
                     #validacion de rango
                     if minHei == None:
                         minHei = self.dataOut.heightList[0]
                     if maxHei == None:
                         maxHei = self.dataOut.heightList[-1]
                     if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
                         print 'minHei: %.2f is out of the heights range'%(minHei)
                         print 'minHei is setting to %.2f'%(self.dataOut.heightList[0])
                         minHei = self.dataOut.heightList[0]
                     if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
                         print 'maxHei: %.2f is out of the heights range'%(maxHei)
                         print 'maxHei is setting to %.2f'%(self.dataOut.heightList[-1])
                         maxHei = self.dataOut.heightList[-1]
                     # validacion de velocidades
                     velrange = self.dataOut.getVelRange(1)
                     if minVel == None:
                         minVel = velrange[0]
                     if maxVel == None:
                         maxVel = velrange[-1]
                     if (minVel < velrange[0]) or (minVel > maxVel):
                         print 'minVel: %.2f is out of the velocity range'%(minVel)
                         print 'minVel is setting to %.2f'%(velrange[0])
                         minVel = velrange[0]
                     if (maxVel > velrange[-1]) or (maxVel < minVel):
                         print 'maxVel: %.2f is out of the velocity range'%(maxVel)
                         print 'maxVel is setting to %.2f'%(velrange[-1])
                         maxVel = velrange[-1]
                     # seleccion de indices para rango
                     minIndex = 0
                     maxIndex = 0
                     heights = self.dataOut.heightList
                     inda = numpy.where(heights >= minHei)
                     indb = numpy.where(heights <= maxHei)
                     try:
                         minIndex = inda[0][0]
                     except:
                         minIndex = 0
                     try:
                         maxIndex = indb[0][-1]
                     except:
                         maxIndex = len(heights)
                     if (minIndex < 0) or (minIndex > maxIndex):
                         raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
                     if (maxIndex >= self.dataOut.nHeights):
                         maxIndex = self.dataOut.nHeights-1
                     # seleccion de indices para velocidades
                     indminvel = numpy.where(velrange >= minVel)
                     indmaxvel = numpy.where(velrange <= maxVel)
                     try:
                         minIndexVel = indminvel[0][0]
                     except:
                         minIndexVel = 0
                     try:
                         maxIndexVel = indmaxvel[0][-1]
                     except:
                         maxIndexVel = len(velrange)
                     #seleccion del espectro
                     data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
                     #estimacion de ruido
                     noise = numpy.zeros(self.dataOut.nChannels)
                     for channel in range(self.dataOut.nChannels):
                         daux = data_spc[channel,:,:]
                         noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
                     self.dataOut.noise_estimation = noise.copy()
                     return 1
             class IncohInt(Operation):
                 __profIndex = 0
                 __withOverapping  = False
                 __byTime = False
                 __initime = None
                 __lastdatatime = None
                 __integrationtime = None
                 __buffer_spc = None
                 __buffer_cspc = None
                 __buffer_dc = None
                 __dataReady = False
                 __timeInterval = None
                 n = None
                 def __init__(self):
                     Operation.__init__(self)
             #         self.isConfig = False
                 def setup(self, n=None, timeInterval=None, overlapping=False):
                     """
                     Set the parameters of the integration class.
                     Inputs:
                         n        :    Number of coherent integrations
                         timeInterval   :    Time of integration. If the parameter "n" is selected this one does not work
                         overlapping    :
                     """
                     self.__initime = None
                     self.__lastdatatime = 0
                     self.__buffer_spc = None
                     self.__buffer_cspc = None
                     self.__buffer_dc = None
                     self.__dataReady = False
                     if n == None and timeInterval == None:
                         raise ValueError, "n or timeInterval should be specified ..."
                     if n != None:
                         self.n = n
                         self.__byTime = False
                     else:
                         self.__integrationtime = timeInterval #if (type(timeInterval)!=integer) -> change this line
                         self.n = 9999
                         self.__byTime = True
                     if overlapping:
                         self.__withOverapping = True
                     else:
                         self.__withOverapping = False
                         self.__buffer_spc = 0
                         self.__buffer_cspc = 0
                         self.__buffer_dc = 0
                     self.__profIndex = 0
                 def putData(self, data_spc, data_cspc, data_dc):
                     """
                     Add a profile to the __buffer_spc and increase in one the __profileIndex
                     """
                     if not self.__withOverapping:
                         self.__buffer_spc += data_spc
                         if data_cspc == None:
                             self.__buffer_cspc = None
                         else:
                             self.__buffer_cspc += data_cspc
                         if data_dc == None:
                             self.__buffer_dc = None
                         else:
                             self.__buffer_dc += data_dc
                         self.__profIndex += 1
                         return
                     #Overlapping data
                     nChannels, nFFTPoints, nHeis = data_spc.shape
                     data_spc = numpy.reshape(data_spc, (1, nChannels, nFFTPoints, nHeis))
                     if data_cspc != None:
                         data_cspc = numpy.reshape(data_cspc, (1, -1, nFFTPoints, nHeis))
                     if data_dc != None:
                         data_dc = numpy.reshape(data_dc, (1, -1, nHeis))
                     #If the buffer is empty then it takes the data value
                     if self.__buffer_spc == None:
                         self.__buffer_spc = data_spc
                         if data_cspc == None:
                             self.__buffer_cspc = None
                         else:
                             self.__buffer_cspc += data_cspc
                         if data_dc == None:
                             self.__buffer_dc = None
                         else:
                             self.__buffer_dc += data_dc
                         self.__profIndex += 1
                         return
                     #If the buffer length is lower than n then stakcing the data value
                     if self.__profIndex < self.n:
                         self.__buffer_spc = numpy.vstack((self.__buffer_spc, data_spc))
                         if data_cspc != None:
                             self.__buffer_cspc = numpy.vstack((self.__buffer_cspc, data_cspc))
                         if data_dc != None:
                             self.__buffer_dc = numpy.vstack((self.__buffer_dc, data_dc))
                         self.__profIndex += 1
                         return
                     #If the buffer length is equal to n then replacing the last buffer value with the data value
                     self.__buffer_spc = numpy.roll(self.__buffer_spc, -1, axis=0)
                     self.__buffer_spc[self.n-1] = data_spc
                     if data_cspc != None:
                         self.__buffer_cspc = numpy.roll(self.__buffer_cspc, -1, axis=0)
                         self.__buffer_cspc[self.n-1] = data_cspc
                     if data_dc != None:
                         self.__buffer_dc = numpy.roll(self.__buffer_dc, -1, axis=0)
                         self.__buffer_dc[self.n-1] = data_dc
                     self.__profIndex = self.n
                     return
                 def pushData(self):
                     """
                     Return the sum of the last profiles and the profiles used in the sum.
                     Affected:
                     self.__profileIndex
                     """
                     data_spc = None
                     data_cspc  = None
                     data_dc = None
                     if not self.__withOverapping:
                         data_spc = self.__buffer_spc
                         data_cspc = self.__buffer_cspc
                         data_dc = self.__buffer_dc
                         n = self.__profIndex
                         self.__buffer_spc = 0
                         self.__buffer_cspc = 0
                         self.__buffer_dc = 0
                         self.__profIndex = 0
                         return data_spc, data_cspc, data_dc, n
                     #Integration with Overlapping
                     data_spc = numpy.sum(self.__buffer_spc, axis=0)
                     if self.__buffer_cspc != None:
                         data_cspc = numpy.sum(self.__buffer_cspc, axis=0)
                     if self.__buffer_dc != None:
                         data_dc = numpy.sum(self.__buffer_dc, axis=0)
                     n = self.__profIndex
                     return data_spc, data_cspc, data_dc, n
                 def byProfiles(self, *args):
                     self.__dataReady = False
                     avgdata_spc = None
                     avgdata_cspc = None
                     avgdata_dc = None
             #         n = None
                     self.putData(*args)
                     if self.__profIndex == self.n:
                         avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
                         self.__dataReady = True
                     return avgdata_spc, avgdata_cspc, avgdata_dc
                 def byTime(self, datatime, *args):
                     self.__dataReady = False
                     avgdata_spc = None
                     avgdata_cspc = None
                     avgdata_dc = None
                     n = None
                     self.putData(*args)
                     if (datatime - self.__initime) >= self.__integrationtime:
                         avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
                         self.n = n
                         self.__dataReady = True
                     return avgdata_spc, avgdata_cspc, avgdata_dc
                 def integrate(self, datatime, *args):
                     if self.__initime == None:
                         self.__initime = datatime
                     if self.__byTime:
                         avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(datatime, *args)
                     else:
                         avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
                     self.__lastdatatime = datatime
                     if avgdata_spc == None:
                         return None, None, None, None
                     avgdatatime = self.__initime
                     try:
                         self.__timeInterval = (self.__lastdatatime - self.__initime)/(self.n - 1)
                     except:
                         self.__timeInterval = self.__lastdatatime - self.__initime
                     deltatime = datatime -self.__lastdatatime
                     if not self.__withOverapping:
                         self.__initime = datatime
                     else:
                         self.__initime += deltatime
                     return avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc
                 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
                     if n==1:
                         dataOut.flagNoData = False
                         return
                     if not self.isConfig:
                         self.setup(n, timeInterval, overlapping)
                         self.isConfig = True
                     avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
                                                                                         dataOut.data_spc,
                                                                                         dataOut.data_cspc,
                                                                                         dataOut.data_dc)
             #        dataOut.timeInterval *= n
                     dataOut.flagNoData = True
                     if self.__dataReady:
                         dataOut.data_spc = avgdata_spc
                         dataOut.data_cspc = avgdata_cspc
                         dataOut.data_dc = avgdata_dc
                         dataOut.nIncohInt *= self.n
                         dataOut.utctime = avgdatatime
                         #dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt * dataOut.nIncohInt * dataOut.nFFTPoints
                         dataOut.timeInterval = self.__timeInterval*self.n
                         dataOut.flagNoData = False

schainpy/model/proc/jroprocessing.py +3 -1

             from jroproc_voltage import *
             from jroproc_spectra import *
             from jroproc_heispectra import *
-            from jroproc_amisr import *
  No newline at end of file
+            from jroproc_amisr import *
+            from jroproc_correlation import *
+            from jroproc_parameters import *
  No newline at end of file

schainpy/test/Meteor_JASMET30Mhz_Beacon.py ~~schainpy/test/Meteor_JASMET_30mhz.py~~ renamed +20 -17

             """
             Se debe verficar que el disco de datos se encuentra montado en el sistema
             """
             import os, sys
             path = os.path.split(os.getcwd())[0]
             sys.path.append(path)
             from controller import *
             desc = "Meteor Experiment Test"
             filename = "meteor20130812.xml"
             controllerObj = Project()
             controllerObj.setup(id = '191', name='meteor_test01', description=desc)
             path = '/home/dsuarez/.gvfs/data on 10.10.20.13/Jasmet50'
+            pathFigure = '/home/jasmet/jasmet30_phase'
+            path = '/home/soporte/Data/JASMET/JASMET_30/2014106'
+            pathFigure = '/home/soporte/workspace/Graficos/JASMET/prueba1'
-            readUnitConfObj = controllerObj.addReadUnit(datatype='Voltage',
+            readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
                                                         path=path,
                                                         startDate='2014/04/15',
                                                         endDate='2014/04/15',
-                                                        startTime='17:00:00',
+                                                        startTime='20:00:00',
                                                         endTime='23:59:59',
                                                         online=0,
-                                                        walk=1)
+                                                        walk=0)
             opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
-            procUnitConfObj0 = controllerObj.addProcUnit(datatype='Voltage', inputId=readUnitConfObj.getId())
+            procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
-            opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
+            # opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
+            #
+            #
+            # opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other')
+            # opObj11.addParameter(name='n', value='2', format='int')
+            # opObj11 = procUnitConfObj0.addOperation(name='VoltageWriter', optype='other')
-            opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other')
+            # opObj11.addParameter(name='path', value='/home/jasmet/jasmet30_abril')
-            opObj11.addParameter(name='n', value='2', format='int')
+            # opObj11.addParameter(name='blocksPerFile', value='100', format='int')
+            # opObj11.addParameter(name='profilesPerBlock', value='200', format='int')
-            opObj11 = procUnitConfObj0.addOperation(name='VoltageWriter', optype='other')
-            opObj11.addParameter(name='path', value='/home/jasmet/jasmet30_abril')
-            opObj11.addParameter(name='blocksPerFile', value='100', format='int')
-            opObj11.addParameter(name='profilesPerBlock', value='200', format='int')
             """
             ########################################### BEACON ##########################################
             """
-            procUnitConfObjBeacon = controllerObj.addProcUnit(datatype='Spectra', inputId=procUnitConfObj0.getId())
+            procUnitConfObjBeacon = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId())
             procUnitConfObjBeacon.addParameter(name='nProfiles', value='200', format='int')
             procUnitConfObjBeacon.addParameter(name='nFFTPoints', value='200', format='int')
             procUnitConfObjBeacon.addParameter(name='pairsList', value='(2,0),(2,1),(2,3),(2,4)', format='pairsList')
             opObj11 = procUnitConfObjBeacon.addOperation(name='IncohInt', optype='other')
             opObj11.addParameter(name='n', value='4', format='int')
             opObj11 = procUnitConfObjBeacon.addOperation(name='getBeaconSignal')
             opObj11 = procUnitConfObjBeacon.addOperation(name='BeaconPhase', optype='other')
             opObj11.addParameter(name='id', value='201', format='int')
             opObj11.addParameter(name='wintitle', value='Beacon Phase', format='str')
             opObj11.addParameter(name='timerange', value='300', format='int')
-            opObj11.addParameter(name='xmin', value='0', format='float')
+            opObj11.addParameter(name='xmin', value='20', format='float')
-            opObj11.addParameter(name='xmax', value='24', format='float')
+            opObj11.addParameter(name='xmax', value='22', format='float')
             opObj11.addParameter(name='ymin', value='-180', format='float')
             opObj11.addParameter(name='ymax', value='180', format='float')
-            opObj11.addParameter(name='figpath', value='/home/jasmet/jasmet30_phase', format='str')
+            opObj11.addParameter(name='figpath', value=pathFigure, format='str')
             print "Escribiendo el archivo XML"
             controllerObj.writeXml(filename)
             print "Leyendo el archivo XML"
             controllerObj.readXml(filename)
             controllerObj.createObjects()
             controllerObj.connectObjects()
             controllerObj.run()

schainpy/test/Meteor_JASMET50Mhz_Beacon.py ~~schainpy/test/Meteor_JASMET_50mhz.py~~ renamed 0 0

NO CONTENT: file renamed from schainpy/test/Meteor_JASMET_50mhz.py to schainpy/test/Meteor_JASMET50Mhz_Beacon.py

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages