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