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

@@ -743,4 +743,225 class Fits:
743	noise = property(getNoise, "I'm the 'nHeights' property.")	743	noise = property(getNoise, "I'm the 'nHeights' property.")
744	datatime = property(getDatatime, "I'm the 'datatime' property")	744	datatime = property(getDatatime, "I'm the 'datatime' property")
745	ltctime = property(getltctime, "I'm the 'ltctime' property")	745	ltctime = property(getltctime, "I'm the 'ltctime' property")
		746
		747	class Correlation(JROData):
		748
		749	noise = None
		750
		751	SNR = None
		752
		753	pairsAutoCorr = None #Pairs of Autocorrelation
		754
		755	#--------------------------------------------------
		756
		757	data_corr = None
		758
		759	data_volt = None
		760
		761	lagT = None # each element value is a profileIndex
		762
		763	lagR = None # each element value is in km
		764
		765	pairsList = None
		766
		767	calculateVelocity = None
		768
		769	nPoints = None
		770
		771	nAvg = None
		772
		773	bufferSize = None
		774
		775	def __init__(self):
		776	'''
		777	Constructor
		778	'''
		779	self.radarControllerHeaderObj = RadarControllerHeader()
		780
		781	self.systemHeaderObj = SystemHeader()
		782
		783	self.type = "Correlation"
		784
		785	self.data = None
		786
		787	self.dtype = None
		788
		789	self.nProfiles = None
		790
		791	self.heightList = None
		792
		793	self.channelList = None
		794
		795	self.flagNoData = True
		796
		797	self.flagTimeBlock = False
		798
		799	self.utctime = None
		800
		801	self.timeZone = None
		802
		803	self.dstFlag = None
		804
		805	self.errorCount = None
		806
		807	self.blocksize = None
		808
		809	self.flagDecodeData = False #asumo q la data no esta decodificada
		810
		811	self.flagDeflipData = False #asumo q la data no esta sin flip
		812
		813	self.pairsList = None
		814
		815	self.nPoints = None
		816
		817	def getLagTRange(self, extrapoints=0):
		818
		819	lagTRange = self.lagT
		820	diff = lagTRange[1] - lagTRange[0]
		821	extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1]
		822	lagTRange = numpy.hstack((lagTRange, extra))
		823
		824	return lagTRange
		825
		826	def getLagRRange(self, extrapoints=0):
		827
		828	return self.lagR
		829
		830	def getPairsList(self):
		831
		832	return self.pairsList
		833
		834	def getCalculateVelocity(self):
		835
		836	return self.calculateVelocity
		837
		838	def getNPoints(self):
		839
		840	return self.nPoints
		841
		842	def getNAvg(self):
		843
		844	return self.nAvg
		845
		846	def getBufferSize(self):
		847
		848	return self.bufferSize
		849
		850	def getPairsAutoCorr(self):
		851	pairsList = self.pairsList
		852	pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan
		853
		854	for l in range(len(pairsList)):
		855	firstChannel = pairsList[l][0]
		856	secondChannel = pairsList[l][1]
		857
		858	#Obteniendo pares de Autocorrelacion
		859	if firstChannel == secondChannel:
		860	pairsAutoCorr[firstChannel] = int(l)
		861
		862	pairsAutoCorr = pairsAutoCorr.astype(int)
		863
		864	return pairsAutoCorr
		865
		866	def getNoise(self, mode = 2):
		867
		868	indR = numpy.where(self.lagR == 0)[0][0]
		869	indT = numpy.where(self.lagT == 0)[0][0]
		870
		871	jspectra0 = self.data_corr[:,:,indR,:]
		872	jspectra = copy.copy(jspectra0)
		873
		874	num_chan = jspectra.shape[0]
		875	num_hei = jspectra.shape[2]
		876
		877	freq_dc = jspectra.shape[1]/2
		878	ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
		879
		880	if ind_vel[0]<0:
		881	ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
		882
		883	if mode == 1:
		884	jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
		885
		886	if mode == 2:
		887
		888	vel = numpy.array([-2,-1,1,2])
		889	xx = numpy.zeros([4,4])
		890
		891	for fil in range(4):
		892	xx[fil,:] = vel[fil]**numpy.asarray(range(4))
		893
		894	xx_inv = numpy.linalg.inv(xx)
		895	xx_aux = xx_inv[0,:]
		896
		897	for ich in range(num_chan):
		898	yy = jspectra[ich,ind_vel,:]
		899	jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
		900
		901	junkid = jspectra[ich,freq_dc,:]<=0
		902	cjunkid = sum(junkid)
		903
		904	if cjunkid.any():
		905	jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
		906
		907	noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
		908
		909	return noise
		910
		911	# pairsList = property(getPairsList, "I'm the 'pairsList' property.")
		912	# nPoints = property(getNPoints, "I'm the 'nPoints' property.")
		913	calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.")
		914	nAvg = property(getNAvg, "I'm the 'nAvg' property.")
		915	bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.")
		916
		917
		918	class Parameters(JROData):
		919
		920	inputUnit = None #Type of data to be processed
		921
		922	operation = None #Type of operation to parametrize
		923
		924	data_param = None #Parameters obtained
		925
		926	data_pre = None #Data Pre Parametrization
