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