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