##// END OF EJS Templates
schain
RSS
Merge branch 'claire' into v3.0-devel
Juan C. Espinoza -
r1275:b8daf6525fed merge
parent child
Show More
Side by Side Unified
Context file:
r1275:b8daf6525fed Loading diff...:
Show file before | Show file after | Hide comments
@@ -23,6 +23,7 from numpy import NaN
23 from scipy.optimize.optimize import OptimizeWarning
23 from scipy.optimize.optimize import OptimizeWarning
24 warnings.filterwarnings('ignore')
24 warnings.filterwarnings('ignore')
25
25
26 import matplotlib.pyplot as plt
26
27
27 SPEED_OF_LIGHT = 299792458
28 SPEED_OF_LIGHT = 299792458
28
29
@@ -807,10 +808,10 class PrecipitationProc(Operation):
807 class FullSpectralAnalysis(Operation):
808 class FullSpectralAnalysis(Operation):
808
809
809 """
810 """
810 Function that implements Full Spectral Analisys technique.
811 Function that implements Full Spectral Analysis technique.
811
812
812 Input:
813 Input:
813 self.dataOut.data_pre : SelfSpectra and CrossSPectra data
814 self.dataOut.data_pre : SelfSpectra and CrossSpectra data
814 self.dataOut.groupList : Pairlist of channels
815 self.dataOut.groupList : Pairlist of channels
815 self.dataOut.ChanDist : Physical distance between receivers
816 self.dataOut.ChanDist : Physical distance between receivers
816
817
@@ -823,15 +824,15 class FullSpectralAnalysis(Operation):
823 Parameters affected: Winds, height range, SNR
824 Parameters affected: Winds, height range, SNR
824
825
825 """
826 """
826 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRlimit=7):
827 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRlimit=7, minheight=None, maxheight=None):
827
828
828 self.indice=int(numpy.random.rand()*1000)
829 self.indice=int(numpy.random.rand()*1000)
829
830
830 spc = dataOut.data_pre[0].copy()
831 spc = dataOut.data_pre[0].copy()
831 cspc = dataOut.data_pre[1]
832 cspc = dataOut.data_pre[1]
832
833
833 """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
834 """Erick: NOTE THE RANGE OF THE PULSE TX MUST BE REMOVED"""
834
835
835 SNRspc = spc.copy()
836 SNRspc = spc.copy()
836 SNRspc[:,:,0:7]= numpy.NaN
837 SNRspc[:,:,0:7]= numpy.NaN
837
838
@@ -842,6 +843,24 class FullSpectralAnalysis(Operation):
842 nProfiles = spc.shape[1]
843 nProfiles = spc.shape[1]
843 nHeights = spc.shape[2]
844 nHeights = spc.shape[2]
844
845
846 # first_height = 0.75 #km (ref: data header 20170822)
847 # resolution_height = 0.075 #km
848 '''
849 finding height range. check this when radar parameters are changed!
850 '''
851 if maxheight is not None:
852 # range_max = math.ceil((maxheight - first_height) / resolution_height) # theoretical
853 range_max = math.ceil(13.26 * maxheight - 3) # empirical, works better
854 else:
855 range_max = nHeights
856 if minheight is not None:
857 # range_min = int((minheight - first_height) / resolution_height) # theoretical
858 range_min = int(13.26 * minheight - 5) # empirical, works better
859 if range_min < 0:
860 range_min = 0
861 else:
862 range_min = 0
863
845 pairsList = dataOut.groupList
864 pairsList = dataOut.groupList
846 if dataOut.ChanDist is not None :
865 if dataOut.ChanDist is not None :
847 ChanDist = dataOut.ChanDist
866 ChanDist = dataOut.ChanDist
@@ -850,15 +869,7 class FullSpectralAnalysis(Operation):
850
869
851 FrecRange = dataOut.spc_range[0]
870 FrecRange = dataOut.spc_range[0]
852
871
853 ySamples=numpy.ones([nChannel,nProfiles])
854 phase=numpy.ones([nChannel,nProfiles])
855 CSPCSamples=numpy.ones([nChannel,nProfiles],dtype=numpy.complex_)
856 coherence=numpy.ones([nChannel,nProfiles])
857 PhaseSlope=numpy.ones(nChannel)
858 PhaseInter=numpy.ones(nChannel)
859 data_SNR=numpy.zeros([nProfiles])
872 data_SNR=numpy.zeros([nProfiles])
860
861 data = dataOut.data_pre
862 noise = dataOut.noise
873 noise = dataOut.noise
863
874
864 dataOut.data_SNR = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0]
875 dataOut.data_SNR = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0]
@@ -871,48 +882,53 class FullSpectralAnalysis(Operation):
871 velocityX=[]
882 velocityX=[]
872 velocityY=[]
883 velocityY=[]
873 velocityV=[]
884 velocityV=[]
874 PhaseLine=[]
885 PhaseLine=[] # unused afterwards
875
886
876 dbSNR = 10*numpy.log10(dataOut.data_SNR)
887 dbSNR = 10*numpy.log10(dataOut.data_SNR)
877 dbSNR = numpy.average(dbSNR,0)
888 dbSNR = numpy.average(dbSNR,0)
878
889
890 '''***********************************************WIND ESTIMATION**************************************'''
891
879 for Height in range(nHeights):
892 for Height in range(nHeights):
880
893
881 [Vzon,Vmer,Vver, GaussCenter, PhaseSlope, FitGaussCSPC]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, dataOut.spc_range, dbSNR[Height], SNRlimit)
894 if Height >= range_min and Height < range_max:
882 PhaseLine = numpy.append(PhaseLine, PhaseSlope)
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
901 else:
902 Vzon,Vmer,Vver = 0., 0., numpy.NaN
903
883
904
884 if abs(Vzon)<100. and abs(Vzon)> 0.:
905 if abs(Vzon) < 100. and abs(Vzon) > 0. and abs(Vmer) < 100. and abs(Vmer) > 0.:
885 velocityX=numpy.append(velocityX, Vzon)#Vmag
906 velocityX=numpy.append(velocityX, Vzon)
886
907 velocityY=numpy.append(velocityY, -Vmer)
908
887 else:
909 else:
888 velocityX=numpy.append(velocityX, numpy.NaN)
910 velocityX=numpy.append(velocityX, numpy.NaN)
889
890 if abs(Vmer)<100. and abs(Vmer) > 0.:
891 velocityY=numpy.append(velocityY, -Vmer)#Vang
892
893 else:
894 velocityY=numpy.append(velocityY, numpy.NaN)
911 velocityY=numpy.append(velocityY, numpy.NaN)
895
912
896 if dbSNR[Height] > SNRlimit:
913 if dbSNR[Height] > SNRlimit:
897 velocityV=numpy.append(velocityV, -Vver)#FirstMoment[Height])
914 velocityV=numpy.append(velocityV, -Vver) # reason for this minus sign -> convention? (taken from Ericks version)
898 else:
915 else:
899 velocityV=numpy.append(velocityV, numpy.NaN)
916 velocityV=numpy.append(velocityV, numpy.NaN)
900
901
917
902
918
903 '''Nota: Cambiar el signo de numpy.array(velocityX) cuando se intente procesar datos de BLTR'''
919 '''Change the numpy.array (velocityX) sign when trying to process BLTR data (Erick)'''
904 data_output[0] = numpy.array(velocityX) #self.moving_average(numpy.array(velocityX) , N=1)
920 data_output[0] = numpy.array(velocityX)
905 data_output[1] = numpy.array(velocityY) #self.moving_average(numpy.array(velocityY) , N=1)
921 data_output[1] = numpy.array(velocityY)
906 data_output[2] = velocityV#FirstMoment
922 data_output[2] = velocityV
907
923
908 xFrec=FrecRange[0:spc.shape[1]]
909
924
910 dataOut.data_output=data_output
925 dataOut.data_output = data_output
911
926
912 return dataOut
927 return dataOut
913
928
914
929
915 def moving_average(self,x, N=2):
930 def moving_average(self,x, N=2):
931 """ convolution for smoothenig data. note that last N-1 values are convolution with zeroes """
916 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
932 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
917
933
918 def gaus(self,xSamples,Amp,Mu,Sigma):
934 def gaus(self,xSamples,Amp,Mu,Sigma):
@@ -921,11 +937,17 class FullSpectralAnalysis(Operation):
921
937
922
938
923 def Moments(self, ySamples, xSamples):
939 def Moments(self, ySamples, xSamples):
924 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
940 '''***
941 Variables corresponding to moments of distribution.
942 Also used as initial coefficients for curve_fit.
943 Vr was corrected. Only a velocity when x is velocity, of course.
944 ***'''
945 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
925 yNorm = ySamples / Pot
946 yNorm = ySamples / Pot
926 Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
947 x_range = (numpy.max(xSamples)-numpy.min(xSamples))
927 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
948 Vr = numpy.nansum( yNorm * xSamples )*x_range/len(xSamples) # Velocidad radial, mu, corrimiento doppler, primer momento
928 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
949 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
950 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
929
951
930 return numpy.array([Pot, Vr, Desv])
952 return numpy.array([Pot, Vr, Desv])
931
953
@@ -1171,7 +1193,366 class FullSpectralAnalysis(Operation):
1171
1193
1172
1194
1173 return Vzon, Vmer, Vver, GaussCenter, PhaseSlope, FitGaussCSPC
1195 return Vzon, Vmer, Vver, GaussCenter, PhaseSlope, FitGaussCSPC
1196
1197
1198
1199 def StopWindEstimation(self, error_code):
1200 '''
1201 the wind calculation and returns zeros
1202 '''
1203 Vzon = 0
1204 Vmer = 0
1205 Vver = numpy.nan
1206 return Vzon, Vmer, Vver, error_code
1207
1208 def AntiAliasing(self, interval, maxstep):
1209 """
1210 function to prevent errors from aliased values when computing phaseslope
1211 """
1212 antialiased = numpy.zeros(len(interval))*0.0
1213 copyinterval = interval.copy()
1214
1215 antialiased[0] = copyinterval[0]
1216
1217 for i in range(1,len(antialiased)):
1218
1219 step = interval[i] - interval[i-1]
1220
1221 if step > maxstep:
1222 copyinterval -= 2*numpy.pi
1223 antialiased[i] = copyinterval[i]
1224
1225 elif step < maxstep*(-1):
1226 copyinterval += 2*numpy.pi
1227 antialiased[i] = copyinterval[i]
1228
1229 else:
1230 antialiased[i] = copyinterval[i].copy()
1231
1232 return antialiased
1233
1234
1235
1236 def TestWindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
1237 """
1238 Function that Calculates Zonal, Meridional and Vertical wind velocities.
1239 Initial Version by E. Bocanegra updated by J. Zibell until Nov. 2019.
1240
1241 Input:
1242 spc, cspc : self spectra and cross spectra data. In Briggs notation something like S_i*(S_i)_conj, (S_j)_conj respectively.
1243 pairsList : Pairlist of channels
1244 ChanDist : array of xi_ij and eta_ij
1245 Height : height at which data is processed
1246 noise : noise in [channels] format for specific height
1247 Abbsisarange : range of the frequencies or velocities
1248 dbSNR, SNRlimit : signal to noise ratio in db, lower limit
1249
1250 Output:
1251 Vzon, Vmer, Vver : wind velocities
1252 error_code : int that states where code is terminated
1253
1254 0 : no error detected
1255 1 : Gaussian of mean spc exceeds widthlimit
1256 2 : no Gaussian of mean spc found
1257 3 : SNR to low or velocity to high -> prec. e.g.
1258 4 : at least one Gaussian of cspc exceeds widthlimit
1259 5 : zero out of three cspc Gaussian fits converged
1260 6 : phase slope fit could not be found
1261 7 : arrays used to fit phase have different length
1262 8 : frequency range is either too short (len <= 5) or very long (> 30% of cspc)
1263
1264 """
1265
1266 error_code = 0
1267
1268
1269 SPC_Samples = numpy.ones([spc.shape[0],spc.shape[1]]) # for normalized spc values for one height
1270 phase = numpy.ones([spc.shape[0],spc.shape[1]]) # phase between channels
1271 CSPC_Samples = numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_) # for normalized cspc values
1272 PhaseSlope = numpy.zeros(spc.shape[0]) # slope of the phases, channelwise
1273 PhaseInter = numpy.ones(spc.shape[0]) # intercept to the slope of the phases, channelwise
1274 xFrec = AbbsisaRange[0][0:spc.shape[1]] # frequency range
1275 xVel = AbbsisaRange[2][0:spc.shape[1]] # velocity range
1276 SPCav = numpy.average(spc, axis=0)-numpy.average(noise) # spc[0]-noise[0]
1277
1278 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)
1279 CSPCmoments = []
1280
1281
1282 '''Getting Eij and Nij'''
1283
1284 Xi01, Xi02, Xi12 = ChanDist[:,0]
1285 Eta01, Eta02, Eta12 = ChanDist[:,1]
1286
1287 # update nov 19
1288 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.
1289
1290 '''************************* SPC is normalized ********************************'''
1291
1292 spc_norm = spc.copy() # need copy() because untouched spc is needed for normalization of cspc below
1293 spc_norm = numpy.where(numpy.isfinite(spc_norm), spc_norm, numpy.NAN)
1294
1295 for i in range(spc.shape[0]):
1296
1297 spc_sub = spc_norm[i,:] - noise[i] # spc not smoothed here or in previous version.
1298
1299 Factor_Norm = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc_sub)) # usually = Freq range / nfft
1300 normalized_spc = spc_sub / (numpy.nansum(numpy.abs(spc_sub)) * Factor_Norm)
1301
1302 xSamples = xFrec # the frequency range is taken
1303 SPC_Samples[i] = normalized_spc # Normalized SPC values are taken
1304
1305 '''********************** FITTING MEAN SPC GAUSSIAN **********************'''
1306
1307 """ the gaussian of the mean: first subtract noise, then normalize. this is legal because
1308 you only fit the curve and don't need the absolute value of height for calculation,
1309 only for estimation of width. for normalization of cross spectra, you need initial,
1310 unnormalized self-spectra With noise.
1311
1312 Technically, you don't even need to normalize the self-spectra, as you only need the
1313 width of the peak. However, it was left this way. Note that the normalization has a flaw:
1314 due to subtraction of the noise, some values are below zero. Raw "spc" values should be
1315 >= 0, as it is the modulus squared of the signals (complex * it's conjugate)
1316 """
1317
1318 SPCMean = numpy.average(SPC_Samples, axis=0)
1319
1320 popt = [1e-10,0,1e-10]
1321 SPCMoments = self.Moments(SPCMean, xSamples)
1322
1323 if dbSNR > SNRlimit and numpy.abs(SPCmoments_vel[1]) < 3:
1324 try:
1325 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.
1326 if popt[2] > widthlimit: # CONDITION
1327 return self.StopWindEstimation(error_code = 1)
1328
1329 FitGauss = self.gaus(xSamples,*popt)
1330
1331 except :#RuntimeError:
1332 return self.StopWindEstimation(error_code = 2)
1333
1334 else:
1335 return self.StopWindEstimation(error_code = 3)
1336
1337
1338
1339 '''***************************** CSPC Normalization *************************
1340 new section:
1341 The Spc spectra are used to normalize the crossspectra. Peaks from precipitation
1342 influence the norm which is not desired. First, a range is identified where the
1343 wind peak is estimated -> sum_wind is sum of those frequencies. Next, the area
1344 around it gets cut off and values replaced by mean determined by the boundary
1345 data -> sum_noise (spc is not normalized here, thats why the noise is important)
1346
1347 The sums are then added and multiplied by range/datapoints, because you need
1348 an integral and not a sum for normalization.
1349
1350 A norm is found according to Briggs 92.
1351 '''
1352
1353 radarWavelength = 0.6741 # meters
1354 count_limit_freq = numpy.abs(popt[1]) + widthlimit # Hz, m/s can be also used if velocity is desired abscissa.
1355 # count_limit_freq = numpy.max(xFrec)
1356
1357 channel_integrals = numpy.zeros(3)
1358
1359 for i in range(spc.shape[0]):
1360 '''
1361 find the point in array corresponding to count_limit frequency.
1362 sum over all frequencies in the range around zero Hz @ math.ceil(N_freq/2)
1363 '''
1364 N_freq = numpy.count_nonzero(~numpy.isnan(spc[i,:]))
1365 count_limit_int = int(math.ceil( count_limit_freq / numpy.max(xFrec) * (N_freq / 2) )) # gives integer point
1366 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.
1367 sum_noise = (numpy.mean(spc[i, :4]) + numpy.mean(spc[i, -6:-2]))/2.0 * (N_freq - 2*count_limit_int)
1368 channel_integrals[i] = (sum_noise + sum_wind) * (2*numpy.max(xFrec) / N_freq)
1369
1370
1371 cross_integrals_peak = numpy.zeros(3)
1372 # cross_integrals_totalrange = numpy.zeros(3)
1373
1374 for i in range(spc.shape[0]):
1375
1376 cspc_norm = cspc[i,:].copy() # cspc not smoothed here or in previous version
1377
1378 chan_index0 = pairsList[i][0]
1379 chan_index1 = pairsList[i][1]
1380
1381 cross_integrals_peak[i] = channel_integrals[chan_index0]*channel_integrals[chan_index1]
1382 normalized_cspc = cspc_norm / numpy.sqrt(cross_integrals_peak[i])
1383 CSPC_Samples[i] = normalized_cspc
1384
1385 ''' Finding cross integrals without subtracting any peaks:'''
1386 # FactorNorm0 = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc[chan_index0,:]))
1387 # FactorNorm1 = 2*numpy.max(xFrec) / numpy.count_nonzero(~numpy.isnan(spc[chan_index1,:]))
1388 # cross_integrals_totalrange[i] = (numpy.nansum(spc[chan_index0,:])) * FactorNorm0 * (numpy.nansum(spc[chan_index1,:])) * FactorNorm1
1389 # normalized_cspc = cspc_norm / numpy.sqrt(cross_integrals_totalrange[i])
1390 # CSPC_Samples[i] = normalized_cspc
1391
1392
1393 phase[i] = numpy.arctan2(CSPC_Samples[i].imag, CSPC_Samples[i].real)
1394
1395
1396 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPC_Samples[0]), xSamples),
1397 self.Moments(numpy.abs(CSPC_Samples[1]), xSamples),
1398 self.Moments(numpy.abs(CSPC_Samples[2]), xSamples)])
1399
1400
1401 '''***Sorting out NaN entries***'''
1402 CSPCMask01 = numpy.abs(CSPC_Samples[0])
1403 CSPCMask02 = numpy.abs(CSPC_Samples[1])
1404 CSPCMask12 = numpy.abs(CSPC_Samples[2])
1405
1406 mask01 = ~numpy.isnan(CSPCMask01)
1407 mask02 = ~numpy.isnan(CSPCMask02)
1408 mask12 = ~numpy.isnan(CSPCMask12)
1409
1410 CSPCMask01 = CSPCMask01[mask01]
1411 CSPCMask02 = CSPCMask02[mask02]
1412 CSPCMask12 = CSPCMask12[mask12]
1413
1174
1414
1415 popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
1416 FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0
1417
1418 '''*******************************FIT GAUSS CSPC************************************'''
1419
1420 try:
1421 popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
1422 if popt01[2] > widthlimit: # CONDITION
1423 return self.StopWindEstimation(error_code = 4)
1424
1425 popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1])
1426 if popt02[2] > widthlimit: # CONDITION
1427 return self.StopWindEstimation(error_code = 4)
1428
1429 popt12,pcov = curve_fit(self.gaus,xSamples[mask12],numpy.abs(CSPCMask12),p0=CSPCmoments[2])
1430 if popt12[2] > widthlimit: # CONDITION
1431 return self.StopWindEstimation(error_code = 4)
1432
1433 FitGauss01 = self.gaus(xSamples, *popt01)
1434 FitGauss02 = self.gaus(xSamples, *popt02)
1435 FitGauss12 = self.gaus(xSamples, *popt12)
1436
1437 except:
1438 return self.StopWindEstimation(error_code = 5)
1439
1440
1441 '''************* Getting Fij ***************'''
1442
1443
1444 #Punto en Eje X de la Gaussiana donde se encuentra el centro -- x-axis point of the gaussian where the center is located
1445 # -> PointGauCenter
1446 GaussCenter = popt[1]
1447 ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()]
1448 PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0]
1449
1450 #Punto e^-1 hubicado en la Gaussiana -- point where e^-1 is located in the gaussian
1451 PeMinus1 = numpy.max(FitGauss) * numpy.exp(-1)
1452 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss"
1453 PointFij = numpy.where(FitGauss==FijClosest)[0][0]
1454
1455 Fij = numpy.abs(xSamples[PointFij] - xSamples[PointGauCenter])
1456
1457 '''********** Taking frequency ranges from mean SPCs **********'''
1458
1459 #GaussCenter = popt[1] #Primer momento 01
1460 GauWidth = popt[2] * 3/2 #Ancho de banda de Gau01 -- Bandwidth of Gau01 TODO why *3/2?
1461 Range = numpy.empty(2)
1462 Range[0] = GaussCenter - GauWidth
1463 Range[1] = GaussCenter + GauWidth
1464 #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)
1465 ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()]
1466 ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()]
1467
1468 PointRangeMin = numpy.where(xSamples==ClosRangeMin)[0][0]
1469 PointRangeMax = numpy.where(xSamples==ClosRangeMax)[0][0]
1470
1471 Range = numpy.array([ PointRangeMin, PointRangeMax ])
1472
1473 FrecRange = xFrec[ Range[0] : Range[1] ]
1474
1475
1476 '''************************** Getting Phase Slope ***************************'''
1477
1478 for i in range(1,3): # Changed to only compute two
1479
1480 if len(FrecRange) > 5 and len(FrecRange) < spc.shape[1] * 0.3:
1481 # PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=1) #used before to smooth phase with N=3
1482 PhaseRange = phase[i,Range[0]:Range[1]].copy()
1483
1484 mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
1485
1486
1487 if len(FrecRange) == len(PhaseRange):
1488
1489 try:
1490 slope, intercept, _, _, _ = stats.linregress(FrecRange[mask], self.AntiAliasing(PhaseRange[mask], 4.5))
1491 PhaseSlope[i] = slope
1492 PhaseInter[i] = intercept
1493
1494 except:
1495 return self.StopWindEstimation(error_code = 6)
1496
1497 else:
1498 return self.StopWindEstimation(error_code = 7)
1499
1500 else:
1501 return self.StopWindEstimation(error_code = 8)
1502
1503
1504
1505 '''*** Constants A-H correspond to the convention as in Briggs and Vincent 1992 ***'''
1506
1507 '''Getting constant C'''
1508 cC=(Fij*numpy.pi)**2
1509
1510 '''****** Getting constants F and G ******'''
1511 MijEijNij = numpy.array([[Xi02,Eta02], [Xi12,Eta12]])
1512 MijResult0 = (-PhaseSlope[1] * cC) / (2*numpy.pi)
1513 MijResult1 = (-PhaseSlope[2] * cC) / (2*numpy.pi)
1514 MijResults = numpy.array([MijResult0,MijResult1])
1515 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1516
1517 '''****** Getting constants A, B and H ******'''
1518 W01 = numpy.nanmax( FitGauss01 )
1519 W02 = numpy.nanmax( FitGauss02 )
1520 W12 = numpy.nanmax( FitGauss12 )
1521
1522 WijResult0 = ((cF * Xi01 + cG * Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi / cC))
1523 WijResult1 = ((cF * Xi02 + cG * Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi / cC))
1524 WijResult2 = ((cF * Xi12 + cG * Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi / cC))
1525
1526 WijResults = numpy.array([WijResult0, WijResult1, WijResult2])
1527
1528 WijEijNij = numpy.array([ [Xi01**2, Eta01**2, 2*Xi01*Eta01] , [Xi02**2, Eta02**2, 2*Xi02*Eta02] , [Xi12**2, Eta12**2, 2*Xi12*Eta12] ])
1529 (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1530
1531 VxVy = numpy.array([[cA,cH],[cH,cB]])
1532 VxVyResults = numpy.array([-cF,-cG])
1533 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1534
1535 Vzon = Vy
1536 Vmer = Vx
1537
1538 # Vmag=numpy.sqrt(Vzon**2+Vmer**2) # unused
1539 # Vang=numpy.arctan2(Vmer,Vzon) # unused
1540
1541
1542 ''' using frequency as abscissa. Due to three channels, the offzenith angle is zero
1543 and Vrad equal to Vver. formula taken from Briggs 92, figure 4.
1544 '''
1545 if numpy.abs( popt[1] ) < 3.5 and len(FrecRange) > 4:
1546 Vver = 0.5 * radarWavelength * popt[1] * 100 # *100 to get cm (/s)
1547 else:
1548 Vver = numpy.NaN
1549
1550 error_code = 0
1551
1552 return Vzon, Vmer, Vver, error_code
1553
1554
1555
1175 class SpectralMoments(Operation):
1556 class SpectralMoments(Operation):
1176
1557
1177 '''
1558 '''
General Comments 0
You need to be logged in to leave comments. Login now