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