746		927
		928	heightRange = None #Heights
		929
		930	abscissaRange = None #Abscissa, can be velocities, lags or time
		931
		932	noise = None #Noise Potency
		933
		934	SNR = None #Signal to Noise Ratio
		935
		936	pairsList = None #List of Pairs for Cross correlations or Cross spectrum
		937
		938	initUtcTime = None #Initial UTC time
		939
		940	paramInterval = None #Time interval to calculate Parameters in seconds
		941
		942	windsInterval = None #Time interval to calculate Winds in seconds
		943
		944	normFactor = None #Normalization Factor
		945
		946	winds = None #Wind estimations
		947
		948	def __init__(self):
		949	'''
		950	Constructor
		951	'''
		952	self.radarControllerHeaderObj = RadarControllerHeader()
		953
		954	self.systemHeaderObj = SystemHeader()
		955
		956	self.type = "Parameters"
		957
		958	def getTimeRange1(self):
		959
		960	datatime = []
		961
		962	datatime.append(self.initUtcTime)
		963	datatime.append(self.initUtcTime + self.windsInterval - 1)
		964
		965	datatime = numpy.array(datatime)
		966
		967	return datatime

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

@@ -572,6 +572,23 class Axes:
572	ylabel=ylabel,	572	ylabel=ylabel,
573	title=title,	573	title=title,
574	colormap=colormap)	574	colormap=colormap)
		575
		576	def polar(self, x, y,
		577	title='', xlabel='',ylabel='',**kwargs):
		578
		579	if self.__firsttime:
		580	self.plot = self.__driver.createPolar(self.ax, x, y, title = title, xlabel = xlabel, ylabel = ylabel)
		581	self.__firsttime = False
		582	self.x_buffer = x
		583	self.y_buffer = y
		584	return
		585
		586	self.x_buffer = numpy.hstack((self.x_buffer,x))
		587	self.y_buffer = numpy.hstack((self.y_buffer,y))
		588	self.__driver.polar(self.plot, self.x_buffer, self.y_buffer, xlabel=xlabel,
		589	ylabel=ylabel,
		590	title=title)
		591
575	def __fillGaps(self, x_buffer, y_buffer, z_buffer):	592	def __fillGaps(self, x_buffer, y_buffer, z_buffer):
576		593
577	deltas = x_buffer[1:] - x_buffer[0:-1]	594	deltas = x_buffer[1:] - x_buffer[0:-1]

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

                     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)

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

@@ -52,14 +52,15 def setTitle(fig, title):
52		52
53	fig.suptitle(title)	53	fig.suptitle(title)
54		54
55	def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan):	55	def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan, polar=False):
56		56
57	matplotlib.pyplot.ioff()	57	matplotlib.pyplot.ioff()
58	matplotlib.pyplot.figure(fig.number)	58	matplotlib.pyplot.figure(fig.number)
59	axes = matplotlib.pyplot.subplot2grid((nrow, ncol),	59	axes = matplotlib.pyplot.subplot2grid((nrow, ncol),
60	(xpos, ypos),	60	(xpos, ypos),
61	colspan=colspan,	61	colspan=colspan,
62	rowspan=rowspan)	62	rowspan=rowspan,
		63	polar=polar)
63		64
64	matplotlib.pyplot.ion()	65	matplotlib.pyplot.ion()
65	return axes	66	return axes
@@ -374,6 +375,49 def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''):
374	for i in range(len(ax.lines)):	375	for i in range(len(ax.lines)):
375	line = ax.lines[i]	376	line = ax.lines[i]
376	line.set_data(x,y[i,:])	377	line.set_data(x,y[i,:])
		378
		379	def createPolar(ax, x, y,
		380	xlabel='', ylabel='', title='', ticksize = 9,
		381	colormap='jet',cblabel='', cbsize="5%",
		382	XAxisAsTime=False):
		383
		384	matplotlib.pyplot.ioff()
		385
		386	ax.plot(x,y,'bo', markersize=5)
		387	# ax.set_rmax(90)
		388	ax.set_ylim(0,90)
		389	ax.set_yticks(numpy.arange(0,90,20))
		390	ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11')
		391	# ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical')
		392	printLabels(ax, xlabel, '', title)
		393	iplot = ax.lines[-1]
		394
		395	if '0.' in matplotlib.__version__[0:2]:
		396	print "The matplotlib version has to be updated to 1.1 or newer"
		397	return iplot
		398
		399	if '1.0.' in matplotlib.__version__[0:4]:
		400	print "The matplotlib version has to be updated to 1.1 or newer"
		401	return iplot
		402
		403	# if grid != None:
		404	# ax.grid(b=True, which='major', axis=grid)
		405
		406	matplotlib.pyplot.tight_layout()
		407
		408	matplotlib.pyplot.ion()
		409
		410
		411	return iplot
		412
		413	def polar(iplot, x, y, xlabel='', ylabel='', title=''):
		414
		415	ax = iplot.get_axes()
		416
		417	# ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11')
		418	printLabels(ax, xlabel, '', title)
		419
		420	set_linedata(ax, x, y, idline=0)
377		421
378	def draw(fig):	422	def draw(fig):
379		423

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

             import numpy
+            import math
             from jroproc_base import ProcessingUnit, Operation
             from model.data.jrodata import Spectra

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

             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.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"

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