##// END OF EJS Templates
schain
RSS
Clean wind estimation FSA
Juan C. Espinoza -
r1276:fe2f3a60be08
parent child
Show More
Side by Side Unified
Context file:
r1276:fe2f3a60be08 Loading diff...:
Show file before | Show file after | Hide comments
@@ -881,8 +881,7 class FullSpectralAnalysis(Operation):
881
881
882 velocityX=[]
882 velocityX=[]
883 velocityY=[]
883 velocityY=[]
884 velocityV=[]
884 velocityV=[]
885 PhaseLine=[] # unused afterwards
886
885
887 dbSNR = 10*numpy.log10(dataOut.data_SNR)
886 dbSNR = 10*numpy.log10(dataOut.data_SNR)
888 dbSNR = numpy.average(dbSNR,0)
887 dbSNR = numpy.average(dbSNR,0)
@@ -892,12 +891,8 class FullSpectralAnalysis(Operation):
892 for Height in range(nHeights):
891 for Height in range(nHeights):
893
892
894 if Height >= range_min and Height < range_max:
893 if Height >= range_min and Height < range_max:
895 # [Vzon,Vmer,Vver, GaussCenter, PhaseSlope, FitGaussCSPC] = self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, dataOut.spc_range, dbSNR[Height], SNRlimit)
894 # error_code unused, yet maybe useful for future analysis.
896
895 [Vzon,Vmer,Vver, error_code] = self.WindEstimation(spc[:,:,Height], cspc[:,:,Height], pairsList, ChanDist, Height, noise, dataOut.spc_range, dbSNR[Height], SNRlimit)
897 # error_code unused, yet maybe useful for future analysis.
898 # Test
899 [Vzon,Vmer,Vver, error_code] = self.TestWindEstimation(spc[:,:,Height], cspc[:,:,Height], pairsList, ChanDist, Height, noise, dataOut.spc_range, dbSNR[Height], SNRlimit)
900
901 else:
896 else:
902 Vzon,Vmer,Vver = 0., 0., numpy.NaN
897 Vzon,Vmer,Vver = 0., 0., numpy.NaN
903
898
@@ -934,8 +929,6 class FullSpectralAnalysis(Operation):
934 def gaus(self,xSamples,Amp,Mu,Sigma):
929 def gaus(self,xSamples,Amp,Mu,Sigma):
935 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
930 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
936
931
937
938
939 def Moments(self, ySamples, xSamples):
932 def Moments(self, ySamples, xSamples):
940 '''***
933 '''***
941 Variables corresponding to moments of distribution.
934 Variables corresponding to moments of distribution.
@@ -950,251 +943,6 class FullSpectralAnalysis(Operation):
950 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
943 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
951
944
952 return numpy.array([Pot, Vr, Desv])
945 return numpy.array([Pot, Vr, Desv])
953
954 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
955
956
957
958 ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
959 phase=numpy.ones([spc.shape[0],spc.shape[1]])
960 CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_)
961 coherence=numpy.ones([spc.shape[0],spc.shape[1]])
962 PhaseSlope=numpy.zeros(spc.shape[0])
963 PhaseInter=numpy.ones(spc.shape[0])
964 xFrec=AbbsisaRange[0][0:spc.shape[1]]
965 xVel =AbbsisaRange[2][0:spc.shape[1]]
966 Vv=numpy.empty(spc.shape[2])*0
967 SPCav = numpy.average(spc, axis=0)-numpy.average(noise) #spc[0]-noise[0]#
968
969 SPCmoments = self.Moments(SPCav[:,Height], xVel )
970 CSPCmoments = []
971 cspcNoise = numpy.empty(3)
972
973 '''Getting Eij and Nij'''
974
975 Xi01=ChanDist[0][0]
976 Eta01=ChanDist[0][1]
977
978 Xi02=ChanDist[1][0]
979 Eta02=ChanDist[1][1]
980
981 Xi12=ChanDist[2][0]
982 Eta12=ChanDist[2][1]
983
984 z = spc.copy()
985 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
986
987 for i in range(spc.shape[0]):
988
989 '''****** Line of Data SPC ******'''
990 zline=z[i,:,Height].copy() - noise[i] # Se resta ruido
991
992 '''****** SPC is normalized ******'''
993 SmoothSPC =self.moving_average(zline.copy(),N=1) # Se suaviza el ruido
994 FactNorm = SmoothSPC/numpy.nansum(SmoothSPC) # SPC Normalizado y suavizado
995
996 xSamples = xFrec # Se toma el rango de frecuncias
997 ySamples[i] = FactNorm # Se toman los valores de SPC normalizado
998
999 for i in range(spc.shape[0]):
1000
1001 '''****** Line of Data CSPC ******'''
1002 cspcLine = ( cspc[i,:,Height].copy())# - noise[i] ) # no! Se resta el ruido
1003 SmoothCSPC =self.moving_average(cspcLine,N=1) # Se suaviza el ruido
1004 cspcNorm = SmoothCSPC/numpy.nansum(SmoothCSPC) # CSPC normalizado y suavizado
1005
1006 '''****** CSPC is normalized with respect to Briggs and Vincent ******'''
1007 chan_index0 = pairsList[i][0]
1008 chan_index1 = pairsList[i][1]
1009
1010 CSPCFactor= numpy.abs(numpy.nansum(ySamples[chan_index0]))**2 * numpy.abs(numpy.nansum(ySamples[chan_index1]))**2
1011 CSPCNorm = cspcNorm / numpy.sqrt(CSPCFactor)
1012
1013 CSPCSamples[i] = CSPCNorm
1014
1015 coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
1016
1017 #coherence[i]= self.moving_average(coherence[i],N=1)
1018
1019 phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
1020
1021 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPCSamples[0]), xSamples),
1022 self.Moments(numpy.abs(CSPCSamples[1]), xSamples),
1023 self.Moments(numpy.abs(CSPCSamples[2]), xSamples)])
1024
1025
1026 popt=[1e-10,0,1e-10]
1027 popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
1028 FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0
1029
1030 CSPCMask01 = numpy.abs(CSPCSamples[0])
1031 CSPCMask02 = numpy.abs(CSPCSamples[1])
1032 CSPCMask12 = numpy.abs(CSPCSamples[2])
1033
1034 mask01 = ~numpy.isnan(CSPCMask01)
1035 mask02 = ~numpy.isnan(CSPCMask02)
1036 mask12 = ~numpy.isnan(CSPCMask12)
1037
1038 #mask = ~numpy.isnan(CSPCMask01)
1039 CSPCMask01 = CSPCMask01[mask01]
1040 CSPCMask02 = CSPCMask02[mask02]
1041 CSPCMask12 = CSPCMask12[mask12]
1042 #CSPCMask01 = numpy.ma.masked_invalid(CSPCMask01)
1043
1044
1045
1046 '''***Fit Gauss CSPC01***'''
1047 if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3 :
1048 try:
1049 popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
1050 popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1])
1051 popt12,pcov = curve_fit(self.gaus,xSamples[mask12],numpy.abs(CSPCMask12),p0=CSPCmoments[2])
1052 FitGauss01 = self.gaus(xSamples,*popt01)
1053 FitGauss02 = self.gaus(xSamples,*popt02)
1054 FitGauss12 = self.gaus(xSamples,*popt12)
1055 except:
1056 FitGauss01=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[0]))
1057 FitGauss02=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[1]))
1058 FitGauss12=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[2]))
1059
1060
1061 CSPCopt = numpy.vstack([popt01,popt02,popt12])
1062
1063 '''****** Getting fij width ******'''
1064
1065 yMean = numpy.average(ySamples, axis=0) # ySamples[0]
1066
1067 '''******* Getting fitting Gaussian *******'''
1068 meanGauss = sum(xSamples*yMean) / len(xSamples) # Mu, velocidad radial (frecuencia)
1069 sigma2 = sum(yMean*(xSamples-meanGauss)**2) / len(xSamples) # Varianza, Ancho espectral (frecuencia)
1070
1071 yMoments = self.Moments(yMean, xSamples)
1072
1073 if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3: # and abs(meanGauss/sigma2) > 0.00001:
1074 try:
1075 popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=yMoments)
1076 FitGauss=self.gaus(xSamples,*popt)
1077
1078 except :#RuntimeError:
1079 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1080
1081
1082 else:
1083 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1084
1085
1086
1087 '''****** Getting Fij ******'''
1088 Fijcspc = CSPCopt[:,2]/2*3
1089
1090
1091 GaussCenter = popt[1] #xFrec[GCpos]
1092 #Punto en Eje X de la Gaussiana donde se encuentra el centro
1093 ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()]
1094 PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0]
1095
1096 #Punto e^-1 hubicado en la Gaussiana
1097 PeMinus1 = numpy.max(FitGauss)* numpy.exp(-1)
1098 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss"
1099 PointFij = numpy.where(FitGauss==FijClosest)[0][0]
1100
1101 if xSamples[PointFij] > xSamples[PointGauCenter]:
1102 Fij = xSamples[PointFij] - xSamples[PointGauCenter]
1103
1104 else:
1105 Fij = xSamples[PointGauCenter] - xSamples[PointFij]
1106
1107
1108 '''****** Taking frequency ranges from SPCs ******'''
1109
1110
1111 #GaussCenter = popt[1] #Primer momento 01
1112 GauWidth = popt[2] *3/2 #Ancho de banda de Gau01
1113 Range = numpy.empty(2)
1114 Range[0] = GaussCenter - GauWidth
1115 Range[1] = GaussCenter + GauWidth
1116 #Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max)
1117 ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()]
1118 ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()]
1119
1120 PointRangeMin = numpy.where(xSamples==ClosRangeMin)[0][0]
1121 PointRangeMax = numpy.where(xSamples==ClosRangeMax)[0][0]
1122
1123 Range=numpy.array([ PointRangeMin, PointRangeMax ])
1124
1125 FrecRange = xFrec[ Range[0] : Range[1] ]
1126 VelRange = xVel[ Range[0] : Range[1] ]
1127
1128
1129 '''****** Getting SCPC Slope ******'''
1130
1131 for i in range(spc.shape[0]):
1132
1133 if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.3:
1134 PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
1135
1136 '''***********************VelRange******************'''
1137
1138 mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
1139
1140 if len(FrecRange) == len(PhaseRange):
1141 try:
1142 slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange[mask], PhaseRange[mask])
1143 PhaseSlope[i]=slope
1144 PhaseInter[i]=intercept
1145 except:
1146 PhaseSlope[i]=0
1147 PhaseInter[i]=0
1148 else:
1149 PhaseSlope[i]=0
1150 PhaseInter[i]=0
1151 else:
1152 PhaseSlope[i]=0
1153 PhaseInter[i]=0
1154
1155
1156 '''Getting constant C'''
1157 cC=(Fij*numpy.pi)**2
1158
1159 '''****** Getting constants F and G ******'''
1160 MijEijNij=numpy.array([[Xi02,Eta02], [Xi12,Eta12]])
1161 MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
1162 MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
1163 MijResults=numpy.array([MijResult0,MijResult1])
1164 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1165
1166 '''****** Getting constants A, B and H ******'''
1167 W01=numpy.nanmax( FitGauss01 ) #numpy.abs(CSPCSamples[0]))
1168 W02=numpy.nanmax( FitGauss02 ) #numpy.abs(CSPCSamples[1]))
1169 W12=numpy.nanmax( FitGauss12 ) #numpy.abs(CSPCSamples[2]))
1170
1171 WijResult0=((cF*Xi01+cG*Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
1172 WijResult1=((cF*Xi02+cG*Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
1173 WijResult2=((cF*Xi12+cG*Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
1174
1175 WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
1176
1177 WijEijNij=numpy.array([ [Xi01**2, Eta01**2, 2*Xi01*Eta01] , [Xi02**2, Eta02**2, 2*Xi02*Eta02] , [Xi12**2, Eta12**2, 2*Xi12*Eta12] ])
1178 (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1179
1180 VxVy=numpy.array([[cA,cH],[cH,cB]])
1181 VxVyResults=numpy.array([-cF,-cG])
1182 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1183
1184 Vzon = Vy
1185 Vmer = Vx
1186 Vmag=numpy.sqrt(Vzon**2+Vmer**2)
1187 Vang=numpy.arctan2(Vmer,Vzon)
1188 if numpy.abs( popt[1] ) < 3.5 and len(FrecRange)>4:
1189 Vver=popt[1]
1190 else:
1191 Vver=numpy.NaN
1192 FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12])
1193
1194
1195 return Vzon, Vmer, Vver, GaussCenter, PhaseSlope, FitGaussCSPC
1196
1197
1198
946
1199 def StopWindEstimation(self, error_code):
947 def StopWindEstimation(self, error_code):
1200 '''
948 '''
@@ -1231,9 +979,7 class FullSpectralAnalysis(Operation):
1231
979
1232 return antialiased
980 return antialiased
1233
981
1234
982 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
1235
1236 def TestWindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
1237 """
983 """
1238 Function that Calculates Zonal, Meridional and Vertical wind velocities.
984 Function that Calculates Zonal, Meridional and Vertical wind velocities.
1239 Initial Version by E. Bocanegra updated by J. Zibell until Nov. 2019.
985 Initial Version by E. Bocanegra updated by J. Zibell until Nov. 2019.
@@ -1552,7 +1298,6 class FullSpectralAnalysis(Operation):
1552 return Vzon, Vmer, Vver, error_code
1298 return Vzon, Vmer, Vver, error_code
1553
1299
1554
1300
1555
1556 class SpectralMoments(Operation):
1301 class SpectralMoments(Operation):
1557
1302
1558 '''
1303 '''
General Comments 0
You need to be logged in to leave comments. Login now