##// END OF EJS Templates
schain
RSS
UPDATE Repository desarroll drifts_schain, antigua version de WindProfiler en comparar 2.py
avaldez -
r1771:ed72d0635cd3
parent child
Show More
Side by Side Unified
Context file:
r1771:ed72d0635cd3 Loading diff...:
Show file before | Show file after | Hide comments
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-cosDirX**2)
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_x**2-theta_y**2)
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 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
154 posy1 = -posx*math.sin(azimuth) + posy*math.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 = -_lambda*phase/(4*math.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 = dist1*tau1/(2*tau2**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 - n**2)/(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#*height*1000/(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(diffX**2 + diffY**2)
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
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1182 lines changed) Show them Hide them
@@ -8,6 +8,7 import re
8 8 import datetime
9 9 import copy
10 10 import sys
11 import os
11 12 import importlib
12 13 import itertools
13 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 19 from .jroproc_base import ProcessingUnit, Operation, MPDecorator
19 20 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
20 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 23 from scipy.optimize import fmin, curve_fit
23 24 from schainpy.utils import log
24 25 import warnings
25 26 from numpy import NaN
26 27 from scipy.optimize.optimize import OptimizeWarning
28 from schainpy.model.io.jroIO_madrigal import MADWriter
27 29 warnings.filterwarnings('ignore')
28 30
29 31 import os
@@ -31,6 +33,12 import csv
31 33 from scipy import signal
32 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 42 SPEED_OF_LIGHT = 299792458
35 43
36 44 '''solving pickling issue'''
@@ -51,7 +59,7 def _unpickle_method(func_name, obj, cls):
51 59 break
52 60 return func.__get__(obj, cls)
53 61
54
62 # @MPDecorator
55 63 class ParametersProc(ProcessingUnit):
56 64
57 65 METHODS = {}
@@ -171,6 +179,8 class ParametersProc(ProcessingUnit):
171 179 self.dataOut.spc_noise = self.dataIn.getNoise()
172 180 self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
173 181 # self.dataOut.normFactor = self.dataIn.normFactor
182 self.dataOut.pairsList = self.dataIn.pairsList
183 self.dataOut.groupList = self.dataIn.pairsList
174 184 self.dataOut.flagNoData = False
175 185
176 186 if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
@@ -224,6 +234,7 class ParametersProc(ProcessingUnit):
224 234 self.__updateObjFromInput()
225 235 self.dataOut.utctimeInit = self.dataIn.utctime
226 236 self.dataOut.paramInterval = self.dataIn.timeInterval
237
227 238 return
228 239
229 240
@@ -285,6 +296,7 class RemoveWideGC(Operation):
285 296 continue
286 297 junk3 = numpy.squeeze(numpy.diff(j1index))
287 298 junk4 = numpy.squeeze(numpy.diff(j2index))
299
288 300 valleyindex = j2index[numpy.where(junk4>1)]
289 301 peakindex = j1index[numpy.where(junk3>1)]
290 302
@@ -294,6 +306,7 class RemoveWideGC(Operation):
294 306 if numpy.size(isvalid) >1 :
295 307 vindex = numpy.argmax(self.spc[ich,gc_values[peakindex[isvalid]],ir])
296 308 isvalid = isvalid[vindex]
309
297 310 # clutter peak
298 311 gcpeak = peakindex[isvalid]
299 312 vl = numpy.where(valleyindex < gcpeak)
@@ -2655,7 +2668,7 class SpectralMoments(Operation):
2655 2668
2656 2669 if (type1 is None): type1 = 0
2657 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 2672 if (dc is None): dc = 0
2660 2673 if (aliasing is None): aliasing = 0
2661 2674 if (oldfd is None): oldfd = 0
@@ -3190,6 +3203,7 class SpectralFitting(Operation):
3190 3203 Operation.__init__(self)
3191 3204 self.i=0
3192 3205 self.isConfig = False
3206 self.aux = 1
3193 3207
3194 3208 def setup(self,nChan,nProf,nHei,nBlocks):
3195 3209 self.__dataReady = False
@@ -3212,6 +3226,7 class SpectralFitting(Operation):
3212 3226 if (wwauto is None): wwauto = 0
3213 3227
3214 3228 if (n0 < 1.e-20): n0 = 1.e-20
3229
3215 3230 freq = oldfreq
3216 3231 vec_power = numpy.zeros(oldspec.shape[1])
3217 3232 vec_fd = numpy.zeros(oldspec.shape[1])
@@ -3620,7 +3635,7 class SpectralFitting(Operation):
3620 3635 junk = tmp
3621 3636 tmp = junk*0.0
3622 3637
3623 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
3638 tmp[list(xpos + (ypos*num_prof))] = junk[list(xpos + (ypos*num_prof))]
3624 3639 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
3625 3640
3626 3641 return array
@@ -3642,6 +3657,146 class SpectralFitting(Operation):
3642 3657 fmom = numpy.sum(doppler*ytemp)/numpy.sum(ytemp)+(index-(npoints/2-1))*numpy.abs(doppler[1]-doppler[0])
3643 3658 smom = numpy.sum(doppler*doppler*ytemp)/numpy.sum(ytemp)
3644 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 3801 def windowing_single_old(self,spc,x,A,B,C,D,nFFTPoints):
3647 3802 '''
@@ -4010,7 +4165,6 class SpectralFitting(Operation):
4010 4165 power = numpy.sum(spectra, axis=1)
4011 4166 nPairs = len(crosspairs)
4012 4167 absc = dataOut.abscissaList[:-1]
4013 print('para escribir h5 ',dataOut.paramInterval)
4014 4168 if not self.isConfig:
4015 4169 self.isConfig = True
4016 4170
@@ -4024,6 +4178,7 class SpectralFitting(Operation):
4024 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 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 4182 dataOut.data_spc = incoh_spectra
4028 4183 dataOut.data_cspc = incoh_cspectra
4029 4184 dataOut.sat_spectra = sat_spectra
@@ -4138,7 +4293,7 class SpectralFitting(Operation):
4138 4293 dataOut.nIncohInt= int(dataOut.nIncohInt)
4139 4294 dataOut.nProfiles = int(dataOut.nProfiles)
4140 4295 dataOut.nCohInt = int(dataOut.nCohInt)
4141 print('sale',dataOut.nProfiles,dataOut.nHeights)
4296 #print('sale',dataOut.nProfiles,dataOut.nHeights)
4142 4297 #dataOut.nFFTPoints=nprofs
4143 4298 #dataOut.normFactor = nprofs
4144 4299 dataOut.channelList = channelList
@@ -4148,7 +4303,7 class SpectralFitting(Operation):
4148 4303 vmax = (300000000/49920000.0/2) / (dataOut.ippSeconds)
4149 4304 #dataOut.ippSeconds=ipp
4150 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 4307 print('Empieza procesamiento offline')
4153 4308 if path != None:
4154 4309 sys.path.append(path)
@@ -4170,6 +4325,8 class SpectralFitting(Operation):
4170 4325 dataOut.data_snr1_i = numpy.zeros((nGroups*2, nHeights))*numpy.nan
4171 4326 # dataOut.smooth_i = numpy.zeros((nGroups*2, nHeights))*numpy.nan
4172 4327
4328 if self.aux == 1:
4329 dataOut.p03last = numpy.zeros(nGroups, 'float32')
4173 4330 for i in range(nGroups):
4174 4331 coord = groupArray[i,:]
4175 4332 #Input data array
@@ -4193,7 +4350,7 class SpectralFitting(Operation):
4193 4350 n1= n1i
4194 4351 my_noises[2*i+0] = n0
4195 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 4354 snrth = 10**(snrth/10.0)
4198 4355 jvelr = numpy.zeros(nHeights, dtype = 'float')
4199 4356 #snr0 = numpy.zeros(nHeights, dtype = 'float')
@@ -4201,7 +4358,9 class SpectralFitting(Operation):
4201 4358 hvalid = [0]
4202 4359
4203 4360 coh2 = abs(dataOut.data_cspc[i,1:nProf,:])**2/(dataOut.data_spc[0+i*2,1:nProf-0,:]*dataOut.data_spc[1+i*2,1:nProf-0,:])
4204
4361
4362 SIGNAL_0 = []
4363 SIGNAL_1 = []
4205 4364 for h in range(nHeights):
4206 4365 smooth = clean_num_aver[i+1,h]
4207 4366 signalpn0 = (dataOut.data_spc[i*2,1:(nProf-0),h])/smooth
@@ -4228,10 +4387,32 class SpectralFitting(Operation):
4228 4387 else: jvelr[h] = absc[0]
4229 4388 if snr0 > 0.1 and snr1 > 0.1: hvalid = numpy.concatenate((hvalid,h), axis=None)
4230 4389 #print(maxp0,absc[maxp0],snr0,jvelr[h])
4390 SIGNAL_0.append(signal0)
4391 SIGNAL_1.append(signal1)
4231 4392
4232 4393 if len(hvalid)> 1: fd0 = numpy.median(jvelr[hvalid[1:]])*-1
4233 else: fd0 = numpy.nan
4394 else: fd0 = numpy.nan
4234 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 4416 for h in range(nHeights):
4236 4417 d = data[:,h]
4237 4418 smooth = clean_num_aver[i+1,h] #dataOut.data_spc[:,1:nProf-0,:]
@@ -4294,6 +4475,16 class SpectralFitting(Operation):
4294 4475 minp = p0*numpy.nan
4295 4476 error0 = numpy.nan
4296 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 4488 else :
4298 4489 data_spc = dataOut.data_spc[coord,:,h]
4299 4490 p0 = numpy.array(self.library.initialValuesFunction(data_spc, constants))
@@ -4308,7 +4499,25 class SpectralFitting(Operation):
4308 4499 dataOut.data_param[i,:,h] = minp
4309 4500 dataOut.data_snr1_i[i*2,h] = numpy.sum(signalpn0/(nProf-1))/n0
4310 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 4521 for ht in range(nHeights-1) :
4313 4522 smooth = coh_num_aver[i+1,ht] #datc[0,ht,0,beam]
4314 4523 dataOut.data_paramC[4*i,ht,1] = smooth
@@ -4393,6 +4602,7 class SpectralFitting(Operation):
4393 4602 sat_fits[1+2*beam,ht] = numpy.sum(signalpn1/(val1_npoints*nProf))/n1
4394 4603
4395 4604 dataOut.sat_fits = sat_fits
4605 self.aux = 0
4396 4606 return dataOut
4397 4607
4398 4608 def __residFunction(self, p, dp, LT, constants):
@@ -4588,7 +4798,6 class WindProfiler(Operation):
4588 4798 return distx, disty, dist, ang
4589 4799 #Calculo de Matrices
4590 4800
4591
4592 4801 def __calculateVelVer(self, phase, lagTRange, _lambda):
4593 4802
4594 4803 Ts = lagTRange[1] - lagTRange[0]
@@ -4715,8 +4924,6 class WindProfiler(Operation):
4715 4924 azim = meteorAux[:, 3]*numpy.pi/180
4716 4925
4717 4926 n = numpy.cos(zen)
4718 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
4719 # l = m*numpy.tan(azim)
4720 4927 l = numpy.sin(zen)*numpy.sin(azim)
4721 4928 m = numpy.sin(zen)*numpy.cos(azim)
4722 4929
@@ -4781,7 +4988,6 class WindProfiler(Operation):
4781 4988
4782 4989 for t in uniqueTime:
4783 4990 metArray1 = metArray[utctime==t,:]
4784 # phaseDerThresh = numpy.pi/4 #reducir Phase thresh
4785 4991 tmet = metArray1[:,1].astype(int)
4786 4992 hmet = metArray1[:,2].astype(int)
4787 4993
@@ -4912,8 +5118,7 class WindProfiler(Operation):
4912 5118
4913 5119 def run(self, dataOut, technique, nHours=1, hmin=70, hmax=110, **kwargs):
4914 5120
4915 param = dataOut.data_param
4916 #if dataOut.abscissaList != None:
5121 param = dataOut.moments ### LA VERSION ANTERIOR trabajaba con param = dataOut.data_param
4917 5122 if numpy.any(dataOut.abscissaList):
4918 5123 absc = dataOut.abscissaList[:-1]
4919 5124 # noise = dataOut.noise
@@ -4933,29 +5138,9 class WindProfiler(Operation):
4933 5138 elif technique == 'SA':
4934 5139
4935 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 5141 kwargs['groupList'] = dataOut.groupList
4956 5142 kwargs['tau'] = dataOut.data_param
4957 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 5144 dataOut.data_output = self.techniqueSA(kwargs)
4960 5145 dataOut.utctimeInit = dataOut.utctime
4961 5146 dataOut.outputInterval = dataOut.timeInterval
@@ -4984,7 +5169,6 class WindProfiler(Operation):
4984 5169 dataOut.outputInterval = nHours*3600
4985 5170
4986 5171 if self.__isConfig == False:
4987 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
4988 5172 #Get Initial LTC time
4989 5173 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
4990 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 5231 if self.__isConfig == False:
5048 5232 dataOut.outputInterval = nMins*60
5049 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
5050 5233 #Get Initial LTC time
5051 5234 initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
5052 5235 minuteAux = initime.minute
@@ -5077,69 +5260,19 class WindProfiler(Operation):
5077 5260 dataOut.flagNoData = False
5078 5261 self.__buffer = None
5079 5262
5080 return
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
5263 return dataOut
5095 5264
5096 n = None
5265 class EWDriftsEstimation(Operation):
5097 5266
5098 5267 def __init__(self):
5099 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-cosDirX**2)
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_x**2-theta_y**2)
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 5270 def __correctValues(self, heiRang, phi, velRadial, SNR):
5137 5271 listPhi = phi.tolist()
5138 5272 maxid = listPhi.index(max(listPhi))
5139 5273 minid = listPhi.index(min(listPhi))
5140 5274
5141 5275 rango = list(range(len(phi)))
5142
5143 5276 heiRang1 = heiRang*math.cos(phi[maxid])
5144 5277 heiRangAux = heiRang*math.cos(phi[minid])
5145 5278 indOut = (heiRang1 < heiRangAux[0]).nonzero()
@@ -5151,13 +5284,18 class WindProfiler(Operation):
5151 5284 for i in rango:
5152 5285 x = heiRang*math.cos(phi[i])
5153 5286 y1 = velRadial[i,:]
5154 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
5155
5287 vali= (numpy.isfinite(y1)==True).nonzero()
5288 y1=y1[vali]
5289 x = x[vali]
5290 f1 = interpolate.interp1d(x,y1,kind = 'cubic',bounds_error=False)
5156 5291 x1 = heiRang1
5157 5292 y11 = f1(x1)
5158
5159 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 5299 y21 = f2(x1)
5162 5300
5163 5301 velRadial1[i,:] = y11
@@ -5165,713 +5303,124 class WindProfiler(Operation):
5165 5303
5166 5304 return heiRang1, velRadial1, SNR1
5167 5305
5168 def __calculateVelUVW(self, A, velRadial):
5169
5170 #Operacion Matricial
5171 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
5172 velUVW[:,:] = numpy.dot(A,velRadial)
5173
5174
5175 return velUVW
5306 def run(self, dataOut, zenith, zenithCorrection,fileDrifts,beam_pos=None):
5307 dataOut.lat = -11.95
5308 dataOut.lon = -76.87
5309 dataOut.spcst = 0.00666
5310 dataOut.pl = 0.0003
5311 dataOut.cbadn = 3
5312 dataOut.inttms = 300
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):
5178 """
5179 Function that implements Doppler Beam Swinging (DBS) technique.
5323 heiRang = dataOut.heightList
5324 velRadial = dataOut.data_param[:,3,:]
5325 velRadialm = dataOut.data_param[:,2:4,:]*-1
5180 5326
5181 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
5182 Direction correction (if necessary), Ranges and SNR
5327 rbufc=dataOut.data_paramC[:,:,0]
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_w0*smooth_w+p_w0C*smooth_wC)/(smooth_w+smooth_wC)
5371 p_w1_a = (p_w1*smooth_w+p_w1C*smooth_wC)/(smooth_w+smooth_wC)
5185 5372
5186 Parameters affected: Winds, height range, SNR
5187 """
5188 velRadial0 = kwargs['velRadial']
5189 heiRang = kwargs['heightList']
5190 SNR0 = kwargs['SNR']
5373 if len(sat_fits) >0 :
5374 p_w0C = p_w0C + sat_fits[0,:]
5375 p_w1C = p_w1C + sat_fits[1,:]
5376
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:
5193 theta_x = numpy.array(kwargs['dirCosx'])
5194 theta_y = numpy.array(kwargs['dirCosy'])
5195 else:
5196 elev = numpy.array(kwargs['elevation'])
5197 azim = numpy.array(kwargs['azimuth'])
5198 theta_x, theta_y = self.__calculateCosDir(elev, azim)
5199 azimuth = kwargs['correctAzimuth']
5200 if 'horizontalOnly' in kwargs:
5201 horizontalOnly = kwargs['horizontalOnly']
5202 else: horizontalOnly = False
5203 if 'correctFactor' in kwargs:
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]
5390 zenith = numpy.array(zenith)
5391 zenith -= zenithCorrection
5392 zenith *= numpy.pi/180
5393 if zenithCorrection != 0 :
5394 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
5395 else :
5396 heiRang1 = heiRang
5397 velRadial1 = velRadial
5398 SNR1 = SNR
5399
5400 alp = zenith[0]
5401 bet = zenith[1]
5214 5402
5215 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
5216 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
5217 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
5403 w_w = velRadial1[0,:]
5404 w_e = velRadial1[1,:]
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
5220 winds = self.__calculateVelUVW(A,velRadial1)
5412 tini=time.localtime(dataOut.utctime)
5221 5413
5222 return winds, heiRang1, SNR1
5223
5224 def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
5225
5226 nPairs = len(pairs_ccf)
5227 posx = numpy.asarray(posx)
5228 posy = numpy.asarray(posy)
5229
5230 #Rotacion Inversa para alinear con el azimuth
5231 if azimuth!= None:
5232 azimuth = azimuth*math.pi/180
5233 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
5234 posy1 = -posx*math.sin(azimuth) + posy*math.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 = -_lambda*phase/(4*math.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 = dist1*tau1/(2*tau2**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#*height*1000/(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(diffX**2 + diffY**2)
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_w0*smooth_w+p_w0C*smooth_wC)/(smooth_w+smooth_wC)
5824 #p_w1 = (p_w1*smooth_w+p_w1C*smooth_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
5414 if tini[3] >= 6 and tini[3] < 18 :
5415 w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w)
5416 w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err)
5417 else:
5418 w_wtmp = numpy.where(numpy.isfinite(w_wC)==True,w_wC,w_w)
5419 w_wtmp = numpy.where(range1 > 200,w_w,w_wtmp)
5420 w_w_errtmp = numpy.where(numpy.isfinite(w_wC)==True,numpy.nan,w_w_err)
5421 w_w_errtmp = numpy.where(range1 > 200,w_w_err,w_w_errtmp)
5422 w_w = w_wtmp
5423 w_w_err = w_w_errtmp
5875 5424
5876 5425 #if my_nbeams == 2:
5877 5426 smooth_eC=ebufc[4,:]
@@ -5891,9 +5440,9 class EWDriftsEstimation(Operation):
5891 5440 smooth_e[val]=0
5892 5441 val = (numpy.isfinite(smooth_eC)==False).nonzero()
5893 5442 smooth_eC[val]=0
5894 #p_e0 = (p_e0*smooth_e+p_e0C*smooth_eC)/(smooth_e+smooth_eC)
5895 #p_e1 = (p_e1*smooth_e+p_e1C*smooth_eC)/(smooth_e+smooth_eC)
5896
5443 p_e0_a = (p_e0*smooth_e+p_e0C*smooth_eC)/(smooth_e+smooth_eC)
5444 p_e1_a = (p_e1*smooth_e+p_e1C*smooth_eC)/(smooth_e+smooth_eC)
5445 #print(w_e,w_eC,p_e0C,sat_fits[2,:])
5897 5446 if len(sat_fits) >0 :
5898 5447 p_e0C = p_e0C + sat_fits[2,:]
5899 5448 p_e1C = p_e1C + sat_fits[3,:]
@@ -5914,6 +5463,22 class EWDriftsEstimation(Operation):
5914 5463
5915 5464 w_err = numpy.sqrt((w_w_err*numpy.sin(bet))**2.+(w_e_err*numpy.sin(alp))**2.)/ numpy.absolute(numpy.cos(alp)*numpy.sin(bet)-numpy.cos(bet)*numpy.sin(alp))
5916 5465 u_err = numpy.sqrt((w_w_err*numpy.cos(bet))**2.+(w_e_err*numpy.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 5483 winds = numpy.vstack((w,u))
5919 5484 dataOut.heightList = heiRang1
@@ -5924,7 +5489,7 class EWDriftsEstimation(Operation):
5924 5489 #print('alt ',range1*numpy.sin(86.1*numpy.pi/180))
5925 5490 #print(numpy.min([dataOut.eldir7,dataOut.eldir8]))
5926 5491 galt = range1*numpy.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 5493 #snr1 = 10*numpy.log10(SNR1[0])
5929 5494 #print(min(snr1), max(snr1))
5930 5495 snr1 = numpy.vstack((p_w0,p_w1,p_e0,p_e1))
@@ -6001,39 +5566,43 class EWDriftsEstimation(Operation):
6001 5566 uA_err[nhei_avg-1] = sigma
6002 5567 uA = uA[6*h_avgs:nhei_avg-0]
6003 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 5571 if my_nbeams == 1: dataOut.drifts_avg = wA
6007 5572 #deltahavg= wA*0.0+deltah
6008 5573 dataOut.range = range1
6009 5574 galtavg = range_aver*numpy.sin(numpy.min([dataOut.elw,dataOut.ele])*numpy.pi/180.)
6010 5575 dataOut.params_avg = numpy.vstack((wA,wA_err,uA,uA_err,range_aver,galtavg,delta_h))
6011
6012 #print('comparando dim de avg ',wA.shape,deltahavg.shape,range_aver.shape)
5576
5577 '''
6013 5578 tini=time.localtime(dataOut.utctime)
6014 5579 datefile= str(tini[0]).zfill(4)+str(tini[1]).zfill(2)+str(tini[2]).zfill(2)
6015 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 5583 f1 = open(nfile,'a')
6018 5584 datedriftavg=str(tini[0])+' '+str(tini[1])+' '+str(tini[2])+' '+str(tini[3])+' '+str(tini[4])
6019 5585 driftavgstr=str(dataOut.drifts_avg)
6020 5586 numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')
6021 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 5589 numpy.savetxt(f1,dataOut.drifts_avg[:,:],fmt='%10.2f')
6023 5590 f1.close()
6024
5591 '''
5592 '''
6025 5593 swfile = fileDrifts+'/jro'+datefile+'drifts_sw.txt'
6026 5594 f1 = open(swfile,'a')
6027 5595 numpy.savetxt(f1,numpy.column_stack([tini[0],tini[1],tini[2],tini[3],tini[4]]),fmt='%4i')
6028 5596 numpy.savetxt(f1,numpy.reshape(heiRang,(1,len(heiRang))),fmt='%10.2f')
6029 5597 numpy.savetxt(f1,dataOut.data_param[:,0,:],fmt='%10.2f')
6030 5598 f1.close()
6031 dataOut.heightListtmp = dataOut.heightList
6032 5599 '''
5600 dataOut.heightListtmp = dataOut.heightList
5601
6033 5602 #Envio data de drifts a mysql
6034 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 5604 mydb = mysql.connector.connect(
6036 host="10.10.110.213",
5605 host="10.10.110.205",
6037 5606 user="user_clima",
6038 5607 password="5D.bh(B2)Y_wRNz9",
6039 5608 database="clima_espacial"
@@ -6056,7 +5625,7 class EWDriftsEstimation(Operation):
6056 5625 mydb.commit()
6057 5626
6058 5627 print(mycursor.rowcount, "record inserted.")
6059 '''
5628
6060 5629 return dataOut
6061 5630
6062 5631 class setHeightDrifts(Operation):
@@ -7317,7 +6886,6 class SMOperations():
7317 6886 chan1Y = chan1
7318 6887 chan2Y = chan2
7319 6888 distances[2:4] = numpy.array([d1,d2])
7320
7321 6889 pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
7322 6890
7323 6891 return pairslist, distances
General Comments 0
You need to be logged in to leave comments. Login now