|
|
@@
-28,7
+28,6
from scipy.optimize import curve_fit
|
|
|
28
|
import warnings
|
|
28
|
import warnings
|
|
|
29
|
from numpy import NaN
|
|
29
|
from numpy import NaN
|
|
|
30
|
from scipy.optimize.optimize import OptimizeWarning
|
|
30
|
from scipy.optimize.optimize import OptimizeWarning
|
|
|
31
|
from IPython.parallel.controller.scheduler import numpy
|
|
|
|
|
|
32
|
warnings.filterwarnings('ignore')
|
|
31
|
warnings.filterwarnings('ignore')
|
|
|
33
|
|
|
32
|
|
|
|
34
|
|
|
33
|
|
|
|
@@
-213,11
+212,12
class SpectralFilters(Operation):
|
|
|
213
|
Operation.__init__(self, **kwargs)
|
|
212
|
Operation.__init__(self, **kwargs)
|
|
|
214
|
self.i=0
|
|
213
|
self.i=0
|
|
|
215
|
|
|
214
|
|
|
|
216
|
def run(self, dataOut, Rain_Velocity_Limit=1.5, Wind_Velocity_Limit=2.5):
|
|
215
|
def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5):
|
|
|
|
|
|
216
|
|
|
|
217
|
|
|
217
|
|
|
|
218
|
#Limite de vientos
|
|
218
|
#Limite de vientos
|
|
|
219
|
LimitR = Rain_Velocity_Limit
|
|
219
|
LimitR = PositiveLimit
|
|
|
220
|
LimitW = Wind_Velocity_Limit
|
|
220
|
LimitN = NegativeLimit
|
|
|
221
|
|
|
221
|
|
|
|
222
|
self.spc = dataOut.data_pre[0].copy()
|
|
222
|
self.spc = dataOut.data_pre[0].copy()
|
|
|
223
|
self.cspc = dataOut.data_pre[1].copy()
|
|
223
|
self.cspc = dataOut.data_pre[1].copy()
|
|
|
@@
-234,29
+234,31
class SpectralFilters(Operation):
|
|
|
234
|
Tmax= 2*numpy.max(dataOut.spc_range[1])
|
|
234
|
Tmax= 2*numpy.max(dataOut.spc_range[1])
|
|
|
235
|
Fmax= 2*numpy.max(dataOut.spc_range[0])
|
|
235
|
Fmax= 2*numpy.max(dataOut.spc_range[0])
|
|
|
236
|
|
|
236
|
|
|
|
237
|
Breaker1R=VelRange[numpy.abs(VelRange-(-LimitR)).argmin()]
|
|
237
|
Breaker1R=VelRange[numpy.abs(VelRange-(-LimitN)).argmin()]
|
|
|
238
|
Breaker1R=numpy.where(VelRange == Breaker1R)
|
|
238
|
Breaker1R=numpy.where(VelRange == Breaker1R)
|
|
|
239
|
|
|
239
|
|
|
|
240
|
Breaker1W=VelRange[numpy.abs(VelRange-(-LimitW)).argmin()]
|
|
240
|
Delta = self.Num_Bin/2 - Breaker1R[0]
|
|
|
241
|
Breaker1W=numpy.where(VelRange == Breaker1W)
|
|
241
|
|
|
|
|
|
|
242
|
#Breaker1W=VelRange[numpy.abs(VelRange-(-LimitW)).argmin()]
|
|
|
|
|
|
243
|
#Breaker1W=numpy.where(VelRange == Breaker1W)
|
|
|
242
|
|
|
244
|
|
|
|
243
|
Breaker2W=VelRange[numpy.abs(VelRange-(LimitW)).argmin()]
|
|
245
|
#Breaker2W=VelRange[numpy.abs(VelRange-(LimitW)).argmin()]
|
|
|
244
|
Breaker2W=numpy.where(VelRange == Breaker2W)
|
|
246
|
#Breaker2W=numpy.where(VelRange == Breaker2W)
|
|
|
245
|
|
|
247
|
|
|
|
246
|
|
|
248
|
|
|
|
247
|
'''Reacomodando SPCrange'''
|
|
249
|
'''Reacomodando SPCrange'''
|
|
|
248
|
|
|
250
|
|
|
|
249
|
VelRange=numpy.roll(VelRange,-Breaker1R[0],axis=0)
|
|
251
|
VelRange=numpy.roll(VelRange,-(self.Num_Bin/2) ,axis=0)
|
|
|
250
|
|
|
252
|
|
|
|
251
|
VelRange[-int(Breaker1R[0]):]+= Vmax
|
|
253
|
VelRange[-(self.Num_Bin/2):]+= Vmax
|
|
|
252
|
|
|
254
|
|
|
|
253
|
FrecRange=numpy.roll(FrecRange,-Breaker1R[0],axis=0)
|
|
255
|
FrecRange=numpy.roll(FrecRange,-(self.Num_Bin/2),axis=0)
|
|
|
254
|
|
|
256
|
|
|
|
255
|
FrecRange[-int(Breaker1R[0]):]+= Fmax
|
|
257
|
FrecRange[-(self.Num_Bin/2):]+= Fmax
|
|
|
256
|
|
|
258
|
|
|
|
257
|
TimeRange=numpy.roll(TimeRange,-Breaker1R[0],axis=0)
|
|
259
|
TimeRange=numpy.roll(TimeRange,-(self.Num_Bin/2),axis=0)
|
|
|
258
|
|
|
260
|
|
|
|
259
|
TimeRange[-int(Breaker1R[0]):]+= Tmax
|
|
261
|
TimeRange[-(self.Num_Bin/2):]+= Tmax
|
|
|
260
|
|
|
262
|
|
|
|
261
|
''' ------------------ '''
|
|
263
|
''' ------------------ '''
|
|
|
262
|
|
|
264
|
|
|
|
@@
-264,19
+266,21
class SpectralFilters(Operation):
|
|
|
264
|
Breaker2R=numpy.where(VelRange == Breaker2R)
|
|
266
|
Breaker2R=numpy.where(VelRange == Breaker2R)
|
|
|
265
|
|
|
267
|
|
|
|
266
|
|
|
268
|
|
|
|
267
|
|
|
269
|
SPCroll = numpy.roll(self.spc,-(self.Num_Bin/2) ,axis=1)
|
|
|
268
|
|
|
|
|
|
|
269
|
SPCroll = numpy.roll(self.spc,-Breaker1R[0],axis=1)
|
|
|
|
|
|
270
|
|
|
270
|
|
|
|
271
|
SPCcut = SPCroll.copy()
|
|
271
|
SPCcut = SPCroll.copy()
|
|
|
272
|
for i in range(self.Num_Chn):
|
|
272
|
for i in range(self.Num_Chn):
|
|
|
|
|
|
273
|
|
|
|
273
|
SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i]
|
|
274
|
SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i]
|
|
|
|
|
|
275
|
SPCcut[i,-int(Delta):,:] = dataOut.noise[i]
|
|
|
|
|
|
276
|
|
|
|
|
|
|
277
|
#self.spc[i, 0:int(Breaker1W[0]) ,:] = dataOut.noise[i]
|
|
|
|
|
|
278
|
#self.spc[i, int(Breaker2W[0]):self.Num_Bin ,:] = dataOut.noise[i]
|
|
|
|
|
|
279
|
|
|
|
|
|
|
280
|
#self.cspc[i, 0:int(Breaker1W[0]) ,:] = dataOut.noise[i]
|
|
|
|
|
|
281
|
#self.cspc[i, int(Breaker2W[0]):self.Num_Bin ,:] = dataOut.noise[i]
|
|
|
274
|
|
|
282
|
|
|
|
275
|
self.spc[i, 0:int(Breaker1W[0]) ,:] = dataOut.noise[i]
|
|
|
|
|
|
276
|
self.spc[i, int(Breaker2W[0]):self.Num_Bin ,:] = dataOut.noise[i]
|
|
|
|
|
|
277
|
|
|
283
|
|
|
|
278
|
self.cspc[i, 0:int(Breaker1W[0]) ,:] = dataOut.noise[i]
|
|
|
|
|
|
279
|
self.cspc[i, int(Breaker2W[0]):self.Num_Bin ,:] = dataOut.noise[i]
|
|
|
|
|
|
280
|
|
|
284
|
|
|
|
281
|
|
|
285
|
|
|
|
282
|
SPC_ch1 = SPCroll
|
|
286
|
SPC_ch1 = SPCroll
|
|
|
@@
-286,7
+290,7
class SpectralFilters(Operation):
|
|
|
286
|
SPCparam = (SPC_ch1, SPC_ch2, self.spc)
|
|
290
|
SPCparam = (SPC_ch1, SPC_ch2, self.spc)
|
|
|
287
|
dataOut.SPCparam = numpy.asarray(SPCparam)
|
|
291
|
dataOut.SPCparam = numpy.asarray(SPCparam)
|
|
|
288
|
|
|
292
|
|
|
|
289
|
dataOut.data_pre= (self.spc , self.cspc)
|
|
293
|
#dataOut.data_pre= (self.spc , self.cspc)
|
|
|
290
|
|
|
294
|
|
|
|
291
|
#dataOut.data_preParam = (self.spc , self.cspc)
|
|
295
|
#dataOut.data_preParam = (self.spc , self.cspc)
|
|
|
292
|
|
|
296
|
|
|
|
@@
-712,8
+716,8
class PrecipitationProc(Operation):
|
|
|
712
|
tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km = 0.93, Altitude=3350):
|
|
716
|
tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km = 0.93, Altitude=3350):
|
|
|
713
|
|
|
717
|
|
|
|
714
|
|
|
718
|
|
|
|
715
|
Velrange = dataOut.spc_range[2]
|
|
719
|
Velrange = dataOut.spcparam_range[2]
|
|
|
716
|
FrecRange = dataOut.spc_range[0]
|
|
720
|
FrecRange = dataOut.spcparam_range[0]
|
|
|
717
|
|
|
721
|
|
|
|
718
|
dV= Velrange[1]-Velrange[0]
|
|
722
|
dV= Velrange[1]-Velrange[0]
|
|
|
719
|
dF= FrecRange[1]-FrecRange[0]
|
|
723
|
dF= FrecRange[1]-FrecRange[0]
|
|
|
@@
-728,7
+732,7
class PrecipitationProc(Operation):
|
|
|
728
|
|
|
732
|
|
|
|
729
|
else:
|
|
733
|
else:
|
|
|
730
|
|
|
734
|
|
|
|
731
|
self.spc = dataOut.SPCparam[1] #dataOut.data_pre[0].copy() #
|
|
735
|
self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() #
|
|
|
732
|
self.Num_Hei = self.spc.shape[2]
|
|
736
|
self.Num_Hei = self.spc.shape[2]
|
|
|
733
|
self.Num_Bin = self.spc.shape[1]
|
|
737
|
self.Num_Bin = self.spc.shape[1]
|
|
|
734
|
self.Num_Chn = self.spc.shape[0]
|
|
738
|
self.Num_Chn = self.spc.shape[0]
|
|
|
@@
-748,18
+752,26
class PrecipitationProc(Operation):
|
|
|
748
|
|
|
752
|
|
|
|
749
|
Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
|
|
753
|
Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
|
|
|
750
|
Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
|
|
754
|
Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
|
|
|
751
|
RadarConstant = 4.1396e+08# Numerator / Denominator
|
|
755
|
RadarConstant = 1/4.1396e+08# Numerator / Denominator #
|
|
|
752
|
print '***'
|
|
756
|
print '***'
|
|
|
753
|
print '*** RadarConstant' , RadarConstant, '****'
|
|
757
|
print '*** RadarConstant' , RadarConstant, '****'
|
|
|
754
|
print '***'
|
|
758
|
print '***'
|
|
|
755
|
''' ============================= '''
|
|
759
|
''' ============================= '''
|
|
|
756
|
|
|
760
|
|
|
|
757
|
SPCmean = numpy.mean(self.spc,0)
|
|
761
|
self.spc[0] = self.spc[0]-dataOut.noise[0]
|
|
|
758
|
ETAf = numpy.zeros([self.Num_Bin,self.Num_Hei])
|
|
762
|
self.spc[1] = self.spc[1]-dataOut.noise[1]
|
|
|
|
|
|
763
|
self.spc[2] = self.spc[2]-dataOut.noise[2]
|
|
|
|
|
|
764
|
|
|
|
|
|
|
765
|
self.spc[ numpy.where(self.spc < 0)] = 0
|
|
|
|
|
|
766
|
|
|
|
|
|
|
767
|
SPCmean = numpy.mean(self.spc,0) - numpy.mean(dataOut.noise)
|
|
|
|
|
|
768
|
SPCmean[ numpy.where(SPCmean < 0)] = 1e-20
|
|
|
|
|
|
769
|
|
|
|
|
|
|
770
|
ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei])
|
|
|
759
|
ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei])
|
|
771
|
ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei])
|
|
|
760
|
ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei])
|
|
772
|
ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei])
|
|
|
761
|
|
|
773
|
|
|
|
762
|
Pr = self.spc[0,:,:]
|
|
774
|
Pr = SPCmean[:,:]
|
|
|
763
|
|
|
775
|
|
|
|
764
|
VelMeteoro = numpy.mean(SPCmean,axis=0)
|
|
776
|
VelMeteoro = numpy.mean(SPCmean,axis=0)
|
|
|
765
|
|
|
777
|
|
|
|
@@
-774,23
+786,31
class PrecipitationProc(Operation):
|
|
|
774
|
Ze = numpy.zeros(self.Num_Hei)
|
|
786
|
Ze = numpy.zeros(self.Num_Hei)
|
|
|
775
|
RR = numpy.zeros(self.Num_Hei)
|
|
787
|
RR = numpy.zeros(self.Num_Hei)
|
|
|
776
|
|
|
788
|
|
|
|
|
|
|
789
|
Range = dataOut.heightList*1000.
|
|
|
777
|
|
|
790
|
|
|
|
778
|
for R in range(self.Num_Hei):
|
|
791
|
for R in range(self.Num_Hei):
|
|
|
779
|
|
|
792
|
|
|
|
780
|
h = R*75 + Altitude #Range from ground to radar pulse altitude
|
|
793
|
h = Range[R] + Altitude #Range from ground to radar pulse altitude
|
|
|
781
|
del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
|
|
794
|
del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
|
|
|
782
|
|
|
795
|
|
|
|
783
|
D_range[:,R] = numpy.log( (9.65 - (Velrange[0:self.Num_Bin] / del_V[R])) / 10.3 ) / -0.6 #Range of Diameter of drops related to velocity
|
|
796
|
D_range[:,R] = numpy.log( (9.65 - (Velrange[0:self.Num_Bin] / del_V[R])) / 10.3 ) / -0.6 #Diameter range [m]x10**-3
|
|
|
784
|
|
|
797
|
|
|
|
785
|
ETAf[:,R] = 1/RadarConstant * Pr[:,R] * (R*0.075)**2 #Reflectivity (ETA)
|
|
798
|
'''NOTA: ETA(n) dn = ETA(f) df
|
|
|
786
|
|
|
799
|
|
|
|
787
|
ETAv[:,R]=ETAf[:,R]*dF/dV
|
|
800
|
dn = 1 Diferencial de muestreo
|
|
|
|
|
|
801
|
df = ETA(n) / ETA(f)
|
|
|
|
|
|
802
|
|
|
|
|
|
|
803
|
'''
|
|
|
|
|
|
804
|
|
|
|
|
|
|
805
|
ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA)
|
|
|
|
|
|
806
|
|
|
|
|
|
|
807
|
ETAv[:,R]=ETAn[:,R]/dV
|
|
|
788
|
|
|
808
|
|
|
|
789
|
ETAd[:,R]=ETAv[:,R]*6.18*exp(-0.6*D_range[:,R])
|
|
809
|
ETAd[:,R]=ETAv[:,R]*6.18*exp(-0.6*D_range[:,R])
|
|
|
790
|
|
|
810
|
|
|
|
791
|
SIGMA[:,R] = numpy.pi**5 / Lambda**4 * Km * D_range[:,R]**6 #Equivalent Section of drops (sigma)
|
|
811
|
SIGMA[:,R] = Km * (D_range[:,R] * 1e-3 )**6 * numpy.pi**5 / Lambda**4 #Equivalent Section of drops (sigma)
|
|
|
792
|
|
|
812
|
|
|
|
793
|
N_dist[:,R] = ETAd[:,R] / SIGMA[:,R]
|
|
813
|
N_dist[:,R] = ETAn[:,R] / SIGMA[:,R]
|
|
|
794
|
|
|
814
|
|
|
|
795
|
DMoments = self.Moments(Pr[:,R], D_range[:,R])
|
|
815
|
DMoments = self.Moments(Pr[:,R], D_range[:,R])
|
|
|
796
|
|
|
816
|
|
|
|
@@
-801,11
+821,11
class PrecipitationProc(Operation):
|
|
|
801
|
popt01[1]= DMoments[1]
|
|
821
|
popt01[1]= DMoments[1]
|
|
|
802
|
D_mean[R]=popt01[1]
|
|
822
|
D_mean[R]=popt01[1]
|
|
|
803
|
|
|
823
|
|
|
|
804
|
Z[R] = numpy.nansum( N_dist[:,R] * D_range[:,R]**6 )
|
|
824
|
Z[R] = numpy.nansum( N_dist[:,R] * (D_range[:,R])**6 )*1e-18
|
|
|
805
|
|
|
825
|
|
|
|
806
|
RR[R] = 6*10**-4.*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate
|
|
826
|
RR[R] = 6*10**-4.*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate
|
|
|
807
|
|
|
827
|
|
|
|
808
|
Ze[R] = (numpy.nansum(ETAd[:,R]) * Lambda**4) / (numpy.pi * Km)
|
|
828
|
Ze[R] = (numpy.nansum( ETAn[:,R]) * Lambda**4) / ( numpy.pi**5 * Km)
|
|
|
809
|
|
|
829
|
|
|
|
810
|
|
|
830
|
|
|
|
811
|
|
|
831
|
|
|
|
@@
-822,12
+842,24
class PrecipitationProc(Operation):
|
|
|
822
|
|
|
842
|
|
|
|
823
|
dataOut.data_param[0]=dBZ
|
|
843
|
dataOut.data_param[0]=dBZ
|
|
|
824
|
dataOut.data_param[1]=dBZe
|
|
844
|
dataOut.data_param[1]=dBZe
|
|
|
825
|
dataOut.data_param[2]=dBRR2
|
|
845
|
dataOut.data_param[2]=RR
|
|
|
826
|
|
|
846
|
|
|
|
827
|
print 'RR SHAPE', dBRR.shape
|
|
847
|
#print 'VELRANGE', Velrange
|
|
|
828
|
#print 'Ze ', dBZe
|
|
848
|
print 'Range', len(Range)
|
|
|
|
|
|
849
|
print 'delv',del_V
|
|
|
|
|
|
850
|
#print 'DRANGE', D_range[:,50]
|
|
|
|
|
|
851
|
print 'NOISE', dataOut.noise[0]
|
|
|
|
|
|
852
|
print 'radarconstant', RadarConstant
|
|
|
|
|
|
853
|
print 'Range', Range
|
|
|
|
|
|
854
|
print 'ETAn SHAPE', ETAn.shape
|
|
|
|
|
|
855
|
print 'ETAn ', numpy.nansum(ETAn, axis=0)
|
|
|
|
|
|
856
|
print 'ETAd ', numpy.nansum(ETAd, axis=0)
|
|
|
|
|
|
857
|
print 'Pr ', numpy.nansum(Pr, axis=0)
|
|
|
|
|
|
858
|
print 'dataOut.SPCparam[1]', numpy.nansum(dataOut.SPCparam[1][0], axis=0)
|
|
|
|
|
|
859
|
print 'Ze ', dBZe
|
|
|
829
|
print 'Z ', dBZ
|
|
860
|
print 'Z ', dBZ
|
|
|
830
|
print 'RR ', dBRR
|
|
861
|
#print 'RR2 ', RR2
|
|
|
|
|
|
862
|
#print 'RR ', RR
|
|
|
831
|
#print 'RR2 ', dBRR2
|
|
863
|
#print 'RR2 ', dBRR2
|
|
|
832
|
#print 'D_mean', D_mean
|
|
864
|
#print 'D_mean', D_mean
|
|
|
833
|
#print 'del_V', del_V
|
|
865
|
#print 'del_V', del_V
|
|
|
@@
-937,9
+969,9
class FullSpectralAnalysis(Operation):
|
|
|
937
|
|
|
969
|
|
|
|
938
|
dataOut.data_SNR = (numpy.mean(spc,axis=1)- noise[0]) / noise[0]
|
|
970
|
dataOut.data_SNR = (numpy.mean(spc,axis=1)- noise[0]) / noise[0]
|
|
|
939
|
|
|
971
|
|
|
|
|
|
|
972
|
dataOut.data_SNR[numpy.where( dataOut.data_SNR <0 )] = 1e-20
|
|
|
940
|
|
|
973
|
|
|
|
941
|
|
|
974
|
|
|
|
942
|
print dataOut.data_SNR.shape
|
|
|
|
|
|
943
|
#FirstMoment = dataOut.moments[0,1,:]#numpy.average(dataOut.data_param[:,1,:],0)
|
|
975
|
#FirstMoment = dataOut.moments[0,1,:]#numpy.average(dataOut.data_param[:,1,:],0)
|
|
|
944
|
#SecondMoment = numpy.average(dataOut.moments[:,2,:],0)
|
|
976
|
#SecondMoment = numpy.average(dataOut.moments[:,2,:],0)
|
|
|
945
|
|
|
977
|
|
|
|
@@
-1023,7
+1055,6
class FullSpectralAnalysis(Operation):
|
|
|
1023
|
def Moments(self, ySamples, xSamples):
|
|
1055
|
def Moments(self, ySamples, xSamples):
|
|
|
1024
|
Pot = numpy.nansum( ySamples ) # Potencia, momento 0
|
|
1056
|
Pot = numpy.nansum( ySamples ) # Potencia, momento 0
|
|
|
1025
|
yNorm = ySamples / Pot
|
|
1057
|
yNorm = ySamples / Pot
|
|
|
1026
|
|
|
|
|
|
|
1027
|
Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
|
|
1058
|
Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
|
|
|
1028
|
Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
|
|
1059
|
Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
|
|
|
1029
|
Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
|
|
1060
|
Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
|
|
|
@@
-1032,9
+1063,7
class FullSpectralAnalysis(Operation):
|
|
|
1032
|
|
|
1063
|
|
|
|
1033
|
def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
|
|
1064
|
def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
|
|
|
1034
|
|
|
1065
|
|
|
|
1035
|
# print ' '
|
|
1066
|
|
|
|
1036
|
# print '######################## Height',Height, (1000 + 75*Height), '##############################'
|
|
|
|
|
|
1037
|
# print ' '
|
|
|
|
|
|
1038
|
|
|
1067
|
|
|
|
1039
|
ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
|
|
1068
|
ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
|
|
|
1040
|
phase=numpy.ones([spc.shape[0],spc.shape[1]])
|
|
1069
|
phase=numpy.ones([spc.shape[0],spc.shape[1]])
|
|
|
@@
-1113,7
+1142,7
class FullSpectralAnalysis(Operation):
|
|
|
1113
|
# print CSPCmoments
|
|
1142
|
# print CSPCmoments
|
|
|
1114
|
# print '#######################'
|
|
1143
|
# print '#######################'
|
|
|
1115
|
|
|
1144
|
|
|
|
1116
|
popt=[1e-10,1e-10,1e-10]
|
|
1145
|
popt=[1e-10,0,1e-10]
|
|
|
1117
|
popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
|
|
1146
|
popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
|
|
|
1118
|
FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0
|
|
1147
|
FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0
|
|
|
1119
|
|
|
1148
|
|
|
|
@@
-1134,7
+1163,7
class FullSpectralAnalysis(Operation):
|
|
|
1134
|
|
|
1163
|
|
|
|
1135
|
|
|
1164
|
|
|
|
1136
|
'''***Fit Gauss CSPC01***'''
|
|
1165
|
'''***Fit Gauss CSPC01***'''
|
|
|
1137
|
if dbSNR > SNRlimit:
|
|
1166
|
if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3 :
|
|
|
1138
|
try:
|
|
1167
|
try:
|
|
|
1139
|
popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
|
|
1168
|
popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
|
|
|
1140
|
popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1])
|
|
1169
|
popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1])
|
|
|
@@
-1160,7
+1189,7
class FullSpectralAnalysis(Operation):
|
|
|
1160
|
|
|
1189
|
|
|
|
1161
|
yMoments = self.Moments(yMean, xSamples)
|
|
1190
|
yMoments = self.Moments(yMean, xSamples)
|
|
|
1162
|
|
|
1191
|
|
|
|
1163
|
if dbSNR > SNRlimit: # and abs(meanGauss/sigma2) > 0.00001:
|
|
1192
|
if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3: # and abs(meanGauss/sigma2) > 0.00001:
|
|
|
1164
|
try:
|
|
1193
|
try:
|
|
|
1165
|
popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=yMoments)
|
|
1194
|
popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=yMoments)
|
|
|
1166
|
FitGauss=self.gaus(xSamples,*popt)
|
|
1195
|
FitGauss=self.gaus(xSamples,*popt)
|
|
|
@@
-1177,7
+1206,7
class FullSpectralAnalysis(Operation):
|
|
|
1177
|
'''****** Getting Fij ******'''
|
|
1206
|
'''****** Getting Fij ******'''
|
|
|
1178
|
Fijcspc = CSPCopt[:,2]/2*3
|
|
1207
|
Fijcspc = CSPCopt[:,2]/2*3
|
|
|
1179
|
|
|
1208
|
|
|
|
1180
|
#GauWidth = (popt[2]/2)*3
|
|
1209
|
|
|
|
1181
|
GaussCenter = popt[1] #xFrec[GCpos]
|
|
1210
|
GaussCenter = popt[1] #xFrec[GCpos]
|
|
|
1182
|
#Punto en Eje X de la Gaussiana donde se encuentra el centro
|
|
1211
|
#Punto en Eje X de la Gaussiana donde se encuentra el centro
|
|
|
1183
|
ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()]
|
|
1212
|
ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()]
|
|
|
@@
-1242,7
+1271,7
class FullSpectralAnalysis(Operation):
|
|
|
1242
|
|
|
1271
|
|
|
|
1243
|
for i in range(spc.shape[0]):
|
|
1272
|
for i in range(spc.shape[0]):
|
|
|
1244
|
|
|
1273
|
|
|
|
1245
|
if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.6:
|
|
1274
|
if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.3:
|
|
|
1246
|
PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
|
|
1275
|
PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
|
|
|
1247
|
|
|
1276
|
|
|
|
1248
|
#print 'Ancho espectral Frecuencias', FrecRange[-1]-FrecRange[0], 'Hz'
|
|
1277
|
#print 'Ancho espectral Frecuencias', FrecRange[-1]-FrecRange[0], 'Hz'
|
|
|
@@
-1310,7
+1339,10
class FullSpectralAnalysis(Operation):
|
|
|
1310
|
Vmer = Vx
|
|
1339
|
Vmer = Vx
|
|
|
1311
|
Vmag=numpy.sqrt(Vzon**2+Vmer**2)
|
|
1340
|
Vmag=numpy.sqrt(Vzon**2+Vmer**2)
|
|
|
1312
|
Vang=numpy.arctan2(Vmer,Vzon)
|
|
1341
|
Vang=numpy.arctan2(Vmer,Vzon)
|
|
|
1313
|
Vver=SPCmoments[1]
|
|
1342
|
if numpy.abs( popt[1] ) < 3.5 and len(FrecRange)>4:
|
|
|
|
|
|
1343
|
Vver=popt[1]
|
|
|
|
|
|
1344
|
else:
|
|
|
|
|
|
1345
|
Vver=numpy.NaN
|
|
|
1314
|
FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12])
|
|
1346
|
FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12])
|
|
|
1315
|
|
|
1347
|
|
|
|
1316
|
|
|
1348
|
|
|
|
@@
-1458,6
+1490,7
class SpectralMoments(Operation):
|
|
|
1458
|
valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
|
|
1490
|
valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
|
|
|
1459
|
power = ( (spec2[valid] - n0) * fwindow[valid] ).sum()
|
|
1491
|
power = ( (spec2[valid] - n0) * fwindow[valid] ).sum()
|
|
|
1460
|
fd = ( (spec2[valid]- n0) * freq[valid] * fwindow[valid] ).sum() / power
|
|
1492
|
fd = ( (spec2[valid]- n0) * freq[valid] * fwindow[valid] ).sum() / power
|
|
|
|
|
|
1493
|
|
|
|
1461
|
w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
|
|
1494
|
w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
|
|
|
1462
|
snr = (spec2.mean()-n0)/n0
|
|
1495
|
snr = (spec2.mean()-n0)/n0
|
|
|
1463
|
|
|
1496
|
|
