##// 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, (1164 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
227
221 self.ich = 0
228 def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5):
222 self.ir = 0
229
223
230
224 def run(self, dataOut, ClutterWidth=2.5):
231 #Limite de vientos
225 # print ('Entering RemoveWideGC ... ')
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))
254
255 valleyindex = j2index[numpy.where(junk4>1)]
256 peakindex = j1index[numpy.where(junk3>1)]
257
258 isvalid = numpy.squeeze(numpy.where(numpy.abs(VelRange[gc_values[peakindex]]) <= 2.5*dv))
259 if numpy.size(isvalid) == 0 :
260 continue
261 if numpy.size(isvalid) >1 :
262 vindex = numpy.argmax(self.spc[ich,gc_values[peakindex[isvalid]],ir])
263 isvalid = isvalid[vindex]
264
265 # clutter peak
266 gcpeak = peakindex[isvalid]
267 vl = numpy.where(valleyindex < gcpeak)
268 if numpy.size(vl) == 0:
269 continue
270 gcvl = valleyindex[vl[0][-1]]
271 vr = numpy.where(valleyindex > gcpeak)
272 if numpy.size(vr) == 0:
273 continue
274 gcvr = valleyindex[vr[0][0]]
275
276 # Removing the clutter
277 interpindex = numpy.array([gc_values[gcvl], gc_values[gcvr]])
278 gcindex = gc_values[gcvl+1:gcvr-1]
279 self.spc_out[ich,gcindex,ir] = numpy.interp(VelRange[gcindex],VelRange[interpindex],self.spc[ich,interpindex,ir])
280
281 dataOut.data_pre[0] = self.spc_out
282 #print ('Leaving RemoveWideGC ... ')
283 return dataOut
241
284
242 VelRange = dataOut.spc_range[2]
285 class SpectralFilters(Operation):
243 TimeRange = dataOut.spc_range[1]
286 ''' This class allows to replace the novalid values with noise for each channel
244 FrecRange = dataOut.spc_range[0]
287 This applies to CLAIRE RADAR
245
246 Vmax= 2*numpy.max(dataOut.spc_range[2])
247 Tmax= 2*numpy.max(dataOut.spc_range[1])
248 Fmax= 2*numpy.max(dataOut.spc_range[0])
249
250 Breaker1R=VelRange[numpy.abs(VelRange-(-LimitN)).argmin()]
251 Breaker1R=numpy.where(VelRange == Breaker1R)
252
253 Delta = self.Num_Bin/2 - Breaker1R[0]
254
255
256 '''Reacomodando SPCrange'''
257
258 VelRange=numpy.roll(VelRange,-(int(self.Num_Bin/2)) ,axis=0)
259
260 VelRange[-(int(self.Num_Bin/2)):]+= Vmax
261
288
262 FrecRange=numpy.roll(FrecRange,-(int(self.Num_Bin/2)),axis=0)
289 PositiveLimit : RightLimit of novalid data
290 NegativeLimit : LeftLimit of novalid data
263
291
264 FrecRange[-(int(self.Num_Bin/2)):]+= Fmax
292 Input:
265
293
266 TimeRange=numpy.roll(TimeRange,-(int(self.Num_Bin/2)),axis=0)
294 self.dataOut.data_pre : SPC and CSPC
295 self.dataOut.spc_range : To select wind and rainfall velocities
267
296
268 TimeRange[-(int(self.Num_Bin/2)):]+= Tmax
297 Affected:
269
298
270 ''' ------------------ '''
299 self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
271
300
272 Breaker2R=VelRange[numpy.abs(VelRange-(LimitR)).argmin()]
301 Written by D. Scipión 29.01.2021
273 Breaker2R=numpy.where(VelRange == Breaker2R)
302 '''
303 def __init__(self):
304 Operation.__init__(self)
305 self.i = 0
306
307 def run(self, dataOut, ):
274
308
309 self.spc = dataOut.data_pre[0].copy()
310 self.Num_Chn = self.spc.shape[0]
311 VelRange = dataOut.spc_range[2]
275
312
276 SPCroll = numpy.roll(self.spc,-(int(self.Num_Bin/2)) ,axis=1)
313 # novalid corresponds to data within the Negative and PositiveLimit
314
277
315
278 SPCcut = SPCroll.copy()
316 # Removing novalid data from the spectra
279 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,135 +338,186 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):
364
396 # print('Entering FitGau')
365 Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
397 # Assigning the variables
366
398 Vrange, ch, wnoise, num_intg, SNRlimit = X
367 SPCparam = []
399 # Noise Limits
368 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
400 noisebl = wnoise * 0.9
369 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
401 noisebh = wnoise * 1.1
370 SPC_ch1[:] = 0#numpy.NaN
402 # Radar Velocity
371 SPC_ch2[:] = 0#numpy.NaN
403 Va = max(Vrange)
372
404 deltav = Vrange[1] - Vrange[0]
373
405 x = numpy.arange(self.Num_Bin)
374
406
407 # print ('stop 0')
408
409 # 5 parameters, 2 Gaussians
410 DGauFitParam = numpy.zeros([5, self.Num_Hei,2])
411 DGauFitParam[:] = numpy.NaN
412
413 # SPCparam = []
414 # SPC_ch1 = numpy.zeros([self.Num_Bin,self.Num_Hei])
415 # SPC_ch2 = numpy.zeros([self.Num_Bin,self.Num_Hei])
416 # SPC_ch1[:] = 0 #numpy.NaN
417 # SPC_ch2[:] = 0 #numpy.NaN
418 # print ('stop 1')
375 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
389
447 # noise per channel.... we might want to use the noise at each range
390 wnoise = noise_ #/ spc_norm_max
448
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
398 minx=numpy.argmin(spc)
457 # print ('stop 2.2')
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
430
493 if powerwidth <= 1:
431 firstpeak=powerlo+powerwidth/10.# first gaussian energy location
494 # print('powerwidth <= 1')
432 secondpeak=powerhi-powerwidth/10.#second gaussian energy location
495 continue
433 midpeak=(firstpeak+secondpeak)/2.
496
434 firstamp=spcs[int(firstpeak)]
497 # print ('stop 6')
435 secondamp=spcs[int(secondpeak)]
498 firstpeak = powerlo + powerwidth/10.# first gaussian energy location
436 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)]
437
504
438 x=numpy.arange( self.Num_Bin )
505 y_data = spc + wnoise
439 y_data=spc+wnoise
440
506
441 ''' single Gaussian '''
507 ''' single Gaussian '''
442 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
508 shift0 = numpy.mod(midpeak+minx, self.Num_Bin )
443 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
509 width0 = powerwidth/4.#Initialization entire power of spectrum divided by 4
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)
449
515 # print ('stop 7.1')
450 chiSq1=lsq1[1];
516 # print (bnds)
451
517
452
518 chiSq1=lsq1[1]
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]
@@ -463,127 +531,142 class GaussianFit(Operation):
463 noise=lsq1[0][4]
531 noise=lsq1[0][4]
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
494 else:
565 else:
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:
503 if a1>0.7*a2:
574 if a1>0.7*a2:
504 choice=1
575 choice = 1
505 else:
576 else:
506 choice=2
577 choice = 2
507 elif gp2>gp1:
578 elif gp2>gp1:
508 if a2>0.7*a1:
579 if a2>0.7*a1:
509 choice=2
580 choice = 2
510 else:
581 else:
511 choice=1
582 choice = 1
512 else:
583 else:
513 choice=numpy.argmax([a1,a2])+1
584 choice = numpy.argmax([a1,a2])+1
514 #else:
585 #else:
515 #choice=argmin([std2a,std2b])+1
586 #choice=argmin([std2a,std2b])+1
516
587
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
520
591 # print ('stop 14')
521 shift0=lsq2[0][0];
592 shift0 = lsq2[0][0]
522 vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
593 vel0 = Vrange[0] + shift0 * deltav
523 shift1=lsq2[0][4];
594 shift1 = lsq2[0][4]
524 vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
595 # vel1=Vrange[0] + shift1 * deltav
525
596
526 max_vel = 1.0
597 # max_vel = 1.0
527
598 # Va = max(Vrange)
599 # deltav = Vrange[1]-Vrange[0]
600 # print ('stop 15')
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
530 shift0=lsq2[0][0]
603 if vel0 > -Va and vel0 < Va : #first peak is in the correct range
531 width0=lsq2[0][1]
604 shift0 = lsq2[0][0]
532 Amplitude0=lsq2[0][2]
605 width0 = lsq2[0][1]
533 p0=lsq2[0][3]
606 Amplitude0 = lsq2[0][2]
534
607 p0 = lsq2[0][3]
535 shift1=lsq2[0][4]
608
536 width1=lsq2[0][5]
609 shift1 = lsq2[0][4]
537 Amplitude1=lsq2[0][6]
610 width1 = lsq2[0][5]
538 p1=lsq2[0][7]
611 Amplitude1 = lsq2[0][6]
539 noise=lsq2[0][8]
612 p1 = lsq2[0][7]
613 noise = lsq2[0][8]
540 else:
614 else:
541 shift1=lsq2[0][0]
615 shift1 = lsq2[0][0]
542 width1=lsq2[0][1]
616 width1 = lsq2[0][1]
543 Amplitude1=lsq2[0][2]
617 Amplitude1 = lsq2[0][2]
544 p1=lsq2[0][3]
618 p1 = lsq2[0][3]
545
619
546 shift0=lsq2[0][4]
620 shift0 = lsq2[0][4]
547 width0=lsq2[0][5]
621 width0 = lsq2[0][5]
548 Amplitude0=lsq2[0][6]
622 Amplitude0 = lsq2[0][6]
549 p0=lsq2[0][7]
623 p0 = lsq2[0][7]
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
656 model0u = amplitude0*numpy.exp(-0.5*abs((x - shift0 - self.Num_Bin)/width0)**power0)
569 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)
570
658 return model0 + model0u + model0d + noise
571 model0d=amplitude0*numpy.exp(-0.5*abs((x-shift0+ self.Num_Bin )/width0)**power0)
572 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
663 model0u = amplitude0*numpy.exp(-0.5*abs((x - shift0 - self.Num_Bin)/width0)**power0)
578 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)
579
665
580 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)
581 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)
582
668 model1d = amplitude1*numpy.exp(-0.5*abs((x - shift1 + self.Num_Bin)/width1)**power1)
583 model1u=amplitude1*numpy.exp(-0.5*abs((x-shift1- self.Num_Bin )/width1)**power1)
669 return model0 + model0u + model0d + model1 + model1u + model1d + noise
584
585 model1d=amplitude1*numpy.exp(-0.5*abs((x-shift1+ self.Num_Bin )/width1)**power1)
586 return model0+model0u+model0d+model1+model1u+model1d+noise
587
670
588 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.
589
672
@@ -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
@@ -670,104 +731,73 class PrecipitationProc(Operation):
670 self.Lambda = Lambda
731 self.Lambda = Lambda
671 self.aL = aL
732 self.aL = aL
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 #
679
743 ExpConstant = 10**(40/10) #Constante Experimental
680 ''' ============================= '''
744
681
745 SignalPower = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
682 self.spc[0] = (self.spc[0]-dataOut.noise[0])
746 for i in range(self.Num_Chn):
683 self.spc[1] = (self.spc[1]-dataOut.noise[1])
747 SignalPower[i,:,:] = self.spc[i,:,:] - dataOut.noise[i]
684 self.spc[2] = (self.spc[2]-dataOut.noise[2])
748 SignalPower[numpy.where(SignalPower < 0)] = 1e-20
685
749
686 self.spc[ numpy.where(self.spc < 0)] = 0
750 SPCmean = numpy.mean(SignalPower, 0)
687
751 Pr = SPCmean[:,:]/dataOut.normFactor
688 SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise))
752
689 SPCmean[ numpy.where(SPCmean < 0)] = 0
753 # Declaring auxiliary variables
690
754 Range = dataOut.heightList*1000. #Range in m
691 ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei])
755 # replicate the heightlist to obtain a matrix [Num_Bin,Num_Hei]
692 ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei])
756 rMtrx = numpy.transpose(numpy.transpose([dataOut.heightList*1000.] * self.Num_Bin))
693 ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei])
757 zMtrx = rMtrx+Altitude
694
758 # replicate the VelRange to obtain a matrix [Num_Bin,Num_Hei]
695 Pr = SPCmean[:,:]
759 VelMtrx = numpy.transpose(numpy.tile(VelRange[:-1], (self.Num_Hei,1)))
696
760
697 VelMeteoro = numpy.mean(SPCmean,axis=0)
761 # height dependence to air density Foote and Du Toit (1969)
698
762 delv_z = 1 + 3.68e-5 * zMtrx + 1.71e-9 * zMtrx**2
699 D_range = numpy.zeros([self.Num_Bin,self.Num_Hei])
763 VMtrx = VelMtrx / delv_z #Normalized velocity
700 SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei])
764 VMtrx[numpy.where(VMtrx> 9.6)] = numpy.NaN
701 N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei])
765 # Diameter is related to the fall speed of falling drops
702 V_mean = numpy.zeros(self.Num_Hei)
766 D_Vz = -1.667 * numpy.log( 0.9369 - 0.097087 * VMtrx ) # D in [mm]
703 del_V = numpy.zeros(self.Num_Hei)
767 # Only valid for D>= 0.16 mm
704 Z = numpy.zeros(self.Num_Hei)
768 D_Vz[numpy.where(D_Vz < 0.16)] = numpy.NaN
705 Ze = numpy.zeros(self.Num_Hei)
769
706 RR = numpy.zeros(self.Num_Hei)
770 #Calculate Radar Reflectivity ETAn
707
771 ETAn = (RadarConstant *ExpConstant) * Pr * rMtrx**2 #Reflectivity (ETA)
708 Range = dataOut.heightList*1000.
772 ETAd = ETAn * 6.18 * exp( -0.6 * D_Vz ) * delv_z
709
773 # Radar Cross Section
710 for R in range(self.Num_Hei):
774 sigmaD = Km2 * (D_Vz * 1e-3 )**6 * numpy.pi**5 / Lambda**4
711
775 # Drop Size Distribution
712 h = Range[R] + Altitude #Range from ground to radar pulse altitude
776 DSD = ETAn / sigmaD
713 del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
777 # Equivalente Reflectivy
714
778 Ze_eqn = numpy.nansum( DSD * D_Vz**6 ,axis=0)
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
779 Ze_org = numpy.nansum(ETAn * Lambda**4, axis=0) / (1e-18*numpy.pi**5 * Km2) # [mm^6 /m^3]
716
780 # RainFall Rate
717 '''NOTA: ETA(n) dn = ETA(f) df
781 RR = 0.0006*numpy.pi * numpy.nansum( D_Vz**3 * DSD * VelMtrx ,0) #mm/hr
718
782
719 dn = 1 Diferencial de muestreo
783 # Censoring the data
720 df = ETA(n) / ETA(f)
784 # Removing data with SNRth < 0dB se debe considerar el SNR por canal
721
785 SNRth = 10**(SNRdBlimit/10) #-30dB
722 '''
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
723
787 W = numpy.nanmean(dataOut.data_dop,0)
724 ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA)
788 W[novalid] = numpy.NaN
725
789 Ze_org[novalid] = numpy.NaN
726 ETAv[:,R]=ETAn[:,R]/dV
790 RR[novalid] = numpy.NaN
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
891 dataOut.data_snr[numpy.where( dataOut.data_snr <0 )] = 1e-20
892
893
907
894 data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN
908 dbSNR = 10*numpy.log10(numpy.average(dataOut.data_snr,0))
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
1028
1008 nChan = spc.shape[0]
1029 SPC_Samples = numpy.ones([spc.shape[0],spc.shape[1]]) # for normalized spc values for one height
1009 nProf = spc.shape[1]
1030 phase = numpy.ones([spc.shape[0],spc.shape[1]]) # phase between channels
1010 nPair = cspc.shape[0]
1031 CSPC_Samples = numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_) # for normalized cspc values
1011
1032 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
1033 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
1034 xFrec = AbbsisaRange[0][0:spc.shape[1]] # frequency range
1014 phase = numpy.zeros([nPair, nProf]) # phase between channels
1035 xVel = AbbsisaRange[2][0:spc.shape[1]] # velocity range
1015 PhaseSlope = numpy.zeros(nPair) # slope of the phases, channelwise
1036 SPCav = numpy.average(spc, axis=0)-numpy.average(noise) # spc[0]-noise[0]
1016 PhaseInter = numpy.zeros(nPair) # intercept to the slope of the phases, channelwise
1037
1017 xFrec = AbbsisaRange[0][:-1] # frequency range
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)
1018 xVel = AbbsisaRange[2][:-1] # velocity range
1039 CSPCmoments = []
1019 xSamples = xFrec # the frequency range is taken
1040
1020 delta_x = xSamples[1] - xSamples[0] # delta_f or delta_x
1021
1022 # only consider velocities with in NegativeLimit and PositiveLimit
1023 if (NegativeLimit is None):
1024 NegativeLimit = numpy.min(xVel)
1025 if (PositiveLimit is None):
1026 PositiveLimit = numpy.max(xVel)
1027 xvalid = numpy.where((xVel > NegativeLimit) & (xVel < PositiveLimit))
1028 xSamples_zoom = xSamples[xvalid]
1041
1029
1042 '''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)
1056 popt = [1e-10,0,1e-10]
1079
1057 # Spectra average
1080 popt = [1e-10,0,1e-10]
1058 SPCMean = numpy.average(SPC_Samples,0)
1081 SPCMoments = self.Moments(SPCMean, xSamples)
1059 # Moments in frequency
1082
1060 SPCMoments = self.Moments(SPCMean[xvalid], xSamples_zoom)
1083 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
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 '''
@@ -1393,13 +1276,13 class SpectralMoments(Operation):
1393 max_spec = aux.max()
1276 max_spec = aux.max()
1394 m = aux.tolist().index(max_spec)
1277 m = aux.tolist().index(max_spec)
1395
1278
1396 #Smooth
1279 # Smooth
1397 if (smooth == 0):
1280 if (smooth == 0):
1398 spec2 = spec
1281 spec2 = spec
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,14 +1308,17 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
1428 power = ((spec2[valid] - n0) * fwindow[valid]).sum()
1311 signal_power = ((spec2[valid] - n0) * fwindow[valid]).mean() # D. Scipión added with correct definition
1312 total_power = (spec2[valid] * fwindow[valid]).mean() # D. Scipión added with correct definition
1313 power = ((spec2[valid] - n0) * fwindow[valid]).sum()
1429 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)
1431 snr = (spec2.mean()-n0)/n0
1316 snr = (spec2.mean()-n0)/n0
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