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