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