##// END OF EJS Templates
27 de Marzo 2018
ebocanegra -
r1154:a08206fb2eae
parent child
Show More
@@ -1,4022 +1,4025
1 1 import numpy
2 2 import math
3 3 from scipy import optimize, interpolate, signal, stats, ndimage
4 4 import scipy
5 5 import re
6 6 import datetime
7 7 import copy
8 8 import sys
9 9 import importlib
10 10 import itertools
11 11 from multiprocessing import Pool, TimeoutError
12 12 from multiprocessing.pool import ThreadPool
13 import copy_reg
14 import cPickle
15 import types
16 from functools import partial
13
14
15
17 16 import time
18 17 #from sklearn.cluster import KMeans
19 18
20 19 import matplotlib.pyplot as plt
21 20
22 21 from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
23 22 from jroproc_base import ProcessingUnit, Operation
24 23 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
25 24 from scipy import asarray as ar,exp
26 25 from scipy.optimize import curve_fit
27 26
28 27 import warnings
29 28 from numpy import NaN
30 29 from scipy.optimize.optimize import OptimizeWarning
31 30 warnings.filterwarnings('ignore')
32 31
33 32
34 33 SPEED_OF_LIGHT = 299792458
35 34
36 35
37 36 '''solving pickling issue'''
38 37
39 38 def _pickle_method(method):
40 39 func_name = method.im_func.__name__
41 40 obj = method.im_self
42 41 cls = method.im_class
43 42 return _unpickle_method, (func_name, obj, cls)
44 43
45 44 def _unpickle_method(func_name, obj, cls):
46 45 for cls in cls.mro():
47 46 try:
48 47 func = cls.__dict__[func_name]
49 48 except KeyError:
50 49 pass
51 50 else:
52 51 break
53 52 return func.__get__(obj, cls)
54 53
55 54 class ParametersProc(ProcessingUnit):
56 55
57 56 nSeconds = None
58 57
59 58 def __init__(self):
60 59 ProcessingUnit.__init__(self)
61 60
62 61 # self.objectDict = {}
63 62 self.buffer = None
64 63 self.firstdatatime = None
65 64 self.profIndex = 0
66 65 self.dataOut = Parameters()
67 66
68 67 def __updateObjFromInput(self):
69 68
70 69 self.dataOut.inputUnit = self.dataIn.type
71 70
72 71 self.dataOut.timeZone = self.dataIn.timeZone
73 72 self.dataOut.dstFlag = self.dataIn.dstFlag
74 73 self.dataOut.errorCount = self.dataIn.errorCount
75 74 self.dataOut.useLocalTime = self.dataIn.useLocalTime
76 75
77 76 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
78 77 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
79 78 self.dataOut.channelList = self.dataIn.channelList
80 79 self.dataOut.heightList = self.dataIn.heightList
81 80 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
82 81 # self.dataOut.nHeights = self.dataIn.nHeights
83 82 # self.dataOut.nChannels = self.dataIn.nChannels
84 83 self.dataOut.nBaud = self.dataIn.nBaud
85 84 self.dataOut.nCode = self.dataIn.nCode
86 85 self.dataOut.code = self.dataIn.code
87 86 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
88 87 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
89 88 # self.dataOut.utctime = self.firstdatatime
90 89 self.dataOut.utctime = self.dataIn.utctime
91 90 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
92 91 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
93 92 self.dataOut.nCohInt = self.dataIn.nCohInt
94 93 # self.dataOut.nIncohInt = 1
95 94 self.dataOut.ippSeconds = self.dataIn.ippSeconds
96 95 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
97 96 self.dataOut.timeInterval1 = self.dataIn.timeInterval
98 97 self.dataOut.heightList = self.dataIn.getHeiRange()
99 98 self.dataOut.frequency = self.dataIn.frequency
100 99 self.dataOut.noise = self.dataIn.noise
101 100
102 101
103 102
104 103 def run(self):
105 104
106 105 #---------------------- Voltage Data ---------------------------
107 106
108 107 if self.dataIn.type == "Voltage":
109 108
110 109 self.__updateObjFromInput()
111 110 self.dataOut.data_pre = self.dataIn.data.copy()
112 111 self.dataOut.flagNoData = False
113 112 self.dataOut.utctimeInit = self.dataIn.utctime
114 113 self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds
115 114 return
116 115
117 116 #---------------------- Spectra Data ---------------------------
118 117
119 118 if self.dataIn.type == "Spectra":
120 119
121 120 self.dataOut.data_pre = (self.dataIn.data_spc , self.dataIn.data_cspc)
122 121 print 'self.dataIn.data_spc', self.dataIn.data_spc.shape
123 122 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
124 123 self.dataOut.spc_noise = self.dataIn.getNoise()
125 124 self.dataOut.spc_range = numpy.asanyarray((self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1) ))
126 125
127 126 self.dataOut.normFactor = self.dataIn.normFactor
128 127 #self.dataOut.outputInterval = self.dataIn.outputInterval
129 128 self.dataOut.groupList = self.dataIn.pairsList
130 129 self.dataOut.flagNoData = False
131 130 #print 'datain chandist ',self.dataIn.ChanDist
132 131 if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
133 132 self.dataOut.ChanDist = self.dataIn.ChanDist
134 133 else: self.dataOut.ChanDist = None
135 134
136 135 print 'datain chandist ',self.dataOut.ChanDist
137 136
138 137 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
139 138 # self.dataOut.VelRange = self.dataIn.VelRange
140 139 #else: self.dataOut.VelRange = None
141 140
142 141 if hasattr(self.dataIn, 'RadarConst'): #Radar Constant
143 142 self.dataOut.RadarConst = self.dataIn.RadarConst
144 143
145 144 if hasattr(self.dataIn, 'NPW'): #NPW
146 145 self.dataOut.NPW = self.dataIn.NPW
147 146
148 147 if hasattr(self.dataIn, 'COFA'): #COFA
149 148 self.dataOut.COFA = self.dataIn.COFA
150 149
151 150
152 151
153 152 #---------------------- Correlation Data ---------------------------
154 153
155 154 if self.dataIn.type == "Correlation":
156 155 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.dataIn.splitFunctions()
157 156
158 157 self.dataOut.data_pre = (self.dataIn.data_cf[acf_ind,:], self.dataIn.data_cf[ccf_ind,:,:])
159 158 self.dataOut.normFactor = (self.dataIn.normFactor[acf_ind,:], self.dataIn.normFactor[ccf_ind,:])
160 159 self.dataOut.groupList = (acf_pairs, ccf_pairs)
161 160
162 161 self.dataOut.abscissaList = self.dataIn.lagRange
163 162 self.dataOut.noise = self.dataIn.noise
164 163 self.dataOut.data_SNR = self.dataIn.SNR
165 164 self.dataOut.flagNoData = False
166 165 self.dataOut.nAvg = self.dataIn.nAvg
167 166
168 167 #---------------------- Parameters Data ---------------------------
169 168
170 169 if self.dataIn.type == "Parameters":
171 170 self.dataOut.copy(self.dataIn)
172 171 self.dataOut.flagNoData = False
173 172
174 173 return True
175 174
176 175 self.__updateObjFromInput()
177 176 self.dataOut.utctimeInit = self.dataIn.utctime
178 177 self.dataOut.paramInterval = self.dataIn.timeInterval
179 178
180 179 return
181 180
182 181
183 182 def target(tups):
184 183
185 184 obj, args = tups
186 185 #print 'TARGETTT', obj, args
187 186 return obj.FitGau(args)
188 187
189 188
190 189 class SpectralFilters(Operation):
191 190
192 191 '''This class allows the Rainfall / Wind Selection for CLAIRE RADAR
193 192
194 193 LimitR : It is the limit in m/s of Rainfall
195 194 LimitW : It is the limit in m/s for Winds
196 195
197 196 Input:
198 197
199 198 self.dataOut.data_pre : SPC and CSPC
200 199 self.dataOut.spc_range : To select wind and rainfall velocities
201 200
202 201 Affected:
203 202
204 203 self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
205 204 self.dataOut.spcparam_range : Used in SpcParamPlot
206 205 self.dataOut.SPCparam : Used in PrecipitationProc
207 206
208 207
209 208 '''
210 209
211 210 def __init__(self, **kwargs):
212 211 Operation.__init__(self, **kwargs)
213 212 self.i=0
214 213
215 214 def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5):
216 215
217 216
218 217 #Limite de vientos
219 218 LimitR = PositiveLimit
220 219 LimitN = NegativeLimit
221 220
222 221 self.spc = dataOut.data_pre[0].copy()
223 222 self.cspc = dataOut.data_pre[1].copy()
224 223
225 224 self.Num_Hei = self.spc.shape[2]
226 225 self.Num_Bin = self.spc.shape[1]
227 226 self.Num_Chn = self.spc.shape[0]
228 227
229 228 VelRange = dataOut.spc_range[2]
230 229 TimeRange = dataOut.spc_range[1]
231 230 FrecRange = dataOut.spc_range[0]
232 231
233 232 Vmax= 2*numpy.max(dataOut.spc_range[2])
234 233 Tmax= 2*numpy.max(dataOut.spc_range[1])
235 234 Fmax= 2*numpy.max(dataOut.spc_range[0])
236 235
237 236 Breaker1R=VelRange[numpy.abs(VelRange-(-LimitN)).argmin()]
238 237 Breaker1R=numpy.where(VelRange == Breaker1R)
239 238
240 239 Delta = self.Num_Bin/2 - Breaker1R[0]
241 240
242 241 #Breaker1W=VelRange[numpy.abs(VelRange-(-LimitW)).argmin()]
243 242 #Breaker1W=numpy.where(VelRange == Breaker1W)
244 243
245 244 #Breaker2W=VelRange[numpy.abs(VelRange-(LimitW)).argmin()]
246 245 #Breaker2W=numpy.where(VelRange == Breaker2W)
247 246
248 247
249 248 '''Reacomodando SPCrange'''
250 249
251 250 VelRange=numpy.roll(VelRange,-(self.Num_Bin/2) ,axis=0)
252 251
253 252 VelRange[-(self.Num_Bin/2):]+= Vmax
254 253
255 254 FrecRange=numpy.roll(FrecRange,-(self.Num_Bin/2),axis=0)
256 255
257 256 FrecRange[-(self.Num_Bin/2):]+= Fmax
258 257
259 258 TimeRange=numpy.roll(TimeRange,-(self.Num_Bin/2),axis=0)
260 259
261 260 TimeRange[-(self.Num_Bin/2):]+= Tmax
262 261
263 262 ''' ------------------ '''
264 263
265 264 Breaker2R=VelRange[numpy.abs(VelRange-(LimitR)).argmin()]
266 265 Breaker2R=numpy.where(VelRange == Breaker2R)
267 266
268 267
269 268 SPCroll = numpy.roll(self.spc,-(self.Num_Bin/2) ,axis=1)
270 269
271 270 SPCcut = SPCroll.copy()
272 271 for i in range(self.Num_Chn):
273 272
274 273 SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i]
275 274 SPCcut[i,-int(Delta):,:] = dataOut.noise[i]
276 275
277 276 SPCcut[i]=SPCcut[i]- dataOut.noise[i]
278 277 SPCcut[ numpy.where( SPCcut<0 ) ] = 1e-20
279 278
280 279 SPCroll[i]=SPCroll[i]-dataOut.noise[i]
281 280 SPCroll[ numpy.where( SPCroll<0 ) ] = 1e-20
282 281
283 282 #self.spc[i, 0:int(Breaker1W[0]) ,:] = dataOut.noise[i]
284 283 #self.spc[i, int(Breaker2W[0]):self.Num_Bin ,:] = dataOut.noise[i]
285 284
286 285 #self.cspc[i, 0:int(Breaker1W[0]) ,:] = dataOut.noise[i]
287 286 #self.cspc[i, int(Breaker2W[0]):self.Num_Bin ,:] = dataOut.noise[i]
288 287
289 288
290 289
291 290
292 291
293 292 SPC_ch1 = SPCroll
294 293
295 294 SPC_ch2 = SPCcut
296 295
297 296 SPCparam = (SPC_ch1, SPC_ch2, self.spc)
298 297 dataOut.SPCparam = numpy.asarray(SPCparam)
299 298
300 299 #dataOut.data_pre= (self.spc , self.cspc)
301 300
302 301 #dataOut.data_preParam = (self.spc , self.cspc)
303 302
304 303 dataOut.spcparam_range=numpy.zeros([self.Num_Chn,self.Num_Bin+1])
305 304
306 305 dataOut.spcparam_range[2]=VelRange
307 306 dataOut.spcparam_range[1]=TimeRange
308 307 dataOut.spcparam_range[0]=FrecRange
309 308
310 309
311 310
312 311
313 312 class GaussianFit(Operation):
314 313
315 314 '''
316 315 Function that fit of one and two generalized gaussians (gg) based
317 316 on the PSD shape across an "power band" identified from a cumsum of
318 317 the measured spectrum - noise.
319 318
320 319 Input:
321 320 self.dataOut.data_pre : SelfSpectra
322 321
323 322 Output:
324 323 self.dataOut.SPCparam : SPC_ch1, SPC_ch2
325 324
326 325 '''
327 326 def __init__(self, **kwargs):
328 327 Operation.__init__(self, **kwargs)
329 328 self.i=0
330 329
331 330
332 331 def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
333 332 """This routine will find a couple of generalized Gaussians to a power spectrum
334 333 input: spc
335 334 output:
336 335 Amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1,noise
337 336 """
338 337
339 338 self.spc = dataOut.data_pre[0].copy()
340 339
341 340
342 341 print 'SelfSpectra Shape', numpy.asarray(self.spc).shape
343 342
344 343
345 344 #plt.figure(50)
346 345 #plt.subplot(121)
347 346 #plt.plot(self.spc,'k',label='spc(66)')
348 347 #plt.plot(xFrec,ySamples[1],'g',label='Ch1')
349 348 #plt.plot(xFrec,ySamples[2],'r',label='Ch2')
350 349 #plt.plot(xFrec,FitGauss,'yo:',label='fit')
351 350 #plt.legend()
352 351 #plt.title('DATOS A ALTURA DE 7500 METROS')
353 352 #plt.show()
354 353
355 354 self.Num_Hei = self.spc.shape[2]
356 355 #self.Num_Bin = len(self.spc)
357 356 self.Num_Bin = self.spc.shape[1]
358 357 self.Num_Chn = self.spc.shape[0]
359 358 Vrange = dataOut.abscissaList
360 359
361 360 GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei])
362 361 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
363 362 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
364 363 SPC_ch1[:] = numpy.NaN
365 364 SPC_ch2[:] = numpy.NaN
366 365
367 366
368 367 start_time = time.time()
369 368
370 369 noise_ = dataOut.spc_noise[0].copy()
371 370
372 371
373 372 pool = Pool(processes=self.Num_Chn)
374 373 args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
375 374 objs = [self for __ in range(self.Num_Chn)]
376 375 attrs = zip(objs, args)
377 376 gauSPC = pool.map(target, attrs)
378 377 dataOut.SPCparam = numpy.asarray(SPCparam)
379 378
380 379
381 380
382 381 print '========================================================'
383 382 print 'total_time: ', time.time()-start_time
384 383
385 384 # re-normalizing spc and noise
386 385 # This part differs from gg1
387 386
388 387
389 388
390 389 ''' Parameters:
391 390 1. Amplitude
392 391 2. Shift
393 392 3. Width
394 393 4. Power
395 394 '''
396 395
397 396
398 397 ###############################################################################
399 398 def FitGau(self, X):
400 399
401 400 Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
402 401 #print 'VARSSSS', ch, pnoise, noise, num_intg
403 402
404 403 #print 'HEIGHTS', self.Num_Hei
405 404
406 405 SPCparam = []
407 406 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
408 407 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
409 408 SPC_ch1[:] = 0#numpy.NaN
410 409 SPC_ch2[:] = 0#numpy.NaN
411 410
412 411
413 412
414 413 for ht in range(self.Num_Hei):
415 414 #print (numpy.asarray(self.spc).shape)
416 415
417 416 #print 'TTTTT', ch , ht
418 417 #print self.spc.shape
419 418
420 419
421 420 spc = numpy.asarray(self.spc)[ch,:,ht]
422 421
423 422 #############################################
424 423 # normalizing spc and noise
425 424 # This part differs from gg1
426 425 spc_norm_max = max(spc)
427 426 #spc = spc / spc_norm_max
428 427 pnoise = pnoise #/ spc_norm_max
429 428 #############################################
430 429
431 430 fatspectra=1.0
432 431
433 432 wnoise = noise_ #/ spc_norm_max
434 433 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
435 434 #if wnoise>1.1*pnoise: # to be tested later
436 435 # wnoise=pnoise
437 436 noisebl=wnoise*0.9;
438 437 noisebh=wnoise*1.1
439 438 spc=spc-wnoise
440 439 # print 'wnoise', noise_[0], spc_norm_max, wnoise
441 440 minx=numpy.argmin(spc)
442 441 #spcs=spc.copy()
443 442 spcs=numpy.roll(spc,-minx)
444 443 cum=numpy.cumsum(spcs)
445 444 tot_noise=wnoise * self.Num_Bin #64;
446 445 #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
447 446 #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
448 447 #snr=tot_signal/tot_noise
449 448 #snr=cum[-1]/tot_noise
450 449 snr = sum(spcs)/tot_noise
451 450 snrdB=10.*numpy.log10(snr)
452 451
453 452 if snrdB < SNRlimit :
454 453 snr = numpy.NaN
455 454 SPC_ch1[:,ht] = 0#numpy.NaN
456 455 SPC_ch1[:,ht] = 0#numpy.NaN
457 456 SPCparam = (SPC_ch1,SPC_ch2)
458 457 continue
459 458 #print 'snr',snrdB #, sum(spcs) , tot_noise
460 459
461 460
462 461
463 462 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
464 463 # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
465 464
466 465 cummax=max(cum);
467 466 epsi=0.08*fatspectra # cumsum to narrow down the energy region
468 467 cumlo=cummax*epsi;
469 468 cumhi=cummax*(1-epsi)
470 469 powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
471 470
472 471
473 472 if len(powerindex) < 1:# case for powerindex 0
474 473 continue
475 474 powerlo=powerindex[0]
476 475 powerhi=powerindex[-1]
477 476 powerwidth=powerhi-powerlo
478 477
479 478 firstpeak=powerlo+powerwidth/10.# first gaussian energy location
480 479 secondpeak=powerhi-powerwidth/10.#second gaussian energy location
481 480 midpeak=(firstpeak+secondpeak)/2.
482 481 firstamp=spcs[int(firstpeak)]
483 482 secondamp=spcs[int(secondpeak)]
484 483 midamp=spcs[int(midpeak)]
485 484
486 485 x=numpy.arange( self.Num_Bin )
487 486 y_data=spc+wnoise
488 487
489 488 ''' single Gaussian '''
490 489 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
491 490 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
492 491 power0=2.
493 492 amplitude0=midamp
494 493 state0=[shift0,width0,amplitude0,power0,wnoise]
495 494 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
496 495 lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
497 496
498 497 chiSq1=lsq1[1];
499 498
500 499
501 500 if fatspectra<1.0 and powerwidth<4:
502 501 choice=0
503 502 Amplitude0=lsq1[0][2]
504 503 shift0=lsq1[0][0]
505 504 width0=lsq1[0][1]
506 505 p0=lsq1[0][3]
507 506 Amplitude1=0.
508 507 shift1=0.
509 508 width1=0.
510 509 p1=0.
511 510 noise=lsq1[0][4]
512 511 #return (numpy.array([shift0,width0,Amplitude0,p0]),
513 512 # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
514 513
515 514 ''' two gaussians '''
516 515 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
517 516 shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
518 517 shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
519 518 width0=powerwidth/6.;
520 519 width1=width0
521 520 power0=2.;
522 521 power1=power0
523 522 amplitude0=firstamp;
524 523 amplitude1=secondamp
525 524 state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
526 525 #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
527 526 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
528 527 #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
529 528
530 529 lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )
531 530
532 531
533 532 chiSq2=lsq2[1];
534 533
535 534
536 535
537 536 oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
538 537
539 538 if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
540 539 if oneG:
541 540 choice=0
542 541 else:
543 542 w1=lsq2[0][1]; w2=lsq2[0][5]
544 543 a1=lsq2[0][2]; a2=lsq2[0][6]
545 544 p1=lsq2[0][3]; p2=lsq2[0][7]
546 545 s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1;
547 546 s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
548 547 gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
549 548
550 549 if gp1>gp2:
551 550 if a1>0.7*a2:
552 551 choice=1
553 552 else:
554 553 choice=2
555 554 elif gp2>gp1:
556 555 if a2>0.7*a1:
557 556 choice=2
558 557 else:
559 558 choice=1
560 559 else:
561 560 choice=numpy.argmax([a1,a2])+1
562 561 #else:
563 562 #choice=argmin([std2a,std2b])+1
564 563
565 564 else: # with low SNR go to the most energetic peak
566 565 choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
567 566
568 567
569 568 shift0=lsq2[0][0];
570 569 vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
571 570 shift1=lsq2[0][4];
572 571 vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
573 572
574 573 max_vel = 1.0
575 574
576 575 #first peak will be 0, second peak will be 1
577 576 if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range
578 577 shift0=lsq2[0][0]
579 578 width0=lsq2[0][1]
580 579 Amplitude0=lsq2[0][2]
581 580 p0=lsq2[0][3]
582 581
583 582 shift1=lsq2[0][4]
584 583 width1=lsq2[0][5]
585 584 Amplitude1=lsq2[0][6]
586 585 p1=lsq2[0][7]
587 586 noise=lsq2[0][8]
588 587 else:
589 588 shift1=lsq2[0][0]
590 589 width1=lsq2[0][1]
591 590 Amplitude1=lsq2[0][2]
592 591 p1=lsq2[0][3]
593 592
594 593 shift0=lsq2[0][4]
595 594 width0=lsq2[0][5]
596 595 Amplitude0=lsq2[0][6]
597 596 p0=lsq2[0][7]
598 597 noise=lsq2[0][8]
599 598
600 599 if Amplitude0<0.05: # in case the peak is noise
601 600 shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN]
602 601 if Amplitude1<0.05:
603 602 shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN]
604 603
605 604
606 605 # if choice==0: # pick the single gaussian fit
607 606 # Amplitude0=lsq1[0][2]
608 607 # shift0=lsq1[0][0]
609 608 # width0=lsq1[0][1]
610 609 # p0=lsq1[0][3]
611 610 # Amplitude1=0.
612 611 # shift1=0.
613 612 # width1=0.
614 613 # p1=0.
615 614 # noise=lsq1[0][4]
616 615 # elif choice==1: # take the first one of the 2 gaussians fitted
617 616 # Amplitude0 = lsq2[0][2]
618 617 # shift0 = lsq2[0][0]
619 618 # width0 = lsq2[0][1]
620 619 # p0 = lsq2[0][3]
621 620 # Amplitude1 = lsq2[0][6] # This is 0 in gg1
622 621 # shift1 = lsq2[0][4] # This is 0 in gg1
623 622 # width1 = lsq2[0][5] # This is 0 in gg1
624 623 # p1 = lsq2[0][7] # This is 0 in gg1
625 624 # noise = lsq2[0][8]
626 625 # else: # the second one
627 626 # Amplitude0 = lsq2[0][6]
628 627 # shift0 = lsq2[0][4]
629 628 # width0 = lsq2[0][5]
630 629 # p0 = lsq2[0][7]
631 630 # Amplitude1 = lsq2[0][2] # This is 0 in gg1
632 631 # shift1 = lsq2[0][0] # This is 0 in gg1
633 632 # width1 = lsq2[0][1] # This is 0 in gg1
634 633 # p1 = lsq2[0][3] # This is 0 in gg1
635 634 # noise = lsq2[0][8]
636 635
637 636 #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0)
638 637 SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
639 638 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
640 639 #print 'SPC_ch1.shape',SPC_ch1.shape
641 640 #print 'SPC_ch2.shape',SPC_ch2.shape
642 641 #dataOut.data_param = SPC_ch1
643 642 SPCparam = (SPC_ch1,SPC_ch2)
644 643 #GauSPC[1] = SPC_ch2
645 644
646 645 # print 'shift0', shift0
647 646 # print 'Amplitude0', Amplitude0
648 647 # print 'width0', width0
649 648 # print 'p0', p0
650 649 # print '========================'
651 650 # print 'shift1', shift1
652 651 # print 'Amplitude1', Amplitude1
653 652 # print 'width1', width1
654 653 # print 'p1', p1
655 654 # print 'noise', noise
656 655 # print 's_noise', wnoise
657 656
658 657 return GauSPC
659 658
660 659 def y_model1(self,x,state):
661 660 shift0,width0,amplitude0,power0,noise=state
662 661 model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
663 662
664 663 model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
665 664
666 665 model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
667 666 return model0+model0u+model0d+noise
668 667
669 668 def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
670 669 shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state
671 670 model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
672 671
673 672 model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
674 673
675 674 model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
676 675 model1=amplitude1*numpy.exp(-0.5*abs((x-shift1)/width1)**power1)
677 676
678 677 model1u=amplitude1*numpy.exp(-0.5*abs((x-shift1- self.Num_Bin )/width1)**power1)
679 678
680 679 model1d=amplitude1*numpy.exp(-0.5*abs((x-shift1+ self.Num_Bin )/width1)**power1)
681 680 return model0+model0u+model0d+model1+model1u+model1d+noise
682 681
683 682 def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is.
684 683
685 684 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model1(x,state)))**2)#/(64-5.) # /(64-5.) can be commented
686 685
687 686 def misfit2(self,state,y_data,x,num_intg):
688 687 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
689 688
690 689
691 690
692 691 class PrecipitationProc(Operation):
693 692
694 693 '''
695 694 Operator that estimates Reflectivity factor (Z), and estimates rainfall Rate (R)
696 695
697 696 Input:
698 697 self.dataOut.data_pre : SelfSpectra
699 698
700 699 Output:
701 700
702 701 self.dataOut.data_output : Reflectivity factor, rainfall Rate
703 702
704 703
705 704 Parameters affected:
706 705 '''
707 706 def gaus(self,xSamples,Amp,Mu,Sigma):
708 707 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
709 708
710 709
711 710
712 711 def Moments(self, ySamples, xSamples):
713 712 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
714 713 yNorm = ySamples / Pot
715 714
716 715 Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
717 716 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
718 717 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
719 718
720 719 return numpy.array([Pot, Vr, Desv])
721 720
722 721 def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118,
723 722 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km = 0.93, Altitude=3350):
724 723
725 724
726 725 Velrange = dataOut.spcparam_range[2]
727 726 FrecRange = dataOut.spcparam_range[0]
728 727
729 728 dV= Velrange[1]-Velrange[0]
730 729 dF= FrecRange[1]-FrecRange[0]
731 730
732 731 if radar == "MIRA35C" :
733 732
734 733 self.spc = dataOut.data_pre[0].copy()
735 734 self.Num_Hei = self.spc.shape[2]
736 735 self.Num_Bin = self.spc.shape[1]
737 736 self.Num_Chn = self.spc.shape[0]
738 737 Ze = self.dBZeMODE2(dataOut)
739 738
740 739 else:
741 740
742 741 self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() #
743 742 self.Num_Hei = self.spc.shape[2]
744 743 self.Num_Bin = self.spc.shape[1]
745 744 self.Num_Chn = self.spc.shape[0]
746 745 print '==================== SPC SHAPE',numpy.shape(self.spc)
747 746
748 747
749 748 ''' Se obtiene la constante del RADAR '''
750 749
751 750 self.Pt = Pt
752 751 self.Gt = Gt
753 752 self.Gr = Gr
754 753 self.Lambda = Lambda
755 754 self.aL = aL
756 755 self.tauW = tauW
757 756 self.ThetaT = ThetaT
758 757 self.ThetaR = ThetaR
759 758
760 759 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
761 760 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
762 RadarConstant = 1e-10 * Numerator / Denominator #
761 RadarConstant = 5e-27 * Numerator / Denominator #
763 762 print '***'
764 763 print '*** RadarConstant' , RadarConstant, '****'
765 764 print '***'
766 765 ''' ============================= '''
767 766
768 767 self.spc[0] = (self.spc[0]-dataOut.noise[0])
769 768 self.spc[1] = (self.spc[1]-dataOut.noise[1])
770 769 self.spc[2] = (self.spc[2]-dataOut.noise[2])
771 770
772 771 self.spc[ numpy.where(self.spc < 0)] = 0
773 772
774 773 SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise))
775 774 SPCmean[ numpy.where(SPCmean < 0)] = 0
776 775
777 776 ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei])
778 777 ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei])
779 778 ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei])
780 779
781 780 Pr = SPCmean[:,:]
782 781
783 782 VelMeteoro = numpy.mean(SPCmean,axis=0)
784 783
785 784 #print '==================== Vel SHAPE',VelMeteoro
786 785
787 786 D_range = numpy.zeros([self.Num_Bin,self.Num_Hei])
788 787 SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei])
789 788 N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei])
790 D_mean = numpy.zeros(self.Num_Hei)
789 V_mean = numpy.zeros(self.Num_Hei)
791 790 del_V = numpy.zeros(self.Num_Hei)
792 791 Z = numpy.zeros(self.Num_Hei)
793 792 Ze = numpy.zeros(self.Num_Hei)
794 793 RR = numpy.zeros(self.Num_Hei)
795 794
796 795 Range = dataOut.heightList*1000.
797 796
798 797 for R in range(self.Num_Hei):
799
798
800 799 h = Range[R] + Altitude #Range from ground to radar pulse altitude
801 800 del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
802 801
803 802 D_range[:,R] = numpy.log( (9.65 - (Velrange[0:self.Num_Bin] / del_V[R])) / 10.3 ) / -0.6 #Diameter range [m]x10**-3
804 803
805 804 '''NOTA: ETA(n) dn = ETA(f) df
806 805
807 806 dn = 1 Diferencial de muestreo
808 807 df = ETA(n) / ETA(f)
809 808
810 809 '''
811 810
812 811 ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA)
813 812
814 813 ETAv[:,R]=ETAn[:,R]/dV
815 814
816 815 ETAd[:,R]=ETAv[:,R]*6.18*exp(-0.6*D_range[:,R])
817 816
818 817 SIGMA[:,R] = Km * (D_range[:,R] * 1e-3 )**6 * numpy.pi**5 / Lambda**4 #Equivalent Section of drops (sigma)
819 818
820 819 N_dist[:,R] = ETAn[:,R] / SIGMA[:,R]
821 820
822 DMoments = self.Moments(Pr[:,R], D_range[:,R])
821 DMoments = self.Moments(Pr[:,R], Velrange[0:self.Num_Bin])
823 822
824 823 try:
825 popt01,pcov = curve_fit(self.gaus, D_range[:,R] , Pr[:,R] , p0=DMoments)
824 popt01,pcov = curve_fit(self.gaus, Velrange[0:self.Num_Bin] , Pr[:,R] , p0=DMoments)
826 825 except:
827 826 popt01=numpy.zeros(3)
828 827 popt01[1]= DMoments[1]
829 D_mean[R]=popt01[1]
828 if popt01[1]<0 or popt01[1]>20:
829 popt01[1]=numpy.NaN
830
831 V_mean[R]=popt01[1]
830 832
831 Z[R] = numpy.nansum( N_dist[:,R] * (D_range[:,R])**6 )*1e-18
833 Z[R] = numpy.nansum( N_dist[:,R] * (D_range[:,R])**6 )#*10**-18
832 834
833 RR[R] = 3.6e-6*1e-9*6*10**-4.*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate
835 RR[R] = 0.0006*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate
834 836
835 Ze[R] = (numpy.nansum( ETAn[:,R]) * Lambda**4) / ( numpy.pi**5 * Km)
837 Ze[R] = (numpy.nansum( ETAn[:,R]) * Lambda**4) / ( 10**-18*numpy.pi**5 * Km)
836 838
837 839
838 840
839 841 RR2 = (Z/200)**(1/1.6)
840 842 dBRR = 10*numpy.log10(RR)
841 843 dBRR2 = 10*numpy.log10(RR2)
842 844
843 845 dBZe = 10*numpy.log10(Ze)
844 846 dBZ = 10*numpy.log10(Z)
845 847
846 848 dataOut.data_output = Z
847 849 dataOut.data_param = numpy.ones([3,self.Num_Hei])
848 850 dataOut.channelList = [0,1,2]
849 851
850 852 dataOut.data_param[0]=dBZ
851 dataOut.data_param[1]=RR2
853 dataOut.data_param[1]=V_mean
852 854 dataOut.data_param[2]=RR
853 855
854 856 #print 'VELRANGE', Velrange
855 print 'Range', len(Range)
856 print 'delv',del_V
857 #print 'DRANGE', D_range[:,50]
858 print 'NOISE', dataOut.noise[0], 10*numpy.log10(dataOut.noise[0])
859 print 'radarconstant', RadarConstant
860 print 'Range', Range
861 # print 'ETAn SHAPE', ETAn.shape
862 # print 'ETAn ', numpy.nansum(ETAn, axis=0)
863 # print 'ETAd ', numpy.nansum(ETAd, axis=0)
864 print 'Pr ', numpy.nansum(Pr, axis=0)
865 print 'dataOut.SPCparam[1]', numpy.nansum(dataOut.SPCparam[1][0], axis=0)
857 #print 'Range', len(Range)
858 #print 'delv',del_V
859 # print 'DRANGE', D_range[:,56]
860 # #print 'NOISE', dataOut.noise[0], 10*numpy.log10(dataOut.noise[0])
861 # print 'radarconstant', RadarConstant
862 # #print 'ETAn SHAPE', ETAn.shape
863 # # print 'ETAn ', numpy.nansum(ETAn, axis=0)
864 # # print 'ETAd ', numpy.nansum(ETAd, axis=0)
865 # print 'Pr ', numpy.nansum(Pr, axis=0)
866 # print 'dataOut.SPCparam[1]', numpy.nansum(dataOut.SPCparam[1][0], axis=0)
866 867 # print 'Ze ', dBZe
867 # print 'Z ', dBZ
868 print 'RR2 ', RR2
868 print 'Z ', dBZ
869 # print 'Ndist',N_dist[:,56]
870 # #print 'RR2 ', RR2
869 871 print 'RR ', RR
872 print 'Vr', V_mean
870 873 #print 'RR2 ', dBRR2
871 874 #print 'D_mean', D_mean
872 875 #print 'del_V', del_V
873 876 #print 'D_range',D_range.shape, D_range[:,30]
874 877 #print 'Velrange', Velrange
875 878 #print 'numpy.nansum( N_dist[:,R]', numpy.nansum( N_dist, axis=0)
876 879 #print 'dataOut.data_param SHAPE', dataOut.data_param.shape
877 880
878 881
879 882 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
880 883
881 884 NPW = dataOut.NPW
882 885 COFA = dataOut.COFA
883 886
884 887 SNR = numpy.array([self.spc[0,:,:] / NPW[0]]) #, self.spc[1,:,:] / NPW[1]])
885 888 RadarConst = dataOut.RadarConst
886 889 #frequency = 34.85*10**9
887 890
888 891 ETA = numpy.zeros(([self.Num_Chn ,self.Num_Hei]))
889 892 data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN
890 893
891 894 ETA = numpy.sum(SNR,1)
892 895 print 'ETA' , ETA
893 896 ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
894 897
895 898 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
896 899
897 900 for r in range(self.Num_Hei):
898 901
899 902 Ze[0,r] = ( ETA[0,r] ) * COFA[0,r][0] * RadarConst * ((r/5000.)**2)
900 903 #Ze[1,r] = ( ETA[1,r] ) * COFA[1,r][0] * RadarConst * ((r/5000.)**2)
901 904
902 905 return Ze
903 906
904 907 # def GetRadarConstant(self):
905 908 #
906 909 # """
907 910 # Constants:
908 911 #
909 912 # Pt: Transmission Power dB 5kW 5000
910 913 # Gt: Transmission Gain dB 24.7 dB 295.1209
911 914 # Gr: Reception Gain dB 18.5 dB 70.7945
912 915 # Lambda: Wavelenght m 0.6741 m 0.6741
913 916 # aL: Attenuation loses dB 4dB 2.5118
914 917 # tauW: Width of transmission pulse s 4us 4e-6
915 918 # ThetaT: Transmission antenna bean angle rad 0.1656317 rad 0.1656317
916 919 # ThetaR: Reception antenna beam angle rad 0.36774087 rad 0.36774087
917 920 #
918 921 # """
919 922 #
920 923 # Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
921 924 # Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
922 925 # RadarConstant = Numerator / Denominator
923 926 #
924 927 # return RadarConstant
925 928
926 929
927 930
928 931 class FullSpectralAnalysis(Operation):
929 932
930 933 """
931 934 Function that implements Full Spectral Analisys technique.
932 935
933 936 Input:
934 937 self.dataOut.data_pre : SelfSpectra and CrossSPectra data
935 938 self.dataOut.groupList : Pairlist of channels
936 939 self.dataOut.ChanDist : Physical distance between receivers
937 940
938 941
939 942 Output:
940 943
941 944 self.dataOut.data_output : Zonal wind, Meridional wind and Vertical wind
942 945
943 946
944 947 Parameters affected: Winds, height range, SNR
945 948
946 949 """
947 950 def run(self, dataOut, E01=None, E02=None, E12=None, N01=None, N02=None, N12=None, SNRlimit=7):
948 951
949 952 self.indice=int(numpy.random.rand()*1000)
950 953
951 954 spc = dataOut.data_pre[0].copy()
952 955 cspc = dataOut.data_pre[1]
953 956
954 957 nChannel = spc.shape[0]
955 958 nProfiles = spc.shape[1]
956 959 nHeights = spc.shape[2]
957 960
958 961 pairsList = dataOut.groupList
959 962 if dataOut.ChanDist is not None :
960 963 ChanDist = dataOut.ChanDist
961 964 else:
962 965 ChanDist = numpy.array([[E01, N01],[E02,N02],[E12,N12]])
963 966
964 967 FrecRange = dataOut.spc_range[0]
965 968
966 969 ySamples=numpy.ones([nChannel,nProfiles])
967 970 phase=numpy.ones([nChannel,nProfiles])
968 971 CSPCSamples=numpy.ones([nChannel,nProfiles],dtype=numpy.complex_)
969 972 coherence=numpy.ones([nChannel,nProfiles])
970 973 PhaseSlope=numpy.ones(nChannel)
971 974 PhaseInter=numpy.ones(nChannel)
972 975 data_SNR=numpy.zeros([nProfiles])
973 976
974 977 data = dataOut.data_pre
975 978 noise = dataOut.noise
976 979
977 980 dataOut.data_SNR = (numpy.mean(spc,axis=1)- noise[0]) / noise[0]
978 981
979 982 dataOut.data_SNR[numpy.where( dataOut.data_SNR <0 )] = 1e-20
980 983
981 984
982 985 #FirstMoment = dataOut.moments[0,1,:]#numpy.average(dataOut.data_param[:,1,:],0)
983 986 #SecondMoment = numpy.average(dataOut.moments[:,2,:],0)
984 987
985 988 #SNRdBMean = []
986 989
987 990 data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN
988 991
989 992 velocityX=[]
990 993 velocityY=[]
991 994 velocityV=[]
992 995 PhaseLine=[]
993 996
994 997 dbSNR = 10*numpy.log10(dataOut.data_SNR)
995 998 dbSNR = numpy.average(dbSNR,0)
996 999
997 1000 for Height in range(nHeights):
998 1001
999 1002 [Vzon,Vmer,Vver, GaussCenter, PhaseSlope, FitGaussCSPC]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, dataOut.spc_range.copy(), dbSNR[Height], SNRlimit)
1000 1003 PhaseLine = numpy.append(PhaseLine, PhaseSlope)
1001 1004
1002 1005 if abs(Vzon)<100. and abs(Vzon)> 0.:
1003 1006 velocityX=numpy.append(velocityX, -Vzon)#Vmag
1004 1007
1005 1008 else:
1006 1009 #print 'Vzon',Vzon
1007 1010 velocityX=numpy.append(velocityX, numpy.NaN)
1008 1011
1009 1012 if abs(Vmer)<100. and abs(Vmer) > 0.:
1010 1013 velocityY=numpy.append(velocityY, -Vmer)#Vang
1011 1014
1012 1015 else:
1013 1016 #print 'Vmer',Vmer
1014 1017 velocityY=numpy.append(velocityY, numpy.NaN)
1015 1018
1016 1019 if dbSNR[Height] > SNRlimit:
1017 1020 velocityV=numpy.append(velocityV, -Vver)#FirstMoment[Height])
1018 1021 else:
1019 1022 velocityV=numpy.append(velocityV, numpy.NaN)
1020 1023 #FirstMoment[Height]= numpy.NaN
1021 1024 # if SNRdBMean[Height] <12:
1022 1025 # FirstMoment[Height] = numpy.NaN
1023 1026 # velocityX[Height] = numpy.NaN
1024 1027 # velocityY[Height] = numpy.NaN
1025 1028
1026 1029
1027 1030
1028 1031 data_output[0] = numpy.array(velocityX) #self.moving_average(numpy.array(velocityX) , N=1)
1029 1032 data_output[1] = numpy.array(velocityY) #self.moving_average(numpy.array(velocityY) , N=1)
1030 1033 data_output[2] = -velocityV#FirstMoment
1031 1034
1032 print 'FirstMoment', data_output[2]
1035 print 'data_output', data_output.shape
1033 1036 #print FirstMoment
1034 1037 # print 'velocityX',numpy.shape(data_output[0])
1035 1038 # print 'velocityX',data_output[0]
1036 1039 # print ' '
1037 1040 # print 'velocityY',numpy.shape(data_output[1])
1038 1041 # print 'velocityY',data_output[1]
1039 1042 # print 'velocityV',data_output[2]
1040 1043 # print 'PhaseLine',PhaseLine
1041 1044 #print numpy.array(velocityY)
1042 1045 #print 'SNR'
1043 1046 #print 10*numpy.log10(dataOut.data_SNR)
1044 1047 #print numpy.shape(10*numpy.log10(dataOut.data_SNR))
1045 1048 print ' '
1046 1049
1047 1050 xFrec=FrecRange[0:spc.shape[1]]
1048 1051
1049 1052 dataOut.data_output=data_output
1050 1053
1051 1054 return
1052 1055
1053 1056
1054 1057 def moving_average(self,x, N=2):
1055 1058 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
1056 1059
1057 1060 def gaus(self,xSamples,Amp,Mu,Sigma):
1058 1061 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
1059 1062
1060 1063
1061
1064
1062 1065 def Moments(self, ySamples, xSamples):
1063 1066 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
1064 1067 yNorm = ySamples / Pot
1065 1068 Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
1066 1069 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
1067 1070 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
1068 1071
1069 1072 return numpy.array([Pot, Vr, Desv])
1070 1073
1071 1074 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
1072 1075
1073 1076
1074 1077
1075 1078 ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
1076 1079 phase=numpy.ones([spc.shape[0],spc.shape[1]])
1077 1080 CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_)
1078 1081 coherence=numpy.ones([spc.shape[0],spc.shape[1]])
1079 1082 PhaseSlope=numpy.zeros(spc.shape[0])
1080 1083 PhaseInter=numpy.ones(spc.shape[0])
1081 1084 xFrec=AbbsisaRange[0][0:spc.shape[1]]
1082 1085 xVel =AbbsisaRange[2][0:spc.shape[1]]
1083 1086 Vv=numpy.empty(spc.shape[2])*0
1084 1087 SPCav = numpy.average(spc, axis=0)-numpy.average(noise) #spc[0]-noise[0]#
1085 1088
1086 1089 SPCmoments = self.Moments(SPCav[:,Height], xVel )
1087 1090 CSPCmoments = []
1088 1091 cspcNoise = numpy.empty(3)
1089 1092
1090 1093 '''Getting Eij and Nij'''
1091 1094
1092 1095 E01=ChanDist[0][0]
1093 1096 N01=ChanDist[0][1]
1094 1097
1095 1098 E02=ChanDist[1][0]
1096 1099 N02=ChanDist[1][1]
1097 1100
1098 1101 E12=ChanDist[2][0]
1099 1102 N12=ChanDist[2][1]
1100 1103
1101 1104 z = spc.copy()
1102 1105 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1103 1106
1104 1107 for i in range(spc.shape[0]):
1105 1108
1106 1109 '''****** Line of Data SPC ******'''
1107 1110 zline=z[i,:,Height].copy() - noise[i] # Se resta ruido
1108 1111
1109 1112 '''****** SPC is normalized ******'''
1110 1113 SmoothSPC =self.moving_average(zline.copy(),N=1) # Se suaviza el ruido
1111 1114 FactNorm = SmoothSPC/numpy.nansum(SmoothSPC) # SPC Normalizado y suavizado
1112 1115
1113 1116 xSamples = xFrec # Se toma el rango de frecuncias
1114 1117 ySamples[i] = FactNorm # Se toman los valores de SPC normalizado
1115 1118
1116 1119 for i in range(spc.shape[0]):
1117 1120
1118 1121 '''****** Line of Data CSPC ******'''
1119 1122 cspcLine = ( cspc[i,:,Height].copy())# - noise[i] ) # no! Se resta el ruido
1120 1123 SmoothCSPC =self.moving_average(cspcLine,N=1) # Se suaviza el ruido
1121 1124 cspcNorm = SmoothCSPC/numpy.nansum(SmoothCSPC) # CSPC normalizado y suavizado
1122 1125
1123 1126 '''****** CSPC is normalized with respect to Briggs and Vincent ******'''
1124 1127 chan_index0 = pairsList[i][0]
1125 1128 chan_index1 = pairsList[i][1]
1126 1129
1127 1130 CSPCFactor= numpy.abs(numpy.nansum(ySamples[chan_index0]))**2 * numpy.abs(numpy.nansum(ySamples[chan_index1]))**2
1128 1131 CSPCNorm = cspcNorm / numpy.sqrt(CSPCFactor)
1129 1132
1130 1133 CSPCSamples[i] = CSPCNorm
1131 1134
1132 1135 coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
1133 1136
1134 1137 #coherence[i]= self.moving_average(coherence[i],N=1)
1135 1138
1136 1139 phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
1137 1140
1138 1141 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPCSamples[0]), xSamples),
1139 1142 self.Moments(numpy.abs(CSPCSamples[1]), xSamples),
1140 1143 self.Moments(numpy.abs(CSPCSamples[2]), xSamples)])
1141 1144
1142 1145 #print '##### SUMA de SPC #####', len(ySamples)
1143 1146 #print numpy.sum(ySamples[0])
1144 1147 #print '##### SUMA de CSPC #####', len(coherence)
1145 1148 #print numpy.sum(numpy.abs(CSPCNorm))
1146 1149 #print numpy.sum(coherence[0])
1147 1150 # print 'len',len(xSamples)
1148 1151 # print 'CSPCmoments', numpy.shape(CSPCmoments)
1149 1152 # print CSPCmoments
1150 1153 # print '#######################'
1151 1154
1152 1155 popt=[1e-10,0,1e-10]
1153 1156 popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
1154 1157 FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0
1155 1158
1156 1159 CSPCMask01 = numpy.abs(CSPCSamples[0])
1157 1160 CSPCMask02 = numpy.abs(CSPCSamples[1])
1158 1161 CSPCMask12 = numpy.abs(CSPCSamples[2])
1159 1162
1160 1163 mask01 = ~numpy.isnan(CSPCMask01)
1161 1164 mask02 = ~numpy.isnan(CSPCMask02)
1162 1165 mask12 = ~numpy.isnan(CSPCMask12)
1163 1166
1164 1167 #mask = ~numpy.isnan(CSPCMask01)
1165 1168 CSPCMask01 = CSPCMask01[mask01]
1166 1169 CSPCMask02 = CSPCMask02[mask02]
1167 1170 CSPCMask12 = CSPCMask12[mask12]
1168 1171 #CSPCMask01 = numpy.ma.masked_invalid(CSPCMask01)
1169 1172
1170 1173
1171 1174
1172 1175 '''***Fit Gauss CSPC01***'''
1173 1176 if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3 :
1174 1177 try:
1175 1178 popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
1176 1179 popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1])
1177 1180 popt12,pcov = curve_fit(self.gaus,xSamples[mask12],numpy.abs(CSPCMask12),p0=CSPCmoments[2])
1178 1181 FitGauss01 = self.gaus(xSamples,*popt01)
1179 1182 FitGauss02 = self.gaus(xSamples,*popt02)
1180 1183 FitGauss12 = self.gaus(xSamples,*popt12)
1181 1184 except:
1182 1185 FitGauss01=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[0]))
1183 1186 FitGauss02=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[1]))
1184 1187 FitGauss12=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[2]))
1185 1188
1186 1189
1187 1190 CSPCopt = numpy.vstack([popt01,popt02,popt12])
1188 1191
1189 1192 '''****** Getting fij width ******'''
1190 1193
1191 1194 yMean = numpy.average(ySamples, axis=0) # ySamples[0]
1192 1195
1193 1196 '''******* Getting fitting Gaussian *******'''
1194 1197 meanGauss = sum(xSamples*yMean) / len(xSamples) # Mu, velocidad radial (frecuencia)
1195 1198 sigma2 = sum(yMean*(xSamples-meanGauss)**2) / len(xSamples) # Varianza, Ancho espectral (frecuencia)
1196 1199
1197 1200 yMoments = self.Moments(yMean, xSamples)
1198 1201
1199 1202 if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3: # and abs(meanGauss/sigma2) > 0.00001:
1200 1203 try:
1201 1204 popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=yMoments)
1202 1205 FitGauss=self.gaus(xSamples,*popt)
1203 1206
1204 1207 except :#RuntimeError:
1205 1208 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1206 1209
1207 1210
1208 1211 else:
1209 1212 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1210 1213
1211 1214
1212 1215
1213 1216 '''****** Getting Fij ******'''
1214 1217 Fijcspc = CSPCopt[:,2]/2*3
1215 1218
1216 1219
1217 1220 GaussCenter = popt[1] #xFrec[GCpos]
1218 1221 #Punto en Eje X de la Gaussiana donde se encuentra el centro
1219 1222 ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()]
1220 1223 PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0]
1221 1224
1222 1225 #Punto e^-1 hubicado en la Gaussiana
1223 1226 PeMinus1 = numpy.max(FitGauss)* numpy.exp(-1)
1224 1227 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss"
1225 1228 PointFij = numpy.where(FitGauss==FijClosest)[0][0]
1226 1229
1227 1230 if xSamples[PointFij] > xSamples[PointGauCenter]:
1228 1231 Fij = xSamples[PointFij] - xSamples[PointGauCenter]
1229 1232
1230 1233 else:
1231 1234 Fij = xSamples[PointGauCenter] - xSamples[PointFij]
1232 1235
1233 1236 # print 'CSPCopt'
1234 1237 # print CSPCopt
1235 1238 # print 'popt'
1236 1239 # print popt
1237 1240 # print '#######################################'
1238 1241 #print 'dataOut.data_param', numpy.shape(data_param)
1239 1242 #print 'dataOut.data_param0', data_param[0,0,Height]
1240 1243 #print 'dataOut.data_param1', data_param[0,1,Height]
1241 1244 #print 'dataOut.data_param2', data_param[0,2,Height]
1242 1245
1243 1246
1244 1247 # print 'yMoments', yMoments
1245 1248 # print 'Moments', SPCmoments
1246 1249 # print 'Fij2 Moment', Fij
1247 1250 # #print 'Fij', Fij, 'popt[2]/2',popt[2]/2
1248 1251 # print 'Fijcspc',Fijcspc
1249 1252 # print '#######################################'
1250 1253
1251 1254
1252 1255 '''****** Taking frequency ranges from SPCs ******'''
1253 1256
1254 1257
1255 1258 #GaussCenter = popt[1] #Primer momento 01
1256 1259 GauWidth = popt[2] *3/2 #Ancho de banda de Gau01
1257 1260 Range = numpy.empty(2)
1258 1261 Range[0] = GaussCenter - GauWidth
1259 1262 Range[1] = GaussCenter + GauWidth
1260 1263 #Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max)
1261 1264 ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()]
1262 1265 ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()]
1263 1266
1264 1267 PointRangeMin = numpy.where(xSamples==ClosRangeMin)[0][0]
1265 1268 PointRangeMax = numpy.where(xSamples==ClosRangeMax)[0][0]
1266 1269
1267 1270 Range=numpy.array([ PointRangeMin, PointRangeMax ])
1268 1271
1269 1272 FrecRange = xFrec[ Range[0] : Range[1] ]
1270 1273 VelRange = xVel[ Range[0] : Range[1] ]
1271 1274
1272 1275
1273 1276 #print 'RANGE: ', Range
1274 1277 #print 'FrecRange', numpy.shape(FrecRange)#,FrecRange
1275 1278 #print 'len: ', len(FrecRange)
1276 1279
1277 1280 '''****** Getting SCPC Slope ******'''
1278 1281
1279 1282 for i in range(spc.shape[0]):
1280 1283
1281 1284 if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.3:
1282 1285 PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
1283 1286
1284 1287 #print 'Ancho espectral Frecuencias', FrecRange[-1]-FrecRange[0], 'Hz'
1285 1288 #print 'Ancho espectral Velocidades', VelRange[-1]-VelRange[0], 'm/s'
1286 1289 #print 'FrecRange', len(FrecRange) , FrecRange
1287 1290 #print 'VelRange', len(VelRange) , VelRange
1288 1291 #print 'PhaseRange', numpy.shape(PhaseRange), PhaseRange
1289 1292 #print ' '
1290 1293
1291 1294 '''***********************VelRange******************'''
1292 1295
1293 1296 mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
1294 1297
1295 1298 if len(FrecRange) == len(PhaseRange):
1296 1299 try:
1297 1300 slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange[mask], PhaseRange[mask])
1298 1301 PhaseSlope[i]=slope
1299 1302 PhaseInter[i]=intercept
1300 1303 except:
1301 1304 PhaseSlope[i]=0
1302 1305 PhaseInter[i]=0
1303 1306 else:
1304 1307 PhaseSlope[i]=0
1305 1308 PhaseInter[i]=0
1306 1309 else:
1307 1310 PhaseSlope[i]=0
1308 1311 PhaseInter[i]=0
1309 1312
1310 1313
1311 1314 '''Getting constant C'''
1312 1315 cC=(Fij*numpy.pi)**2
1313 1316
1314 1317 '''****** Getting constants F and G ******'''
1315 1318 MijEijNij=numpy.array([[E02,N02], [E12,N12]])
1316 1319 MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
1317 1320 MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
1318 1321 MijResults=numpy.array([MijResult0,MijResult1])
1319 1322 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1320 1323
1321 1324 '''****** Getting constants A, B and H ******'''
1322 1325 W01=numpy.nanmax( FitGauss01 ) #numpy.abs(CSPCSamples[0]))
1323 1326 W02=numpy.nanmax( FitGauss02 ) #numpy.abs(CSPCSamples[1]))
1324 1327 W12=numpy.nanmax( FitGauss12 ) #numpy.abs(CSPCSamples[2]))
1325 1328
1326 1329 WijResult0=((cF*E01+cG*N01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
1327 1330 WijResult1=((cF*E02+cG*N02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
1328 1331 WijResult2=((cF*E12+cG*N12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
1329 1332
1330 1333 WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
1331 1334
1332 1335 WijEijNij=numpy.array([ [E01**2, N01**2, 2*E01*N01] , [E02**2, N02**2, 2*E02*N02] , [E12**2, N12**2, 2*E12*N12] ])
1333 1336 (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1334 1337
1335 1338 VxVy=numpy.array([[cA,cH],[cH,cB]])
1336 1339 VxVyResults=numpy.array([-cF,-cG])
1337 1340 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1338 1341
1339 1342 #print 'MijResults, cC, PhaseSlope', MijResults, cC, PhaseSlope
1340 1343 #print 'W01,02,12', W01, W02, W12
1341 1344 #print 'WijResult0,1,2',WijResult0, WijResult1, WijResult2, 'Results', WijResults
1342 1345 #print 'cA,cB,cH, cF, cG', cA, cB, cH, cF, cG
1343 1346 #print 'VxVy', VxVyResults
1344 1347 #print '###########################****************************************'
1345 1348 Vzon = Vy
1346 1349 Vmer = Vx
1347 1350 Vmag=numpy.sqrt(Vzon**2+Vmer**2)
1348 1351 Vang=numpy.arctan2(Vmer,Vzon)
1349 1352 if numpy.abs( popt[1] ) < 3.5 and len(FrecRange)>4:
1350 1353 Vver=popt[1]
1351 1354 else:
1352 1355 Vver=numpy.NaN
1353 1356 FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12])
1354 1357
1355 1358
1356 1359 # ''' Ploteo por altura '''
1357 1360 # if Height == 28:
1358 1361 # for i in range(3):
1359 1362 # #print 'FASE', numpy.shape(phase), y[25]
1360 1363 # #print numpy.shape(coherence)
1361 1364 # fig = plt.figure(10+self.indice)
1362 1365 # #plt.plot( x[0:256],coherence[:,25] )
1363 1366 # #cohAv = numpy.average(coherence[i],1)
1364 1367 # Pendiente = FrecRange * PhaseSlope[i]
1365 1368 # plt.plot( FrecRange, Pendiente)
1366 1369 # plt.plot( xFrec,phase[i])
1367 1370 #
1368 1371 # CSPCmean = numpy.mean(numpy.abs(CSPCSamples),0)
1369 1372 # #plt.plot(xFrec, FitGauss01)
1370 1373 # #plt.plot(xFrec, CSPCmean)
1371 1374 # #plt.plot(xFrec, numpy.abs(CSPCSamples[0]))
1372 1375 # #plt.plot(xFrec, FitGauss)
1373 1376 # #plt.plot(xFrec, yMean)
1374 1377 # #plt.plot(xFrec, numpy.abs(coherence[0]))
1375 1378 #
1376 1379 # #plt.axis([-12, 12, 15, 50])
1377 1380 # #plt.title("%s" %( '%s %s, Channel %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S") , i)))
1378 1381 # plt.ylabel('Desfase [rad]')
1379 1382 # #plt.ylabel('CSPC normalizado')
1380 1383 # plt.xlabel('Frec range [Hz]')
1381 1384
1382 1385 #fig.savefig('/home/erick/Documents/Pics/to{}.png'.format(self.indice))
1383 1386
1384 1387 # plt.show()
1385 1388 # self.indice=self.indice+1
1386 1389
1387 1390
1388 1391
1389 1392
1390 1393
1391 1394 # print 'vzon y vmer', Vzon, Vmer
1392 1395 return Vzon, Vmer, Vver, GaussCenter, PhaseSlope, FitGaussCSPC
1393 1396
1394 1397 class SpectralMoments(Operation):
1395 1398
1396 1399 '''
1397 1400 Function SpectralMoments()
1398 1401
1399 1402 Calculates moments (power, mean, standard deviation) and SNR of the signal
1400 1403
1401 1404 Type of dataIn: Spectra
1402 1405
1403 1406 Configuration Parameters:
1404 1407
1405 1408 dirCosx : Cosine director in X axis
1406 1409 dirCosy : Cosine director in Y axis
1407 1410
1408 1411 elevation :
1409 1412 azimuth :
1410 1413
1411 1414 Input:
1412 1415 channelList : simple channel list to select e.g. [2,3,7]
1413 1416 self.dataOut.data_pre : Spectral data
1414 1417 self.dataOut.abscissaList : List of frequencies
1415 1418 self.dataOut.noise : Noise level per channel
1416 1419
1417 1420 Affected:
1418 1421 self.dataOut.moments : Parameters per channel
1419 1422 self.dataOut.data_SNR : SNR per channel
1420 1423
1421 1424 '''
1422 1425
1423 1426 def run(self, dataOut):
1424 1427
1425 1428 #dataOut.data_pre = dataOut.data_pre[0]
1426 1429 data = dataOut.data_pre[0]
1427 1430 absc = dataOut.abscissaList[:-1]
1428 1431 noise = dataOut.noise
1429 1432 nChannel = data.shape[0]
1430 1433 data_param = numpy.zeros((nChannel, 4, data.shape[2]))
1431 1434
1432 1435 for ind in range(nChannel):
1433 1436 data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] )
1434 1437
1435 1438 dataOut.moments = data_param[:,1:,:]
1436 1439 dataOut.data_SNR = data_param[:,0]
1437 1440 return
1438 1441
1439 1442 def __calculateMoments(self, oldspec, oldfreq, n0,
1440 1443 nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
1441 1444
1442 1445 if (nicoh == None): nicoh = 1
1443 1446 if (graph == None): graph = 0
1444 1447 if (smooth == None): smooth = 0
1445 1448 elif (self.smooth < 3): smooth = 0
1446 1449
1447 1450 if (type1 == None): type1 = 0
1448 1451 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
1449 1452 if (snrth == None): snrth = -3
1450 1453 if (dc == None): dc = 0
1451 1454 if (aliasing == None): aliasing = 0
1452 1455 if (oldfd == None): oldfd = 0
1453 1456 if (wwauto == None): wwauto = 0
1454 1457
1455 1458 if (n0 < 1.e-20): n0 = 1.e-20
1456 1459
1457 1460 freq = oldfreq
1458 1461 vec_power = numpy.zeros(oldspec.shape[1])
1459 1462 vec_fd = numpy.zeros(oldspec.shape[1])
1460 1463 vec_w = numpy.zeros(oldspec.shape[1])
1461 1464 vec_snr = numpy.zeros(oldspec.shape[1])
1462 1465
1463 1466 oldspec = numpy.ma.masked_invalid(oldspec)
1464 1467
1465 1468 for ind in range(oldspec.shape[1]):
1466 1469
1467 1470 spec = oldspec[:,ind]
1468 1471 aux = spec*fwindow
1469 1472 max_spec = aux.max()
1470 1473 m = list(aux).index(max_spec)
1471 1474
1472 1475 #Smooth
1473 1476 if (smooth == 0): spec2 = spec
1474 1477 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
1475 1478
1476 1479 # Calculo de Momentos
1477 1480 bb = spec2[range(m,spec2.size)]
1478 1481 bb = (bb<n0).nonzero()
1479 1482 bb = bb[0]
1480 1483
1481 1484 ss = spec2[range(0,m + 1)]
1482 1485 ss = (ss<n0).nonzero()
1483 1486 ss = ss[0]
1484 1487
1485 1488 if (bb.size == 0):
1486 1489 bb0 = spec.size - 1 - m
1487 1490 else:
1488 1491 bb0 = bb[0] - 1
1489 1492 if (bb0 < 0):
1490 1493 bb0 = 0
1491 1494
1492 1495 if (ss.size == 0): ss1 = 1
1493 1496 else: ss1 = max(ss) + 1
1494 1497
1495 1498 if (ss1 > m): ss1 = m
1496 1499
1497 1500 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
1498 1501 power = ( (spec2[valid] - n0) * fwindow[valid] ).sum()
1499 1502 fd = ( (spec2[valid]- n0) * freq[valid] * fwindow[valid] ).sum() / power
1500 1503
1501 1504 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
1502 1505 snr = (spec2.mean()-n0)/n0
1503 1506
1504 1507 if (snr < 1.e-20) :
1505 1508 snr = 1.e-20
1506 1509
1507 1510 vec_power[ind] = power
1508 1511 vec_fd[ind] = fd
1509 1512 vec_w[ind] = w
1510 1513 vec_snr[ind] = snr
1511 1514
1512 1515 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
1513 1516 return moments
1514 1517
1515 1518 #------------------ Get SA Parameters --------------------------
1516 1519
1517 1520 def GetSAParameters(self):
1518 1521 #SA en frecuencia
1519 1522 pairslist = self.dataOut.groupList
1520 1523 num_pairs = len(pairslist)
1521 1524
1522 1525 vel = self.dataOut.abscissaList
1523 1526 spectra = self.dataOut.data_pre
1524 1527 cspectra = self.dataIn.data_cspc
1525 1528 delta_v = vel[1] - vel[0]
1526 1529
1527 1530 #Calculating the power spectrum
1528 1531 spc_pow = numpy.sum(spectra, 3)*delta_v
1529 1532 #Normalizing Spectra
1530 1533 norm_spectra = spectra/spc_pow
1531 1534 #Calculating the norm_spectra at peak
1532 1535 max_spectra = numpy.max(norm_spectra, 3)
1533 1536
1534 1537 #Normalizing Cross Spectra
1535 1538 norm_cspectra = numpy.zeros(cspectra.shape)
1536 1539
1537 1540 for i in range(num_chan):
1538 1541 norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
1539 1542
1540 1543 max_cspectra = numpy.max(norm_cspectra,2)
1541 1544 max_cspectra_index = numpy.argmax(norm_cspectra, 2)
1542 1545
1543 1546 for i in range(num_pairs):
1544 1547 cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
1545 1548 #------------------- Get Lags ----------------------------------
1546 1549
1547 1550 class SALags(Operation):
1548 1551 '''
1549 1552 Function GetMoments()
1550 1553
1551 1554 Input:
1552 1555 self.dataOut.data_pre
1553 1556 self.dataOut.abscissaList
1554 1557 self.dataOut.noise
1555 1558 self.dataOut.normFactor
1556 1559 self.dataOut.data_SNR
1557 1560 self.dataOut.groupList
1558 1561 self.dataOut.nChannels
1559 1562
1560 1563 Affected:
1561 1564 self.dataOut.data_param
1562 1565
1563 1566 '''
1564 1567 def run(self, dataOut):
1565 1568 data_acf = dataOut.data_pre[0]
1566 1569 data_ccf = dataOut.data_pre[1]
1567 1570 normFactor_acf = dataOut.normFactor[0]
1568 1571 normFactor_ccf = dataOut.normFactor[1]
1569 1572 pairs_acf = dataOut.groupList[0]
1570 1573 pairs_ccf = dataOut.groupList[1]
1571 1574
1572 1575 nHeights = dataOut.nHeights
1573 1576 absc = dataOut.abscissaList
1574 1577 noise = dataOut.noise
1575 1578 SNR = dataOut.data_SNR
1576 1579 nChannels = dataOut.nChannels
1577 1580 # pairsList = dataOut.groupList
1578 1581 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
1579 1582
1580 1583 for l in range(len(pairs_acf)):
1581 1584 data_acf[l,:,:] = data_acf[l,:,:]/normFactor_acf[l,:]
1582 1585
1583 1586 for l in range(len(pairs_ccf)):
1584 1587 data_ccf[l,:,:] = data_ccf[l,:,:]/normFactor_ccf[l,:]
1585 1588
1586 1589 dataOut.data_param = numpy.zeros((len(pairs_ccf)*2 + 1, nHeights))
1587 1590 dataOut.data_param[:-1,:] = self.__calculateTaus(data_acf, data_ccf, absc)
1588 1591 dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc)
1589 1592 return
1590 1593
1591 1594 # def __getPairsAutoCorr(self, pairsList, nChannels):
1592 1595 #
1593 1596 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1594 1597 #
1595 1598 # for l in range(len(pairsList)):
1596 1599 # firstChannel = pairsList[l][0]
1597 1600 # secondChannel = pairsList[l][1]
1598 1601 #
1599 1602 # #Obteniendo pares de Autocorrelacion
1600 1603 # if firstChannel == secondChannel:
1601 1604 # pairsAutoCorr[firstChannel] = int(l)
1602 1605 #
1603 1606 # pairsAutoCorr = pairsAutoCorr.astype(int)
1604 1607 #
1605 1608 # pairsCrossCorr = range(len(pairsList))
1606 1609 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1607 1610 #
1608 1611 # return pairsAutoCorr, pairsCrossCorr
1609 1612
1610 1613 def __calculateTaus(self, data_acf, data_ccf, lagRange):
1611 1614
1612 1615 lag0 = data_acf.shape[1]/2
1613 1616 #Funcion de Autocorrelacion
1614 1617 mean_acf = stats.nanmean(data_acf, axis = 0)
1615 1618
1616 1619 #Obtencion Indice de TauCross
1617 1620 ind_ccf = data_ccf.argmax(axis = 1)
1618 1621 #Obtencion Indice de TauAuto
1619 1622 ind_acf = numpy.zeros(ind_ccf.shape,dtype = 'int')
1620 1623 ccf_lag0 = data_ccf[:,lag0,:]
1621 1624
1622 1625 for i in range(ccf_lag0.shape[0]):
1623 1626 ind_acf[i,:] = numpy.abs(mean_acf - ccf_lag0[i,:]).argmin(axis = 0)
1624 1627
1625 1628 #Obtencion de TauCross y TauAuto
1626 1629 tau_ccf = lagRange[ind_ccf]
1627 1630 tau_acf = lagRange[ind_acf]
1628 1631
1629 1632 Nan1, Nan2 = numpy.where(tau_ccf == lagRange[0])
1630 1633
1631 1634 tau_ccf[Nan1,Nan2] = numpy.nan
1632 1635 tau_acf[Nan1,Nan2] = numpy.nan
1633 1636 tau = numpy.vstack((tau_ccf,tau_acf))
1634 1637
1635 1638 return tau
1636 1639
1637 1640 def __calculateLag1Phase(self, data, lagTRange):
1638 1641 data1 = stats.nanmean(data, axis = 0)
1639 1642 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
1640 1643
1641 1644 phase = numpy.angle(data1[lag1,:])
1642 1645
1643 1646 return phase
1644 1647
1645 1648 class SpectralFitting(Operation):
1646 1649 '''
1647 1650 Function GetMoments()
1648 1651
1649 1652 Input:
1650 1653 Output:
1651 1654 Variables modified:
1652 1655 '''
1653 1656
1654 1657 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
1655 1658
1656 1659
1657 1660 if path != None:
1658 1661 sys.path.append(path)
1659 1662 self.dataOut.library = importlib.import_module(file)
1660 1663
1661 1664 #To be inserted as a parameter
1662 1665 groupArray = numpy.array(groupList)
1663 1666 # groupArray = numpy.array([[0,1],[2,3]])
1664 1667 self.dataOut.groupList = groupArray
1665 1668
1666 1669 nGroups = groupArray.shape[0]
1667 1670 nChannels = self.dataIn.nChannels
1668 1671 nHeights=self.dataIn.heightList.size
1669 1672
1670 1673 #Parameters Array
1671 1674 self.dataOut.data_param = None
1672 1675
1673 1676 #Set constants
1674 1677 constants = self.dataOut.library.setConstants(self.dataIn)
1675 1678 self.dataOut.constants = constants
1676 1679 M = self.dataIn.normFactor
1677 1680 N = self.dataIn.nFFTPoints
1678 1681 ippSeconds = self.dataIn.ippSeconds
1679 1682 K = self.dataIn.nIncohInt
1680 1683 pairsArray = numpy.array(self.dataIn.pairsList)
1681 1684
1682 1685 #List of possible combinations
1683 1686 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
1684 1687 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
1685 1688
1686 1689 if getSNR:
1687 1690 listChannels = groupArray.reshape((groupArray.size))
1688 1691 listChannels.sort()
1689 1692 noise = self.dataIn.getNoise()
1690 1693 self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
1691 1694
1692 1695 for i in range(nGroups):
1693 1696 coord = groupArray[i,:]
1694 1697
1695 1698 #Input data array
1696 1699 data = self.dataIn.data_spc[coord,:,:]/(M*N)
1697 1700 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
1698 1701
1699 1702 #Cross Spectra data array for Covariance Matrixes
1700 1703 ind = 0
1701 1704 for pairs in listComb:
1702 1705 pairsSel = numpy.array([coord[x],coord[y]])
1703 1706 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
1704 1707 ind += 1
1705 1708 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
1706 1709 dataCross = dataCross**2/K
1707 1710
1708 1711 for h in range(nHeights):
1709 1712 # print self.dataOut.heightList[h]
1710 1713
1711 1714 #Input
1712 1715 d = data[:,h]
1713 1716
1714 1717 #Covariance Matrix
1715 1718 D = numpy.diag(d**2/K)
1716 1719 ind = 0
1717 1720 for pairs in listComb:
1718 1721 #Coordinates in Covariance Matrix
1719 1722 x = pairs[0]
1720 1723 y = pairs[1]
1721 1724 #Channel Index
1722 1725 S12 = dataCross[ind,:,h]
1723 1726 D12 = numpy.diag(S12)
1724 1727 #Completing Covariance Matrix with Cross Spectras
1725 1728 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
1726 1729 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
1727 1730 ind += 1
1728 1731 Dinv=numpy.linalg.inv(D)
1729 1732 L=numpy.linalg.cholesky(Dinv)
1730 1733 LT=L.T
1731 1734
1732 1735 dp = numpy.dot(LT,d)
1733 1736
1734 1737 #Initial values
1735 1738 data_spc = self.dataIn.data_spc[coord,:,h]
1736 1739
1737 1740 if (h>0)and(error1[3]<5):
1738 1741 p0 = self.dataOut.data_param[i,:,h-1]
1739 1742 else:
1740 1743 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
1741 1744
1742 1745 try:
1743 1746 #Least Squares
1744 1747 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
1745 1748 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
1746 1749 #Chi square error
1747 1750 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
1748 1751 #Error with Jacobian
1749 1752 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
1750 1753 except:
1751 1754 minp = p0*numpy.nan
1752 1755 error0 = numpy.nan
1753 1756 error1 = p0*numpy.nan
1754 1757
1755 1758 #Save
1756 1759 if self.dataOut.data_param == None:
1757 1760 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
1758 1761 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
1759 1762
1760 1763 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
1761 1764 self.dataOut.data_param[i,:,h] = minp
1762 1765 return
1763 1766
1764 1767 def __residFunction(self, p, dp, LT, constants):
1765 1768
1766 1769 fm = self.dataOut.library.modelFunction(p, constants)
1767 1770 fmp=numpy.dot(LT,fm)
1768 1771
1769 1772 return dp-fmp
1770 1773
1771 1774 def __getSNR(self, z, noise):
1772 1775
1773 1776 avg = numpy.average(z, axis=1)
1774 1777 SNR = (avg.T-noise)/noise
1775 1778 SNR = SNR.T
1776 1779 return SNR
1777 1780
1778 1781 def __chisq(p,chindex,hindex):
1779 1782 #similar to Resid but calculates CHI**2
1780 1783 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
1781 1784 dp=numpy.dot(LT,d)
1782 1785 fmp=numpy.dot(LT,fm)
1783 1786 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
1784 1787 return chisq
1785 1788
1786 1789 class WindProfiler(Operation):
1787 1790
1788 1791 __isConfig = False
1789 1792
1790 1793 __initime = None
1791 1794 __lastdatatime = None
1792 1795 __integrationtime = None
1793 1796
1794 1797 __buffer = None
1795 1798
1796 1799 __dataReady = False
1797 1800
1798 1801 __firstdata = None
1799 1802
1800 1803 n = None
1801 1804
1802 1805 def __init__(self, **kwargs):
1803 1806 Operation.__init__(self, **kwargs)
1804 1807
1805 1808 def __calculateCosDir(self, elev, azim):
1806 1809 zen = (90 - elev)*numpy.pi/180
1807 1810 azim = azim*numpy.pi/180
1808 1811 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
1809 1812 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
1810 1813
1811 1814 signX = numpy.sign(numpy.cos(azim))
1812 1815 signY = numpy.sign(numpy.sin(azim))
1813 1816
1814 1817 cosDirX = numpy.copysign(cosDirX, signX)
1815 1818 cosDirY = numpy.copysign(cosDirY, signY)
1816 1819 return cosDirX, cosDirY
1817 1820
1818 1821 def __calculateAngles(self, theta_x, theta_y, azimuth):
1819 1822
1820 1823 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
1821 1824 zenith_arr = numpy.arccos(dir_cosw)
1822 1825 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
1823 1826
1824 1827 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
1825 1828 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
1826 1829
1827 1830 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
1828 1831
1829 1832 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
1830 1833
1831 1834 #
1832 1835 if horOnly:
1833 1836 A = numpy.c_[dir_cosu,dir_cosv]
1834 1837 else:
1835 1838 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
1836 1839 A = numpy.asmatrix(A)
1837 1840 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
1838 1841
1839 1842 return A1
1840 1843
1841 1844 def __correctValues(self, heiRang, phi, velRadial, SNR):
1842 1845 listPhi = phi.tolist()
1843 1846 maxid = listPhi.index(max(listPhi))
1844 1847 minid = listPhi.index(min(listPhi))
1845 1848
1846 1849 rango = range(len(phi))
1847 1850 # rango = numpy.delete(rango,maxid)
1848 1851
1849 1852 heiRang1 = heiRang*math.cos(phi[maxid])
1850 1853 heiRangAux = heiRang*math.cos(phi[minid])
1851 1854 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1852 1855 heiRang1 = numpy.delete(heiRang1,indOut)
1853 1856
1854 1857 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1855 1858 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1856 1859
1857 1860 for i in rango:
1858 1861 x = heiRang*math.cos(phi[i])
1859 1862 y1 = velRadial[i,:]
1860 1863 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1861 1864
1862 1865 x1 = heiRang1
1863 1866 y11 = f1(x1)
1864 1867
1865 1868 y2 = SNR[i,:]
1866 1869 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1867 1870 y21 = f2(x1)
1868 1871
1869 1872 velRadial1[i,:] = y11
1870 1873 SNR1[i,:] = y21
1871 1874
1872 1875 return heiRang1, velRadial1, SNR1
1873 1876
1874 1877 def __calculateVelUVW(self, A, velRadial):
1875 1878
1876 1879 #Operacion Matricial
1877 1880 # velUVW = numpy.zeros((velRadial.shape[1],3))
1878 1881 # for ind in range(velRadial.shape[1]):
1879 1882 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
1880 1883 # velUVW = velUVW.transpose()
1881 1884 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
1882 1885 velUVW[:,:] = numpy.dot(A,velRadial)
1883 1886
1884 1887
1885 1888 return velUVW
1886 1889
1887 1890 # def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
1888 1891
1889 1892 def techniqueDBS(self, kwargs):
1890 1893 """
1891 1894 Function that implements Doppler Beam Swinging (DBS) technique.
1892 1895
1893 1896 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1894 1897 Direction correction (if necessary), Ranges and SNR
1895 1898
1896 1899 Output: Winds estimation (Zonal, Meridional and Vertical)
1897 1900
1898 1901 Parameters affected: Winds, height range, SNR
1899 1902 """
1900 1903 velRadial0 = kwargs['velRadial']
1901 1904 heiRang = kwargs['heightList']
1902 1905 SNR0 = kwargs['SNR']
1903 1906
1904 1907 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
1905 1908 theta_x = numpy.array(kwargs['dirCosx'])
1906 1909 theta_y = numpy.array(kwargs['dirCosy'])
1907 1910 else:
1908 1911 elev = numpy.array(kwargs['elevation'])
1909 1912 azim = numpy.array(kwargs['azimuth'])
1910 1913 theta_x, theta_y = self.__calculateCosDir(elev, azim)
1911 1914 azimuth = kwargs['correctAzimuth']
1912 1915 if kwargs.has_key('horizontalOnly'):
1913 1916 horizontalOnly = kwargs['horizontalOnly']
1914 1917 else: horizontalOnly = False
1915 1918 if kwargs.has_key('correctFactor'):
1916 1919 correctFactor = kwargs['correctFactor']
1917 1920 else: correctFactor = 1
1918 1921 if kwargs.has_key('channelList'):
1919 1922 channelList = kwargs['channelList']
1920 1923 if len(channelList) == 2:
1921 1924 horizontalOnly = True
1922 1925 arrayChannel = numpy.array(channelList)
1923 1926 param = param[arrayChannel,:,:]
1924 1927 theta_x = theta_x[arrayChannel]
1925 1928 theta_y = theta_y[arrayChannel]
1926 1929
1927 1930 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
1928 1931 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
1929 1932 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
1930 1933
1931 1934 #Calculo de Componentes de la velocidad con DBS
1932 1935 winds = self.__calculateVelUVW(A,velRadial1)
1933 1936
1934 1937 return winds, heiRang1, SNR1
1935 1938
1936 1939 def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
1937 1940
1938 1941 nPairs = len(pairs_ccf)
1939 1942 posx = numpy.asarray(posx)
1940 1943 posy = numpy.asarray(posy)
1941 1944
1942 1945 #Rotacion Inversa para alinear con el azimuth
1943 1946 if azimuth!= None:
1944 1947 azimuth = azimuth*math.pi/180
1945 1948 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
1946 1949 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
1947 1950 else:
1948 1951 posx1 = posx
1949 1952 posy1 = posy
1950 1953
1951 1954 #Calculo de Distancias
1952 1955 distx = numpy.zeros(nPairs)
1953 1956 disty = numpy.zeros(nPairs)
1954 1957 dist = numpy.zeros(nPairs)
1955 1958 ang = numpy.zeros(nPairs)
1956 1959
1957 1960 for i in range(nPairs):
1958 1961 distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
1959 1962 disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
1960 1963 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
1961 1964 ang[i] = numpy.arctan2(disty[i],distx[i])
1962 1965
1963 1966 return distx, disty, dist, ang
1964 1967 #Calculo de Matrices
1965 1968 # nPairs = len(pairs)
1966 1969 # ang1 = numpy.zeros((nPairs, 2, 1))
1967 1970 # dist1 = numpy.zeros((nPairs, 2, 1))
1968 1971 #
1969 1972 # for j in range(nPairs):
1970 1973 # dist1[j,0,0] = dist[pairs[j][0]]
1971 1974 # dist1[j,1,0] = dist[pairs[j][1]]
1972 1975 # ang1[j,0,0] = ang[pairs[j][0]]
1973 1976 # ang1[j,1,0] = ang[pairs[j][1]]
1974 1977 #
1975 1978 # return distx,disty, dist1,ang1
1976 1979
1977 1980
1978 1981 def __calculateVelVer(self, phase, lagTRange, _lambda):
1979 1982
1980 1983 Ts = lagTRange[1] - lagTRange[0]
1981 1984 velW = -_lambda*phase/(4*math.pi*Ts)
1982 1985
1983 1986 return velW
1984 1987
1985 1988 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
1986 1989 nPairs = tau1.shape[0]
1987 1990 nHeights = tau1.shape[1]
1988 1991 vel = numpy.zeros((nPairs,3,nHeights))
1989 1992 dist1 = numpy.reshape(dist, (dist.size,1))
1990 1993
1991 1994 angCos = numpy.cos(ang)
1992 1995 angSin = numpy.sin(ang)
1993 1996
1994 1997 vel0 = dist1*tau1/(2*tau2**2)
1995 1998 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
1996 1999 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
1997 2000
1998 2001 ind = numpy.where(numpy.isinf(vel))
1999 2002 vel[ind] = numpy.nan
2000 2003
2001 2004 return vel
2002 2005
2003 2006 # def __getPairsAutoCorr(self, pairsList, nChannels):
2004 2007 #
2005 2008 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
2006 2009 #
2007 2010 # for l in range(len(pairsList)):
2008 2011 # firstChannel = pairsList[l][0]
2009 2012 # secondChannel = pairsList[l][1]
2010 2013 #
2011 2014 # #Obteniendo pares de Autocorrelacion
2012 2015 # if firstChannel == secondChannel:
2013 2016 # pairsAutoCorr[firstChannel] = int(l)
2014 2017 #
2015 2018 # pairsAutoCorr = pairsAutoCorr.astype(int)
2016 2019 #
2017 2020 # pairsCrossCorr = range(len(pairsList))
2018 2021 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
2019 2022 #
2020 2023 # return pairsAutoCorr, pairsCrossCorr
2021 2024
2022 2025 # def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
2023 2026 def techniqueSA(self, kwargs):
2024 2027
2025 2028 """
2026 2029 Function that implements Spaced Antenna (SA) technique.
2027 2030
2028 2031 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
2029 2032 Direction correction (if necessary), Ranges and SNR
2030 2033
2031 2034 Output: Winds estimation (Zonal, Meridional and Vertical)
2032 2035
2033 2036 Parameters affected: Winds
2034 2037 """
2035 2038 position_x = kwargs['positionX']
2036 2039 position_y = kwargs['positionY']
2037 2040 azimuth = kwargs['azimuth']
2038 2041
2039 2042 if kwargs.has_key('correctFactor'):
2040 2043 correctFactor = kwargs['correctFactor']
2041 2044 else:
2042 2045 correctFactor = 1
2043 2046
2044 2047 groupList = kwargs['groupList']
2045 2048 pairs_ccf = groupList[1]
2046 2049 tau = kwargs['tau']
2047 2050 _lambda = kwargs['_lambda']
2048 2051
2049 2052 #Cross Correlation pairs obtained
2050 2053 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels)
2051 2054 # pairsArray = numpy.array(pairsList)[pairsCrossCorr]
2052 2055 # pairsSelArray = numpy.array(pairsSelected)
2053 2056 # pairs = []
2054 2057 #
2055 2058 # #Wind estimation pairs obtained
2056 2059 # for i in range(pairsSelArray.shape[0]/2):
2057 2060 # ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
2058 2061 # ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
2059 2062 # pairs.append((ind1,ind2))
2060 2063
2061 2064 indtau = tau.shape[0]/2
2062 2065 tau1 = tau[:indtau,:]
2063 2066 tau2 = tau[indtau:-1,:]
2064 2067 # tau1 = tau1[pairs,:]
2065 2068 # tau2 = tau2[pairs,:]
2066 2069 phase1 = tau[-1,:]
2067 2070
2068 2071 #---------------------------------------------------------------------
2069 2072 #Metodo Directo
2070 2073 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
2071 2074 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
2072 2075 winds = stats.nanmean(winds, axis=0)
2073 2076 #---------------------------------------------------------------------
2074 2077 #Metodo General
2075 2078 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
2076 2079 # #Calculo Coeficientes de Funcion de Correlacion
2077 2080 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
2078 2081 # #Calculo de Velocidades
2079 2082 # winds = self.calculateVelUV(F,G,A,B,H)
2080 2083
2081 2084 #---------------------------------------------------------------------
2082 2085 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
2083 2086 winds = correctFactor*winds
2084 2087 return winds
2085 2088
2086 2089 def __checkTime(self, currentTime, paramInterval, outputInterval):
2087 2090
2088 2091 dataTime = currentTime + paramInterval
2089 2092 deltaTime = dataTime - self.__initime
2090 2093
2091 2094 if deltaTime >= outputInterval or deltaTime < 0:
2092 2095 self.__dataReady = True
2093 2096 return
2094 2097
2095 2098 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
2096 2099 '''
2097 2100 Function that implements winds estimation technique with detected meteors.
2098 2101
2099 2102 Input: Detected meteors, Minimum meteor quantity to wind estimation
2100 2103
2101 2104 Output: Winds estimation (Zonal and Meridional)
2102 2105
2103 2106 Parameters affected: Winds
2104 2107 '''
2105 2108 # print arrayMeteor.shape
2106 2109 #Settings
2107 2110 nInt = (heightMax - heightMin)/2
2108 2111 # print nInt
2109 2112 nInt = int(nInt)
2110 2113 # print nInt
2111 2114 winds = numpy.zeros((2,nInt))*numpy.nan
2112 2115
2113 2116 #Filter errors
2114 2117 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
2115 2118 finalMeteor = arrayMeteor[error,:]
2116 2119
2117 2120 #Meteor Histogram
2118 2121 finalHeights = finalMeteor[:,2]
2119 2122 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
2120 2123 nMeteorsPerI = hist[0]
2121 2124 heightPerI = hist[1]
2122 2125
2123 2126 #Sort of meteors
2124 2127 indSort = finalHeights.argsort()
2125 2128 finalMeteor2 = finalMeteor[indSort,:]
2126 2129
2127 2130 # Calculating winds
2128 2131 ind1 = 0
2129 2132 ind2 = 0
2130 2133
2131 2134 for i in range(nInt):
2132 2135 nMet = nMeteorsPerI[i]
2133 2136 ind1 = ind2
2134 2137 ind2 = ind1 + nMet
2135 2138
2136 2139 meteorAux = finalMeteor2[ind1:ind2,:]
2137 2140
2138 2141 if meteorAux.shape[0] >= meteorThresh:
2139 2142 vel = meteorAux[:, 6]
2140 2143 zen = meteorAux[:, 4]*numpy.pi/180
2141 2144 azim = meteorAux[:, 3]*numpy.pi/180
2142 2145
2143 2146 n = numpy.cos(zen)
2144 2147 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
2145 2148 # l = m*numpy.tan(azim)
2146 2149 l = numpy.sin(zen)*numpy.sin(azim)
2147 2150 m = numpy.sin(zen)*numpy.cos(azim)
2148 2151
2149 2152 A = numpy.vstack((l, m)).transpose()
2150 2153 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
2151 2154 windsAux = numpy.dot(A1, vel)
2152 2155
2153 2156 winds[0,i] = windsAux[0]
2154 2157 winds[1,i] = windsAux[1]
2155 2158
2156 2159 return winds, heightPerI[:-1]
2157 2160
2158 2161 def techniqueNSM_SA(self, **kwargs):
2159 2162 metArray = kwargs['metArray']
2160 2163 heightList = kwargs['heightList']
2161 2164 timeList = kwargs['timeList']
2162 2165
2163 2166 rx_location = kwargs['rx_location']
2164 2167 groupList = kwargs['groupList']
2165 2168 azimuth = kwargs['azimuth']
2166 2169 dfactor = kwargs['dfactor']
2167 2170 k = kwargs['k']
2168 2171
2169 2172 azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
2170 2173 d = dist*dfactor
2171 2174 #Phase calculation
2172 2175 metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
2173 2176
2174 2177 metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
2175 2178
2176 2179 velEst = numpy.zeros((heightList.size,2))*numpy.nan
2177 2180 azimuth1 = azimuth1*numpy.pi/180
2178 2181
2179 2182 for i in range(heightList.size):
2180 2183 h = heightList[i]
2181 2184 indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
2182 2185 metHeight = metArray1[indH,:]
2183 2186 if metHeight.shape[0] >= 2:
2184 2187 velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities
2185 2188 iazim = metHeight[:,1].astype(int)
2186 2189 azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths
2187 2190 A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
2188 2191 A = numpy.asmatrix(A)
2189 2192 A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()
2190 2193 velHor = numpy.dot(A1,velAux)
2191 2194
2192 2195 velEst[i,:] = numpy.squeeze(velHor)
2193 2196 return velEst
2194 2197
2195 2198 def __getPhaseSlope(self, metArray, heightList, timeList):
2196 2199 meteorList = []
2197 2200 #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
2198 2201 #Putting back together the meteor matrix
2199 2202 utctime = metArray[:,0]
2200 2203 uniqueTime = numpy.unique(utctime)
2201 2204
2202 2205 phaseDerThresh = 0.5
2203 2206 ippSeconds = timeList[1] - timeList[0]
2204 2207 sec = numpy.where(timeList>1)[0][0]
2205 2208 nPairs = metArray.shape[1] - 6
2206 2209 nHeights = len(heightList)
2207 2210
2208 2211 for t in uniqueTime:
2209 2212 metArray1 = metArray[utctime==t,:]
2210 2213 # phaseDerThresh = numpy.pi/4 #reducir Phase thresh
2211 2214 tmet = metArray1[:,1].astype(int)
2212 2215 hmet = metArray1[:,2].astype(int)
2213 2216
2214 2217 metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
2215 2218 metPhase[:,:] = numpy.nan
2216 2219 metPhase[:,hmet,tmet] = metArray1[:,6:].T
2217 2220
2218 2221 #Delete short trails
2219 2222 metBool = ~numpy.isnan(metPhase[0,:,:])
2220 2223 heightVect = numpy.sum(metBool, axis = 1)
2221 2224 metBool[heightVect<sec,:] = False
2222 2225 metPhase[:,heightVect<sec,:] = numpy.nan
2223 2226
2224 2227 #Derivative
2225 2228 metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
2226 2229 phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
2227 2230 metPhase[phDerAux] = numpy.nan
2228 2231
2229 2232 #--------------------------METEOR DETECTION -----------------------------------------
2230 2233 indMet = numpy.where(numpy.any(metBool,axis=1))[0]
2231 2234
2232 2235 for p in numpy.arange(nPairs):
2233 2236 phase = metPhase[p,:,:]
2234 2237 phDer = metDer[p,:,:]
2235 2238
2236 2239 for h in indMet:
2237 2240 height = heightList[h]
2238 2241 phase1 = phase[h,:] #82
2239 2242 phDer1 = phDer[h,:]
2240 2243
2241 2244 phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap
2242 2245
2243 2246 indValid = numpy.where(~numpy.isnan(phase1))[0]
2244 2247 initMet = indValid[0]
2245 2248 endMet = 0
2246 2249
2247 2250 for i in range(len(indValid)-1):
2248 2251
2249 2252 #Time difference
2250 2253 inow = indValid[i]
2251 2254 inext = indValid[i+1]
2252 2255 idiff = inext - inow
2253 2256 #Phase difference
2254 2257 phDiff = numpy.abs(phase1[inext] - phase1[inow])
2255 2258
2256 2259 if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor
2257 2260 sizeTrail = inow - initMet + 1
2258 2261 if sizeTrail>3*sec: #Too short meteors
2259 2262 x = numpy.arange(initMet,inow+1)*ippSeconds
2260 2263 y = phase1[initMet:inow+1]
2261 2264 ynnan = ~numpy.isnan(y)
2262 2265 x = x[ynnan]
2263 2266 y = y[ynnan]
2264 2267 slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
2265 2268 ylin = x*slope + intercept
2266 2269 rsq = r_value**2
2267 2270 if rsq > 0.5:
2268 2271 vel = slope#*height*1000/(k*d)
2269 2272 estAux = numpy.array([utctime,p,height, vel, rsq])
2270 2273 meteorList.append(estAux)
2271 2274 initMet = inext
2272 2275 metArray2 = numpy.array(meteorList)
2273 2276
2274 2277 return metArray2
2275 2278
2276 2279 def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
2277 2280
2278 2281 azimuth1 = numpy.zeros(len(pairslist))
2279 2282 dist = numpy.zeros(len(pairslist))
2280 2283
2281 2284 for i in range(len(rx_location)):
2282 2285 ch0 = pairslist[i][0]
2283 2286 ch1 = pairslist[i][1]
2284 2287
2285 2288 diffX = rx_location[ch0][0] - rx_location[ch1][0]
2286 2289 diffY = rx_location[ch0][1] - rx_location[ch1][1]
2287 2290 azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
2288 2291 dist[i] = numpy.sqrt(diffX**2 + diffY**2)
2289 2292
2290 2293 azimuth1 -= azimuth0
2291 2294 return azimuth1, dist
2292 2295
2293 2296 def techniqueNSM_DBS(self, **kwargs):
2294 2297 metArray = kwargs['metArray']
2295 2298 heightList = kwargs['heightList']
2296 2299 timeList = kwargs['timeList']
2297 2300 zenithList = kwargs['zenithList']
2298 2301 nChan = numpy.max(cmet) + 1
2299 2302 nHeights = len(heightList)
2300 2303
2301 2304 utctime = metArray[:,0]
2302 2305 cmet = metArray[:,1]
2303 2306 hmet = metArray1[:,3].astype(int)
2304 2307 h1met = heightList[hmet]*zenithList[cmet]
2305 2308 vmet = metArray1[:,5]
2306 2309
2307 2310 for i in range(nHeights - 1):
2308 2311 hmin = heightList[i]
2309 2312 hmax = heightList[i + 1]
2310 2313
2311 2314 vthisH = vmet[(h1met>=hmin) & (h1met<hmax)]
2312 2315
2313 2316
2314 2317
2315 2318 return data_output
2316 2319
2317 2320 def run(self, dataOut, technique, positionY, positionX, azimuth, **kwargs):
2318 2321
2319 2322 param = dataOut.data_param
2320 2323 if dataOut.abscissaList != None:
2321 2324 absc = dataOut.abscissaList[:-1]
2322 2325 noise = dataOut.noise
2323 2326 heightList = dataOut.heightList
2324 2327 SNR = dataOut.data_SNR
2325 2328
2326 2329 if technique == 'DBS':
2327 2330
2328 2331 kwargs['velRadial'] = param[:,1,:] #Radial velocity
2329 2332 kwargs['heightList'] = heightList
2330 2333 kwargs['SNR'] = SNR
2331 2334
2332 2335 dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(kwargs) #DBS Function
2333 2336 dataOut.utctimeInit = dataOut.utctime
2334 2337 dataOut.outputInterval = dataOut.paramInterval
2335 2338
2336 2339 elif technique == 'SA':
2337 2340
2338 2341 #Parameters
2339 2342 # position_x = kwargs['positionX']
2340 2343 # position_y = kwargs['positionY']
2341 2344 # azimuth = kwargs['azimuth']
2342 2345 #
2343 2346 # if kwargs.has_key('crosspairsList'):
2344 2347 # pairs = kwargs['crosspairsList']
2345 2348 # else:
2346 2349 # pairs = None
2347 2350 #
2348 2351 # if kwargs.has_key('correctFactor'):
2349 2352 # correctFactor = kwargs['correctFactor']
2350 2353 # else:
2351 2354 # correctFactor = 1
2352 2355
2353 2356 # tau = dataOut.data_param
2354 2357 # _lambda = dataOut.C/dataOut.frequency
2355 2358 # pairsList = dataOut.groupList
2356 2359 # nChannels = dataOut.nChannels
2357 2360
2358 2361 kwargs['groupList'] = dataOut.groupList
2359 2362 kwargs['tau'] = dataOut.data_param
2360 2363 kwargs['_lambda'] = dataOut.C/dataOut.frequency
2361 2364 # dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
2362 2365 dataOut.data_output = self.techniqueSA(kwargs)
2363 2366 dataOut.utctimeInit = dataOut.utctime
2364 2367 dataOut.outputInterval = dataOut.timeInterval
2365 2368
2366 2369 elif technique == 'Meteors':
2367 2370 dataOut.flagNoData = True
2368 2371 self.__dataReady = False
2369 2372
2370 2373 if kwargs.has_key('nHours'):
2371 2374 nHours = kwargs['nHours']
2372 2375 else:
2373 2376 nHours = 1
2374 2377
2375 2378 if kwargs.has_key('meteorsPerBin'):
2376 2379 meteorThresh = kwargs['meteorsPerBin']
2377 2380 else:
2378 2381 meteorThresh = 6
2379 2382
2380 2383 if kwargs.has_key('hmin'):
2381 2384 hmin = kwargs['hmin']
2382 2385 else: hmin = 70
2383 2386 if kwargs.has_key('hmax'):
2384 2387 hmax = kwargs['hmax']
2385 2388 else: hmax = 110
2386 2389
2387 2390 dataOut.outputInterval = nHours*3600
2388 2391
2389 2392 if self.__isConfig == False:
2390 2393 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
2391 2394 #Get Initial LTC time
2392 2395 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
2393 2396 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
2394 2397
2395 2398 self.__isConfig = True
2396 2399
2397 2400 if self.__buffer == None:
2398 2401 self.__buffer = dataOut.data_param
2399 2402 self.__firstdata = copy.copy(dataOut)
2400 2403
2401 2404 else:
2402 2405 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
2403 2406
2404 2407 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
2405 2408
2406 2409 if self.__dataReady:
2407 2410 dataOut.utctimeInit = self.__initime
2408 2411
2409 2412 self.__initime += dataOut.outputInterval #to erase time offset
2410 2413
2411 2414 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
2412 2415 dataOut.flagNoData = False
2413 2416 self.__buffer = None
2414 2417
2415 2418 elif technique == 'Meteors1':
2416 2419 dataOut.flagNoData = True
2417 2420 self.__dataReady = False
2418 2421
2419 2422 if kwargs.has_key('nMins'):
2420 2423 nMins = kwargs['nMins']
2421 2424 else: nMins = 20
2422 2425 if kwargs.has_key('rx_location'):
2423 2426 rx_location = kwargs['rx_location']
2424 2427 else: rx_location = [(0,1),(1,1),(1,0)]
2425 2428 if kwargs.has_key('azimuth'):
2426 2429 azimuth = kwargs['azimuth']
2427 2430 else: azimuth = 51
2428 2431 if kwargs.has_key('dfactor'):
2429 2432 dfactor = kwargs['dfactor']
2430 2433 if kwargs.has_key('mode'):
2431 2434 mode = kwargs['mode']
2432 2435 else: mode = 'SA'
2433 2436
2434 2437 #Borrar luego esto
2435 2438 if dataOut.groupList == None:
2436 2439 dataOut.groupList = [(0,1),(0,2),(1,2)]
2437 2440 groupList = dataOut.groupList
2438 2441 C = 3e8
2439 2442 freq = 50e6
2440 2443 lamb = C/freq
2441 2444 k = 2*numpy.pi/lamb
2442 2445
2443 2446 timeList = dataOut.abscissaList
2444 2447 heightList = dataOut.heightList
2445 2448
2446 2449 if self.__isConfig == False:
2447 2450 dataOut.outputInterval = nMins*60
2448 2451 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
2449 2452 #Get Initial LTC time
2450 2453 initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
2451 2454 minuteAux = initime.minute
2452 2455 minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
2453 2456 self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
2454 2457
2455 2458 self.__isConfig = True
2456 2459
2457 2460 if self.__buffer == None:
2458 2461 self.__buffer = dataOut.data_param
2459 2462 self.__firstdata = copy.copy(dataOut)
2460 2463
2461 2464 else:
2462 2465 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
2463 2466
2464 2467 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
2465 2468
2466 2469 if self.__dataReady:
2467 2470 dataOut.utctimeInit = self.__initime
2468 2471 self.__initime += dataOut.outputInterval #to erase time offset
2469 2472
2470 2473 metArray = self.__buffer
2471 2474 if mode == 'SA':
2472 2475 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
2473 2476 elif mode == 'DBS':
2474 2477 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList)
2475 2478 dataOut.data_output = dataOut.data_output.T
2476 2479 dataOut.flagNoData = False
2477 2480 self.__buffer = None
2478 2481
2479 2482 return
2480 2483
2481 2484 class EWDriftsEstimation(Operation):
2482 2485
2483 2486 def __init__(self):
2484 2487 Operation.__init__(self)
2485 2488
2486 2489 def __correctValues(self, heiRang, phi, velRadial, SNR):
2487 2490 listPhi = phi.tolist()
2488 2491 maxid = listPhi.index(max(listPhi))
2489 2492 minid = listPhi.index(min(listPhi))
2490 2493
2491 2494 rango = range(len(phi))
2492 2495 # rango = numpy.delete(rango,maxid)
2493 2496
2494 2497 heiRang1 = heiRang*math.cos(phi[maxid])
2495 2498 heiRangAux = heiRang*math.cos(phi[minid])
2496 2499 indOut = (heiRang1 < heiRangAux[0]).nonzero()
2497 2500 heiRang1 = numpy.delete(heiRang1,indOut)
2498 2501
2499 2502 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
2500 2503 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
2501 2504
2502 2505 for i in rango:
2503 2506 x = heiRang*math.cos(phi[i])
2504 2507 y1 = velRadial[i,:]
2505 2508 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
2506 2509
2507 2510 x1 = heiRang1
2508 2511 y11 = f1(x1)
2509 2512
2510 2513 y2 = SNR[i,:]
2511 2514 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
2512 2515 y21 = f2(x1)
2513 2516
2514 2517 velRadial1[i,:] = y11
2515 2518 SNR1[i,:] = y21
2516 2519
2517 2520 return heiRang1, velRadial1, SNR1
2518 2521
2519 2522 def run(self, dataOut, zenith, zenithCorrection):
2520 2523 heiRang = dataOut.heightList
2521 2524 velRadial = dataOut.data_param[:,3,:]
2522 2525 SNR = dataOut.data_SNR
2523 2526
2524 2527 zenith = numpy.array(zenith)
2525 2528 zenith -= zenithCorrection
2526 2529 zenith *= numpy.pi/180
2527 2530
2528 2531 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
2529 2532
2530 2533 alp = zenith[0]
2531 2534 bet = zenith[1]
2532 2535
2533 2536 w_w = velRadial1[0,:]
2534 2537 w_e = velRadial1[1,:]
2535 2538
2536 2539 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
2537 2540 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
2538 2541
2539 2542 winds = numpy.vstack((u,w))
2540 2543
2541 2544 dataOut.heightList = heiRang1
2542 2545 dataOut.data_output = winds
2543 2546 dataOut.data_SNR = SNR1
2544 2547
2545 2548 dataOut.utctimeInit = dataOut.utctime
2546 2549 dataOut.outputInterval = dataOut.timeInterval
2547 2550 return
2548 2551
2549 2552 #--------------- Non Specular Meteor ----------------
2550 2553
2551 2554 class NonSpecularMeteorDetection(Operation):
2552 2555
2553 2556 def run(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):
2554 2557 data_acf = self.dataOut.data_pre[0]
2555 2558 data_ccf = self.dataOut.data_pre[1]
2556 2559
2557 2560 lamb = self.dataOut.C/self.dataOut.frequency
2558 2561 tSamp = self.dataOut.ippSeconds*self.dataOut.nCohInt
2559 2562 paramInterval = self.dataOut.paramInterval
2560 2563
2561 2564 nChannels = data_acf.shape[0]
2562 2565 nLags = data_acf.shape[1]
2563 2566 nProfiles = data_acf.shape[2]
2564 2567 nHeights = self.dataOut.nHeights
2565 2568 nCohInt = self.dataOut.nCohInt
2566 2569 sec = numpy.round(nProfiles/self.dataOut.paramInterval)
2567 2570 heightList = self.dataOut.heightList
2568 2571 ippSeconds = self.dataOut.ippSeconds*self.dataOut.nCohInt*self.dataOut.nAvg
2569 2572 utctime = self.dataOut.utctime
2570 2573
2571 2574 self.dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds)
2572 2575
2573 2576 #------------------------ SNR --------------------------------------
2574 2577 power = data_acf[:,0,:,:].real
2575 2578 noise = numpy.zeros(nChannels)
2576 2579 SNR = numpy.zeros(power.shape)
2577 2580 for i in range(nChannels):
2578 2581 noise[i] = hildebrand_sekhon(power[i,:], nCohInt)
2579 2582 SNR[i] = (power[i]-noise[i])/noise[i]
2580 2583 SNRm = numpy.nanmean(SNR, axis = 0)
2581 2584 SNRdB = 10*numpy.log10(SNR)
2582 2585
2583 2586 if mode == 'SA':
2584 2587 nPairs = data_ccf.shape[0]
2585 2588 #---------------------- Coherence and Phase --------------------------
2586 2589 phase = numpy.zeros(data_ccf[:,0,:,:].shape)
2587 2590 # phase1 = numpy.copy(phase)
2588 2591 coh1 = numpy.zeros(data_ccf[:,0,:,:].shape)
2589 2592
2590 2593 for p in range(nPairs):
2591 2594 ch0 = self.dataOut.groupList[p][0]
2592 2595 ch1 = self.dataOut.groupList[p][1]
2593 2596 ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:])
2594 2597 phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter
2595 2598 # phase1[p,:,:] = numpy.angle(ccf) #median filter
2596 2599 coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter
2597 2600 # coh1[p,:,:] = numpy.abs(ccf) #median filter
2598 2601 coh = numpy.nanmax(coh1, axis = 0)
2599 2602 # struc = numpy.ones((5,1))
2600 2603 # coh = ndimage.morphology.grey_dilation(coh, size=(10,1))
2601 2604 #---------------------- Radial Velocity ----------------------------
2602 2605 phaseAux = numpy.mean(numpy.angle(data_acf[:,1,:,:]), axis = 0)
2603 2606 velRad = phaseAux*lamb/(4*numpy.pi*tSamp)
2604 2607
2605 2608 if allData:
2606 2609 boolMetFin = ~numpy.isnan(SNRm)
2607 2610 # coh[:-1,:] = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
2608 2611 else:
2609 2612 #------------------------ Meteor mask ---------------------------------
2610 2613 # #SNR mask
2611 2614 # boolMet = (SNRdB>SNRthresh)#|(~numpy.isnan(SNRdB))
2612 2615 #
2613 2616 # #Erase small objects
2614 2617 # boolMet1 = self.__erase_small(boolMet, 2*sec, 5)
2615 2618 #
2616 2619 # auxEEJ = numpy.sum(boolMet1,axis=0)
2617 2620 # indOver = auxEEJ>nProfiles*0.8 #Use this later
2618 2621 # indEEJ = numpy.where(indOver)[0]
2619 2622 # indNEEJ = numpy.where(~indOver)[0]
2620 2623 #
2621 2624 # boolMetFin = boolMet1
2622 2625 #
2623 2626 # if indEEJ.size > 0:
2624 2627 # boolMet1[:,indEEJ] = False #Erase heights with EEJ
2625 2628 #
2626 2629 # boolMet2 = coh > cohThresh
2627 2630 # boolMet2 = self.__erase_small(boolMet2, 2*sec,5)
2628 2631 #
2629 2632 # #Final Meteor mask
2630 2633 # boolMetFin = boolMet1|boolMet2
2631 2634
2632 2635 #Coherence mask
2633 2636 boolMet1 = coh > 0.75
2634 2637 struc = numpy.ones((30,1))
2635 2638 boolMet1 = ndimage.morphology.binary_dilation(boolMet1, structure=struc)
2636 2639
2637 2640 #Derivative mask
2638 2641 derPhase = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
2639 2642 boolMet2 = derPhase < 0.2
2640 2643 # boolMet2 = ndimage.morphology.binary_opening(boolMet2)
2641 2644 # boolMet2 = ndimage.morphology.binary_closing(boolMet2, structure = numpy.ones((10,1)))
2642 2645 boolMet2 = ndimage.median_filter(boolMet2,size=5)
2643 2646 boolMet2 = numpy.vstack((boolMet2,numpy.full((1,nHeights), True, dtype=bool)))
2644 2647 # #Final mask
2645 2648 # boolMetFin = boolMet2
2646 2649 boolMetFin = boolMet1&boolMet2
2647 2650 # boolMetFin = ndimage.morphology.binary_dilation(boolMetFin)
2648 2651 #Creating data_param
2649 2652 coordMet = numpy.where(boolMetFin)
2650 2653
2651 2654 tmet = coordMet[0]
2652 2655 hmet = coordMet[1]
2653 2656
2654 2657 data_param = numpy.zeros((tmet.size, 6 + nPairs))
2655 2658 data_param[:,0] = utctime
2656 2659 data_param[:,1] = tmet
2657 2660 data_param[:,2] = hmet
2658 2661 data_param[:,3] = SNRm[tmet,hmet]
2659 2662 data_param[:,4] = velRad[tmet,hmet]
2660 2663 data_param[:,5] = coh[tmet,hmet]
2661 2664 data_param[:,6:] = phase[:,tmet,hmet].T
2662 2665
2663 2666 elif mode == 'DBS':
2664 2667 self.dataOut.groupList = numpy.arange(nChannels)
2665 2668
2666 2669 #Radial Velocities
2667 2670 # phase = numpy.angle(data_acf[:,1,:,:])
2668 2671 phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1))
2669 2672 velRad = phase*lamb/(4*numpy.pi*tSamp)
2670 2673
2671 2674 #Spectral width
2672 2675 acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1))
2673 2676 acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1))
2674 2677
2675 2678 spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2))
2676 2679 # velRad = ndimage.median_filter(velRad, size = (1,5,1))
2677 2680 if allData:
2678 2681 boolMetFin = ~numpy.isnan(SNRdB)
2679 2682 else:
2680 2683 #SNR
2681 2684 boolMet1 = (SNRdB>SNRthresh) #SNR mask
2682 2685 boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5))
2683 2686
2684 2687 #Radial velocity
2685 2688 boolMet2 = numpy.abs(velRad) < 30
2686 2689 boolMet2 = ndimage.median_filter(boolMet2, (1,5,5))
2687 2690
2688 2691 #Spectral Width
2689 2692 boolMet3 = spcWidth < 30
2690 2693 boolMet3 = ndimage.median_filter(boolMet3, (1,5,5))
2691 2694 # boolMetFin = self.__erase_small(boolMet1, 10,5)
2692 2695 boolMetFin = boolMet1&boolMet2&boolMet3
2693 2696
2694 2697 #Creating data_param
2695 2698 coordMet = numpy.where(boolMetFin)
2696 2699
2697 2700 cmet = coordMet[0]
2698 2701 tmet = coordMet[1]
2699 2702 hmet = coordMet[2]
2700 2703
2701 2704 data_param = numpy.zeros((tmet.size, 7))
2702 2705 data_param[:,0] = utctime
2703 2706 data_param[:,1] = cmet
2704 2707 data_param[:,2] = tmet
2705 2708 data_param[:,3] = hmet
2706 2709 data_param[:,4] = SNR[cmet,tmet,hmet].T
2707 2710 data_param[:,5] = velRad[cmet,tmet,hmet].T
2708 2711 data_param[:,6] = spcWidth[cmet,tmet,hmet].T
2709 2712
2710 2713 # self.dataOut.data_param = data_int
2711 2714 if len(data_param) == 0:
2712 2715 self.dataOut.flagNoData = True
2713 2716 else:
2714 2717 self.dataOut.data_param = data_param
2715 2718
2716 2719 def __erase_small(self, binArray, threshX, threshY):
2717 2720 labarray, numfeat = ndimage.measurements.label(binArray)
2718 2721 binArray1 = numpy.copy(binArray)
2719 2722
2720 2723 for i in range(1,numfeat + 1):
2721 2724 auxBin = (labarray==i)
2722 2725 auxSize = auxBin.sum()
2723 2726
2724 2727 x,y = numpy.where(auxBin)
2725 2728 widthX = x.max() - x.min()
2726 2729 widthY = y.max() - y.min()
2727 2730
2728 2731 #width X: 3 seg -> 12.5*3
2729 2732 #width Y:
2730 2733
2731 2734 if (auxSize < 50) or (widthX < threshX) or (widthY < threshY):
2732 2735 binArray1[auxBin] = False
2733 2736
2734 2737 return binArray1
2735 2738
2736 2739 #--------------- Specular Meteor ----------------
2737 2740
2738 2741 class SMDetection(Operation):
2739 2742 '''
2740 2743 Function DetectMeteors()
2741 2744 Project developed with paper:
2742 2745 HOLDSWORTH ET AL. 2004
2743 2746
2744 2747 Input:
2745 2748 self.dataOut.data_pre
2746 2749
2747 2750 centerReceiverIndex: From the channels, which is the center receiver
2748 2751
2749 2752 hei_ref: Height reference for the Beacon signal extraction
2750 2753 tauindex:
2751 2754 predefinedPhaseShifts: Predefined phase offset for the voltge signals
2752 2755
2753 2756 cohDetection: Whether to user Coherent detection or not
2754 2757 cohDet_timeStep: Coherent Detection calculation time step
2755 2758 cohDet_thresh: Coherent Detection phase threshold to correct phases
2756 2759
2757 2760 noise_timeStep: Noise calculation time step
2758 2761 noise_multiple: Noise multiple to define signal threshold
2759 2762
2760 2763 multDet_timeLimit: Multiple Detection Removal time limit in seconds
2761 2764 multDet_rangeLimit: Multiple Detection Removal range limit in km
2762 2765
2763 2766 phaseThresh: Maximum phase difference between receiver to be consider a meteor
2764 2767 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
2765 2768
2766 2769 hmin: Minimum Height of the meteor to use it in the further wind estimations
2767 2770 hmax: Maximum Height of the meteor to use it in the further wind estimations
2768 2771 azimuth: Azimuth angle correction
2769 2772
2770 2773 Affected:
2771 2774 self.dataOut.data_param
2772 2775
2773 2776 Rejection Criteria (Errors):
2774 2777 0: No error; analysis OK
2775 2778 1: SNR < SNR threshold
2776 2779 2: angle of arrival (AOA) ambiguously determined
2777 2780 3: AOA estimate not feasible
2778 2781 4: Large difference in AOAs obtained from different antenna baselines
2779 2782 5: echo at start or end of time series
2780 2783 6: echo less than 5 examples long; too short for analysis
2781 2784 7: echo rise exceeds 0.3s
2782 2785 8: echo decay time less than twice rise time
2783 2786 9: large power level before echo
2784 2787 10: large power level after echo
2785 2788 11: poor fit to amplitude for estimation of decay time
2786 2789 12: poor fit to CCF phase variation for estimation of radial drift velocity
2787 2790 13: height unresolvable echo: not valid height within 70 to 110 km
2788 2791 14: height ambiguous echo: more then one possible height within 70 to 110 km
2789 2792 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
2790 2793 16: oscilatory echo, indicating event most likely not an underdense echo
2791 2794
2792 2795 17: phase difference in meteor Reestimation
2793 2796
2794 2797 Data Storage:
2795 2798 Meteors for Wind Estimation (8):
2796 2799 Utc Time | Range Height
2797 2800 Azimuth Zenith errorCosDir
2798 2801 VelRad errorVelRad
2799 2802 Phase0 Phase1 Phase2 Phase3
2800 2803 TypeError
2801 2804
2802 2805 '''
2803 2806
2804 2807 def run(self, dataOut, hei_ref = None, tauindex = 0,
2805 2808 phaseOffsets = None,
2806 2809 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
2807 2810 noise_timeStep = 4, noise_multiple = 4,
2808 2811 multDet_timeLimit = 1, multDet_rangeLimit = 3,
2809 2812 phaseThresh = 20, SNRThresh = 5,
2810 2813 hmin = 50, hmax=150, azimuth = 0,
2811 2814 channelPositions = None) :
2812 2815
2813 2816
2814 2817 #Getting Pairslist
2815 2818 if channelPositions == None:
2816 2819 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2817 2820 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2818 2821 meteorOps = SMOperations()
2819 2822 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2820 2823 heiRang = dataOut.getHeiRange()
2821 2824 #Get Beacon signal - No Beacon signal anymore
2822 2825 # newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
2823 2826 #
2824 2827 # if hei_ref != None:
2825 2828 # newheis = numpy.where(self.dataOut.heightList>hei_ref)
2826 2829 #
2827 2830
2828 2831
2829 2832 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
2830 2833 # see if the user put in pre defined phase shifts
2831 2834 voltsPShift = dataOut.data_pre.copy()
2832 2835
2833 2836 # if predefinedPhaseShifts != None:
2834 2837 # hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
2835 2838 #
2836 2839 # # elif beaconPhaseShifts:
2837 2840 # # #get hardware phase shifts using beacon signal
2838 2841 # # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
2839 2842 # # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
2840 2843 #
2841 2844 # else:
2842 2845 # hardwarePhaseShifts = numpy.zeros(5)
2843 2846 #
2844 2847 # voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
2845 2848 # for i in range(self.dataOut.data_pre.shape[0]):
2846 2849 # voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
2847 2850
2848 2851 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
2849 2852
2850 2853 #Remove DC
2851 2854 voltsDC = numpy.mean(voltsPShift,1)
2852 2855 voltsDC = numpy.mean(voltsDC,1)
2853 2856 for i in range(voltsDC.shape[0]):
2854 2857 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
2855 2858
2856 2859 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
2857 2860 # voltsPShift = voltsPShift[:,:,:newheis[0][0]]
2858 2861
2859 2862 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
2860 2863 #Coherent Detection
2861 2864 if cohDetection:
2862 2865 #use coherent detection to get the net power
2863 2866 cohDet_thresh = cohDet_thresh*numpy.pi/180
2864 2867 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, dataOut.timeInterval, pairslist0, cohDet_thresh)
2865 2868
2866 2869 #Non-coherent detection!
2867 2870 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
2868 2871 #********** END OF COH/NON-COH POWER CALCULATION**********************
2869 2872
2870 2873 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
2871 2874 #Get noise
2872 2875 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, dataOut.timeInterval)
2873 2876 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
2874 2877 #Get signal threshold
2875 2878 signalThresh = noise_multiple*noise
2876 2879 #Meteor echoes detection
2877 2880 listMeteors = self.__findMeteors(powerNet, signalThresh)
2878 2881 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
2879 2882
2880 2883 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
2881 2884 #Parameters
2882 2885 heiRange = dataOut.getHeiRange()
2883 2886 rangeInterval = heiRange[1] - heiRange[0]
2884 2887 rangeLimit = multDet_rangeLimit/rangeInterval
2885 2888 timeLimit = multDet_timeLimit/dataOut.timeInterval
2886 2889 #Multiple detection removals
2887 2890 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
2888 2891 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
2889 2892
2890 2893 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
2891 2894 #Parameters
2892 2895 phaseThresh = phaseThresh*numpy.pi/180
2893 2896 thresh = [phaseThresh, noise_multiple, SNRThresh]
2894 2897 #Meteor reestimation (Errors N 1, 6, 12, 17)
2895 2898 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, dataOut.timeInterval, dataOut.frequency)
2896 2899 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
2897 2900 #Estimation of decay times (Errors N 7, 8, 11)
2898 2901 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, dataOut.timeInterval, dataOut.frequency)
2899 2902 #******************* END OF METEOR REESTIMATION *******************
2900 2903
2901 2904 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
2902 2905 #Calculating Radial Velocity (Error N 15)
2903 2906 radialStdThresh = 10
2904 2907 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, dataOut.timeInterval)
2905 2908
2906 2909 if len(listMeteors4) > 0:
2907 2910 #Setting New Array
2908 2911 date = dataOut.utctime
2909 2912 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
2910 2913
2911 2914 #Correcting phase offset
2912 2915 if phaseOffsets != None:
2913 2916 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
2914 2917 arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
2915 2918
2916 2919 #Second Pairslist
2917 2920 pairsList = []
2918 2921 pairx = (0,1)
2919 2922 pairy = (2,3)
2920 2923 pairsList.append(pairx)
2921 2924 pairsList.append(pairy)
2922 2925
2923 2926 jph = numpy.array([0,0,0,0])
2924 2927 h = (hmin,hmax)
2925 2928 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
2926 2929
2927 2930 # #Calculate AOA (Error N 3, 4)
2928 2931 # #JONES ET AL. 1998
2929 2932 # error = arrayParameters[:,-1]
2930 2933 # AOAthresh = numpy.pi/8
2931 2934 # phases = -arrayParameters[:,9:13]
2932 2935 # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
2933 2936 #
2934 2937 # #Calculate Heights (Error N 13 and 14)
2935 2938 # error = arrayParameters[:,-1]
2936 2939 # Ranges = arrayParameters[:,2]
2937 2940 # zenith = arrayParameters[:,5]
2938 2941 # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
2939 2942 # error = arrayParameters[:,-1]
2940 2943 #********************* END OF PARAMETERS CALCULATION **************************
2941 2944
2942 2945 #***************************+ PASS DATA TO NEXT STEP **********************
2943 2946 # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
2944 2947 dataOut.data_param = arrayParameters
2945 2948
2946 2949 if arrayParameters == None:
2947 2950 dataOut.flagNoData = True
2948 2951 else:
2949 2952 dataOut.flagNoData = True
2950 2953
2951 2954 return
2952 2955
2953 2956 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
2954 2957
2955 2958 minIndex = min(newheis[0])
2956 2959 maxIndex = max(newheis[0])
2957 2960
2958 2961 voltage = voltage0[:,:,minIndex:maxIndex+1]
2959 2962 nLength = voltage.shape[1]/n
2960 2963 nMin = 0
2961 2964 nMax = 0
2962 2965 phaseOffset = numpy.zeros((len(pairslist),n))
2963 2966
2964 2967 for i in range(n):
2965 2968 nMax += nLength
2966 2969 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
2967 2970 phaseCCF = numpy.mean(phaseCCF, axis = 2)
2968 2971 phaseOffset[:,i] = phaseCCF.transpose()
2969 2972 nMin = nMax
2970 2973 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
2971 2974
2972 2975 #Remove Outliers
2973 2976 factor = 2
2974 2977 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
2975 2978 dw = numpy.std(wt,axis = 1)
2976 2979 dw = dw.reshape((dw.size,1))
2977 2980 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
2978 2981 phaseOffset[ind] = numpy.nan
2979 2982 phaseOffset = stats.nanmean(phaseOffset, axis=1)
2980 2983
2981 2984 return phaseOffset
2982 2985
2983 2986 def __shiftPhase(self, data, phaseShift):
2984 2987 #this will shift the phase of a complex number
2985 2988 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
2986 2989 return dataShifted
2987 2990
2988 2991 def __estimatePhaseDifference(self, array, pairslist):
2989 2992 nChannel = array.shape[0]
2990 2993 nHeights = array.shape[2]
2991 2994 numPairs = len(pairslist)
2992 2995 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
2993 2996 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
2994 2997
2995 2998 #Correct phases
2996 2999 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
2997 3000 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
2998 3001
2999 3002 if indDer[0].shape[0] > 0:
3000 3003 for i in range(indDer[0].shape[0]):
3001 3004 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
3002 3005 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
3003 3006
3004 3007 # for j in range(numSides):
3005 3008 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
3006 3009 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
3007 3010 #
3008 3011 #Linear
3009 3012 phaseInt = numpy.zeros((numPairs,1))
3010 3013 angAllCCF = phaseCCF[:,[0,1,3,4],0]
3011 3014 for j in range(numPairs):
3012 3015 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
3013 3016 phaseInt[j] = fit[1]
3014 3017 #Phase Differences
3015 3018 phaseDiff = phaseInt - phaseCCF[:,2,:]
3016 3019 phaseArrival = phaseInt.reshape(phaseInt.size)
3017 3020
3018 3021 #Dealias
3019 3022 phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))
3020 3023 # indAlias = numpy.where(phaseArrival > numpy.pi)
3021 3024 # phaseArrival[indAlias] -= 2*numpy.pi
3022 3025 # indAlias = numpy.where(phaseArrival < -numpy.pi)
3023 3026 # phaseArrival[indAlias] += 2*numpy.pi
3024 3027
3025 3028 return phaseDiff, phaseArrival
3026 3029
3027 3030 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
3028 3031 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
3029 3032 #find the phase shifts of each channel over 1 second intervals
3030 3033 #only look at ranges below the beacon signal
3031 3034 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
3032 3035 numBlocks = int(volts.shape[1]/numProfPerBlock)
3033 3036 numHeights = volts.shape[2]
3034 3037 nChannel = volts.shape[0]
3035 3038 voltsCohDet = volts.copy()
3036 3039
3037 3040 pairsarray = numpy.array(pairslist)
3038 3041 indSides = pairsarray[:,1]
3039 3042 # indSides = numpy.array(range(nChannel))
3040 3043 # indSides = numpy.delete(indSides, indCenter)
3041 3044 #
3042 3045 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
3043 3046 listBlocks = numpy.array_split(volts, numBlocks, 1)
3044 3047
3045 3048 startInd = 0
3046 3049 endInd = 0
3047 3050
3048 3051 for i in range(numBlocks):
3049 3052 startInd = endInd
3050 3053 endInd = endInd + listBlocks[i].shape[1]
3051 3054
3052 3055 arrayBlock = listBlocks[i]
3053 3056 # arrayBlockCenter = listCenter[i]
3054 3057
3055 3058 #Estimate the Phase Difference
3056 3059 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
3057 3060 #Phase Difference RMS
3058 3061 arrayPhaseRMS = numpy.abs(phaseDiff)
3059 3062 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
3060 3063 indPhase = numpy.where(phaseRMSaux==4)
3061 3064 #Shifting
3062 3065 if indPhase[0].shape[0] > 0:
3063 3066 for j in range(indSides.size):
3064 3067 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
3065 3068 voltsCohDet[:,startInd:endInd,:] = arrayBlock
3066 3069
3067 3070 return voltsCohDet
3068 3071
3069 3072 def __calculateCCF(self, volts, pairslist ,laglist):
3070 3073
3071 3074 nHeights = volts.shape[2]
3072 3075 nPoints = volts.shape[1]
3073 3076 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
3074 3077
3075 3078 for i in range(len(pairslist)):
3076 3079 volts1 = volts[pairslist[i][0]]
3077 3080 volts2 = volts[pairslist[i][1]]
3078 3081
3079 3082 for t in range(len(laglist)):
3080 3083 idxT = laglist[t]
3081 3084 if idxT >= 0:
3082 3085 vStacked = numpy.vstack((volts2[idxT:,:],
3083 3086 numpy.zeros((idxT, nHeights),dtype='complex')))
3084 3087 else:
3085 3088 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
3086 3089 volts2[:(nPoints + idxT),:]))
3087 3090 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
3088 3091
3089 3092 vStacked = None
3090 3093 return voltsCCF
3091 3094
3092 3095 def __getNoise(self, power, timeSegment, timeInterval):
3093 3096 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
3094 3097 numBlocks = int(power.shape[0]/numProfPerBlock)
3095 3098 numHeights = power.shape[1]
3096 3099
3097 3100 listPower = numpy.array_split(power, numBlocks, 0)
3098 3101 noise = numpy.zeros((power.shape[0], power.shape[1]))
3099 3102 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
3100 3103
3101 3104 startInd = 0
3102 3105 endInd = 0
3103 3106
3104 3107 for i in range(numBlocks): #split por canal
3105 3108 startInd = endInd
3106 3109 endInd = endInd + listPower[i].shape[0]
3107 3110
3108 3111 arrayBlock = listPower[i]
3109 3112 noiseAux = numpy.mean(arrayBlock, 0)
3110 3113 # noiseAux = numpy.median(noiseAux)
3111 3114 # noiseAux = numpy.mean(arrayBlock)
3112 3115 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
3113 3116
3114 3117 noiseAux1 = numpy.mean(arrayBlock)
3115 3118 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
3116 3119
3117 3120 return noise, noise1
3118 3121
3119 3122 def __findMeteors(self, power, thresh):
3120 3123 nProf = power.shape[0]
3121 3124 nHeights = power.shape[1]
3122 3125 listMeteors = []
3123 3126
3124 3127 for i in range(nHeights):
3125 3128 powerAux = power[:,i]
3126 3129 threshAux = thresh[:,i]
3127 3130
3128 3131 indUPthresh = numpy.where(powerAux > threshAux)[0]
3129 3132 indDNthresh = numpy.where(powerAux <= threshAux)[0]
3130 3133
3131 3134 j = 0
3132 3135
3133 3136 while (j < indUPthresh.size - 2):
3134 3137 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
3135 3138 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
3136 3139 indDNthresh = indDNthresh[indDNAux]
3137 3140
3138 3141 if (indDNthresh.size > 0):
3139 3142 indEnd = indDNthresh[0] - 1
3140 3143 indInit = indUPthresh[j]
3141 3144
3142 3145 meteor = powerAux[indInit:indEnd + 1]
3143 3146 indPeak = meteor.argmax() + indInit
3144 3147 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
3145 3148
3146 3149 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
3147 3150 j = numpy.where(indUPthresh == indEnd)[0] + 1
3148 3151 else: j+=1
3149 3152 else: j+=1
3150 3153
3151 3154 return listMeteors
3152 3155
3153 3156 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
3154 3157
3155 3158 arrayMeteors = numpy.asarray(listMeteors)
3156 3159 listMeteors1 = []
3157 3160
3158 3161 while arrayMeteors.shape[0] > 0:
3159 3162 FLAs = arrayMeteors[:,4]
3160 3163 maxFLA = FLAs.argmax()
3161 3164 listMeteors1.append(arrayMeteors[maxFLA,:])
3162 3165
3163 3166 MeteorInitTime = arrayMeteors[maxFLA,1]
3164 3167 MeteorEndTime = arrayMeteors[maxFLA,3]
3165 3168 MeteorHeight = arrayMeteors[maxFLA,0]
3166 3169
3167 3170 #Check neighborhood
3168 3171 maxHeightIndex = MeteorHeight + rangeLimit
3169 3172 minHeightIndex = MeteorHeight - rangeLimit
3170 3173 minTimeIndex = MeteorInitTime - timeLimit
3171 3174 maxTimeIndex = MeteorEndTime + timeLimit
3172 3175
3173 3176 #Check Heights
3174 3177 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
3175 3178 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
3176 3179 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
3177 3180
3178 3181 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
3179 3182
3180 3183 return listMeteors1
3181 3184
3182 3185 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
3183 3186 numHeights = volts.shape[2]
3184 3187 nChannel = volts.shape[0]
3185 3188
3186 3189 thresholdPhase = thresh[0]
3187 3190 thresholdNoise = thresh[1]
3188 3191 thresholdDB = float(thresh[2])
3189 3192
3190 3193 thresholdDB1 = 10**(thresholdDB/10)
3191 3194 pairsarray = numpy.array(pairslist)
3192 3195 indSides = pairsarray[:,1]
3193 3196
3194 3197 pairslist1 = list(pairslist)
3195 3198 pairslist1.append((0,1))
3196 3199 pairslist1.append((3,4))
3197 3200
3198 3201 listMeteors1 = []
3199 3202 listPowerSeries = []
3200 3203 listVoltageSeries = []
3201 3204 #volts has the war data
3202 3205
3203 3206 if frequency == 30e6:
3204 3207 timeLag = 45*10**-3
3205 3208 else:
3206 3209 timeLag = 15*10**-3
3207 3210 lag = numpy.ceil(timeLag/timeInterval)
3208 3211
3209 3212 for i in range(len(listMeteors)):
3210 3213
3211 3214 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
3212 3215 meteorAux = numpy.zeros(16)
3213 3216
3214 3217 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
3215 3218 mHeight = listMeteors[i][0]
3216 3219 mStart = listMeteors[i][1]
3217 3220 mPeak = listMeteors[i][2]
3218 3221 mEnd = listMeteors[i][3]
3219 3222
3220 3223 #get the volt data between the start and end times of the meteor
3221 3224 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
3222 3225 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
3223 3226
3224 3227 #3.6. Phase Difference estimation
3225 3228 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
3226 3229
3227 3230 #3.7. Phase difference removal & meteor start, peak and end times reestimated
3228 3231 #meteorVolts0.- all Channels, all Profiles
3229 3232 meteorVolts0 = volts[:,:,mHeight]
3230 3233 meteorThresh = noise[:,mHeight]*thresholdNoise
3231 3234 meteorNoise = noise[:,mHeight]
3232 3235 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
3233 3236 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
3234 3237
3235 3238 #Times reestimation
3236 3239 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
3237 3240 if mStart1.size > 0:
3238 3241 mStart1 = mStart1[-1] + 1
3239 3242
3240 3243 else:
3241 3244 mStart1 = mPeak
3242 3245
3243 3246 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
3244 3247 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
3245 3248 if mEndDecayTime1.size == 0:
3246 3249 mEndDecayTime1 = powerNet0.size
3247 3250 else:
3248 3251 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
3249 3252 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
3250 3253
3251 3254 #meteorVolts1.- all Channels, from start to end
3252 3255 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
3253 3256 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
3254 3257 if meteorVolts2.shape[1] == 0:
3255 3258 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
3256 3259 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
3257 3260 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
3258 3261 ##################### END PARAMETERS REESTIMATION #########################
3259 3262
3260 3263 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
3261 3264 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
3262 3265 if meteorVolts2.shape[1] > 0:
3263 3266 #Phase Difference re-estimation
3264 3267 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
3265 3268 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
3266 3269 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
3267 3270 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
3268 3271 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
3269 3272
3270 3273 #Phase Difference RMS
3271 3274 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
3272 3275 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
3273 3276 #Data from Meteor
3274 3277 mPeak1 = powerNet1.argmax() + mStart1
3275 3278 mPeakPower1 = powerNet1.max()
3276 3279 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
3277 3280 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
3278 3281 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
3279 3282 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
3280 3283 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
3281 3284 #Vectorize
3282 3285 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
3283 3286 meteorAux[7:11] = phaseDiffint[0:4]
3284 3287
3285 3288 #Rejection Criterions
3286 3289 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
3287 3290 meteorAux[-1] = 17
3288 3291 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
3289 3292 meteorAux[-1] = 1
3290 3293
3291 3294
3292 3295 else:
3293 3296 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
3294 3297 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
3295 3298 PowerSeries = 0
3296 3299
3297 3300 listMeteors1.append(meteorAux)
3298 3301 listPowerSeries.append(PowerSeries)
3299 3302 listVoltageSeries.append(meteorVolts1)
3300 3303
3301 3304 return listMeteors1, listPowerSeries, listVoltageSeries
3302 3305
3303 3306 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
3304 3307
3305 3308 threshError = 10
3306 3309 #Depending if it is 30 or 50 MHz
3307 3310 if frequency == 30e6:
3308 3311 timeLag = 45*10**-3
3309 3312 else:
3310 3313 timeLag = 15*10**-3
3311 3314 lag = numpy.ceil(timeLag/timeInterval)
3312 3315
3313 3316 listMeteors1 = []
3314 3317
3315 3318 for i in range(len(listMeteors)):
3316 3319 meteorPower = listPower[i]
3317 3320 meteorAux = listMeteors[i]
3318 3321
3319 3322 if meteorAux[-1] == 0:
3320 3323
3321 3324 try:
3322 3325 indmax = meteorPower.argmax()
3323 3326 indlag = indmax + lag
3324 3327
3325 3328 y = meteorPower[indlag:]
3326 3329 x = numpy.arange(0, y.size)*timeLag
3327 3330
3328 3331 #first guess
3329 3332 a = y[0]
3330 3333 tau = timeLag
3331 3334 #exponential fit
3332 3335 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
3333 3336 y1 = self.__exponential_function(x, *popt)
3334 3337 #error estimation
3335 3338 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
3336 3339
3337 3340 decayTime = popt[1]
3338 3341 riseTime = indmax*timeInterval
3339 3342 meteorAux[11:13] = [decayTime, error]
3340 3343
3341 3344 #Table items 7, 8 and 11
3342 3345 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
3343 3346 meteorAux[-1] = 7
3344 3347 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
3345 3348 meteorAux[-1] = 8
3346 3349 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
3347 3350 meteorAux[-1] = 11
3348 3351
3349 3352
3350 3353 except:
3351 3354 meteorAux[-1] = 11
3352 3355
3353 3356
3354 3357 listMeteors1.append(meteorAux)
3355 3358
3356 3359 return listMeteors1
3357 3360
3358 3361 #Exponential Function
3359 3362
3360 3363 def __exponential_function(self, x, a, tau):
3361 3364 y = a*numpy.exp(-x/tau)
3362 3365 return y
3363 3366
3364 3367 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
3365 3368
3366 3369 pairslist1 = list(pairslist)
3367 3370 pairslist1.append((0,1))
3368 3371 pairslist1.append((3,4))
3369 3372 numPairs = len(pairslist1)
3370 3373 #Time Lag
3371 3374 timeLag = 45*10**-3
3372 3375 c = 3e8
3373 3376 lag = numpy.ceil(timeLag/timeInterval)
3374 3377 freq = 30e6
3375 3378
3376 3379 listMeteors1 = []
3377 3380
3378 3381 for i in range(len(listMeteors)):
3379 3382 meteorAux = listMeteors[i]
3380 3383 if meteorAux[-1] == 0:
3381 3384 mStart = listMeteors[i][1]
3382 3385 mPeak = listMeteors[i][2]
3383 3386 mLag = mPeak - mStart + lag
3384 3387
3385 3388 #get the volt data between the start and end times of the meteor
3386 3389 meteorVolts = listVolts[i]
3387 3390 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
3388 3391
3389 3392 #Get CCF
3390 3393 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
3391 3394
3392 3395 #Method 2
3393 3396 slopes = numpy.zeros(numPairs)
3394 3397 time = numpy.array([-2,-1,1,2])*timeInterval
3395 3398 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
3396 3399
3397 3400 #Correct phases
3398 3401 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
3399 3402 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
3400 3403
3401 3404 if indDer[0].shape[0] > 0:
3402 3405 for i in range(indDer[0].shape[0]):
3403 3406 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
3404 3407 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
3405 3408
3406 3409 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
3407 3410 for j in range(numPairs):
3408 3411 fit = stats.linregress(time, angAllCCF[j,:])
3409 3412 slopes[j] = fit[0]
3410 3413
3411 3414 #Remove Outlier
3412 3415 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
3413 3416 # slopes = numpy.delete(slopes,indOut)
3414 3417 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
3415 3418 # slopes = numpy.delete(slopes,indOut)
3416 3419
3417 3420 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
3418 3421 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
3419 3422 meteorAux[-2] = radialError
3420 3423 meteorAux[-3] = radialVelocity
3421 3424
3422 3425 #Setting Error
3423 3426 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
3424 3427 if numpy.abs(radialVelocity) > 200:
3425 3428 meteorAux[-1] = 15
3426 3429 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
3427 3430 elif radialError > radialStdThresh:
3428 3431 meteorAux[-1] = 12
3429 3432
3430 3433 listMeteors1.append(meteorAux)
3431 3434 return listMeteors1
3432 3435
3433 3436 def __setNewArrays(self, listMeteors, date, heiRang):
3434 3437
3435 3438 #New arrays
3436 3439 arrayMeteors = numpy.array(listMeteors)
3437 3440 arrayParameters = numpy.zeros((len(listMeteors), 13))
3438 3441
3439 3442 #Date inclusion
3440 3443 # date = re.findall(r'\((.*?)\)', date)
3441 3444 # date = date[0].split(',')
3442 3445 # date = map(int, date)
3443 3446 #
3444 3447 # if len(date)<6:
3445 3448 # date.append(0)
3446 3449 #
3447 3450 # date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
3448 3451 # arrayDate = numpy.tile(date, (len(listMeteors), 1))
3449 3452 arrayDate = numpy.tile(date, (len(listMeteors)))
3450 3453
3451 3454 #Meteor array
3452 3455 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
3453 3456 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
3454 3457
3455 3458 #Parameters Array
3456 3459 arrayParameters[:,0] = arrayDate #Date
3457 3460 arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range
3458 3461 arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error
3459 3462 arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases
3460 3463 arrayParameters[:,-1] = arrayMeteors[:,-1] #Error
3461 3464
3462 3465
3463 3466 return arrayParameters
3464 3467
3465 3468 class CorrectSMPhases(Operation):
3466 3469
3467 3470 def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):
3468 3471
3469 3472 arrayParameters = dataOut.data_param
3470 3473 pairsList = []
3471 3474 pairx = (0,1)
3472 3475 pairy = (2,3)
3473 3476 pairsList.append(pairx)
3474 3477 pairsList.append(pairy)
3475 3478 jph = numpy.zeros(4)
3476 3479
3477 3480 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
3478 3481 # arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
3479 3482 arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
3480 3483
3481 3484 meteorOps = SMOperations()
3482 3485 if channelPositions == None:
3483 3486 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
3484 3487 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
3485 3488
3486 3489 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
3487 3490 h = (hmin,hmax)
3488 3491
3489 3492 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
3490 3493
3491 3494 dataOut.data_param = arrayParameters
3492 3495 return
3493 3496
3494 3497 class SMPhaseCalibration(Operation):
3495 3498
3496 3499 __buffer = None
3497 3500
3498 3501 __initime = None
3499 3502
3500 3503 __dataReady = False
3501 3504
3502 3505 __isConfig = False
3503 3506
3504 3507 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
3505 3508
3506 3509 dataTime = currentTime + paramInterval
3507 3510 deltaTime = dataTime - initTime
3508 3511
3509 3512 if deltaTime >= outputInterval or deltaTime < 0:
3510 3513 return True
3511 3514
3512 3515 return False
3513 3516
3514 3517 def __getGammas(self, pairs, d, phases):
3515 3518 gammas = numpy.zeros(2)
3516 3519
3517 3520 for i in range(len(pairs)):
3518 3521
3519 3522 pairi = pairs[i]
3520 3523
3521 3524 phip3 = phases[:,pairi[1]]
3522 3525 d3 = d[pairi[1]]
3523 3526 phip2 = phases[:,pairi[0]]
3524 3527 d2 = d[pairi[0]]
3525 3528 #Calculating gamma
3526 3529 # jdcos = alp1/(k*d1)
3527 3530 # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
3528 3531 jgamma = -phip2*d3/d2 - phip3
3529 3532 jgamma = numpy.angle(numpy.exp(1j*jgamma))
3530 3533 # jgamma[jgamma>numpy.pi] -= 2*numpy.pi
3531 3534 # jgamma[jgamma<-numpy.pi] += 2*numpy.pi
3532 3535
3533 3536 #Revised distribution
3534 3537 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
3535 3538
3536 3539 #Histogram
3537 3540 nBins = 64.0
3538 3541 rmin = -0.5*numpy.pi
3539 3542 rmax = 0.5*numpy.pi
3540 3543 phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
3541 3544
3542 3545 meteorsY = phaseHisto[0]
3543 3546 phasesX = phaseHisto[1][:-1]
3544 3547 width = phasesX[1] - phasesX[0]
3545 3548 phasesX += width/2
3546 3549
3547 3550 #Gaussian aproximation
3548 3551 bpeak = meteorsY.argmax()
3549 3552 peak = meteorsY.max()
3550 3553 jmin = bpeak - 5
3551 3554 jmax = bpeak + 5 + 1
3552 3555
3553 3556 if jmin<0:
3554 3557 jmin = 0
3555 3558 jmax = 6
3556 3559 elif jmax > meteorsY.size:
3557 3560 jmin = meteorsY.size - 6
3558 3561 jmax = meteorsY.size
3559 3562
3560 3563 x0 = numpy.array([peak,bpeak,50])
3561 3564 coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
3562 3565
3563 3566 #Gammas
3564 3567 gammas[i] = coeff[0][1]
3565 3568
3566 3569 return gammas
3567 3570
3568 3571 def __residualFunction(self, coeffs, y, t):
3569 3572
3570 3573 return y - self.__gauss_function(t, coeffs)
3571 3574
3572 3575 def __gauss_function(self, t, coeffs):
3573 3576
3574 3577 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
3575 3578
3576 3579 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
3577 3580 meteorOps = SMOperations()
3578 3581 nchan = 4
3579 3582 pairx = pairsList[0]
3580 3583 pairy = pairsList[1]
3581 3584 center_xangle = 0
3582 3585 center_yangle = 0
3583 3586 range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
3584 3587 ntimes = len(range_angle)
3585 3588
3586 3589 nstepsx = 20.0
3587 3590 nstepsy = 20.0
3588 3591
3589 3592 for iz in range(ntimes):
3590 3593 min_xangle = -range_angle[iz]/2 + center_xangle
3591 3594 max_xangle = range_angle[iz]/2 + center_xangle
3592 3595 min_yangle = -range_angle[iz]/2 + center_yangle
3593 3596 max_yangle = range_angle[iz]/2 + center_yangle
3594 3597
3595 3598 inc_x = (max_xangle-min_xangle)/nstepsx
3596 3599 inc_y = (max_yangle-min_yangle)/nstepsy
3597 3600
3598 3601 alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
3599 3602 alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
3600 3603 penalty = numpy.zeros((nstepsx,nstepsy))
3601 3604 jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
3602 3605 jph = numpy.zeros(nchan)
3603 3606
3604 3607 # Iterations looking for the offset
3605 3608 for iy in range(int(nstepsy)):
3606 3609 for ix in range(int(nstepsx)):
3607 3610 jph[pairy[1]] = alpha_y[iy]
3608 3611 jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
3609 3612
3610 3613 jph[pairx[1]] = alpha_x[ix]
3611 3614 jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
3612 3615
3613 3616 jph_array[:,ix,iy] = jph
3614 3617
3615 3618 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
3616 3619 error = meteorsArray1[:,-1]
3617 3620 ind1 = numpy.where(error==0)[0]
3618 3621 penalty[ix,iy] = ind1.size
3619 3622
3620 3623 i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
3621 3624 phOffset = jph_array[:,i,j]
3622 3625
3623 3626 center_xangle = phOffset[pairx[1]]
3624 3627 center_yangle = phOffset[pairy[1]]
3625 3628
3626 3629 phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
3627 3630 phOffset = phOffset*180/numpy.pi
3628 3631 return phOffset
3629 3632
3630 3633
3631 3634 def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1):
3632 3635
3633 3636 dataOut.flagNoData = True
3634 3637 self.__dataReady = False
3635 3638 dataOut.outputInterval = nHours*3600
3636 3639
3637 3640 if self.__isConfig == False:
3638 3641 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
3639 3642 #Get Initial LTC time
3640 3643 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
3641 3644 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
3642 3645
3643 3646 self.__isConfig = True
3644 3647
3645 3648 if self.__buffer == None:
3646 3649 self.__buffer = dataOut.data_param.copy()
3647 3650
3648 3651 else:
3649 3652 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
3650 3653
3651 3654 self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
3652 3655
3653 3656 if self.__dataReady:
3654 3657 dataOut.utctimeInit = self.__initime
3655 3658 self.__initime += dataOut.outputInterval #to erase time offset
3656 3659
3657 3660 freq = dataOut.frequency
3658 3661 c = dataOut.C #m/s
3659 3662 lamb = c/freq
3660 3663 k = 2*numpy.pi/lamb
3661 3664 azimuth = 0
3662 3665 h = (hmin, hmax)
3663 3666 pairs = ((0,1),(2,3))
3664 3667
3665 3668 if channelPositions == None:
3666 3669 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
3667 3670 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
3668 3671 meteorOps = SMOperations()
3669 3672 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
3670 3673
3671 3674 # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
3672 3675
3673 3676 meteorsArray = self.__buffer
3674 3677 error = meteorsArray[:,-1]
3675 3678 boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
3676 3679 ind1 = numpy.where(boolError)[0]
3677 3680 meteorsArray = meteorsArray[ind1,:]
3678 3681 meteorsArray[:,-1] = 0
3679 3682 phases = meteorsArray[:,8:12]
3680 3683
3681 3684 #Calculate Gammas
3682 3685 gammas = self.__getGammas(pairs, distances, phases)
3683 3686 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
3684 3687 #Calculate Phases
3685 3688 phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
3686 3689 phasesOff = phasesOff.reshape((1,phasesOff.size))
3687 3690 dataOut.data_output = -phasesOff
3688 3691 dataOut.flagNoData = False
3689 3692 self.__buffer = None
3690 3693
3691 3694
3692 3695 return
3693 3696
3694 3697 class SMOperations():
3695 3698
3696 3699 def __init__(self):
3697 3700
3698 3701 return
3699 3702
3700 3703 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph):
3701 3704
3702 3705 arrayParameters = arrayParameters0.copy()
3703 3706 hmin = h[0]
3704 3707 hmax = h[1]
3705 3708
3706 3709 #Calculate AOA (Error N 3, 4)
3707 3710 #JONES ET AL. 1998
3708 3711 AOAthresh = numpy.pi/8
3709 3712 error = arrayParameters[:,-1]
3710 3713 phases = -arrayParameters[:,8:12] + jph
3711 3714 # phases = numpy.unwrap(phases)
3712 3715 arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth)
3713 3716
3714 3717 #Calculate Heights (Error N 13 and 14)
3715 3718 error = arrayParameters[:,-1]
3716 3719 Ranges = arrayParameters[:,1]
3717 3720 zenith = arrayParameters[:,4]
3718 3721 arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
3719 3722
3720 3723 #----------------------- Get Final data ------------------------------------
3721 3724 # error = arrayParameters[:,-1]
3722 3725 # ind1 = numpy.where(error==0)[0]
3723 3726 # arrayParameters = arrayParameters[ind1,:]
3724 3727
3725 3728 return arrayParameters
3726 3729
3727 3730 def __getAOA(self, phases, pairsList, directions, error, AOAthresh, azimuth):
3728 3731
3729 3732 arrayAOA = numpy.zeros((phases.shape[0],3))
3730 3733 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions)
3731 3734
3732 3735 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
3733 3736 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
3734 3737 arrayAOA[:,2] = cosDirError
3735 3738
3736 3739 azimuthAngle = arrayAOA[:,0]
3737 3740 zenithAngle = arrayAOA[:,1]
3738 3741
3739 3742 #Setting Error
3740 3743 indError = numpy.where(numpy.logical_or(error == 3, error == 4))[0]
3741 3744 error[indError] = 0
3742 3745 #Number 3: AOA not fesible
3743 3746 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
3744 3747 error[indInvalid] = 3
3745 3748 #Number 4: Large difference in AOAs obtained from different antenna baselines
3746 3749 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
3747 3750 error[indInvalid] = 4
3748 3751 return arrayAOA, error
3749 3752
3750 3753 def __getDirectionCosines(self, arrayPhase, pairsList, distances):
3751 3754
3752 3755 #Initializing some variables
3753 3756 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
3754 3757 ang_aux = ang_aux.reshape(1,ang_aux.size)
3755 3758
3756 3759 cosdir = numpy.zeros((arrayPhase.shape[0],2))
3757 3760 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
3758 3761
3759 3762
3760 3763 for i in range(2):
3761 3764 ph0 = arrayPhase[:,pairsList[i][0]]
3762 3765 ph1 = arrayPhase[:,pairsList[i][1]]
3763 3766 d0 = distances[pairsList[i][0]]
3764 3767 d1 = distances[pairsList[i][1]]
3765 3768
3766 3769 ph0_aux = ph0 + ph1
3767 3770 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
3768 3771 # ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi
3769 3772 # ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi
3770 3773 #First Estimation
3771 3774 cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
3772 3775
3773 3776 #Most-Accurate Second Estimation
3774 3777 phi1_aux = ph0 - ph1
3775 3778 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
3776 3779 #Direction Cosine 1
3777 3780 cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1))
3778 3781
3779 3782 #Searching the correct Direction Cosine
3780 3783 cosdir0_aux = cosdir0[:,i]
3781 3784 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
3782 3785 #Minimum Distance
3783 3786 cosDiff = (cosdir1 - cosdir0_aux)**2
3784 3787 indcos = cosDiff.argmin(axis = 1)
3785 3788 #Saving Value obtained
3786 3789 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
3787 3790
3788 3791 return cosdir0, cosdir
3789 3792
3790 3793 def __calculateAOA(self, cosdir, azimuth):
3791 3794 cosdirX = cosdir[:,0]
3792 3795 cosdirY = cosdir[:,1]
3793 3796
3794 3797 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
3795 3798 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth#0 deg north, 90 deg east
3796 3799 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
3797 3800
3798 3801 return angles
3799 3802
3800 3803 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
3801 3804
3802 3805 Ramb = 375 #Ramb = c/(2*PRF)
3803 3806 Re = 6371 #Earth Radius
3804 3807 heights = numpy.zeros(Ranges.shape)
3805 3808
3806 3809 R_aux = numpy.array([0,1,2])*Ramb
3807 3810 R_aux = R_aux.reshape(1,R_aux.size)
3808 3811
3809 3812 Ranges = Ranges.reshape(Ranges.size,1)
3810 3813
3811 3814 Ri = Ranges + R_aux
3812 3815 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
3813 3816
3814 3817 #Check if there is a height between 70 and 110 km
3815 3818 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
3816 3819 ind_h = numpy.where(h_bool == 1)[0]
3817 3820
3818 3821 hCorr = hi[ind_h, :]
3819 3822 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
3820 3823
3821 3824 hCorr = hi[ind_hCorr]
3822 3825 heights[ind_h] = hCorr
3823 3826
3824 3827 #Setting Error
3825 3828 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
3826 3829 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
3827 3830 indError = numpy.where(numpy.logical_or(error == 13, error == 14))[0]
3828 3831 error[indError] = 0
3829 3832 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
3830 3833 error[indInvalid2] = 14
3831 3834 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
3832 3835 error[indInvalid1] = 13
3833 3836
3834 3837 return heights, error
3835 3838
3836 3839 def getPhasePairs(self, channelPositions):
3837 3840 chanPos = numpy.array(channelPositions)
3838 3841 listOper = list(itertools.combinations(range(5),2))
3839 3842
3840 3843 distances = numpy.zeros(4)
3841 3844 axisX = []
3842 3845 axisY = []
3843 3846 distX = numpy.zeros(3)
3844 3847 distY = numpy.zeros(3)
3845 3848 ix = 0
3846 3849 iy = 0
3847 3850
3848 3851 pairX = numpy.zeros((2,2))
3849 3852 pairY = numpy.zeros((2,2))
3850 3853
3851 3854 for i in range(len(listOper)):
3852 3855 pairi = listOper[i]
3853 3856
3854 3857 posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:])
3855 3858
3856 3859 if posDif[0] == 0:
3857 3860 axisY.append(pairi)
3858 3861 distY[iy] = posDif[1]
3859 3862 iy += 1
3860 3863 elif posDif[1] == 0:
3861 3864 axisX.append(pairi)
3862 3865 distX[ix] = posDif[0]
3863 3866 ix += 1
3864 3867
3865 3868 for i in range(2):
3866 3869 if i==0:
3867 3870 dist0 = distX
3868 3871 axis0 = axisX
3869 3872 else:
3870 3873 dist0 = distY
3871 3874 axis0 = axisY
3872 3875
3873 3876 side = numpy.argsort(dist0)[:-1]
3874 3877 axis0 = numpy.array(axis0)[side,:]
3875 3878 chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0])
3876 3879 axis1 = numpy.unique(numpy.reshape(axis0,4))
3877 3880 side = axis1[axis1 != chanC]
3878 3881 diff1 = chanPos[chanC,i] - chanPos[side[0],i]
3879 3882 diff2 = chanPos[chanC,i] - chanPos[side[1],i]
3880 3883 if diff1<0:
3881 3884 chan2 = side[0]
3882 3885 d2 = numpy.abs(diff1)
3883 3886 chan1 = side[1]
3884 3887 d1 = numpy.abs(diff2)
3885 3888 else:
3886 3889 chan2 = side[1]
3887 3890 d2 = numpy.abs(diff2)
3888 3891 chan1 = side[0]
3889 3892 d1 = numpy.abs(diff1)
3890 3893
3891 3894 if i==0:
3892 3895 chanCX = chanC
3893 3896 chan1X = chan1
3894 3897 chan2X = chan2
3895 3898 distances[0:2] = numpy.array([d1,d2])
3896 3899 else:
3897 3900 chanCY = chanC
3898 3901 chan1Y = chan1
3899 3902 chan2Y = chan2
3900 3903 distances[2:4] = numpy.array([d1,d2])
3901 3904 # axisXsides = numpy.reshape(axisX[ix,:],4)
3902 3905 #
3903 3906 # channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0])
3904 3907 # channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0])
3905 3908 #
3906 3909 # ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0]
3907 3910 # ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0]
3908 3911 # channel25X = int(pairX[0,ind25X])
3909 3912 # channel20X = int(pairX[1,ind20X])
3910 3913 # ind25Y = numpy.where(pairY[0,:] != channelCentY)[0][0]
3911 3914 # ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0]
3912 3915 # channel25Y = int(pairY[0,ind25Y])
3913 3916 # channel20Y = int(pairY[1,ind20Y])
3914 3917
3915 3918 # pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
3916 3919 pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
3917 3920
3918 3921 return pairslist, distances
3919 3922 # def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
3920 3923 #
3921 3924 # arrayAOA = numpy.zeros((phases.shape[0],3))
3922 3925 # cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
3923 3926 #
3924 3927 # arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
3925 3928 # cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
3926 3929 # arrayAOA[:,2] = cosDirError
3927 3930 #
3928 3931 # azimuthAngle = arrayAOA[:,0]
3929 3932 # zenithAngle = arrayAOA[:,1]
3930 3933 #
3931 3934 # #Setting Error
3932 3935 # #Number 3: AOA not fesible
3933 3936 # indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
3934 3937 # error[indInvalid] = 3
3935 3938 # #Number 4: Large difference in AOAs obtained from different antenna baselines
3936 3939 # indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
3937 3940 # error[indInvalid] = 4
3938 3941 # return arrayAOA, error
3939 3942 #
3940 3943 # def __getDirectionCosines(self, arrayPhase, pairsList):
3941 3944 #
3942 3945 # #Initializing some variables
3943 3946 # ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
3944 3947 # ang_aux = ang_aux.reshape(1,ang_aux.size)
3945 3948 #
3946 3949 # cosdir = numpy.zeros((arrayPhase.shape[0],2))
3947 3950 # cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
3948 3951 #
3949 3952 #
3950 3953 # for i in range(2):
3951 3954 # #First Estimation
3952 3955 # phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
3953 3956 # #Dealias
3954 3957 # indcsi = numpy.where(phi0_aux > numpy.pi)
3955 3958 # phi0_aux[indcsi] -= 2*numpy.pi
3956 3959 # indcsi = numpy.where(phi0_aux < -numpy.pi)
3957 3960 # phi0_aux[indcsi] += 2*numpy.pi
3958 3961 # #Direction Cosine 0
3959 3962 # cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
3960 3963 #
3961 3964 # #Most-Accurate Second Estimation
3962 3965 # phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
3963 3966 # phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
3964 3967 # #Direction Cosine 1
3965 3968 # cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
3966 3969 #
3967 3970 # #Searching the correct Direction Cosine
3968 3971 # cosdir0_aux = cosdir0[:,i]
3969 3972 # cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
3970 3973 # #Minimum Distance
3971 3974 # cosDiff = (cosdir1 - cosdir0_aux)**2
3972 3975 # indcos = cosDiff.argmin(axis = 1)
3973 3976 # #Saving Value obtained
3974 3977 # cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
3975 3978 #
3976 3979 # return cosdir0, cosdir
3977 3980 #
3978 3981 # def __calculateAOA(self, cosdir, azimuth):
3979 3982 # cosdirX = cosdir[:,0]
3980 3983 # cosdirY = cosdir[:,1]
3981 3984 #
3982 3985 # zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
3983 3986 # azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
3984 3987 # angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
3985 3988 #
3986 3989 # return angles
3987 3990 #
3988 3991 # def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
3989 3992 #
3990 3993 # Ramb = 375 #Ramb = c/(2*PRF)
3991 3994 # Re = 6371 #Earth Radius
3992 3995 # heights = numpy.zeros(Ranges.shape)
3993 3996 #
3994 3997 # R_aux = numpy.array([0,1,2])*Ramb
3995 3998 # R_aux = R_aux.reshape(1,R_aux.size)
3996 3999 #
3997 4000 # Ranges = Ranges.reshape(Ranges.size,1)
3998 4001 #
3999 4002 # Ri = Ranges + R_aux
4000 4003 # hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
4001 4004 #
4002 4005 # #Check if there is a height between 70 and 110 km
4003 4006 # h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
4004 4007 # ind_h = numpy.where(h_bool == 1)[0]
4005 4008 #
4006 4009 # hCorr = hi[ind_h, :]
4007 4010 # ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
4008 4011 #
4009 4012 # hCorr = hi[ind_hCorr]
4010 4013 # heights[ind_h] = hCorr
4011 4014 #
4012 4015 # #Setting Error
4013 4016 # #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
4014 4017 # #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
4015 4018 #
4016 4019 # indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
4017 4020 # error[indInvalid2] = 14
4018 4021 # indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
4019 4022 # error[indInvalid1] = 13
4020 4023 #
4021 4024 # return heights, error
4022 4025 No newline at end of file
@@ -1,1 +1,1
1 <Project description="Segundo Test" id="191" name="test01"><ReadUnit datatype="SpectraReader" id="1911" inputId="0" name="SpectraReader"><Operation id="19111" name="run" priority="1" type="self"><Parameter format="str" id="191111" name="datatype" value="SpectraReader" /><Parameter format="str" id="191112" name="path" value="/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra" /><Parameter format="date" id="191113" name="startDate" value="2018/02/01" /><Parameter format="date" id="191114" name="endDate" value="2018/02/01" /><Parameter format="time" id="191115" name="startTime" value="17:00:00" /><Parameter format="time" id="191116" name="endTime" value="20:00:00" /><Parameter format="int" id="191118" name="online" value="0" /><Parameter format="int" id="191119" name="walk" value="1" /></Operation><Operation id="19112" name="printInfo" priority="2" type="self" /><Operation id="19113" name="printNumberOfBlock" priority="3" type="self" /></ReadUnit><ProcUnit datatype="Parameters" id="1913" inputId="1912" name="ParametersProc"><Operation id="19131" name="run" priority="1" type="self" /><Operation id="19132" name="SpectralFilters" priority="2" type="other"><Parameter format="float" id="191321" name="PositiveLimit" value="1.5" /><Parameter format="float" id="191322" name="NegativeLimit" value="3.5" /></Operation><Operation id="19133" name="FullSpectralAnalysis" priority="3" type="other"><Parameter format="float" id="191331" name="SNRlimit" value="-16" /><Parameter format="float" id="191332" name="E01" value="1.500" /><Parameter format="float" id="191333" name="E02" value="1.500" /><Parameter format="float" id="191334" name="E12" value="0" /><Parameter format="float" id="191335" name="N01" value="0.875" /><Parameter format="float" id="191336" name="N02" value="-0.875" /><Parameter format="float" id="191337" name="N12" value="-1.750" /></Operation><Operation id="19134" name="WindProfilerPlot" priority="4" type="other"><Parameter format="int" id="191341" name="id" value="4" /><Parameter format="str" id="191342" name="wintitle" value="Wind Profiler" /><Parameter format="float" id="191343" name="xmin" value="17" /><Parameter format="float" id="191344" name="xmax" value="20" /><Parameter format="float" id="191345" name="ymin" value="0" /><Parameter format="int" id="191346" name="ymax" value="11" /><Parameter format="float" id="191347" name="zmin" value="-20" /><Parameter format="float" id="191348" name="zmax" value="20" /><Parameter format="float" id="191349" name="SNRmin" value="-20" /><Parameter format="float" id="191350" name="SNRmax" value="20" /><Parameter format="float" id="191351" name="zmin_ver" value="-200" /><Parameter format="float" id="191352" name="zmax_ver" value="200" /><Parameter format="float" id="191353" name="SNRthresh" value="-20" /><Parameter format="int" id="191354" name="save" value="1" /><Parameter format="str" id="191355" name="figpath" value="/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra" /></Operation><Operation id="19135" name="PrecipitationProc" priority="5" type="other" /><Operation id="19136" name="ParametersPlot" priority="6" type="other"><Parameter format="int" id="191361" name="id" value="10" /><Parameter format="str" id="191362" name="wintitle" value="First_gg" /><Parameter format="str" id="191363" name="colormap" value="ocean_r" /><Parameter format="int" id="191364" name="zmin" value="50" /><Parameter format="int" id="191365" name="zmax" value="100" /><Parameter format="int" id="191366" name="ymin" value="0" /><Parameter format="int" id="191367" name="ymax" value="11" /><Parameter format="int" id="191368" name="xmin" value="17" /><Parameter format="int" id="191369" name="xmax" value="20" /><Parameter format="int" id="191370" name="save" value="1" /><Parameter format="str" id="191371" name="figpath" value="/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra" /></Operation><Operation id="19137" name="SpcParamPlot" priority="7" type="other"><Parameter format="int" id="191371" name="id" value="21" /><Parameter format="str" id="191372" name="wintitle" value="Primer eco removido" /><Parameter format="str" id="191373" name="xaxis" value="velocity" /><Parameter format="int" id="191374" name="showprofile" value="1" /><Parameter format="int" id="191375" name="zmin" value="10" /><Parameter format="int" id="191376" name="zmax" value="40" /><Parameter format="int" id="191377" name="Selector" value="1" /></Operation><Operation id="19138" name="ParamWriter" priority="8" type="other"><Parameter format="str" id="191381" name="path" value="/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdatatest/test1024" /><Parameter format="int" id="191382" name="blocksPerFile" value="100" /><Parameter format="list" id="191383" name="metadataList" value="heightList,timeZone,paramInterval" /><Parameter format="list" id="191384" name="dataList" value="data_output,data_SNR,utctime,utctimeInit" /></Operation></ProcUnit><ProcUnit datatype="SpectraProc" id="1912" inputId="1911" name="SpectraProc"><Operation id="19121" name="run" priority="1" type="self" /><Operation id="19122" name="setRadarFrequency" priority="2" type="self"><Parameter format="float" id="191221" name="frequency" value="445.09e6" /></Operation></ProcUnit></Project> No newline at end of file
1 <Project description="Segundo Test" id="191" name="test01"><ReadUnit datatype="SpectraReader" id="1911" inputId="0" name="SpectraReader"><Operation id="19111" name="run" priority="1" type="self"><Parameter format="str" id="191111" name="datatype" value="SpectraReader" /><Parameter format="str" id="191112" name="path" value="/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra" /><Parameter format="date" id="191113" name="startDate" value="2018/02/01" /><Parameter format="date" id="191114" name="endDate" value="2018/02/01" /><Parameter format="time" id="191115" name="startTime" value="17:00:00" /><Parameter format="time" id="191116" name="endTime" value="20:00:00" /><Parameter format="int" id="191118" name="online" value="0" /><Parameter format="int" id="191119" name="walk" value="1" /></Operation><Operation id="19112" name="printInfo" priority="2" type="self" /><Operation id="19113" name="printNumberOfBlock" priority="3" type="self" /></ReadUnit><ProcUnit datatype="Parameters" id="1913" inputId="1912" name="ParametersProc"><Operation id="19131" name="run" priority="1" type="self" /><Operation id="19132" name="SpectralFilters" priority="2" type="other"><Parameter format="float" id="191321" name="PositiveLimit" value="1.5" /><Parameter format="float" id="191322" name="NegativeLimit" value="3.5" /></Operation><Operation id="19133" name="PrecipitationProc" priority="3" type="other" /><Operation id="19134" name="ParametersPlot" priority="4" type="other"><Parameter format="int" id="191341" name="id" value="10" /><Parameter format="str" id="191342" name="wintitle" value="First_gg" /><Parameter format="str" id="191343" name="colormap" value="ocean_r" /><Parameter format="int" id="191344" name="zmin" value="00" /><Parameter format="int" id="191345" name="zmax" value="40" /><Parameter format="int" id="191346" name="ymin" value="0" /><Parameter format="int" id="191347" name="ymax" value="11" /><Parameter format="int" id="191348" name="xmin" value="17" /><Parameter format="int" id="191349" name="xmax" value="20" /><Parameter format="int" id="191350" name="save" value="1" /><Parameter format="str" id="191351" name="figpath" value="/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra" /></Operation><Operation id="19135" name="ParamWriter" priority="5" type="other"><Parameter format="str" id="191351" name="path" value="/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdatatest/test1024" /><Parameter format="int" id="191352" name="blocksPerFile" value="100" /><Parameter format="list" id="191353" name="metadataList" value="heightList,timeZone,paramInterval" /><Parameter format="list" id="191354" name="dataList" value="data_output,data_SNR,utctime,utctimeInit" /></Operation></ProcUnit><ProcUnit datatype="SpectraProc" id="1912" inputId="1911" name="SpectraProc"><Operation id="19121" name="run" priority="1" type="self" /><Operation id="19122" name="setRadarFrequency" priority="2" type="self"><Parameter format="float" id="191221" name="frequency" value="445.09e6" /></Operation></ProcUnit></Project> No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now