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