##// END OF EJS Templates
schain-jc
RSS
Merge with Claire_proc
Juan C. Espinoza -
r1188:b1c63639a885 merge
parent child
Show More
Side by Side Unified
Context file:
r1188:b1c63639a885 Loading diff...:
Show file before | Show file after | Hide comments
@@ -112,5 +112,3 schainpy/scripts/
112 .vscode
112 .vscode
113 trash
113 trash
114 *.log
114 *.log
115 schainpy/scripts/testDigitalRF.py
116 schainpy/scripts/testDigitalRFWriter.py
Show file before | Show file after | Hide comments
@@ -190,8 +190,7 class JROData(GenericData):
190 profileIndex = None
190 profileIndex = None
191 error = None
191 error = None
192 data = None
192 data = None
193 data_plt = None
193 nmodes = None
194
195
194
196 def __str__(self):
195 def __str__(self):
197
196
@@ -462,6 +461,7 class Spectra(JROData):
462 ippFactor = None
461 ippFactor = None
463 profileIndex = 0
462 profileIndex = 0
464 plotting = "spectra"
463 plotting = "spectra"
464
465 def __init__(self):
465 def __init__(self):
466 '''
466 '''
467 Constructor
467 Constructor
@@ -554,9 +554,12 class Spectra(JROData):
554
554
555 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
555 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
556 velrange = deltav * (numpy.arange(self.nFFTPoints +
556 velrange = deltav * (numpy.arange(self.nFFTPoints +
557 extrapoints) - self.nFFTPoints / 2.) # - deltav/2
557 extrapoints) - self.nFFTPoints / 2.)
558
558
559 return velrange
559 if self.nmodes:
560 return velrange/self.nmodes
561 else:
562 return velrange
560
563
561 def getNPairs(self):
564 def getNPairs(self):
562
565
Show file before | Show file after | Hide comments
@@ -7,14 +7,190 from .plotting_codes import *
7 from schainpy.model.proc.jroproc_base import MPDecorator
7 from schainpy.model.proc.jroproc_base import MPDecorator
8 from schainpy.utils import log
8 from schainpy.utils import log
9
9
10 class FitGauPlot_(Figure):
10 class ParamLine_(Figure):
11
12 isConfig = None
13
14 def __init__(self):
15
16 self.isConfig = False
17 self.WIDTH = 300
18 self.HEIGHT = 200
19 self.counter_imagwr = 0
20
21 def getSubplots(self):
22
23 nrow = self.nplots
24 ncol = 3
25 return nrow, ncol
26
27 def setup(self, id, nplots, wintitle, show):
28
29 self.nplots = nplots
30
31 self.createFigure(id=id,
32 wintitle=wintitle,
33 show=show)
34
35 nrow,ncol = self.getSubplots()
36 colspan = 3
37 rowspan = 1
38
39 for i in range(nplots):
40 self.addAxes(nrow, ncol, i, 0, colspan, rowspan)
41
42 def plot_iq(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax):
43 yreal = y[channelIndexList,:].real
44 yimag = y[channelIndexList,:].imag
45
46 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
47 xlabel = "Range (Km)"
48 ylabel = "Intensity - IQ"
49
50 if not self.isConfig:
51 nplots = len(channelIndexList)
52
53 self.setup(id=id,
54 nplots=nplots,
55 wintitle='',
56 show=show)
57
58 if xmin == None: xmin = numpy.nanmin(x)
59 if xmax == None: xmax = numpy.nanmax(x)
60 if ymin == None: ymin = min(numpy.nanmin(yreal),numpy.nanmin(yimag))
61 if ymax == None: ymax = max(numpy.nanmax(yreal),numpy.nanmax(yimag))
62
63 self.isConfig = True
64
65 self.setWinTitle(title)
66
67 for i in range(len(self.axesList)):
68 title = "Channel %d" %(i)
69 axes = self.axesList[i]
70
71 axes.pline(x, yreal[i,:],
72 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
73 xlabel=xlabel, ylabel=ylabel, title=title)
74
75 axes.addpline(x, yimag[i,:], idline=1, color="red", linestyle="solid", lw=2)
76
77 def plot_power(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax):
78 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
79 yreal = y.real
80
81 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
82 xlabel = "Range (Km)"
83 ylabel = "Intensity"
84
85 if not self.isConfig:
86 nplots = len(channelIndexList)
87
88 self.setup(id=id,
89 nplots=nplots,
90 wintitle='',
91 show=show)
92
93 if xmin == None: xmin = numpy.nanmin(x)
94 if xmax == None: xmax = numpy.nanmax(x)
95 if ymin == None: ymin = numpy.nanmin(yreal)
96 if ymax == None: ymax = numpy.nanmax(yreal)
97
98 self.isConfig = True
99
100 self.setWinTitle(title)
101
102 for i in range(len(self.axesList)):
103 title = "Channel %d" %(i)
104 axes = self.axesList[i]
105 ychannel = yreal[i,:]
106 axes.pline(x, ychannel,
107 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
108 xlabel=xlabel, ylabel=ylabel, title=title)
109
110
111 def run(self, dataOut, id, wintitle="", channelList=None,
112 xmin=None, xmax=None, ymin=None, ymax=None, save=False,
113 figpath='./', figfile=None, show=True, wr_period=1,
114 ftp=False, server=None, folder=None, username=None, password=None):
115
116 """
117
118 Input:
119 dataOut :
120 id :
121 wintitle :
122 channelList :
123 xmin : None,
124 xmax : None,
125 ymin : None,
126 ymax : None,
127 """
128
129 if channelList == None:
130 channelIndexList = dataOut.channelIndexList
131 else:
132 channelIndexList = []
133 for channel in channelList:
134 if channel not in dataOut.channelList:
135 raise ValueError("Channel %d is not in dataOut.channelList" % channel)
136 channelIndexList.append(dataOut.channelList.index(channel))
137
138 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
139
140 y = dataOut.RR
141
142 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
143 xlabel = "Range (Km)"
144 ylabel = "Intensity"
145
146 if not self.isConfig:
147 nplots = len(channelIndexList)
148
149 self.setup(id=id,
150 nplots=nplots,
151 wintitle='',
152 show=show)
153
154 if xmin == None: xmin = numpy.nanmin(x)
155 if xmax == None: xmax = numpy.nanmax(x)
156 if ymin == None: ymin = numpy.nanmin(y)
157 if ymax == None: ymax = numpy.nanmax(y)
158
159 self.isConfig = True
160
161 self.setWinTitle(title)
162
163 for i in range(len(self.axesList)):
164 title = "Channel %d" %(i)
165 axes = self.axesList[i]
166 ychannel = y[i,:]
167 axes.pline(x, ychannel,
168 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
169 xlabel=xlabel, ylabel=ylabel, title=title)
170
171
172 self.draw()
173
174 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + "_" + str(dataOut.profileIndex)
175 figfile = self.getFilename(name = str_datetime)
176
177 self.save(figpath=figpath,
178 figfile=figfile,
179 save=save,
180 ftp=ftp,
181 wr_period=wr_period,
182 thisDatetime=thisDatetime)
183
184
185
186 class SpcParamPlot_(Figure):
11
187
12 isConfig = None
188 isConfig = None
13 __nsubplots = None
189 __nsubplots = None
14
190
15 WIDTHPROF = None
191 WIDTHPROF = None
16 HEIGHTPROF = None
192 HEIGHTPROF = None
17 PREFIX = 'fitgau'
193 PREFIX = 'SpcParam'
18
194
19 def __init__(self, **kwargs):
195 def __init__(self, **kwargs):
20 Figure.__init__(self, **kwargs)
196 Figure.__init__(self, **kwargs)
@@ -83,7 +259,7 class FitGauPlot_(Figure):
83 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
259 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
84 server=None, folder=None, username=None, password=None,
260 server=None, folder=None, username=None, password=None,
85 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
261 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False,
86 xaxis="frequency", colormap='jet', normFactor=None , GauSelector = 1):
262 xaxis="frequency", colormap='jet', normFactor=None , Selector = 0):
87
263
88 """
264 """
89
265
@@ -119,23 +295,22 class FitGauPlot_(Figure):
119 # else:
295 # else:
120 # factor = normFactor
296 # factor = normFactor
121 if xaxis == "frequency":
297 if xaxis == "frequency":
122 x = dataOut.spc_range[0]
298 x = dataOut.spcparam_range[0]
123 xlabel = "Frequency (kHz)"
299 xlabel = "Frequency (kHz)"
124
300
125 elif xaxis == "time":
301 elif xaxis == "time":
126 x = dataOut.spc_range[1]
302 x = dataOut.spcparam_range[1]
127 xlabel = "Time (ms)"
303 xlabel = "Time (ms)"
128
304
129 else:
305 else:
130 x = dataOut.spc_range[2]
306 x = dataOut.spcparam_range[2]
131 xlabel = "Velocity (m/s)"
307 xlabel = "Velocity (m/s)"
132
308
133 ylabel = "Range (Km)"
309 ylabel = "Range (km)"
134
310
135 y = dataOut.getHeiRange()
311 y = dataOut.getHeiRange()
136
312
137 z = dataOut.GauSPC[:,GauSelector,:,:] #GauSelector] #dataOut.data_spc/factor
313 z = dataOut.SPCparam[Selector] /1966080.0#/ dataOut.normFactor#GauSelector] #dataOut.data_spc/factor
138 print('GausSPC', z[0,32,10:40])
139 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
314 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
140 zdB = 10*numpy.log10(z)
315 zdB = 10*numpy.log10(z)
141
316
@@ -657,7 +832,7 class WindProfilerPlot_(Figure):
657 # tmax = None
832 # tmax = None
658
833
659 x = dataOut.getTimeRange1(dataOut.paramInterval)
834 x = dataOut.getTimeRange1(dataOut.paramInterval)
660 y = dataOut.heightList
835 y = dataOut.heightList
661 z = dataOut.data_output.copy()
836 z = dataOut.data_output.copy()
662 nplots = z.shape[0] #Number of wind dimensions estimated
837 nplots = z.shape[0] #Number of wind dimensions estimated
663 nplotsw = nplots
838 nplotsw = nplots
@@ -666,13 +841,14 class WindProfilerPlot_(Figure):
666 #If there is a SNR function defined
841 #If there is a SNR function defined
667 if dataOut.data_SNR is not None:
842 if dataOut.data_SNR is not None:
668 nplots += 1
843 nplots += 1
669 SNR = dataOut.data_SNR
844 SNR = dataOut.data_SNR[0]
670 SNRavg = numpy.average(SNR, axis=0)
845 SNRavg = SNR#numpy.average(SNR, axis=0)
671
846
672 SNRdB = 10*numpy.log10(SNR)
847 SNRdB = 10*numpy.log10(SNR)
673 SNRavgdB = 10*numpy.log10(SNRavg)
848 SNRavgdB = 10*numpy.log10(SNRavg)
674
849
675 if SNRthresh == None: SNRthresh = -5.0
850 if SNRthresh == None:
851 SNRthresh = -5.0
676 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
852 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
677
853
678 for i in range(nplotsw):
854 for i in range(nplotsw):
@@ -741,8 +917,7 class WindProfilerPlot_(Figure):
741 axes = self.axesList[i*self.__nsubplots]
917 axes = self.axesList[i*self.__nsubplots]
742
918
743 z1 = z[i,:].reshape((1,-1))*windFactor[i]
919 z1 = z[i,:].reshape((1,-1))*windFactor[i]
744 #z1=numpy.ma.masked_where(z1==0.,z1)
920
745
746 axes.pcolorbuffer(x, y, z1,
921 axes.pcolorbuffer(x, y, z1,
747 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
922 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
748 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
923 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
@@ -792,8 +967,8 class ParametersPlot_(Figure):
792 self.isConfig = False
967 self.isConfig = False
793 self.__nsubplots = 1
968 self.__nsubplots = 1
794
969
795 self.WIDTH = 800
970 self.WIDTH = 300
796 self.HEIGHT = 180
971 self.HEIGHT = 550
797 self.WIDTHPROF = 120
972 self.WIDTHPROF = 120
798 self.HEIGHTPROF = 0
973 self.HEIGHTPROF = 0
799 self.counter_imagwr = 0
974 self.counter_imagwr = 0
@@ -905,7 +1080,7 class ParametersPlot_(Figure):
905 # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
1080 # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
906 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
1081 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
907 xlabel = ""
1082 xlabel = ""
908 ylabel = "Range (Km)"
1083 ylabel = "Range (km)"
909
1084
910 update_figfile = False
1085 update_figfile = False
911
1086
@@ -949,24 +1124,81 class ParametersPlot_(Figure):
949
1124
950 self.setWinTitle(title)
1125 self.setWinTitle(title)
951
1126
952 for i in range(self.nchan):
1127 # for i in range(self.nchan):
953 index = channelIndexList[i]
1128 # index = channelIndexList[i]
954 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1129 # title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
955 axes = self.axesList[i*self.plotFact]
1130 # axes = self.axesList[i*self.plotFact]
956 z1 = z[i,:].reshape((1,-1))
1131 # z1 = z[i,:].reshape((1,-1))
957 axes.pcolorbuffer(x, y, z1,
1132 # axes.pcolorbuffer(x, y, z1,
958 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
1133 # xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
959 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1134 # xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
960 ticksize=9, cblabel='', cbsize="1%",colormap=colormap)
1135 # ticksize=9, cblabel='', cbsize="1%",colormap=colormap)
961
1136 #
962 if showSNR:
1137 # if showSNR:
963 title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1138 # title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
964 axes = self.axesList[i*self.plotFact + 1]
1139 # axes = self.axesList[i*self.plotFact + 1]
965 SNRdB1 = SNRdB[i,:].reshape((1,-1))
1140 # SNRdB1 = SNRdB[i,:].reshape((1,-1))
966 axes.pcolorbuffer(x, y, SNRdB1,
1141 # axes.pcolorbuffer(x, y, SNRdB1,
967 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1142 # xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
968 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1143 # xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
969 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
1144 # ticksize=9, cblabel='', cbsize="1%",colormap='jet')
1145
1146 i=0
1147 index = channelIndexList[i]
1148 title = "Factor de reflectividad Z [dBZ]"
1149 axes = self.axesList[i*self.plotFact]
1150 z1 = z[i,:].reshape((1,-1))
1151 axes.pcolorbuffer(x, y, z1,
1152 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
1153 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1154 ticksize=9, cblabel='', cbsize="1%",colormap=colormap)
1155
1156 if showSNR:
1157 title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1158 axes = self.axesList[i*self.plotFact + 1]
1159 SNRdB1 = SNRdB[i,:].reshape((1,-1))
1160 axes.pcolorbuffer(x, y, SNRdB1,
1161 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1162 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1163 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
1164
1165 i=1
1166 index = channelIndexList[i]
1167 title = "Velocidad vertical Doppler [m/s]"
1168 axes = self.axesList[i*self.plotFact]
1169 z1 = z[i,:].reshape((1,-1))
1170 axes.pcolorbuffer(x, y, z1,
1171 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=-10, zmax=10,
1172 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1173 ticksize=9, cblabel='', cbsize="1%",colormap='seismic_r')
1174
1175 if showSNR:
1176 title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1177 axes = self.axesList[i*self.plotFact + 1]
1178 SNRdB1 = SNRdB[i,:].reshape((1,-1))
1179 axes.pcolorbuffer(x, y, SNRdB1,
1180 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1181 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1182 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
1183
1184 i=2
1185 index = channelIndexList[i]
1186 title = "Intensidad de lluvia [mm/h]"
1187 axes = self.axesList[i*self.plotFact]
1188 z1 = z[i,:].reshape((1,-1))
1189 axes.pcolorbuffer(x, y, z1,
1190 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=0, zmax=40,
1191 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1192 ticksize=9, cblabel='', cbsize="1%",colormap='ocean_r')
1193
1194 if showSNR:
1195 title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1196 axes = self.axesList[i*self.plotFact + 1]
1197 SNRdB1 = SNRdB[i,:].reshape((1,-1))
1198 axes.pcolorbuffer(x, y, SNRdB1,
1199 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1200 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1201 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
970
1202
971
1203
972 self.draw()
1204 self.draw()
@@ -1067,9 +1299,8 class Parameters1Plot_(Figure):
1067 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1299 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1068 server=None, folder=None, username=None, password=None,
1300 server=None, folder=None, username=None, password=None,
1069 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1301 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1070 #print inspect.getargspec(self.run).args
1071 """
1072
1302
1303 """
1073 Input:
1304 Input:
1074 dataOut :
1305 dataOut :
1075 id :
1306 id :
Show file before | Show file after | Hide comments
@@ -42,6 +42,8 class SpectraPlot_(Figure):
42
42
43 self.__xfilter_ena = False
43 self.__xfilter_ena = False
44 self.__yfilter_ena = False
44 self.__yfilter_ena = False
45
46 self.indice=1
45
47
46 def getSubplots(self):
48 def getSubplots(self):
47
49
@@ -139,10 +141,9 class SpectraPlot_(Figure):
139 x = dataOut.getVelRange(1)
141 x = dataOut.getVelRange(1)
140 xlabel = "Velocity (m/s)"
142 xlabel = "Velocity (m/s)"
141
143
142 ylabel = "Range (Km)"
144 ylabel = "Range (km)"
143
145
144 y = dataOut.getHeiRange()
146 y = dataOut.getHeiRange()
145
146 z = dataOut.data_spc/factor
147 z = dataOut.data_spc/factor
147 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
148 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
148 zdB = 10*numpy.log10(z)
149 zdB = 10*numpy.log10(z)
@@ -155,6 +156,7 class SpectraPlot_(Figure):
155
156
156 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
157 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
157 title = wintitle + " Spectra"
158 title = wintitle + " Spectra"
159
158 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
160 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
159 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
161 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
160
162
@@ -223,6 +225,7 class SpectraPlot_(Figure):
223 ftp=ftp,
225 ftp=ftp,
224 wr_period=wr_period,
226 wr_period=wr_period,
225 thisDatetime=thisDatetime)
227 thisDatetime=thisDatetime)
228
226
229
227 return dataOut
230 return dataOut
228 @MPDecorator
231 @MPDecorator
@@ -252,6 +255,8 class CrossSpectraPlot_(Figure):
252 self.EXP_CODE = None
255 self.EXP_CODE = None
253 self.SUB_EXP_CODE = None
256 self.SUB_EXP_CODE = None
254 self.PLOT_POS = None
257 self.PLOT_POS = None
258
259 self.indice=0
255
260
256 def getSubplots(self):
261 def getSubplots(self):
257
262
@@ -396,6 +401,7 class CrossSpectraPlot_(Figure):
396 self.isConfig = True
401 self.isConfig = True
397
402
398 self.setWinTitle(title)
403 self.setWinTitle(title)
404
399
405
400 for i in range(self.nplots):
406 for i in range(self.nplots):
401 pair = dataOut.pairsList[pairsIndexList[i]]
407 pair = dataOut.pairsList[pairsIndexList[i]]
@@ -420,7 +426,7 class CrossSpectraPlot_(Figure):
420 xlabel=xlabel, ylabel=ylabel, title=title,
426 xlabel=xlabel, ylabel=ylabel, title=title,
421 ticksize=9, colormap=power_cmap, cblabel='')
427 ticksize=9, colormap=power_cmap, cblabel='')
422
428
423 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:])
429 coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:] / numpy.sqrt( dataOut.data_spc[chan_index0,:,:]*dataOut.data_spc[chan_index1,:,:] )
424 coherence = numpy.abs(coherenceComplex)
430 coherence = numpy.abs(coherenceComplex)
425 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
431 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
426 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
432 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
@@ -439,8 +445,6 class CrossSpectraPlot_(Figure):
439 xlabel=xlabel, ylabel=ylabel, title=title,
445 xlabel=xlabel, ylabel=ylabel, title=title,
440 ticksize=9, colormap=phase_cmap, cblabel='')
446 ticksize=9, colormap=phase_cmap, cblabel='')
441
447
442
443
444 self.draw()
448 self.draw()
445
449
446 self.save(figpath=figpath,
450 self.save(figpath=figpath,
@@ -470,7 +474,7 class RTIPlot_(Figure):
470 self.__nsubplots = 1
474 self.__nsubplots = 1
471
475
472 self.WIDTH = 800
476 self.WIDTH = 800
473 self.HEIGHT = 180
477 self.HEIGHT = 250
474 self.WIDTHPROF = 120
478 self.WIDTHPROF = 120
475 self.HEIGHTPROF = 0
479 self.HEIGHTPROF = 0
476 self.counter_imagwr = 0
480 self.counter_imagwr = 0
@@ -1497,9 +1501,6 class BeaconPhase_(Figure):
1497 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1501 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1498 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1502 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1499
1503
1500 #print "Phase %d%d" %(pair[0], pair[1])
1501 #print phase[dataOut.beacon_heiIndexList]
1502
1503 if dataOut.beacon_heiIndexList:
1504 if dataOut.beacon_heiIndexList:
1504 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1505 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1505 else:
1506 else:
Show file before | Show file after | Hide comments
@@ -434,7 +434,6 def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel=''
434 def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''):
434 def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''):
435
435
436 ax = iplot.axes
436 ax = iplot.axes
437
438 printLabels(ax, xlabel, ylabel, title)
437 printLabels(ax, xlabel, ylabel, title)
439
438
440 for i in range(len(ax.lines)):
439 for i in range(len(ax.lines)):
Show file before | Show file after | Hide comments
@@ -111,7 +111,6 class BLTRParamReader(JRODataReader, ProcessingUnit):
111 timezone=0,
111 timezone=0,
112 status_value=0,
112 status_value=0,
113 **kwargs):
113 **kwargs):
114
115 self.path = path
114 self.path = path
116 self.startDate = startDate
115 self.startDate = startDate
117 self.endDate = endDate
116 self.endDate = endDate
Show file before | Show file after | Hide comments
@@ -1815,7 +1815,7 class JRODataWriter(JRODataIO):
1815
1815
1816 return 1
1816 return 1
1817
1817
1818 def run(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs):
1818 def run(self, dataOut, path, blocksPerFile=100, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs):
1819
1819
1820 if not(self.isConfig):
1820 if not(self.isConfig):
1821
1821
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (607 lines changed) Show them Hide them
@@ -63,9 +63,6 class Header(object):
63 if attr:
63 if attr:
64 message += "%s = %s" % ("size", attr) + "\n"
64 message += "%s = %s" % ("size", attr) + "\n"
65
65
66 # print message
67
68
69 FILE_STRUCTURE = numpy.dtype([ # HEADER 48bytes
66 FILE_STRUCTURE = numpy.dtype([ # HEADER 48bytes
70 ('FileMgcNumber', '<u4'), # 0x23020100
67 ('FileMgcNumber', '<u4'), # 0x23020100
71 # No Of FDT data records in this file (0 or more)
68 # No Of FDT data records in this file (0 or more)
@@ -94,29 +91,6 class FileHeaderBLTR(Header):
94
91
95 header = numpy.fromfile(startFp, FILE_STRUCTURE, 1)
92 header = numpy.fromfile(startFp, FILE_STRUCTURE, 1)
96
93
97 print(' ')
98 print('puntero file header', startFp.tell())
99 print(' ')
100
101 ''' numpy.fromfile(file, dtype, count, sep='')
102 file : file or str
103 Open file object or filename.
104
105 dtype : data-type
106 Data type of the returned array. For binary files, it is used to determine
107 the size and byte-order of the items in the file.
108
109 count : int
110 Number of items to read. -1 means all items (i.e., the complete file).
111
112 sep : str
113 Separator between items if file is a text file. Empty ("") separator means
114 the file should be treated as binary. Spaces (" ") in the separator match zero
115 or more whitespace characters. A separator consisting only of spaces must match
116 at least one whitespace.
117
118 '''
119
120 self.FileMgcNumber = hex(header['FileMgcNumber'][0])
94 self.FileMgcNumber = hex(header['FileMgcNumber'][0])
121 # No Of FDT data records in this file (0 or more)
95 # No Of FDT data records in this file (0 or more)
122 self.nFDTdataRecors = int(header['nFDTdataRecors'][0])
96 self.nFDTdataRecors = int(header['nFDTdataRecors'][0])
@@ -124,8 +98,6 class FileHeaderBLTR(Header):
124 self.OffsetStartHeader = int(header['OffsetStartHeader'][0])
98 self.OffsetStartHeader = int(header['OffsetStartHeader'][0])
125 self.SiteName = str(header['SiteName'][0])
99 self.SiteName = str(header['SiteName'][0])
126
100
127 # print 'Numero de bloques', self.nFDTdataRecors
128
129 if self.size < 48:
101 if self.size < 48:
130 return 0
102 return 0
131
103
@@ -316,36 +288,10 class RecordHeaderBLTR(Header):
316 self.OffsetStartHeader = 48
288 self.OffsetStartHeader = 48
317
289
318 def RHread(self, fp):
290 def RHread(self, fp):
319 # print fp
320 # startFp = open('/home/erick/Documents/Data/huancayo.20161019.22.fdt',"rb") #The method tell() returns the current position of the file read/write pointer within the file.
321 # The method tell() returns the current position of the file read/write pointer within the file.
322 startFp = open(fp, "rb")
291 startFp = open(fp, "rb")
323 # RecCounter=0
324 # Off2StartNxtRec=811248
325 OffRHeader = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
292 OffRHeader = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
326 print(' ')
327 print('puntero Record Header', startFp.tell())
328 print(' ')
329
330 startFp.seek(OffRHeader, os.SEEK_SET)
293 startFp.seek(OffRHeader, os.SEEK_SET)
331
332 print(' ')
333 print('puntero Record Header con seek', startFp.tell())
334 print(' ')
335
336 # print 'Posicion del bloque: ',OffRHeader
337
338 header = numpy.fromfile(startFp, RECORD_STRUCTURE, 1)
294 header = numpy.fromfile(startFp, RECORD_STRUCTURE, 1)
339
340 print(' ')
341 print('puntero Record Header con seek', startFp.tell())
342 print(' ')
343
344 print(' ')
345 #
346 # print 'puntero Record Header despues de seek', header.tell()
347 print(' ')
348
349 self.RecMgcNumber = hex(header['RecMgcNumber'][0]) # 0x23030001
295 self.RecMgcNumber = hex(header['RecMgcNumber'][0]) # 0x23030001
350 self.RecCounter = int(header['RecCounter'][0])
296 self.RecCounter = int(header['RecCounter'][0])
351 self.Off2StartNxtRec = int(header['Off2StartNxtRec'][0])
297 self.Off2StartNxtRec = int(header['Off2StartNxtRec'][0])
@@ -397,52 +343,9 class RecordHeaderBLTR(Header):
397
343
398 self.RHsize = 180 + 20 * self.nChannels
344 self.RHsize = 180 + 20 * self.nChannels
399 self.Datasize = self.nProfiles * self.nChannels * self.nHeights * 2 * 4
345 self.Datasize = self.nProfiles * self.nChannels * self.nHeights * 2 * 4
400 # print 'Datasize',self.Datasize
401 endFp = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
346 endFp = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
402
347
403 print('==============================================')
348
404 print('RecMgcNumber ', self.RecMgcNumber)
405 print('RecCounter ', self.RecCounter)
406 print('Off2StartNxtRec ', self.Off2StartNxtRec)
407 print('Off2StartData ', self.Off2StartData)
408 print('Range Resolution ', self.SampResolution)
409 print('First Height ', self.StartRangeSamp)
410 print('PRF (Hz) ', self.PRFhz)
411 print('Heights (K) ', self.nHeights)
412 print('Channels (N) ', self.nChannels)
413 print('Profiles (J) ', self.nProfiles)
414 print('iCoh ', self.nCohInt)
415 print('iInCoh ', self.nIncohInt)
416 print('BeamAngleAzim ', self.BeamAngleAzim)
417 print('BeamAngleZen ', self.BeamAngleZen)
418
419 # print 'ModoEnUso ',self.DualModeIndex
420 # print 'UtcTime ',self.nUtime
421 # print 'MiliSec ',self.nMilisec
422 # print 'Exp TagName ',self.ExpTagName
423 # print 'Exp Comment ',self.ExpComment
424 # print 'FFT Window Index ',self.FFTwindowingInd
425 # print 'N Dig. Channels ',self.nDigChannels
426 print('Size de bloque ', self.RHsize)
427 print('DataSize ', self.Datasize)
428 print('BeamAngleAzim ', self.BeamAngleAzim)
429 # print 'AntennaCoord0 ',self.AntennaCoord0
430 # print 'AntennaAngl0 ',self.AntennaAngl0
431 # print 'AntennaCoord1 ',self.AntennaCoord1
432 # print 'AntennaAngl1 ',self.AntennaAngl1
433 # print 'AntennaCoord2 ',self.AntennaCoord2
434 # print 'AntennaAngl2 ',self.AntennaAngl2
435 print('RecPhaseCalibr0 ', self.RecPhaseCalibr0)
436 print('RecPhaseCalibr1 ', self.RecPhaseCalibr1)
437 print('RecPhaseCalibr2 ', self.RecPhaseCalibr2)
438 print('RecAmpCalibr0 ', self.RecAmpCalibr0)
439 print('RecAmpCalibr1 ', self.RecAmpCalibr1)
440 print('RecAmpCalibr2 ', self.RecAmpCalibr2)
441 print('ReceiverGaindB0 ', self.ReceiverGaindB0)
442 print('ReceiverGaindB1 ', self.ReceiverGaindB1)
443 print('ReceiverGaindB2 ', self.ReceiverGaindB2)
444 print('==============================================')
445
446 if OffRHeader > endFp:
349 if OffRHeader > endFp:
447 sys.stderr.write(
350 sys.stderr.write(
448 "Warning %s: Size value read from System Header is lower than it has to be\n" % fp)
351 "Warning %s: Size value read from System Header is lower than it has to be\n" % fp)
@@ -537,9 +440,6 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
537 FileList.append(IndexFile)
440 FileList.append(IndexFile)
538 nFiles += 1
441 nFiles += 1
539
442
540 # print 'Files2Read'
541 # print 'Existen '+str(nFiles)+' archivos .fdt'
542
543 self.filenameList = FileList # List of files from least to largest by names
443 self.filenameList = FileList # List of files from least to largest by names
544
444
545 def run(self, **kwargs):
445 def run(self, **kwargs):
@@ -553,7 +453,6 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
553 self.isConfig = True
453 self.isConfig = True
554
454
555 self.getData()
455 self.getData()
556 # print 'running'
557
456
558 def setup(self, path=None,
457 def setup(self, path=None,
559 startDate=None,
458 startDate=None,
@@ -590,22 +489,19 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
590
489
591 if self.flagNoMoreFiles:
490 if self.flagNoMoreFiles:
592 self.dataOut.flagNoData = True
491 self.dataOut.flagNoData = True
593 print('NoData se vuelve true')
594 return 0
492 return 0
595
493
596 self.fp = self.path
494 self.fp = self.path
597 self.Files2Read(self.fp)
495 self.Files2Read(self.fp)
598 self.readFile(self.fp)
496 self.readFile(self.fp)
599 self.dataOut.data_spc = self.data_spc
497 self.dataOut.data_spc = self.data_spc
600 self.dataOut.data_cspc = self.data_cspc
498 self.dataOut.data_cspc =self.data_cspc
601 self.dataOut.data_output = self.data_output
499 self.dataOut.data_output=self.data_output
602
500
603 print('self.dataOut.data_output', shape(self.dataOut.data_output))
501 return self.dataOut.data_spc
604
502
605 # self.removeDC()
503
606 return self.dataOut.data_spc
504 def readFile(self,fp):
607
608 def readFile(self, fp):
609 '''
505 '''
610 You must indicate if you are reading in Online or Offline mode and load the
506 You must indicate if you are reading in Online or Offline mode and load the
611 The parameters for this file reading mode.
507 The parameters for this file reading mode.
@@ -615,23 +511,18 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
615 1. Get the BLTR FileHeader.
511 1. Get the BLTR FileHeader.
616 2. Start reading the first block.
512 2. Start reading the first block.
617 '''
513 '''
618
514
619 # The address of the folder is generated the name of the .fdt file that will be read
620 print("File: ", self.fileSelector + 1)
621
622 if self.fileSelector < len(self.filenameList):
515 if self.fileSelector < len(self.filenameList):
623
516
624 self.fpFile = str(fp) + '/' + \
517 self.fpFile = str(fp) + '/' + \
625 str(self.filenameList[self.fileSelector])
518 str(self.filenameList[self.fileSelector])
626 # print self.fpFile
627 fheader = FileHeaderBLTR()
519 fheader = FileHeaderBLTR()
628 fheader.FHread(self.fpFile) # Bltr FileHeader Reading
520 fheader.FHread(self.fpFile) # Bltr FileHeader Reading
629 self.nFDTdataRecors = fheader.nFDTdataRecors
521 self.nFDTdataRecors = fheader.nFDTdataRecors
630
522
631 self.readBlock() # Block reading
523 self.readBlock() # Block reading
632 else:
524 else:
633 print('readFile FlagNoData becomes true')
525 self.flagNoMoreFiles=True
634 self.flagNoMoreFiles = True
635 self.dataOut.flagNoData = True
526 self.dataOut.flagNoData = True
636 return 0
527 return 0
637
528
@@ -658,12 +549,11 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
658 2. Fill the buffer with the current block number.
549 2. Fill the buffer with the current block number.
659
550
660 '''
551 '''
661
552
662 if self.BlockCounter < self.nFDTdataRecors - 2:
553 if self.BlockCounter < self.nFDTdataRecors-1:
663 print(self.nFDTdataRecors, 'CONDICION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
554 if self.ReadMode==1:
664 if self.ReadMode == 1:
555 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter+1)
665 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter + 1)
556 elif self.ReadMode==0:
666 elif self.ReadMode == 0:
667 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter)
557 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter)
668
558
669 rheader.RHread(self.fpFile) # Bltr FileHeader Reading
559 rheader.RHread(self.fpFile) # Bltr FileHeader Reading
@@ -683,31 +573,26 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
683
573
684 self.nRdPairs = len(self.dataOut.pairsList)
574 self.nRdPairs = len(self.dataOut.pairsList)
685 self.dataOut.nRdPairs = self.nRdPairs
575 self.dataOut.nRdPairs = self.nRdPairs
686
576 self.__firstHeigth=rheader.StartRangeSamp
687 self.__firstHeigth = rheader.StartRangeSamp
577 self.__deltaHeigth=rheader.SampResolution
688 self.__deltaHeigth = rheader.SampResolution
578 self.dataOut.heightList= self.__firstHeigth + numpy.array(range(self.nHeights))*self.__deltaHeigth
689 self.dataOut.heightList = self.__firstHeigth + \
579 self.dataOut.channelList = range(self.nChannels)
690 numpy.array(list(range(self.nHeights))) * self.__deltaHeigth
580 self.dataOut.nProfiles=rheader.nProfiles
691 self.dataOut.channelList = list(range(self.nChannels))
581 self.dataOut.nIncohInt=rheader.nIncohInt
692 self.dataOut.nProfiles = rheader.nProfiles
582 self.dataOut.nCohInt=rheader.nCohInt
693 self.dataOut.nIncohInt = rheader.nIncohInt
583 self.dataOut.ippSeconds= 1/float(rheader.PRFhz)
694 self.dataOut.nCohInt = rheader.nCohInt
584 self.dataOut.PRF=rheader.PRFhz
695 self.dataOut.ippSeconds = 1 / float(rheader.PRFhz)
585 self.dataOut.nFFTPoints=rheader.nProfiles
696 self.dataOut.PRF = rheader.PRFhz
586 self.dataOut.utctime=rheader.nUtime
697 self.dataOut.nFFTPoints = rheader.nProfiles
587 self.dataOut.timeZone=0
698 self.dataOut.utctime = rheader.nUtime
588 self.dataOut.normFactor= self.dataOut.nProfiles*self.dataOut.nIncohInt*self.dataOut.nCohInt
699 self.dataOut.timeZone = 0
589 self.dataOut.outputInterval= self.dataOut.ippSeconds * self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles
700 self.dataOut.normFactor = self.dataOut.nProfiles * \
590
701 self.dataOut.nIncohInt * self.dataOut.nCohInt
591 self.data_output=numpy.ones([3,rheader.nHeights])*numpy.NaN
702 self.dataOut.outputInterval = self.dataOut.ippSeconds * \
592 self.dataOut.velocityX=[]
703 self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles
593 self.dataOut.velocityY=[]
704
594 self.dataOut.velocityV=[]
705 self.data_output = numpy.ones([3, rheader.nHeights]) * numpy.NaN
595
706 print('self.data_output', shape(self.data_output))
707 self.dataOut.velocityX = []
708 self.dataOut.velocityY = []
709 self.dataOut.velocityV = []
710
711 '''Block Reading, the Block Data is received and Reshape is used to give it
596 '''Block Reading, the Block Data is received and Reshape is used to give it
712 shape.
597 shape.
713 '''
598 '''
@@ -734,18 +619,17 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
734 y = rho * numpy.sin(phi)
619 y = rho * numpy.sin(phi)
735 return(x, y)
620 return(x, y)
736
621
737 if self.DualModeIndex == self.ReadMode:
622 if self.DualModeIndex==self.ReadMode:
738
623
739 self.data_fft = numpy.fromfile(
624 self.data_fft = numpy.fromfile( startDATA, [('complex','<c8')],self.nProfiles*self.nChannels*self.nHeights )
740 startDATA, [('complex', '<c8')], self.nProfiles * self.nChannels * self.nHeights)
625 self.data_fft = numpy.empty(101376)
741
626
742 self.data_fft = self.data_fft.astype(numpy.dtype('complex'))
627 self.data_fft=self.data_fft.astype(numpy.dtype('complex'))
743
628
744 self.data_block = numpy.reshape(
629 self.data_block=numpy.reshape(self.data_fft,(self.nHeights, self.nChannels, self.nProfiles ))
745 self.data_fft, (self.nHeights, self.nChannels, self.nProfiles))
630
746
631 self.data_block = numpy.transpose(self.data_block, (1,2,0))
747 self.data_block = numpy.transpose(self.data_block, (1, 2, 0))
632
748
749 copy = self.data_block.copy()
633 copy = self.data_block.copy()
750 spc = copy * numpy.conjugate(copy)
634 spc = copy * numpy.conjugate(copy)
751
635
@@ -756,18 +640,8 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
756
640
757 z = self.data_spc.copy() # /factor
641 z = self.data_spc.copy() # /factor
758 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
642 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
759 #zdB = 10*numpy.log10(z)
643 self.dataOut.data_spc=self.data_spc
760 print(' ')
644 self.noise = self.dataOut.getNoise(ymin_index=80, ymax_index=132)#/factor
761 print('Z: ')
762 print(shape(z))
763 print(' ')
764 print(' ')
765
766 self.dataOut.data_spc = self.data_spc
767
768 self.noise = self.dataOut.getNoise(
769 ymin_index=80, ymax_index=132) # /factor
770 #noisedB = 10*numpy.log10(self.noise)
771
645
772 ySamples = numpy.ones([3, self.nProfiles])
646 ySamples = numpy.ones([3, self.nProfiles])
773 phase = numpy.ones([3, self.nProfiles])
647 phase = numpy.ones([3, self.nProfiles])
@@ -778,20 +652,16 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
778 PhaseInter = numpy.ones(3)
652 PhaseInter = numpy.ones(3)
779
653
780 '''****** Getting CrossSpectra ******'''
654 '''****** Getting CrossSpectra ******'''
781 cspc = self.data_block.copy()
655 cspc=self.data_block.copy()
782 self.data_cspc = self.data_block.copy()
656 self.data_cspc=self.data_block.copy()
783
657
784 xFrec = self.getVelRange(1)
658 xFrec=self.getVelRange(1)
785 VelRange = self.getVelRange(1)
659 VelRange=self.getVelRange(1)
786 self.dataOut.VelRange = VelRange
660 self.dataOut.VelRange=VelRange
787 # print ' '
661
788 # print ' '
662
789 # print 'xFrec',xFrec
663 for i in range(self.nRdPairs):
790 # print ' '
664
791 # print ' '
792 # Height=35
793 for i in range(self.nRdPairs):
794
795 chan_index0 = self.dataOut.pairsList[i][0]
665 chan_index0 = self.dataOut.pairsList[i][0]
796 chan_index1 = self.dataOut.pairsList[i][1]
666 chan_index1 = self.dataOut.pairsList[i][1]
797
667
@@ -820,361 +690,8 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
820
690
821 self.dataOut.ChanDist = self.ChanDist
691 self.dataOut.ChanDist = self.ChanDist
822
692
823
693 self.BlockCounter+=2
824 # for Height in range(self.nHeights):
694
825 #
826 # for i in range(self.nRdPairs):
827 #
828 # '''****** Line of Data SPC ******'''
829 # zline=z[i,:,Height]
830 #
831 # '''****** DC is removed ******'''
832 # DC=Find(zline,numpy.amax(zline))
833 # zline[DC]=(zline[DC-1]+zline[DC+1])/2
834 #
835 #
836 # '''****** SPC is normalized ******'''
837 # FactNorm= zline.copy() / numpy.sum(zline.copy())
838 # FactNorm= FactNorm/numpy.sum(FactNorm)
839 #
840 # SmoothSPC=moving_average(FactNorm,N=3)
841 #
842 # xSamples = ar(range(len(SmoothSPC)))
843 # ySamples[i] = SmoothSPC-self.noise[i]
844 #
845 # for i in range(self.nRdPairs):
846 #
847 # '''****** Line of Data CSPC ******'''
848 # cspcLine=self.data_cspc[i,:,Height].copy()
849 #
850 #
851 #
852 # '''****** CSPC is normalized ******'''
853 # chan_index0 = self.dataOut.pairsList[i][0]
854 # chan_index1 = self.dataOut.pairsList[i][1]
855 # CSPCFactor= numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])
856 #
857 #
858 # CSPCNorm= cspcLine.copy() / numpy.sqrt(CSPCFactor)
859 #
860 #
861 # CSPCSamples[i] = CSPCNorm-self.noise[i]
862 # coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
863 #
864 # '''****** DC is removed ******'''
865 # DC=Find(coherence[i],numpy.amax(coherence[i]))
866 # coherence[i][DC]=(coherence[i][DC-1]+coherence[i][DC+1])/2
867 # coherence[i]= moving_average(coherence[i],N=2)
868 #
869 # phase[i] = moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
870 #
871 #
872 # '''****** Getting fij width ******'''
873 #
874 # yMean=[]
875 # yMean2=[]
876 #
877 # for j in range(len(ySamples[1])):
878 # yMean=numpy.append(yMean,numpy.average([ySamples[0,j],ySamples[1,j],ySamples[2,j]]))
879 #
880 # '''******* Getting fitting Gaussian ******'''
881 # meanGauss=sum(xSamples*yMean) / len(xSamples)
882 # sigma=sum(yMean*(xSamples-meanGauss)**2) / len(xSamples)
883 # #print 'Height',Height,'SNR', meanGauss/sigma**2
884 #
885 # if (abs(meanGauss/sigma**2) > 0.0001) :
886 #
887 # try:
888 # popt,pcov = curve_fit(gaus,xSamples,yMean,p0=[1,meanGauss,sigma])
889 #
890 # if numpy.amax(popt)>numpy.amax(yMean)*0.3:
891 # FitGauss=gaus(xSamples,*popt)
892 #
893 # else:
894 # FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
895 # print 'Verificador: Dentro', Height
896 # except RuntimeError:
897 #
898 # try:
899 # for j in range(len(ySamples[1])):
900 # yMean2=numpy.append(yMean2,numpy.average([ySamples[1,j],ySamples[2,j]]))
901 # popt,pcov = curve_fit(gaus,xSamples,yMean2,p0=[1,meanGauss,sigma])
902 # FitGauss=gaus(xSamples,*popt)
903 # print 'Verificador: Exepcion1', Height
904 # except RuntimeError:
905 #
906 # try:
907 # popt,pcov = curve_fit(gaus,xSamples,ySamples[1],p0=[1,meanGauss,sigma])
908 # FitGauss=gaus(xSamples,*popt)
909 # print 'Verificador: Exepcion2', Height
910 # except RuntimeError:
911 # FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
912 # print 'Verificador: Exepcion3', Height
913 # else:
914 # FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
915 # #print 'Verificador: Fuera', Height
916 #
917 #
918 #
919 # Maximun=numpy.amax(yMean)
920 # eMinus1=Maximun*numpy.exp(-1)
921 #
922 # HWpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1)))
923 # HalfWidth= xFrec[HWpos]
924 # GCpos=Find(FitGauss, numpy.amax(FitGauss))
925 # Vpos=Find(FactNorm, numpy.amax(FactNorm))
926 # #Vpos=numpy.sum(FactNorm)/len(FactNorm)
927 # #Vpos=Find(FactNorm, min(FactNorm, key=lambda value:abs(value- numpy.mean(FactNorm) )))
928 # #print 'GCpos',GCpos, numpy.amax(FitGauss), 'HWpos',HWpos
929 # '''****** Getting Fij ******'''
930 #
931 # GaussCenter=xFrec[GCpos]
932 # if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0):
933 # Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001
934 # else:
935 # Fij=abs(GaussCenter-HalfWidth)+0.0000001
936 #
937 # '''****** Getting Frecuency range of significant data ******'''
938 #
939 # Rangpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10)))
940 #
941 # if Rangpos<GCpos:
942 # Range=numpy.array([Rangpos,2*GCpos-Rangpos])
943 # else:
944 # Range=numpy.array([2*GCpos-Rangpos,Rangpos])
945 #
946 # FrecRange=xFrec[Range[0]:Range[1]]
947 #
948 # #print 'FrecRange', FrecRange
949 # '''****** Getting SCPC Slope ******'''
950 #
951 # for i in range(self.nRdPairs):
952 #
953 # if len(FrecRange)>5 and len(FrecRange)<self.nProfiles*0.5:
954 # PhaseRange=moving_average(phase[i,Range[0]:Range[1]],N=3)
955 #
956 # slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange)
957 # PhaseSlope[i]=slope
958 # PhaseInter[i]=intercept
959 # else:
960 # PhaseSlope[i]=0
961 # PhaseInter[i]=0
962 #
963 # # plt.figure(i+15)
964 # # plt.title('FASE ( CH%s*CH%s )' %(self.dataOut.pairsList[i][0],self.dataOut.pairsList[i][1]))
965 # # plt.xlabel('Frecuencia (KHz)')
966 # # plt.ylabel('Magnitud')
967 # # #plt.subplot(311+i)
968 # # plt.plot(FrecRange,PhaseRange,'b')
969 # # plt.plot(FrecRange,FrecRange*PhaseSlope[i]+PhaseInter[i],'r')
970 #
971 # #plt.axis([-0.6, 0.2, -3.2, 3.2])
972 #
973 #
974 # '''Getting constant C'''
975 # cC=(Fij*numpy.pi)**2
976 #
977 # # '''Getting Eij and Nij'''
978 # # (AntennaX0,AntennaY0)=pol2cart(rheader.AntennaCoord0, rheader.AntennaAngl0*numpy.pi/180)
979 # # (AntennaX1,AntennaY1)=pol2cart(rheader.AntennaCoord1, rheader.AntennaAngl1*numpy.pi/180)
980 # # (AntennaX2,AntennaY2)=pol2cart(rheader.AntennaCoord2, rheader.AntennaAngl2*numpy.pi/180)
981 # #
982 # # E01=AntennaX0-AntennaX1
983 # # N01=AntennaY0-AntennaY1
984 # #
985 # # E02=AntennaX0-AntennaX2
986 # # N02=AntennaY0-AntennaY2
987 # #
988 # # E12=AntennaX1-AntennaX2
989 # # N12=AntennaY1-AntennaY2
990 #
991 # '''****** Getting constants F and G ******'''
992 # MijEijNij=numpy.array([[E02,N02], [E12,N12]])
993 # MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
994 # MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
995 # MijResults=numpy.array([MijResult0,MijResult1])
996 # (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
997 #
998 # '''****** Getting constants A, B and H ******'''
999 # W01=numpy.amax(coherence[0])
1000 # W02=numpy.amax(coherence[1])
1001 # W12=numpy.amax(coherence[2])
1002 #
1003 # WijResult0=((cF*E01+cG*N01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
1004 # WijResult1=((cF*E02+cG*N02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
1005 # WijResult2=((cF*E12+cG*N12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
1006 #
1007 # WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
1008 #
1009 # WijEijNij=numpy.array([ [E01**2, N01**2, 2*E01*N01] , [E02**2, N02**2, 2*E02*N02] , [E12**2, N12**2, 2*E12*N12] ])
1010 # (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1011 #
1012 # VxVy=numpy.array([[cA,cH],[cH,cB]])
1013 #
1014 # VxVyResults=numpy.array([-cF,-cG])
1015 # (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1016 # Vzon = Vy
1017 # Vmer = Vx
1018 # Vmag=numpy.sqrt(Vzon**2+Vmer**2)
1019 # Vang=numpy.arctan2(Vmer,Vzon)
1020 #
1021 # if abs(Vy)<100 and abs(Vy)> 0.:
1022 # self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, Vzon) #Vmag
1023 # #print 'Vmag',Vmag
1024 # else:
1025 # self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, NaN)
1026 #
1027 # if abs(Vx)<100 and abs(Vx) > 0.:
1028 # self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, Vmer) #Vang
1029 # #print 'Vang',Vang
1030 # else:
1031 # self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, NaN)
1032 #
1033 # if abs(GaussCenter)<2:
1034 # self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, xFrec[Vpos])
1035 #
1036 # else:
1037 # self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, NaN)
1038 #
1039 #
1040 # # print '********************************************'
1041 # # print 'HalfWidth ', HalfWidth
1042 # # print 'Maximun ', Maximun
1043 # # print 'eMinus1 ', eMinus1
1044 # # print 'Rangpos ', Rangpos
1045 # # print 'GaussCenter ',GaussCenter
1046 # # print 'E01 ',E01
1047 # # print 'N01 ',N01
1048 # # print 'E02 ',E02
1049 # # print 'N02 ',N02
1050 # # print 'E12 ',E12
1051 # # print 'N12 ',N12
1052 # #print 'self.dataOut.velocityX ', self.dataOut.velocityX
1053 # # print 'Fij ', Fij
1054 # # print 'cC ', cC
1055 # # print 'cF ', cF
1056 # # print 'cG ', cG
1057 # # print 'cA ', cA
1058 # # print 'cB ', cB
1059 # # print 'cH ', cH
1060 # # print 'Vx ', Vx
1061 # # print 'Vy ', Vy
1062 # # print 'Vmag ', Vmag
1063 # # print 'Vang ', Vang*180/numpy.pi
1064 # # print 'PhaseSlope ',PhaseSlope[0]
1065 # # print 'PhaseSlope ',PhaseSlope[1]
1066 # # print 'PhaseSlope ',PhaseSlope[2]
1067 # # print '********************************************'
1068 # #print 'data_output',shape(self.dataOut.velocityX), shape(self.dataOut.velocityY)
1069 #
1070 # #print 'self.dataOut.velocityX', len(self.dataOut.velocityX)
1071 # #print 'self.dataOut.velocityY', len(self.dataOut.velocityY)
1072 # #print 'self.dataOut.velocityV', self.dataOut.velocityV
1073 #
1074 # self.data_output[0]=numpy.array(self.dataOut.velocityX)
1075 # self.data_output[1]=numpy.array(self.dataOut.velocityY)
1076 # self.data_output[2]=numpy.array(self.dataOut.velocityV)
1077 #
1078 # prin= self.data_output[0][~numpy.isnan(self.data_output[0])]
1079 # print ' '
1080 # print 'VmagAverage',numpy.mean(prin)
1081 # print ' '
1082 # # plt.figure(5)
1083 # # plt.subplot(211)
1084 # # plt.plot(self.dataOut.velocityX,'yo:')
1085 # # plt.subplot(212)
1086 # # plt.plot(self.dataOut.velocityY,'yo:')
1087 #
1088 # # plt.figure(1)
1089 # # # plt.subplot(121)
1090 # # # plt.plot(xFrec,ySamples[0],'k',label='Ch0')
1091 # # # plt.plot(xFrec,ySamples[1],'g',label='Ch1')
1092 # # # plt.plot(xFrec,ySamples[2],'r',label='Ch2')
1093 # # # plt.plot(xFrec,FitGauss,'yo:',label='fit')
1094 # # # plt.legend()
1095 # # plt.title('DATOS A ALTURA DE 2850 METROS')
1096 # #
1097 # # plt.xlabel('Frecuencia (KHz)')
1098 # # plt.ylabel('Magnitud')
1099 # # # plt.subplot(122)
1100 # # # plt.title('Fit for Time Constant')
1101 # # #plt.plot(xFrec,zline)
1102 # # #plt.plot(xFrec,SmoothSPC,'g')
1103 # # plt.plot(xFrec,FactNorm)
1104 # # plt.axis([-4, 4, 0, 0.15])
1105 # # # plt.xlabel('SelfSpectra KHz')
1106 # #
1107 # # plt.figure(10)
1108 # # # plt.subplot(121)
1109 # # plt.plot(xFrec,ySamples[0],'b',label='Ch0')
1110 # # plt.plot(xFrec,ySamples[1],'y',label='Ch1')
1111 # # plt.plot(xFrec,ySamples[2],'r',label='Ch2')
1112 # # # plt.plot(xFrec,FitGauss,'yo:',label='fit')
1113 # # plt.legend()
1114 # # plt.title('SELFSPECTRA EN CANALES')
1115 # #
1116 # # plt.xlabel('Frecuencia (KHz)')
1117 # # plt.ylabel('Magnitud')
1118 # # # plt.subplot(122)
1119 # # # plt.title('Fit for Time Constant')
1120 # # #plt.plot(xFrec,zline)
1121 # # #plt.plot(xFrec,SmoothSPC,'g')
1122 # # # plt.plot(xFrec,FactNorm)
1123 # # # plt.axis([-4, 4, 0, 0.15])
1124 # # # plt.xlabel('SelfSpectra KHz')
1125 # #
1126 # # plt.figure(9)
1127 # #
1128 # #
1129 # # plt.title('DATOS SUAVIZADOS')
1130 # # plt.xlabel('Frecuencia (KHz)')
1131 # # plt.ylabel('Magnitud')
1132 # # plt.plot(xFrec,SmoothSPC,'g')
1133 # #
1134 # # #plt.plot(xFrec,FactNorm)
1135 # # plt.axis([-4, 4, 0, 0.15])
1136 # # # plt.xlabel('SelfSpectra KHz')
1137 # # #
1138 # # plt.figure(2)
1139 # # # #plt.subplot(121)
1140 # # plt.plot(xFrec,yMean,'r',label='Mean SelfSpectra')
1141 # # plt.plot(xFrec,FitGauss,'yo:',label='Ajuste Gaussiano')
1142 # # # plt.plot(xFrec[Rangpos],FitGauss[Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.1)))],'bo')
1143 # # # #plt.plot(xFrec,phase)
1144 # # # plt.xlabel('Suavizado, promediado KHz')
1145 # # plt.title('SELFSPECTRA PROMEDIADO')
1146 # # # #plt.subplot(122)
1147 # # # #plt.plot(xSamples,zline)
1148 # # plt.xlabel('Frecuencia (KHz)')
1149 # # plt.ylabel('Magnitud')
1150 # # plt.legend()
1151 # # #
1152 # # # plt.figure(3)
1153 # # # plt.subplot(311)
1154 # # # #plt.plot(xFrec,phase[0])
1155 # # # plt.plot(xFrec,phase[0],'g')
1156 # # # plt.subplot(312)
1157 # # # plt.plot(xFrec,phase[1],'g')
1158 # # # plt.subplot(313)
1159 # # # plt.plot(xFrec,phase[2],'g')
1160 # # # #plt.plot(xFrec,phase[2])
1161 # # #
1162 # # # plt.figure(4)
1163 # # #
1164 # # # plt.plot(xSamples,coherence[0],'b')
1165 # # # plt.plot(xSamples,coherence[1],'r')
1166 # # # plt.plot(xSamples,coherence[2],'g')
1167 # # plt.show()
1168 # # #
1169 # # # plt.clf()
1170 # # # plt.cla()
1171 # # # plt.close()
1172 #
1173 # print ' '
1174
1175 self.BlockCounter += 2
1176
1177 else:
695 else:
1178 self.fileSelector += 1
696 self.fileSelector+=1
1179 self.BlockCounter = 0
697 self.BlockCounter=0
1180 print("Next File") No newline at end of file
Show file before | Show file after | Hide comments
@@ -179,9 +179,6 class ParamReader(JRODataReader,ProcessingUnit):
179 print("[Reading] %d file(s) was(were) found in time range: %s - %s" %(len(filenameList), startTime, endTime))
179 print("[Reading] %d file(s) was(were) found in time range: %s - %s" %(len(filenameList), startTime, endTime))
180 print()
180 print()
181
181
182 # for i in range(len(filenameList)):
183 # print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
184
185 self.filenameList = filenameList
182 self.filenameList = filenameList
186 self.datetimeList = datetimeList
183 self.datetimeList = datetimeList
187
184
@@ -504,20 +501,11 class ParamReader(JRODataReader,ProcessingUnit):
504
501
505 def getData(self):
502 def getData(self):
506
503
507 # if self.flagNoMoreFiles:
508 # self.dataOut.flagNoData = True
509 # print 'Process finished'
510 # return 0
511 #
512 if self.blockIndex==self.blocksPerFile:
504 if self.blockIndex==self.blocksPerFile:
513 if not( self.__setNextFileOffline() ):
505 if not( self.__setNextFileOffline() ):
514 self.dataOut.flagNoData = True
506 self.dataOut.flagNoData = True
515 return 0
507 return 0
516
508
517 # if self.datablock == None: # setear esta condicion cuando no hayan datos por leers
518 # self.dataOut.flagNoData = True
519 # return 0
520 # self.__readData()
521 self.__setDataOut()
509 self.__setDataOut()
522 self.dataOut.flagNoData = False
510 self.dataOut.flagNoData = False
523
511
@@ -637,7 +625,10 class ParamWriter(Operation):
637 dsDict['variable'] = self.dataList[i]
625 dsDict['variable'] = self.dataList[i]
638 #--------------------- Conditionals ------------------------
626 #--------------------- Conditionals ------------------------
639 #There is no data
627 #There is no data
628
629
640 if dataAux is None:
630 if dataAux is None:
631
641 return 0
632 return 0
642
633
643 #Not array, just a number
634 #Not array, just a number
@@ -821,7 +812,7 class ParamWriter(Operation):
821 return False
812 return False
822
813
823 def setNextFile(self):
814 def setNextFile(self):
824
815
825 ext = self.ext
816 ext = self.ext
826 path = self.path
817 path = self.path
827 setFile = self.setFile
818 setFile = self.setFile
@@ -1095,7 +1086,6 class ParamWriter(Operation):
1095 return
1086 return
1096
1087
1097 self.isConfig = True
1088 self.isConfig = True
1098 # self.putMetadata()
1099 self.setNextFile()
1089 self.setNextFile()
1100
1090
1101 self.putData()
1091 self.putData()
Show file before | Show file after | Hide comments
@@ -413,9 +413,7 class SpectraWriter(JRODataWriter, Operation):
413
413
414 data_dc = None
414 data_dc = None
415
415
416 # dataOut = None
416 def __init__(self):
417
418 def __init__(self):#, **kwargs):
419 """
417 """
420 Inicializador de la clase SpectraWriter para la escritura de datos de espectros.
418 Inicializador de la clase SpectraWriter para la escritura de datos de espectros.
421
419
@@ -429,9 +427,7 class SpectraWriter(JRODataWriter, Operation):
429 Return: None
427 Return: None
430 """
428 """
431
429
432 Operation.__init__(self)#, **kwargs)
430 Operation.__init__(self)
433
434 #self.isConfig = False
435
431
436 self.nTotalBlocks = 0
432 self.nTotalBlocks = 0
437
433
@@ -496,7 +492,7 class SpectraWriter(JRODataWriter, Operation):
496
492
497
493
498 def writeBlock(self):
494 def writeBlock(self):
499 """
495 """processingHeaderObj
500 Escribe el buffer en el file designado
496 Escribe el buffer en el file designado
501
497
502 Affected:
498 Affected:
@@ -519,8 +515,10 class SpectraWriter(JRODataWriter, Operation):
519 data.tofile(self.fp)
515 data.tofile(self.fp)
520
516
521 if self.data_cspc is not None:
517 if self.data_cspc is not None:
522 data = numpy.zeros( self.shape_cspc_Buffer, self.dtype )
518
523 cspc = numpy.transpose( self.data_cspc, (0,2,1) )
519 cspc = numpy.transpose( self.data_cspc, (0,2,1) )
520 #data = numpy.zeros( numpy.shape(cspc), self.dtype )
521 #print 'data.shape', self.shape_cspc_Buffer
524 if not self.processingHeaderObj.shif_fft:
522 if not self.processingHeaderObj.shif_fft:
525 cspc = numpy.roll( cspc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones
523 cspc = numpy.roll( cspc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones
526 data['real'] = cspc.real
524 data['real'] = cspc.real
@@ -529,8 +527,9 class SpectraWriter(JRODataWriter, Operation):
529 data.tofile(self.fp)
527 data.tofile(self.fp)
530
528
531 if self.data_dc is not None:
529 if self.data_dc is not None:
532 data = numpy.zeros( self.shape_dc_Buffer, self.dtype )
530
533 dc = self.data_dc
531 dc = self.data_dc
532 data = numpy.zeros( numpy.shape(dc), self.dtype )
534 data['real'] = dc.real
533 data['real'] = dc.real
535 data['imag'] = dc.imag
534 data['imag'] = dc.imag
536 data = data.reshape((-1))
535 data = data.reshape((-1))
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (1209 lines changed) Show them Hide them
@@ -10,11 +10,7 import importlib
10 import itertools
10 import itertools
11 from multiprocessing import Pool, TimeoutError
11 from multiprocessing import Pool, TimeoutError
12 from multiprocessing.pool import ThreadPool
12 from multiprocessing.pool import ThreadPool
13 import types
14 from functools import partial
15 import time
13 import time
16 #from sklearn.cluster import KMeans
17
18
14
19 from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
15 from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
20 from .jroproc_base import ProcessingUnit, Operation, MPDecorator
16 from .jroproc_base import ProcessingUnit, Operation, MPDecorator
@@ -128,6 +124,7 class ParametersProc(ProcessingUnit):
128 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
124 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
129 self.dataOut.spc_noise = self.dataIn.getNoise()
125 self.dataOut.spc_noise = self.dataIn.getNoise()
130 self.dataOut.spc_range = (self.dataIn.getFreqRange(1)/1000. , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
126 self.dataOut.spc_range = (self.dataIn.getFreqRange(1)/1000. , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
127 # self.dataOut.normFactor = self.dataIn.normFactor
131 self.dataOut.pairsList = self.dataIn.pairsList
128 self.dataOut.pairsList = self.dataIn.pairsList
132 self.dataOut.groupList = self.dataIn.pairsList
129 self.dataOut.groupList = self.dataIn.pairsList
133 self.dataOut.flagNoData = False
130 self.dataOut.flagNoData = False
@@ -136,9 +133,9 class ParametersProc(ProcessingUnit):
136 self.dataOut.ChanDist = self.dataIn.ChanDist
133 self.dataOut.ChanDist = self.dataIn.ChanDist
137 else: self.dataOut.ChanDist = None
134 else: self.dataOut.ChanDist = None
138
135
139 if hasattr(self.dataIn, 'VelRange'): #Velocities range
136 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
140 self.dataOut.VelRange = self.dataIn.VelRange
137 # self.dataOut.VelRange = self.dataIn.VelRange
141 else: self.dataOut.VelRange = None
138 #else: self.dataOut.VelRange = None
142
139
143 if hasattr(self.dataIn, 'RadarConst'): #Radar Constant
140 if hasattr(self.dataIn, 'RadarConst'): #Radar Constant
144 self.dataOut.RadarConst = self.dataIn.RadarConst
141 self.dataOut.RadarConst = self.dataIn.RadarConst
@@ -184,9 +181,112 class ParametersProc(ProcessingUnit):
184 def target(tups):
181 def target(tups):
185
182
186 obj, args = tups
183 obj, args = tups
187 #print 'TARGETTT', obj, args
184
188 return obj.FitGau(args)
185 return obj.FitGau(args)
189
186
187
188 class SpectralFilters(Operation):
189
190 '''This class allows the Rainfall / Wind Selection for CLAIRE RADAR
191
192 LimitR : It is the limit in m/s of Rainfall
193 LimitW : It is the limit in m/s for Winds
194
195 Input:
196
197 self.dataOut.data_pre : SPC and CSPC
198 self.dataOut.spc_range : To select wind and rainfall velocities
199
200 Affected:
201
202 self.dataOut.data_pre : It is used for the new SPC and CSPC ranges of wind
203 self.dataOut.spcparam_range : Used in SpcParamPlot
204 self.dataOut.SPCparam : Used in PrecipitationProc
205
206
207 '''
208
209 def __init__(self):
210 Operation.__init__(self)
211 self.i=0
212
213 def run(self, dataOut, PositiveLimit=1.5, NegativeLimit=2.5):
214
215
216 #Limite de vientos
217 LimitR = PositiveLimit
218 LimitN = NegativeLimit
219
220 self.spc = dataOut.data_pre[0].copy()
221 self.cspc = dataOut.data_pre[1].copy()
222
223 self.Num_Hei = self.spc.shape[2]
224 self.Num_Bin = self.spc.shape[1]
225 self.Num_Chn = self.spc.shape[0]
226
227 VelRange = dataOut.spc_range[2]
228 TimeRange = dataOut.spc_range[1]
229 FrecRange = dataOut.spc_range[0]
230
231 Vmax= 2*numpy.max(dataOut.spc_range[2])
232 Tmax= 2*numpy.max(dataOut.spc_range[1])
233 Fmax= 2*numpy.max(dataOut.spc_range[0])
234
235 Breaker1R=VelRange[numpy.abs(VelRange-(-LimitN)).argmin()]
236 Breaker1R=numpy.where(VelRange == Breaker1R)
237
238 Delta = self.Num_Bin/2 - Breaker1R[0]
239
240
241 '''Reacomodando SPCrange'''
242
243 VelRange=numpy.roll(VelRange,-(self.Num_Bin/2) ,axis=0)
244
245 VelRange[-(self.Num_Bin/2):]+= Vmax
246
247 FrecRange=numpy.roll(FrecRange,-(self.Num_Bin/2),axis=0)
248
249 FrecRange[-(self.Num_Bin/2):]+= Fmax
250
251 TimeRange=numpy.roll(TimeRange,-(self.Num_Bin/2),axis=0)
252
253 TimeRange[-(self.Num_Bin/2):]+= Tmax
254
255 ''' ------------------ '''
256
257 Breaker2R=VelRange[numpy.abs(VelRange-(LimitR)).argmin()]
258 Breaker2R=numpy.where(VelRange == Breaker2R)
259
260
261 SPCroll = numpy.roll(self.spc,-(self.Num_Bin/2) ,axis=1)
262
263 SPCcut = SPCroll.copy()
264 for i in range(self.Num_Chn):
265
266 SPCcut[i,0:int(Breaker2R[0]),:] = dataOut.noise[i]
267 SPCcut[i,-int(Delta):,:] = dataOut.noise[i]
268
269 SPCcut[i]=SPCcut[i]- dataOut.noise[i]
270 SPCcut[ numpy.where( SPCcut<0 ) ] = 1e-20
271
272 SPCroll[i]=SPCroll[i]-dataOut.noise[i]
273 SPCroll[ numpy.where( SPCroll<0 ) ] = 1e-20
274
275 SPC_ch1 = SPCroll
276
277 SPC_ch2 = SPCcut
278
279 SPCparam = (SPC_ch1, SPC_ch2, self.spc)
280 dataOut.SPCparam = numpy.asarray(SPCparam)
281
282
283 dataOut.spcparam_range=numpy.zeros([self.Num_Chn,self.Num_Bin+1])
284
285 dataOut.spcparam_range[2]=VelRange
286 dataOut.spcparam_range[1]=TimeRange
287 dataOut.spcparam_range[0]=FrecRange
288
289
190 class GaussianFit(Operation):
290 class GaussianFit(Operation):
191
291
192 '''
292 '''
@@ -198,15 +298,15 class GaussianFit(Operation):
198 self.dataOut.data_pre : SelfSpectra
298 self.dataOut.data_pre : SelfSpectra
199
299
200 Output:
300 Output:
201 self.dataOut.GauSPC : SPC_ch1, SPC_ch2
301 self.dataOut.SPCparam : SPC_ch1, SPC_ch2
202
302
203 '''
303 '''
204 def __init__(self, **kwargs):
304 def __init__(self):
205 Operation.__init__(self, **kwargs)
305 Operation.__init__(self)
206 self.i=0
306 self.i=0
207
307
208
308
209 def run(self, dataOut, num_intg=7, pnoise=1., vel_arr=None, SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
309 def run(self, dataOut, num_intg=7, pnoise=1., SNRlimit=-9): #num_intg: Incoherent integrations, pnoise: Noise, vel_arr: range of velocities, similar to the ftt points
210 """This routine will find a couple of generalized Gaussians to a power spectrum
310 """This routine will find a couple of generalized Gaussians to a power spectrum
211 input: spc
311 input: spc
212 output:
312 output:
@@ -214,31 +314,12 class GaussianFit(Operation):
214 """
314 """
215
315
216 self.spc = dataOut.data_pre[0].copy()
316 self.spc = dataOut.data_pre[0].copy()
217
218
219 print('SelfSpectra Shape', numpy.asarray(self.spc).shape)
220
221
222 #plt.figure(50)
223 #plt.subplot(121)
224 #plt.plot(self.spc,'k',label='spc(66)')
225 #plt.plot(xFrec,ySamples[1],'g',label='Ch1')
226 #plt.plot(xFrec,ySamples[2],'r',label='Ch2')
227 #plt.plot(xFrec,FitGauss,'yo:',label='fit')
228 #plt.legend()
229 #plt.title('DATOS A ALTURA DE 7500 METROS')
230 #plt.show()
231
232 self.Num_Hei = self.spc.shape[2]
317 self.Num_Hei = self.spc.shape[2]
233 #self.Num_Bin = len(self.spc)
234 self.Num_Bin = self.spc.shape[1]
318 self.Num_Bin = self.spc.shape[1]
235 self.Num_Chn = self.spc.shape[0]
319 self.Num_Chn = self.spc.shape[0]
236
237 Vrange = dataOut.abscissaList
320 Vrange = dataOut.abscissaList
238
321
239 #print 'self.spc2', numpy.asarray(self.spc).shape
322 GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei])
240
241 GauSPC = numpy.empty([2,self.Num_Bin,self.Num_Hei])
242 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
323 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
243 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
324 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
244 SPC_ch1[:] = numpy.NaN
325 SPC_ch1[:] = numpy.NaN
@@ -250,272 +331,12 class GaussianFit(Operation):
250 noise_ = dataOut.spc_noise[0].copy()
331 noise_ = dataOut.spc_noise[0].copy()
251
332
252
333
253
254 pool = Pool(processes=self.Num_Chn)
334 pool = Pool(processes=self.Num_Chn)
255 args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
335 args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
256 objs = [self for __ in range(self.Num_Chn)]
336 objs = [self for __ in range(self.Num_Chn)]
257 attrs = list(zip(objs, args))
337 attrs = list(zip(objs, args))
258 gauSPC = pool.map(target, attrs)
338 gauSPC = pool.map(target, attrs)
259 dataOut.GauSPC = numpy.asarray(gauSPC)
339 dataOut.SPCparam = numpy.asarray(SPCparam)
260 # ret = []
261 # for n in range(self.Num_Chn):
262 # self.FitGau(args[n])
263 # dataOut.GauSPC = ret
264
265
266
267 # for ch in range(self.Num_Chn):
268 #
269 # for ht in range(self.Num_Hei):
270 # #print (numpy.asarray(self.spc).shape)
271 # spc = numpy.asarray(self.spc)[ch,:,ht]
272 #
273 # #############################################
274 # # normalizing spc and noise
275 # # This part differs from gg1
276 # spc_norm_max = max(spc)
277 # spc = spc / spc_norm_max
278 # pnoise = pnoise / spc_norm_max
279 # #############################################
280 #
281 # if abs(vel_arr[0])<15.0: # this switch is for spectra collected with different length IPP's
282 # fatspectra=1.0
283 # else:
284 # fatspectra=0.5
285 #
286 # wnoise = noise_ / spc_norm_max
287 # #print 'wnoise', noise_, dataOut.spc_noise[0], wnoise
288 # #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
289 # #if wnoise>1.1*pnoise: # to be tested later
290 # # wnoise=pnoise
291 # noisebl=wnoise*0.9; noisebh=wnoise*1.1
292 # spc=spc-wnoise
293 #
294 # minx=numpy.argmin(spc)
295 # spcs=numpy.roll(spc,-minx)
296 # cum=numpy.cumsum(spcs)
297 # tot_noise=wnoise * self.Num_Bin #64;
298 # #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
299 # #snr=tot_signal/tot_noise
300 # #snr=cum[-1]/tot_noise
301 #
302 # #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
303 #
304 # snr = sum(spcs)/tot_noise
305 # snrdB=10.*numpy.log10(snr)
306 #
307 # #if snrdB < -9 :
308 # # snrdB = numpy.NaN
309 # # continue
310 #
311 # #print 'snr',snrdB # , sum(spcs) , tot_noise
312 #
313 #
314 # #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
315 # # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
316 #
317 # cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region
318 # cumlo=cummax*epsi;
319 # cumhi=cummax*(1-epsi)
320 # powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
321 #
322 # #if len(powerindex)==1:
323 # ##return [numpy.mod(powerindex[0]+minx,64),None,None,None,],[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
324 # #return [numpy.mod(powerindex[0]+minx, self.Num_Bin ),None,None,None,],[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
325 # #elif len(powerindex)<4*fatspectra:
326 # #return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
327 #
328 # if len(powerindex) < 1:# case for powerindex 0
329 # continue
330 # powerlo=powerindex[0]
331 # powerhi=powerindex[-1]
332 # powerwidth=powerhi-powerlo
333 #
334 # firstpeak=powerlo+powerwidth/10.# first gaussian energy location
335 # secondpeak=powerhi-powerwidth/10.#second gaussian energy location
336 # midpeak=(firstpeak+secondpeak)/2.
337 # firstamp=spcs[int(firstpeak)]
338 # secondamp=spcs[int(secondpeak)]
339 # midamp=spcs[int(midpeak)]
340 # #x=numpy.spc.shape[1]
341 #
342 # #x=numpy.arange(64)
343 # x=numpy.arange( self.Num_Bin )
344 # y_data=spc+wnoise
345 #
346 # # single gaussian
347 # #shift0=numpy.mod(midpeak+minx,64)
348 # shift0=numpy.mod(midpeak+minx, self.Num_Bin )
349 # width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
350 # power0=2.
351 # amplitude0=midamp
352 # state0=[shift0,width0,amplitude0,power0,wnoise]
353 # #bnds=((0,63),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
354 # bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(noisebl,noisebh))
355 # #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth),(0,None),(0.5,3.),(0.1,0.5))
356 # # bnds = range of fft, power width, amplitude, power, noise
357 # lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
358 #
359 # chiSq1=lsq1[1];
360 # jack1= self.y_jacobian1(x,lsq1[0])
361 #
362 #
363 # try:
364 # sigmas1=numpy.sqrt(chiSq1*numpy.diag(numpy.linalg.inv(numpy.dot(jack1.T,jack1))))
365 # except:
366 # std1=32.; sigmas1=numpy.ones(5)
367 # else:
368 # std1=sigmas1[0]
369 #
370 #
371 # if fatspectra<1.0 and powerwidth<4:
372 # choice=0
373 # Amplitude0=lsq1[0][2]
374 # shift0=lsq1[0][0]
375 # width0=lsq1[0][1]
376 # p0=lsq1[0][3]
377 # Amplitude1=0.
378 # shift1=0.
379 # width1=0.
380 # p1=0.
381 # noise=lsq1[0][4]
382 # #return (numpy.array([shift0,width0,Amplitude0,p0]),
383 # # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
384 #
385 # # two gaussians
386 # #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
387 # shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
388 # shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
389 # width0=powerwidth/6.;
390 # width1=width0
391 # power0=2.;
392 # power1=power0
393 # amplitude0=firstamp;
394 # amplitude1=secondamp
395 # state0=[shift0,width0,amplitude0,power0,shift1,width1,amplitude1,power1,wnoise]
396 # #bnds=((0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(0,63),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
397 # bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
398 # #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
399 #
400 # lsq2=fmin_l_bfgs_b(self.misfit2,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
401 #
402 #
403 # chiSq2=lsq2[1]; jack2=self.y_jacobian2(x,lsq2[0])
404 #
405 #
406 # try:
407 # sigmas2=numpy.sqrt(chiSq2*numpy.diag(numpy.linalg.inv(numpy.dot(jack2.T,jack2))))
408 # except:
409 # std2a=32.; std2b=32.; sigmas2=numpy.ones(9)
410 # else:
411 # std2a=sigmas2[0]; std2b=sigmas2[4]
412 #
413 #
414 #
415 # oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
416 #
417 # if snrdB>-9: # when SNR is strong pick the peak with least shift (LOS velocity) error
418 # if oneG:
419 # choice=0
420 # else:
421 # w1=lsq2[0][1]; w2=lsq2[0][5]
422 # a1=lsq2[0][2]; a2=lsq2[0][6]
423 # p1=lsq2[0][3]; p2=lsq2[0][7]
424 # s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1; s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
425 # gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
426 #
427 # if gp1>gp2:
428 # if a1>0.7*a2:
429 # choice=1
430 # else:
431 # choice=2
432 # elif gp2>gp1:
433 # if a2>0.7*a1:
434 # choice=2
435 # else:
436 # choice=1
437 # else:
438 # choice=numpy.argmax([a1,a2])+1
439 # #else:
440 # #choice=argmin([std2a,std2b])+1
441 #
442 # else: # with low SNR go to the most energetic peak
443 # choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
444 #
445 # #print 'choice',choice
446 #
447 # if choice==0: # pick the single gaussian fit
448 # Amplitude0=lsq1[0][2]
449 # shift0=lsq1[0][0]
450 # width0=lsq1[0][1]
451 # p0=lsq1[0][3]
452 # Amplitude1=0.
453 # shift1=0.
454 # width1=0.
455 # p1=0.
456 # noise=lsq1[0][4]
457 # elif choice==1: # take the first one of the 2 gaussians fitted
458 # Amplitude0 = lsq2[0][2]
459 # shift0 = lsq2[0][0]
460 # width0 = lsq2[0][1]
461 # p0 = lsq2[0][3]
462 # Amplitude1 = lsq2[0][6] # This is 0 in gg1
463 # shift1 = lsq2[0][4] # This is 0 in gg1
464 # width1 = lsq2[0][5] # This is 0 in gg1
465 # p1 = lsq2[0][7] # This is 0 in gg1
466 # noise = lsq2[0][8]
467 # else: # the second one
468 # Amplitude0 = lsq2[0][6]
469 # shift0 = lsq2[0][4]
470 # width0 = lsq2[0][5]
471 # p0 = lsq2[0][7]
472 # Amplitude1 = lsq2[0][2] # This is 0 in gg1
473 # shift1 = lsq2[0][0] # This is 0 in gg1
474 # width1 = lsq2[0][1] # This is 0 in gg1
475 # p1 = lsq2[0][3] # This is 0 in gg1
476 # noise = lsq2[0][8]
477 #
478 # #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0)
479 # SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
480 # SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
481 # #print 'SPC_ch1.shape',SPC_ch1.shape
482 # #print 'SPC_ch2.shape',SPC_ch2.shape
483 # #dataOut.data_param = SPC_ch1
484 # GauSPC[0] = SPC_ch1
485 # GauSPC[1] = SPC_ch2
486
487 # #plt.gcf().clear()
488 # plt.figure(50+self.i)
489 # self.i=self.i+1
490 # #plt.subplot(121)
491 # plt.plot(self.spc,'k')#,label='spc(66)')
492 # plt.plot(SPC_ch1[ch,ht],'b')#,label='gg1')
493 # #plt.plot(SPC_ch2,'r')#,label='gg2')
494 # #plt.plot(xFrec,ySamples[1],'g',label='Ch1')
495 # #plt.plot(xFrec,ySamples[2],'r',label='Ch2')
496 # #plt.plot(xFrec,FitGauss,'yo:',label='fit')
497 # plt.legend()
498 # plt.title('DATOS A ALTURA DE 7500 METROS')
499 # plt.show()
500 # print 'shift0', shift0
501 # print 'Amplitude0', Amplitude0
502 # print 'width0', width0
503 # print 'p0', p0
504 # print '========================'
505 # print 'shift1', shift1
506 # print 'Amplitude1', Amplitude1
507 # print 'width1', width1
508 # print 'p1', p1
509 # print 'noise', noise
510 # print 's_noise', wnoise
511
512 print('========================================================')
513 print('total_time: ', time.time()-start_time)
514
515 # re-normalizing spc and noise
516 # This part differs from gg1
517
518
519
340
520 ''' Parameters:
341 ''' Parameters:
521 1. Amplitude
342 1. Amplitude
@@ -524,16 +345,11 class GaussianFit(Operation):
524 4. Power
345 4. Power
525 '''
346 '''
526
347
527
528 ###############################################################################
529 def FitGau(self, X):
348 def FitGau(self, X):
530
349
531 Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
350 Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
532 #print 'VARSSSS', ch, pnoise, noise, num_intg
351
533
352 SPCparam = []
534 #print 'HEIGHTS', self.Num_Hei
535
536 GauSPC = []
537 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
353 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
538 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
354 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
539 SPC_ch1[:] = 0#numpy.NaN
355 SPC_ch1[:] = 0#numpy.NaN
@@ -542,10 +358,6 class GaussianFit(Operation):
542
358
543
359
544 for ht in range(self.Num_Hei):
360 for ht in range(self.Num_Hei):
545 #print (numpy.asarray(self.spc).shape)
546
547 #print 'TTTTT', ch , ht
548 #print self.spc.shape
549
361
550
362
551 spc = numpy.asarray(self.spc)[ch,:,ht]
363 spc = numpy.asarray(self.spc)[ch,:,ht]
@@ -554,27 +366,26 class GaussianFit(Operation):
554 # normalizing spc and noise
366 # normalizing spc and noise
555 # This part differs from gg1
367 # This part differs from gg1
556 spc_norm_max = max(spc)
368 spc_norm_max = max(spc)
557 spc = spc / spc_norm_max
369 #spc = spc / spc_norm_max
558 pnoise = pnoise / spc_norm_max
370 pnoise = pnoise #/ spc_norm_max
559 #############################################
371 #############################################
560
372
561 fatspectra=1.0
373 fatspectra=1.0
562
374
563 wnoise = noise_ / spc_norm_max
375 wnoise = noise_ #/ spc_norm_max
564 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
376 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
565 #if wnoise>1.1*pnoise: # to be tested later
377 #if wnoise>1.1*pnoise: # to be tested later
566 # wnoise=pnoise
378 # wnoise=pnoise
567 noisebl=wnoise*0.9; noisebh=wnoise*1.1
379 noisebl=wnoise*0.9;
380 noisebh=wnoise*1.1
568 spc=spc-wnoise
381 spc=spc-wnoise
569 # print 'wnoise', noise_[0], spc_norm_max, wnoise
382
570 minx=numpy.argmin(spc)
383 minx=numpy.argmin(spc)
384 #spcs=spc.copy()
571 spcs=numpy.roll(spc,-minx)
385 spcs=numpy.roll(spc,-minx)
572 cum=numpy.cumsum(spcs)
386 cum=numpy.cumsum(spcs)
573 tot_noise=wnoise * self.Num_Bin #64;
387 tot_noise=wnoise * self.Num_Bin #64;
574 #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
388
575 #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
576 #snr=tot_signal/tot_noise
577 #snr=cum[-1]/tot_noise
578 snr = sum(spcs)/tot_noise
389 snr = sum(spcs)/tot_noise
579 snrdB=10.*numpy.log10(snr)
390 snrdB=10.*numpy.log10(snr)
580
391
@@ -582,16 +393,15 class GaussianFit(Operation):
582 snr = numpy.NaN
393 snr = numpy.NaN
583 SPC_ch1[:,ht] = 0#numpy.NaN
394 SPC_ch1[:,ht] = 0#numpy.NaN
584 SPC_ch1[:,ht] = 0#numpy.NaN
395 SPC_ch1[:,ht] = 0#numpy.NaN
585 GauSPC = (SPC_ch1,SPC_ch2)
396 SPCparam = (SPC_ch1,SPC_ch2)
586 continue
397 continue
587 #print 'snr',snrdB #, sum(spcs) , tot_noise
588
589
398
590
399
591 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
400 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
592 # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
401 # return [None,]*4,[None,]*4,None,snrdB,None,None,[None,]*5,[None,]*9,None
593
402
594 cummax=max(cum); epsi=0.08*fatspectra # cumsum to narrow down the energy region
403 cummax=max(cum);
404 epsi=0.08*fatspectra # cumsum to narrow down the energy region
595 cumlo=cummax*epsi;
405 cumlo=cummax*epsi;
596 cumhi=cummax*(1-epsi)
406 cumhi=cummax*(1-epsi)
597 powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
407 powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
@@ -613,7 +423,7 class GaussianFit(Operation):
613 x=numpy.arange( self.Num_Bin )
423 x=numpy.arange( self.Num_Bin )
614 y_data=spc+wnoise
424 y_data=spc+wnoise
615
425
616 # single gaussian
426 ''' single Gaussian '''
617 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
427 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
618 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
428 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
619 power0=2.
429 power0=2.
@@ -623,16 +433,7 class GaussianFit(Operation):
623 lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
433 lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
624
434
625 chiSq1=lsq1[1];
435 chiSq1=lsq1[1];
626 jack1= self.y_jacobian1(x,lsq1[0])
436
627
628
629 try:
630 sigmas1=numpy.sqrt(chiSq1*numpy.diag(numpy.linalg.inv(numpy.dot(jack1.T,jack1))))
631 except:
632 std1=32.; sigmas1=numpy.ones(5)
633 else:
634 std1=sigmas1[0]
635
636
437
637 if fatspectra<1.0 and powerwidth<4:
438 if fatspectra<1.0 and powerwidth<4:
638 choice=0
439 choice=0
@@ -648,7 +449,7 class GaussianFit(Operation):
648 #return (numpy.array([shift0,width0,Amplitude0,p0]),
449 #return (numpy.array([shift0,width0,Amplitude0,p0]),
649 # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
450 # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
650
451
651 # two gaussians
452 ''' two gaussians '''
652 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
453 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
653 shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
454 shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
654 shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
455 shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
@@ -663,24 +464,16 class GaussianFit(Operation):
663 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
464 bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(noisebl,noisebh))
664 #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
465 #bnds=(( 0,(self.Num_Bin-1) ),(1,powerwidth/2.),(0,None),(0.5,3.),( 0,(self.Num_Bin-1)),(1,powerwidth/2.),(0,None),(0.5,3.),(0.1,0.5))
665
466
666 lsq2=fmin_l_bfgs_b(self.misfit2,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
467 lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )
667
468
668
469
669 chiSq2=lsq2[1]; jack2=self.y_jacobian2(x,lsq2[0])
470 chiSq2=lsq2[1];
471
670
472
671
473
672 try:
673 sigmas2=numpy.sqrt(chiSq2*numpy.diag(numpy.linalg.inv(numpy.dot(jack2.T,jack2))))
674 except:
675 std2a=32.; std2b=32.; sigmas2=numpy.ones(9)
676 else:
677 std2a=sigmas2[0]; std2b=sigmas2[4]
678
679
680
681 oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
474 oneG=(chiSq1<5 and chiSq1/chiSq2<2.0) and (abs(lsq2[0][0]-lsq2[0][4])<(lsq2[0][1]+lsq2[0][5])/3. or abs(lsq2[0][0]-lsq2[0][4])<10)
682
475
683 if snrdB>-6: # when SNR is strong pick the peak with least shift (LOS velocity) error
476 if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
684 if oneG:
477 if oneG:
685 choice=0
478 choice=0
686 else:
479 else:
@@ -690,7 +483,7 class GaussianFit(Operation):
690 s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1;
483 s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1;
691 s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
484 s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
692 gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
485 gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
693
486
694 if gp1>gp2:
487 if gp1>gp2:
695 if a1>0.7*a2:
488 if a1>0.7*a2:
696 choice=1
489 choice=1
@@ -710,13 +503,15 class GaussianFit(Operation):
710 choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
503 choice=numpy.argmax([lsq1[0][2]*lsq1[0][1],lsq2[0][2]*lsq2[0][1],lsq2[0][6]*lsq2[0][5]])
711
504
712
505
713 shift0=lsq2[0][0]; vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
506 shift0=lsq2[0][0];
714 shift1=lsq2[0][4]; vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
507 vel0=Vrange[0] + shift0*(Vrange[1]-Vrange[0])
508 shift1=lsq2[0][4];
509 vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
715
510
716 max_vel = 20
511 max_vel = 1.0
717
512
718 #first peak will be 0, second peak will be 1
513 #first peak will be 0, second peak will be 1
719 if vel0 > 0 and vel0 < max_vel : #first peak is in the correct range
514 if vel0 > -1.0 and vel0 < max_vel : #first peak is in the correct range
720 shift0=lsq2[0][0]
515 shift0=lsq2[0][0]
721 width0=lsq2[0][1]
516 width0=lsq2[0][1]
722 Amplitude0=lsq2[0][2]
517 Amplitude0=lsq2[0][2]
@@ -739,120 +534,18 class GaussianFit(Operation):
739 p0=lsq2[0][7]
534 p0=lsq2[0][7]
740 noise=lsq2[0][8]
535 noise=lsq2[0][8]
741
536
742 if Amplitude0<0.1: # in case the peak is noise
537 if Amplitude0<0.05: # in case the peak is noise
743 shift0,width0,Amplitude0,p0 = 4*[numpy.NaN]
538 shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN]
744 if Amplitude1<0.1:
539 if Amplitude1<0.05:
745 shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
540 shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN]
746
541
747
542
748 # if choice==0: # pick the single gaussian fit
749 # Amplitude0=lsq1[0][2]
750 # shift0=lsq1[0][0]
751 # width0=lsq1[0][1]
752 # p0=lsq1[0][3]
753 # Amplitude1=0.
754 # shift1=0.
755 # width1=0.
756 # p1=0.
757 # noise=lsq1[0][4]
758 # elif choice==1: # take the first one of the 2 gaussians fitted
759 # Amplitude0 = lsq2[0][2]
760 # shift0 = lsq2[0][0]
761 # width0 = lsq2[0][1]
762 # p0 = lsq2[0][3]
763 # Amplitude1 = lsq2[0][6] # This is 0 in gg1
764 # shift1 = lsq2[0][4] # This is 0 in gg1
765 # width1 = lsq2[0][5] # This is 0 in gg1
766 # p1 = lsq2[0][7] # This is 0 in gg1
767 # noise = lsq2[0][8]
768 # else: # the second one
769 # Amplitude0 = lsq2[0][6]
770 # shift0 = lsq2[0][4]
771 # width0 = lsq2[0][5]
772 # p0 = lsq2[0][7]
773 # Amplitude1 = lsq2[0][2] # This is 0 in gg1
774 # shift1 = lsq2[0][0] # This is 0 in gg1
775 # width1 = lsq2[0][1] # This is 0 in gg1
776 # p1 = lsq2[0][3] # This is 0 in gg1
777 # noise = lsq2[0][8]
778
779 #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0)
780 SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
543 SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
781 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
544 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
782 #print 'SPC_ch1.shape',SPC_ch1.shape
545 SPCparam = (SPC_ch1,SPC_ch2)
783 #print 'SPC_ch2.shape',SPC_ch2.shape
784 #dataOut.data_param = SPC_ch1
785 GauSPC = (SPC_ch1,SPC_ch2)
786 #GauSPC[1] = SPC_ch2
787
788 # print 'shift0', shift0
789 # print 'Amplitude0', Amplitude0
790 # print 'width0', width0
791 # print 'p0', p0
792 # print '========================'
793 # print 'shift1', shift1
794 # print 'Amplitude1', Amplitude1
795 # print 'width1', width1
796 # print 'p1', p1
797 # print 'noise', noise
798 # print 's_noise', wnoise
799
546
800 return GauSPC
801
802
803 def y_jacobian1(self,x,state): # This function is for further analysis of generalized Gaussians, it is not too importan for the signal discrimination.
804 y_model=self.y_model1(x,state)
805 s0,w0,a0,p0,n=state
806 e0=((x-s0)/w0)**2;
807
808 e0u=((x-s0-self.Num_Bin)/w0)**2;
809
810 e0d=((x-s0+self.Num_Bin)/w0)**2
811 m0=numpy.exp(-0.5*e0**(p0/2.));
812 m0u=numpy.exp(-0.5*e0u**(p0/2.));
813 m0d=numpy.exp(-0.5*e0d**(p0/2.))
814 JA=m0+m0u+m0d
815 JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d)
816
817 JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)
818
819 JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2
820 jack1=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,1./y_model])
821 return jack1.T
822
823 def y_jacobian2(self,x,state):
824 y_model=self.y_model2(x,state)
825 s0,w0,a0,p0,s1,w1,a1,p1,n=state
826 e0=((x-s0)/w0)**2;
827
828 e0u=((x-s0- self.Num_Bin )/w0)**2;
829
830 e0d=((x-s0+ self.Num_Bin )/w0)**2
831 e1=((x-s1)/w1)**2;
832
547
833 e1u=((x-s1- self.Num_Bin )/w1)**2;
548 return GauSPC
834
835 e1d=((x-s1+ self.Num_Bin )/w1)**2
836 m0=numpy.exp(-0.5*e0**(p0/2.));
837 m0u=numpy.exp(-0.5*e0u**(p0/2.));
838 m0d=numpy.exp(-0.5*e0d**(p0/2.))
839 m1=numpy.exp(-0.5*e1**(p1/2.));
840 m1u=numpy.exp(-0.5*e1u**(p1/2.));
841 m1d=numpy.exp(-0.5*e1d**(p1/2.))
842 JA=m0+m0u+m0d
843 JA1=m1+m1u+m1d
844 JP=(-1/4.)*a0*m0*e0**(p0/2.)*numpy.log(e0)+(-1/4.)*a0*m0u*e0u**(p0/2.)*numpy.log(e0u)+(-1/4.)*a0*m0d*e0d**(p0/2.)*numpy.log(e0d)
845 JP1=(-1/4.)*a1*m1*e1**(p1/2.)*numpy.log(e1)+(-1/4.)*a1*m1u*e1u**(p1/2.)*numpy.log(e1u)+(-1/4.)*a1*m1d*e1d**(p1/2.)*numpy.log(e1d)
846
847 JS=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)
848
849 JS1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1)
850
851 JW=(p0/w0/2.)*a0*m0*e0**(p0/2.-1)*((x-s0)/w0)**2+(p0/w0/2.)*a0*m0u*e0u**(p0/2.-1)*((x-s0- self.Num_Bin )/w0)**2+(p0/w0/2.)*a0*m0d*e0d**(p0/2.-1)*((x-s0+ self.Num_Bin )/w0)**2
852
853 JW1=(p1/w1/2.)*a1*m1*e1**(p1/2.-1)*((x-s1)/w1)**2+(p1/w1/2.)*a1*m1u*e1u**(p1/2.-1)*((x-s1- self.Num_Bin )/w1)**2+(p1/w1/2.)*a1*m1d*e1d**(p1/2.-1)*((x-s1+ self.Num_Bin )/w1)**2
854 jack2=numpy.sqrt(7)*numpy.array([JS/y_model,JW/y_model,JA/y_model,JP/y_model,JS1/y_model,JW1/y_model,JA1/y_model,JP1/y_model,1./y_model])
855 return jack2.T
856
549
857 def y_model1(self,x,state):
550 def y_model1(self,x,state):
858 shift0,width0,amplitude0,power0,noise=state
551 shift0,width0,amplitude0,power0,noise=state
@@ -884,6 +577,7 class GaussianFit(Operation):
884 def misfit2(self,state,y_data,x,num_intg):
577 def misfit2(self,state,y_data,x,num_intg):
885 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
578 return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.)
886
579
580
887
581
888 class PrecipitationProc(Operation):
582 class PrecipitationProc(Operation):
889
583
@@ -900,24 +594,61 class PrecipitationProc(Operation):
900
594
901 Parameters affected:
595 Parameters affected:
902 '''
596 '''
903
904
597
905 def run(self, dataOut, radar=None, Pt=None, Gt=None, Gr=None, Lambda=None, aL=None,
598 def __init__(self):
906 tauW=None, ThetaT=None, ThetaR=None, Km = 0.93, Altitude=None):
599 Operation.__init__(self)
600 self.i=0
601
602
603 def gaus(self,xSamples,Amp,Mu,Sigma):
604 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
605
606
607
608 def Moments(self, ySamples, xSamples):
609 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
610 yNorm = ySamples / Pot
611
612 Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
613 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
614 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
615
616 return numpy.array([Pot, Vr, Desv])
617
618 def run(self, dataOut, radar=None, Pt=5000, Gt=295.1209, Gr=70.7945, Lambda=0.6741, aL=2.5118,
619 tauW=4e-06, ThetaT=0.1656317, ThetaR=0.36774087, Km = 0.93, Altitude=3350):
907
620
908 self.spc = dataOut.data_pre[0].copy()
909 self.Num_Hei = self.spc.shape[2]
910 self.Num_Bin = self.spc.shape[1]
911 self.Num_Chn = self.spc.shape[0]
912
621
913 Velrange = dataOut.abscissaList
622 Velrange = dataOut.spcparam_range[2]
623 FrecRange = dataOut.spcparam_range[0]
624
625 dV= Velrange[1]-Velrange[0]
626 dF= FrecRange[1]-FrecRange[0]
914
627
915 if radar == "MIRA35C" :
628 if radar == "MIRA35C" :
916
629
630 self.spc = dataOut.data_pre[0].copy()
631 self.Num_Hei = self.spc.shape[2]
632 self.Num_Bin = self.spc.shape[1]
633 self.Num_Chn = self.spc.shape[0]
917 Ze = self.dBZeMODE2(dataOut)
634 Ze = self.dBZeMODE2(dataOut)
918
635
919 else:
636 else:
920
637
638 self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() #
639
640 """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
641
642 self.spc[:,:,0:7]= numpy.NaN
643
644 """##########################################"""
645
646 self.Num_Hei = self.spc.shape[2]
647 self.Num_Bin = self.spc.shape[1]
648 self.Num_Chn = self.spc.shape[0]
649
650 ''' Se obtiene la constante del RADAR '''
651
921 self.Pt = Pt
652 self.Pt = Pt
922 self.Gt = Gt
653 self.Gt = Gt
923 self.Gr = Gr
654 self.Gr = Gr
@@ -927,48 +658,101 class PrecipitationProc(Operation):
927 self.ThetaT = ThetaT
658 self.ThetaT = ThetaT
928 self.ThetaR = ThetaR
659 self.ThetaR = ThetaR
929
660
930 RadarConstant = GetRadarConstant()
661 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
931 SPCmean = numpy.mean(self.spc,0)
662 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
932 ETA = numpy.zeros(self.Num_Hei)
663 RadarConstant = 5e-26 * Numerator / Denominator #
933 Pr = numpy.sum(SPCmean,0)
934
664
935 #for R in range(self.Num_Hei):
665 ''' ============================= '''
936 # ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
666
937
667 self.spc[0] = (self.spc[0]-dataOut.noise[0])
938 D_range = numpy.zeros(self.Num_Hei)
668 self.spc[1] = (self.spc[1]-dataOut.noise[1])
939 EqSec = numpy.zeros(self.Num_Hei)
669 self.spc[2] = (self.spc[2]-dataOut.noise[2])
670
671 self.spc[ numpy.where(self.spc < 0)] = 0
672
673 SPCmean = (numpy.mean(self.spc,0) - numpy.mean(dataOut.noise))
674 SPCmean[ numpy.where(SPCmean < 0)] = 0
675
676 ETAn = numpy.zeros([self.Num_Bin,self.Num_Hei])
677 ETAv = numpy.zeros([self.Num_Bin,self.Num_Hei])
678 ETAd = numpy.zeros([self.Num_Bin,self.Num_Hei])
679
680 Pr = SPCmean[:,:]
681
682 VelMeteoro = numpy.mean(SPCmean,axis=0)
683
684 D_range = numpy.zeros([self.Num_Bin,self.Num_Hei])
685 SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei])
686 N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei])
687 V_mean = numpy.zeros(self.Num_Hei)
940 del_V = numpy.zeros(self.Num_Hei)
688 del_V = numpy.zeros(self.Num_Hei)
689 Z = numpy.zeros(self.Num_Hei)
690 Ze = numpy.zeros(self.Num_Hei)
691 RR = numpy.zeros(self.Num_Hei)
692
693 Range = dataOut.heightList*1000.
941
694
942 for R in range(self.Num_Hei):
695 for R in range(self.Num_Hei):
943 ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
944
696
945 h = R + Altitude #Range from ground to radar pulse altitude
697 h = Range[R] + Altitude #Range from ground to radar pulse altitude
946 del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
698 del_V[R] = 1 + 3.68 * 10**-5 * h + 1.71 * 10**-9 * h**2 #Density change correction for velocity
947
699
948 D_range[R] = numpy.log( (9.65 - (Velrange[R]/del_V[R])) / 10.3 ) / -0.6 #Range of Diameter of drops related to velocity
700 D_range[:,R] = numpy.log( (9.65 - (Velrange[0:self.Num_Bin] / del_V[R])) / 10.3 ) / -0.6 #Diameter range [m]x10**-3
949 SIGMA[R] = numpy.pi**5 / Lambda**4 * Km * D_range[R]**6 #Equivalent Section of drops (sigma)
701
702 '''NOTA: ETA(n) dn = ETA(f) df
703
704 dn = 1 Diferencial de muestreo
705 df = ETA(n) / ETA(f)
706
707 '''
708
709 ETAn[:,R] = RadarConstant * Pr[:,R] * (Range[R] )**2 #Reflectivity (ETA)
710
711 ETAv[:,R]=ETAn[:,R]/dV
712
713 ETAd[:,R]=ETAv[:,R]*6.18*exp(-0.6*D_range[:,R])
714
715 SIGMA[:,R] = Km * (D_range[:,R] * 1e-3 )**6 * numpy.pi**5 / Lambda**4 #Equivalent Section of drops (sigma)
716
717 N_dist[:,R] = ETAn[:,R] / SIGMA[:,R]
718
719 DMoments = self.Moments(Pr[:,R], Velrange[0:self.Num_Bin])
720
721 try:
722 popt01,pcov = curve_fit(self.gaus, Velrange[0:self.Num_Bin] , Pr[:,R] , p0=DMoments)
723 except:
724 popt01=numpy.zeros(3)
725 popt01[1]= DMoments[1]
726
727 if popt01[1]<0 or popt01[1]>20:
728 popt01[1]=numpy.NaN
729
730
731 V_mean[R]=popt01[1]
732
733 Z[R] = numpy.nansum( N_dist[:,R] * (D_range[:,R])**6 )#*10**-18
734
735 RR[R] = 0.0006*numpy.pi * numpy.nansum( D_range[:,R]**3 * N_dist[:,R] * Velrange[0:self.Num_Bin] ) #Rainfall rate
736
737 Ze[R] = (numpy.nansum( ETAn[:,R]) * Lambda**4) / ( 10**-18*numpy.pi**5 * Km)
950
738
951 N_dist[R] = ETA[R] / SIGMA[R]
952
953 Ze = (ETA * Lambda**4) / (numpy.pi * Km)
954 Z = numpy.sum( N_dist * D_range**6 )
955 RR = 6*10**-4*numpy.pi * numpy.sum( D_range**3 * N_dist * Velrange ) #Rainfall rate
956
739
957
740
958 RR = (Ze/200)**(1/1.6)
741 RR2 = (Z/200)**(1/1.6)
959 dBRR = 10*numpy.log10(RR)
742 dBRR = 10*numpy.log10(RR)
743 dBRR2 = 10*numpy.log10(RR2)
960
744
961 dBZe = 10*numpy.log10(Ze)
745 dBZe = 10*numpy.log10(Ze)
962 dataOut.data_output = Ze
746 dBZ = 10*numpy.log10(Z)
963 dataOut.data_param = numpy.ones([2,self.Num_Hei])
747
964 dataOut.channelList = [0,1]
748 dataOut.data_output = RR[8]
965 print('channelList', dataOut.channelList)
749 dataOut.data_param = numpy.ones([3,self.Num_Hei])
966 dataOut.data_param[0]=dBZe
750 dataOut.channelList = [0,1,2]
967 dataOut.data_param[1]=dBRR
751
968 print('RR SHAPE', dBRR.shape)
752 dataOut.data_param[0]=dBZ
969 print('Ze SHAPE', dBZe.shape)
753 dataOut.data_param[1]=V_mean
970 print('dataOut.data_param SHAPE', dataOut.data_param.shape)
754 dataOut.data_param[2]=RR
971
755
972
756
973 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
757 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
974
758
@@ -983,7 +767,7 class PrecipitationProc(Operation):
983 data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN
767 data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN
984
768
985 ETA = numpy.sum(SNR,1)
769 ETA = numpy.sum(SNR,1)
986 print('ETA' , ETA)
770
987 ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
771 ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
988
772
989 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
773 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
@@ -995,26 +779,27 class PrecipitationProc(Operation):
995
779
996 return Ze
780 return Ze
997
781
998 def GetRadarConstant(self):
782 # def GetRadarConstant(self):
999
783 #
1000 """
784 # """
1001 Constants:
785 # Constants:
1002
786 #
1003 Pt: Transmission Power dB
787 # Pt: Transmission Power dB 5kW 5000
1004 Gt: Transmission Gain dB
788 # Gt: Transmission Gain dB 24.7 dB 295.1209
1005 Gr: Reception Gain dB
789 # Gr: Reception Gain dB 18.5 dB 70.7945
1006 Lambda: Wavelenght m
790 # Lambda: Wavelenght m 0.6741 m 0.6741
1007 aL: Attenuation loses dB
791 # aL: Attenuation loses dB 4dB 2.5118
1008 tauW: Width of transmission pulse s
792 # tauW: Width of transmission pulse s 4us 4e-6
1009 ThetaT: Transmission antenna bean angle rad
793 # ThetaT: Transmission antenna bean angle rad 0.1656317 rad 0.1656317
1010 ThetaR: Reception antenna beam angle rad
794 # ThetaR: Reception antenna beam angle rad 0.36774087 rad 0.36774087
1011
795 #
1012 """
796 # """
1013 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
797 #
1014 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
798 # Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
1015 RadarConstant = Numerator / Denominator
799 # Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
1016
800 # RadarConstant = Numerator / Denominator
1017 return RadarConstant
801 #
802 # return RadarConstant
1018
803
1019
804
1020
805
@@ -1037,10 +822,20 class FullSpectralAnalysis(Operation):
1037 Parameters affected: Winds, height range, SNR
822 Parameters affected: Winds, height range, SNR
1038
823
1039 """
824 """
1040 def run(self, dataOut, E01=None, E02=None, E12=None, N01=None, N02=None, N12=None, SNRlimit=7):
825 def run(self, dataOut, Xi01=None, Xi02=None, Xi12=None, Eta01=None, Eta02=None, Eta12=None, SNRlimit=7):
826
827 self.indice=int(numpy.random.rand()*1000)
1041
828
1042 spc = dataOut.data_pre[0].copy()
829 spc = dataOut.data_pre[0].copy()
1043 cspc = dataOut.data_pre[1].copy()
830 cspc = dataOut.data_pre[1]
831
832 """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
833
834 SNRspc = spc.copy()
835 SNRspc[:,:,0:7]= numpy.NaN
836
837 """##########################################"""
838
1044
839
1045 nChannel = spc.shape[0]
840 nChannel = spc.shape[0]
1046 nProfiles = spc.shape[1]
841 nProfiles = spc.shape[1]
@@ -1050,14 +845,9 class FullSpectralAnalysis(Operation):
1050 if dataOut.ChanDist is not None :
845 if dataOut.ChanDist is not None :
1051 ChanDist = dataOut.ChanDist
846 ChanDist = dataOut.ChanDist
1052 else:
847 else:
1053 ChanDist = numpy.array([[E01, N01],[E02,N02],[E12,N12]])
848 ChanDist = numpy.array([[Xi01, Eta01],[Xi02,Eta02],[Xi12,Eta12]])
1054
1055 #print 'ChanDist', ChanDist
1056
849
1057 if dataOut.VelRange is not None:
850 FrecRange = dataOut.spc_range[0]
1058 VelRange= dataOut.VelRange
1059 else:
1060 VelRange= dataOut.abscissaList
1061
851
1062 ySamples=numpy.ones([nChannel,nProfiles])
852 ySamples=numpy.ones([nChannel,nProfiles])
1063 phase=numpy.ones([nChannel,nProfiles])
853 phase=numpy.ones([nChannel,nProfiles])
@@ -1065,113 +855,108 class FullSpectralAnalysis(Operation):
1065 coherence=numpy.ones([nChannel,nProfiles])
855 coherence=numpy.ones([nChannel,nProfiles])
1066 PhaseSlope=numpy.ones(nChannel)
856 PhaseSlope=numpy.ones(nChannel)
1067 PhaseInter=numpy.ones(nChannel)
857 PhaseInter=numpy.ones(nChannel)
1068 dataSNR = dataOut.data_SNR
858 data_SNR=numpy.zeros([nProfiles])
1069
1070
1071
859
1072 data = dataOut.data_pre
860 data = dataOut.data_pre
1073 noise = dataOut.noise
861 noise = dataOut.noise
1074 print('noise',noise)
1075 #SNRdB = 10*numpy.log10(dataOut.data_SNR)
1076
862
1077 FirstMoment = numpy.average(dataOut.data_param[:,1,:],0)
863 dataOut.data_SNR = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0]
1078 #SNRdBMean = []
1079
1080
864
1081 #for j in range(nHeights):
865 dataOut.data_SNR[numpy.where( dataOut.data_SNR <0 )] = 1e-20
1082 # FirstMoment = numpy.append(FirstMoment,numpy.mean([dataOut.data_param[0,1,j],dataOut.data_param[1,1,j],dataOut.data_param[2,1,j]]))
866
1083 # SNRdBMean = numpy.append(SNRdBMean,numpy.mean([SNRdB[0,j],SNRdB[1,j],SNRdB[2,j]]))
867
1084
868 data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN
1085 data_output=numpy.ones([3,spc.shape[2]])*numpy.NaN
1086
869
1087 velocityX=[]
870 velocityX=[]
1088 velocityY=[]
871 velocityY=[]
1089 velocityV=[]
872 velocityV=[]
873 PhaseLine=[]
1090
874
1091 dbSNR = 10*numpy.log10(dataSNR)
875 dbSNR = 10*numpy.log10(dataOut.data_SNR)
1092 dbSNR = numpy.average(dbSNR,0)
876 dbSNR = numpy.average(dbSNR,0)
877
1093 for Height in range(nHeights):
878 for Height in range(nHeights):
1094
879
1095 [Vzon,Vmer,Vver, GaussCenter]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR[Height], SNRlimit)
880 [Vzon,Vmer,Vver, GaussCenter, PhaseSlope, FitGaussCSPC]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, dataOut.spc_range.copy(), dbSNR[Height], SNRlimit)
881 PhaseLine = numpy.append(PhaseLine, PhaseSlope)
1096
882
1097 if abs(Vzon)<100. and abs(Vzon)> 0.:
883 if abs(Vzon)<100. and abs(Vzon)> 0.:
1098 velocityX=numpy.append(velocityX, Vzon)#Vmag
884 velocityX=numpy.append(velocityX, Vzon)#Vmag
1099
885
1100 else:
886 else:
1101 print('Vzon',Vzon)
1102 velocityX=numpy.append(velocityX, numpy.NaN)
887 velocityX=numpy.append(velocityX, numpy.NaN)
1103
888
1104 if abs(Vmer)<100. and abs(Vmer) > 0.:
889 if abs(Vmer)<100. and abs(Vmer) > 0.:
1105 velocityY=numpy.append(velocityY, Vmer)#Vang
890 velocityY=numpy.append(velocityY, -Vmer)#Vang
1106
891
1107 else:
892 else:
1108 print('Vmer',Vmer)
1109 velocityY=numpy.append(velocityY, numpy.NaN)
893 velocityY=numpy.append(velocityY, numpy.NaN)
1110
894
1111 if dbSNR[Height] > SNRlimit:
895 if dbSNR[Height] > SNRlimit:
1112 velocityV=numpy.append(velocityV, FirstMoment[Height])
896 velocityV=numpy.append(velocityV, -Vver)#FirstMoment[Height])
1113 else:
897 else:
1114 velocityV=numpy.append(velocityV, numpy.NaN)
898 velocityV=numpy.append(velocityV, numpy.NaN)
1115 #FirstMoment[Height]= numpy.NaN
899
1116 # if SNRdBMean[Height] <12:
900
1117 # FirstMoment[Height] = numpy.NaN
1118 # velocityX[Height] = numpy.NaN
1119 # velocityY[Height] = numpy.NaN
1120
1121
1122 data_output[0]=numpy.array(velocityX)
1123 data_output[1]=numpy.array(velocityY)
1124 data_output[2]=-velocityV#FirstMoment
1125
1126 print(' ')
1127 #print 'FirstMoment'
1128 #print FirstMoment
1129 print('velocityX',data_output[0])
1130 print(' ')
1131 print('velocityY',data_output[1])
1132 #print numpy.array(velocityY)
1133 print(' ')
1134 #print 'SNR'
1135 #print 10*numpy.log10(dataOut.data_SNR)
1136 #print numpy.shape(10*numpy.log10(dataOut.data_SNR))
1137 print(' ')
1138
901
902 '''Nota: Cambiar el signo de numpy.array(velocityX) cuando se intente procesar datos de BLTR'''
903 data_output[0] = numpy.array(velocityX) #self.moving_average(numpy.array(velocityX) , N=1)
904 data_output[1] = numpy.array(velocityY) #self.moving_average(numpy.array(velocityY) , N=1)
905 data_output[2] = velocityV#FirstMoment
906
907 xFrec=FrecRange[0:spc.shape[1]]
1139
908
1140 dataOut.data_output=data_output
909 dataOut.data_output=data_output
910
1141 return
911 return
1142
912
1143
913
1144 def moving_average(self,x, N=2):
914 def moving_average(self,x, N=2):
1145 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
915 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
1146
916
1147 def gaus(self,xSamples,a,x0,sigma):
917 def gaus(self,xSamples,Amp,Mu,Sigma):
1148 return a*numpy.exp(-(xSamples-x0)**2/(2*sigma**2))
918 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
919
920
1149
921
1150 def Find(self,x,value):
922 def Moments(self, ySamples, xSamples):
1151 for index in range(len(x)):
923 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
1152 if x[index]==value:
924 yNorm = ySamples / Pot
1153 return index
925 Vr = numpy.nansum( yNorm * xSamples ) # Velocidad radial, mu, corrimiento doppler, primer momento
926 Sigma2 = abs(numpy.nansum( yNorm * ( xSamples - Vr )**2 )) # Segundo Momento
927 Desv = Sigma2**0.5 # Desv. Estandar, Ancho espectral
928
929 return numpy.array([Pot, Vr, Desv])
1154
930
1155 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR, SNRlimit):
931 def WindEstimation(self, spc, cspc, pairsList, ChanDist, Height, noise, AbbsisaRange, dbSNR, SNRlimit):
932
933
1156
934
1157 ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
935 ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
1158 phase=numpy.ones([spc.shape[0],spc.shape[1]])
936 phase=numpy.ones([spc.shape[0],spc.shape[1]])
1159 CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_)
937 CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_)
1160 coherence=numpy.ones([spc.shape[0],spc.shape[1]])
938 coherence=numpy.ones([spc.shape[0],spc.shape[1]])
1161 PhaseSlope=numpy.ones(spc.shape[0])
939 PhaseSlope=numpy.zeros(spc.shape[0])
1162 PhaseInter=numpy.ones(spc.shape[0])
940 PhaseInter=numpy.ones(spc.shape[0])
1163 xFrec=VelRange
941 xFrec=AbbsisaRange[0][0:spc.shape[1]]
942 xVel =AbbsisaRange[2][0:spc.shape[1]]
943 Vv=numpy.empty(spc.shape[2])*0
944 SPCav = numpy.average(spc, axis=0)-numpy.average(noise) #spc[0]-noise[0]#
945
946 SPCmoments = self.Moments(SPCav[:,Height], xVel )
947 CSPCmoments = []
948 cspcNoise = numpy.empty(3)
1164
949
1165 '''Getting Eij and Nij'''
950 '''Getting Eij and Nij'''
1166
951
1167 E01=ChanDist[0][0]
952 Xi01=ChanDist[0][0]
1168 N01=ChanDist[0][1]
953 Eta01=ChanDist[0][1]
1169
954
1170 E02=ChanDist[1][0]
955 Xi02=ChanDist[1][0]
1171 N02=ChanDist[1][1]
956 Eta02=ChanDist[1][1]
1172
957
1173 E12=ChanDist[2][0]
958 Xi12=ChanDist[2][0]
1174 N12=ChanDist[2][1]
959 Eta12=ChanDist[2][1]
1175
960
1176 z = spc.copy()
961 z = spc.copy()
1177 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
962 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
@@ -1179,176 +964,197 class FullSpectralAnalysis(Operation):
1179 for i in range(spc.shape[0]):
964 for i in range(spc.shape[0]):
1180
965
1181 '''****** Line of Data SPC ******'''
966 '''****** Line of Data SPC ******'''
1182 zline=z[i,:,Height]
967 zline=z[i,:,Height].copy() - noise[i] # Se resta ruido
1183
968
1184 '''****** SPC is normalized ******'''
969 '''****** SPC is normalized ******'''
1185 FactNorm= (zline.copy()-noise[i]) / numpy.sum(zline.copy())
970 SmoothSPC =self.moving_average(zline.copy(),N=1) # Se suaviza el ruido
1186 FactNorm= FactNorm/numpy.sum(FactNorm)
971 FactNorm = SmoothSPC/numpy.nansum(SmoothSPC) # SPC Normalizado y suavizado
1187
972
1188 SmoothSPC=self.moving_average(FactNorm,N=3)
973 xSamples = xFrec # Se toma el rango de frecuncias
1189
974 ySamples[i] = FactNorm # Se toman los valores de SPC normalizado
1190 xSamples = ar(list(range(len(SmoothSPC))))
1191 ySamples[i] = SmoothSPC
1192
1193 #dbSNR=10*numpy.log10(dataSNR)
1194 print(' ')
1195 print(' ')
1196 print(' ')
1197
1198 #print 'dataSNR', dbSNR.shape, dbSNR[0,40:120]
1199 print('SmoothSPC', SmoothSPC.shape, SmoothSPC[0:20])
1200 print('noise',noise)
1201 print('zline',zline.shape, zline[0:20])
1202 print('FactNorm',FactNorm.shape, FactNorm[0:20])
1203 print('FactNorm suma', numpy.sum(FactNorm))
1204
975
1205 for i in range(spc.shape[0]):
976 for i in range(spc.shape[0]):
1206
977
1207 '''****** Line of Data CSPC ******'''
978 '''****** Line of Data CSPC ******'''
1208 cspcLine=cspc[i,:,Height].copy()
979 cspcLine = ( cspc[i,:,Height].copy())# - noise[i] ) # no! Se resta el ruido
980 SmoothCSPC =self.moving_average(cspcLine,N=1) # Se suaviza el ruido
981 cspcNorm = SmoothCSPC/numpy.nansum(SmoothCSPC) # CSPC normalizado y suavizado
1209
982
1210 '''****** CSPC is normalized ******'''
983 '''****** CSPC is normalized with respect to Briggs and Vincent ******'''
1211 chan_index0 = pairsList[i][0]
984 chan_index0 = pairsList[i][0]
1212 chan_index1 = pairsList[i][1]
985 chan_index1 = pairsList[i][1]
1213 CSPCFactor= abs(numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])) #
1214
986
1215 CSPCNorm = (cspcLine.copy() -noise[i]) / numpy.sqrt(CSPCFactor)
987 CSPCFactor= numpy.abs(numpy.nansum(ySamples[chan_index0]))**2 * numpy.abs(numpy.nansum(ySamples[chan_index1]))**2
988 CSPCNorm = cspcNorm / numpy.sqrt(CSPCFactor)
1216
989
1217 CSPCSamples[i] = CSPCNorm
990 CSPCSamples[i] = CSPCNorm
991
1218 coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
992 coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
1219
993
1220 coherence[i]= self.moving_average(coherence[i],N=2)
994 #coherence[i]= self.moving_average(coherence[i],N=1)
1221
995
1222 phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
996 phase[i] = self.moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
1223
997
1224 print('cspcLine', cspcLine.shape, cspcLine[0:20])
998 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPCSamples[0]), xSamples),
1225 print('CSPCFactor', CSPCFactor)#, CSPCFactor[0:20]
999 self.Moments(numpy.abs(CSPCSamples[1]), xSamples),
1226 print(numpy.sum(ySamples[chan_index0]), numpy.sum(ySamples[chan_index1]), -noise[i])
1000 self.Moments(numpy.abs(CSPCSamples[2]), xSamples)])
1227 print('CSPCNorm', CSPCNorm.shape, CSPCNorm[0:20])
1228 print('CSPCNorm suma', numpy.sum(CSPCNorm))
1229 print('CSPCSamples', CSPCSamples.shape, CSPCSamples[0,0:20])
1230
1001
1231 '''****** Getting fij width ******'''
1002
1003 popt=[1e-10,0,1e-10]
1004 popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
1005 FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0
1006
1007 CSPCMask01 = numpy.abs(CSPCSamples[0])
1008 CSPCMask02 = numpy.abs(CSPCSamples[1])
1009 CSPCMask12 = numpy.abs(CSPCSamples[2])
1010
1011 mask01 = ~numpy.isnan(CSPCMask01)
1012 mask02 = ~numpy.isnan(CSPCMask02)
1013 mask12 = ~numpy.isnan(CSPCMask12)
1232
1014
1233 yMean=[]
1015 #mask = ~numpy.isnan(CSPCMask01)
1234 yMean2=[]
1016 CSPCMask01 = CSPCMask01[mask01]
1017 CSPCMask02 = CSPCMask02[mask02]
1018 CSPCMask12 = CSPCMask12[mask12]
1019 #CSPCMask01 = numpy.ma.masked_invalid(CSPCMask01)
1235
1020
1236 for j in range(len(ySamples[1])):
1237 yMean=numpy.append(yMean,numpy.mean([ySamples[0,j],ySamples[1,j],ySamples[2,j]]))
1238
1021
1239 '''******* Getting fitting Gaussian ******'''
1240 meanGauss=sum(xSamples*yMean) / len(xSamples)
1241 sigma=sum(yMean*(xSamples-meanGauss)**2) / len(xSamples)
1242
1022
1243 print('****************************')
1023 '''***Fit Gauss CSPC01***'''
1244 print('len(xSamples): ',len(xSamples))
1024 if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3 :
1245 print('yMean: ', yMean.shape, yMean[0:20])
1025 try:
1246 print('ySamples', ySamples.shape, ySamples[0,0:20])
1026 popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
1247 print('xSamples: ',xSamples.shape, xSamples[0:20])
1027 popt02,pcov = curve_fit(self.gaus,xSamples[mask02],numpy.abs(CSPCMask02),p0=CSPCmoments[1])
1028 popt12,pcov = curve_fit(self.gaus,xSamples[mask12],numpy.abs(CSPCMask12),p0=CSPCmoments[2])
1029 FitGauss01 = self.gaus(xSamples,*popt01)
1030 FitGauss02 = self.gaus(xSamples,*popt02)
1031 FitGauss12 = self.gaus(xSamples,*popt12)
1032 except:
1033 FitGauss01=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[0]))
1034 FitGauss02=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[1]))
1035 FitGauss12=numpy.ones(len(xSamples))*numpy.mean(numpy.abs(CSPCSamples[2]))
1036
1037
1038 CSPCopt = numpy.vstack([popt01,popt02,popt12])
1039
1040 '''****** Getting fij width ******'''
1041
1042 yMean = numpy.average(ySamples, axis=0) # ySamples[0]
1043
1044 '''******* Getting fitting Gaussian *******'''
1045 meanGauss = sum(xSamples*yMean) / len(xSamples) # Mu, velocidad radial (frecuencia)
1046 sigma2 = sum(yMean*(xSamples-meanGauss)**2) / len(xSamples) # Varianza, Ancho espectral (frecuencia)
1248
1047
1249 print('meanGauss',meanGauss)
1048 yMoments = self.Moments(yMean, xSamples)
1250 print('sigma',sigma)
1251
1049
1252 #if (abs(meanGauss/sigma**2) > 0.0001) : #0.000000001):
1050 if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3: # and abs(meanGauss/sigma2) > 0.00001:
1253 if dbSNR > SNRlimit :
1051 try:
1254 try:
1052 popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=yMoments)
1255 popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=[1,meanGauss,sigma])
1053 FitGauss=self.gaus(xSamples,*popt)
1256
1054
1257 if numpy.amax(popt)>numpy.amax(yMean)*0.3:
1258 FitGauss=self.gaus(xSamples,*popt)
1259
1260 else:
1261 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1262 print('Verificador: Dentro', Height)
1263 except :#RuntimeError:
1055 except :#RuntimeError:
1264 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1056 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1265
1057
1266
1058
1267 else:
1059 else:
1268 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1060 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1269
1061
1270 Maximun=numpy.amax(yMean)
1271 eMinus1=Maximun*numpy.exp(-1)#*0.8
1272
1062
1273 HWpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1)))
1274 HalfWidth= xFrec[HWpos]
1275 GCpos=self.Find(FitGauss, numpy.amax(FitGauss))
1276 Vpos=self.Find(FactNorm, numpy.amax(FactNorm))
1277
1278 #Vpos=FirstMoment[]
1279
1063
1280 '''****** Getting Fij ******'''
1064 '''****** Getting Fij ******'''
1065 Fijcspc = CSPCopt[:,2]/2*3
1066
1067
1068 GaussCenter = popt[1] #xFrec[GCpos]
1069 #Punto en Eje X de la Gaussiana donde se encuentra el centro
1070 ClosestCenter = xSamples[numpy.abs(xSamples-GaussCenter).argmin()]
1071 PointGauCenter = numpy.where(xSamples==ClosestCenter)[0][0]
1072
1073 #Punto e^-1 hubicado en la Gaussiana
1074 PeMinus1 = numpy.max(FitGauss)* numpy.exp(-1)
1075 FijClosest = FitGauss[numpy.abs(FitGauss-PeMinus1).argmin()] # El punto mas cercano a "Peminus1" dentro de "FitGauss"
1076 PointFij = numpy.where(FitGauss==FijClosest)[0][0]
1281
1077
1282 GaussCenter=xFrec[GCpos]
1078 if xSamples[PointFij] > xSamples[PointGauCenter]:
1283 if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0):
1079 Fij = xSamples[PointFij] - xSamples[PointGauCenter]
1284 Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001
1080
1285 else:
1081 else:
1286 Fij=abs(GaussCenter-HalfWidth)+0.0000001
1082 Fij = xSamples[PointGauCenter] - xSamples[PointFij]
1083
1084
1085 '''****** Taking frequency ranges from SPCs ******'''
1287
1086
1288 '''****** Getting Frecuency range of significant data ******'''
1289
1087
1290 Rangpos=self.Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10)))
1088 #GaussCenter = popt[1] #Primer momento 01
1089 GauWidth = popt[2] *3/2 #Ancho de banda de Gau01
1090 Range = numpy.empty(2)
1091 Range[0] = GaussCenter - GauWidth
1092 Range[1] = GaussCenter + GauWidth
1093 #Punto en Eje X de la Gaussiana donde se encuentra ancho de banda (min:max)
1094 ClosRangeMin = xSamples[numpy.abs(xSamples-Range[0]).argmin()]
1095 ClosRangeMax = xSamples[numpy.abs(xSamples-Range[1]).argmin()]
1291
1096
1292 if Rangpos<GCpos:
1097 PointRangeMin = numpy.where(xSamples==ClosRangeMin)[0][0]
1293 Range=numpy.array([Rangpos,2*GCpos-Rangpos])
1098 PointRangeMax = numpy.where(xSamples==ClosRangeMax)[0][0]
1294 elif Rangpos< ( len(xFrec)- len(xFrec)*0.1):
1099
1295 Range=numpy.array([2*GCpos-Rangpos,Rangpos])
1100 Range=numpy.array([ PointRangeMin, PointRangeMax ])
1296 else:
1101
1297 Range = numpy.array([0,0])
1102 FrecRange = xFrec[ Range[0] : Range[1] ]
1103 VelRange = xVel[ Range[0] : Range[1] ]
1298
1104
1299 print(' ')
1300 print('GCpos',GCpos, ( len(xFrec)- len(xFrec)*0.1))
1301 print('Rangpos',Rangpos)
1302 print('RANGE: ', Range)
1303 FrecRange=xFrec[Range[0]:Range[1]]
1304
1105
1305 '''****** Getting SCPC Slope ******'''
1106 '''****** Getting SCPC Slope ******'''
1306
1107
1307 for i in range(spc.shape[0]):
1108 for i in range(spc.shape[0]):
1308
1109
1309 if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.5:
1110 if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.3:
1310 PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
1111 PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
1112
1113 '''***********************VelRange******************'''
1114
1115 mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
1311
1116
1312 print('FrecRange', len(FrecRange) , FrecRange)
1313 print('PhaseRange', len(PhaseRange), PhaseRange)
1314 print(' ')
1315 if len(FrecRange) == len(PhaseRange):
1117 if len(FrecRange) == len(PhaseRange):
1316 slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange)
1118 try:
1317 PhaseSlope[i]=slope
1119 slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange[mask], PhaseRange[mask])
1318 PhaseInter[i]=intercept
1120 PhaseSlope[i]=slope
1121 PhaseInter[i]=intercept
1122 except:
1123 PhaseSlope[i]=0
1124 PhaseInter[i]=0
1319 else:
1125 else:
1320 PhaseSlope[i]=0
1126 PhaseSlope[i]=0
1321 PhaseInter[i]=0
1127 PhaseInter[i]=0
1322 else:
1128 else:
1323 PhaseSlope[i]=0
1129 PhaseSlope[i]=0
1324 PhaseInter[i]=0
1130 PhaseInter[i]=0
1325
1131
1132
1326 '''Getting constant C'''
1133 '''Getting constant C'''
1327 cC=(Fij*numpy.pi)**2
1134 cC=(Fij*numpy.pi)**2
1328
1135
1329 '''****** Getting constants F and G ******'''
1136 '''****** Getting constants F and G ******'''
1330 MijEijNij=numpy.array([[E02,N02], [E12,N12]])
1137 MijEijNij=numpy.array([[Xi02,Eta02], [Xi12,Eta12]])
1331 MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
1138 MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
1332 MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
1139 MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
1333 MijResults=numpy.array([MijResult0,MijResult1])
1140 MijResults=numpy.array([MijResult0,MijResult1])
1334 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1141 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1335
1142
1336 '''****** Getting constants A, B and H ******'''
1143 '''****** Getting constants A, B and H ******'''
1337 W01=numpy.amax(coherence[0])
1144 W01=numpy.nanmax( FitGauss01 ) #numpy.abs(CSPCSamples[0]))
1338 W02=numpy.amax(coherence[1])
1145 W02=numpy.nanmax( FitGauss02 ) #numpy.abs(CSPCSamples[1]))
1339 W12=numpy.amax(coherence[2])
1146 W12=numpy.nanmax( FitGauss12 ) #numpy.abs(CSPCSamples[2]))
1340
1147
1341 WijResult0=((cF*E01+cG*N01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
1148 WijResult0=((cF*Xi01+cG*Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
1342 WijResult1=((cF*E02+cG*N02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
1149 WijResult1=((cF*Xi02+cG*Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
1343 WijResult2=((cF*E12+cG*N12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
1150 WijResult2=((cF*Xi12+cG*Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
1344
1151
1345 WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
1152 WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
1346
1153
1347 WijEijNij=numpy.array([ [E01**2, N01**2, 2*E01*N01] , [E02**2, N02**2, 2*E02*N02] , [E12**2, N12**2, 2*E12*N12] ])
1154 WijEijNij=numpy.array([ [Xi01**2, Eta01**2, 2*Xi01*Eta01] , [Xi02**2, Eta02**2, 2*Xi02*Eta02] , [Xi12**2, Eta12**2, 2*Xi12*Eta12] ])
1348 (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1155 (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1349
1156
1350 VxVy=numpy.array([[cA,cH],[cH,cB]])
1157 VxVy=numpy.array([[cA,cH],[cH,cB]])
1351
1352 VxVyResults=numpy.array([-cF,-cG])
1158 VxVyResults=numpy.array([-cF,-cG])
1353 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1159 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1354
1160
@@ -1356,10 +1162,15 class FullSpectralAnalysis(Operation):
1356 Vmer = Vx
1162 Vmer = Vx
1357 Vmag=numpy.sqrt(Vzon**2+Vmer**2)
1163 Vmag=numpy.sqrt(Vzon**2+Vmer**2)
1358 Vang=numpy.arctan2(Vmer,Vzon)
1164 Vang=numpy.arctan2(Vmer,Vzon)
1359 Vver=xFrec[Vpos]
1165 if numpy.abs( popt[1] ) < 3.5 and len(FrecRange)>4:
1360 print('vzon y vmer', Vzon, Vmer)
1166 Vver=popt[1]
1361 return Vzon, Vmer, Vver, GaussCenter
1167 else:
1362
1168 Vver=numpy.NaN
1169 FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12])
1170
1171
1172 return Vzon, Vmer, Vver, GaussCenter, PhaseSlope, FitGaussCSPC
1173
1363 class SpectralMoments(Operation):
1174 class SpectralMoments(Operation):
1364
1175
1365 '''
1176 '''
@@ -1384,7 +1195,7 class SpectralMoments(Operation):
1384 self.dataOut.noise : Noise level per channel
1195 self.dataOut.noise : Noise level per channel
1385
1196
1386 Affected:
1197 Affected:
1387 self.dataOut.data_param : Parameters per channel
1198 self.dataOut.moments : Parameters per channel
1388 self.dataOut.data_SNR : SNR per channel
1199 self.dataOut.data_SNR : SNR per channel
1389
1200
1390 '''
1201 '''
@@ -1401,7 +1212,7 class SpectralMoments(Operation):
1401 for ind in range(nChannel):
1212 for ind in range(nChannel):
1402 data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] )
1213 data_param[ind,:,:] = self.__calculateMoments( data[ind,:,:] , absc , noise[ind] )
1403
1214
1404 dataOut.data_param = data_param[:,1:,:]
1215 dataOut.moments = data_param[:,1:,:]
1405 dataOut.data_SNR = data_param[:,0]
1216 dataOut.data_SNR = data_param[:,0]
1406 dataOut.data_DOP = data_param[:,1]
1217 dataOut.data_DOP = data_param[:,1]
1407 dataOut.data_MEAN = data_param[:,2]
1218 dataOut.data_MEAN = data_param[:,2]
@@ -1431,6 +1242,8 class SpectralMoments(Operation):
1431 vec_fd = numpy.zeros(oldspec.shape[1])
1242 vec_fd = numpy.zeros(oldspec.shape[1])
1432 vec_w = numpy.zeros(oldspec.shape[1])
1243 vec_w = numpy.zeros(oldspec.shape[1])
1433 vec_snr = numpy.zeros(oldspec.shape[1])
1244 vec_snr = numpy.zeros(oldspec.shape[1])
1245
1246 oldspec = numpy.ma.masked_invalid(oldspec)
1434
1247
1435 for ind in range(oldspec.shape[1]):
1248 for ind in range(oldspec.shape[1]):
1436
1249
@@ -1469,7 +1282,7 class SpectralMoments(Operation):
1469 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
1282 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
1470 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
1283 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
1471 snr = (spec2.mean()-n0)/n0
1284 snr = (spec2.mean()-n0)/n0
1472
1285
1473 if (snr < 1.e-20) :
1286 if (snr < 1.e-20) :
1474 snr = 1.e-20
1287 snr = 1.e-20
1475
1288
@@ -1477,7 +1290,7 class SpectralMoments(Operation):
1477 vec_fd[ind] = fd
1290 vec_fd[ind] = fd
1478 vec_w[ind] = w
1291 vec_w[ind] = w
1479 vec_snr[ind] = snr
1292 vec_snr[ind] = snr
1480
1293
1481 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
1294 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
1482 return moments
1295 return moments
1483
1296
@@ -1675,7 +1488,6 class SpectralFitting(Operation):
1675 dataCross = dataCross**2/K
1488 dataCross = dataCross**2/K
1676
1489
1677 for h in range(nHeights):
1490 for h in range(nHeights):
1678 # print self.dataOut.heightList[h]
1679
1491
1680 #Input
1492 #Input
1681 d = data[:,h]
1493 d = data[:,h]
@@ -1734,7 +1546,7 class SpectralFitting(Operation):
1734
1546
1735 fm = self.dataOut.library.modelFunction(p, constants)
1547 fm = self.dataOut.library.modelFunction(p, constants)
1736 fmp=numpy.dot(LT,fm)
1548 fmp=numpy.dot(LT,fm)
1737
1549
1738 return dp-fmp
1550 return dp-fmp
1739
1551
1740 def __getSNR(self, z, noise):
1552 def __getSNR(self, z, noise):
@@ -1768,8 +1580,8 class WindProfiler(Operation):
1768
1580
1769 n = None
1581 n = None
1770
1582
1771 def __init__(self, **kwargs):
1583 def __init__(self):
1772 Operation.__init__(self, **kwargs)
1584 Operation.__init__(self)
1773
1585
1774 def __calculateCosDir(self, elev, azim):
1586 def __calculateCosDir(self, elev, azim):
1775 zen = (90 - elev)*numpy.pi/180
1587 zen = (90 - elev)*numpy.pi/180
@@ -2071,12 +1883,9 class WindProfiler(Operation):
2071
1883
2072 Parameters affected: Winds
1884 Parameters affected: Winds
2073 '''
1885 '''
2074 # print arrayMeteor.shape
2075 #Settings
1886 #Settings
2076 nInt = (heightMax - heightMin)/2
1887 nInt = (heightMax - heightMin)/2
2077 # print nInt
2078 nInt = int(nInt)
1888 nInt = int(nInt)
2079 # print nInt
2080 winds = numpy.zeros((2,nInt))*numpy.nan
1889 winds = numpy.zeros((2,nInt))*numpy.nan
2081
1890
2082 #Filter errors
1891 #Filter errors
@@ -2475,8 +2284,8 class WindProfiler(Operation):
2475
2284
2476 class EWDriftsEstimation(Operation):
2285 class EWDriftsEstimation(Operation):
2477
2286
2478 def __init__(self, **kwargs):
2287 def __init__(self):
2479 Operation.__init__(self, **kwargs)
2288 Operation.__init__(self)
2480
2289
2481 def __correctValues(self, heiRang, phi, velRadial, SNR):
2290 def __correctValues(self, heiRang, phi, velRadial, SNR):
2482 listPhi = phi.tolist()
2291 listPhi = phi.tolist()
Show file before | Show file after | Hide comments
@@ -159,9 +159,7 class SpectraProc(ProcessingUnit):
159 dtype='complex')
159 dtype='complex')
160
160
161 if self.dataIn.flagDataAsBlock:
161 if self.dataIn.flagDataAsBlock:
162 # data dimension: [nChannels, nProfiles, nSamples]
163 nVoltProfiles = self.dataIn.data.shape[1]
162 nVoltProfiles = self.dataIn.data.shape[1]
164 # nVoltProfiles = self.dataIn.nProfiles
165
163
166 if nVoltProfiles == nProfiles:
164 if nVoltProfiles == nProfiles:
167 self.buffer = self.dataIn.data.copy()
165 self.buffer = self.dataIn.data.copy()
@@ -299,7 +297,57 class SpectraProc(ProcessingUnit):
299 self.__selectPairsByChannel(self.dataOut.channelList)
297 self.__selectPairsByChannel(self.dataOut.channelList)
300
298
301 return 1
299 return 1
300
301
302 def selectFFTs(self, minFFT, maxFFT ):
303 """
304 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
305 minFFT<= FFT <= maxFFT
306 """
307
308 if (minFFT > maxFFT):
309 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
310
311 if (minFFT < self.dataOut.getFreqRange()[0]):
312 minFFT = self.dataOut.getFreqRange()[0]
313
314 if (maxFFT > self.dataOut.getFreqRange()[-1]):
315 maxFFT = self.dataOut.getFreqRange()[-1]
316
317 minIndex = 0
318 maxIndex = 0
319 FFTs = self.dataOut.getFreqRange()
320
321 inda = numpy.where(FFTs >= minFFT)
322 indb = numpy.where(FFTs <= maxFFT)
323
324 try:
325 minIndex = inda[0][0]
326 except:
327 minIndex = 0
328
329 try:
330 maxIndex = indb[0][-1]
331 except:
332 maxIndex = len(FFTs)
333
334 self.selectFFTsByIndex(minIndex, maxIndex)
302
335
336 return 1
337
338
339 def setH0(self, h0, deltaHeight = None):
340
341 if not deltaHeight:
342 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
343
344 nHeights = self.dataOut.nHeights
345
346 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
347
348 self.dataOut.heightList = newHeiRange
349
350
303 def selectHeights(self, minHei, maxHei):
351 def selectHeights(self, minHei, maxHei):
304 """
352 """
305 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
353 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
@@ -316,9 +364,9 class SpectraProc(ProcessingUnit):
316 1 si el metodo se ejecuto con exito caso contrario devuelve 0
364 1 si el metodo se ejecuto con exito caso contrario devuelve 0
317 """
365 """
318
366
367
319 if (minHei > maxHei):
368 if (minHei > maxHei):
320 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (
369 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei))
321 minHei, maxHei))
322
370
323 if (minHei < self.dataOut.heightList[0]):
371 if (minHei < self.dataOut.heightList[0]):
324 minHei = self.dataOut.heightList[0]
372 minHei = self.dataOut.heightList[0]
@@ -344,6 +392,7 class SpectraProc(ProcessingUnit):
344 maxIndex = len(heights)
392 maxIndex = len(heights)
345
393
346 self.selectHeightsByIndex(minIndex, maxIndex)
394 self.selectHeightsByIndex(minIndex, maxIndex)
395
347
396
348 return 1
397 return 1
349
398
@@ -389,6 +438,40 class SpectraProc(ProcessingUnit):
389
438
390 return 1
439 return 1
391
440
441 def selectFFTsByIndex(self, minIndex, maxIndex):
442 """
443
444 """
445
446 if (minIndex < 0) or (minIndex > maxIndex):
447 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
448
449 if (maxIndex >= self.dataOut.nProfiles):
450 maxIndex = self.dataOut.nProfiles-1
451
452 #Spectra
453 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
454
455 data_cspc = None
456 if self.dataOut.data_cspc is not None:
457 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
458
459 data_dc = None
460 if self.dataOut.data_dc is not None:
461 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
462
463 self.dataOut.data_spc = data_spc
464 self.dataOut.data_cspc = data_cspc
465 self.dataOut.data_dc = data_dc
466
467 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
468 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
469 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
470
471 return 1
472
473
474
392 def selectHeightsByIndex(self, minIndex, maxIndex):
475 def selectHeightsByIndex(self, minIndex, maxIndex):
393 """
476 """
394 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
477 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
@@ -494,7 +577,32 class SpectraProc(ProcessingUnit):
494
577
495 return 1
578 return 1
496
579
497 def removeInterference(self, interf=2, hei_interf=None, nhei_interf=None, offhei_interf=None):
580 def removeInterference2(self):
581
582 cspc = self.dataOut.data_cspc
583 spc = self.dataOut.data_spc
584 Heights = numpy.arange(cspc.shape[2])
585 realCspc = numpy.abs(cspc)
586
587 for i in range(cspc.shape[0]):
588 LinePower= numpy.sum(realCspc[i], axis=0)
589 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
590 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
591 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
592 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
593 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
594
595
596 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
597 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
598 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
599 cspc[i,InterferenceRange,:] = numpy.NaN
600
601
602
603 self.dataOut.data_cspc = cspc
604
605 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
498
606
499 jspectra = self.dataOut.data_spc
607 jspectra = self.dataOut.data_spc
500 jcspectra = self.dataOut.data_cspc
608 jcspectra = self.dataOut.data_cspc
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
NO CONTENT: file was removed
General Comments 0
You need to be logged in to leave comments. Login now