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