##// 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 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 133 self.dataOut.data_spc = self.dataIn.data_spc
134 134 self.dataOut.data_cspc = self.dataIn.data_cspc
135 135 self.dataOut.nProfiles = self.dataIn.nProfiles
@@ -199,13 +199,11 def target(tups):
199 199
200 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):
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
206 ClutterWidth : Width to look for the clutter peak
209 207
210 208 Input:
211 209
@@ -215,91 +213,103 class SpectralFilters(Operation):
215 213 Affected:
216 214
217 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 219 def __init__(self):
225 220 Operation.__init__(self)
226 self.i=0
227
228 def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5):
229
230
231 #Limite de vientos
232 LimitR = PositiveLimit
233 LimitN = NegativeLimit
221 self.i = 0
222 self.ich = 0
223 self.ir = 0
224
225 def run(self, dataOut, ClutterWidth=2.5):
234 226
235 227 self.spc = dataOut.data_pre[0].copy()
236 self.cspc = dataOut.data_pre[1].copy()
237
238 self.Num_Hei = self.spc.shape[2]
239 self.Num_Bin = self.spc.shape[1]
228 self.spc_out = dataOut.data_pre[0].copy()
240 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]
243 TimeRange = dataOut.spc_range[1]
244 FrecRange = dataOut.spc_range[0]
245
246 Vmax= 2*numpy.max(dataOut.spc_range[2])
247 Tmax= 2*numpy.max(dataOut.spc_range[1])
248 Fmax= 2*numpy.max(dataOut.spc_range[0])
249
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
278 class SpectralFilters(Operation):
279 ''' This class allows to replace the novalid values with noise for each channel
280 This applies to CLAIRE RADAR
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()]
273 Breaker2R=numpy.where(VelRange == Breaker2R)
294 Written by D. Scipión 29.01.2021
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 310 for i in range(self.Num_Chn):
280
281 SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i]
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
311 self.spc[i,novalid,:] = dataOut.noise[i]
312 dataOut.data_pre[0] = self.spc
303 313 return dataOut
304 314
305 315 class GaussianFit(Operation):
@@ -367,8 +377,8 class GaussianFit(Operation):
367 377 SPCparam = []
368 378 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
369 379 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
370 SPC_ch1[:] = 0#numpy.NaN
371 SPC_ch2[:] = 0#numpy.NaN
380 SPC_ch1[:] = 0 #numpy.NaN
381 SPC_ch2[:] = 0 #numpy.NaN
372 382
373 383
374 384
@@ -614,31 +624,8 class PrecipitationProc(Operation):
614 624 Operation.__init__(self)
615 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 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):
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]
628 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km2 = 0.93, Altitude=3350):
642 629
643 630 if radar == "MIRA35C" :
644 631
@@ -650,18 +637,17 class PrecipitationProc(Operation):
650 637
651 638 else:
652 639
653 self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() #
654
655 """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
640 self.spc = dataOut.data_pre[0].copy()
656 641
642 #NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX
657 643 self.spc[:,:,0:7]= numpy.NaN
658 644
659 """##########################################"""
660
661 645 self.Num_Hei = self.spc.shape[2]
662 646 self.Num_Bin = self.spc.shape[1]
663 647 self.Num_Chn = self.spc.shape[0]
664 648
649 VelRange = dataOut.spc_range[2]
650
665 651 ''' Se obtiene la constante del RADAR '''
666 652
667 653 self.Pt = Pt
@@ -670,104 +656,72 class PrecipitationProc(Operation):
670 656 self.Lambda = Lambda
671 657 self.aL = aL
672 658 self.tauW = tauW
673 self.ThetaT = ThetaT
659 self.ThetaT = ThetaT
674 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 665 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
677 666 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
678 667 RadarConstant = 10e-26 * Numerator / Denominator #
679
680 ''' ============================= '''
681
682 self.spc[0] = (self.spc[0]-dataOut.noise[0])
683 self.spc[1] = (self.spc[1]-dataOut.noise[1])
684 self.spc[2] = (self.spc[2]-dataOut.noise[2])
685
686 self.spc[ numpy.where(self.spc < 0)] = 0
687
688 SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise))
689 SPCmean[ numpy.where(SPCmean < 0)] = 0
690
691 ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei])
692 ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei])
693 ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei])
694
695 Pr = SPCmean[:,:]
696
697 VelMeteoro = numpy.mean(SPCmean,axis=0)
698
699 D_range = numpy.zeros([self.Num_Bin,self.Num_Hei])
700 SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei])
701 N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei])
702 V_mean = numpy.zeros(self.Num_Hei)
703 del_V = numpy.zeros(self.Num_Hei)
704 Z = numpy.zeros(self.Num_Hei)
705 Ze = numpy.zeros(self.Num_Hei)
706 RR = numpy.zeros(self.Num_Hei)
707
708 Range = dataOut.heightList*1000.
709
710 for R in range(self.Num_Hei):
711
712 h = Range[R] + Altitude #Range from ground to radar pulse altitude
713 del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
714
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
716
717 '''NOTA: ETA(n) dn = ETA(f) df
718
719 dn = 1 Diferencial de muestreo
720 df = ETA(n) / ETA(f)
721
722 '''
723
724 ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA)
725
726 ETAv[:,R]=ETAn[:,R]/dV
727
728 ETAd[:,R]=ETAv[:,R]*6.18*exp(-0.6*D_range[:,R])
729
730 SIGMA[:,R] = Km * (D_range[:,R] * 1e-3 )**6 * numpy.pi**5 / Lambda**4 #Equivalent Section of drops (sigma)
731
732 N_dist[:,R] = ETAn[:,R] / SIGMA[:,R]
733
734 DMoments = self.Moments(Pr[:,R], Velrange[0:self.Num_Bin])
735
736 try:
737 popt01,pcov = curve_fit(self.gaus, Velrange[0:self.Num_Bin] , Pr[:,R] , p0=DMoments)
738 except:
739 popt01=numpy.zeros(3)
740 popt01[1]= DMoments[1]
741
742 if popt01[1]<0 or popt01[1]>20:
743 popt01[1]=numpy.NaN
744
745
746 V_mean[R]=popt01[1]
747
748 Z[R] = numpy.nansum( N_dist[:,R] * (D_range[:,R])**6 )#*10**-18
749
750 RR[R] = 0.0006*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate
751
752 Ze[R] = (numpy.nansum( ETAn[:,R]) * Lambda**4) / ( 10**-18*numpy.pi**5 * Km)
753
754
755
756 RR2 = (Z/200)**(1/1.6)
757 dBRR = 10*numpy.log10(RR)
758 dBRR2 = 10*numpy.log10(RR2)
759
760 dBZe = 10*numpy.log10(Ze)
761 dBZ = 10*numpy.log10(Z)
668 ExpConstant = 10**(40/10) #Constante Experimental
669
670 SignalPower = numpy.zeros([self.Num_Chn,self.Num_Bin,self.Num_Hei])
671 for i in range(self.Num_Chn):
672 SignalPower[i,:,:] = self.spc[i,:,:] - dataOut.noise[i]
673 SignalPower[numpy.where(SignalPower < 0)] = 1e-20
674
675 SPCmean = numpy.mean(SignalPower, 0)
676
677 Pr = SPCmean[:,:]/dataOut.normFactor
678
679 # Declaring auxiliary variables
680 Range = dataOut.heightList*1000. #Range in m
681 # replicate the heightlist to obtain a matrix [Num_Bin,Num_Hei]
682 rMtrx = numpy.transpose(numpy.transpose([dataOut.heightList*1000.] * self.Num_Bin))
683 zMtrx = rMtrx+Altitude
684 # replicate the VelRange to obtain a matrix [Num_Bin,Num_Hei]
685 VelMtrx = numpy.transpose(numpy.tile(VelRange[:-1], (self.Num_Hei,1)))
686
687 # height dependence to air density Foote and Du Toit (1969)
688 delv_z = 1 + 3.68e-5 * zMtrx + 1.71e-9 * zMtrx**2
689 VMtrx = VelMtrx / delv_z #Normalized velocity
690 VMtrx[numpy.where(VMtrx> 9.6)] = numpy.NaN
691 # Diameter is related to the fall speed of falling drops
692 D_Vz = -1.667 * numpy.log( 0.9369 - 0.097087 * VMtrx ) # D in [mm]
693 # Only valid for D>= 0.16 mm
694 D_Vz[numpy.where(D_Vz < 0.16)] = numpy.NaN
695
696 #Calculate Radar Reflectivity ETAn
697 ETAn = (RadarConstant *ExpConstant) * Pr * rMtrx**2 #Reflectivity (ETA)
698 ETAd = ETAn * 6.18 * exp( -0.6 * D_Vz ) * delv_z
699 # Radar Cross Section
700 sigmaD = Km2 * (D_Vz * 1e-3 )**6 * numpy.pi**5 / Lambda**4
701 # Drop Size Distribution
702 DSD = ETAn / sigmaD
703 # Equivalente Reflectivy
704 Ze_eqn = numpy.nansum( DSD * D_Vz**6 ,axis=0)
705 Ze_org = numpy.nansum(ETAn * Lambda**4, axis=0) / (1e-18*numpy.pi**5 * Km2) # [mm^6 /m^3]
706 # RainFall Rate
707 RR = 0.0006*numpy.pi * numpy.nansum( D_Vz**3 * DSD * VelMtrx ,0) #mm/hr
708
709 # Censoring the data
710 # Removing data with SNRth < 0dB se debe considerar el SNR por canal
711 SNRth = 10**(-30/10) #-20dB
712 novalid = numpy.where((dataOut.data_SNR[0,:] <SNRth) | (dataOut.data_SNR[1,:] <SNRth) | (dataOut.data_SNR[2,:] <SNRth)) # AND condition. Maybe OR condition better
713 W = numpy.nanmean(dataOut.data_DOP,0)
714 W[novalid] = numpy.NaN
715 Ze_org[novalid] = numpy.NaN
716 RR[novalid] = numpy.NaN
762 717
763 718 dataOut.data_output = RR[8]
764 719 dataOut.data_param = numpy.ones([3,self.Num_Hei])
765 720 dataOut.channelList = [0,1,2]
766
767 dataOut.data_param[0]=dBZ
768 dataOut.data_param[1]=V_mean
721
722 dataOut.data_param[0]=10*numpy.log10(Ze_org)
723 dataOut.data_param[1]=W
769 724 dataOut.data_param[2]=RR
770
771 725 return dataOut
772 726
773 727 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
@@ -784,7 +738,7 class PrecipitationProc(Operation):
784 738
785 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 743 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
790 744
@@ -832,7 +786,7 class FullSpectralAnalysis(Operation):
832 786
833 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 792 Parameters affected: Winds, height range, SNR
@@ -848,7 +802,7 class FullSpectralAnalysis(Operation):
848 802 """Erick: NOTE THE RANGE OF THE PULSE TX MUST BE REMOVED"""
849 803
850 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 874 velocityY=numpy.append(velocityY, numpy.NaN)
921 875
922 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 878 else:
925 879 velocityV=numpy.append(velocityV, numpy.NaN)
926 880
@@ -944,19 +898,12 class FullSpectralAnalysis(Operation):
944 898 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
945 899
946 900 def Moments(self, ySamples, xSamples):
947 '''***
948 Variables corresponding to moments of distribution.
949 Also used as initial coefficients for curve_fit.
950 Vr was corrected. Only a velocity when x is velocity, of course.
951 ***'''
952 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
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])
901 Power = numpy.nanmean(ySamples) # Power, 0th Moment
902 yNorm = ySamples / numpy.nansum(ySamples)
903 RadVel = numpy.nansum(xSamples * yNorm) # Radial Velocity, 1st Moment
904 Sigma2 = abs(numpy.nansum(yNorm * (xSamples - RadVel)**2)) # Spectral Width, 2nd Moment
905 StdDev = Sigma2**0.5 # Desv. Estandar, Ancho espectral
906 return numpy.array([Power,RadVel,StdDev])
960 907
961 908 def StopWindEstimation(self, error_code):
962 909 '''
@@ -1033,12 +980,13 class FullSpectralAnalysis(Operation):
1033 980 PhaseInter = numpy.ones(spc.shape[0]) # intercept to the slope of the phases, channelwise
1034 981 xFrec = AbbsisaRange[0][0:spc.shape[1]] # frequency range
1035 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 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 988 CSPCmoments = []
1040 989
1041
1042 990 '''Getting Eij and Nij'''
1043 991
1044 992 Xi01, Xi02, Xi12 = ChanDist[:,0]
@@ -1052,10 +1000,13 class FullSpectralAnalysis(Operation):
1052 1000 spc_norm = spc.copy() # need copy() because untouched spc is needed for normalization of cspc below
1053 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 1006 for i in range(spc.shape[0]):
1056 1007
1057 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 1010 Factor_Norm = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc_sub)) # usually = Freq range / nfft
1060 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 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 1067 count_limit_freq = numpy.abs(popt[1]) + widthlimit # Hz, m/s can be also used if velocity is desired abscissa.
1115 1068 # count_limit_freq = numpy.max(xFrec)
1116 1069
@@ -1393,13 +1346,13 class SpectralMoments(Operation):
1393 1346 max_spec = aux.max()
1394 1347 m = aux.tolist().index(max_spec)
1395 1348
1396 #Smooth
1349 # Smooth
1397 1350 if (smooth == 0):
1398 1351 spec2 = spec
1399 1352 else:
1400 1353 spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
1401 1354
1402 # Calculo de Momentos
1355 # Moments Estimation
1403 1356 bb = spec2[numpy.arange(m,spec2.size)]
1404 1357 bb = (bb<n0).nonzero()
1405 1358 bb = bb[0]
@@ -1425,14 +1378,17 class SpectralMoments(Operation):
1425 1378
1426 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 1384 fd = ((spec2[valid]- n0)*freq[valid] * fwindow[valid]).sum() / power
1430 1385 w = numpy.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum() / power)
1431 1386 snr = (spec2.mean()-n0)/n0
1432 1387 if (snr < 1.e-20) :
1433 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 1392 vec_fd[ind] = fd
1437 1393 vec_w[ind] = w
1438 1394 vec_snr[ind] = snr
General Comments 0
You need to be logged in to leave comments. Login now