The requested changes are too big and content was truncated. Show full diff
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 |
1 | NO CONTENT: modified file |
|
NO CONTENT: modified file | ||
The requested commit or file is too big and content was truncated. Show full diff |
General Comments 0
You need to be logged in to leave comments.
Login now