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