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

		@@ -744,3 +744,224 class Fits:
744	744	datatime = property(getDatatime, "I'm the 'datatime' property")
745	745	ltctime = property(getltctime, "I'm the 'ltctime' property")
746	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
	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

                                              ylabel=ylabel,
                                              title=title,
                                              colormap=colormap)
+                 def polar(self, x, y,
+                            title='', xlabel='',ylabel='',**kwargs):
+                     if self.__firsttime:
+                         self.plot = self.__driver.createPolar(self.ax, x, y, title = title, xlabel = xlabel, ylabel = ylabel)
+                         self.__firsttime = False
+                         self.x_buffer = x
+                         self.y_buffer = y
+                         return
+                     self.x_buffer = numpy.hstack((self.x_buffer,x))
+                     self.y_buffer = numpy.hstack((self.y_buffer,y))
+                     self.__driver.polar(self.plot, self.x_buffer, self.y_buffer, xlabel=xlabel,
+                                                         ylabel=ylabel,
+                                                         title=title)
                  def __fillGaps(self, x_buffer, y_buffer, z_buffer):
                      deltas = x_buffer[1:] - x_buffer[0:-1]

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

                  fig.suptitle(title)
-             def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan):
+             def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan, polar=False):
                  matplotlib.pyplot.ioff()
                  matplotlib.pyplot.figure(fig.number)
                  axes = matplotlib.pyplot.subplot2grid((nrow, ncol),
                                                      (xpos, ypos),
                                                      colspan=colspan,
-                                                     rowspan=rowspan)
+                                                     rowspan=rowspan,
+                                                     polar=polar)
                  matplotlib.pyplot.ion()
                  return axes
                      line = ax.lines[i]
                      line.set_data(x,y[i,:])
+             def createPolar(ax, x, y,
+                             xlabel='', ylabel='', title='', ticksize = 9,
+                              colormap='jet',cblabel='', cbsize="5%",
+                              XAxisAsTime=False):
+                 matplotlib.pyplot.ioff()
+                 ax.plot(x,y,'bo', markersize=5)
+             #     ax.set_rmax(90)
+                 ax.set_ylim(0,90)
+                 ax.set_yticks(numpy.arange(0,90,20))
+                 ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11')
+             #     ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical')
+                 printLabels(ax, xlabel, '', title)
+                 iplot = ax.lines[-1]
+                 if '0.' in matplotlib.__version__[0:2]:
+                     print "The matplotlib version has to be updated to 1.1 or newer"
+                     return iplot
+                 if '1.0.' in matplotlib.__version__[0:4]:
+                     print "The matplotlib version has to be updated to 1.1 or newer"
+                     return iplot
+             #     if grid != None:
+             #         ax.grid(b=True, which='major', axis=grid)
+                 matplotlib.pyplot.tight_layout()
+                 matplotlib.pyplot.ion()
+                 return iplot
+             def polar(iplot, x, y, xlabel='', ylabel='', title=''):
+                 ax = iplot.get_axes()
+             #     ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11')
+                 printLabels(ax, xlabel, '', title)
+                 set_linedata(ax, x, y, idline=0)
              def draw(fig):
                  if type(fig) == 'int':

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 +2 0

              from jroproc_spectra import *
              from jroproc_heispectra import *
              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='CohInt', optype='other')
-             opObj11.addParameter(name='n', value='2', 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')
+             # 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='xmax', value='24', format='float')
+             opObj11.addParameter(name='xmin', value='20', 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