##// END OF EJS Templates
schain
RSS
Update SpectralParameters
Danny Scipión -
r1345:fa1ffba10c02
parent child
Show More
Side by Side Unified
Context file:
r1345:fa1ffba10c02 Loading diff...:
Show file before | Show file after | Hide comments
@@ -129,7 +129,7 class ParametersProc(ProcessingUnit):
129
129
130 if self.dataIn.type == "Spectra":
130 if self.dataIn.type == "Spectra":
131
131
132 self.dataOut.data_pre = (self.dataIn.data_spc, self.dataIn.data_cspc)
132 self.dataOut.data_pre = [self.dataIn.data_spc, self.dataIn.data_cspc]
133 self.dataOut.data_spc = self.dataIn.data_spc
133 self.dataOut.data_spc = self.dataIn.data_spc
134 self.dataOut.data_cspc = self.dataIn.data_cspc
134 self.dataOut.data_cspc = self.dataIn.data_cspc
135 self.dataOut.nProfiles = self.dataIn.nProfiles
135 self.dataOut.nProfiles = self.dataIn.nProfiles
@@ -199,13 +199,11 def target(tups):
199
199
200 return obj.FitGau(args)
200 return obj.FitGau(args)
201
201
202 class RemoveWideGC(Operation):
203 ''' This class remove the wide clutter and replace it with a simple interpolation points
204 This mainly applies to CLAIRE radar
202
205
203 class SpectralFilters(Operation):
206 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
207
210 Input:
208 Input:
211
209
@@ -215,91 +213,103 class SpectralFilters(Operation):
215 Affected:
213 Affected:
216
214
217 self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
215 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
216
217 Written by D. Scipión 25.02.2021
222 '''
218 '''
223
224 def __init__(self):
219 def __init__(self):
225 Operation.__init__(self)
220 Operation.__init__(self)
226 self.i=0
221 self.i = 0
227
222 self.ich = 0
228 def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5):
223 self.ir = 0
229
224
230
225 def run(self, dataOut, ClutterWidth=2.5):
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[1]
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.where(numpy.abs(VelRange) <= ClutterWidth)
236
237 # Removing novalid data from the spectra
238 for ich in range(self.Num_Chn):
239 # Estimate the noise at aech range
240
241 for ir in range(self.Num_Hei):
242 # Estimate the noise at aech range
243 HSn = hildebrand_sekhon(self.spc[ich,:,ir],dataOut.nIncohInt)
244 # Removing the noise floor at each range
245 novalid = numpy.where(self.spc[ich,:,ir] < HSn)
246 self.spc[novalid,ir] = HSn
247
248 junk = [HSn, numpy.transpose(self.spc[ich,gc_values,ir]), HSn]
249 j1index = numpy.where(numpy.diff[junk]>0)
250 j2index = numpy.where(numpy.diff[junk]<0)
251 junk3 = numpy.diff(j1index)
252 junk4 = numpy.diff(j2index)
253
254 valleyindex = j2index[junk4>1]
255 peakindex = j1index[junk3>1]
256
257 # Identify the clutter (close to zero)
258 isvalid = numpy.where(where.abs(VelRange[gc_values[peakindex]]) <= 2.5*dv)
259 # if isempty(isvalid)
260 # continue
261 # if numel(isvalid)>1
262 # [~,vindex]= numpy.max(spc[gc_values[peakindex[isvalid]],ir])
263 # isvalid = isvalid[vindex]
264 # clutter peak
265 gcpeak = peakindex[isvalid]
266 # left and right index of the clutter
267 gcvl = valleyindex[numpy.where(valleyindex < gcpeak, 1, 'last' )]
268 gcvr = valleyindex[numpy.where(valleyindex > gcpeak, 1)]
269
270 # Removing the clutter
271 interpindex = [gc_values(gcvl), gc_values(gcvr)]
272 gcindex = gc_values[gcvl+1:gcvr-1]
273 self.spc_out[ich,gcindex,ir] = numpy.interp(VelRange[gcindex],VelRange[interpindex],self.spc[ich,interpindex,ir])
274
275 dataOut.data_pre[0] = self.spc_out
276 return dataOut
241
277
242 VelRange = dataOut.spc_range[2]
278 class SpectralFilters(Operation):
243 TimeRange = dataOut.spc_range[1]
279 ''' This class allows to replace the novalid values with noise for each channel
244 FrecRange = dataOut.spc_range[0]
280 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
281
262 FrecRange=numpy.roll(FrecRange,-(int(self.Num_Bin/2)),axis=0)
282 PositiveLimit : RightLimit of novalid data
283 NegativeLimit : LeftLimit of novalid data
263
284
264 FrecRange[-(int(self.Num_Bin/2)):]+= Fmax
285 Input:
265
286
266 TimeRange=numpy.roll(TimeRange,-(int(self.Num_Bin/2)),axis=0)
287 self.dataOut.data_pre : SPC and CSPC
288 self.dataOut.spc_range : To select wind and rainfall velocities
267
289
268 TimeRange[-(int(self.Num_Bin/2)):]+= Tmax
290 Affected:
269
291
270 ''' ------------------ '''
292 self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
271
293
272 Breaker2R=VelRange[numpy.abs(VelRange-(LimitR)).argmin()]
294 Written by D. Scipión 29.01.2021
273 Breaker2R=numpy.where(VelRange == Breaker2R)
295 '''
296 def __init__(self):
297 Operation.__init__(self)
298 self.i = 0
299
300 def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=-1.5):
274
301
302 self.spc = dataOut.data_pre[0].copy()
303 self.Num_Chn = self.spc.shape[0]
304 VelRange = dataOut.spc_range[2]
275
305
276 SPCroll = numpy.roll(self.spc,-(int(self.Num_Bin/2)) ,axis=1)
306 # novalid corresponds to data within the Negative and PositiveLimit
307 novalid = numpy.where((VelRange[:-1] >= NegativeLimit) & (VelRange[:-1] <= PositiveLimit))
277
308
278 SPCcut = SPCroll.copy()
309 # Removing novalid data from the spectra
279 for i in range(self.Num_Chn):
310 for i in range(self.Num_Chn):
280
311 self.spc[i,novalid,:] = dataOut.noise[i]
281 SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i]
312 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
313 return dataOut
304
314
305 class GaussianFit(Operation):
315 class GaussianFit(Operation):
@@ -367,8 +377,8 class GaussianFit(Operation):
367 SPCparam = []
377 SPCparam = []
368 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
378 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
369 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
379 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
370 SPC_ch1[:] = 0#numpy.NaN
380 SPC_ch1[:] = 0 #numpy.NaN
371 SPC_ch2[:] = 0#numpy.NaN
381 SPC_ch2[:] = 0 #numpy.NaN
372
382
373
383
374
384
@@ -614,31 +624,8 class PrecipitationProc(Operation):
614 Operation.__init__(self)
624 Operation.__init__(self)
615 self.i=0
625 self.i=0
616
626
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,
627 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):
628 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km2 = 0.93, Altitude=3350):
635
636
637 Velrange = dataOut.spcparam_range[2]
638 FrecRange = dataOut.spcparam_range[0]
639
640 dV= Velrange[1]-Velrange[0]
641 dF= FrecRange[1]-FrecRange[0]
642
629
643 if radar == "MIRA35C" :
630 if radar == "MIRA35C" :
644
631
@@ -650,18 +637,17 class PrecipitationProc(Operation):
650
637
651 else:
638 else:
652
639
653 self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() #
640 self.spc = dataOut.data_pre[0].copy()
654
655 """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
656
641
642 #NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX
657 self.spc[:,:,0:7]= numpy.NaN
643 self.spc[:,:,0:7]= numpy.NaN
658
644
659 """##########################################"""
660
661 self.Num_Hei = self.spc.shape[2]
645 self.Num_Hei = self.spc.shape[2]
662 self.Num_Bin = self.spc.shape[1]
646 self.Num_Bin = self.spc.shape[1]
663 self.Num_Chn = self.spc.shape[0]
647 self.Num_Chn = self.spc.shape[0]
664
648
649 VelRange = dataOut.spc_range[2]
650
665 ''' Se obtiene la constante del RADAR '''
651 ''' Se obtiene la constante del RADAR '''
666
652
667 self.Pt = Pt
653 self.Pt = Pt
@@ -670,104 +656,72 class PrecipitationProc(Operation):
670 self.Lambda = Lambda
656 self.Lambda = Lambda
671 self.aL = aL
657 self.aL = aL
672 self.tauW = tauW
658 self.tauW = tauW
673 self.ThetaT = ThetaT
659 self.ThetaT = ThetaT
674 self.ThetaR = ThetaR
660 self.ThetaR = ThetaR
661 self.GSys = 10**(36.63/10) # Ganancia de los LNA 36.63 dB
662 self.lt = 10**(1.67/10) # Perdida en cables Tx 1.67 dB
663 self.lr = 10**(5.73/10) # Perdida en cables Rx 5.73 dB
675
664
676 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
665 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)
666 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
678 RadarConstant = 10e-26 * Numerator / Denominator #
667 RadarConstant = 10e-26 * Numerator / Denominator #
679
668 ExpConstant = 10**(40/10) #Constante Experimental
680 ''' ============================= '''
669
681
670 SignalPower = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
682 self.spc[0] = (self.spc[0]-dataOut.noise[0])
671 for i in range(self.Num_Chn):
683 self.spc[1] = (self.spc[1]-dataOut.noise[1])
672 SignalPower[i,:,:] = self.spc[i,:,:] - dataOut.noise[i]
684 self.spc[2] = (self.spc[2]-dataOut.noise[2])
673 SignalPower[numpy.where(SignalPower < 0)] = 1e-20
685
674
686 self.spc[ numpy.where(self.spc < 0)] = 0
675 SPCmean = numpy.mean(SignalPower, 0)
687
676
688 SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise))
677 Pr = SPCmean[:,:]/dataOut.normFactor
689 SPCmean[ numpy.where(SPCmean < 0)] = 0
678
690
679 # Declaring auxiliary variables
691 ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei])
680 Range = dataOut.heightList*1000. #Range in m
692 ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei])
681 # replicate the heightlist to obtain a matrix [Num_Bin,Num_Hei]
693 ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei])
682 rMtrx = numpy.transpose(numpy.transpose([dataOut.heightList*1000.] * self.Num_Bin))
694
683 zMtrx = rMtrx+Altitude
695 Pr = SPCmean[:,:]
684 # replicate the VelRange to obtain a matrix [Num_Bin,Num_Hei]
696
685 VelMtrx = numpy.transpose(numpy.tile(VelRange[:-1], (self.Num_Hei,1)))
697 VelMeteoro = numpy.mean(SPCmean,axis=0)
686
698
687 # height dependence to air density Foote and Du Toit (1969)
699 D_range = numpy.zeros([self.Num_Bin,self.Num_Hei])
688 delv_z = 1 + 3.68e-5 * zMtrx + 1.71e-9 * zMtrx**2
700 SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei])
689 VMtrx = VelMtrx / delv_z #Normalized velocity
701 N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei])
690 VMtrx[numpy.where(VMtrx> 9.6)] = numpy.NaN
702 V_mean = numpy.zeros(self.Num_Hei)
691 # Diameter is related to the fall speed of falling drops
703 del_V = numpy.zeros(self.Num_Hei)
692 D_Vz = -1.667 * numpy.log( 0.9369 - 0.097087 * VMtrx ) # D in [mm]
704 Z = numpy.zeros(self.Num_Hei)
693 # Only valid for D>= 0.16 mm
705 Ze = numpy.zeros(self.Num_Hei)
694 D_Vz[numpy.where(D_Vz < 0.16)] = numpy.NaN
706 RR = numpy.zeros(self.Num_Hei)
695
707
696 #Calculate Radar Reflectivity ETAn
708 Range = dataOut.heightList*1000.
697 ETAn = (RadarConstant *ExpConstant) * Pr * rMtrx**2 #Reflectivity (ETA)
709
698 ETAd = ETAn * 6.18 * exp( -0.6 * D_Vz ) * delv_z
710 for R in range(self.Num_Hei):
699 # Radar Cross Section
711
700 sigmaD = Km2 * (D_Vz * 1e-3 )**6 * numpy.pi**5 / Lambda**4
712 h = Range[R] + Altitude #Range from ground to radar pulse altitude
701 # Drop Size Distribution
713 del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
702 DSD = ETAn / sigmaD
714
703 # Equivalente Reflectivy
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
704 Ze_eqn = numpy.nansum( DSD * D_Vz**6 ,axis=0)
716
705 Ze_org = numpy.nansum(ETAn * Lambda**4, axis=0) / (1e-18*numpy.pi**5 * Km2) # [mm^6 /m^3]
717 '''NOTA: ETA(n) dn = ETA(f) df
706 # RainFall Rate
718
707 RR = 0.0006*numpy.pi * numpy.nansum( D_Vz**3 * DSD * VelMtrx ,0) #mm/hr
719 dn = 1 Diferencial de muestreo
708
720 df = ETA(n) / ETA(f)
709 # Censoring the data
721
710 # Removing data with SNRth < 0dB se debe considerar el SNR por canal
722 '''
711 SNRth = 10**(-30/10) #-20dB
723
712 novalid = numpy.where((dataOut.data_SNR[0,:] <SNRth) | (dataOut.data_SNR[1,:] <SNRth) | (dataOut.data_SNR[2,:] <SNRth)) # AND condition. Maybe OR condition better
724 ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA)
713 W = numpy.nanmean(dataOut.data_DOP,0)
725
714 W[novalid] = numpy.NaN
726 ETAv[:,R]=ETAn[:,R]/dV
715 Ze_org[novalid] = numpy.NaN
727
716 RR[novalid] = numpy.NaN
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
717
763 dataOut.data_output = RR[8]
718 dataOut.data_output = RR[8]
764 dataOut.data_param = numpy.ones([3,self.Num_Hei])
719 dataOut.data_param = numpy.ones([3,self.Num_Hei])
765 dataOut.channelList = [0,1,2]
720 dataOut.channelList = [0,1,2]
766
721
767 dataOut.data_param[0]=dBZ
722 dataOut.data_param[0]=10*numpy.log10(Ze_org)
768 dataOut.data_param[1]=V_mean
723 dataOut.data_param[1]=W
769 dataOut.data_param[2]=RR
724 dataOut.data_param[2]=RR
770
771 return dataOut
725 return dataOut
772
726
773 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
727 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
@@ -784,7 +738,7 class PrecipitationProc(Operation):
784
738
785 ETA = numpy.sum(SNR,1)
739 ETA = numpy.sum(SNR,1)
786
740
787 ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
741 ETA = numpy.where(ETA != 0. , ETA, numpy.NaN)
788
742
789 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
743 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
790
744
@@ -832,7 +786,7 class FullSpectralAnalysis(Operation):
832
786
833 Output:
787 Output:
834
788
835 self.dataOut.data_output : Zonal wind, Meridional wind and Vertical wind
789 self.dataOut.data_output : Zonal wind, Meridional wind, and Vertical wind
836
790
837
791
838 Parameters affected: Winds, height range, SNR
792 Parameters affected: Winds, height range, SNR
@@ -848,7 +802,7 class FullSpectralAnalysis(Operation):
848 """Erick: NOTE THE RANGE OF THE PULSE TX MUST BE REMOVED"""
802 """Erick: NOTE THE RANGE OF THE PULSE TX MUST BE REMOVED"""
849
803
850 SNRspc = spc.copy()
804 SNRspc = spc.copy()
851 SNRspc[:,:,0:7]= numpy.NaN
805 SNRspc[:,:,0:7]= numpy.NaN # D. Scipión... the cleaning should not be hardwired in the code... it needs to be flexible... NEEDS TO BE REMOVED
852
806
853 """##########################################"""
807 """##########################################"""
854
808
@@ -920,7 +874,7 class FullSpectralAnalysis(Operation):
920 velocityY=numpy.append(velocityY, numpy.NaN)
874 velocityY=numpy.append(velocityY, numpy.NaN)
921
875
922 if dbSNR[Height] > SNRlimit:
876 if dbSNR[Height] > SNRlimit:
923 velocityV=numpy.append(velocityV, -Vver) # reason for this minus sign -> convention? (taken from Ericks version)
877 velocityV=numpy.append(velocityV, -Vver) # reason for this minus sign -> convention? (taken from Ericks version) D.S. yes!
924 else:
878 else:
925 velocityV=numpy.append(velocityV, numpy.NaN)
879 velocityV=numpy.append(velocityV, numpy.NaN)
926
880
@@ -944,19 +898,12 class FullSpectralAnalysis(Operation):
944 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
898 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
945
899
946 def Moments(self, ySamples, xSamples):
900 def Moments(self, ySamples, xSamples):
947 '''***
901 Power = numpy.nanmean(ySamples) # Power, 0th Moment
948 Variables corresponding to moments of distribution.
902 yNorm = ySamples / numpy.nansum(ySamples)
949 Also used as initial coefficients for curve_fit.
903 RadVel = numpy.nansum(xSamples * yNorm) # Radial Velocity, 1st Moment
950 Vr was corrected. Only a velocity when x is velocity, of course.
904 Sigma2 = abs(numpy.nansum(yNorm * (xSamples - RadVel)**2)) # Spectral Width, 2nd Moment
951 ***'''
905 StdDev = Sigma2**0.5 # Desv. Estandar, Ancho espectral
952 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
906 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
907
961 def StopWindEstimation(self, error_code):
908 def StopWindEstimation(self, error_code):
962 '''
909 '''
@@ -1033,12 +980,13 class FullSpectralAnalysis(Operation):
1033 PhaseInter = numpy.ones(spc.shape[0]) # intercept to the slope of the phases, channelwise
980 PhaseInter = numpy.ones(spc.shape[0]) # intercept to the slope of the phases, channelwise
1034 xFrec = AbbsisaRange[0][0:spc.shape[1]] # frequency range
981 xFrec = AbbsisaRange[0][0:spc.shape[1]] # frequency range
1035 xVel = AbbsisaRange[2][0:spc.shape[1]] # velocity range
982 xVel = AbbsisaRange[2][0:spc.shape[1]] # velocity range
1036 SPCav = numpy.average(spc, axis=0)-numpy.average(noise) # spc[0]-noise[0]
983 SPCav = numpy.average(spc, axis=0)-numpy.average(noise) # spc[0]-noise[0] %D.S. why??? I suggest only spc....
1037
984
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)
985 SPCmoments_vel = self.Moments(SPCav, xVel ) # SPCmoments_vel[1] corresponds to vertical velocity and is used to determine if signal corresponds to wind (if .. <3)
986 # D.S. I suggest to each moment to be calculated independently, because the signal level y/o interferences are not the same in all channels and
987 # signal or SNR seems to be contaminated
1039 CSPCmoments = []
988 CSPCmoments = []
1040
989
1041
1042 '''Getting Eij and Nij'''
990 '''Getting Eij and Nij'''
1043
991
1044 Xi01, Xi02, Xi12 = ChanDist[:,0]
992 Xi01, Xi02, Xi12 = ChanDist[:,0]
@@ -1052,10 +1000,13 class FullSpectralAnalysis(Operation):
1052 spc_norm = spc.copy() # need copy() because untouched spc is needed for normalization of cspc below
1000 spc_norm = spc.copy() # need copy() because untouched spc is needed for normalization of cspc below
1053 spc_norm = numpy.where(numpy.isfinite(spc_norm), spc_norm, numpy.NAN)
1001 spc_norm = numpy.where(numpy.isfinite(spc_norm), spc_norm, numpy.NAN)
1054
1002
1003 # D. Scipión: It is necessary to define DeltaF or DeltaV... it is wrong to use Factor_Norm. It's constant... not a variable
1004
1005 # For each channel
1055 for i in range(spc.shape[0]):
1006 for i in range(spc.shape[0]):
1056
1007
1057 spc_sub = spc_norm[i,:] - noise[i] # spc not smoothed here or in previous version.
1008 spc_sub = spc_norm[i,:] - noise[i] # spc not smoothed here or in previous version.
1058
1009 # D. Scipión: Factor_Norm has to be replaced by DeltaF or DeltaV - It's a constant
1059 Factor_Norm = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc_sub)) # usually = Freq range / nfft
1010 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)
1011 normalized_spc = spc_sub / (numpy.nansum(numpy.abs(spc_sub)) * Factor_Norm)
1061
1012
@@ -1111,6 +1062,8 class FullSpectralAnalysis(Operation):
1111 '''
1062 '''
1112
1063
1113 radarWavelength = 0.6741 # meters
1064 radarWavelength = 0.6741 # meters
1065 # D.S. This does not need to hardwired... It has to be in function of the radar frequency
1066
1114 count_limit_freq = numpy.abs(popt[1]) + widthlimit # Hz, m/s can be also used if velocity is desired abscissa.
1067 count_limit_freq = numpy.abs(popt[1]) + widthlimit # Hz, m/s can be also used if velocity is desired abscissa.
1115 # count_limit_freq = numpy.max(xFrec)
1068 # count_limit_freq = numpy.max(xFrec)
1116
1069
@@ -1393,13 +1346,13 class SpectralMoments(Operation):
1393 max_spec = aux.max()
1346 max_spec = aux.max()
1394 m = aux.tolist().index(max_spec)
1347 m = aux.tolist().index(max_spec)
1395
1348
1396 #Smooth
1349 # Smooth
1397 if (smooth == 0):
1350 if (smooth == 0):
1398 spec2 = spec
1351 spec2 = spec
1399 else:
1352 else:
1400 spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
1353 spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
1401
1354
1402 # Calculo de Momentos
1355 # Moments Estimation
1403 bb = spec2[numpy.arange(m,spec2.size)]
1356 bb = spec2[numpy.arange(m,spec2.size)]
1404 bb = (bb<n0).nonzero()
1357 bb = (bb<n0).nonzero()
1405 bb = bb[0]
1358 bb = bb[0]
@@ -1425,14 +1378,17 class SpectralMoments(Operation):
1425
1378
1426 valid = numpy.arange(int(m + bb0 - ss1 + 1)) + ss1
1379 valid = numpy.arange(int(m + bb0 - ss1 + 1)) + ss1
1427
1380
1428 power = ((spec2[valid] - n0) * fwindow[valid]).sum()
1381 signal_power = ((spec2[valid] - n0) * fwindow[valid]).mean() # D. Scipión added with correct definition
1382 total_power = (spec2[valid] * fwindow[valid]).mean() # D. Scipión added with correct definition
1383 power = ((spec2[valid] - n0) * fwindow[valid]).sum()
1429 fd = ((spec2[valid]- n0)*freq[valid] * fwindow[valid]).sum() / power
1384 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)
1385 w = numpy.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum() / power)
1431 snr = (spec2.mean()-n0)/n0
1386 snr = (spec2.mean()-n0)/n0
1432 if (snr < 1.e-20) :
1387 if (snr < 1.e-20) :
1433 snr = 1.e-20
1388 snr = 1.e-20
1434
1389
1435 vec_power[ind] = power
1390 # vec_power[ind] = power #D. Scipión replaced with the line below
1391 vec_power[ind] = total_power
1436 vec_fd[ind] = fd
1392 vec_fd[ind] = fd
1437 vec_w[ind] = w
1393 vec_w[ind] = w
1438 vec_snr[ind] = snr
1394 vec_snr[ind] = snr
General Comments 0
You need to be logged in to leave comments. Login now