##// END OF EJS Templates
schain
RSS
Actualización del procesamiento de CLAIRE
Danny Scipión -
r1356:218999bf0aa4
parent child
Show More
Side by Side Unified
Context file:
r1356:218999bf0aa4 Loading diff...:
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (880 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):
@@ -223,56 +222,64 class RemoveWideGC(Operation):
223 222 self.ir = 0
224 223
225 224 def run(self, dataOut, ClutterWidth=2.5):
225 # print ('Entering RemoveWideGC ... ')
226 226
227 227 self.spc = dataOut.data_pre[0].copy()
228 228 self.spc_out = dataOut.data_pre[0].copy()
229 229 self.Num_Chn = self.spc.shape[0]
230 self.Num_Hei = self.spc.shape[1]
230 self.Num_Hei = self.spc.shape[2]
231 231 VelRange = dataOut.spc_range[2][:-1]
232 232 dv = VelRange[1]-VelRange[0]
233 233
234 234 # Find the velocities that corresponds to zero
235 gc_values = numpy.where(numpy.abs(VelRange) <= ClutterWidth)
235 gc_values = numpy.squeeze(numpy.where(numpy.abs(VelRange) <= ClutterWidth))
236 236
237 237 # Removing novalid data from the spectra
238 for ich in range(self.Num_Chn):
239 # Estimate the noise at aech range
240
241 for ir in range(self.Num_Hei):
242 # Estimate the noise at aech range
238 for ich in range(self.Num_Chn) :
239 for ir in range(self.Num_Hei) :
240 # Estimate the noise at each range
243 241 HSn = hildebrand_sekhon(self.spc[ich,:,ir],dataOut.nIncohInt)
242
244 243 # Removing the noise floor at each range
245 244 novalid = numpy.where(self.spc[ich,:,ir] < HSn)
246 self.spc[novalid,ir] = HSn
247
248 junk = [HSn, numpy.transpose(self.spc[ich,gc_values,ir]), HSn]
249 j1index = numpy.where(numpy.diff[junk]>0)
250 j2index = numpy.where(numpy.diff[junk]<0)
251 junk3 = numpy.diff(j1index)
252 junk4 = numpy.diff(j2index)
253
254 valleyindex = j2index[junk4>1]
255 peakindex = j1index[junk3>1]
256
257 # Identify the clutter (close to zero)
258 isvalid = numpy.where(where.abs(VelRange[gc_values[peakindex]]) <= 2.5*dv)
259 # if isempty(isvalid)
260 # continue
261 # if numel(isvalid)>1
262 # [~,vindex]= numpy.max(spc[gc_values[peakindex[isvalid]],ir])
263 # isvalid = isvalid[vindex]
245 self.spc[ich,novalid,ir] = HSn
246
247 junk = numpy.append(numpy.insert(numpy.squeeze(self.spc[ich,gc_values,ir]),0,HSn),HSn)
248 j1index = numpy.squeeze(numpy.where(numpy.diff(junk)>0))
249 j2index = numpy.squeeze(numpy.where(numpy.diff(junk)<0))
250 if ((numpy.size(j1index)<=1) | (numpy.size(j2index)<=1)) :
251 continue
252 junk3 = numpy.squeeze(numpy.diff(j1index))
253 junk4 = numpy.squeeze(numpy.diff(j2index))
254
255 valleyindex = j2index[numpy.where(junk4>1)]
256 peakindex = j1index[numpy.where(junk3>1)]
257
258 isvalid = numpy.squeeze(numpy.where(numpy.abs(VelRange[gc_values[peakindex]]) <= 2.5*dv))
259 if numpy.size(isvalid) == 0 :
260 continue
261 if numpy.size(isvalid) >1 :
262 vindex = numpy.argmax(self.spc[ich,gc_values[peakindex[isvalid]],ir])
263 isvalid = isvalid[vindex]
264
264 265 # clutter peak
265 266 gcpeak = peakindex[isvalid]
266 # left and right index of the clutter
267 gcvl = valleyindex[numpy.where(valleyindex < gcpeak, 1, 'last' )]
268 gcvr = valleyindex[numpy.where(valleyindex > gcpeak, 1)]
267 vl = numpy.where(valleyindex < gcpeak)
268 if numpy.size(vl) == 0:
269 continue
270 gcvl = valleyindex[vl[0][-1]]
271 vr = numpy.where(valleyindex > gcpeak)
272 if numpy.size(vr) == 0:
273 continue
274 gcvr = valleyindex[vr[0][0]]
269 275
270 276 # Removing the clutter
271 interpindex = [gc_values(gcvl), gc_values(gcvr)]
277 interpindex = numpy.array([gc_values[gcvl], gc_values[gcvr]])
272 278 gcindex = gc_values[gcvl+1:gcvr-1]
273 279 self.spc_out[ich,gcindex,ir] = numpy.interp(VelRange[gcindex],VelRange[interpindex],self.spc[ich,interpindex,ir])
274 280
275 281 dataOut.data_pre[0] = self.spc_out
282 #print ('Leaving RemoveWideGC ... ')
276 283 return dataOut
277 284
278 285 class SpectralFilters(Operation):
@@ -297,14 +304,14 class SpectralFilters(Operation):
297 304 Operation.__init__(self)
298 305 self.i = 0
299 306
300 def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=-1.5):
307 def run(self, dataOut, ):
301 308
302 309 self.spc = dataOut.data_pre[0].copy()
303 310 self.Num_Chn = self.spc.shape[0]
304 311 VelRange = dataOut.spc_range[2]
305 312
306 313 # novalid corresponds to data within the Negative and PositiveLimit
307 novalid = numpy.where((VelRange[:-1] >= NegativeLimit) & (VelRange[:-1] <= PositiveLimit))
314
308 315
309 316 # Removing novalid data from the spectra
310 317 for i in range(self.Num_Chn):
@@ -331,135 +338,186 class GaussianFit(Operation):
331 338 self.i=0
332 339
333 340
334 def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
341 # def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
342 def run(self, dataOut, SNRdBlimit=-9, method='generalized'):
335 343 """This routine will find a couple of generalized Gaussians to a power spectrum
344 methods: generalized, squared
336 345 input: spc
337 346 output:
338 Amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1,noise
347 noise, amplitude0,shift0,width0,p0,Amplitude1,shift1,width1,p1
339 348 """
340
349 print ('Entering ',method,' double Gaussian fit')
341 350 self.spc = dataOut.data_pre[0].copy()
342 351 self.Num_Hei = self.spc.shape[2]
343 352 self.Num_Bin = self.spc.shape[1]
344 353 self.Num_Chn = self.spc.shape[0]
345 Vrange = dataOut.abscissaList
346
347 GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei])
348 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
349 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
350 SPC_ch1[:] = numpy.NaN
351 SPC_ch2[:] = numpy.NaN
352
353 354
354 355 start_time = time.time()
355 356
356 noise_ = dataOut.spc_noise[0].copy()
357
358
359 357 pool = Pool(processes=self.Num_Chn)
360 args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
358 args = [(dataOut.spc_range[2], ich, dataOut.spc_noise[ich], dataOut.nIncohInt, SNRdBlimit) for ich in range(self.Num_Chn)]
361 359 objs = [self for __ in range(self.Num_Chn)]
362 360 attrs = list(zip(objs, args))
363 gauSPC = pool.map(target, attrs)
364 dataOut.SPCparam = numpy.asarray(SPCparam)
365
366 ''' Parameters:
367 1. Amplitude
368 2. Shift
369 3. Width
370 4. Power
371 '''
361 DGauFitParam = pool.map(target, attrs)
362 # Parameters:
363 # 0. Noise, 1. Amplitude, 2. Shift, 3. Width 4. Power
364 dataOut.DGauFitParams = numpy.asarray(DGauFitParam)
365
366 # Double Gaussian Curves
367 gau0 = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
368 gau0[:] = numpy.NaN
369 gau1 = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
370 gau1[:] = numpy.NaN
371 x_mtr = numpy.transpose(numpy.tile(dataOut.getVelRange(1)[:-1], (self.Num_Hei,1)))
372 for iCh in range(self.Num_Chn):
373 N0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][0,:,0]] * self.Num_Bin))
374 N1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][0,:,1]] * self.Num_Bin))
375 A0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][1,:,0]] * self.Num_Bin))
376 A1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][1,:,1]] * self.Num_Bin))
377 v0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][2,:,0]] * self.Num_Bin))
378 v1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][2,:,1]] * self.Num_Bin))
379 s0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][3,:,0]] * self.Num_Bin))
380 s1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][3,:,1]] * self.Num_Bin))
381 if method == 'genealized':
382 p0 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][4,:,0]] * self.Num_Bin))
383 p1 = numpy.transpose(numpy.transpose([dataOut.DGauFitParams[iCh][4,:,1]] * self.Num_Bin))
384 elif method == 'squared':
385 p0 = 2.
386 p1 = 2.
387 gau0[iCh] = A0*numpy.exp(-0.5*numpy.abs((x_mtr-v0)/s0)**p0)+N0
388 gau1[iCh] = A1*numpy.exp(-0.5*numpy.abs((x_mtr-v1)/s1)**p1)+N1
389 dataOut.GaussFit0 = gau0
390 dataOut.GaussFit1 = gau1
391
392 print('Leaving ',method ,' double Gaussian fit')
393 return dataOut
372 394
373 395 def FitGau(self, X):
374
375 Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
376
377 SPCparam = []
378 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
379 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
380 SPC_ch1[:] = 0 #numpy.NaN
381 SPC_ch2[:] = 0 #numpy.NaN
382
383
384
396 # print('Entering FitGau')
397 # Assigning the variables
398 Vrange, ch, wnoise, num_intg, SNRlimit = X
399 # Noise Limits
400 noisebl = wnoise * 0.9
401 noisebh = wnoise * 1.1
402 # Radar Velocity
403 Va = max(Vrange)
404 deltav = Vrange[1] - Vrange[0]
405 x = numpy.arange(self.Num_Bin)
406
407 # print ('stop 0')
408
409 # 5 parameters, 2 Gaussians
410 DGauFitParam = numpy.zeros([5, self.Num_Hei,2])
411 DGauFitParam[:] = numpy.NaN
412
413 # SPCparam = []
414 # SPC_ch1 = numpy.zeros([self.Num_Bin,self.Num_Hei])
415 # SPC_ch2 = numpy.zeros([self.Num_Bin,self.Num_Hei])
416 # SPC_ch1[:] = 0 #numpy.NaN
417 # SPC_ch2[:] = 0 #numpy.NaN
418 # print ('stop 1')
385 419 for ht in range(self.Num_Hei):
386
387
420 # print (ht)
421 # print ('stop 2')
422 # Spectra at each range
388 423 spc = numpy.asarray(self.spc)[ch,:,ht]
424 snr = ( spc.mean() - wnoise ) / wnoise
425 snrdB = 10.*numpy.log10(snr)
389 426
427 #print ('stop 3')
428 if snrdB < SNRlimit :
429 # snr = numpy.NaN
430 # SPC_ch1[:,ht] = 0#numpy.NaN
431 # SPC_ch1[:,ht] = 0#numpy.NaN
432 # SPCparam = (SPC_ch1,SPC_ch2)
433 # print ('SNR less than SNRth')
434 continue
435 # wnoise = hildebrand_sekhon(spc,num_intg)
436 # print ('stop 2.01')
390 437 #############################################
391 438 # normalizing spc and noise
392 439 # This part differs from gg1
393 spc_norm_max = max(spc)
440 # spc_norm_max = max(spc) #commented by D. Scipión 19.03.2021
394 441 #spc = spc / spc_norm_max
395 pnoise = pnoise #/ spc_norm_max
442 # pnoise = pnoise #/ spc_norm_max #commented by D. Scipión 19.03.2021
396 443 #############################################
397 444
445 # print ('stop 2.1')
398 446 fatspectra=1.0
399
400 wnoise = noise_ #/ spc_norm_max
447 # noise per channel.... we might want to use the noise at each range
448
449 # wnoise = noise_ #/ spc_norm_max #commented by D. Scipión 19.03.2021
401 450 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
402 451 #if wnoise>1.1*pnoise: # to be tested later
403 452 # wnoise=pnoise
404 noisebl=wnoise*0.9;
405 noisebh=wnoise*1.1
406 spc=spc-wnoise
453 # noisebl = wnoise*0.9
454 # noisebh = wnoise*1.1
455 spc = spc - wnoise # signal
407 456
408 minx=numpy.argmin(spc)
457 # print ('stop 2.2')
458 minx = numpy.argmin(spc)
409 459 #spcs=spc.copy()
410 spcs=numpy.roll(spc,-minx)
411 cum=numpy.cumsum(spcs)
412 tot_noise=wnoise * self.Num_Bin #64;
413
414 snr = sum(spcs)/tot_noise
415 snrdB=10.*numpy.log10(snr)
416
417 if snrdB < SNRlimit :
418 snr = numpy.NaN
419 SPC_ch1[:,ht] = 0#numpy.NaN
420 SPC_ch1[:,ht] = 0#numpy.NaN
421 SPCparam = (SPC_ch1,SPC_ch2)
422 continue
460 spcs = numpy.roll(spc,-minx)
461 cum = numpy.cumsum(spcs)
462 # tot_noise = wnoise * self.Num_Bin #64;
463
464 # print ('stop 2.3')
465 # snr = sum(spcs) / tot_noise
466 # snrdB = 10.*numpy.log10(snr)
467 #print ('stop 3')
468 # if snrdB < SNRlimit :
469 # snr = numpy.NaN
470 # SPC_ch1[:,ht] = 0#numpy.NaN
471 # SPC_ch1[:,ht] = 0#numpy.NaN
472 # SPCparam = (SPC_ch1,SPC_ch2)
473 # print ('SNR less than SNRth')
474 # continue
423 475
424 476
425 477 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
426 478 # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
427
428 cummax=max(cum);
429 epsi=0.08*fatspectra # cumsum to narrow down the energy region
430 cumlo=cummax*epsi;
431 cumhi=cummax*(1-epsi)
432 powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
433
434
479 # print ('stop 4')
480 cummax = max(cum)
481 epsi = 0.08 * fatspectra # cumsum to narrow down the energy region
482 cumlo = cummax * epsi
483 cumhi = cummax * (1-epsi)
484 powerindex = numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
485
486 # print ('stop 5')
435 487 if len(powerindex) < 1:# case for powerindex 0
488 # print ('powerindex < 1')
436 489 continue
437 powerlo=powerindex[0]
438 powerhi=powerindex[-1]
439 powerwidth=powerhi-powerlo
440
441 firstpeak=powerlo+powerwidth/10.# first gaussian energy location
442 secondpeak=powerhi-powerwidth/10.#second gaussian energy location
443 midpeak=(firstpeak+secondpeak)/2.
444 firstamp=spcs[int(firstpeak)]
445 secondamp=spcs[int(secondpeak)]
446 midamp=spcs[int(midpeak)]
490 powerlo = powerindex[0]
491 powerhi = powerindex[-1]
492 powerwidth = powerhi-powerlo
493 if powerwidth <= 1:
494 # print('powerwidth <= 1')
495 continue
496
497 # print ('stop 6')
498 firstpeak = powerlo + powerwidth/10.# first gaussian energy location
499 secondpeak = powerhi - powerwidth/10. #second gaussian energy location
500 midpeak = (firstpeak + secondpeak)/2.
501 firstamp = spcs[int(firstpeak)]
502 secondamp = spcs[int(secondpeak)]
503 midamp = spcs[int(midpeak)]
447 504
448 x=numpy.arange( self.Num_Bin )
449 y_data=spc+wnoise
505 y_data = spc + wnoise
450 506
451 507 ''' single Gaussian '''
452 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
453 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
454 power0=2.
455 amplitude0=midamp
456 state0=[shift0,width0,amplitude0,power0,wnoise]
457 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
458 lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
459
460 chiSq1=lsq1[1];
461
462
508 shift0 = numpy.mod(midpeak+minx, self.Num_Bin )
509 width0 = powerwidth/4.#Initialization entire power of spectrum divided by 4
510 power0 = 2.
511 amplitude0 = midamp
512 state0 = [shift0,width0,amplitude0,power0,wnoise]
513 bnds = ((0,self.Num_Bin-1),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
514 lsq1 = fmin_l_bfgs_b(self.misfit1, state0, args=(y_data,x,num_intg), bounds=bnds, approx_grad=True)
515 # print ('stop 7.1')
516 # print (bnds)
517
518 chiSq1=lsq1[1]
519
520 # print ('stop 8')
463 521 if fatspectra<1.0 and powerwidth<4:
464 522 choice=0
465 523 Amplitude0=lsq1[0][2]
@@ -473,127 +531,142 class GaussianFit(Operation):
473 531 noise=lsq1[0][4]
474 532 #return (numpy.array([shift0,width0,Amplitude0,p0]),
475 533 # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
476
477 ''' two gaussians '''
534
535 # print ('stop 9')
536 ''' two Gaussians '''
478 537 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
479 shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
480 shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
481 width0=powerwidth/6.;
482 width1=width0
483 power0=2.;
484 power1=power0
485 amplitude0=firstamp;
486 amplitude1=secondamp
487 state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
538 shift0 = numpy.mod(firstpeak+minx, self.Num_Bin )
539 shift1 = numpy.mod(secondpeak+minx, self.Num_Bin )
540 width0 = powerwidth/6.
541 width1 = width0
542 power0 = 2.
543 power1 = power0
544 amplitude0 = firstamp
545 amplitude1 = secondamp
546 state0 = [shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
488 547 #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
489 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
548 bnds=((0,self.Num_Bin-1),(1,powerwidth/2.),(0,None),(0.5,3.),(0,self.Num_Bin-1),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
490 549 #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
491 550
551 # print ('stop 10')
492 552 lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )
493 553
554 # print ('stop 11')
555 chiSq2 = lsq2[1]
494 556
495 chiSq2=lsq2[1];
496
497
557 # print ('stop 12')
498 558
499 oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
559 oneG = (chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
500 560
561 # print ('stop 13')
501 562 if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
502 563 if oneG:
503 choice=0
564 choice = 0
504 565 else:
505 w1=lsq2[0][1]; w2=lsq2[0][5]
506 a1=lsq2[0][2]; a2=lsq2[0][6]
507 p1=lsq2[0][3]; p2=lsq2[0][7]
508 s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1;
509 s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
510 gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
566 w1 = lsq2[0][1]; w2 = lsq2[0][5]
567 a1 = lsq2[0][2]; a2 = lsq2[0][6]
568 p1 = lsq2[0][3]; p2 = lsq2[0][7]
569 s1 = (2**(1+1./p1))*scipy.special.gamma(1./p1)/p1
570 s2 = (2**(1+1./p2))*scipy.special.gamma(1./p2)/p2
571 gp1 = a1*w1*s1; gp2 = a2*w2*s2 # power content of each ggaussian with proper p scaling
511 572
512 573 if gp1>gp2:
513 574 if a1>0.7*a2:
514 choice=1
575 choice = 1
515 576 else:
516 choice=2
577 choice = 2
517 578 elif gp2>gp1:
518 579 if a2>0.7*a1:
519 choice=2
580 choice = 2
520 581 else:
521 choice=1
582 choice = 1
522 583 else:
523 choice=numpy.argmax([a1,a2])+1
584 choice = numpy.argmax([a1,a2])+1
524 585 #else:
525 586 #choice=argmin([std2a,std2b])+1
526 587
527 588 else: # with low SNR go to the most energetic peak
528 choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
529
530
531 shift0=lsq2[0][0];
532 vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
533 shift1=lsq2[0][4];
534 vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
535
536 max_vel = 1.0
537
589 choice = numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
590
591 # print ('stop 14')
592 shift0 = lsq2[0][0]
593 vel0 = Vrange[0] + shift0 * deltav
594 shift1 = lsq2[0][4]
595 # vel1=Vrange[0] + shift1 * deltav
596
597 # max_vel = 1.0
598 # Va = max(Vrange)
599 # deltav = Vrange[1]-Vrange[0]
600 # print ('stop 15')
538 601 #first peak will be 0, second peak will be 1
539 if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range
540 shift0=lsq2[0][0]
541 width0=lsq2[0][1]
542 Amplitude0=lsq2[0][2]
543 p0=lsq2[0][3]
544
545 shift1=lsq2[0][4]
546 width1=lsq2[0][5]
547 Amplitude1=lsq2[0][6]
548 p1=lsq2[0][7]
549 noise=lsq2[0][8]
602 # if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range # Commented by D.Scipión 19.03.2021
603 if vel0 > -Va and vel0 < Va : #first peak is in the correct range
604 shift0 = lsq2[0][0]
605 width0 = lsq2[0][1]
606 Amplitude0 = lsq2[0][2]
607 p0 = lsq2[0][3]
608
609 shift1 = lsq2[0][4]
610 width1 = lsq2[0][5]
611 Amplitude1 = lsq2[0][6]
612 p1 = lsq2[0][7]
613 noise = lsq2[0][8]
550 614 else:
551 shift1=lsq2[0][0]
552 width1=lsq2[0][1]
553 Amplitude1=lsq2[0][2]
554 p1=lsq2[0][3]
615 shift1 = lsq2[0][0]
616 width1 = lsq2[0][1]
617 Amplitude1 = lsq2[0][2]
618 p1 = lsq2[0][3]
555 619
556 shift0=lsq2[0][4]
557 width0=lsq2[0][5]
558 Amplitude0=lsq2[0][6]
559 p0=lsq2[0][7]
560 noise=lsq2[0][8]
620 shift0 = lsq2[0][4]
621 width0 = lsq2[0][5]
622 Amplitude0 = lsq2[0][6]
623 p0 = lsq2[0][7]
624 noise = lsq2[0][8]
561 625
562 626 if Amplitude0<0.05: # in case the peak is noise
563 shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN]
627 shift0,width0,Amplitude0,p0 = 4*[numpy.NaN]
564 628 if Amplitude1<0.05:
565 shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN]
566
567
568 SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
569 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
570 SPCparam = (SPC_ch1,SPC_ch2)
571
572
573 return GauSPC
629 shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
630
631 # print ('stop 16 ')
632 # SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0)/width0)**p0)
633 # SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1)/width1)**p1)
634 # SPCparam = (SPC_ch1,SPC_ch2)
635
636 DGauFitParam[0,ht,0] = noise
637 DGauFitParam[0,ht,1] = noise
638 DGauFitParam[1,ht,0] = Amplitude0
639 DGauFitParam[1,ht,1] = Amplitude1
640 DGauFitParam[2,ht,0] = Vrange[0] + shift0 * deltav
641 DGauFitParam[2,ht,1] = Vrange[0] + shift1 * deltav
642 DGauFitParam[3,ht,0] = width0 * deltav
643 DGauFitParam[3,ht,1] = width1 * deltav
644 DGauFitParam[4,ht,0] = p0
645 DGauFitParam[4,ht,1] = p1
646
647 # print (DGauFitParam.shape)
648 # print ('Leaving FitGau')
649 return DGauFitParam
650 # return SPCparam
651 # return GauSPC
574 652
575 653 def y_model1(self,x,state):
576 shift0,width0,amplitude0,power0,noise=state
577 model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
578
579 model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
580
581 model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
582 return model0+model0u+model0d+noise
654 shift0, width0, amplitude0, power0, noise = state
655 model0 = amplitude0*numpy.exp(-0.5*abs((x - shift0)/width0)**power0)
656 model0u = amplitude0*numpy.exp(-0.5*abs((x - shift0 - self.Num_Bin)/width0)**power0)
657 model0d = amplitude0*numpy.exp(-0.5*abs((x - shift0 + self.Num_Bin)/width0)**power0)
658 return model0 + model0u + model0d + noise
583 659
584 660 def y_model2(self,x,state): #Equation for two generalized Gaussians with Nyquist
585 shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,noise=state
586 model0=amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
587
588 model0u=amplitude0*numpy.exp(-0.5*abs((x-shift0- self.Num_Bin )/width0)**power0)
589
590 model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
591 model1=amplitude1*numpy.exp(-0.5*abs((x-shift1)/width1)**power1)
592
593 model1u=amplitude1*numpy.exp(-0.5*abs((x-shift1- self.Num_Bin )/width1)**power1)
594
595 model1d=amplitude1*numpy.exp(-0.5*abs((x-shift1+ self.Num_Bin )/width1)**power1)
596 return model0+model0u+model0d+model1+model1u+model1d+noise
661 shift0, width0, amplitude0, power0, shift1, width1, amplitude1, power1, noise = state
662 model0 = amplitude0*numpy.exp(-0.5*abs((x-shift0)/width0)**power0)
663 model0u = amplitude0*numpy.exp(-0.5*abs((x - shift0 - self.Num_Bin)/width0)**power0)
664 model0d = amplitude0*numpy.exp(-0.5*abs((x - shift0 + self.Num_Bin)/width0)**power0)
665
666 model1 = amplitude1*numpy.exp(-0.5*abs((x - shift1)/width1)**power1)
667 model1u = amplitude1*numpy.exp(-0.5*abs((x - shift1 - self.Num_Bin)/width1)**power1)
668 model1d = amplitude1*numpy.exp(-0.5*abs((x - shift1 + self.Num_Bin)/width1)**power1)
669 return model0 + model0u + model0d + model1 + model1u + model1d + noise
597 670
598 671 def misfit1(self,state,y_data,x,num_intg): # This function compares how close real data is with the model data, the close it is, the better it is.
599 672
@@ -625,7 +698,9 class PrecipitationProc(Operation):
625 698 self.i=0
626 699
627 700 def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118,
628 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km2 = 0.93, Altitude=3350):
701 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km2 = 0.93, Altitude=3350,SNRdBlimit=-30):
702
703 # print ('Entering PrecepitationProc ... ')
629 704
630 705 if radar == "MIRA35C" :
631 706
@@ -673,7 +748,6 class PrecipitationProc(Operation):
673 748 SignalPower[numpy.where(SignalPower < 0)] = 1e-20
674 749
675 750 SPCmean = numpy.mean(SignalPower, 0)
676
677 751 Pr = SPCmean[:,:]/dataOut.normFactor
678 752
679 753 # Declaring auxiliary variables
@@ -708,9 +782,9 class PrecipitationProc(Operation):
708 782
709 783 # Censoring the data
710 784 # Removing data with SNRth < 0dB se debe considerar el SNR por canal
711 SNRth = 10**(-30/10) #-20dB
712 novalid = numpy.where((dataOut.data_SNR[0,:] <SNRth) | (dataOut.data_SNR[1,:] <SNRth) | (dataOut.data_SNR[2,:] <SNRth)) # AND condition. Maybe OR condition better
713 W = numpy.nanmean(dataOut.data_DOP,0)
785 SNRth = 10**(SNRdBlimit/10) #-30dB
786 novalid = numpy.where((dataOut.data_snr[0,:] <SNRth) | (dataOut.data_snr[1,:] <SNRth) | (dataOut.data_snr[2,:] <SNRth)) # AND condition. Maybe OR condition better
787 W = numpy.nanmean(dataOut.data_dop,0)
714 788 W[novalid] = numpy.NaN
715 789 Ze_org[novalid] = numpy.NaN
716 790 RR[novalid] = numpy.NaN
@@ -720,8 +794,10 class PrecipitationProc(Operation):
720 794 dataOut.channelList = [0,1,2]
721 795
722 796 dataOut.data_param[0]=10*numpy.log10(Ze_org)
723 dataOut.data_param[1]=W
797 dataOut.data_param[1]=-W
724 798 dataOut.data_param[2]=RR
799
800 # print ('Leaving PrecepitationProc ... ')
725 801 return dataOut
726 802
727 803 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
@@ -792,23 +868,11 class FullSpectralAnalysis(Operation):
792 868 Parameters affected: Winds, height range, SNR
793 869
794 870 """
795 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRlimit=7, minheight=None, maxheight=None):
796
797 self.indice=int(numpy.random.rand()*1000)
871 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRdBlimit=-30,
872 minheight=None, maxheight=None, NegativeLimit=None, PositiveLimit=None):
798 873
799 874 spc = dataOut.data_pre[0].copy()
800 875 cspc = dataOut.data_pre[1]
801
802 """Erick: NOTE THE RANGE OF THE PULSE TX MUST BE REMOVED"""
803
804 SNRspc = spc.copy()
805 SNRspc[:,:,0:7]= numpy.NaN # D. Scipión... the cleaning should not be hardwired in the code... it needs to be flexible... NEEDS TO BE REMOVED
806
807 """##########################################"""
808
809
810 nChannel = spc.shape[0]
811 nProfiles = spc.shape[1]
812 876 nHeights = spc.shape[2]
813 877
814 878 # first_height = 0.75 #km (ref: data header 20170822)
@@ -835,112 +899,81 class FullSpectralAnalysis(Operation):
835 899 else:
836 900 ChanDist = numpy.array([[Xi01, Eta01],[Xi02,Eta02],[Xi12,Eta12]])
837 901
838 FrecRange = dataOut.spc_range[0]
839
840 data_SNR=numpy.zeros([nProfiles])
841 noise = dataOut.noise
842
843 dataOut.data_snr = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0]
844
845 dataOut.data_snr[numpy.where( dataOut.data_snr <0 )] = 1e-20
846
847
848 data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN
849
850 velocityX=[]
851 velocityY=[]
852 velocityV=[]
902 # 4 variables: zonal, meridional, vertical, and average SNR
903 data_param = numpy.zeros([4,nHeights]) * numpy.NaN
904 velocityX = numpy.zeros([nHeights]) * numpy.NaN
905 velocityY = numpy.zeros([nHeights]) * numpy.NaN
906 velocityZ = numpy.zeros([nHeights]) * numpy.NaN
853 907
854 dbSNR = 10*numpy.log10(dataOut.data_snr)
855 dbSNR = numpy.average(dbSNR,0)
908 dbSNR = 10*numpy.log10(numpy.average(dataOut.data_snr,0))
856 909
857 910 '''***********************************************WIND ESTIMATION**************************************'''
858
859 911 for Height in range(nHeights):
860 912
861 913 if Height >= range_min and Height < range_max:
862 # error_code unused, yet maybe useful for future analysis.
863 [Vzon,Vmer,Vver, error_code] = self.WindEstimation(spc[:,:,Height], cspc[:,:,Height], pairsList, ChanDist, Height, noise, dataOut.spc_range, dbSNR[Height], SNRlimit)
864 else:
865 Vzon,Vmer,Vver = 0., 0., numpy.NaN
866
867
868 if abs(Vzon) < 100. and abs(Vzon) > 0. and abs(Vmer) < 100. and abs(Vmer) > 0.:
869 velocityX=numpy.append(velocityX, Vzon)
870 velocityY=numpy.append(velocityY, -Vmer)
871
872 else:
873 velocityX=numpy.append(velocityX, numpy.NaN)
874 velocityY=numpy.append(velocityY, numpy.NaN)
875
876 if dbSNR[Height] > SNRlimit:
877 velocityV=numpy.append(velocityV, -Vver) # reason for this minus sign -> convention? (taken from Ericks version) D.S. yes!
878 else:
879 velocityV=numpy.append(velocityV, numpy.NaN)
880
881
882 '''Change the numpy.array (velocityX) sign when trying to process BLTR data (Erick)'''
883 data_output[0] = numpy.array(velocityX)
884 data_output[1] = numpy.array(velocityY)
885 data_output[2] = velocityV
886
887
888 dataOut.data_output = data_output
889
914 # error_code will be useful in future analysis
915 [Vzon,Vmer,Vver, error_code] = self.WindEstimation(spc[:,:,Height], cspc[:,:,Height], pairsList,
916 ChanDist, Height, dataOut.noise, dataOut.spc_range, dbSNR[Height], SNRdBlimit, NegativeLimit, PositiveLimit,dataOut.frequency)
917
918 if abs(Vzon) < 100. and abs(Vmer) < 100.:
919 velocityX[Height] = Vzon
920 velocityY[Height] = -Vmer
921 velocityZ[Height] = Vver
922
923 # Censoring data with SNR threshold
924 dbSNR [dbSNR < SNRdBlimit] = numpy.NaN
925
926 data_param[0] = velocityX
927 data_param[1] = velocityY
928 data_param[2] = velocityZ
929 data_param[3] = dbSNR
930 dataOut.data_param = data_param
890 931 return dataOut
891 932
892
893 933 def moving_average(self,x, N=2):
894 934 """ convolution for smoothenig data. note that last N-1 values are convolution with zeroes """
895 935 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
896 936
897 937 def gaus(self,xSamples,Amp,Mu,Sigma):
898 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
938 return Amp * numpy.exp(-0.5*((xSamples - Mu)/Sigma)**2)
899 939
900 940 def Moments(self, ySamples, xSamples):
901 941 Power = numpy.nanmean(ySamples) # Power, 0th Moment
902 942 yNorm = ySamples / numpy.nansum(ySamples)
903 943 RadVel = numpy.nansum(xSamples * yNorm) # Radial Velocity, 1st Moment
904 Sigma2 = abs(numpy.nansum(yNorm * (xSamples - RadVel)**2)) # Spectral Width, 2nd Moment
905 StdDev = Sigma2**0.5 # Desv. Estandar, Ancho espectral
944 Sigma2 = numpy.nansum(yNorm * (xSamples - RadVel)**2) # Spectral Width, 2nd Moment
945 StdDev = numpy.sqrt(numpy.abs(Sigma2)) # Desv. Estandar, Ancho espectral
906 946 return numpy.array([Power,RadVel,StdDev])
907 947
908 948 def StopWindEstimation(self, error_code):
909 '''
910 the wind calculation and returns zeros
911 '''
912 Vzon = 0
913 Vmer = 0
914 Vver = numpy.nan
949 Vzon = numpy.NaN
950 Vmer = numpy.NaN
951 Vver = numpy.NaN
915 952 return Vzon, Vmer, Vver, error_code
916 953
917 954 def AntiAliasing(self, interval, maxstep):
918 955 """
919 956 function to prevent errors from aliased values when computing phaseslope
920 957 """
921 antialiased = numpy.zeros(len(interval))*0.0
958 antialiased = numpy.zeros(len(interval))
922 959 copyinterval = interval.copy()
923 960
924 961 antialiased[0] = copyinterval[0]
925 962
926 963 for i in range(1,len(antialiased)):
927
928 964 step = interval[i] - interval[i-1]
929
930 965 if step > maxstep:
931 966 copyinterval -= 2*numpy.pi
932 967 antialiased[i] = copyinterval[i]
933
934 968 elif step < maxstep*(-1):
935 969 copyinterval += 2*numpy.pi
936 970 antialiased[i] = copyinterval[i]
937
938 971 else:
939 972 antialiased[i] = copyinterval[i].copy()
940 973
941 974 return antialiased
942 975
943 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
976 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit, NegativeLimit, PositiveLimit, radfreq):
944 977 """
945 978 Function that Calculates Zonal, Meridional and Vertical wind velocities.
946 979 Initial Version by E. Bocanegra updated by J. Zibell until Nov. 2019.
@@ -972,46 +1005,40 class FullSpectralAnalysis(Operation):
972 1005
973 1006 error_code = 0
974 1007
975
976 SPC_Samples = numpy.ones([spc.shape[0],spc.shape[1]]) # for normalized spc values for one height
977 phase = numpy.ones([spc.shape[0],spc.shape[1]]) # phase between channels
978 CSPC_Samples = numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_) # for normalized cspc values
979 PhaseSlope = numpy.zeros(spc.shape[0]) # slope of the phases, channelwise
980 PhaseInter = numpy.ones(spc.shape[0]) # intercept to the slope of the phases, channelwise
981 xFrec = AbbsisaRange[0][0:spc.shape[1]] # frequency range
982 xVel = AbbsisaRange[2][0:spc.shape[1]] # velocity range
983 SPCav = numpy.average(spc, axis=0)-numpy.average(noise) # spc[0]-noise[0] %D.S. why??? I suggest only spc....
984
985 SPCmoments_vel = self.Moments(SPCav, xVel ) # SPCmoments_vel[1] corresponds to vertical velocity and is used to determine if signal corresponds to wind (if .. <3)
986 # D.S. I suggest to each moment to be calculated independently, because the signal level y/o interferences are not the same in all channels and
987 # signal or SNR seems to be contaminated
988 CSPCmoments = []
1008 nChan = spc.shape[0]
1009 nProf = spc.shape[1]
1010 nPair = cspc.shape[0]
1011
1012 SPC_Samples = numpy.zeros([nChan, nProf]) # for normalized spc values for one height
1013 CSPC_Samples = numpy.zeros([nPair, nProf], dtype=numpy.complex_) # for normalized cspc values
1014 phase = numpy.zeros([nPair, nProf]) # phase between channels
1015 PhaseSlope = numpy.zeros(nPair) # slope of the phases, channelwise
1016 PhaseInter = numpy.zeros(nPair) # intercept to the slope of the phases, channelwise
1017 xFrec = AbbsisaRange[0][:-1] # frequency range
1018 xVel = AbbsisaRange[2][:-1] # velocity range
1019 xSamples = xFrec # the frequency range is taken
1020 delta_x = xSamples[1] - xSamples[0] # delta_f or delta_x
1021
1022 # only consider velocities with in NegativeLimit and PositiveLimit
1023 if (NegativeLimit is None):
1024 NegativeLimit = numpy.min(xVel)
1025 if (PositiveLimit is None):
1026 PositiveLimit = numpy.max(xVel)
1027 xvalid = numpy.where((xVel > NegativeLimit) & (xVel < PositiveLimit))
1028 xSamples_zoom = xSamples[xvalid]
989 1029
990 1030 '''Getting Eij and Nij'''
991
992 1031 Xi01, Xi02, Xi12 = ChanDist[:,0]
993 1032 Eta01, Eta02, Eta12 = ChanDist[:,1]
994 1033
995 # update nov 19
996 widthlimit = 7 # maximum width in Hz of the gaussian, empirically determined. Anything above 10 is unrealistic, often values between 1 and 5 correspond to proper fits.
997
1034 # spwd limit - updated by D. Scipión 30.03.2021
1035 widthlimit = 10
998 1036 '''************************* SPC is normalized ********************************'''
999
1000 spc_norm = spc.copy() # need copy() because untouched spc is needed for normalization of cspc below
1001 spc_norm = numpy.where(numpy.isfinite(spc_norm), spc_norm, numpy.NAN)
1002
1003 # D. Scipión: It is necessary to define DeltaF or DeltaV... it is wrong to use Factor_Norm. It's constant... not a variable
1004
1037 spc_norm = spc.copy()
1005 1038 # For each channel
1006 for i in range(spc.shape[0]):
1007
1008 spc_sub = spc_norm[i,:] - noise[i] # spc not smoothed here or in previous version.
1009 # D. Scipión: Factor_Norm has to be replaced by DeltaF or DeltaV - It's a constant
1010 Factor_Norm = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc_sub)) # usually = Freq range / nfft
1011 normalized_spc = spc_sub / (numpy.nansum(numpy.abs(spc_sub)) * Factor_Norm)
1012
1013 xSamples = xFrec # the frequency range is taken
1014 SPC_Samples[i] = normalized_spc # Normalized SPC values are taken
1039 for i in range(nChan):
1040 spc_sub = spc_norm[i,:] - noise[i] # only the signal power
1041 SPC_Samples[i] = spc_sub / (numpy.nansum(spc_sub) * delta_x)
1015 1042
1016 1043 '''********************** FITTING MEAN SPC GAUSSIAN **********************'''
1017 1044
@@ -1025,30 +1052,26 class FullSpectralAnalysis(Operation):
1025 1052 due to subtraction of the noise, some values are below zero. Raw "spc" values should be
1026 1053 >= 0, as it is the modulus squared of the signals (complex * it's conjugate)
1027 1054 """
1028
1029 SPCMean = numpy.average(SPC_Samples, axis=0)
1030
1031 popt = [1e-10,0,1e-10]
1032 SPCMoments = self.Moments(SPCMean, xSamples)
1033
1034 if dbSNR > SNRlimit and numpy.abs(SPCmoments_vel[1]) < 3:
1055 # initial conditions
1056 popt = [1e-10,0,1e-10]
1057 # Spectra average
1058 SPCMean = numpy.average(SPC_Samples,0)
1059 # Moments in frequency
1060 SPCMoments = self.Moments(SPCMean[xvalid], xSamples_zoom)
1061
1062 # Gauss Fit SPC in frequency domain
1063 if dbSNR > SNRlimit: # only if SNR > SNRth
1035 1064 try:
1036 popt,pcov = curve_fit(self.gaus,xSamples,SPCMean,p0=SPCMoments)#, bounds=(-numpy.inf, [numpy.inf, numpy.inf, 10])). Setting bounds does not make the code faster but only keeps the fit from finding the minimum.
1037 if popt[2] > widthlimit: # CONDITION
1065 popt,pcov = curve_fit(self.gaus,xSamples_zoom,SPCMean[xvalid],p0=SPCMoments)
1066 if popt[2] <= 0 or popt[2] > widthlimit: # CONDITION
1038 1067 return self.StopWindEstimation(error_code = 1)
1039
1040 FitGauss = self.gaus(xSamples,*popt)
1041
1068 FitGauss = self.gaus(xSamples_zoom,*popt)
1042 1069 except :#RuntimeError:
1043 1070 return self.StopWindEstimation(error_code = 2)
1044
1045 1071 else:
1046 1072 return self.StopWindEstimation(error_code = 3)
1047 1073
1048
1049
1050 1074 '''***************************** CSPC Normalization *************************
1051 new section:
1052 1075 The Spc spectra are used to normalize the crossspectra. Peaks from precipitation
1053 1076 influence the norm which is not desired. First, a range is identified where the
1054 1077 wind peak is estimated -> sum_wind is sum of those frequencies. Next, the area
@@ -1060,161 +1083,82 class FullSpectralAnalysis(Operation):
1060 1083
1061 1084 A norm is found according to Briggs 92.
1062 1085 '''
1063
1064 radarWavelength = 0.6741 # meters
1065 # D.S. This does not need to hardwired... It has to be in function of the radar frequency
1066
1067 count_limit_freq = numpy.abs(popt[1]) + widthlimit # Hz, m/s can be also used if velocity is desired abscissa.
1068 # count_limit_freq = numpy.max(xFrec)
1069
1070 channel_integrals = numpy.zeros(3)
1071
1072 for i in range(spc.shape[0]):
1073 '''
1074 find the point in array corresponding to count_limit frequency.
1075 sum over all frequencies in the range around zero Hz @ math.ceil(N_freq/2)
1076 '''
1077 N_freq = numpy.count_nonzero(~numpy.isnan(spc[i,:]))
1078 count_limit_int = int(math.ceil( count_limit_freq / numpy.max(xFrec) * (N_freq / 2) )) # gives integer point
1079 sum_wind = numpy.nansum( spc[i, (math.ceil(N_freq/2) - count_limit_int) : (math.ceil(N_freq / 2) + count_limit_int)] ) #N_freq/2 is where frequency (velocity) is zero, i.e. middle of spectrum.
1080 sum_noise = (numpy.mean(spc[i, :4]) + numpy.mean(spc[i, -6:-2]))/2.0 * (N_freq - 2*count_limit_int)
1081 channel_integrals[i] = (sum_noise + sum_wind) * (2*numpy.max(xFrec) / N_freq)
1082
1083
1084 cross_integrals_peak = numpy.zeros(3)
1085 # cross_integrals_totalrange = numpy.zeros(3)
1086
1087 for i in range(spc.shape[0]):
1088
1089 cspc_norm = cspc[i,:].copy() # cspc not smoothed here or in previous version
1090
1086 # for each pair
1087 for i in range(nPair):
1088 cspc_norm = cspc[i,:].copy()
1091 1089 chan_index0 = pairsList[i][0]
1092 1090 chan_index1 = pairsList[i][1]
1093
1094 cross_integrals_peak[i] = channel_integrals[chan_index0]*channel_integrals[chan_index1]
1095 normalized_cspc = cspc_norm / numpy.sqrt(cross_integrals_peak[i])
1096 CSPC_Samples[i] = normalized_cspc
1097
1098 ''' Finding cross integrals without subtracting any peaks:'''
1099 # FactorNorm0 = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc[chan_index0,:]))
1100 # FactorNorm1 = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc[chan_index1,:]))
1101 # cross_integrals_totalrange[i] = (numpy.nansum(spc[chan_index0,:])) * FactorNorm0 * (numpy.nansum(spc[chan_index1,:])) * FactorNorm1
1102 # normalized_cspc = cspc_norm / numpy.sqrt(cross_integrals_totalrange[i])
1103 # CSPC_Samples[i] = normalized_cspc
1104
1105
1091 CSPC_Samples[i] = cspc_norm / (numpy.sqrt(numpy.nansum(spc_norm[chan_index0])*numpy.nansum(spc_norm[chan_index1])) * delta_x)
1106 1092 phase[i] = numpy.arctan2(CSPC_Samples[i].imag, CSPC_Samples[i].real)
1107 1093
1094 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPC_Samples[0,xvalid]), xSamples_zoom),
1095 self.Moments(numpy.abs(CSPC_Samples[1,xvalid]), xSamples_zoom),
1096 self.Moments(numpy.abs(CSPC_Samples[2,xvalid]), xSamples_zoom)])
1108 1097
1109 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPC_Samples[0]), xSamples),
1110 self.Moments(numpy.abs(CSPC_Samples[1]), xSamples),
1111 self.Moments(numpy.abs(CSPC_Samples[2]), xSamples)])
1112
1113
1114 '''***Sorting out NaN entries***'''
1115 CSPCMask01 = numpy.abs(CSPC_Samples[0])
1116 CSPCMask02 = numpy.abs(CSPC_Samples[1])
1117 CSPCMask12 = numpy.abs(CSPC_Samples[2])
1118
1119 mask01 = ~numpy.isnan(CSPCMask01)
1120 mask02 = ~numpy.isnan(CSPCMask02)
1121 mask12 = ~numpy.isnan(CSPCMask12)
1122
1123 CSPCMask01 = CSPCMask01[mask01]
1124 CSPCMask02 = CSPCMask02[mask02]
1125 CSPCMask12 = CSPCMask12[mask12]
1126
1127
1128 popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
1129 FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0
1098 popt01, popt02, popt12 = [1e-10,0,1e-10], [1e-10,0,1e-10] ,[1e-10,0,1e-10]
1099 FitGauss01, FitGauss02, FitGauss12 = numpy.zeros(len(xSamples)), numpy.zeros(len(xSamples)), numpy.zeros(len(xSamples))
1130 1100
1131 1101 '''*******************************FIT GAUSS CSPC************************************'''
1132
1133 1102 try:
1134 popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
1103 popt01,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[0][xvalid]),p0=CSPCmoments[0])
1135 1104 if popt01[2] > widthlimit: # CONDITION
1136 1105 return self.StopWindEstimation(error_code = 4)
1137
1138 popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1])
1106 popt02,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[1][xvalid]),p0=CSPCmoments[1])
1139 1107 if popt02[2] > widthlimit: # CONDITION
1140 1108 return self.StopWindEstimation(error_code = 4)
1141
1142 popt12,pcov = curve_fit(self.gaus,xSamples[mask12],numpy.abs(CSPCMask12),p0=CSPCmoments[2])
1109 popt12,pcov = curve_fit(self.gaus,xSamples_zoom,numpy.abs(CSPC_Samples[2][xvalid]),p0=CSPCmoments[2])
1143 1110 if popt12[2] > widthlimit: # CONDITION
1144 1111 return self.StopWindEstimation(error_code = 4)
1145 1112
1146 FitGauss01 = self.gaus(xSamples, *popt01)
1147 FitGauss02 = self.gaus(xSamples, *popt02)
1148 FitGauss12 = self.gaus(xSamples, *popt12)
1149
1113 FitGauss01 = self.gaus(xSamples_zoom, *popt01)
1114 FitGauss02 = self.gaus(xSamples_zoom, *popt02)
1115 FitGauss12 = self.gaus(xSamples_zoom, *popt12)
1150 1116 except:
1151 1117 return self.StopWindEstimation(error_code = 5)
1152 1118
1153 1119
1154 1120 '''************* Getting Fij ***************'''
1155
1156
1157 #Punto en Eje X de la Gaussiana donde se encuentra el centro -- x-axis point of the gaussian where the center is located
1158 # -> PointGauCenter
1121 # x-axis point of the gaussian where the center is located from GaussFit of spectra
1159 1122 GaussCenter = popt[1]
1160 ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()]
1161 PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0]
1123 ClosestCenter = xSamples_zoom[numpy.abs(xSamples_zoom-GaussCenter).argmin()]
1124 PointGauCenter = numpy.where(xSamples_zoom==ClosestCenter)[0][0]
1162 1125
1163 #Punto e^-1 hubicado en la Gaussiana -- point where e^-1 is located in the gaussian
1126 # Point where e^-1 is located in the gaussian
1164 1127 PeMinus1 = numpy.max(FitGauss) * numpy.exp(-1)
1165 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss"
1128 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # The closest point to"Peminus1" in "FitGauss"
1166 1129 PointFij = numpy.where(FitGauss==FijClosest)[0][0]
1167
1168 Fij = numpy.abs(xSamples[PointFij] - xSamples[PointGauCenter])
1130 Fij = numpy.abs(xSamples_zoom[PointFij] - xSamples_zoom[PointGauCenter])
1169 1131
1170 1132 '''********** Taking frequency ranges from mean SPCs **********'''
1171
1172 #GaussCenter = popt[1] #Primer momento 01
1173 GauWidth = popt[2] * 3/2 #Ancho de banda de Gau01 -- Bandwidth of Gau01 TODO why *3/2?
1133 GauWidth = popt[2] * 3/2 # Bandwidth of Gau01
1174 1134 Range = numpy.empty(2)
1175 1135 Range[0] = GaussCenter - GauWidth
1176 1136 Range[1] = GaussCenter + GauWidth
1177 #Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max) -- Point in x-axis where the bandwidth is located (min:max)
1178 ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()]
1179 ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()]
1180
1181 PointRangeMin = numpy.where(xSamples==ClosRangeMin)[0][0]
1182 PointRangeMax = numpy.where(xSamples==ClosRangeMax)[0][0]
1183
1137 # Point in x-axis where the bandwidth is located (min:max)
1138 ClosRangeMin = xSamples_zoom[numpy.abs(xSamples_zoom-Range[0]).argmin()]
1139 ClosRangeMax = xSamples_zoom[numpy.abs(xSamples_zoom-Range[1]).argmin()]
1140 PointRangeMin = numpy.where(xSamples_zoom==ClosRangeMin)[0][0]
1141 PointRangeMax = numpy.where(xSamples_zoom==ClosRangeMax)[0][0]
1184 1142 Range = numpy.array([ PointRangeMin, PointRangeMax ])
1185
1186 FrecRange = xFrec[ Range[0] : Range[1] ]
1187
1143 FrecRange = xSamples_zoom[ Range[0] : Range[1] ]
1188 1144
1189 1145 '''************************** Getting Phase Slope ***************************'''
1190
1191 for i in range(1,3): # Changed to only compute two
1192
1193 if len(FrecRange) > 5 and len(FrecRange) < spc.shape[1] * 0.3:
1194 # PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=1) #used before to smooth phase with N=3
1195 PhaseRange = phase[i,Range[0]:Range[1]].copy()
1196
1146 for i in range(nPair):
1147 if len(FrecRange) > 5:
1148 PhaseRange = phase[i, xvalid[0][Range[0]:Range[1]]].copy()
1197 1149 mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
1198
1199
1200 1150 if len(FrecRange) == len(PhaseRange):
1201
1202 1151 try:
1203 1152 slope, intercept, _, _, _ = stats.linregress(FrecRange[mask], self.AntiAliasing(PhaseRange[mask], 4.5))
1204 1153 PhaseSlope[i] = slope
1205 1154 PhaseInter[i] = intercept
1206
1207 1155 except:
1208 1156 return self.StopWindEstimation(error_code = 6)
1209
1210 1157 else:
1211 1158 return self.StopWindEstimation(error_code = 7)
1212
1213 1159 else:
1214 1160 return self.StopWindEstimation(error_code = 8)
1215 1161
1216
1217
1218 1162 '''*** Constants A-H correspond to the convention as in Briggs and Vincent 1992 ***'''
1219 1163
1220 1164 '''Getting constant C'''
@@ -1222,9 +1166,12 class FullSpectralAnalysis(Operation):
1222 1166
1223 1167 '''****** Getting constants F and G ******'''
1224 1168 MijEijNij = numpy.array([[Xi02,Eta02], [Xi12,Eta12]])
1225 MijResult0 = (-PhaseSlope[1] * cC) / (2*numpy.pi)
1226 MijResult1 = (-PhaseSlope[2] * cC) / (2*numpy.pi)
1227 MijResults = numpy.array([MijResult0,MijResult1])
1169 # MijEijNij = numpy.array([[Xi01,Eta01], [Xi02,Eta02], [Xi12,Eta12]])
1170 # MijResult0 = (-PhaseSlope[0] * cC) / (2*numpy.pi)
1171 MijResult1 = (-PhaseSlope[1] * cC) / (2*numpy.pi)
1172 MijResult2 = (-PhaseSlope[2] * cC) / (2*numpy.pi)
1173 # MijResults = numpy.array([MijResult0, MijResult1, MijResult2])
1174 MijResults = numpy.array([MijResult1, MijResult2])
1228 1175 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1229 1176
1230 1177 '''****** Getting constants A, B and H ******'''
@@ -1232,39 +1179,22 class FullSpectralAnalysis(Operation):
1232 1179 W02 = numpy.nanmax( FitGauss02 )
1233 1180 W12 = numpy.nanmax( FitGauss12 )
1234 1181
1235 WijResult0 = ((cF * Xi01 + cG * Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi / cC))
1236 WijResult1 = ((cF * Xi02 + cG * Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi / cC))
1237 WijResult2 = ((cF * Xi12 + cG * Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi / cC))
1238
1239 WijResults = numpy.array([WijResult0, WijResult1, WijResult2])
1182 WijResult01 = ((cF * Xi01 + cG * Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi / cC))
1183 WijResult02 = ((cF * Xi02 + cG * Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi / cC))
1184 WijResult12 = ((cF * Xi12 + cG * Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi / cC))
1185 WijResults = numpy.array([WijResult01, WijResult02, WijResult12])
1240 1186
1241 1187 WijEijNij = numpy.array([ [Xi01**2, Eta01**2, 2*Xi01*Eta01] , [Xi02**2, Eta02**2, 2*Xi02*Eta02] , [Xi12**2, Eta12**2, 2*Xi12*Eta12] ])
1242 1188 (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1243 1189
1244 1190 VxVy = numpy.array([[cA,cH],[cH,cB]])
1245 1191 VxVyResults = numpy.array([-cF,-cG])
1246 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1247
1248 Vzon = Vy
1249 Vmer = Vx
1250
1251 # Vmag=numpy.sqrt(Vzon**2+Vmer**2) # unused
1252 # Vang=numpy.arctan2(Vmer,Vzon) # unused
1253
1254
1255 ''' using frequency as abscissa. Due to three channels, the offzenith angle is zero
1256 and Vrad equal to Vver. formula taken from Briggs 92, figure 4.
1257 '''
1258 if numpy.abs( popt[1] ) < 3.5 and len(FrecRange) > 4:
1259 Vver = 0.5 * radarWavelength * popt[1] * 100 # *100 to get cm (/s)
1260 else:
1261 Vver = numpy.NaN
1262
1192 (Vmer,Vzon) = numpy.linalg.solve(VxVy, VxVyResults)
1193 Vver = -SPCMoments[1]*SPEED_OF_LIGHT/(2*radfreq)
1263 1194 error_code = 0
1264 1195
1265 1196 return Vzon, Vmer, Vver, error_code
1266 1197
1267
1268 1198 class SpectralMoments(Operation):
1269 1199
1270 1200 '''
General Comments 0
You need to be logged in to leave comments. Login now