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