##// END OF EJS Templates
schain
RSS
Merge branch 'master' of http://jro-dev.igp.gob.pe/rhodecode/schain
jespinoza -
r1361:1dcf28d5b762 merge
parent child
Show More
Side by Side Unified
Context file:
r1361:1dcf28d5b762 Loading diff...:
Show file before | Show file after | Hide comments
@@ -68,6 +68,7 class BLTRParametersProc(ProcessingUnit):
68 68 SNRavgdB = 10*numpy.log10(SNRavg)
69 69 self.dataOut.data_snr_avg_db = SNRavgdB.reshape(1, *SNRavgdB.shape)
70 70
71 # Censoring Data
71 72 if snr_threshold is not None:
72 73 for i in range(3):
73 74 self.dataOut.data_param[i][SNRavgdB <= snr_threshold] = numpy.nan
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1164 lines changed) Show them Hide them
@@ -27,7 +27,6 import matplotlib.pyplot as plt
27 27
28 28 SPEED_OF_LIGHT = 299792458
29 29
30
31 30 '''solving pickling issue'''
32 31
33 32 def _pickle_method(method):
@@ -129,7 +128,7 class ParametersProc(ProcessingUnit):
129 128
130 129 if self.dataIn.type == "Spectra":
131 130
132 self.dataOut.data_pre = (self.dataIn.data_spc, self.dataIn.data_cspc)
131 self.dataOut.data_pre = [self.dataIn.data_spc, self.dataIn.data_cspc]
133 132 self.dataOut.data_spc = self.dataIn.data_spc
134 133 self.dataOut.data_cspc = self.dataIn.data_cspc
135 134 self.dataOut.nProfiles = self.dataIn.nProfiles
@@ -199,13 +198,11 def target(tups):
199 198
200 199 return obj.FitGau(args)
201 200
201 class RemoveWideGC(Operation):
202 ''' This class remove the wide clutter and replace it with a simple interpolation points
203 This mainly applies to CLAIRE radar
202 204
203 class SpectralFilters(Operation):
204
205 '''This class allows the Rainfall / Wind Selection for CLAIRE RADAR
206
207 LimitR : It is the limit in m/s of Rainfall
208 LimitW : It is the limit in m/s for Winds
205 ClutterWidth : Width to look for the clutter peak
209 206
210 207 Input:
211 208
@@ -215,91 +212,111 class SpectralFilters(Operation):
215 212 Affected:
216 213
217 214 self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
218 self.dataOut.spcparam_range : Used in SpcParamPlot
219 self.dataOut.SPCparam : Used in PrecipitationProc
220
221 215
216 Written by D. Scipión 25.02.2021
222 217 '''
223
224 218 def __init__(self):
225 219 Operation.__init__(self)
226 self.i=0
227
228 def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5):
229
230
231 #Limite de vientos
232 LimitR = PositiveLimit
233 LimitN = NegativeLimit
220 self.i = 0
221 self.ich = 0
222 self.ir = 0
223
224 def run(self, dataOut, ClutterWidth=2.5):
225 # print ('Entering RemoveWideGC ... ')
234 226
235 227 self.spc = dataOut.data_pre[0].copy()
236 self.cspc = dataOut.data_pre[1].copy()
237
238 self.Num_Hei = self.spc.shape[2]
239 self.Num_Bin = self.spc.shape[1]
228 self.spc_out = dataOut.data_pre[0].copy()
240 229 self.Num_Chn = self.spc.shape[0]
230 self.Num_Hei = self.spc.shape[2]
231 VelRange = dataOut.spc_range[2][:-1]
232 dv = VelRange[1]-VelRange[0]
233
234 # Find the velocities that corresponds to zero
235 gc_values = numpy.squeeze(numpy.where(numpy.abs(VelRange) <= ClutterWidth))
236
237 # Removing novalid data from the spectra
238 for ich in range(self.Num_Chn) :
239 for ir in range(self.Num_Hei) :
240 # Estimate the noise at each range
241 HSn = hildebrand_sekhon(self.spc[ich,:,ir],dataOut.nIncohInt)
242
243 # Removing the noise floor at each range
244 novalid = numpy.where(self.spc[ich,:,ir] < HSn)
245 self.spc[ich,novalid,ir] = HSn
246
247 junk = numpy.append(numpy.insert(numpy.squeeze(self.spc[ich,gc_values,ir]),0,HSn),HSn)
248 j1index = numpy.squeeze(numpy.where(numpy.diff(junk)>0))
249 j2index = numpy.squeeze(numpy.where(numpy.diff(junk)<0))
250 if ((numpy.size(j1index)<=1) | (numpy.size(j2index)<=1)) :
251 continue
252 junk3 = numpy.squeeze(numpy.diff(j1index))
253 junk4 = numpy.squeeze(numpy.diff(j2index))
254
255 valleyindex = j2index[numpy.where(junk4>1)]
256 peakindex = j1index[numpy.where(junk3>1)]
257
258 isvalid = numpy.squeeze(numpy.where(numpy.abs(VelRange[gc_values[peakindex]]) <= 2.5*dv))
259 if numpy.size(isvalid) == 0 :
260 continue
261 if numpy.size(isvalid) >1 :
262 vindex = numpy.argmax(self.spc[ich,gc_values[peakindex[isvalid]],ir])
263 isvalid = isvalid[vindex]
264
265 # clutter peak
266 gcpeak = peakindex[isvalid]
267 vl = numpy.where(valleyindex < gcpeak)
268 if numpy.size(vl) == 0:
269 continue
270 gcvl = valleyindex[vl[0][-1]]
271 vr = numpy.where(valleyindex > gcpeak)
272 if numpy.size(vr) == 0:
273 continue
274 gcvr = valleyindex[vr[0][0]]
275
276 # Removing the clutter
277 interpindex = numpy.array([gc_values[gcvl], gc_values[gcvr]])
278 gcindex = gc_values[gcvl+1:gcvr-1]
279 self.spc_out[ich,gcindex,ir] = numpy.interp(VelRange[gcindex],VelRange[interpindex],self.spc[ich,interpindex,ir])
280
281 dataOut.data_pre[0] = self.spc_out
282 #print ('Leaving RemoveWideGC ... ')
283 return dataOut
241 284
242 VelRange = dataOut.spc_range[2]
243 TimeRange = dataOut.spc_range[1]
244 FrecRange = dataOut.spc_range[0]
245
246 Vmax= 2*numpy.max(dataOut.spc_range[2])
247 Tmax= 2*numpy.max(dataOut.spc_range[1])
248 Fmax= 2*numpy.max(dataOut.spc_range[0])
249
250 Breaker1R=VelRange[numpy.abs(VelRange-(-LimitN)).argmin()]
251 Breaker1R=numpy.where(VelRange == Breaker1R)
252
253 Delta = self.Num_Bin/2 - Breaker1R[0]
254
255
256 '''Reacomodando SPCrange'''
257
258 VelRange=numpy.roll(VelRange,-(int(self.Num_Bin/2)) ,axis=0)
259
260 VelRange[-(int(self.Num_Bin/2)):]+= Vmax
285 class SpectralFilters(Operation):
286 ''' This class allows to replace the novalid values with noise for each channel
287 This applies to CLAIRE RADAR
261 288
262 FrecRange=numpy.roll(FrecRange,-(int(self.Num_Bin/2)),axis=0)
289 PositiveLimit : RightLimit of novalid data
290 NegativeLimit : LeftLimit of novalid data
263 291
264 FrecRange[-(int(self.Num_Bin/2)):]+= Fmax
292 Input:
265 293
266 TimeRange=numpy.roll(TimeRange,-(int(self.Num_Bin/2)),axis=0)
294 self.dataOut.data_pre : SPC and CSPC
295 self.dataOut.spc_range : To select wind and rainfall velocities
267 296
268 TimeRange[-(int(self.Num_Bin/2)):]+= Tmax
297 Affected:
269 298
270 ''' ------------------ '''
299 self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
271 300
272 Breaker2R=VelRange[numpy.abs(VelRange-(LimitR)).argmin()]
273 Breaker2R=numpy.where(VelRange == Breaker2R)
301 Written by D. Scipión 29.01.2021
302 '''
303 def __init__(self):
304 Operation.__init__(self)
305 self.i = 0
306
307 def run(self, dataOut, ):
274 308
309 self.spc = dataOut.data_pre[0].copy()
310 self.Num_Chn = self.spc.shape[0]
311 VelRange = dataOut.spc_range[2]
275 312
276 SPCroll = numpy.roll(self.spc,-(int(self.Num_Bin/2)) ,axis=1)
313 # novalid corresponds to data within the Negative and PositiveLimit
314
277 315
278 SPCcut = SPCroll.copy()
316 # Removing novalid data from the spectra
279 317 for i in range(self.Num_Chn):
280
281 SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i]
282 SPCcut[i,-int(Delta):,:] = dataOut.noise[i]
283
284 SPCcut[i]=SPCcut[i]- dataOut.noise[i]
285 SPCcut[ numpy.where( SPCcut<0 ) ] = 1e-20
286
287 SPCroll[i]=SPCroll[i]-dataOut.noise[i]
288 SPCroll[ numpy.where( SPCroll<0 ) ] = 1e-20
289
290 SPC_ch1 = SPCroll
291
292 SPC_ch2 = SPCcut
293
294 SPCparam = (SPC_ch1, SPC_ch2, self.spc)
295 dataOut.SPCparam = numpy.asarray(SPCparam)
296
297
298 dataOut.spcparam_range=numpy.zeros([self.Num_Chn,self.Num_Bin+1])
299
300 dataOut.spcparam_range[2]=VelRange
301 dataOut.spcparam_range[1]=TimeRange
302 dataOut.spcparam_range[0]=FrecRange
318 self.spc[i,novalid,:] = dataOut.noise[i]
319 dataOut.data_pre[0] = self.spc
303 320 return dataOut
304 321
305 322 class GaussianFit(Operation):
@@ -321,135 +338,186 class GaussianFit(Operation):
321 338 self.i=0
322 339
323 340
324 def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
341 # def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
342 def run(self, dataOut, SNRdBlimit=-9, method='generalized'):
325 343 """This routine will find a couple of generalized Gaussians to a power spectrum
344 methods: generalized, squared
326 345 input: spc
327 346 output:
328 Amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1,noise
347 noise, amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1
329 348 """
330
349 print ('Entering ',method,' double Gaussian fit')
331 350 self.spc = dataOut.data_pre[0].copy()
332 351 self.Num_Hei = self.spc.shape[2]
333 352 self.Num_Bin = self.spc.shape[1]
334 353 self.Num_Chn = self.spc.shape[0]
335 Vrange = dataOut.abscissaList
336
337 GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei])
338 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
339 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
340 SPC_ch1[:] = numpy.NaN
341 SPC_ch2[:] = numpy.NaN
342
343 354
344 355 start_time = time.time()
345 356
346 noise_ = dataOut.spc_noise[0].copy()
347
348
349 357 pool = Pool(processes=self.Num_Chn)
350 args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
358 args = [(dataOut.spc_range[2], ich, dataOut.spc_noise[ich], dataOut.nIncohInt, SNRdBlimit) for ich in range(self.Num_Chn)]
351 359 objs = [self for __ in range(self.Num_Chn)]
352 360 attrs = list(zip(objs, args))
353 gauSPC = pool.map(target, attrs)
354 dataOut.SPCparam = numpy.asarray(SPCparam)
355
356 ''' Parameters:
357 1. Amplitude
358 2. Shift
359 3. Width
360 4. Power
361 '''
361 DGauFitParam = pool.map(target, attrs)
362 # Parameters:
363 # 0. Noise, 1. Amplitude, 2. Shift, 3. Width 4. Power
364 dataOut.DGauFitParams = numpy.asarray(DGauFitParam)
365
366 # Double Gaussian Curves
367 gau0 = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
368 gau0[:] = numpy.NaN
369 gau1 = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
370 gau1[:] = numpy.NaN
371 x_mtr = numpy.transpose(numpy.tile(dataOut.getVelRange(1)[:-1], (self.Num_Hei,1)))
372 for iCh in range(self.Num_Chn):
373 N0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][0,:,0]] * self.Num_Bin))
374 N1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][0,:,1]] * self.Num_Bin))
375 A0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][1,:,0]] * self.Num_Bin))
376 A1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][1,:,1]] * self.Num_Bin))
377 v0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][2,:,0]] * self.Num_Bin))
378 v1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][2,:,1]] * self.Num_Bin))
379 s0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][3,:,0]] * self.Num_Bin))
380 s1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][3,:,1]] * self.Num_Bin))
381 if method == 'genealized':
382 p0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][4,:,0]] * self.Num_Bin))
383 p1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][4,:,1]] * self.Num_Bin))
384 elif method == 'squared':
385 p0 = 2.
386 p1 = 2.
387 gau0[iCh] = A0*numpy.exp(-0.5*numpy.abs((x_mtr-v0)/s0)**p0)+N0
388 gau1[iCh] = A1*numpy.exp(-0.5*numpy.abs((x_mtr-v1)/s1)**p1)+N1
389 dataOut.GaussFit0 = gau0
390 dataOut.GaussFit1 = gau1
391
392 print('Leaving ',method ,' double Gaussian fit')
393 return dataOut
362 394
363 395 def FitGau(self, X):
364
365 Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
366
367 SPCparam = []
368 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
369 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
370 SPC_ch1[:] = 0#numpy.NaN
371 SPC_ch2[:] = 0#numpy.NaN
372
373
374
396 # print('Entering FitGau')
397 # Assigning the variables
398 Vrange, ch, wnoise, num_intg, SNRlimit = X
399 # Noise Limits
400 noisebl = wnoise * 0.9
401 noisebh = wnoise * 1.1
402 # Radar Velocity
403 Va = max(Vrange)
404 deltav = Vrange[1] - Vrange[0]
405 x = numpy.arange(self.Num_Bin)
406
407 # print ('stop 0')
408
409 # 5 parameters, 2 Gaussians
410 DGauFitParam = numpy.zeros([5, self.Num_Hei,2])
411 DGauFitParam[:] = numpy.NaN
412
413 # SPCparam = []
414 # SPC_ch1 = numpy.zeros([self.Num_Bin,self.Num_Hei])
415 # SPC_ch2 = numpy.zeros([self.Num_Bin,self.Num_Hei])
416 # SPC_ch1[:] = 0 #numpy.NaN
417 # SPC_ch2[:] = 0 #numpy.NaN
418 # print ('stop 1')
375 419 for ht in range(self.Num_Hei):
376
377
420 # print (ht)
421 # print ('stop 2')
422 # Spectra at each range
378 423 spc = numpy.asarray(self.spc)[ch,:,ht]
424 snr = ( spc.mean() - wnoise ) / wnoise
425 snrdB = 10.*numpy.log10(snr)
379 426
427 #print ('stop 3')
428 if snrdB < SNRlimit :
429 # snr = numpy.NaN
430 # SPC_ch1[:,ht] = 0#numpy.NaN
431 # SPC_ch1[:,ht] = 0#numpy.NaN
432 # SPCparam = (SPC_ch1,SPC_ch2)
433 # print ('SNR less than SNRth')
434 continue
435 # wnoise = hildebrand_sekhon(spc,num_intg)
436 # print ('stop 2.01')
380 437 #############################################
381 438 # normalizing spc and noise
382 439 # This part differs from gg1
383 spc_norm_max = max(spc)
440 # spc_norm_max = max(spc) #commented by D. Scipión 19.03.2021
384 441 #spc = spc / spc_norm_max
385 pnoise = pnoise #/ spc_norm_max
442 # pnoise = pnoise #/ spc_norm_max #commented by D. Scipión 19.03.2021
386 443 #############################################
387 444
445 # print ('stop 2.1')
388 446 fatspectra=1.0
389
390 wnoise = noise_ #/ spc_norm_max
447 # noise per channel.... we might want to use the noise at each range
448
449 # wnoise = noise_ #/ spc_norm_max #commented by D. Scipión 19.03.2021
391 450 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
392 451 #if wnoise>1.1*pnoise: # to be tested later
393 452 # wnoise=pnoise
394 noisebl=wnoise*0.9;
395 noisebh=wnoise*1.1
396 spc=spc-wnoise
453 # noisebl = wnoise*0.9
454 # noisebh = wnoise*1.1
455 spc = spc - wnoise # signal
397 456
398 minx=numpy.argmin(spc)
457 # print ('stop 2.2')
458 minx = numpy.argmin(spc)
399 459 #spcs=spc.copy()
400 spcs=numpy.roll(spc,-minx)
401 cum=numpy.cumsum(spcs)
402 tot_noise=wnoise * self.Num_Bin #64;
403
404 snr = sum(spcs)/tot_noise
405 snrdB=10.*numpy.log10(snr)
406
407 if snrdB < SNRlimit :
408 snr = numpy.NaN
409 SPC_ch1[:,ht] = 0#numpy.NaN
410 SPC_ch1[:,ht] = 0#numpy.NaN
411 SPCparam = (SPC_ch1,SPC_ch2)
412 continue
460 spcs = numpy.roll(spc,-minx)
461 cum = numpy.cumsum(spcs)
462 # tot_noise = wnoise * self.Num_Bin #64;
463
464 # print ('stop 2.3')
465 # snr = sum(spcs) / tot_noise
466 # snrdB = 10.*numpy.log10(snr)
467 #print ('stop 3')
468 # if snrdB < SNRlimit :
469 # snr = numpy.NaN
470 # SPC_ch1[:,ht] = 0#numpy.NaN
471 # SPC_ch1[:,ht] = 0#numpy.NaN
472 # SPCparam = (SPC_ch1,SPC_ch2)
473 # print ('SNR less than SNRth')
474 # continue
413 475
414 476
415 477 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
416 478 # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
417
418 cummax=max(cum);
419 epsi=0.08*fatspectra # cumsum to narrow down the energy region
420 cumlo=cummax*epsi;
421 cumhi=cummax*(1-epsi)
422 powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
423
424
479 # print ('stop 4')
480 cummax = max(cum)
481 epsi = 0.08 * fatspectra # cumsum to narrow down the energy region
482 cumlo = cummax * epsi
483 cumhi = cummax * (1-epsi)
484 powerindex = numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
485
486 # print ('stop 5')
425 487 if len(powerindex) < 1:# case for powerindex 0
488 # print ('powerindex < 1')
426 489 continue
427 powerlo=powerindex[0]
428 powerhi=powerindex[-1]
429 powerwidth=powerhi-powerlo
430
431 firstpeak=powerlo+powerwidth/10.# first gaussian energy location
432 secondpeak=powerhi-powerwidth/10.#second gaussian energy location
433 midpeak=(firstpeak+secondpeak)/2.
434 firstamp=spcs[int(firstpeak)]
435 secondamp=spcs[int(secondpeak)]
436 midamp=spcs[int(midpeak)]
490 powerlo = powerindex[0]
491 powerhi = powerindex[-1]
492 powerwidth = powerhi-powerlo
493 if powerwidth <= 1:
494 # print('powerwidth <= 1')
495 continue
496
497 # print ('stop 6')
498 firstpeak = powerlo + powerwidth/10.# first gaussian energy location
499 secondpeak = powerhi - powerwidth/10. #second gaussian energy location
500 midpeak = (firstpeak + secondpeak)/2.
501 firstamp = spcs[int(firstpeak)]
502 secondamp = spcs[int(secondpeak)]
503 midamp = spcs[int(midpeak)]
437 504
438 x=numpy.arange( self.Num_Bin )
439 y_data=spc+wnoise
505 y_data = spc + wnoise
440 506
441 507 ''' single Gaussian '''
442 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
443 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
444 power0=2.
445 amplitude0=midamp
446 state0=[shift0,width0,amplitude0,power0,wnoise]
447 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
448 lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
449
450 chiSq1=lsq1[1];
451
452
508 shift0 = numpy.mod(midpeak+minx, self.Num_Bin )
509 width0 = powerwidth/4.#Initialization entire power of spectrum divided by 4
510 power0 = 2.
511 amplitude0 = midamp
512 state0 = [shift0,width0,amplitude0,power0,wnoise]
513 bnds = ((0,self.Num_Bin-1),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
514 lsq1 = fmin_l_bfgs_b(self.misfit1, state0, args=(y_data,x,num_intg), bounds=bnds, approx_grad=True)
515 # print ('stop 7.1')
516 # print (bnds)
517
518 chiSq1=lsq1[1]
519
520 # print ('stop 8')
453 521 if fatspectra<1.0 and powerwidth<4:
454 522 choice=0
455 523 Amplitude0=lsq1[0][2]
@@ -463,127 +531,142 class GaussianFit(Operation):
463 531 noise=lsq1[0][4]
464 532 #return (numpy.array([shift0,width0,Amplitude0,p0]),
465 533 # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
466
467 ''' two gaussians '''
534
535 # print ('stop 9')
536 ''' two Gaussians '''
468 537 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
469 shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
470 shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
471 width0=powerwidth/6.;
472 width1=width0
473 power0=2.;
474 power1=power0
475 amplitude0=firstamp;
476 amplitude1=secondamp
477 state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
538 shift0 = numpy.mod(firstpeak+minx, self.Num_Bin )
539 shift1 = numpy.mod(secondpeak+minx, self.Num_Bin )
540 width0 = powerwidth/6.
541 width1 = width0
542 power0 = 2.
543 power1 = power0
544 amplitude0 = firstamp
545 amplitude1 = secondamp
546 state0 = [shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
478 547 #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
479 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
548 bnds=((0,self.Num_Bin-1),(1,powerwidth/2.),(0,None),(0.5,3.),(0,self.Num_Bin-1),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
480 549 #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
481 550
551 # print ('stop 10')
482 552 lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )
483 553
554 # print ('stop 11')
555 chiSq2 = lsq2[1]
484 556
485 chiSq2=lsq2[1];
486
487
557 # print ('stop 12')
488 558
489 oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
559 oneG = (chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
490 560
561 # print ('stop 13')
491 562 if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
492 563 if oneG:
493 choice=0
564 choice = 0
494 565 else:
495 w1=lsq2[0][1]; w2=lsq2[0][5]
496 a1=lsq2[0][2]; a2=lsq2[0][6]
497 p1=lsq2[0][3]; p2=lsq2[0][7]
498 s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1;
499 s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
500 gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
566 w1 = lsq2[0][1]; w2 = lsq2[0][5]
567 a1 = lsq2[0][2]; a2 = lsq2[0][6]
568 p1 = lsq2[0][3]; p2 = lsq2[0][7]
569 s1 = (2**(1+1./p1))*scipy.special.gamma(1./p1)/p1
570 s2 = (2**(1+1./p2))*scipy.special.gamma(1./p2)/p2
571 gp1 = a1*w1*s1; gp2 = a2*w2*s2 # power content of each ggaussian with proper p scaling
501 572
502 573 if gp1>gp2:
503 574 if a1>0.7*a2:
504 choice=1
575 choice = 1
505 576 else:
506 choice=2
577 choice = 2
507 578 elif gp2>gp1:
508 579 if a2>0.7*a1:
509 choice=2
580 choice = 2
510 581 else:
511 choice=1
582 choice = 1
512 583 else:
513 choice=numpy.argmax([a1,a2])+1
584 choice = numpy.argmax([a1,a2])+1
514 585 #else:
515 586 #choice=argmin([std2a,std2b])+1
516 587
517 588 else: # with low SNR go to the most energetic peak
518 choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
519
520
521 shift0=lsq2[0][0];
522 vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
523 shift1=lsq2[0][4];
524 vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
525
526 max_vel = 1.0
527
589 choice = numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
590
591 # print ('stop 14')
592 shift0 = lsq2[0][0]
593 vel0 = Vrange[0] + shift0 * deltav
594 shift1 = lsq2[0][4]
595 # vel1=Vrange[0] + shift1 * deltav
596
597 # max_vel = 1.0
598 # Va = max(Vrange)
599 # deltav = Vrange[1]-Vrange[0]
600 # print ('stop 15')
528 601 #first peak will be 0, second peak will be 1
529 if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range
530 shift0=lsq2[0][0]
531 width0=lsq2[0][1]
532 Amplitude0=lsq2[0][2]
533 p0=lsq2[0][3]
534
535 shift1=lsq2[0][4]
536 width1=lsq2[0][5]
537 Amplitude1=lsq2[0][6]
538 p1=lsq2[0][7]
539 noise=lsq2[0][8]
602 # if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range # Commented by D.Scipión 19.03.2021
603 if vel0 > -Va and vel0 < Va : #first peak is in the correct range
604 shift0 = lsq2[0][0]
605 width0 = lsq2[0][1]
606 Amplitude0 = lsq2[0][2]
607 p0 = lsq2[0][3]
608
609 shift1 = lsq2[0][4]
610 width1 = lsq2[0][5]
611 Amplitude1 = lsq2[0][6]
612 p1 = lsq2[0][7]
613 noise = lsq2[0][8]
540 614 else:
541 shift1=lsq2[0][0]
542 width1=lsq2[0][1]
543 Amplitude1=lsq2[0][2]
544 p1=lsq2[0][3]
615 shift1 = lsq2[0][0]
616 width1 = lsq2[0][1]
617 Amplitude1 = lsq2[0][2]
618 p1 = lsq2[0][3]
545 619
546 shift0=lsq2[0][4]
547 width0=lsq2[0][5]
548 Amplitude0=lsq2[0][6]
549 p0=lsq2[0][7]
550 noise=lsq2[0][8]
620 shift0 = lsq2[0][4]
621 width0 = lsq2[0][5]
622 Amplitude0 = lsq2[0][6]
623 p0 = lsq2[0][7]
624 noise = lsq2[0][8]
551 625
552 626 if Amplitude0<0.05: # in case the peak is noise
553 shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN]
627 shift0,width0,Amplitude0,p0 = 4*[numpy.NaN]
554 628 if Amplitude1<0.05:
555 shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN]
556
557
558 SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
559 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
560 SPCparam = (SPC_ch1,SPC_ch2)
561
562
563 return GauSPC
629 shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
630
631 # print ('stop 16 ')
632 # SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0)/width0)**p0)
633 # SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1)/width1)**p1)
634 # SPCparam = (SPC_ch1,SPC_ch2)
635
636 DGauFitParam[0,ht,0] = noise
637 DGauFitParam[0,ht,1] = noise
638 DGauFitParam[1,ht,0] = Amplitude0
639 DGauFitParam[1,ht,1] = Amplitude1
640 DGauFitParam[2,ht,0] = Vrange[0] + shift0 * deltav
641 DGauFitParam[2,ht,1] = Vrange[0] + shift1 * deltav
642 DGauFitParam[3,ht,0] = width0 * deltav
643 DGauFitParam[3,ht,1] = width1 * deltav
644 DGauFitParam[4,ht,0] = p0
645 DGauFitParam[4,ht,1] = p1
646
647 # print (DGauFitParam.shape)
648 # print ('Leaving FitGau')
649 return DGauFitParam
650 # return SPCparam
651 # return GauSPC
564 652
565 653 def y_model1(self,x,state):
566 shift0,width0,amplitude0,power0,noise=state
567 model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
568
569 model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
570
571 model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
572 return model0+model0u+model0d+noise
654 shift0, width0, amplitude0, power0, noise = state
655 model0 = amplitude0*numpy.exp(-0.5*abs((x - shift0)/width0)**power0)
656 model0u = amplitude0*numpy.exp(-0.5*abs((x - shift0 - self.Num_Bin)/width0)**power0)
657 model0d = amplitude0*numpy.exp(-0.5*abs((x - shift0 + self.Num_Bin)/width0)**power0)
658 return model0 + model0u + model0d + noise
573 659
574 660 def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
575 shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state
576 model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
577
578 model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
579
580 model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
581 model1=amplitude1*numpy.exp(-0.5*abs((x-shift1)/width1)**power1)
582
583 model1u=amplitude1*numpy.exp(-0.5*abs((x-shift1- self.Num_Bin )/width1)**power1)
584
585 model1d=amplitude1*numpy.exp(-0.5*abs((x-shift1+ self.Num_Bin )/width1)**power1)
586 return model0+model0u+model0d+model1+model1u+model1d+noise
661 shift0, width0, amplitude0, power0, shift1, width1, amplitude1, power1, noise = state
662 model0 = amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
663 model0u = amplitude0*numpy.exp(-0.5*abs((x - shift0 - self.Num_Bin)/width0)**power0)
664 model0d = amplitude0*numpy.exp(-0.5*abs((x - shift0 + self.Num_Bin)/width0)**power0)
665
666 model1 = amplitude1*numpy.exp(-0.5*abs((x - shift1)/width1)**power1)
667 model1u = amplitude1*numpy.exp(-0.5*abs((x - shift1 - self.Num_Bin)/width1)**power1)
668 model1d = amplitude1*numpy.exp(-0.5*abs((x - shift1 + self.Num_Bin)/width1)**power1)
669 return model0 + model0u + model0d + model1 + model1u + model1d + noise
587 670
588 671 def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is.
589 672
@@ -614,31 +697,10 class PrecipitationProc(Operation):
614 697 Operation.__init__(self)
615 698 self.i=0
616 699
617
618 def gaus(self,xSamples,Amp,Mu,Sigma):
619 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
620
621
622
623 def Moments(self, ySamples, xSamples):
624 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
625 yNorm = ySamples / Pot
626
627 Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
628 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
629 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
630
631 return numpy.array([Pot, Vr, Desv])
632
633 700 def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118,
634 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km = 0.93, Altitude=3350):
635
701 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km2 = 0.93, Altitude=3350,SNRdBlimit=-30):
636 702
637 Velrange = dataOut.spcparam_range[2]
638 FrecRange = dataOut.spcparam_range[0]
639
640 dV= Velrange[1]-Velrange[0]
641 dF= FrecRange[1]-FrecRange[0]
703 # print ('Entering PrecepitationProc ... ')
642 704
643 705 if radar == "MIRA35C" :
644 706
@@ -650,18 +712,17 class PrecipitationProc(Operation):
650 712
651 713 else:
652 714
653 self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() #
654
655 """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
715 self.spc = dataOut.data_pre[0].copy()
656 716
717 #NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX
657 718 self.spc[:,:,0:7]= numpy.NaN
658 719
659 """##########################################"""
660
661 720 self.Num_Hei = self.spc.shape[2]
662 721 self.Num_Bin = self.spc.shape[1]
663 722 self.Num_Chn = self.spc.shape[0]
664 723
724 VelRange = dataOut.spc_range[2]
725
665 726 ''' Se obtiene la constante del RADAR '''
666 727
667 728 self.Pt = Pt
@@ -670,104 +731,73 class PrecipitationProc(Operation):
670 731 self.Lambda = Lambda
671 732 self.aL = aL
672 733 self.tauW = tauW
673 self.ThetaT = ThetaT
734 self.ThetaT = ThetaT
674 735 self.ThetaR = ThetaR
736 self.GSys = 10**(36.63/10) # Ganancia de los LNA 36.63 dB
737 self.lt = 10**(1.67/10) # Perdida en cables Tx 1.67 dB
738 self.lr = 10**(5.73/10) # Perdida en cables Rx 5.73 dB
675 739
676 740 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
677 741 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
678 742 RadarConstant = 10e-26 * Numerator / Denominator #
679
680 ''' ============================= '''
681
682 self.spc[0] = (self.spc[0]-dataOut.noise[0])
683 self.spc[1] = (self.spc[1]-dataOut.noise[1])
684 self.spc[2] = (self.spc[2]-dataOut.noise[2])
685
686 self.spc[ numpy.where(self.spc < 0)] = 0
687
688 SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise))
689 SPCmean[ numpy.where(SPCmean < 0)] = 0
690
691 ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei])
692 ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei])
693 ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei])
694
695 Pr = SPCmean[:,:]
696
697 VelMeteoro = numpy.mean(SPCmean,axis=0)
698
699 D_range = numpy.zeros([self.Num_Bin,self.Num_Hei])
700 SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei])
701 N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei])
702 V_mean = numpy.zeros(self.Num_Hei)
703 del_V = numpy.zeros(self.Num_Hei)
704 Z = numpy.zeros(self.Num_Hei)
705 Ze = numpy.zeros(self.Num_Hei)
706 RR = numpy.zeros(self.Num_Hei)
707
708 Range = dataOut.heightList*1000.
709
710 for R in range(self.Num_Hei):
711
712 h = Range[R] + Altitude #Range from ground to radar pulse altitude
713 del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
714
715 D_range[:,R] = numpy.log( (9.65 - (Velrange[0:self.Num_Bin] / del_V[R])) / 10.3 ) / -0.6 #Diameter range [m]x10**-3
716
717 '''NOTA: ETA(n) dn = ETA(f) df
718
719 dn = 1 Diferencial de muestreo
720 df = ETA(n) / ETA(f)
721
722 '''
723
724 ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA)
725
726 ETAv[:,R]=ETAn[:,R]/dV
727
728 ETAd[:,R]=ETAv[:,R]*6.18*exp(-0.6*D_range[:,R])
729
730 SIGMA[:,R] = Km * (D_range[:,R] * 1e-3 )**6 * numpy.pi**5 / Lambda**4 #Equivalent Section of drops (sigma)
731
732 N_dist[:,R] = ETAn[:,R] / SIGMA[:,R]
733
734 DMoments = self.Moments(Pr[:,R], Velrange[0:self.Num_Bin])
735
736 try:
737 popt01,pcov = curve_fit(self.gaus, Velrange[0:self.Num_Bin] , Pr[:,R] , p0=DMoments)
738 except:
739 popt01=numpy.zeros(3)
740 popt01[1]= DMoments[1]
741
742 if popt01[1]<0 or popt01[1]>20:
743 popt01[1]=numpy.NaN
744
745
746 V_mean[R]=popt01[1]
747
748 Z[R] = numpy.nansum( N_dist[:,R] * (D_range[:,R])**6 )#*10**-18
749
750 RR[R] = 0.0006*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate
751
752 Ze[R] = (numpy.nansum( ETAn[:,R]) * Lambda**4) / ( 10**-18*numpy.pi**5 * Km)
753
754
755
756 RR2 = (Z/200)**(1/1.6)
757 dBRR = 10*numpy.log10(RR)
758 dBRR2 = 10*numpy.log10(RR2)
759
760 dBZe = 10*numpy.log10(Ze)
761 dBZ = 10*numpy.log10(Z)
743 ExpConstant = 10**(40/10) #Constante Experimental
744
745 SignalPower = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
746 for i in range(self.Num_Chn):
747 SignalPower[i,:,:] = self.spc[i,:,:] - dataOut.noise[i]
748 SignalPower[numpy.where(SignalPower < 0)] = 1e-20
749
750 SPCmean = numpy.mean(SignalPower, 0)
751 Pr = SPCmean[:,:]/dataOut.normFactor
752
753 # Declaring auxiliary variables
754 Range = dataOut.heightList*1000. #Range in m
755 # replicate the heightlist to obtain a matrix [Num_Bin,Num_Hei]
756 rMtrx = numpy.transpose(numpy.transpose([dataOut.heightList*1000.] * self.Num_Bin))
757 zMtrx = rMtrx+Altitude
758 # replicate the VelRange to obtain a matrix [Num_Bin,Num_Hei]
759 VelMtrx = numpy.transpose(numpy.tile(VelRange[:-1], (self.Num_Hei,1)))
760
761 # height dependence to air density Foote and Du Toit (1969)
762 delv_z = 1 + 3.68e-5 * zMtrx + 1.71e-9 * zMtrx**2
763 VMtrx = VelMtrx / delv_z #Normalized velocity
764 VMtrx[numpy.where(VMtrx> 9.6)] = numpy.NaN
765 # Diameter is related to the fall speed of falling drops
766 D_Vz = -1.667 * numpy.log( 0.9369 - 0.097087 * VMtrx ) # D in [mm]
767 # Only valid for D>= 0.16 mm
768 D_Vz[numpy.where(D_Vz < 0.16)] = numpy.NaN
769
770 #Calculate Radar Reflectivity ETAn
771 ETAn = (RadarConstant *ExpConstant) * Pr * rMtrx**2 #Reflectivity (ETA)
772 ETAd = ETAn * 6.18 * exp( -0.6 * D_Vz ) * delv_z
773 # Radar Cross Section
774 sigmaD = Km2 * (D_Vz * 1e-3 )**6 * numpy.pi**5 / Lambda**4
775 # Drop Size Distribution
776 DSD = ETAn / sigmaD
777 # Equivalente Reflectivy
778 Ze_eqn = numpy.nansum( DSD * D_Vz**6 ,axis=0)
779 Ze_org = numpy.nansum(ETAn * Lambda**4, axis=0) / (1e-18*numpy.pi**5 * Km2) # [mm^6 /m^3]
780 # RainFall Rate
781 RR = 0.0006*numpy.pi * numpy.nansum( D_Vz**3 * DSD * VelMtrx ,0) #mm/hr
782
783 # Censoring the data
784 # Removing data with SNRth < 0dB se debe considerar el SNR por canal
785 SNRth = 10**(SNRdBlimit/10) #-30dB
786 novalid = numpy.where((dataOut.data_snr[0,:] <SNRth) | (dataOut.data_snr[1,:] <SNRth) | (dataOut.data_snr[2,:] <SNRth)) # AND condition. Maybe OR condition better
787 W = numpy.nanmean(dataOut.data_dop,0)
788 W[novalid] = numpy.NaN
789 Ze_org[novalid] = numpy.NaN
790 RR[novalid] = numpy.NaN
762 791
763 792 dataOut.data_output = RR[8]
764 793 dataOut.data_param = numpy.ones([3,self.Num_Hei])
765 794 dataOut.channelList = [0,1,2]
766
767 dataOut.data_param[0]=dBZ
768 dataOut.data_param[1]=V_mean
795
796 dataOut.data_param[0]=10*numpy.log10(Ze_org)
797 dataOut.data_param[1]=-W
769 798 dataOut.data_param[2]=RR
770 799
800 # print ('Leaving PrecepitationProc ... ')
771 801 return dataOut
772 802
773 803 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
@@ -784,7 +814,7 class PrecipitationProc(Operation):
784 814
785 815 ETA = numpy.sum(SNR,1)
786 816
787 ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
817 ETA = numpy.where(ETA != 0. , ETA, numpy.NaN)
788 818
789 819 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
790 820
@@ -832,29 +862,17 class FullSpectralAnalysis(Operation):
832 862
833 863 Output:
834 864
835 self.dataOut.data_output : Zonal wind, Meridional wind and Vertical wind
865 self.dataOut.data_output : Zonal wind, Meridional wind, and Vertical wind
836 866
837 867
838 868 Parameters affected: Winds, height range, SNR
839 869
840 870 """
841 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRlimit=7, minheight=None, maxheight=None):
842
843 self.indice=int(numpy.random.rand()*1000)
871 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRdBlimit=-30,
872 minheight=None, maxheight=None, NegativeLimit=None, PositiveLimit=None):
844 873
845 874 spc = dataOut.data_pre[0].copy()
846 875 cspc = dataOut.data_pre[1]
847
848 """Erick: NOTE THE RANGE OF THE PULSE TX MUST BE REMOVED"""
849
850 SNRspc = spc.copy()
851 SNRspc[:,:,0:7]= numpy.NaN
852
853 """##########################################"""
854
855
856 nChannel = spc.shape[0]
857 nProfiles = spc.shape[1]
858 876 nHeights = spc.shape[2]
859 877
860 878 # first_height = 0.75 #km (ref: data header 20170822)
@@ -881,119 +899,81 class FullSpectralAnalysis(Operation):
881 899 else:
882 900 ChanDist = numpy.array([[Xi01, Eta01],[Xi02,Eta02],[Xi12,Eta12]])
883 901
884 FrecRange = dataOut.spc_range[0]
885
886 data_SNR=numpy.zeros([nProfiles])
887 noise = dataOut.noise
888
889 dataOut.data_snr = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0]
890
891 dataOut.data_snr[numpy.where( dataOut.data_snr <0 )] = 1e-20
892
902 # 4 variables: zonal, meridional, vertical, and average SNR
903 data_param = numpy.zeros([4,nHeights]) * numpy.NaN
904 velocityX = numpy.zeros([nHeights]) * numpy.NaN
905 velocityY = numpy.zeros([nHeights]) * numpy.NaN
906 velocityZ = numpy.zeros([nHeights]) * numpy.NaN
893 907
894 data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN
895
896 velocityX=[]
897 velocityY=[]
898 velocityV=[]
899
900 dbSNR = 10*numpy.log10(dataOut.data_snr)
901 dbSNR = numpy.average(dbSNR,0)
908 dbSNR = 10*numpy.log10(numpy.average(dataOut.data_snr,0))
902 909
903 910 '''***********************************************WIND ESTIMATION**************************************'''
904
905 911 for Height in range(nHeights):
906 912
907 913 if Height >= range_min and Height < range_max:
908 # error_code unused, yet maybe useful for future analysis.
909 [Vzon,Vmer,Vver, error_code] = self.WindEstimation(spc[:,:,Height], cspc[:,:,Height], pairsList, ChanDist, Height, noise, dataOut.spc_range, dbSNR[Height], SNRlimit)
910 else:
911 Vzon,Vmer,Vver = 0., 0., numpy.NaN
912
913
914 if abs(Vzon) < 100. and abs(Vzon) > 0. and abs(Vmer) < 100. and abs(Vmer) > 0.:
915 velocityX=numpy.append(velocityX, Vzon)
916 velocityY=numpy.append(velocityY, -Vmer)
917
918 else:
919 velocityX=numpy.append(velocityX, numpy.NaN)
920 velocityY=numpy.append(velocityY, numpy.NaN)
921
922 if dbSNR[Height] > SNRlimit:
923 velocityV=numpy.append(velocityV, -Vver) # reason for this minus sign -> convention? (taken from Ericks version)
924 else:
925 velocityV=numpy.append(velocityV, numpy.NaN)
926
927
928 '''Change the numpy.array (velocityX) sign when trying to process BLTR data (Erick)'''
929 data_output[0] = numpy.array(velocityX)
930 data_output[1] = numpy.array(velocityY)
931 data_output[2] = velocityV
932
933
934 dataOut.data_output = data_output
935
914 # error_code will be useful in future analysis
915 [Vzon,Vmer,Vver, error_code] = self.WindEstimation(spc[:,:,Height], cspc[:,:,Height], pairsList,
916 ChanDist, Height, dataOut.noise, dataOut.spc_range, dbSNR[Height], SNRdBlimit, NegativeLimit, PositiveLimit,dataOut.frequency)
917
918 if abs(Vzon) < 100. and abs(Vmer) < 100.:
919 velocityX[Height] = Vzon
920 velocityY[Height] = -Vmer
921 velocityZ[Height] = Vver
922
923 # Censoring data with SNR threshold
924 dbSNR [dbSNR < SNRdBlimit] = numpy.NaN
925
926 data_param[0] = velocityX
927 data_param[1] = velocityY
928 data_param[2] = velocityZ
929 data_param[3] = dbSNR
930 dataOut.data_param = data_param
936 931 return dataOut
937 932
938
939 933 def moving_average(self,x, N=2):
940 934 """ convolution for smoothenig data. note that last N-1 values are convolution with zeroes """
941 935 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
942 936
943 937 def gaus(self,xSamples,Amp,Mu,Sigma):
944 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
938 return Amp * numpy.exp(-0.5*((xSamples - Mu)/Sigma)**2)
945 939
946 940 def Moments(self, ySamples, xSamples):
947 '''***
948 Variables corresponding to moments of distribution.
949 Also used as initial coefficients for curve_fit.
950 Vr was corrected. Only a velocity when x is velocity, of course.
951 ***'''
952 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
953 yNorm = ySamples / Pot
954 x_range = (numpy.max(xSamples)-numpy.min(xSamples))
955 Vr = numpy.nansum( yNorm * xSamples )*x_range/len(xSamples) # Velocidad radial, mu, corrimiento doppler, primer momento
956 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
957 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
958
959 return numpy.array([Pot, Vr, Desv])
941 Power = numpy.nanmean(ySamples) # Power, 0th Moment
942 yNorm = ySamples / numpy.nansum(ySamples)
943 RadVel = numpy.nansum(xSamples * yNorm) # Radial Velocity, 1st Moment
944 Sigma2 = numpy.nansum(yNorm * (xSamples - RadVel)**2) # Spectral Width, 2nd Moment
945 StdDev = numpy.sqrt(numpy.abs(Sigma2)) # Desv. Estandar, Ancho espectral
946 return numpy.array([Power,RadVel,StdDev])
960 947
961 948 def StopWindEstimation(self, error_code):
962 '''
963 the wind calculation and returns zeros
964 '''
965 Vzon = 0
966 Vmer = 0
967 Vver = numpy.nan
949 Vzon = numpy.NaN
950 Vmer = numpy.NaN
951 Vver = numpy.NaN
968 952 return Vzon, Vmer, Vver, error_code
969 953
970 954 def AntiAliasing(self, interval, maxstep):
971 955 """
972 956 function to prevent errors from aliased values when computing phaseslope
973 957 """
974 antialiased = numpy.zeros(len(interval))*0.0
958 antialiased = numpy.zeros(len(interval))
975 959 copyinterval = interval.copy()
976 960
977 961 antialiased[0] = copyinterval[0]
978 962
979 963 for i in range(1,len(antialiased)):
980
981 964 step = interval[i] - interval[i-1]
982
983 965 if step > maxstep:
984 966 copyinterval -= 2*numpy.pi
985 967 antialiased[i] = copyinterval[i]
986
987 968 elif step < maxstep*(-1):
988 969 copyinterval += 2*numpy.pi
989 970 antialiased[i] = copyinterval[i]
990
991 971 else:
992 972 antialiased[i] = copyinterval[i].copy()
993 973
994 974 return antialiased
995 975
996 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
976 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit, NegativeLimit, PositiveLimit, radfreq):
997 977 """
998 978 Function that Calculates Zonal, Meridional and Vertical wind velocities.
999 979 Initial Version by E. Bocanegra updated by J. Zibell until Nov. 2019.
@@ -1025,42 +1005,40 class FullSpectralAnalysis(Operation):
1025 1005
1026 1006 error_code = 0
1027 1007
1028
1029 SPC_Samples = numpy.ones([spc.shape[0],spc.shape[1]]) # for normalized spc values for one height
1030 phase = numpy.ones([spc.shape[0],spc.shape[1]]) # phase between channels
1031 CSPC_Samples = numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_) # for normalized cspc values
1032 PhaseSlope = numpy.zeros(spc.shape[0]) # slope of the phases, channelwise
1033 PhaseInter = numpy.ones(spc.shape[0]) # intercept to the slope of the phases, channelwise
1034 xFrec = AbbsisaRange[0][0:spc.shape[1]] # frequency range
1035 xVel = AbbsisaRange[2][0:spc.shape[1]] # velocity range
1036 SPCav = numpy.average(spc, axis=0)-numpy.average(noise) # spc[0]-noise[0]
1037
1038 SPCmoments_vel = self.Moments(SPCav, xVel ) # SPCmoments_vel[1] corresponds to vertical velocity and is used to determine if signal corresponds to wind (if .. <3)
1039 CSPCmoments = []
1040
1008 nChan = spc.shape[0]
1009 nProf = spc.shape[1]
1010 nPair = cspc.shape[0]
1011
1012 SPC_Samples = numpy.zeros([nChan, nProf]) # for normalized spc values for one height
1013 CSPC_Samples = numpy.zeros([nPair, nProf], dtype=numpy.complex_) # for normalized cspc values
1014 phase = numpy.zeros([nPair, nProf]) # phase between channels
1015 PhaseSlope = numpy.zeros(nPair) # slope of the phases, channelwise
1016 PhaseInter = numpy.zeros(nPair) # intercept to the slope of the phases, channelwise
1017 xFrec = AbbsisaRange[0][:-1] # frequency range
1018 xVel = AbbsisaRange[2][:-1] # velocity range
1019 xSamples = xFrec # the frequency range is taken
1020 delta_x = xSamples[1] - xSamples[0] # delta_f or delta_x
1021
1022 # only consider velocities with in NegativeLimit and PositiveLimit
1023 if (NegativeLimit is None):
1024 NegativeLimit = numpy.min(xVel)
1025 if (PositiveLimit is None):
1026 PositiveLimit = numpy.max(xVel)
1027 xvalid = numpy.where((xVel > NegativeLimit) & (xVel < PositiveLimit))
1028 xSamples_zoom = xSamples[xvalid]
1041 1029
1042 1030 '''Getting Eij and Nij'''
1043
1044 1031 Xi01, Xi02, Xi12 = ChanDist[:,0]
1045 1032 Eta01, Eta02, Eta12 = ChanDist[:,1]
1046 1033
1047 # update nov 19
1048 widthlimit = 7 # maximum width in Hz of the gaussian, empirically determined. Anything above 10 is unrealistic, often values between 1 and 5 correspond to proper fits.
1049
1034 # spwd limit - updated by D. Scipión 30.03.2021
1035 widthlimit = 10
1050 1036 '''************************* SPC is normalized ********************************'''
1051
1052 spc_norm = spc.copy() # need copy() because untouched spc is needed for normalization of cspc below
1053 spc_norm = numpy.where(numpy.isfinite(spc_norm), spc_norm, numpy.NAN)
1054
1055 for i in range(spc.shape[0]):
1056
1057 spc_sub = spc_norm[i,:] - noise[i] # spc not smoothed here or in previous version.
1058
1059 Factor_Norm = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc_sub)) # usually = Freq range / nfft
1060 normalized_spc = spc_sub / (numpy.nansum(numpy.abs(spc_sub)) * Factor_Norm)
1061
1062 xSamples = xFrec # the frequency range is taken
1063 SPC_Samples[i] = normalized_spc # Normalized SPC values are taken
1037 spc_norm = spc.copy()
1038 # For each channel
1039 for i in range(nChan):
1040 spc_sub = spc_norm[i,:] - noise[i] # only the signal power
1041 SPC_Samples[i] = spc_sub / (numpy.nansum(spc_sub) * delta_x)
1064 1042
1065 1043 '''********************** FITTING MEAN SPC GAUSSIAN **********************'''
1066 1044
@@ -1074,30 +1052,26 class FullSpectralAnalysis(Operation):
1074 1052 due to subtraction of the noise, some values are below zero. Raw "spc" values should be
1075 1053 >= 0, as it is the modulus squared of the signals (complex * it's conjugate)
1076 1054 """
1077
1078 SPCMean = numpy.average(SPC_Samples, axis=0)
1079
1080 popt = [1e-10,0,1e-10]
1081 SPCMoments = self.Moments(SPCMean, xSamples)
1082
1083 if dbSNR > SNRlimit and numpy.abs(SPCmoments_vel[1]) < 3:
1055 # initial conditions
1056 popt = [1e-10,0,1e-10]
1057 # Spectra average
1058 SPCMean = numpy.average(SPC_Samples,0)
1059 # Moments in frequency
1060 SPCMoments = self.Moments(SPCMean[xvalid], xSamples_zoom)
1061
1062 # Gauss Fit SPC in frequency domain
1063 if dbSNR > SNRlimit: # only if SNR > SNRth
1084 1064 try:
1085 popt,pcov = curve_fit(self.gaus,xSamples,SPCMean,p0=SPCMoments)#, bounds=(-numpy.inf, [numpy.inf, numpy.inf, 10])). Setting bounds does not make the code faster but only keeps the fit from finding the minimum.
1086 if popt[2] > widthlimit: # CONDITION
1065 popt,pcov = curve_fit(self.gaus,xSamples_zoom,SPCMean[xvalid],p0=SPCMoments)
1066 if popt[2] <= 0 or popt[2] > widthlimit: # CONDITION
1087 1067 return self.StopWindEstimation(error_code = 1)
1088
1089 FitGauss = self.gaus(xSamples,*popt)
1090
1068 FitGauss = self.gaus(xSamples_zoom,*popt)
1091 1069 except :#RuntimeError:
1092 1070 return self.StopWindEstimation(error_code = 2)
1093
1094 1071 else:
1095 1072 return self.StopWindEstimation(error_code = 3)
1096 1073
1097
1098
1099 1074 '''***************************** CSPC Normalization *************************
1100 new section:
1101 1075 The Spc spectra are used to normalize the crossspectra. Peaks from precipitation
1102 1076 influence the norm which is not desired. First, a range is identified where the
1103 1077 wind peak is estimated -> sum_wind is sum of those frequencies. Next, the area
@@ -1109,159 +1083,82 class FullSpectralAnalysis(Operation):
1109 1083
1110 1084 A norm is found according to Briggs 92.
1111 1085 '''
1112
1113 radarWavelength = 0.6741 # meters
1114 count_limit_freq = numpy.abs(popt[1]) + widthlimit # Hz, m/s can be also used if velocity is desired abscissa.
1115 # count_limit_freq = numpy.max(xFrec)
1116
1117 channel_integrals = numpy.zeros(3)
1118
1119 for i in range(spc.shape[0]):
1120 '''
1121 find the point in array corresponding to count_limit frequency.
1122 sum over all frequencies in the range around zero Hz @ math.ceil(N_freq/2)
1123 '''
1124 N_freq = numpy.count_nonzero(~numpy.isnan(spc[i,:]))
1125 count_limit_int = int(math.ceil( count_limit_freq / numpy.max(xFrec) * (N_freq / 2) )) # gives integer point
1126 sum_wind = numpy.nansum( spc[i, (math.ceil(N_freq/2) - count_limit_int) : (math.ceil(N_freq / 2) + count_limit_int)] ) #N_freq/2 is where frequency (velocity) is zero, i.e. middle of spectrum.
1127 sum_noise = (numpy.mean(spc[i, :4]) + numpy.mean(spc[i, -6:-2]))/2.0 * (N_freq - 2*count_limit_int)
1128 channel_integrals[i] = (sum_noise + sum_wind) * (2*numpy.max(xFrec) / N_freq)
1129
1130
1131 cross_integrals_peak = numpy.zeros(3)
1132 # cross_integrals_totalrange = numpy.zeros(3)
1133
1134 for i in range(spc.shape[0]):
1135
1136 cspc_norm = cspc[i,:].copy() # cspc not smoothed here or in previous version
1137
1086 # for each pair
1087 for i in range(nPair):
1088 cspc_norm = cspc[i,:].copy()
1138 1089 chan_index0 = pairsList[i][0]
1139 1090 chan_index1 = pairsList[i][1]
1140
1141 cross_integrals_peak[i] = channel_integrals[chan_index0]*channel_integrals[chan_index1]
1142 normalized_cspc = cspc_norm / numpy.sqrt(cross_integrals_peak[i])
1143 CSPC_Samples[i] = normalized_cspc
1144
1145 ''' Finding cross integrals without subtracting any peaks:'''
1146 # FactorNorm0 = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc[chan_index0,:]))
1147 # FactorNorm1 = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc[chan_index1,:]))
1148 # cross_integrals_totalrange[i] = (numpy.nansum(spc[chan_index0,:])) * FactorNorm0 * (numpy.nansum(spc[chan_index1,:])) * FactorNorm1
1149 # normalized_cspc = cspc_norm / numpy.sqrt(cross_integrals_totalrange[i])
1150 # CSPC_Samples[i] = normalized_cspc
1151
1152
1091 CSPC_Samples[i] = cspc_norm / (numpy.sqrt(numpy.nansum(spc_norm[chan_index0])*numpy.nansum(spc_norm[chan_index1])) * delta_x)
1153 1092 phase[i] = numpy.arctan2(CSPC_Samples[i].imag, CSPC_Samples[i].real)
1154 1093
1094 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPC_Samples[0,xvalid]), xSamples_zoom),
1095 self.Moments(numpy.abs(CSPC_Samples[1,xvalid]), xSamples_zoom),
1096 self.Moments(numpy.abs(CSPC_Samples[2,xvalid]), xSamples_zoom)])
1155 1097
1156 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPC_Samples[0]), xSamples),
1157 self.Moments(numpy.abs(CSPC_Samples[1]), xSamples),
1158 self.Moments(numpy.abs(CSPC_Samples[2]), xSamples)])
1159
1160
1161 '''***Sorting out NaN entries***'''
1162 CSPCMask01 = numpy.abs(CSPC_Samples[0])
1163 CSPCMask02 = numpy.abs(CSPC_Samples[1])
1164 CSPCMask12 = numpy.abs(CSPC_Samples[2])
1165
1166 mask01 = ~numpy.isnan(CSPCMask01)
1167 mask02 = ~numpy.isnan(CSPCMask02)
1168 mask12 = ~numpy.isnan(CSPCMask12)
1169
1170 CSPCMask01 = CSPCMask01[mask01]
1171 CSPCMask02 = CSPCMask02[mask02]
1172 CSPCMask12 = CSPCMask12[mask12]
1173
1174
1175 popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
1176 FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0
1098 popt01, popt02, popt12 = [1e-10,0,1e-10], [1e-10,0,1e-10] ,[1e-10,0,1e-10]
1099 FitGauss01, FitGauss02, FitGauss12 = numpy.zeros(len(xSamples)), numpy.zeros(len(xSamples)), numpy.zeros(len(xSamples))
1177 1100
1178 1101 '''*******************************FIT GAUSS CSPC************************************'''
1179
1180 1102 try:
1181 popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
1103 popt01,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[0][xvalid]),p0=CSPCmoments[0])
1182 1104 if popt01[2] > widthlimit: # CONDITION
1183 1105 return self.StopWindEstimation(error_code = 4)
1184
1185 popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1])
1106 popt02,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[1][xvalid]),p0=CSPCmoments[1])
1186 1107 if popt02[2] > widthlimit: # CONDITION
1187 1108 return self.StopWindEstimation(error_code = 4)
1188
1189 popt12,pcov = curve_fit(self.gaus,xSamples[mask12],numpy.abs(CSPCMask12),p0=CSPCmoments[2])
1109 popt12,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[2][xvalid]),p0=CSPCmoments[2])
1190 1110 if popt12[2] > widthlimit: # CONDITION
1191 1111 return self.StopWindEstimation(error_code = 4)
1192 1112
1193 FitGauss01 = self.gaus(xSamples, *popt01)
1194 FitGauss02 = self.gaus(xSamples, *popt02)
1195 FitGauss12 = self.gaus(xSamples, *popt12)
1196
1113 FitGauss01 = self.gaus(xSamples_zoom, *popt01)
1114 FitGauss02 = self.gaus(xSamples_zoom, *popt02)
1115 FitGauss12 = self.gaus(xSamples_zoom, *popt12)
1197 1116 except:
1198 1117 return self.StopWindEstimation(error_code = 5)
1199 1118
1200 1119
1201 1120 '''************* Getting Fij ***************'''
1202
1203
1204 #Punto en Eje X de la Gaussiana donde se encuentra el centro -- x-axis point of the gaussian where the center is located
1205 # -> PointGauCenter
1121 # x-axis point of the gaussian where the center is located from GaussFit of spectra
1206 1122 GaussCenter = popt[1]
1207 ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()]
1208 PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0]
1123 ClosestCenter = xSamples_zoom[numpy.abs(xSamples_zoom-GaussCenter).argmin()]
1124 PointGauCenter = numpy.where(xSamples_zoom==ClosestCenter)[0][0]
1209 1125
1210 #Punto e^-1 hubicado en la Gaussiana -- point where e^-1 is located in the gaussian
1126 # Point where e^-1 is located in the gaussian
1211 1127 PeMinus1 = numpy.max(FitGauss) * numpy.exp(-1)
1212 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss"
1128 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # The closest point to"Peminus1" in "FitGauss"
1213 1129 PointFij = numpy.where(FitGauss==FijClosest)[0][0]
1214
1215 Fij = numpy.abs(xSamples[PointFij] - xSamples[PointGauCenter])
1130 Fij = numpy.abs(xSamples_zoom[PointFij] - xSamples_zoom[PointGauCenter])
1216 1131
1217 1132 '''********** Taking frequency ranges from mean SPCs **********'''
1218
1219 #GaussCenter = popt[1] #Primer momento 01
1220 GauWidth = popt[2] * 3/2 #Ancho de banda de Gau01 -- Bandwidth of Gau01 TODO why *3/2?
1133 GauWidth = popt[2] * 3/2 # Bandwidth of Gau01
1221 1134 Range = numpy.empty(2)
1222 1135 Range[0] = GaussCenter - GauWidth
1223 1136 Range[1] = GaussCenter + GauWidth
1224 #Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max) -- Point in x-axis where the bandwidth is located (min:max)
1225 ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()]
1226 ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()]
1227
1228 PointRangeMin = numpy.where(xSamples==ClosRangeMin)[0][0]
1229 PointRangeMax = numpy.where(xSamples==ClosRangeMax)[0][0]
1230
1137 # Point in x-axis where the bandwidth is located (min:max)
1138 ClosRangeMin = xSamples_zoom[numpy.abs(xSamples_zoom-Range[0]).argmin()]
1139 ClosRangeMax = xSamples_zoom[numpy.abs(xSamples_zoom-Range[1]).argmin()]
1140 PointRangeMin = numpy.where(xSamples_zoom==ClosRangeMin)[0][0]
1141 PointRangeMax = numpy.where(xSamples_zoom==ClosRangeMax)[0][0]
1231 1142 Range = numpy.array([ PointRangeMin, PointRangeMax ])
1232
1233 FrecRange = xFrec[ Range[0] : Range[1] ]
1234
1143 FrecRange = xSamples_zoom[ Range[0] : Range[1] ]
1235 1144
1236 1145 '''************************** Getting Phase Slope ***************************'''
1237
1238 for i in range(1,3): # Changed to only compute two
1239
1240 if len(FrecRange) > 5 and len(FrecRange) < spc.shape[1] * 0.3:
1241 # PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=1) #used before to smooth phase with N=3
1242 PhaseRange = phase[i,Range[0]:Range[1]].copy()
1243
1146 for i in range(nPair):
1147 if len(FrecRange) > 5:
1148 PhaseRange = phase[i, xvalid[0][Range[0]:Range[1]]].copy()
1244 1149 mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
1245
1246
1247 1150 if len(FrecRange) == len(PhaseRange):
1248
1249 1151 try:
1250 1152 slope, intercept, _, _, _ = stats.linregress(FrecRange[mask], self.AntiAliasing(PhaseRange[mask], 4.5))
1251 1153 PhaseSlope[i] = slope
1252 1154 PhaseInter[i] = intercept
1253
1254 1155 except:
1255 1156 return self.StopWindEstimation(error_code = 6)
1256
1257 1157 else:
1258 1158 return self.StopWindEstimation(error_code = 7)
1259
1260 1159 else:
1261 1160 return self.StopWindEstimation(error_code = 8)
1262 1161
1263
1264
1265 1162 '''*** Constants A-H correspond to the convention as in Briggs and Vincent 1992 ***'''
1266 1163
1267 1164 '''Getting constant C'''
@@ -1269,9 +1166,12 class FullSpectralAnalysis(Operation):
1269 1166
1270 1167 '''****** Getting constants F and G ******'''
1271 1168 MijEijNij = numpy.array([[Xi02,Eta02], [Xi12,Eta12]])
1272 MijResult0 = (-PhaseSlope[1] * cC) / (2*numpy.pi)
1273 MijResult1 = (-PhaseSlope[2] * cC) / (2*numpy.pi)
1274 MijResults = numpy.array([MijResult0,MijResult1])
1169 # MijEijNij = numpy.array([[Xi01,Eta01], [Xi02,Eta02], [Xi12,Eta12]])
1170 # MijResult0 = (-PhaseSlope[0] * cC) / (2*numpy.pi)
1171 MijResult1 = (-PhaseSlope[1] * cC) / (2*numpy.pi)
1172 MijResult2 = (-PhaseSlope[2] * cC) / (2*numpy.pi)
1173 # MijResults = numpy.array([MijResult0, MijResult1, MijResult2])
1174 MijResults = numpy.array([MijResult1, MijResult2])
1275 1175 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1276 1176
1277 1177 '''****** Getting constants A, B and H ******'''
@@ -1279,39 +1179,22 class FullSpectralAnalysis(Operation):
1279 1179 W02 = numpy.nanmax( FitGauss02 )
1280 1180 W12 = numpy.nanmax( FitGauss12 )
1281 1181
1282 WijResult0 = ((cF * Xi01 + cG * Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi / cC))
1283 WijResult1 = ((cF * Xi02 + cG * Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi / cC))
1284 WijResult2 = ((cF * Xi12 + cG * Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi / cC))
1285
1286 WijResults = numpy.array([WijResult0, WijResult1, WijResult2])
1182 WijResult01 = ((cF * Xi01 + cG * Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi / cC))
1183 WijResult02 = ((cF * Xi02 + cG * Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi / cC))
1184 WijResult12 = ((cF * Xi12 + cG * Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi / cC))
1185 WijResults = numpy.array([WijResult01, WijResult02, WijResult12])
1287 1186
1288 1187 WijEijNij = numpy.array([ [Xi01**2, Eta01**2, 2*Xi01*Eta01] , [Xi02**2, Eta02**2, 2*Xi02*Eta02] , [Xi12**2, Eta12**2, 2*Xi12*Eta12] ])
1289 1188 (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1290 1189
1291 1190 VxVy = numpy.array([[cA,cH],[cH,cB]])
1292 1191 VxVyResults = numpy.array([-cF,-cG])
1293 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1294
1295 Vzon = Vy
1296 Vmer = Vx
1297
1298 # Vmag=numpy.sqrt(Vzon**2+Vmer**2) # unused
1299 # Vang=numpy.arctan2(Vmer,Vzon) # unused
1300
1301
1302 ''' using frequency as abscissa. Due to three channels, the offzenith angle is zero
1303 and Vrad equal to Vver. formula taken from Briggs 92, figure 4.
1304 '''
1305 if numpy.abs( popt[1] ) < 3.5 and len(FrecRange) > 4:
1306 Vver = 0.5 * radarWavelength * popt[1] * 100 # *100 to get cm (/s)
1307 else:
1308 Vver = numpy.NaN
1309
1192 (Vmer,Vzon) = numpy.linalg.solve(VxVy, VxVyResults)
1193 Vver = -SPCMoments[1]*SPEED_OF_LIGHT/(2*radfreq)
1310 1194 error_code = 0
1311 1195
1312 1196 return Vzon, Vmer, Vver, error_code
1313 1197
1314
1315 1198 class SpectralMoments(Operation):
1316 1199
1317 1200 '''
@@ -1393,13 +1276,13 class SpectralMoments(Operation):
1393 1276 max_spec = aux.max()
1394 1277 m = aux.tolist().index(max_spec)
1395 1278
1396 #Smooth
1279 # Smooth
1397 1280 if (smooth == 0):
1398 1281 spec2 = spec
1399 1282 else:
1400 1283 spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
1401 1284
1402 # Calculo de Momentos
1285 # Moments Estimation
1403 1286 bb = spec2[numpy.arange(m,spec2.size)]
1404 1287 bb = (bb<n0).nonzero()
1405 1288 bb = bb[0]
@@ -1425,14 +1308,17 class SpectralMoments(Operation):
1425 1308
1426 1309 valid = numpy.arange(int(m + bb0 - ss1 + 1)) + ss1
1427 1310
1428 power = ((spec2[valid] - n0) * fwindow[valid]).sum()
1311 signal_power = ((spec2[valid] - n0) * fwindow[valid]).mean() # D. Scipión added with correct definition
1312 total_power = (spec2[valid] * fwindow[valid]).mean() # D. Scipión added with correct definition
1313 power = ((spec2[valid] - n0) * fwindow[valid]).sum()
1429 1314 fd = ((spec2[valid]- n0)*freq[valid] * fwindow[valid]).sum() / power
1430 1315 w = numpy.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum() / power)
1431 1316 snr = (spec2.mean()-n0)/n0
1432 1317 if (snr < 1.e-20) :
1433 1318 snr = 1.e-20
1434 1319
1435 vec_power[ind] = power
1320 # vec_power[ind] = power #D. Scipión replaced with the line below
1321 vec_power[ind] = total_power
1436 1322 vec_fd[ind] = fd
1437 1323 vec_w[ind] = w
1438 1324 vec_snr[ind] = snr
General Comments 0
You need to be logged in to leave comments. Login now