schain Commit - r1771:ed72d0635cd3 · Jicamarca Repository

UPDATE Repository desarroll drifts_schain, antigua version de WindProfiler en comparar 2.py

avaldez -

r1771:ed72d0635cd3

parent child

Context file:

r1771:ed72d0635cd3

schainpy/scripts/comparar2.py created 644 +663 0

This diff has been collapsed as it changes many lines, (663 lines changed) Show them Hide them
	@@ -0,0 +1,663
		1	#primer windprofiler
		2	class WindProfiler(Operation):
		3
		4	__isConfig = False
		5
		6	__initime = None
		7	__lastdatatime = None
		8	__integrationtime = None
		9
		10	__buffer = None
		11
		12	__dataReady = False
		13
		14	__firstdata = None
		15
		16	n = None
		17
		18	def __init__(self):
		19	Operation.__init__(self)
		20
		21	def __calculateCosDir(self, elev, azim):
		22	zen = (90 - elev)*numpy.pi/180
		23	azim = azim*numpy.pi/180
		24	cosDirX = numpy.sqrt((1-numpy.cos(zen)2)/((1+numpy.tan(azim)2)))
		25	cosDirY = numpy.sqrt(1-numpy.cos(zen)2-cosDirX2)
		26
		27	signX = numpy.sign(numpy.cos(azim))
		28	signY = numpy.sign(numpy.sin(azim))
		29
		30	cosDirX = numpy.copysign(cosDirX, signX)
		31	cosDirY = numpy.copysign(cosDirY, signY)
		32	return cosDirX, cosDirY
		33
		34	def __calculateAngles(self, theta_x, theta_y, azimuth):
		35
		36	dir_cosw = numpy.sqrt(1-theta_x2-theta_y2)
		37	zenith_arr = numpy.arccos(dir_cosw)
		38	azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
		39
		40	dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
		41	dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
		42
		43	return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
		44
		45	def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
		46
		47	if horOnly:
		48	A = numpy.c_[dir_cosu,dir_cosv]
		49	else:
		50	A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
		51	A = numpy.asmatrix(A)
		52	A1 = numpy.linalg.inv(A.transpose()A)A.transpose()
		53
		54	return A1
		55
		56	def __correctValues(self, heiRang, phi, velRadial, SNR):
		57	listPhi = phi.tolist()
		58	maxid = listPhi.index(max(listPhi))
		59	minid = listPhi.index(min(listPhi))
		60
		61	rango = list(range(len(phi)))
		62
		63	heiRang1 = heiRang*math.cos(phi[maxid])
		64	heiRangAux = heiRang*math.cos(phi[minid])
		65	indOut = (heiRang1 < heiRangAux[0]).nonzero()
		66	heiRang1 = numpy.delete(heiRang1,indOut)
		67
		68	velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
		69	SNR1 = numpy.zeros([len(phi),len(heiRang1)])
		70
		71	for i in rango:
		72	x = heiRang*math.cos(phi[i])
		73	y1 = velRadial[i,:]
		74	f1 = interpolate.interp1d(x,y1,kind = 'cubic')
		75
		76	x1 = heiRang1
		77	y11 = f1(x1)
		78
		79	y2 = SNR[i,:]
		80	f2 = interpolate.interp1d(x,y2,kind = 'cubic')
		81	y21 = f2(x1)
		82
		83	velRadial1[i,:] = y11
		84	SNR1[i,:] = y21
		85
		86	return heiRang1, velRadial1, SNR1
		87
		88	def __calculateVelUVW(self, A, velRadial):
		89
		90	#Operacion Matricial
		91	velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
		92	velUVW[:,:] = numpy.dot(A,velRadial)
		93
		94
		95	return velUVW
		96
		97	def techniqueDBS(self, kwargs):
		98	"""
		99	Function that implements Doppler Beam Swinging (DBS) technique.
		100
		101	Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
		102	Direction correction (if necessary), Ranges and SNR
		103
		104	Output: Winds estimation (Zonal, Meridional and Vertical)
		105
		106	Parameters affected: Winds, height range, SNR
		107	"""
		108	velRadial0 = kwargs['velRadial']
		109	heiRang = kwargs['heightList']
		110	SNR0 = kwargs['SNR']
		111
		112	if 'dirCosx' in kwargs and 'dirCosy' in kwargs:
		113	theta_x = numpy.array(kwargs['dirCosx'])
		114	theta_y = numpy.array(kwargs['dirCosy'])
		115	else:
		116	elev = numpy.array(kwargs['elevation'])
		117	azim = numpy.array(kwargs['azimuth'])
		118	theta_x, theta_y = self.__calculateCosDir(elev, azim)
		119	azimuth = kwargs['correctAzimuth']
		120	if 'horizontalOnly' in kwargs:
		121	horizontalOnly = kwargs['horizontalOnly']
		122	else: horizontalOnly = False
		123	if 'correctFactor' in kwargs:
		124	correctFactor = kwargs['correctFactor']
		125	else: correctFactor = 1
		126	if 'channelList' in kwargs:
		127	channelList = kwargs['channelList']
		128	if len(channelList) == 2:
		129	horizontalOnly = True
		130	arrayChannel = numpy.array(channelList)
		131	param = param[arrayChannel,:,:]
		132	theta_x = theta_x[arrayChannel]
		133	theta_y = theta_y[arrayChannel]
		134
		135	azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
		136	heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
		137	A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
		138
		139	#Calculo de Componentes de la velocidad con DBS
		140	winds = self.__calculateVelUVW(A,velRadial1)
		141
		142	return winds, heiRang1, SNR1
		143
		144	def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
		145
		146	nPairs = len(pairs_ccf)
		147	posx = numpy.asarray(posx)
		148	posy = numpy.asarray(posy)
		149
		150	#Rotacion Inversa para alinear con el azimuth
		151	if azimuth!= None:
		152	azimuth = azimuth*math.pi/180
		153	posx1 = posxmath.cos(azimuth) + posymath.sin(azimuth)
		154	posy1 = -posxmath.sin(azimuth) + posymath.cos(azimuth)
		155	else:
		156	posx1 = posx
		157	posy1 = posy
		158
		159	#Calculo de Distancias
		160	distx = numpy.zeros(nPairs)
		161	disty = numpy.zeros(nPairs)
		162	dist = numpy.zeros(nPairs)
		163	ang = numpy.zeros(nPairs)
		164
		165	for i in range(nPairs):
		166	distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
		167	disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
		168	dist[i] = numpy.sqrt(distx[i]2 + disty[i]2)
		169	ang[i] = numpy.arctan2(disty[i],distx[i])
		170
		171	return distx, disty, dist, ang
		172	#Calculo de Matrices
		173
		174
		175	def __calculateVelVer(self, phase, lagTRange, _lambda):
		176
		177	Ts = lagTRange[1] - lagTRange[0]
		178	velW = -_lambdaphase/(4math.pi*Ts)
		179
		180	return velW
		181
		182	def __calculateVelHorDir(self, dist, tau1, tau2, ang):
		183	nPairs = tau1.shape[0]
		184	nHeights = tau1.shape[1]
		185	vel = numpy.zeros((nPairs,3,nHeights))
		186	dist1 = numpy.reshape(dist, (dist.size,1))
		187
		188	angCos = numpy.cos(ang)
		189	angSin = numpy.sin(ang)
		190
		191	vel0 = dist1tau1/(2tau2**2)
		192	vel[:,0,:] = (vel0*angCos).sum(axis = 1)
		193	vel[:,1,:] = (vel0*angSin).sum(axis = 1)
		194
		195	ind = numpy.where(numpy.isinf(vel))
		196	vel[ind] = numpy.nan
		197
		198	return vel
		199
		200	def techniqueSA(self, kwargs):
		201
		202	"""
		203	Function that implements Spaced Antenna (SA) technique.
		204
		205	Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
		206	Direction correction (if necessary), Ranges and SNR
		207
		208	Output: Winds estimation (Zonal, Meridional and Vertical)
		209
		210	Parameters affected: Winds
		211	"""
		212	position_x = kwargs['positionX']
		213	position_y = kwargs['positionY']
		214	azimuth = kwargs['azimuth']
		215
		216	if 'correctFactor' in kwargs:
		217	correctFactor = kwargs['correctFactor']
		218	else:
		219	correctFactor = 1
		220
		221	groupList = kwargs['groupList']
		222	pairs_ccf = groupList[1]
		223	tau = kwargs['tau']
		224	_lambda = kwargs['_lambda']
		225
		226	#Cross Correlation pairs obtained
		227
		228	indtau = tau.shape[0]/2
		229	tau1 = tau[:indtau,:]
		230	tau2 = tau[indtau:-1,:]
		231	phase1 = tau[-1,:]
		232
		233	#---------------------------------------------------------------------
		234	#Metodo Directo
		235	distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
		236	winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
		237	winds = stats.nanmean(winds, axis=0)
		238	#---------------------------------------------------------------------
		239	#Metodo General
		240
		241	#---------------------------------------------------------------------
		242	winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
		243	winds = correctFactor*winds
		244	return winds
		245
		246	def __checkTime(self, currentTime, paramInterval, outputInterval):
		247
		248	dataTime = currentTime + paramInterval
		249	deltaTime = dataTime - self.__initime
		250
		251	if deltaTime >= outputInterval or deltaTime < 0:
		252	self.__dataReady = True
		253	return
		254
		255	def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
		256	'''
		257	Function that implements winds estimation technique with detected meteors.
		258
		259	Input: Detected meteors, Minimum meteor quantity to wind estimation
		260
		261	Output: Winds estimation (Zonal and Meridional)
		262
		263	Parameters affected: Winds
		264	'''
		265	#Settings
		266	nInt = (heightMax - heightMin)/2
		267	nInt = int(nInt)
		268	winds = numpy.zeros((2,nInt))*numpy.nan
		269
		270	#Filter errors
		271	error = numpy.where(arrayMeteor[:,-1] == 0)[0]
		272	finalMeteor = arrayMeteor[error,:]
		273
		274	#Meteor Histogram
		275	finalHeights = finalMeteor[:,2]
		276	hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
		277	nMeteorsPerI = hist[0]
		278	heightPerI = hist[1]
		279
		280	#Sort of meteors
		281	indSort = finalHeights.argsort()
		282	finalMeteor2 = finalMeteor[indSort,:]
		283
		284	# Calculating winds
		285	ind1 = 0
		286	ind2 = 0
		287
		288	for i in range(nInt):
		289	nMet = nMeteorsPerI[i]
		290	ind1 = ind2
		291	ind2 = ind1 + nMet
		292
		293	meteorAux = finalMeteor2[ind1:ind2,:]
		294
		295	if meteorAux.shape[0] >= meteorThresh:
		296	vel = meteorAux[:, 6]
		297	zen = meteorAux[:, 4]*numpy.pi/180
		298	azim = meteorAux[:, 3]*numpy.pi/180
		299
		300	n = numpy.cos(zen)
		301	# m = (1 - n2)/(1 - numpy.tan(azim)2)
		302	# l = m*numpy.tan(azim)
		303	l = numpy.sin(zen)*numpy.sin(azim)
		304	m = numpy.sin(zen)*numpy.cos(azim)
		305
		306	A = numpy.vstack((l, m)).transpose()
		307	A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
		308	windsAux = numpy.dot(A1, vel)
		309
		310	winds[0,i] = windsAux[0]
		311	winds[1,i] = windsAux[1]
		312
		313	return winds, heightPerI[:-1]
		314
		315	def techniqueNSM_SA(self, **kwargs):
		316	metArray = kwargs['metArray']
		317	heightList = kwargs['heightList']
		318	timeList = kwargs['timeList']
		319
		320	rx_location = kwargs['rx_location']
		321	groupList = kwargs['groupList']
		322	azimuth = kwargs['azimuth']
		323	dfactor = kwargs['dfactor']
		324	k = kwargs['k']
		325
		326	azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
		327	d = dist*dfactor
		328	#Phase calculation
		329	metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
		330
		331	metArray1[:,-2] = metArray1[:,-2]metArray1[:,2]1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
		332
		333	velEst = numpy.zeros((heightList.size,2))*numpy.nan
		334	azimuth1 = azimuth1*numpy.pi/180
		335
		336	for i in range(heightList.size):
		337	h = heightList[i]
		338	indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
		339	metHeight = metArray1[indH,:]
		340	if metHeight.shape[0] >= 2:
		341	velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities
		342	iazim = metHeight[:,1].astype(int)
		343	azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths
		344	A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
		345	A = numpy.asmatrix(A)
		346	A1 = numpy.linalg.pinv(A.transpose()A)A.transpose()
		347	velHor = numpy.dot(A1,velAux)
		348
		349	velEst[i,:] = numpy.squeeze(velHor)
		350	return velEst
		351
		352	def __getPhaseSlope(self, metArray, heightList, timeList):
		353	meteorList = []
		354	#utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
		355	#Putting back together the meteor matrix
		356	utctime = metArray[:,0]
		357	uniqueTime = numpy.unique(utctime)
		358
		359	phaseDerThresh = 0.5
		360	ippSeconds = timeList[1] - timeList[0]
		361	sec = numpy.where(timeList>1)[0][0]
		362	nPairs = metArray.shape[1] - 6
		363	nHeights = len(heightList)
		364
		365	for t in uniqueTime:
		366	metArray1 = metArray[utctime==t,:]
		367	# phaseDerThresh = numpy.pi/4 #reducir Phase thresh
		368	tmet = metArray1[:,1].astype(int)
		369	hmet = metArray1[:,2].astype(int)
		370
		371	metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
		372	metPhase[:,:] = numpy.nan
		373	metPhase[:,hmet,tmet] = metArray1[:,6:].T
		374
		375	#Delete short trails
		376	metBool = ~numpy.isnan(metPhase[0,:,:])
		377	heightVect = numpy.sum(metBool, axis = 1)
		378	metBool[heightVect<sec,:] = False
		379	metPhase[:,heightVect<sec,:] = numpy.nan
		380
		381	#Derivative
		382	metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
		383	phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
		384	metPhase[phDerAux] = numpy.nan
		385
		386	#--------------------------METEOR DETECTION -----------------------------------------
		387	indMet = numpy.where(numpy.any(metBool,axis=1))[0]
		388
		389	for p in numpy.arange(nPairs):
		390	phase = metPhase[p,:,:]
		391	phDer = metDer[p,:,:]
		392
		393	for h in indMet:
		394	height = heightList[h]
		395	phase1 = phase[h,:] #82
		396	phDer1 = phDer[h,:]
		397
		398	phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap
		399
		400	indValid = numpy.where(~numpy.isnan(phase1))[0]
		401	initMet = indValid[0]
		402	endMet = 0
		403
		404	for i in range(len(indValid)-1):
		405
		406	#Time difference
		407	inow = indValid[i]
		408	inext = indValid[i+1]
		409	idiff = inext - inow
		410	#Phase difference
		411	phDiff = numpy.abs(phase1[inext] - phase1[inow])
		412
		413	if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor
		414	sizeTrail = inow - initMet + 1
		415	if sizeTrail>3*sec: #Too short meteors
		416	x = numpy.arange(initMet,inow+1)*ippSeconds
		417	y = phase1[initMet:inow+1]
		418	ynnan = ~numpy.isnan(y)
		419	x = x[ynnan]
		420	y = y[ynnan]
		421	slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
		422	ylin = x*slope + intercept
		423	rsq = r_value**2
		424	if rsq > 0.5:
		425	vel = slope#height1000/(k*d)
		426	estAux = numpy.array([utctime,p,height, vel, rsq])
		427	meteorList.append(estAux)
		428	initMet = inext
		429	metArray2 = numpy.array(meteorList)
		430
		431	return metArray2
		432
		433	def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
		434
		435	azimuth1 = numpy.zeros(len(pairslist))
		436	dist = numpy.zeros(len(pairslist))
		437
		438	for i in range(len(rx_location)):
		439	ch0 = pairslist[i][0]
		440	ch1 = pairslist[i][1]
		441
		442	diffX = rx_location[ch0][0] - rx_location[ch1][0]
		443	diffY = rx_location[ch0][1] - rx_location[ch1][1]
		444	azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
		445	dist[i] = numpy.sqrt(diffX2 + diffY2)
		446
		447	azimuth1 -= azimuth0
		448	return azimuth1, dist
		449
		450	def techniqueNSM_DBS(self, **kwargs):
		451	metArray = kwargs['metArray']
		452	heightList = kwargs['heightList']
		453	timeList = kwargs['timeList']
		454	azimuth = kwargs['azimuth']
		455	theta_x = numpy.array(kwargs['theta_x'])
		456	theta_y = numpy.array(kwargs['theta_y'])
		457
		458	utctime = metArray[:,0]
		459	cmet = metArray[:,1].astype(int)
		460	hmet = metArray[:,3].astype(int)
		461	SNRmet = metArray[:,4]
		462	vmet = metArray[:,5]
		463	spcmet = metArray[:,6]
		464
		465	nChan = numpy.max(cmet) + 1
		466	nHeights = len(heightList)
		467
		468	azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
		469	hmet = heightList[hmet]
		470	h1met = hmet*numpy.cos(zenith_arr[cmet]) #Corrected heights
		471
		472	velEst = numpy.zeros((heightList.size,2))*numpy.nan
		473
		474	for i in range(nHeights - 1):
		475	hmin = heightList[i]
		476	hmax = heightList[i + 1]
		477
		478	thisH = (h1met>=hmin) & (h1met<hmax) & (cmet!=2) & (SNRmet>8) & (vmet<50) & (spcmet<10)
		479	indthisH = numpy.where(thisH)
		480
		481	if numpy.size(indthisH) > 3:
		482
		483	vel_aux = vmet[thisH]
		484	chan_aux = cmet[thisH]
		485	cosu_aux = dir_cosu[chan_aux]
		486	cosv_aux = dir_cosv[chan_aux]
		487	cosw_aux = dir_cosw[chan_aux]
		488
		489	nch = numpy.size(numpy.unique(chan_aux))
		490	if nch > 1:
		491	A = self.__calculateMatA(cosu_aux, cosv_aux, cosw_aux, True)
		492	velEst[i,:] = numpy.dot(A,vel_aux)
		493
		494	return velEst
		495
		496	def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
		497
		498	param = dataOut.data_param
		499	#if dataOut.abscissaList != None:
		500	if numpy.any(dataOut.abscissaList):
		501	absc = dataOut.abscissaList[:-1]
		502	# noise = dataOut.noise
		503	heightList = dataOut.heightList
		504	SNR = dataOut.data_snr
		505
		506	if technique == 'DBS':
		507
		508	kwargs['velRadial'] = param[:,1,:] #Radial velocity
		509	kwargs['heightList'] = heightList
		510	kwargs['SNR'] = SNR
		511
		512	dataOut.data_output, dataOut.heightList, dataOut.data_snr = self.techniqueDBS(kwargs) #DBS Function
		513	dataOut.utctimeInit = dataOut.utctime
		514	dataOut.outputInterval = dataOut.paramInterval
		515
		516	elif technique == 'SA':
		517
		518	#Parameters
		519	# position_x = kwargs['positionX']
		520	# position_y = kwargs['positionY']
		521	# azimuth = kwargs['azimuth']
		522	#
		523	# if kwargs.has_key('crosspairsList'):
		524	# pairs = kwargs['crosspairsList']
		525	# else:
		526	# pairs = None
		527	#
		528	# if kwargs.has_key('correctFactor'):
		529	# correctFactor = kwargs['correctFactor']
		530	# else:
		531	# correctFactor = 1
		532
		533	# tau = dataOut.data_param
		534	# _lambda = dataOut.C/dataOut.frequency
		535	# pairsList = dataOut.groupList
		536	# nChannels = dataOut.nChannels
		537
		538	kwargs['groupList'] = dataOut.groupList
		539	kwargs['tau'] = dataOut.data_param
		540	kwargs['_lambda'] = dataOut.C/dataOut.frequency
		541	# dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
		542	dataOut.data_output = self.techniqueSA(kwargs)
		543	dataOut.utctimeInit = dataOut.utctime
		544	dataOut.outputInterval = dataOut.timeInterval
		545
		546	elif technique == 'Meteors':
		547	dataOut.flagNoData = True
		548	self.__dataReady = False
		549
		550	if 'nHours' in kwargs:
		551	nHours = kwargs['nHours']
		552	else:
		553	nHours = 1
		554
		555	if 'meteorsPerBin' in kwargs:
		556	meteorThresh = kwargs['meteorsPerBin']
		557	else:
		558	meteorThresh = 6
		559
		560	if 'hmin' in kwargs:
		561	hmin = kwargs['hmin']
		562	else: hmin = 70
		563	if 'hmax' in kwargs:
		564	hmax = kwargs['hmax']
		565	else: hmax = 110
		566
		567	dataOut.outputInterval = nHours*3600
		568
		569	if self.__isConfig == False:
		570	# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
		571	#Get Initial LTC time
		572	self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
		573	self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
		574
		575	self.__isConfig = True
		576
		577	if self.__buffer is None:
		578	self.__buffer = dataOut.data_param
		579	self.__firstdata = copy.copy(dataOut)
		580
		581	else:
		582	self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
		583
		584	self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
		585
		586	if self.__dataReady:
		587	dataOut.utctimeInit = self.__initime
		588
		589	self.__initime += dataOut.outputInterval #to erase time offset
		590
		591	dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
		592	dataOut.flagNoData = False
		593	self.__buffer = None
		594
		595	elif technique == 'Meteors1':
		596	dataOut.flagNoData = True
		597	self.__dataReady = False
		598
		599	if 'nMins' in kwargs:
		600	nMins = kwargs['nMins']
		601	else: nMins = 20
		602	if 'rx_location' in kwargs:
		603	rx_location = kwargs['rx_location']
		604	else: rx_location = [(0,1),(1,1),(1,0)]
		605	if 'azimuth' in kwargs:
		606	azimuth = kwargs['azimuth']
		607	else: azimuth = 51.06
		608	if 'dfactor' in kwargs:
		609	dfactor = kwargs['dfactor']
		610	if 'mode' in kwargs:
		611	mode = kwargs['mode']
		612	if 'theta_x' in kwargs:
		613	theta_x = kwargs['theta_x']
		614	if 'theta_y' in kwargs:
		615	theta_y = kwargs['theta_y']
		616	else: mode = 'SA'
		617
		618	#Borrar luego esto
		619	if dataOut.groupList is None:
		620	dataOut.groupList = [(0,1),(0,2),(1,2)]
		621	groupList = dataOut.groupList
		622	C = 3e8
		623	freq = 50e6
		624	lamb = C/freq
		625	k = 2*numpy.pi/lamb
		626
		627	timeList = dataOut.abscissaList
		628	heightList = dataOut.heightList
		629
		630	if self.__isConfig == False:
		631	dataOut.outputInterval = nMins*60
		632	# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
		633	#Get Initial LTC time
		634	initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
		635	minuteAux = initime.minute
		636	minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
		637	self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
		638
		639	self.__isConfig = True
		640
		641	if self.__buffer is None:
		642	self.__buffer = dataOut.data_param
		643	self.__firstdata = copy.copy(dataOut)
		644
		645	else:
		646	self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
		647
		648	self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
		649
		650	if self.__dataReady:
		651	dataOut.utctimeInit = self.__initime
		652	self.__initime += dataOut.outputInterval #to erase time offset
		653
		654	metArray = self.__buffer
		655	if mode == 'SA':
		656	dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
		657	elif mode == 'DBS':
		658	dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList, azimuth=azimuth, theta_x=theta_x, theta_y=theta_y)
		659	dataOut.data_output = dataOut.data_output.T
		660	dataOut.flagNoData = False
		661	self.__buffer = None
		662
		663	return

