##// 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 882 velocityX=[]
883 883 velocityY=[]
884 velocityV=[]
885 PhaseLine=[] # unused afterwards
884 velocityV=[]
886 885
887 886 dbSNR = 10*numpy.log10(dataOut.data_SNR)
888 887 dbSNR = numpy.average(dbSNR,0)
@@ -892,12 +891,8 class FullSpectralAnalysis(Operation):
892 891 for Height in range(nHeights):
893 892
894 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)
896
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
894 # error_code unused, yet maybe useful for future analysis.
895 [Vzon,Vmer,Vver, error_code] = self.WindEstimation(spc[:,:,Height], cspc[:,:,Height], pairsList, ChanDist, Height, noise, dataOut.spc_range, dbSNR[Height], SNRlimit)
901 896 else:
902 897 Vzon,Vmer,Vver = 0., 0., numpy.NaN
903 898
@@ -934,8 +929,6 class FullSpectralAnalysis(Operation):
934 929 def gaus(self,xSamples,Amp,Mu,Sigma):
935 930 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
936 931
937
938
939 932 def Moments(self, ySamples, xSamples):
940 933 '''***
941 934 Variables corresponding to moments of distribution.
@@ -950,251 +943,6 class FullSpectralAnalysis(Operation):
950 943 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
951 944
952 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 947 def StopWindEstimation(self, error_code):
1200 948 '''
@@ -1231,9 +979,7 class FullSpectralAnalysis(Operation):
1231 979
1232 980 return antialiased
1233 981
1234
1235
1236 def TestWindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
982 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
1237 983 """
1238 984 Function that Calculates Zonal, Meridional and Vertical wind velocities.
1239 985 Initial Version by E. Bocanegra updated by J. Zibell until Nov. 2019.
@@ -1552,7 +1298,6 class FullSpectralAnalysis(Operation):
1552 1298 return Vzon, Vmer, Vver, error_code
1553 1299
1554 1300
1555
1556 1301 class SpectralMoments(Operation):
1557 1302
1558 1303 '''
General Comments 0
You need to be logged in to leave comments. Login now