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