schainpy/model/proc/jroproc_parameters.py +375 -807

This diff has been collapsed as it changes many lines, (1182 lines changed) Show them Hide them
	@@ -8,6 +8,7 import re
	8	import datetime	8	import datetime
	9	import copy	9	import copy
	10	import sys	10	import sys
			11	import os
	11	import importlib	12	import importlib
	12	import itertools	13	import itertools
	13	from multiprocessing import Pool, TimeoutError	14	from multiprocessing import Pool, TimeoutError
	@@ -18,12 +19,13 from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papamete
	18	from .jroproc_base import ProcessingUnit, Operation, MPDecorator	19	from .jroproc_base import ProcessingUnit, Operation, MPDecorator
	19	from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon	20	from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
	20	from schainpy.model.data.jrodata import Spectra	21	from schainpy.model.data.jrodata import Spectra
	21	#from scipy import asarray as ar,exp	22	from scipy import asarray as ar,exp
	22	from scipy.optimize import fmin, curve_fit	23	from scipy.optimize import fmin, curve_fit
	23	from schainpy.utils import log	24	from schainpy.utils import log
	24	import warnings	25	import warnings
	25	from numpy import NaN	26	from numpy import NaN
	26	from scipy.optimize.optimize import OptimizeWarning	27	from scipy.optimize.optimize import OptimizeWarning
			28	from schainpy.model.io.jroIO_madrigal import MADWriter
	27	warnings.filterwarnings('ignore')	29	warnings.filterwarnings('ignore')
	28		30
	29	import os	31	import os
	@@ -31,6 +33,12 import csv
	31	from scipy import signal	33	from scipy import signal
	32	import matplotlib.pyplot as plt	34	import matplotlib.pyplot as plt
	33		35
			36	import json
			37	import mysql.connector
			38	from scipy.signal import savgol_filter
			39	from scipy.signal import argrelextrema
			40	from scipy.signal import find_peaks
			41
	34	SPEED_OF_LIGHT = 299792458	42	SPEED_OF_LIGHT = 299792458
	35		43
	36	'''solving pickling issue'''	44	'''solving pickling issue'''
	@@ -51,7 +59,7 def _unpickle_method(func_name, obj, cls):
	51	break	59	break
	52	return func.__get__(obj, cls)	60	return func.__get__(obj, cls)
	53		61
	54		62	# @MPDecorator
	55	class ParametersProc(ProcessingUnit):	63	class ParametersProc(ProcessingUnit):
	56		64
	57	METHODS = {}	65	METHODS = {}
	@@ -171,6 +179,8 class ParametersProc(ProcessingUnit):
	171	self.dataOut.spc_noise = self.dataIn.getNoise()	179	self.dataOut.spc_noise = self.dataIn.getNoise()
	172	self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))	180	self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
	173	# self.dataOut.normFactor = self.dataIn.normFactor	181	# self.dataOut.normFactor = self.dataIn.normFactor
			182	self.dataOut.pairsList = self.dataIn.pairsList
			183	self.dataOut.groupList = self.dataIn.pairsList
	174	self.dataOut.flagNoData = False	184	self.dataOut.flagNoData = False
	175		185
	176	if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels	186	if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
	@@ -224,6 +234,7 class ParametersProc(ProcessingUnit):
	224	self.__updateObjFromInput()	234	self.__updateObjFromInput()
	225	self.dataOut.utctimeInit = self.dataIn.utctime	235	self.dataOut.utctimeInit = self.dataIn.utctime
	226	self.dataOut.paramInterval = self.dataIn.timeInterval	236	self.dataOut.paramInterval = self.dataIn.timeInterval
			237
	227	return	238	return
	228		239
	229		240
	@@ -285,6 +296,7 class RemoveWideGC(Operation):
	285	continue	296	continue
	286	junk3 = numpy.squeeze(numpy.diff(j1index))	297	junk3 = numpy.squeeze(numpy.diff(j1index))
	287	junk4 = numpy.squeeze(numpy.diff(j2index))	298	junk4 = numpy.squeeze(numpy.diff(j2index))
			299
	288	valleyindex = j2index[numpy.where(junk4>1)]	300	valleyindex = j2index[numpy.where(junk4>1)]
	289	peakindex = j1index[numpy.where(junk3>1)]	301	peakindex = j1index[numpy.where(junk3>1)]
	290		302
	@@ -294,6 +306,7 class RemoveWideGC(Operation):
	294	if numpy.size(isvalid) >1 :	306	if numpy.size(isvalid) >1 :
	295	vindex = numpy.argmax(self.spc[ich,gc_values[peakindex[isvalid]],ir])	307	vindex = numpy.argmax(self.spc[ich,gc_values[peakindex[isvalid]],ir])
	296	isvalid = isvalid[vindex]	308	isvalid = isvalid[vindex]
			309
	297	# clutter peak	310	# clutter peak
	298	gcpeak = peakindex[isvalid]	311	gcpeak = peakindex[isvalid]
	299	vl = numpy.where(valleyindex < gcpeak)	312	vl = numpy.where(valleyindex < gcpeak)
	@@ -2655,7 +2668,7 class SpectralMoments(Operation):
	2655		2668
	2656	if (type1 is None): type1 = 0	2669	if (type1 is None): type1 = 0
	2657	if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1	2670	if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
	2658	if (snrth is None): snrth = -3 ~~#-20.0~~	2671	if (snrth is None): snrth = -3
	2659	if (dc is None): dc = 0	2672	if (dc is None): dc = 0
	2660	if (aliasing is None): aliasing = 0	2673	if (aliasing is None): aliasing = 0
	2661	if (oldfd is None): oldfd = 0	2674	if (oldfd is None): oldfd = 0
	@@ -3190,6 +3203,7 class SpectralFitting(Operation):
	3190	Operation.__init__(self)	3203	Operation.__init__(self)
	3191	self.i=0	3204	self.i=0
	3192	self.isConfig = False	3205	self.isConfig = False
			3206	self.aux = 1
	3193		3207
	3194	def setup(self,nChan,nProf,nHei,nBlocks):	3208	def setup(self,nChan,nProf,nHei,nBlocks):
	3195	self.__dataReady = False	3209	self.__dataReady = False
	@@ -3212,6 +3226,7 class SpectralFitting(Operation):
	3212	if (wwauto is None): wwauto = 0	3226	if (wwauto is None): wwauto = 0
	3213		3227
	3214	if (n0 < 1.e-20): n0 = 1.e-20	3228	if (n0 < 1.e-20): n0 = 1.e-20
			3229
	3215	freq = oldfreq	3230	freq = oldfreq
	3216	vec_power = numpy.zeros(oldspec.shape[1])	3231	vec_power = numpy.zeros(oldspec.shape[1])
	3217	vec_fd = numpy.zeros(oldspec.shape[1])	3232	vec_fd = numpy.zeros(oldspec.shape[1])
	@@ -3620,7 +3635,7 class SpectralFitting(Operation):
	3620	junk = tmp	3635	junk = tmp
	3621	tmp = junk*0.0	3636	tmp = junk*0.0
	3622		3637
	3623	tmp[list(xpos + (yposnum_~~hei~~))] = junk[list(xpos + (yposnum_~~hei~~))]	3638	tmp[list(xpos + (yposnum_prof))] = junk[list(xpos + (yposnum_prof))]
	3624	array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))	3639	array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
	3625		3640
	3626	return array	3641	return array
	@@ -3642,6 +3657,146 class SpectralFitting(Operation):
	3642	fmom = numpy.sum(dopplerytemp)/numpy.sum(ytemp)+(index-(npoints/2-1))numpy.abs(doppler[1]-doppler[0])	3657	fmom = numpy.sum(dopplerytemp)/numpy.sum(ytemp)+(index-(npoints/2-1))numpy.abs(doppler[1]-doppler[0])
	3643	smom = numpy.sum(dopplerdopplerytemp)/numpy.sum(ytemp)	3658	smom = numpy.sum(dopplerdopplerytemp)/numpy.sum(ytemp)
	3644	return [fmom,numpy.sqrt(smom)]	3659	return [fmom,numpy.sqrt(smom)]
			3660
			3661	def find_max_excluding_neighbors_and_borders(self,matrix):
			3662	n, m = matrix.shape # Dimensions of the matrix
			3663	if n < 3 or m < 3:
			3664	raise ValueError("Matrix too small to exclude borders and neighbors of NaN.")
			3665
			3666	# Create a mask for the borders
			3667	border_mask = numpy.zeros_like(matrix, dtype=bool)
			3668	border_mask[0, :] = True # Top border
			3669	border_mask[-1, :] = True # Bottom border
			3670	border_mask[:, 0] = True # Left border
			3671	border_mask[:, -1] = True # Right border
			3672
			3673	# Exclude elements within 3 indices from the border
			3674	three_index_border_mask = numpy.zeros_like(matrix, dtype=bool)
			3675	three_index_border_mask[0:3, :] = True # Top 3 rows
			3676	three_index_border_mask[-3:, :] = True # Bottom 3 rows
			3677	three_index_border_mask[:, 0:3] = True # Left 3 columns
			3678	three_index_border_mask[:, -3:] = True # Right 3 columns
			3679
			3680	# Create a mask for neighbors of NaNs
			3681	nan_mask = numpy.isnan(matrix)
			3682	nan_neighbors_mask = numpy.zeros_like(matrix, dtype=bool)
			3683	for i in range(1, n - 1):
			3684	for j in range(1, m - 1):
			3685	if nan_mask[i-1:i+2, j-1:j+2].any(): # Check 3x3 region around (i, j)
			3686	nan_neighbors_mask[i, j] = True
			3687
			3688	# Combine all masks (borders, 3 indices from border, and NaN neighbors)
			3689	combined_mask = border_mask \| three_index_border_mask \| nan_neighbors_mask
			3690
			3691	# Exclude these positions from the matrix
			3692	valid_matrix = numpy.where(combined_mask, numpy.nan, matrix)
			3693
			3694	# Find the maximum value index in the valid matrix
			3695	if numpy.isnan(valid_matrix).all():
			3696	return None # No valid maximum
			3697	max_index = numpy.unravel_index(numpy.nanargmax(valid_matrix), matrix.shape)
			3698	return max_index, valid_matrix
			3699
			3700	def construct_three_phase_path(self,maximas_2, max_index):
			3701	"""
			3702	Constructs a path through maximas_2 starting from the middle point (max_index)
			3703	and going downwards, passing through the starting point, and then going upwards.
			3704
			3705	Parameters:
			3706	- maximas_2 (numpy.ndarray): 2D array of points.
			3707	- max_index (tuple): Starting position as (row_index, starting_value).
			3708
			3709	Returns:
			3710	- list: A path of points in the three phases: down, through the start, then up.
			3711	"""
			3712	row_index, starting_value = max_index
			3713	path = []
			3714	# Initialize the path with the starting value
			3715	path_o = [starting_value]
			3716	path_up = []
			3717	path_dw = []
			3718	threshold = 4
			3719	flag_broken = 0
			3720	print(row_index,maximas_2.shape[0])
			3721	# Phase 1: Going downwards (rows below the starting point)
			3722	current_value = starting_value
			3723	for i in range(row_index + 1, maximas_2.shape[0]):
			3724
			3725	if flag_broken == 1:
			3726	path_up.append(numpy.nan)
			3727	continue
			3728
			3729	distances = numpy.abs(maximas_2[i] - current_value)
			3730	closest_index = numpy.argmin(distances)
			3731
			3732	if numpy.nanmin(distances) > threshold:
			3733	path_up.append(numpy.nan)
			3734	flag_broken = 1
			3735	continue
			3736
			3737	current_value = maximas_2[i, closest_index]
			3738	path_up.append(current_value)
			3739
			3740	flag_broken = 0
			3741
			3742	# Phase 2: Going upwards (rows above the starting point)
			3743	current_value = starting_value
			3744	for i in range(row_index - 1, -1, -1):
			3745	if flag_broken == 1:
			3746	path_dw.append(numpy.nan)
			3747	continue
			3748	distances = numpy.abs(maximas_2[i] - current_value)
			3749	closest_index = numpy.argmin(distances)
			3750
			3751	if numpy.nanmin(distances) > threshold:
			3752	path_dw.append(numpy.nan)
			3753	flag_broken = 1
			3754	continue
			3755
			3756	current_value = maximas_2[i, closest_index]
			3757	path_dw.append(current_value)
			3758	path.append(path_dw[::-1])
			3759	path.append(path_o)
			3760	path.append(path_up)
			3761	return numpy.concatenate(path)
			3762
			3763	def Vr_correction(self,dataset):
			3764
			3765	dataset = savgol_filter(dataset, window_length=15, polyorder=3, axis=1)
			3766	dataset = savgol_filter(dataset, window_length=15, polyorder=3, axis=1)
			3767	dataset = savgol_filter(dataset, window_length=7, polyorder=3, axis=0)
			3768	dataset[:12,:] = numpy.nan
			3769	dataset[45:, :] = numpy.nan
			3770	max_index, _ = self.find_max_excluding_neighbors_and_borders(dataset)
			3771	max_index = numpy.array(max_index)
			3772	print(max_index)
			3773	maximas = argrelextrema(dataset, numpy.greater, axis=1)
			3774	heighs = numpy.unique(maximas[0])
			3775	#grouped_arrays = {value: maximas[0][numpy.where(maximas[0] == value)] for value in numpy.unique(maximas[0])}
			3776	grouped_arrays = [maximas[0][numpy.where(maximas[0] == value)] for value in numpy.unique(maximas[0])]
			3777
			3778	maximas_2_list = []
			3779	num_maximas = 3
			3780	for i, h in enumerate(heighs):
			3781	n = len(grouped_arrays[i])
			3782	idx = numpy.where(maximas[0] == h)[0]
			3783	maxs_ = numpy.argsort(dataset[h, maximas[1][idx]])[::-1][:num_maximas]
			3784	maxs = maximas[1][idx][maxs_]
			3785	heigh = numpy.ones(len(maxs)) * h
			3786	arr = numpy.array([heigh, maxs])
			3787	if n < num_maximas:
			3788	arr = numpy.hstack((arr, numpy.full((arr.shape[0], num_maximas - arr.shape[1]), numpy.nan)))
			3789	maximas_2_list.append(arr)
			3790
			3791	maximas_2 = numpy.array(maximas_2_list)
			3792	_maximas_2_ = maximas_2[:, 1]
			3793	#
			3794	max_index[0] = numpy.where(heighs == max_index[0])[0][0]
			3795	#
			3796	path = self.construct_three_phase_path(_maximas_2_, max_index)
			3797	return int(numpy.nanmedian(path))
			3798
			3799
	3645		3800
	3646	def windowing_single_old(self,spc,x,A,B,C,D,nFFTPoints):	3801	def windowing_single_old(self,spc,x,A,B,C,D,nFFTPoints):
	3647	'''	3802	'''
	@@ -4010,7 +4165,6 class SpectralFitting(Operation):
	4010	power = numpy.sum(spectra, axis=1)	4165	power = numpy.sum(spectra, axis=1)
	4011	nPairs = len(crosspairs)	4166	nPairs = len(crosspairs)
	4012	absc = dataOut.abscissaList[:-1]	4167	absc = dataOut.abscissaList[:-1]
	4013	print('para escribir h5 ',dataOut.paramInterval)
	4014	if not self.isConfig:	4168	if not self.isConfig:
	4015	self.isConfig = True	4169	self.isConfig = True
	4016		4170
	@@ -4024,6 +4178,7 class SpectralFitting(Operation):
	4024	my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver = self.__DiffCoherent(jspectra, jcspectra, dataOut, noise, snrth,coh_th, hei_th)	4178	my_incoh_spectra ,my_incoh_cspectra,my_incoh_aver,my_coh_aver, incoh_spectra, coh_spectra, incoh_cspectra, coh_cspectra, incoh_aver, coh_aver = self.__DiffCoherent(jspectra, jcspectra, dataOut, noise, snrth,coh_th, hei_th)
	4025		4179
	4026	clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)	4180	clean_coh_spectra, clean_coh_cspectra, clean_coh_aver = self.__CleanCoherent(snrth, coh_spectra, coh_cspectra, coh_aver, dataOut, noise,1,index)
			4181
	4027	dataOut.data_spc = incoh_spectra	4182	dataOut.data_spc = incoh_spectra
	4028	dataOut.data_cspc = incoh_cspectra	4183	dataOut.data_cspc = incoh_cspectra
	4029	dataOut.sat_spectra = sat_spectra	4184	dataOut.sat_spectra = sat_spectra
	@@ -4138,7 +4293,7 class SpectralFitting(Operation):
	4138	dataOut.nIncohInt= int(dataOut.nIncohInt)	4293	dataOut.nIncohInt= int(dataOut.nIncohInt)
	4139	dataOut.nProfiles = int(dataOut.nProfiles)	4294	dataOut.nProfiles = int(dataOut.nProfiles)
	4140	dataOut.nCohInt = int(dataOut.nCohInt)	4295	dataOut.nCohInt = int(dataOut.nCohInt)
	4141	print('sale',dataOut.nProfiles,dataOut.nHeights)	4296	#print('sale',dataOut.nProfiles,dataOut.nHeights)
	4142	#dataOut.nFFTPoints=nprofs	4297	#dataOut.nFFTPoints=nprofs
	4143	#dataOut.normFactor = nprofs	4298	#dataOut.normFactor = nprofs
	4144	dataOut.channelList = channelList	4299	dataOut.channelList = channelList
	@@ -4148,7 +4303,7 class SpectralFitting(Operation):
	4148	vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds)	4303	vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds)
	4149	#dataOut.ippSeconds=ipp	4304	#dataOut.ippSeconds=ipp
	4150	absc = vmax*( numpy.arange(nProf,dtype='float')-nProf/2.)/nProf	4305	absc = vmax*( numpy.arange(nProf,dtype='float')-nProf/2.)/nProf
	4151	print('sale 2',dataOut.ippSeconds,M,N)	4306	#print('sale 2',dataOut.ippSeconds,M,N)
	4152	print('Empieza procesamiento offline')	4307	print('Empieza procesamiento offline')
	4153	if path != None:	4308	if path != None:
	4154	sys.path.append(path)	4309	sys.path.append(path)
	@@ -4170,6 +4325,8 class SpectralFitting(Operation):
	4170	dataOut.data_snr1_i = numpy.zeros((nGroups2, nHeights))numpy.nan	4325	dataOut.data_snr1_i = numpy.zeros((nGroups2, nHeights))numpy.nan
	4171	# dataOut.smooth_i = numpy.zeros((nGroups2, nHeights))numpy.nan	4326	# dataOut.smooth_i = numpy.zeros((nGroups2, nHeights))numpy.nan
	4172		4327
			4328	if self.aux == 1:
			4329	dataOut.p03last = numpy.zeros(nGroups, 'float32')
	4173	for i in range(nGroups):	4330	for i in range(nGroups):
	4174	coord = groupArray[i,:]	4331	coord = groupArray[i,:]
	4175	#Input data array	4332	#Input data array
	@@ -4193,7 +4350,7 class SpectralFitting(Operation):
	4193	n1= n1i	4350	n1= n1i
	4194	my_noises[2*i+0] = n0	4351	my_noises[2*i+0] = n0
	4195	my_noises[2*i+1] = n1	4352	my_noises[2*i+1] = n1
	4196	snrth = -13 #-13.0 # -4 -16 -25	4353	snrth = -12.0 #-11.0 # -4 -16 -25
	4197	snrth = 10**(snrth/10.0)	4354	snrth = 10**(snrth/10.0)
	4198	jvelr = numpy.zeros(nHeights, dtype = 'float')	4355	jvelr = numpy.zeros(nHeights, dtype = 'float')
	4199	#snr0 = numpy.zeros(nHeights, dtype = 'float')	4356	#snr0 = numpy.zeros(nHeights, dtype = 'float')
	@@ -4201,7 +4358,9 class SpectralFitting(Operation):
	4201	hvalid = [0]	4358	hvalid = [0]
	4202		4359
	4203	coh2 = abs(dataOut.data_cspc[i,1:nProf,:])*2/(dataOut.data_spc[0+i2,1:nProf-0,:]dataOut.data_spc[1+i2,1:nProf-0,:])	4360	coh2 = abs(dataOut.data_cspc[i,1:nProf,:])*2/(dataOut.data_spc[0+i2,1:nProf-0,:]dataOut.data_spc[1+i2,1:nProf-0,:])
	4204		4361
			4362	SIGNAL_0 = []
			4363	SIGNAL_1 = []
	4205	for h in range(nHeights):	4364	for h in range(nHeights):
	4206	smooth = clean_num_aver[i+1,h]	4365	smooth = clean_num_aver[i+1,h]
	4207	signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth	4366	signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
	@@ -4228,10 +4387,32 class SpectralFitting(Operation):
	4228	else: jvelr[h] = absc[0]	4387	else: jvelr[h] = absc[0]
	4229	if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)	4388	if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)
	4230	#print(maxp0,absc[maxp0],snr0,jvelr[h])	4389	#print(maxp0,absc[maxp0],snr0,jvelr[h])
			4390	SIGNAL_0.append(signal0)
			4391	SIGNAL_1.append(signal1)
	4231		4392
	4232	if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1	4393	if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1
	4233	else: fd0 = numpy.nan	4394	else: fd0 = numpy.nan
	4234	print(fd0)	4395	print(fd0)
			4396	fd0n = [fd0]
			4397	SIGNAL_0 = numpy.array(SIGNAL_0)
			4398	SIGNAL_1 = numpy.array(SIGNAL_1)
			4399	### ---- Signal correction
			4400
			4401	diff = numpy.abs(dataOut.p03last[i] - fd0)
			4402	if (filec != None) and (self.aux == 0) and (diff > 6):
			4403	print('hace correccion de vr')
			4404	try:
			4405	maxs_0 = self.Vr_correction(SIGNAL_0)
			4406	maxs_1 = self.Vr_correction(SIGNAL_1)
			4407	fd0_aux = -1*(absc[maxs_0] + absc[maxs_1]) / 2
			4408	if numpy.abs(dataOut.p03last[i] - fd0_aux ) > diff: print("Wrong correction omitted", fd0_aux)
			4409	else: fd0 = fd0_aux; print("Changed fd0: ", fd0)
			4410	except Exception as e:
			4411	print("Error in correction processes skiped: ", e)
			4412
			4413	dataOut.p03last[i] = fd0
			4414	### ---
			4415	fd0n = [fd0]
	4235	for h in range(nHeights):	4416	for h in range(nHeights):
	4236	d = data[:,h]	4417	d = data[:,h]
	4237	smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]	4418	smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]
	@@ -4294,6 +4475,16 class SpectralFitting(Operation):
	4294	minp = p0*numpy.nan	4475	minp = p0*numpy.nan
	4295	error0 = numpy.nan	4476	error0 = numpy.nan
	4296	error1 = p0*numpy.nan	4477	error1 = p0*numpy.nan
			4478	# s_sq = (self.__residFunction(minp,dp,LT,constants)).sum()/(len(dp)-len(p0))
			4479	# covp = covp*s_sq
			4480	# error = []
			4481	# for ip in range(len(minp)):
			4482	# try:
			4483	# error.append(numpy.absolute(covp[ip][ip])**0.5)
			4484	# except:
			4485	# error.append( 0.00 )
			4486	#if i==1 and h==11 and index == 139: print(p0, minp,data_spc)
			4487	fd0n = numpy.concatenate((fd0n,minp[3]), axis=None)
	4297	else :	4488	else :
	4298	data_spc = dataOut.data_spc[coord,:,h]	4489	data_spc = dataOut.data_spc[coord,:,h]
	4299	p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))	4490	p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))
	@@ -4308,7 +4499,25 class SpectralFitting(Operation):
	4308	dataOut.data_param[i,:,h] = minp	4499	dataOut.data_param[i,:,h] = minp
	4309	dataOut.data_snr1_i[i*2,h] = numpy.sum(signalpn0/(nProf-1))/n0	4500	dataOut.data_snr1_i[i*2,h] = numpy.sum(signalpn0/(nProf-1))/n0
	4310	dataOut.data_snr1_i[i*2+1,h] = numpy.sum(signalpn1/(nProf-1))/n1	4501	dataOut.data_snr1_i[i*2+1,h] = numpy.sum(signalpn1/(nProf-1))/n1
	4311		4502	#dataOut.smooth_i[i*2,h] = clean_num_aver[i+1,h]
			4503	#print(fd0,dataOut.data_param[i,3,h])
			4504	#print(fd0,dataOut.data_param[i,3,:])
			4505	fd0n = fd0n[1:]
			4506	dvr = abs(dataOut.data_param[i,3,:]-numpy.nanmedian(fd0n))
			4507	indbad=numpy.where(dvr > 25)
			4508	esf = 0
			4509	if filec != None:
			4510	esf = self.weightf.esffit(esf,tini.tm_year,tini.tm_yday,index)
			4511	#print(indbad,'media ',numpy.nanmedian(fd0n), esf)
			4512
			4513	#plt.plot(abs(dataOut.data_param[i,3,:]-numpy.nanmedian(fd0n)),'o')
			4514	#plt.plot(dataOut.data_param[i,3,:])
			4515	#plt.ylim(-200,200)
			4516	if esf == 0 :
			4517	dataOut.data_param[i,3,indbad] = numpy.nanmedian(fd0n)
			4518	print('corrige')
			4519	#plt.plot(dataOut.data_param[i,3,:])
			4520	#plt.show()
	4312	for ht in range(nHeights-1) :	4521	for ht in range(nHeights-1) :
	4313	smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]	4522	smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
	4314	dataOut.data_paramC[4*i,ht,1] = smooth	4523	dataOut.data_paramC[4*i,ht,1] = smooth
	@@ -4393,6 +4602,7 class SpectralFitting(Operation):
	4393	sat_fits[1+2beam,ht] = numpy.sum(signalpn1/(val1_npointsnProf))/n1	4602	sat_fits[1+2beam,ht] = numpy.sum(signalpn1/(val1_npointsnProf))/n1
	4394		4603
	4395	dataOut.sat_fits = sat_fits	4604	dataOut.sat_fits = sat_fits
			4605	self.aux = 0
	4396	return dataOut	4606	return dataOut
	4397		4607
	4398	def __residFunction(self, p, dp, LT, constants):	4608	def __residFunction(self, p, dp, LT, constants):
	@@ -4588,7 +4798,6 class WindProfiler(Operation):
	4588	return distx, disty, dist, ang	4798	return distx, disty, dist, ang
	4589	#Calculo de Matrices	4799	#Calculo de Matrices
	4590		4800
	4591
	4592	def __calculateVelVer(self, phase, lagTRange, _lambda):	4801	def __calculateVelVer(self, phase, lagTRange, _lambda):
	4593		4802
	4594	Ts = lagTRange[1] - lagTRange[0]	4803	Ts = lagTRange[1] - lagTRange[0]
	@@ -4715,8 +4924,6 class WindProfiler(Operation):
	4715	azim = meteorAux[:, 3]*numpy.pi/180	4924	azim = meteorAux[:, 3]*numpy.pi/180
	4716		4925
	4717	n = numpy.cos(zen)	4926	n = numpy.cos(zen)
	4718	# m = (1 - n2)/(1 - numpy.tan(azim)2)
	4719	# l = m*numpy.tan(azim)
	4720	l = numpy.sin(zen)*numpy.sin(azim)	4927	l = numpy.sin(zen)*numpy.sin(azim)
	4721	m = numpy.sin(zen)*numpy.cos(azim)	4928	m = numpy.sin(zen)*numpy.cos(azim)
	4722		4929
	@@ -4781,7 +4988,6 class WindProfiler(Operation):
	4781		4988
	4782	for t in uniqueTime:	4989	for t in uniqueTime:
	4783	metArray1 = metArray[utctime==t,:]	4990	metArray1 = metArray[utctime==t,:]
	4784	# phaseDerThresh = numpy.pi/4 #reducir Phase thresh
	4785	tmet = metArray1[:,1].astype(int)	4991	tmet = metArray1[:,1].astype(int)
	4786	hmet = metArray1[:,2].astype(int)	4992	hmet = metArray1[:,2].astype(int)
	4787		4993
	@@ -4912,8 +5118,7 class WindProfiler(Operation):
	4912		5118
	4913	def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):	5119	def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
	4914		5120
	4915	param = dataOut.data_param	5121	param = dataOut.moments ### LA VERSION ANTERIOR trabajaba con param = dataOut.data_param
	4916	#if dataOut.abscissaList != None:
	4917	if numpy.any(dataOut.abscissaList):	5122	if numpy.any(dataOut.abscissaList):
	4918	absc = dataOut.abscissaList[:-1]	5123	absc = dataOut.abscissaList[:-1]
	4919	# noise = dataOut.noise	5124	# noise = dataOut.noise
	@@ -4933,29 +5138,9 class WindProfiler(Operation):
	4933	elif technique == 'SA':	5138	elif technique == 'SA':
	4934		5139
	4935	#Parameters	5140	#Parameters
	4936	# position_x = kwargs['positionX']
	4937	# position_y = kwargs['positionY']
	4938	# azimuth = kwargs['azimuth']
	4939	#
	4940	# if kwargs.has_key('crosspairsList'):
	4941	# pairs = kwargs['crosspairsList']
	4942	# else:
	4943	# pairs = None
	4944	#
	4945	# if kwargs.has_key('correctFactor'):
	4946	# correctFactor = kwargs['correctFactor']
	4947	# else:
	4948	# correctFactor = 1
	4949
	4950	# tau = dataOut.data_param
	4951	# _lambda = dataOut.C/dataOut.frequency
	4952	# pairsList = dataOut.groupList
	4953	# nChannels = dataOut.nChannels
	4954
	4955	kwargs['groupList'] = dataOut.groupList	5141	kwargs['groupList'] = dataOut.groupList
	4956	kwargs['tau'] = dataOut.data_param	5142	kwargs['tau'] = dataOut.data_param
	4957	kwargs['_lambda'] = dataOut.C/dataOut.frequency	5143	kwargs['_lambda'] = dataOut.C/dataOut.frequency
	4958	# dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
	4959	dataOut.data_output = self.techniqueSA(kwargs)	5144	dataOut.data_output = self.techniqueSA(kwargs)
	4960	dataOut.utctimeInit = dataOut.utctime	5145	dataOut.utctimeInit = dataOut.utctime
	4961	dataOut.outputInterval = dataOut.timeInterval	5146	dataOut.outputInterval = dataOut.timeInterval
	@@ -4984,7 +5169,6 class WindProfiler(Operation):
	4984	dataOut.outputInterval = nHours*3600	5169	dataOut.outputInterval = nHours*3600
	4985		5170
	4986	if self.__isConfig == False:	5171	if self.__isConfig == False:
	4987	# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
	4988	#Get Initial LTC time	5172	#Get Initial LTC time
	4989	self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)	5173	self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
	4990	self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()	5174	self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
	@@ -5046,7 +5230,6 class WindProfiler(Operation):
	5046		5230
	5047	if self.__isConfig == False:	5231	if self.__isConfig == False:
	5048	dataOut.outputInterval = nMins*60	5232	dataOut.outputInterval = nMins*60
	5049	# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
	5050	#Get Initial LTC time	5233	#Get Initial LTC time
	5051	initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)	5234	initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
	5052	minuteAux = initime.minute	5235	minuteAux = initime.minute
	@@ -5077,69 +5260,19 class WindProfiler(Operation):
	5077	dataOut.flagNoData = False	5260	dataOut.flagNoData = False
	5078	self.__buffer = None	5261	self.__buffer = None
	5079		5262
	5080	return	5263	return dataOut
	5081
	5082	class WindProfiler(Operation):
	5083
	5084	__isConfig = False
	5085
	5086	__initime = None
	5087	__lastdatatime = None
	5088	__integrationtime = None
	5089
	5090	__buffer = None
	5091
	5092	__dataReady = False
	5093
	5094	__firstdata = None
	5095		5264
	5096	n = None	5265	class EWDriftsEstimation(Operation):
	5097		5266
	5098	def __init__(self):	5267	def __init__(self):
	5099	Operation.__init__(self)	5268	Operation.__init__(self)
	5100		5269
	5101	def __calculateCosDir(self, elev, azim):
	5102	zen = (90 - elev)*numpy.pi/180
	5103	azim = azim*numpy.pi/180
	5104	cosDirX = numpy.sqrt((1-numpy.cos(zen)2)/((1+numpy.tan(azim)2)))
	5105	cosDirY = numpy.sqrt(1-numpy.cos(zen)2-cosDirX2)
	5106
	5107	signX = numpy.sign(numpy.cos(azim))
	5108	signY = numpy.sign(numpy.sin(azim))
	5109
	5110	cosDirX = numpy.copysign(cosDirX, signX)
	5111	cosDirY = numpy.copysign(cosDirY, signY)
	5112	return cosDirX, cosDirY
	5113
	5114	def __calculateAngles(self, theta_x, theta_y, azimuth):
	5115
	5116	dir_cosw = numpy.sqrt(1-theta_x2-theta_y2)
	5117	zenith_arr = numpy.arccos(dir_cosw)
	5118	azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
	5119
	5120	dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
	5121	dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
	5122
	5123	return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
	5124
	5125	def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
	5126
	5127	if horOnly:
	5128	A = numpy.c_[dir_cosu,dir_cosv]
	5129	else:
	5130	A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
	5131	A = numpy.asmatrix(A)
	5132	A1 = numpy.linalg.inv(A.transpose()A)A.transpose()
	5133
	5134	return A1
	5135
	5136	def __correctValues(self, heiRang, phi, velRadial, SNR):	5270	def __correctValues(self, heiRang, phi, velRadial, SNR):
	5137	listPhi = phi.tolist()	5271	listPhi = phi.tolist()
	5138	maxid = listPhi.index(max(listPhi))	5272	maxid = listPhi.index(max(listPhi))
	5139	minid = listPhi.index(min(listPhi))	5273	minid = listPhi.index(min(listPhi))
	5140		5274
	5141	rango = list(range(len(phi)))	5275	rango = list(range(len(phi)))
	5142
	5143	heiRang1 = heiRang*math.cos(phi[maxid])	5276	heiRang1 = heiRang*math.cos(phi[maxid])
	5144	heiRangAux = heiRang*math.cos(phi[minid])	5277	heiRangAux = heiRang*math.cos(phi[minid])
	5145	indOut = (heiRang1 < heiRangAux[0]).nonzero()	5278	indOut = (heiRang1 < heiRangAux[0]).nonzero()
	@@ -5151,13 +5284,18 class WindProfiler(Operation):
	5151	for i in rango:	5284	for i in rango:
	5152	x = heiRang*math.cos(phi[i])	5285	x = heiRang*math.cos(phi[i])
	5153	y1 = velRadial[i,:]	5286	y1 = velRadial[i,:]
	5154	f1 = interpolate.interp1d(x,y1,kind = 'cubic')	5287	vali= (numpy.isfinite(y1)==True).nonzero()
	5155		5288	y1=y1[vali]
			5289	x = x[vali]
			5290	f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False)
	5156	x1 = heiRang1	5291	x1 = heiRang1
	5157	y11 = f1(x1)	5292	y11 = f1(x1)
	5158
	5159	y2 = SNR[i,:]	5293	y2 = SNR[i,:]
	5160	f2 = interpolate.interp1d(x,y2,kind = 'cubic')	5294	x = heiRang*math.cos(phi[i])
			5295	vali= (y2 != -1).nonzero()
			5296	y2 = y2[vali]
			5297	x = x[vali]
			5298	f2 = interpolate.interp1d(x,y2,kind = 'cubic',bounds_error=False)
	5161	y21 = f2(x1)	5299	y21 = f2(x1)
	5162		5300
	5163	velRadial1[i,:] = y11	5301	velRadial1[i,:] = y11
	@@ -5165,713 +5303,124 class WindProfiler(Operation):
	5165		5303
	5166	return heiRang1, velRadial1, SNR1	5304	return heiRang1, velRadial1, SNR1
	5167		5305
	5168	def __calculateVelUVW(self, A, velRadial):	5306	def run(self, dataOut, zenith, zenithCorrection,fileDrifts,beam_pos=None):
	5169		5307	dataOut.lat = -11.95
	5170	#Operacion Matricial	5308	dataOut.lon = -76.87
	5171	velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))	5309	dataOut.spcst = 0.00666
	5172	velUVW[:,:] = numpy.dot(A,velRadial)	5310	dataOut.pl = 0.0003
	5173		5311	dataOut.cbadn = 3
	5174		5312	dataOut.inttms = 300
	5175	return velUVW	5313	if numpy.any(beam_pos) :
			5314	dataOut.azw = beam_pos[0] # -115.687
			5315	dataOut.elw = beam_pos[1] # 86.1095
			5316	dataOut.aze = beam_pos[2] # 130.052
			5317	dataOut.ele = beam_pos[3] # 87.6558
			5318	dataOut.jro14 = numpy.log10(dataOut.spc_noise[0]/dataOut.normFactor)
			5319	dataOut.jro15 = numpy.log10(dataOut.spc_noise[1]/dataOut.normFactor)
			5320	dataOut.jro16 = numpy.log10(dataOut.spc_noise[2]/dataOut.normFactor)
			5321	dataOut.nwlos = numpy.log10(dataOut.spc_noise[3]/dataOut.normFactor)
	5176		5322
	5177	def techniqueDBS(self, kwargs):	5323	heiRang = dataOut.heightList
	5178	"""	5324	velRadial = dataOut.data_param[:,3,:]
	5179	Function that implements Doppler Beam Swinging (DBS) technique.	5325	velRadialm = dataOut.data_param[:,2:4,:]*-1
	5180		5326
	5181	Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,	5327	rbufc=dataOut.data_paramC[:,:,0]
	5182	Direction correction (if necessary), Ranges and SNR	5328	ebufc=dataOut.data_paramC[:,:,1]
			5329	SNR = dataOut.data_snr1_i
			5330	rbufi = dataOut.data_snr1_i
			5331	velRerr = dataOut.data_error[:,4,:]
			5332	range1 = dataOut.heightList
			5333	nhei = len(range1)
			5334
			5335	sat_fits = dataOut.sat_fits
			5336
			5337	channels = dataOut.channelList
			5338	nChan = len(channels)
			5339	my_nbeams = nChan/2
			5340	if my_nbeams == 2:
			5341	moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]],[velRadialm[1,:]],[velRadialm[1,:]]))
			5342	else :
			5343	moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]]))
			5344	dataOut.moments=moments
			5345	#Incoherent
			5346	smooth_w = dataOut.clean_num_aver[0,:]
			5347	chisq_w = dataOut.data_error[0,0,:]
			5348	p_w0 = rbufi[0,:]
			5349	p_w1 = rbufi[1,:]
	5183		5350
	5184	Output: Winds estimation (Zonal, Meridional and Vertical)	5351	# Coherent
			5352	smooth_wC = ebufc[0,:]
			5353	p_w0C = rbufc[0,:]
			5354	p_w1C = rbufc[1,:]
			5355	w_wC = rbufc[2,:]-1 #radial_sign(radial EQ 1)
			5356	t_wC = rbufc[3,:]
			5357	val = (numpy.isfinite(p_w0)==False).nonzero()
			5358	p_w0[val]=0
			5359	val = (numpy.isfinite(p_w1)==False).nonzero()
			5360	p_w1[val]=0
			5361	val = (numpy.isfinite(p_w0C)==False).nonzero()
			5362	p_w0C[val]=0
			5363	val = (numpy.isfinite(p_w1C)==False).nonzero()
			5364	p_w1C[val]=0
			5365	val = (numpy.isfinite(smooth_w)==False).nonzero()
			5366	smooth_w[val]=0
			5367	val = (numpy.isfinite(smooth_wC)==False).nonzero()
			5368	smooth_wC[val]=0
			5369
			5370	p_w0_a = (p_w0smooth_w+p_w0Csmooth_wC)/(smooth_w+smooth_wC)
			5371	p_w1_a = (p_w1smooth_w+p_w1Csmooth_wC)/(smooth_w+smooth_wC)
	5185		5372
	5186	Parameters affected: Winds, height range, SNR	5373	if len(sat_fits) >0 :
	5187	"""	5374	p_w0C = p_w0C + sat_fits[0,:]
	5188	velRadial0 = kwargs['velRadial']	5375	p_w1C = p_w1C + sat_fits[1,:]
	5189	heiRang = kwargs['heightList']	5376
	5190	SNR0 = kwargs['SNR']	5377	if my_nbeams == 1:
			5378	w = velRadial[0,:]
			5379	winds = velRadial.copy()
			5380	w_err = velRerr[0,:]
			5381	werrtmp = numpy.where(numpy.isfinite(w)==True,w_err,numpy.nan)
			5382	w_err = werrtmp.copy()
			5383	u = w*numpy.nan
			5384	u_err = w_err*numpy.nan
			5385	p_e0 = p_w0*numpy.nan
			5386	p_e1 = p_w1*numpy.nan
			5387	#snr1 = 10*numpy.log10(SNR[0])
			5388	if my_nbeams == 2:
	5191		5389
	5192	if 'dirCosx' in kwargs and 'dirCosy' in kwargs:	5390	zenith = numpy.array(zenith)
	5193	theta_x = numpy.array(kwargs['dirCosx'])	5391	zenith -= zenithCorrection
	5194	theta_y = numpy.array(kwargs['dirCosy'])	5392	zenith *= numpy.pi/180
	5195	else:	5393	if zenithCorrection != 0 :
	5196	elev = numpy.array(kwargs['elevation'])	5394	heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
	5197	azim = numpy.array(kwargs['azimuth'])	5395	else :
	5198	theta_x, theta_y = self.__calculateCosDir(elev, azim)	5396	heiRang1 = heiRang
	5199	azimuth = kwargs['correctAzimuth']	5397	velRadial1 = velRadial
	5200	if 'horizontalOnly' in kwargs:	5398	SNR1 = SNR
	5201	horizontalOnly = kwargs['horizontalOnly']	5399
	5202	else: horizontalOnly = False	5400	alp = zenith[0]
	5203	if 'correctFactor' in kwargs:	5401	bet = zenith[1]
	5204	correctFactor = kwargs['correctFactor']
	5205	else: correctFactor = 1
	5206	if 'channelList' in kwargs:
	5207	channelList = kwargs['channelList']
	5208	if len(channelList) == 2:
	5209	horizontalOnly = True
	5210	arrayChannel = numpy.array(channelList)
	5211	param = param[arrayChannel,:,:]
	5212	theta_x = theta_x[arrayChannel]
	5213	theta_y = theta_y[arrayChannel]
	5214		5402
	5215	azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)	5403	w_w = velRadial1[0,:]
	5216	heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)	5404	w_e = velRadial1[1,:]
	5217	A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)	5405	w_w_err = velRerr[0,:]
			5406	w_e_err = velRerr[1,:]
			5407	smooth_e = dataOut.clean_num_aver[2,:]
			5408	chisq_e = dataOut.data_error[1,0,:]
			5409	p_e0 = rbufi[2,:]
			5410	p_e1 = rbufi[3,:]
	5218		5411
	5219	#Calculo de Componentes de la velocidad con DBS	5412	tini=time.localtime(dataOut.utctime)
	5220	winds = self.__calculateVelUVW(A,velRadial1)
	5221		5413
	5222	return winds, heiRang1, SNR1	5414	if tini[3] >= 6 and tini[3] < 18 :
	5223		5415	w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w)
	5224	def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):	5416	w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err)
	5225		5417	else:
	5226	nPairs = len(pairs_ccf)	5418	w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w)
	5227	posx = numpy.asarray(posx)	5419	w_wtmp = numpy.where(range1 > 200,w_w,w_wtmp)
	5228	posy = numpy.asarray(posy)	5420	w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err)
	5229		5421	w_w_errtmp = numpy.where(range1 > 200,w_w_err,w_w_errtmp)
	5230	#Rotacion Inversa para alinear con el azimuth	5422	w_w = w_wtmp
	5231	if azimuth!= None:	5423	w_w_err = w_w_errtmp
	5232	azimuth = azimuth*math.pi/180
	5233	posx1 = posxmath.cos(azimuth) + posymath.sin(azimuth)
	5234	posy1 = -posxmath.sin(azimuth) + posymath.cos(azimuth)
	5235	else:
	5236	posx1 = posx
	5237	posy1 = posy
	5238
	5239	#Calculo de Distancias
	5240	distx = numpy.zeros(nPairs)
	5241	disty = numpy.zeros(nPairs)
	5242	dist = numpy.zeros(nPairs)
	5243	ang = numpy.zeros(nPairs)
	5244
	5245	for i in range(nPairs):
	5246	distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
	5247	disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
	5248	dist[i] = numpy.sqrt(distx[i]2 + disty[i]2)
	5249	ang[i] = numpy.arctan2(disty[i],distx[i])
	5250
	5251	return distx, disty, dist, ang
	5252	#Calculo de Matrices
	5253
	5254	def __calculateVelVer(self, phase, lagTRange, _lambda):
	5255
	5256	Ts = lagTRange[1] - lagTRange[0]
	5257	velW = -_lambdaphase/(4math.pi*Ts)
	5258
	5259	return velW
	5260
	5261	def __calculateVelHorDir(self, dist, tau1, tau2, ang):
	5262	nPairs = tau1.shape[0]
	5263	nHeights = tau1.shape[1]
	5264	vel = numpy.zeros((nPairs,3,nHeights))
	5265	dist1 = numpy.reshape(dist, (dist.size,1))
	5266
	5267	angCos = numpy.cos(ang)
	5268	angSin = numpy.sin(ang)
	5269
	5270	vel0 = dist1tau1/(2tau2**2)
	5271	vel[:,0,:] = (vel0*angCos).sum(axis = 1)
	5272	vel[:,1,:] = (vel0*angSin).sum(axis = 1)
	5273
	5274	ind = numpy.where(numpy.isinf(vel))
	5275	vel[ind] = numpy.nan
	5276
	5277	return vel
	5278
	5279	def techniqueSA(self, kwargs):
	5280
	5281	"""
	5282	Function that implements Spaced Antenna (SA) technique.
	5283
	5284	Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
	5285	Direction correction (if necessary), Ranges and SNR
	5286
	5287	Output: Winds estimation (Zonal, Meridional and Vertical)
	5288
	5289	Parameters affected: Winds
	5290	"""
	5291	position_x = kwargs['positionX']
	5292	position_y = kwargs['positionY']
	5293	azimuth = kwargs['azimuth']
	5294
	5295	if 'correctFactor' in kwargs:
	5296	correctFactor = kwargs['correctFactor']
	5297	else:
	5298	correctFactor = 1
	5299
	5300	groupList = kwargs['groupList']
	5301	pairs_ccf = groupList[1]
	5302	tau = kwargs['tau']
	5303	_lambda = kwargs['_lambda']
	5304
	5305	#Cross Correlation pairs obtained
	5306
	5307	indtau = tau.shape[0]/2
	5308	tau1 = tau[:indtau,:]
	5309	tau2 = tau[indtau:-1,:]
	5310	phase1 = tau[-1,:]
	5311
	5312	#---------------------------------------------------------------------
	5313	#Metodo Directo
	5314	distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
	5315	winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
	5316	winds = stats.nanmean(winds, axis=0)
	5317	#---------------------------------------------------------------------
	5318	#Metodo General
	5319
	5320	#---------------------------------------------------------------------
	5321	winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
	5322	winds = correctFactor*winds
	5323	return winds
	5324
	5325	def __checkTime(self, currentTime, paramInterval, outputInterval):
	5326
	5327	dataTime = currentTime + paramInterval
	5328	deltaTime = dataTime - self.__initime
	5329
	5330	if deltaTime >= outputInterval or deltaTime < 0:
	5331	self.__dataReady = True
	5332	return
	5333
	5334	def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
	5335	'''
	5336	Function that implements winds estimation technique with detected meteors.
	5337
	5338	Input: Detected meteors, Minimum meteor quantity to wind estimation
	5339
	5340	Output: Winds estimation (Zonal and Meridional)
	5341
	5342	Parameters affected: Winds
	5343	'''
	5344	#Settings
	5345	nInt = (heightMax - heightMin)/2
	5346	nInt = int(nInt)
	5347	winds = numpy.zeros((2,nInt))*numpy.nan
	5348
	5349	#Filter errors
	5350	error = numpy.where(arrayMeteor[:,-1] == 0)[0]
	5351	finalMeteor = arrayMeteor[error,:]
	5352
	5353	#Meteor Histogram
	5354	finalHeights = finalMeteor[:,2]
	5355	hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
	5356	nMeteorsPerI = hist[0]
	5357	heightPerI = hist[1]
	5358
	5359	#Sort of meteors
	5360	indSort = finalHeights.argsort()
	5361	finalMeteor2 = finalMeteor[indSort,:]
	5362
	5363	# Calculating winds
	5364	ind1 = 0
	5365	ind2 = 0
	5366
	5367	for i in range(nInt):
	5368	nMet = nMeteorsPerI[i]
	5369	ind1 = ind2
	5370	ind2 = ind1 + nMet
	5371
	5372	meteorAux = finalMeteor2[ind1:ind2,:]
	5373
	5374	if meteorAux.shape[0] >= meteorThresh:
	5375	vel = meteorAux[:, 6]
	5376	zen = meteorAux[:, 4]*numpy.pi/180
	5377	azim = meteorAux[:, 3]*numpy.pi/180
	5378
	5379	n = numpy.cos(zen)
	5380	l = numpy.sin(zen)*numpy.sin(azim)
	5381	m = numpy.sin(zen)*numpy.cos(azim)
	5382
	5383	A = numpy.vstack((l, m)).transpose()
	5384	A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
	5385	windsAux = numpy.dot(A1, vel)
	5386
	5387	winds[0,i] = windsAux[0]
	5388	winds[1,i] = windsAux[1]
	5389
	5390	return winds, heightPerI[:-1]
	5391
	5392	def techniqueNSM_SA(self, **kwargs):
	5393	metArray = kwargs['metArray']
	5394	heightList = kwargs['heightList']
	5395	timeList = kwargs['timeList']
	5396
	5397	rx_location = kwargs['rx_location']
	5398	groupList = kwargs['groupList']
	5399	azimuth = kwargs['azimuth']
	5400	dfactor = kwargs['dfactor']
	5401	k = kwargs['k']
	5402
	5403	azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
	5404	d = dist*dfactor
	5405	#Phase calculation
	5406	metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
	5407
	5408	metArray1[:,-2] = metArray1[:,-2]metArray1[:,2]1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
	5409
	5410	velEst = numpy.zeros((heightList.size,2))*numpy.nan
	5411	azimuth1 = azimuth1*numpy.pi/180
	5412
	5413	for i in range(heightList.size):
	5414	h = heightList[i]
	5415	indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
	5416	metHeight = metArray1[indH,:]
	5417	if metHeight.shape[0] >= 2:
	5418	velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities
	5419	iazim = metHeight[:,1].astype(int)
	5420	azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths
	5421	A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
	5422	A = numpy.asmatrix(A)
	5423	A1 = numpy.linalg.pinv(A.transpose()A)A.transpose()
	5424	velHor = numpy.dot(A1,velAux)
	5425
	5426	velEst[i,:] = numpy.squeeze(velHor)
	5427	return velEst
	5428
	5429	def __getPhaseSlope(self, metArray, heightList, timeList):
	5430	meteorList = []
	5431	#utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
	5432	#Putting back together the meteor matrix
	5433	utctime = metArray[:,0]
	5434	uniqueTime = numpy.unique(utctime)
	5435
	5436	phaseDerThresh = 0.5
	5437	ippSeconds = timeList[1] - timeList[0]
	5438	sec = numpy.where(timeList>1)[0][0]
	5439	nPairs = metArray.shape[1] - 6
	5440	nHeights = len(heightList)
	5441
	5442	for t in uniqueTime:
	5443	metArray1 = metArray[utctime==t,:]
	5444	tmet = metArray1[:,1].astype(int)
	5445	hmet = metArray1[:,2].astype(int)
	5446
	5447	metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
	5448	metPhase[:,:] = numpy.nan
	5449	metPhase[:,hmet,tmet] = metArray1[:,6:].T
	5450
	5451	#Delete short trails
	5452	metBool = ~numpy.isnan(metPhase[0,:,:])
	5453	heightVect = numpy.sum(metBool, axis = 1)
	5454	metBool[heightVect<sec,:] = False
	5455	metPhase[:,heightVect<sec,:] = numpy.nan
	5456
	5457	#Derivative
	5458	metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
	5459	phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
	5460	metPhase[phDerAux] = numpy.nan
	5461
	5462	#--------------------------METEOR DETECTION -----------------------------------------
	5463	indMet = numpy.where(numpy.any(metBool,axis=1))[0]
	5464
	5465	for p in numpy.arange(nPairs):
	5466	phase = metPhase[p,:,:]
	5467	phDer = metDer[p,:,:]
	5468
	5469	for h in indMet:
	5470	height = heightList[h]
	5471	phase1 = phase[h,:] #82
	5472	phDer1 = phDer[h,:]
	5473
	5474	phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap
	5475
	5476	indValid = numpy.where(~numpy.isnan(phase1))[0]
	5477	initMet = indValid[0]
	5478	endMet = 0
	5479
	5480	for i in range(len(indValid)-1):
	5481
	5482	#Time difference
	5483	inow = indValid[i]
	5484	inext = indValid[i+1]
	5485	idiff = inext - inow
	5486	#Phase difference
	5487	phDiff = numpy.abs(phase1[inext] - phase1[inow])
	5488
	5489	if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor
	5490	sizeTrail = inow - initMet + 1
	5491	if sizeTrail>3*sec: #Too short meteors
	5492	x = numpy.arange(initMet,inow+1)*ippSeconds
	5493	y = phase1[initMet:inow+1]
	5494	ynnan = ~numpy.isnan(y)
	5495	x = x[ynnan]
	5496	y = y[ynnan]
	5497	slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
	5498	ylin = x*slope + intercept
	5499	rsq = r_value**2
	5500	if rsq > 0.5:
	5501	vel = slope#height1000/(k*d)
	5502	estAux = numpy.array([utctime,p,height, vel, rsq])
	5503	meteorList.append(estAux)
	5504	initMet = inext
	5505	metArray2 = numpy.array(meteorList)
	5506
	5507	return metArray2
	5508
	5509	def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
	5510
	5511	azimuth1 = numpy.zeros(len(pairslist))
	5512	dist = numpy.zeros(len(pairslist))
	5513
	5514	for i in range(len(rx_location)):
	5515	ch0 = pairslist[i][0]
	5516	ch1 = pairslist[i][1]
	5517
	5518	diffX = rx_location[ch0][0] - rx_location[ch1][0]
	5519	diffY = rx_location[ch0][1] - rx_location[ch1][1]
	5520	azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
	5521	dist[i] = numpy.sqrt(diffX2 + diffY2)
	5522
	5523	azimuth1 -= azimuth0
	5524	return azimuth1, dist
	5525
	5526	def techniqueNSM_DBS(self, **kwargs):
	5527	metArray = kwargs['metArray']
	5528	heightList = kwargs['heightList']
	5529	timeList = kwargs['timeList']
	5530	azimuth = kwargs['azimuth']
	5531	theta_x = numpy.array(kwargs['theta_x'])
	5532	theta_y = numpy.array(kwargs['theta_y'])
	5533
	5534	utctime = metArray[:,0]
	5535	cmet = metArray[:,1].astype(int)
	5536	hmet = metArray[:,3].astype(int)
	5537	SNRmet = metArray[:,4]
	5538	vmet = metArray[:,5]
	5539	spcmet = metArray[:,6]
	5540
	5541	nChan = numpy.max(cmet) + 1
	5542	nHeights = len(heightList)
	5543
	5544	azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
	5545	hmet = heightList[hmet]
	5546	h1met = hmet*numpy.cos(zenith_arr[cmet]) #Corrected heights
	5547
	5548	velEst = numpy.zeros((heightList.size,2))*numpy.nan
	5549
	5550	for i in range(nHeights - 1):
	5551	hmin = heightList[i]
	5552	hmax = heightList[i + 1]
	5553
	5554	thisH = (h1met>=hmin) & (h1met<hmax) & (cmet!=2) & (SNRmet>8) & (vmet<50) & (spcmet<10)
	5555	indthisH = numpy.where(thisH)
	5556
	5557	if numpy.size(indthisH) > 3:
	5558
	5559	vel_aux = vmet[thisH]
	5560	chan_aux = cmet[thisH]
	5561	cosu_aux = dir_cosu[chan_aux]
	5562	cosv_aux = dir_cosv[chan_aux]
	5563	cosw_aux = dir_cosw[chan_aux]
	5564
	5565	nch = numpy.size(numpy.unique(chan_aux))
	5566	if nch > 1:
	5567	A = self.__calculateMatA(cosu_aux, cosv_aux, cosw_aux, True)
	5568	velEst[i,:] = numpy.dot(A,vel_aux)
	5569
	5570	return velEst
	5571
	5572	def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
	5573
	5574	param = dataOut.moments
	5575	if numpy.any(dataOut.abscissaList):
	5576	absc = dataOut.abscissaList[:-1]
	5577	# noise = dataOut.noise
	5578	heightList = dataOut.heightList
	5579	SNR = dataOut.data_snr
	5580
	5581	if technique == 'DBS':
	5582
	5583	kwargs['velRadial'] = param[:,1,:] #Radial velocity
	5584	kwargs['heightList'] = heightList
	5585	kwargs['SNR'] = SNR
	5586
	5587	dataOut.data_output, dataOut.heightList, dataOut.data_snr = self.techniqueDBS(kwargs) #DBS Function
	5588	dataOut.utctimeInit = dataOut.utctime
	5589	dataOut.outputInterval = dataOut.paramInterval
	5590
	5591	elif technique == 'SA':
	5592
	5593	#Parameters
	5594	kwargs['groupList'] = dataOut.groupList
	5595	kwargs['tau'] = dataOut.data_param
	5596	kwargs['_lambda'] = dataOut.C/dataOut.frequency
	5597	dataOut.data_output = self.techniqueSA(kwargs)
	5598	dataOut.utctimeInit = dataOut.utctime
	5599	dataOut.outputInterval = dataOut.timeInterval
	5600
	5601	elif technique == 'Meteors':
	5602	dataOut.flagNoData = True
	5603	self.__dataReady = False
	5604
	5605	if 'nHours' in kwargs:
	5606	nHours = kwargs['nHours']
	5607	else:
	5608	nHours = 1
	5609
	5610	if 'meteorsPerBin' in kwargs:
	5611	meteorThresh = kwargs['meteorsPerBin']
	5612	else:
	5613	meteorThresh = 6
	5614
	5615	if 'hmin' in kwargs:
	5616	hmin = kwargs['hmin']
	5617	else: hmin = 70
	5618	if 'hmax' in kwargs:
	5619	hmax = kwargs['hmax']
	5620	else: hmax = 110
	5621
	5622	dataOut.outputInterval = nHours*3600
	5623
	5624	if self.__isConfig == False:
	5625	#Get Initial LTC time
	5626	self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
	5627	self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
	5628
	5629	self.__isConfig = True
	5630
	5631	if self.__buffer is None:
	5632	self.__buffer = dataOut.data_param
	5633	self.__firstdata = copy.copy(dataOut)
	5634
	5635	else:
	5636	self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
	5637
	5638	self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
	5639
	5640	if self.__dataReady:
	5641	dataOut.utctimeInit = self.__initime
	5642
	5643	self.__initime += dataOut.outputInterval #to erase time offset
	5644
	5645	dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
	5646	dataOut.flagNoData = False
	5647	self.__buffer = None
	5648
	5649	elif technique == 'Meteors1':
	5650	dataOut.flagNoData = True
	5651	self.__dataReady = False
	5652
	5653	if 'nMins' in kwargs:
	5654	nMins = kwargs['nMins']
	5655	else: nMins = 20
	5656	if 'rx_location' in kwargs:
	5657	rx_location = kwargs['rx_location']
	5658	else: rx_location = [(0,1),(1,1),(1,0)]
	5659	if 'azimuth' in kwargs:
	5660	azimuth = kwargs['azimuth']
	5661	else: azimuth = 51.06
	5662	if 'dfactor' in kwargs:
	5663	dfactor = kwargs['dfactor']
	5664	if 'mode' in kwargs:
	5665	mode = kwargs['mode']
	5666	if 'theta_x' in kwargs:
	5667	theta_x = kwargs['theta_x']
	5668	if 'theta_y' in kwargs:
	5669	theta_y = kwargs['theta_y']
	5670	else: mode = 'SA'
	5671
	5672	#Borrar luego esto
	5673	if dataOut.groupList is None:
	5674	dataOut.groupList = [(0,1),(0,2),(1,2)]
	5675	groupList = dataOut.groupList
	5676	C = 3e8
	5677	freq = 50e6
	5678	lamb = C/freq
	5679	k = 2*numpy.pi/lamb
	5680
	5681	timeList = dataOut.abscissaList
	5682	heightList = dataOut.heightList
	5683
	5684	if self.__isConfig == False:
	5685	dataOut.outputInterval = nMins*60
	5686	#Get Initial LTC time
	5687	initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
	5688	minuteAux = initime.minute
	5689	minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
	5690	self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
	5691
	5692	self.__isConfig = True
	5693
	5694	if self.__buffer is None:
	5695	self.__buffer = dataOut.data_param
	5696	self.__firstdata = copy.copy(dataOut)
	5697
	5698	else:
	5699	self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
	5700
	5701	self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
	5702
	5703	if self.__dataReady:
	5704	dataOut.utctimeInit = self.__initime
	5705	self.__initime += dataOut.outputInterval #to erase time offset
	5706
	5707	metArray = self.__buffer
	5708	if mode == 'SA':
	5709	dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
	5710	elif mode == 'DBS':
	5711	dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList, azimuth=azimuth, theta_x=theta_x, theta_y=theta_y)
	5712	dataOut.data_output = dataOut.data_output.T
	5713	dataOut.flagNoData = False
	5714	self.__buffer = None
	5715
	5716	return dataOut
	5717
	5718	class EWDriftsEstimation(Operation):
	5719
	5720	def __init__(self):
	5721	Operation.__init__(self)
	5722
	5723	def __correctValues(self, heiRang, phi, velRadial, SNR):
	5724	listPhi = phi.tolist()
	5725	maxid = listPhi.index(max(listPhi))
	5726	minid = listPhi.index(min(listPhi))
	5727
	5728	rango = list(range(len(phi)))
	5729	heiRang1 = heiRang*math.cos(phi[maxid])
	5730	heiRangAux = heiRang*math.cos(phi[minid])
	5731	indOut = (heiRang1 < heiRangAux[0]).nonzero()
	5732	heiRang1 = numpy.delete(heiRang1,indOut)
	5733
	5734	velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
	5735	SNR1 = numpy.zeros([len(phi),len(heiRang1)])
	5736
	5737	for i in rango:
	5738	x = heiRang*math.cos(phi[i])
	5739	y1 = velRadial[i,:]
	5740	vali= (numpy.isfinite(y1)==True).nonzero()
	5741	y1=y1[vali]
	5742	x = x[vali]
	5743	f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False)
	5744	x1 = heiRang1
	5745	y11 = f1(x1)
	5746	y2 = SNR[i,:]
	5747	x = heiRang*math.cos(phi[i])
	5748	vali= (y2 != -1).nonzero()
	5749	y2 = y2[vali]
	5750	x = x[vali]
	5751	f2 = interpolate.interp1d(x,y2,kind = 'cubic',bounds_error=False)
	5752	y21 = f2(x1)
	5753
	5754	velRadial1[i,:] = y11
	5755	SNR1[i,:] = y21
	5756
	5757	return heiRang1, velRadial1, SNR1
	5758
	5759	def run(self, dataOut, zenith, zenithCorrection,fileDrifts):
	5760
	5761	dataOut.lat = -11.95
	5762	dataOut.lon = -76.87
	5763	dataOut.spcst = 0.00666
	5764	dataOut.pl = 0.0003
	5765	dataOut.cbadn = 3
	5766	dataOut.inttms = 300
	5767	dataOut.azw = -115.687
	5768	dataOut.elw = 86.1095
	5769	dataOut.aze = 130.052
	5770	dataOut.ele = 87.6558
	5771	dataOut.jro14 = numpy.log10(dataOut.spc_noise[0]/dataOut.normFactor)
	5772	dataOut.jro15 = numpy.log10(dataOut.spc_noise[1]/dataOut.normFactor)
	5773	dataOut.jro16 = numpy.log10(dataOut.spc_noise[2]/dataOut.normFactor)
	5774	dataOut.nwlos = numpy.log10(dataOut.spc_noise[3]/dataOut.normFactor)
	5775
	5776	heiRang = dataOut.heightList
	5777	velRadial = dataOut.data_param[:,3,:]
	5778	velRadialm = dataOut.data_param[:,2:4,:]*-1
	5779
	5780	rbufc=dataOut.data_paramC[:,:,0]
	5781	ebufc=dataOut.data_paramC[:,:,1]
	5782	SNR = dataOut.data_snr1_i
	5783	rbufi = dataOut.data_snr1_i
	5784	velRerr = dataOut.data_error[:,4,:]
	5785	range1 = dataOut.heightList
	5786	nhei = len(range1)
	5787
	5788	sat_fits = dataOut.sat_fits
	5789
	5790	channels = dataOut.channelList
	5791	nChan = len(channels)
	5792	my_nbeams = nChan/2
	5793	if my_nbeams == 2:
	5794	moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]],[velRadialm[1,:]],[velRadialm[1,:]]))
	5795	else :
	5796	moments=numpy.vstack(([velRadialm[0,:]],[velRadialm[0,:]]))
	5797	dataOut.moments=moments
	5798	#Incoherent
	5799	smooth_w = dataOut.clean_num_aver[0,:]
	5800	chisq_w = dataOut.data_error[0,0,:]
	5801	p_w0 = rbufi[0,:]
	5802	p_w1 = rbufi[1,:]
	5803
	5804	# Coherent
	5805	smooth_wC = ebufc[0,:]
	5806	p_w0C = rbufc[0,:]
	5807	p_w1C = rbufc[1,:]
	5808	w_wC = rbufc[2,:]-1 #radial_sign(radial EQ 1)
	5809	t_wC = rbufc[3,:]
	5810	val = (numpy.isfinite(p_w0)==False).nonzero()
	5811	p_w0[val]=0
	5812	val = (numpy.isfinite(p_w1)==False).nonzero()
	5813	p_w1[val]=0
	5814	val = (numpy.isfinite(p_w0C)==False).nonzero()
	5815	p_w0C[val]=0
	5816	val = (numpy.isfinite(p_w1C)==False).nonzero()
	5817	p_w1C[val]=0
	5818	val = (numpy.isfinite(smooth_w)==False).nonzero()
	5819	smooth_w[val]=0
	5820	val = (numpy.isfinite(smooth_wC)==False).nonzero()
	5821	smooth_wC[val]=0
	5822
	5823	#p_w0 = (p_w0smooth_w+p_w0Csmooth_wC)/(smooth_w+smooth_wC)
	5824	#p_w1 = (p_w1smooth_w+p_w1Csmooth_wC)/(smooth_w+smooth_wC)
	5825
	5826	if len(sat_fits) >0 :
	5827	p_w0C = p_w0C + sat_fits[0,:]
	5828	p_w1C = p_w1C + sat_fits[1,:]
	5829
	5830	if my_nbeams == 1:
	5831	w = velRadial[0,:]
	5832	winds = velRadial.copy()
	5833	w_err = velRerr[0,:]
	5834	u = w*numpy.nan
	5835	u_err = w_err*numpy.nan
	5836	p_e0 = p_w0*numpy.nan
	5837	p_e1 = p_w1*numpy.nan
	5838	#snr1 = 10*numpy.log10(SNR[0])
	5839	if my_nbeams == 2:
	5840
	5841	zenith = numpy.array(zenith)
	5842	zenith -= zenithCorrection
	5843	zenith *= numpy.pi/180
	5844	if zenithCorrection != 0 :
	5845	heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
	5846	else :
	5847	heiRang1 = heiRang
	5848	velRadial1 = velRadial
	5849	SNR1 = SNR
	5850
	5851	alp = zenith[0]
	5852	bet = zenith[1]
	5853
	5854	w_w = velRadial1[0,:]
	5855	w_e = velRadial1[1,:]
	5856	w_w_err = velRerr[0,:]
	5857	w_e_err = velRerr[1,:]
	5858	smooth_e = dataOut.clean_num_aver[2,:]
	5859	chisq_e = dataOut.data_error[1,0,:]
	5860	p_e0 = rbufi[2,:]
	5861	p_e1 = rbufi[3,:]
	5862
	5863	tini=time.localtime(dataOut.utctime)
	5864
	5865	if tini[3] >= 6 and tini[3] < 18 :
	5866	w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w)
	5867	w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err)
	5868	else:
	5869	w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w)
	5870	w_wtmp = numpy.where(range1 > 200,w_w,w_wtmp)
	5871	w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err)
	5872	w_w_errtmp = numpy.where(range1 > 200,w_w_err,w_w_errtmp)
	5873	w_w = w_wtmp
	5874	w_w_err = w_w_errtmp
	5875		5424
	5876	#if my_nbeams == 2:	5425	#if my_nbeams == 2:
	5877	smooth_eC=ebufc[4,:]	5426	smooth_eC=ebufc[4,:]
	@@ -5891,9 +5440,9 class EWDriftsEstimation(Operation):
	5891	smooth_e[val]=0	5440	smooth_e[val]=0
	5892	val = (numpy.isfinite(smooth_eC)==False).nonzero()	5441	val = (numpy.isfinite(smooth_eC)==False).nonzero()
	5893	smooth_eC[val]=0	5442	smooth_eC[val]=0
	5894	#p_e0 = (p_e0smooth_e+p_e0Csmooth_eC)/(smooth_e+smooth_eC)	5443	p_e0_a = (p_e0smooth_e+p_e0Csmooth_eC)/(smooth_e+smooth_eC)
	5895	#p_e1 = (p_e1smooth_e+p_e1Csmooth_eC)/(smooth_e+smooth_eC)	5444	p_e1_a = (p_e1smooth_e+p_e1Csmooth_eC)/(smooth_e+smooth_eC)
	5896		5445	#print(w_e,w_eC,p_e0C,sat_fits[2,:])
	5897	if len(sat_fits) >0 :	5446	if len(sat_fits) >0 :
	5898	p_e0C = p_e0C + sat_fits[2,:]	5447	p_e0C = p_e0C + sat_fits[2,:]
	5899	p_e1C = p_e1C + sat_fits[3,:]	5448	p_e1C = p_e1C + sat_fits[3,:]
	@@ -5914,6 +5463,22 class EWDriftsEstimation(Operation):
	5914		5463
	5915	w_err = numpy.sqrt((w_w_errnumpy.sin(bet))2.+(w_e_errnumpy.sin(alp))*2.)/ numpy.absolute(numpy.cos(alp)numpy.sin(bet)-numpy.cos(bet)*numpy.sin(alp))	5464	w_err = numpy.sqrt((w_w_errnumpy.sin(bet))2.+(w_e_errnumpy.sin(alp))*2.)/ numpy.absolute(numpy.cos(alp)numpy.sin(bet)-numpy.cos(bet)*numpy.sin(alp))
	5916	u_err = numpy.sqrt((w_w_errnumpy.cos(bet))2.+(w_e_errnumpy.cos(alp))*2.)/ numpy.absolute(numpy.cos(alp)numpy.sin(bet)-numpy.cos(bet)*numpy.sin(alp))	5465	u_err = numpy.sqrt((w_w_errnumpy.cos(bet))2.+(w_e_errnumpy.cos(alp))*2.)/ numpy.absolute(numpy.cos(alp)numpy.sin(bet)-numpy.cos(bet)*numpy.sin(alp))
			5466
			5467	#wtmp = numpy.where(w < -100 ,numpy.nan,w)
			5468	wtmp = numpy.where(abs(w) > 100 ,numpy.nan,w)
			5469	werrtmp = numpy.where(abs(w_err) > 100 ,numpy.nan,w_err)
			5470	werrtmp = numpy.where(numpy.isfinite(wtmp)==True,werrtmp,numpy.nan)
			5471	#wtmp = numpy.where(numpy.isfinite(werrtmp)==True,wtmp,numpy.nan)
			5472
			5473	#utmp = numpy.where(u < -500,numpy.nan,u)
			5474	utmp = numpy.where(abs(u) > 500 ,numpy.nan,u)
			5475	uerrtmp = numpy.where(abs(u_err) > 500 ,numpy.nan,u_err)
			5476	uerrtmp = numpy.where(numpy.isfinite(utmp)==True,uerrtmp,numpy.nan)
			5477	#utmp = numpy.where(numpy.isfinite(uerrtmp)==True,utmp,numpy.nan)
			5478	w= wtmp.copy()
			5479	u= utmp.copy()
			5480	w_err= werrtmp.copy()
			5481	u_err= uerrtmp.copy()
	5917		5482
	5918	winds = numpy.vstack((w,u))	5483	winds = numpy.vstack((w,u))
	5919	dataOut.heightList = heiRang1	5484	dataOut.heightList = heiRang1
	@@ -5924,7 +5489,7 class EWDriftsEstimation(Operation):
	5924	#print('alt ',range1numpy.sin(86.1numpy.pi/180))	5489	#print('alt ',range1numpy.sin(86.1numpy.pi/180))
	5925	#print(numpy.min([dataOut.eldir7,dataOut.eldir8]))	5490	#print(numpy.min([dataOut.eldir7,dataOut.eldir8]))
	5926	galt = range1numpy.sin(numpy.min([dataOut.elw,dataOut.ele])numpy.pi/180.)	5491	galt = range1numpy.sin(numpy.min([dataOut.elw,dataOut.ele])numpy.pi/180.)
	5927	dataOut.params = numpy.vstack((range1,galt,w,w_err,u,u_err,w_w,w_w_err,w_e,w_e_err,numpy.log10(p_w0),numpy.log10(p_w0C),numpy.log10(p_w1),numpy.log10(p_w1C),numpy.log10(p_e0),numpy.log10(p_e0C),numpy.log10(p_e1),numpy.log10(p_e1C),chisq_w,chisq_e))	5492	dataOut.params = numpy.vstack((range1,galt,w,w_err,u,u_err,w_w,w_w_err,w_e,w_e_err,numpy.log10(p_w0),numpy.log10(p_w0C),numpy.log10(p_w1),numpy.log10(p_w1C),numpy.log10(p_e0),numpy.log10(p_e0C),numpy.log10(p_e1),numpy.log10(p_e1C),chisq_w,chisq_e,numpy.log10(p_w0_a),numpy.log10(p_w1_a),numpy.log10(p_e0_a),numpy.log10(p_e1_a)))
	5928	#snr1 = 10*numpy.log10(SNR1[0])	5493	#snr1 = 10*numpy.log10(SNR1[0])
	5929	#print(min(snr1), max(snr1))	5494	#print(min(snr1), max(snr1))
	5930	snr1 = numpy.vstack((p_w0,p_w1,p_e0,p_e1))	5495	snr1 = numpy.vstack((p_w0,p_w1,p_e0,p_e1))
	@@ -6001,39 +5566,43 class EWDriftsEstimation(Operation):
	6001	uA_err[nhei_avg-1] = sigma	5566	uA_err[nhei_avg-1] = sigma
	6002	uA = uA[6*h_avgs:nhei_avg-0]	5567	uA = uA[6*h_avgs:nhei_avg-0]
	6003	uA_err = uA_err[6*h_avgs:nhei_avg-0]	5568	uA_err = uA_err[6*h_avgs:nhei_avg-0]
	6004	dataOut.drifts_avg = numpy.vstack((wA,uA))	5569	dataOut.drifts_avg = numpy.vstack((wA,uA))
	6005		5570
	6006	if my_nbeams == 1: dataOut.drifts_avg = wA	5571	if my_nbeams == 1: dataOut.drifts_avg = wA
	6007	#deltahavg= wA*0.0+deltah	5572	#deltahavg= wA*0.0+deltah
	6008	dataOut.range = range1	5573	dataOut.range = range1
	6009	galtavg = range_avernumpy.sin(numpy.min([dataOut.elw,dataOut.ele])numpy.pi/180.)	5574	galtavg = range_avernumpy.sin(numpy.min([dataOut.elw,dataOut.ele])numpy.pi/180.)
	6010	dataOut.params_avg = numpy.vstack((wA,wA_err,uA,uA_err,range_aver,galtavg,delta_h))	5575	dataOut.params_avg = numpy.vstack((wA,wA_err,uA,uA_err,range_aver,galtavg,delta_h))
	6011		5576
	6012	#print('comparando dim de avg ',wA.shape,deltahavg.shape,range_aver.shape)	5577	'''
	6013	tini=time.localtime(dataOut.utctime)	5578	tini=time.localtime(dataOut.utctime)
	6014	datefile= str(tini[0]).zfill(4)+str(tini[1]).zfill(2)+str(tini[2]).zfill(2)	5579	datefile= str(tini[0]).zfill(4)+str(tini[1]).zfill(2)+str(tini[2]).zfill(2)
	6015	nfile = fileDrifts+'/jro'+datefile+'drifts_sch3.txt'	5580	nfile = fileDrifts+'/jro'+datefile+'drifts_sch3.txt'
	6016		5581	#nfile = '/home/pcondor/Database/ewdriftsschain2019/jro'+datefile+'drifts_sch3.txt'
			5582	#print(len(dataOut.drifts_avg),dataOut.drifts_avg.shape)
	6017	f1 = open(nfile,'a')	5583	f1 = open(nfile,'a')
	6018	datedriftavg=str(tini[0])+' '+str(tini[1])+' '+str(tini[2])+' '+str(tini[3])+' '+str(tini[4])	5584	datedriftavg=str(tini[0])+' '+str(tini[1])+' '+str(tini[2])+' '+str(tini[3])+' '+str(tini[4])
	6019	driftavgstr=str(dataOut.drifts_avg)	5585	driftavgstr=str(dataOut.drifts_avg)
	6020	numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')	5586	numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')
	6021	numpy.savetxt(f1,numpy.reshape(range_aver,(1,len(range_aver))) ,fmt='%10.2f')	5587	numpy.savetxt(f1,numpy.reshape(range_aver,(1,len(range_aver))) ,fmt='%10.2f')
			5588	#numpy.savetxt(f1,numpy.reshape(dataOut.drifts_avg,(7,len(dataOut.drifts_avg))) ,fmt='%10.2f')
	6022	numpy.savetxt(f1,dataOut.drifts_avg[:,:],fmt='%10.2f')	5589	numpy.savetxt(f1,dataOut.drifts_avg[:,:],fmt='%10.2f')
	6023	f1.close()	5590	f1.close()
	6024		5591	'''
			5592	'''
	6025	swfile = fileDrifts+'/jro'+datefile+'drifts_sw.txt'	5593	swfile = fileDrifts+'/jro'+datefile+'drifts_sw.txt'
	6026	f1 = open(swfile,'a')	5594	f1 = open(swfile,'a')
	6027	numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')	5595	numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')
	6028	numpy.savetxt(f1,numpy.reshape(heiRang,(1,len(heiRang))),fmt='%10.2f')	5596	numpy.savetxt(f1,numpy.reshape(heiRang,(1,len(heiRang))),fmt='%10.2f')
	6029	numpy.savetxt(f1,dataOut.data_param[:,0,:],fmt='%10.2f')	5597	numpy.savetxt(f1,dataOut.data_param[:,0,:],fmt='%10.2f')
	6030	f1.close()	5598	f1.close()
	6031	dataOut.heightListtmp = dataOut.heightList
	6032	'''	5599	'''
			5600	dataOut.heightListtmp = dataOut.heightList
			5601
	6033	#Envio data de drifts a mysql	5602	#Envio data de drifts a mysql
	6034	fechad = str(tini[0]).zfill(4)+'-'+str(tini[1]).zfill(2)+'-'+str(tini[2]).zfill(2)+' '+str(tini[3]).zfill(2)+':'+str(tini[4]).zfill(2)+':'+str(0).zfill(2)	5603	fechad = str(tini[0]).zfill(4)+'-'+str(tini[1]).zfill(2)+'-'+str(tini[2]).zfill(2)+' '+str(tini[3]).zfill(2)+':'+str(tini[4]).zfill(2)+':'+str(0).zfill(2)
	6035	mydb = mysql.connector.connect(	5604	mydb = mysql.connector.connect(
	6036	host="10.10.110.213",	5605	host="10.10.110.205",
	6037	user="user_clima",	5606	user="user_clima",
	6038	password="5D.bh(B2)Y_wRNz9",	5607	password="5D.bh(B2)Y_wRNz9",
	6039	database="clima_espacial"	5608	database="clima_espacial"
	@@ -6056,7 +5625,7 class EWDriftsEstimation(Operation):
	6056	mydb.commit()	5625	mydb.commit()
	6057		5626
	6058	print(mycursor.rowcount, "record inserted.")	5627	print(mycursor.rowcount, "record inserted.")
	6059	'''	5628
	6060	return dataOut	5629	return dataOut
	6061		5630
	6062	class setHeightDrifts(Operation):	5631	class setHeightDrifts(Operation):
	@@ -7317,7 +6886,6 class SMOperations():
	7317	chan1Y = chan1	6886	chan1Y = chan1
	7318	chan2Y = chan2	6887	chan2Y = chan2
	7319	distances[2:4] = numpy.array([d1,d2])	6888	distances[2:4] = numpy.array([d1,d2])
	7320
	7321	pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]	6889	pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
	7322		6890
	7323	return pairslist, distances	6891	return pairslist, distances

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