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