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