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