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