##// 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 112 .vscode
113 113 trash
114 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 190 profileIndex = None
191 191 error = None
192 192 data = None
193 data_plt = None
194
193 nmodes = None
195 194
196 195 def __str__(self):
197 196
@@ -462,6 +461,7 class Spectra(JROData):
462 461 ippFactor = None
463 462 profileIndex = 0
464 463 plotting = "spectra"
464
465 465 def __init__(self):
466 466 '''
467 467 Constructor
@@ -554,9 +554,12 class Spectra(JROData):
554 554
555 555 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
556 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 564 def getNPairs(self):
562 565
Show file before | Show file after | Hide comments
@@ -7,14 +7,190 from .plotting_codes import *
7 7 from schainpy.model.proc.jroproc_base import MPDecorator
8 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 188 isConfig = None
13 189 __nsubplots = None
14 190
15 191 WIDTHPROF = None
16 192 HEIGHTPROF = None
17 PREFIX = 'fitgau'
193 PREFIX = 'SpcParam'
18 194
19 195 def __init__(self, **kwargs):
20 196 Figure.__init__(self, **kwargs)
@@ -83,7 +259,7 class FitGauPlot_(Figure):
83 259 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
84 260 server=None, folder=None, username=None, password=None,
85 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 295 # else:
120 296 # factor = normFactor
121 297 if xaxis == "frequency":
122 x = dataOut.spc_range[0]
298 x = dataOut.spcparam_range[0]
123 299 xlabel = "Frequency (kHz)"
124 300
125 301 elif xaxis == "time":
126 x = dataOut.spc_range[1]
302 x = dataOut.spcparam_range[1]
127 303 xlabel = "Time (ms)"
128 304
129 else:
130 x = dataOut.spc_range[2]
305 else:
306 x = dataOut.spcparam_range[2]
131 307 xlabel = "Velocity (m/s)"
132 308
133 ylabel = "Range (Km)"
309 ylabel = "Range (km)"
134 310
135 311 y = dataOut.getHeiRange()
136 312
137 z = dataOut.GauSPC[:,GauSelector,:,:] #GauSelector] #dataOut.data_spc/factor
138 print('GausSPC', z[0,32,10:40])
313 z = dataOut.SPCparam[Selector] /1966080.0#/ dataOut.normFactor#GauSelector] #dataOut.data_spc/factor
139 314 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
140 315 zdB = 10*numpy.log10(z)
141 316
@@ -657,7 +832,7 class WindProfilerPlot_(Figure):
657 832 # tmax = None
658 833
659 834 x = dataOut.getTimeRange1(dataOut.paramInterval)
660 y = dataOut.heightList
835 y = dataOut.heightList
661 836 z = dataOut.data_output.copy()
662 837 nplots = z.shape[0] #Number of wind dimensions estimated
663 838 nplotsw = nplots
@@ -666,13 +841,14 class WindProfilerPlot_(Figure):
666 841 #If there is a SNR function defined
667 842 if dataOut.data_SNR is not None:
668 843 nplots += 1
669 SNR = dataOut.data_SNR
670 SNRavg = numpy.average(SNR, axis=0)
844 SNR = dataOut.data_SNR[0]
845 SNRavg = SNR#numpy.average(SNR, axis=0)
671 846
672 847 SNRdB = 10*numpy.log10(SNR)
673 848 SNRavgdB = 10*numpy.log10(SNRavg)
674 849
675 if SNRthresh == None: SNRthresh = -5.0
850 if SNRthresh == None:
851 SNRthresh = -5.0
676 852 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
677 853
678 854 for i in range(nplotsw):
@@ -741,8 +917,7 class WindProfilerPlot_(Figure):
741 917 axes = self.axesList[i*self.__nsubplots]
742 918
743 919 z1 = z[i,:].reshape((1,-1))*windFactor[i]
744 #z1=numpy.ma.masked_where(z1==0.,z1)
745
920
746 921 axes.pcolorbuffer(x, y, z1,
747 922 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
748 923 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
@@ -792,8 +967,8 class ParametersPlot_(Figure):
792 967 self.isConfig = False
793 968 self.__nsubplots = 1
794 969
795 self.WIDTH = 800
796 self.HEIGHT = 180
970 self.WIDTH = 300
971 self.HEIGHT = 550
797 972 self.WIDTHPROF = 120
798 973 self.HEIGHTPROF = 0
799 974 self.counter_imagwr = 0
@@ -905,7 +1080,7 class ParametersPlot_(Figure):
905 1080 # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
906 1081 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
907 1082 xlabel = ""
908 ylabel = "Range (Km)"
1083 ylabel = "Range (km)"
909 1084
910 1085 update_figfile = False
911 1086
@@ -949,24 +1124,81 class ParametersPlot_(Figure):
949 1124
950 1125 self.setWinTitle(title)
951 1126
952 for i in range(self.nchan):
953 index = channelIndexList[i]
954 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
955 axes = self.axesList[i*self.plotFact]
956 z1 = z[i,:].reshape((1,-1))
957 axes.pcolorbuffer(x, y, z1,
958 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,
960 ticksize=9, cblabel='', cbsize="1%",colormap=colormap)
961
962 if showSNR:
963 title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
964 axes = self.axesList[i*self.plotFact + 1]
965 SNRdB1 = SNRdB[i,:].reshape((1,-1))
966 axes.pcolorbuffer(x, y, SNRdB1,
967 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,
969 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
1127 # for i in range(self.nchan):
1128 # index = channelIndexList[i]
1129 # title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1130 # axes = self.axesList[i*self.plotFact]
1131 # z1 = z[i,:].reshape((1,-1))
1132 # axes.pcolorbuffer(x, y, z1,
1133 # xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
1134 # xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1135 # ticksize=9, cblabel='', cbsize="1%",colormap=colormap)
1136 #
1137 # if showSNR:
1138 # title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1139 # axes = self.axesList[i*self.plotFact + 1]
1140 # SNRdB1 = SNRdB[i,:].reshape((1,-1))
1141 # axes.pcolorbuffer(x, y, SNRdB1,
1142 # xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1143 # xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
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 1204 self.draw()
@@ -1067,9 +1299,8 class Parameters1Plot_(Figure):
1067 1299 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1068 1300 server=None, folder=None, username=None, password=None,
1069 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 1304 Input:
1074 1305 dataOut :
1075 1306 id :
Show file before | Show file after | Hide comments
@@ -42,6 +42,8 class SpectraPlot_(Figure):
42 42
43 43 self.__xfilter_ena = False
44 44 self.__yfilter_ena = False
45
46 self.indice=1
45 47
46 48 def getSubplots(self):
47 49
@@ -139,10 +141,9 class SpectraPlot_(Figure):
139 141 x = dataOut.getVelRange(1)
140 142 xlabel = "Velocity (m/s)"
141 143
142 ylabel = "Range (Km)"
144 ylabel = "Range (km)"
143 145
144 146 y = dataOut.getHeiRange()
145
146 147 z = dataOut.data_spc/factor
147 148 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
148 149 zdB = 10*numpy.log10(z)
@@ -155,6 +156,7 class SpectraPlot_(Figure):
155 156
156 157 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
157 158 title = wintitle + " Spectra"
159
158 160 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
159 161 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
160 162
@@ -223,6 +225,7 class SpectraPlot_(Figure):
223 225 ftp=ftp,
224 226 wr_period=wr_period,
225 227 thisDatetime=thisDatetime)
228
226 229
227 230 return dataOut
228 231 @MPDecorator
@@ -252,6 +255,8 class CrossSpectraPlot_(Figure):
252 255 self.EXP_CODE = None
253 256 self.SUB_EXP_CODE = None
254 257 self.PLOT_POS = None
258
259 self.indice=0
255 260
256 261 def getSubplots(self):
257 262
@@ -396,6 +401,7 class CrossSpectraPlot_(Figure):
396 401 self.isConfig = True
397 402
398 403 self.setWinTitle(title)
404
399 405
400 406 for i in range(self.nplots):
401 407 pair = dataOut.pairsList[pairsIndexList[i]]
@@ -420,7 +426,7 class CrossSpectraPlot_(Figure):
420 426 xlabel=xlabel, ylabel=ylabel, title=title,
421 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 430 coherence = numpy.abs(coherenceComplex)
425 431 # phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi
426 432 phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi
@@ -439,8 +445,6 class CrossSpectraPlot_(Figure):
439 445 xlabel=xlabel, ylabel=ylabel, title=title,
440 446 ticksize=9, colormap=phase_cmap, cblabel='')
441 447
442
443
444 448 self.draw()
445 449
446 450 self.save(figpath=figpath,
@@ -470,7 +474,7 class RTIPlot_(Figure):
470 474 self.__nsubplots = 1
471 475
472 476 self.WIDTH = 800
473 self.HEIGHT = 180
477 self.HEIGHT = 250
474 478 self.WIDTHPROF = 120
475 479 self.HEIGHTPROF = 0
476 480 self.counter_imagwr = 0
@@ -1497,9 +1501,6 class BeaconPhase_(Figure):
1497 1501 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1498 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 1504 if dataOut.beacon_heiIndexList:
1504 1505 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1505 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 434 def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''):
435 435
436 436 ax = iplot.axes
437
438 437 printLabels(ax, xlabel, ylabel, title)
439 438
440 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 111 timezone=0,
112 112 status_value=0,
113 113 **kwargs):
114
115 114 self.path = path
116 115 self.startDate = startDate
117 116 self.endDate = endDate
Show file before | Show file after | Hide comments
@@ -1815,7 +1815,7 class JRODataWriter(JRODataIO):
1815 1815
1816 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 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 63 if attr:
64 64 message += "%s = %s" % ("size", attr) + "\n"
65 65
66 # print message
67
68
69 66 FILE_STRUCTURE = numpy.dtype([ # HEADER 48bytes
70 67 ('FileMgcNumber', '<u4'), # 0x23020100
71 68 # No Of FDT data records in this file (0 or more)
@@ -94,29 +91,6 class FileHeaderBLTR(Header):
94 91
95 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 94 self.FileMgcNumber = hex(header['FileMgcNumber'][0])
121 95 # No Of FDT data records in this file (0 or more)
122 96 self.nFDTdataRecors = int(header['nFDTdataRecors'][0])
@@ -124,8 +98,6 class FileHeaderBLTR(Header):
124 98 self.OffsetStartHeader = int(header['OffsetStartHeader'][0])
125 99 self.SiteName = str(header['SiteName'][0])
126 100
127 # print 'Numero de bloques', self.nFDTdataRecors
128
129 101 if self.size < 48:
130 102 return 0
131 103
@@ -316,36 +288,10 class RecordHeaderBLTR(Header):
316 288 self.OffsetStartHeader = 48
317 289
318 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 291 startFp = open(fp, "rb")
323 # RecCounter=0
324 # Off2StartNxtRec=811248
325 292 OffRHeader = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
326 print(' ')
327 print('puntero Record Header', startFp.tell())
328 print(' ')
329
330 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 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 295 self.RecMgcNumber = hex(header['RecMgcNumber'][0]) # 0x23030001
350 296 self.RecCounter = int(header['RecCounter'][0])
351 297 self.Off2StartNxtRec = int(header['Off2StartNxtRec'][0])
@@ -397,52 +343,9 class RecordHeaderBLTR(Header):
397 343
398 344 self.RHsize = 180 + 20 * self.nChannels
399 345 self.Datasize = self.nProfiles * self.nChannels * self.nHeights * 2 * 4
400 # print 'Datasize',self.Datasize
401 346 endFp = self.OffsetStartHeader + self.RecCounter * self.Off2StartNxtRec
402 347
403 print('==============================================')
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
348
446 349 if OffRHeader > endFp:
447 350 sys.stderr.write(
448 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 440 FileList.append(IndexFile)
538 441 nFiles += 1
539 442
540 # print 'Files2Read'
541 # print 'Existen '+str(nFiles)+' archivos .fdt'
542
543 443 self.filenameList = FileList # List of files from least to largest by names
544 444
545 445 def run(self, **kwargs):
@@ -553,7 +453,6 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
553 453 self.isConfig = True
554 454
555 455 self.getData()
556 # print 'running'
557 456
558 457 def setup(self, path=None,
559 458 startDate=None,
@@ -590,22 +489,19 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
590 489
591 490 if self.flagNoMoreFiles:
592 491 self.dataOut.flagNoData = True
593 print('NoData se vuelve true')
594 492 return 0
595 493
596 494 self.fp = self.path
597 495 self.Files2Read(self.fp)
598 496 self.readFile(self.fp)
599 497 self.dataOut.data_spc = self.data_spc
600 self.dataOut.data_cspc = self.data_cspc
601 self.dataOut.data_output = self.data_output
602
603 print('self.dataOut.data_output', shape(self.dataOut.data_output))
604
605 # self.removeDC()
606 return self.dataOut.data_spc
607
608 def readFile(self, fp):
498 self.dataOut.data_cspc =self.data_cspc
499 self.dataOut.data_output=self.data_output
500
501 return self.dataOut.data_spc
502
503
504 def readFile(self,fp):
609 505 '''
610 506 You must indicate if you are reading in Online or Offline mode and load the
611 507 The parameters for this file reading mode.
@@ -615,23 +511,18 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
615 511 1. Get the BLTR FileHeader.
616 512 2. Start reading the first block.
617 513 '''
618
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
514
622 515 if self.fileSelector < len(self.filenameList):
623 516
624 517 self.fpFile = str(fp) + '/' + \
625 518 str(self.filenameList[self.fileSelector])
626 # print self.fpFile
627 519 fheader = FileHeaderBLTR()
628 520 fheader.FHread(self.fpFile) # Bltr FileHeader Reading
629 521 self.nFDTdataRecors = fheader.nFDTdataRecors
630 522
631 523 self.readBlock() # Block reading
632 524 else:
633 print('readFile FlagNoData becomes true')
634 self.flagNoMoreFiles = True
525 self.flagNoMoreFiles=True
635 526 self.dataOut.flagNoData = True
636 527 return 0
637 528
@@ -658,12 +549,11 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
658 549 2. Fill the buffer with the current block number.
659 550
660 551 '''
661
662 if self.BlockCounter < self.nFDTdataRecors - 2:
663 print(self.nFDTdataRecors, 'CONDICION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
664 if self.ReadMode == 1:
665 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter + 1)
666 elif self.ReadMode == 0:
552
553 if self.BlockCounter < self.nFDTdataRecors-1:
554 if self.ReadMode==1:
555 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter+1)
556 elif self.ReadMode==0:
667 557 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter)
668 558
669 559 rheader.RHread(self.fpFile) # Bltr FileHeader Reading
@@ -683,31 +573,26 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
683 573
684 574 self.nRdPairs = len(self.dataOut.pairsList)
685 575 self.dataOut.nRdPairs = self.nRdPairs
686
687 self.__firstHeigth = rheader.StartRangeSamp
688 self.__deltaHeigth = rheader.SampResolution
689 self.dataOut.heightList = self.__firstHeigth + \
690 numpy.array(list(range(self.nHeights))) * self.__deltaHeigth
691 self.dataOut.channelList = list(range(self.nChannels))
692 self.dataOut.nProfiles = rheader.nProfiles
693 self.dataOut.nIncohInt = rheader.nIncohInt
694 self.dataOut.nCohInt = rheader.nCohInt
695 self.dataOut.ippSeconds = 1 / float(rheader.PRFhz)
696 self.dataOut.PRF = rheader.PRFhz
697 self.dataOut.nFFTPoints = rheader.nProfiles
698 self.dataOut.utctime = rheader.nUtime
699 self.dataOut.timeZone = 0
700 self.dataOut.normFactor = self.dataOut.nProfiles * \
701 self.dataOut.nIncohInt * self.dataOut.nCohInt
702 self.dataOut.outputInterval = self.dataOut.ippSeconds * \
703 self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles
704
705 self.data_output = numpy.ones([3, rheader.nHeights]) * numpy.NaN
706 print('self.data_output', shape(self.data_output))
707 self.dataOut.velocityX = []
708 self.dataOut.velocityY = []
709 self.dataOut.velocityV = []
710
576 self.__firstHeigth=rheader.StartRangeSamp
577 self.__deltaHeigth=rheader.SampResolution
578 self.dataOut.heightList= self.__firstHeigth + numpy.array(range(self.nHeights))*self.__deltaHeigth
579 self.dataOut.channelList = range(self.nChannels)
580 self.dataOut.nProfiles=rheader.nProfiles
581 self.dataOut.nIncohInt=rheader.nIncohInt
582 self.dataOut.nCohInt=rheader.nCohInt
583 self.dataOut.ippSeconds= 1/float(rheader.PRFhz)
584 self.dataOut.PRF=rheader.PRFhz
585 self.dataOut.nFFTPoints=rheader.nProfiles
586 self.dataOut.utctime=rheader.nUtime
587 self.dataOut.timeZone=0
588 self.dataOut.normFactor= self.dataOut.nProfiles*self.dataOut.nIncohInt*self.dataOut.nCohInt
589 self.dataOut.outputInterval= self.dataOut.ippSeconds * self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles
590
591 self.data_output=numpy.ones([3,rheader.nHeights])*numpy.NaN
592 self.dataOut.velocityX=[]
593 self.dataOut.velocityY=[]
594 self.dataOut.velocityV=[]
595
711 596 '''Block Reading, the Block Data is received and Reshape is used to give it
712 597 shape.
713 598 '''
@@ -734,18 +619,17 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
734 619 y = rho * numpy.sin(phi)
735 620 return(x, y)
736 621
737 if self.DualModeIndex == self.ReadMode:
738
739 self.data_fft = numpy.fromfile(
740 startDATA, [('complex', '<c8')], self.nProfiles * self.nChannels * self.nHeights)
741
742 self.data_fft = self.data_fft.astype(numpy.dtype('complex'))
743
744 self.data_block = numpy.reshape(
745 self.data_fft, (self.nHeights, self.nChannels, self.nProfiles))
746
747 self.data_block = numpy.transpose(self.data_block, (1, 2, 0))
748
622 if self.DualModeIndex==self.ReadMode:
623
624 self.data_fft = numpy.fromfile( startDATA, [('complex','<c8')],self.nProfiles*self.nChannels*self.nHeights )
625 self.data_fft = numpy.empty(101376)
626
627 self.data_fft=self.data_fft.astype(numpy.dtype('complex'))
628
629 self.data_block=numpy.reshape(self.data_fft,(self.nHeights, self.nChannels, self.nProfiles ))
630
631 self.data_block = numpy.transpose(self.data_block, (1,2,0))
632
749 633 copy = self.data_block.copy()
750 634 spc = copy * numpy.conjugate(copy)
751 635
@@ -756,18 +640,8 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
756 640
757 641 z = self.data_spc.copy() # /factor
758 642 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
759 #zdB = 10*numpy.log10(z)
760 print(' ')
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)
643 self.dataOut.data_spc=self.data_spc
644 self.noise = self.dataOut.getNoise(ymin_index=80, ymax_index=132)#/factor
771 645
772 646 ySamples = numpy.ones([3, self.nProfiles])
773 647 phase = numpy.ones([3, self.nProfiles])
@@ -778,20 +652,16 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
778 652 PhaseInter = numpy.ones(3)
779 653
780 654 '''****** Getting CrossSpectra ******'''
781 cspc = self.data_block.copy()
782 self.data_cspc = self.data_block.copy()
783
784 xFrec = self.getVelRange(1)
785 VelRange = self.getVelRange(1)
786 self.dataOut.VelRange = VelRange
787 # print ' '
788 # print ' '
789 # print 'xFrec',xFrec
790 # print ' '
791 # print ' '
792 # Height=35
793 for i in range(self.nRdPairs):
794
655 cspc=self.data_block.copy()
656 self.data_cspc=self.data_block.copy()
657
658 xFrec=self.getVelRange(1)
659 VelRange=self.getVelRange(1)
660 self.dataOut.VelRange=VelRange
661
662
663 for i in range(self.nRdPairs):
664
795 665 chan_index0 = self.dataOut.pairsList[i][0]
796 666 chan_index1 = self.dataOut.pairsList[i][1]
797 667
@@ -820,361 +690,8 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
820 690
821 691 self.dataOut.ChanDist = self.ChanDist
822 692
823
824 # for Height in range(self.nHeights):
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
693 self.BlockCounter+=2
694
1177 695 else:
1178 self.fileSelector += 1
1179 self.BlockCounter = 0
1180 print("Next File") No newline at end of file
696 self.fileSelector+=1
697 self.BlockCounter=0
Show file before | Show file after | Hide comments
@@ -179,9 +179,6 class ParamReader(JRODataReader,ProcessingUnit):
179 179 print("[Reading] %d file(s) was(were) found in time range: %s - %s" %(len(filenameList), startTime, endTime))
180 180 print()
181 181
182 # for i in range(len(filenameList)):
183 # print "[Reading] %s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
184
185 182 self.filenameList = filenameList
186 183 self.datetimeList = datetimeList
187 184
@@ -504,20 +501,11 class ParamReader(JRODataReader,ProcessingUnit):
504 501
505 502 def getData(self):
506 503
507 # if self.flagNoMoreFiles:
508 # self.dataOut.flagNoData = True
509 # print 'Process finished'
510 # return 0
511 #
512 504 if self.blockIndex==self.blocksPerFile:
513 505 if not( self.__setNextFileOffline() ):
514 506 self.dataOut.flagNoData = True
515 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 509 self.__setDataOut()
522 510 self.dataOut.flagNoData = False
523 511
@@ -637,7 +625,10 class ParamWriter(Operation):
637 625 dsDict['variable'] = self.dataList[i]
638 626 #--------------------- Conditionals ------------------------
639 627 #There is no data
628
629
640 630 if dataAux is None:
631
641 632 return 0
642 633
643 634 #Not array, just a number
@@ -821,7 +812,7 class ParamWriter(Operation):
821 812 return False
822 813
823 814 def setNextFile(self):
824
815
825 816 ext = self.ext
826 817 path = self.path
827 818 setFile = self.setFile
@@ -1095,7 +1086,6 class ParamWriter(Operation):
1095 1086 return
1096 1087
1097 1088 self.isConfig = True
1098 # self.putMetadata()
1099 1089 self.setNextFile()
1100 1090
1101 1091 self.putData()
Show file before | Show file after | Hide comments
@@ -413,9 +413,7 class SpectraWriter(JRODataWriter, Operation):
413 413
414 414 data_dc = None
415 415
416 # dataOut = None
417
418 def __init__(self):#, **kwargs):
416 def __init__(self):
419 417 """
420 418 Inicializador de la clase SpectraWriter para la escritura de datos de espectros.
421 419
@@ -429,9 +427,7 class SpectraWriter(JRODataWriter, Operation):
429 427 Return: None
430 428 """
431 429
432 Operation.__init__(self)#, **kwargs)
433
434 #self.isConfig = False
430 Operation.__init__(self)
435 431
436 432 self.nTotalBlocks = 0
437 433
@@ -496,7 +492,7 class SpectraWriter(JRODataWriter, Operation):
496 492
497 493
498 494 def writeBlock(self):
499 """
495 """processingHeaderObj
500 496 Escribe el buffer en el file designado
501 497
502 498 Affected:
@@ -519,8 +515,10 class SpectraWriter(JRODataWriter, Operation):
519 515 data.tofile(self.fp)
520 516
521 517 if self.data_cspc is not None:
522 data = numpy.zeros( self.shape_cspc_Buffer, self.dtype )
518
523 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 522 if not self.processingHeaderObj.shif_fft:
525 523 cspc = numpy.roll( cspc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones
526 524 data['real'] = cspc.real
@@ -529,8 +527,9 class SpectraWriter(JRODataWriter, Operation):
529 527 data.tofile(self.fp)
530 528
531 529 if self.data_dc is not None:
532 data = numpy.zeros( self.shape_dc_Buffer, self.dtype )
530
533 531 dc = self.data_dc
532 data = numpy.zeros( numpy.shape(dc), self.dtype )
534 533 data['real'] = dc.real
535 534 data['imag'] = dc.imag
536 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 10 import itertools
11 11 from multiprocessing import Pool, TimeoutError
12 12 from multiprocessing.pool import ThreadPool
13 import types
14 from functools import partial
15 13 import time
16 #from sklearn.cluster import KMeans
17
18 14
19 15 from scipy.optimize import fmin_l_bfgs_b #optimize with bounds on state papameters
20 16 from .jroproc_base import ProcessingUnit, Operation, MPDecorator
@@ -128,6 +124,7 class ParametersProc(ProcessingUnit):
128 124 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
129 125 self.dataOut.spc_noise = self.dataIn.getNoise()
130 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 128 self.dataOut.pairsList = self.dataIn.pairsList
132 129 self.dataOut.groupList = self.dataIn.pairsList
133 130 self.dataOut.flagNoData = False
@@ -136,9 +133,9 class ParametersProc(ProcessingUnit):
136 133 self.dataOut.ChanDist = self.dataIn.ChanDist
137 134 else: self.dataOut.ChanDist = None
138 135
139 if hasattr(self.dataIn, 'VelRange'): #Velocities range
140 self.dataOut.VelRange = self.dataIn.VelRange
141 else: self.dataOut.VelRange = None
136 #if hasattr(self.dataIn, 'VelRange'): #Velocities range
137 # self.dataOut.VelRange = self.dataIn.VelRange
138 #else: self.dataOut.VelRange = None
142 139
143 140 if hasattr(self.dataIn, 'RadarConst'): #Radar Constant
144 141 self.dataOut.RadarConst = self.dataIn.RadarConst
@@ -184,9 +181,112 class ParametersProc(ProcessingUnit):
184 181 def target(tups):
185 182
186 183 obj, args = tups
187 #print 'TARGETTT', obj, args
184
188 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 290 class GaussianFit(Operation):
191 291
192 292 '''
@@ -198,15 +298,15 class GaussianFit(Operation):
198 298 self.dataOut.data_pre : SelfSpectra
199 299
200 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):
205 Operation.__init__(self, **kwargs)
304 def __init__(self):
305 Operation.__init__(self)
206 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 310 """This routine will find a couple of generalized Gaussians to a power spectrum
211 311 input: spc
212 312 output:
@@ -214,31 +314,12 class GaussianFit(Operation):
214 314 """
215 315
216 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 317 self.Num_Hei = self.spc.shape[2]
233 #self.Num_Bin = len(self.spc)
234 318 self.Num_Bin = self.spc.shape[1]
235 319 self.Num_Chn = self.spc.shape[0]
236
237 320 Vrange = dataOut.abscissaList
238 321
239 #print 'self.spc2', numpy.asarray(self.spc).shape
240
241 GauSPC = numpy.empty([2,self.Num_Bin,self.Num_Hei])
322 GauSPC = numpy.empty([self.Num_Chn,self.Num_Bin,self.Num_Hei])
242 323 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
243 324 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
244 325 SPC_ch1[:] = numpy.NaN
@@ -250,272 +331,12 class GaussianFit(Operation):
250 331 noise_ = dataOut.spc_noise[0].copy()
251 332
252 333
253
254 334 pool = Pool(processes=self.Num_Chn)
255 335 args = [(Vrange, Ch, pnoise, noise_, num_intg, SNRlimit) for Ch in range(self.Num_Chn)]
256 336 objs = [self for __ in range(self.Num_Chn)]
257 337 attrs = list(zip(objs, args))
258 338 gauSPC = pool.map(target, attrs)
259 dataOut.GauSPC = numpy.asarray(gauSPC)
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
339 dataOut.SPCparam = numpy.asarray(SPCparam)
519 340
520 341 ''' Parameters:
521 342 1. Amplitude
@@ -524,16 +345,11 class GaussianFit(Operation):
524 345 4. Power
525 346 '''
526 347
527
528 ###############################################################################
529 348 def FitGau(self, X):
530 349
531 350 Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
532 #print 'VARSSSS', ch, pnoise, noise, num_intg
533
534 #print 'HEIGHTS', self.Num_Hei
535
536 GauSPC = []
351
352 SPCparam = []
537 353 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
538 354 SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei])
539 355 SPC_ch1[:] = 0#numpy.NaN
@@ -542,10 +358,6 class GaussianFit(Operation):
542 358
543 359
544 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 363 spc = numpy.asarray(self.spc)[ch,:,ht]
@@ -554,27 +366,26 class GaussianFit(Operation):
554 366 # normalizing spc and noise
555 367 # This part differs from gg1
556 368 spc_norm_max = max(spc)
557 spc = spc / spc_norm_max
558 pnoise = pnoise / spc_norm_max
369 #spc = spc / spc_norm_max
370 pnoise = pnoise #/ spc_norm_max
559 371 #############################################
560
372
561 373 fatspectra=1.0
562 374
563 wnoise = noise_ / spc_norm_max
375 wnoise = noise_ #/ spc_norm_max
564 376 #wnoise,stdv,i_max,index =enoise(spc,num_intg) #noise estimate using Hildebrand Sekhon, only wnoise is used
565 377 #if wnoise>1.1*pnoise: # to be tested later
566 378 # wnoise=pnoise
567 noisebl=wnoise*0.9; noisebh=wnoise*1.1
379 noisebl=wnoise*0.9;
380 noisebh=wnoise*1.1
568 381 spc=spc-wnoise
569 # print 'wnoise', noise_[0], spc_norm_max, wnoise
382
570 383 minx=numpy.argmin(spc)
384 #spcs=spc.copy()
571 385 spcs=numpy.roll(spc,-minx)
572 386 cum=numpy.cumsum(spcs)
573 387 tot_noise=wnoise * self.Num_Bin #64;
574 #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
575 #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
576 #snr=tot_signal/tot_noise
577 #snr=cum[-1]/tot_noise
388
578 389 snr = sum(spcs)/tot_noise
579 390 snrdB=10.*numpy.log10(snr)
580 391
@@ -582,16 +393,15 class GaussianFit(Operation):
582 393 snr = numpy.NaN
583 394 SPC_ch1[:,ht] = 0#numpy.NaN
584 395 SPC_ch1[:,ht] = 0#numpy.NaN
585 GauSPC = (SPC_ch1,SPC_ch2)
396 SPCparam = (SPC_ch1,SPC_ch2)
586 397 continue
587 #print 'snr',snrdB #, sum(spcs) , tot_noise
588
589 398
590 399
591 400 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
592 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 405 cumlo=cummax*epsi;
596 406 cumhi=cummax*(1-epsi)
597 407 powerindex=numpy.array(numpy.where(numpy.logical_and(cum>cumlo, cum<cumhi))[0])
@@ -613,7 +423,7 class GaussianFit(Operation):
613 423 x=numpy.arange( self.Num_Bin )
614 424 y_data=spc+wnoise
615 425
616 # single gaussian
426 ''' single Gaussian '''
617 427 shift0=numpy.mod(midpeak+minx, self.Num_Bin )
618 428 width0=powerwidth/4.#Initialization entire power of spectrum divided by 4
619 429 power0=2.
@@ -623,16 +433,7 class GaussianFit(Operation):
623 433 lsq1=fmin_l_bfgs_b(self.misfit1,state0,args=(y_data,x,num_intg),bounds=bnds,approx_grad=True)
624 434
625 435 chiSq1=lsq1[1];
626 jack1= self.y_jacobian1(x,lsq1[0])
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
436
636 437
637 438 if fatspectra<1.0 and powerwidth<4:
638 439 choice=0
@@ -648,7 +449,7 class GaussianFit(Operation):
648 449 #return (numpy.array([shift0,width0,Amplitude0,p0]),
649 450 # numpy.array([shift1,width1,Amplitude1,p1]),noise,snrdB,chiSq1,6.,sigmas1,[None,]*9,choice)
650 451
651 # two gaussians
452 ''' two gaussians '''
652 453 #shift0=numpy.mod(firstpeak+minx,64); shift1=numpy.mod(secondpeak+minx,64)
653 454 shift0=numpy.mod(firstpeak+minx, self.Num_Bin );
654 455 shift1=numpy.mod(secondpeak+minx, self.Num_Bin )
@@ -663,24 +464,16 class GaussianFit(Operation):
663 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 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)
667
668
669 chiSq2=lsq2[1]; jack2=self.y_jacobian2(x,lsq2[0])
467 lsq2 = fmin_l_bfgs_b( self.misfit2 , state0 , args=(y_data,x,num_intg) , bounds=bnds , approx_grad=True )
468
469
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 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
683 if snrdB>-6: # when SNR is strong pick the peak with least shift (LOS velocity) error
475
476 if snrdB>-12: # when SNR is strong pick the peak with least shift (LOS velocity) error
684 477 if oneG:
685 478 choice=0
686 479 else:
@@ -690,7 +483,7 class GaussianFit(Operation):
690 483 s1=(2**(1+1./p1))*scipy.special.gamma(1./p1)/p1;
691 484 s2=(2**(1+1./p2))*scipy.special.gamma(1./p2)/p2;
692 485 gp1=a1*w1*s1; gp2=a2*w2*s2 # power content of each ggaussian with proper p scaling
693
486
694 487 if gp1>gp2:
695 488 if a1>0.7*a2:
696 489 choice=1
@@ -710,13 +503,15 class GaussianFit(Operation):
710 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])
714 shift1=lsq2[0][4]; vel1=Vrange[0] + shift1*(Vrange[1]-Vrange[0])
506 shift0=lsq2[0][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 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 515 shift0=lsq2[0][0]
721 516 width0=lsq2[0][1]
722 517 Amplitude0=lsq2[0][2]
@@ -739,120 +534,18 class GaussianFit(Operation):
739 534 p0=lsq2[0][7]
740 535 noise=lsq2[0][8]
741 536
742 if Amplitude0<0.1: # in case the peak is noise
743 shift0,width0,Amplitude0,p0 = 4*[numpy.NaN]
744 if Amplitude1<0.1:
745 shift1,width1,Amplitude1,p1 = 4*[numpy.NaN]
746
747
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)
537 if Amplitude0<0.05: # in case the peak is noise
538 shift0,width0,Amplitude0,p0 = [0,0,0,0]#4*[numpy.NaN]
539 if Amplitude1<0.05:
540 shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN]
541
542
780 543 SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
781 544 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
782 #print 'SPC_ch1.shape',SPC_ch1.shape
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
545 SPCparam = (SPC_ch1,SPC_ch2)
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;
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
548 return GauSPC
856 549
857 550 def y_model1(self,x,state):
858 551 shift0,width0,amplitude0,power0,noise=state
@@ -884,6 +577,7 class GaussianFit(Operation):
884 577 def misfit2(self,state,y_data,x,num_intg):
885 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 582 class PrecipitationProc(Operation):
889 583
@@ -900,24 +594,61 class PrecipitationProc(Operation):
900 594
901 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,
906 tauW=None, ThetaT=None, ThetaR=None, Km = 0.93, Altitude=None):
598 def __init__(self):
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 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 634 Ze = self.dBZeMODE2(dataOut)
918 635
919 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 652 self.Pt = Pt
922 653 self.Gt = Gt
923 654 self.Gr = Gr
@@ -927,48 +658,101 class PrecipitationProc(Operation):
927 658 self.ThetaT = ThetaT
928 659 self.ThetaR = ThetaR
929 660
930 RadarConstant = GetRadarConstant()
931 SPCmean = numpy.mean(self.spc,0)
932 ETA = numpy.zeros(self.Num_Hei)
933 Pr = numpy.sum(SPCmean,0)
661 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
662 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
663 RadarConstant = 5e-26 * Numerator / Denominator #
934 664
935 #for R in range(self.Num_Hei):
936 # ETA[R] = RadarConstant * Pr[R] * R**2 #Reflectivity (ETA)
937
938 D_range = numpy.zeros(self.Num_Hei)
939 EqSec = numpy.zeros(self.Num_Hei)
665 ''' ============================= '''
666
667 self.spc[0] = (self.spc[0]-dataOut.noise[0])
668 self.spc[1] = (self.spc[1]-dataOut.noise[1])
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 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 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 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
949 SIGMA[R] = numpy.pi**5 / Lambda**4 * Km * D_range[R]**6 #Equivalent Section of drops (sigma)
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
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 742 dBRR = 10*numpy.log10(RR)
743 dBRR2 = 10*numpy.log10(RR2)
960 744
961 745 dBZe = 10*numpy.log10(Ze)
962 dataOut.data_output = Ze
963 dataOut.data_param = numpy.ones([2,self.Num_Hei])
964 dataOut.channelList = [0,1]
965 print('channelList', dataOut.channelList)
966 dataOut.data_param[0]=dBZe
967 dataOut.data_param[1]=dBRR
968 print('RR SHAPE', dBRR.shape)
969 print('Ze SHAPE', dBZe.shape)
970 print('dataOut.data_param SHAPE', dataOut.data_param.shape)
971
746 dBZ = 10*numpy.log10(Z)
747
748 dataOut.data_output = RR[8]
749 dataOut.data_param = numpy.ones([3,self.Num_Hei])
750 dataOut.channelList = [0,1,2]
751
752 dataOut.data_param[0]=dBZ
753 dataOut.data_param[1]=V_mean
754 dataOut.data_param[2]=RR
755
972 756
973 757 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
974 758
@@ -983,7 +767,7 class PrecipitationProc(Operation):
983 767 data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN
984 768
985 769 ETA = numpy.sum(SNR,1)
986 print('ETA' , ETA)
770
987 771 ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
988 772
989 773 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
@@ -995,26 +779,27 class PrecipitationProc(Operation):
995 779
996 780 return Ze
997 781
998 def GetRadarConstant(self):
999
1000 """
1001 Constants:
1002
1003 Pt: Transmission Power dB
1004 Gt: Transmission Gain dB
1005 Gr: Reception Gain dB
1006 Lambda: Wavelenght m
1007 aL: Attenuation loses dB
1008 tauW: Width of transmission pulse s
1009 ThetaT: Transmission antenna bean angle rad
1010 ThetaR: Reception antenna beam angle rad
1011
1012 """
1013 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
1014 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
1015 RadarConstant = Numerator / Denominator
1016
1017 return RadarConstant
782 # def GetRadarConstant(self):
783 #
784 # """
785 # Constants:
786 #
787 # Pt: Transmission Power dB 5kW 5000
788 # Gt: Transmission Gain dB 24.7 dB 295.1209
789 # Gr: Reception Gain dB 18.5 dB 70.7945
790 # Lambda: Wavelenght m 0.6741 m 0.6741
791 # aL: Attenuation loses dB 4dB 2.5118
792 # tauW: Width of transmission pulse s 4us 4e-6
793 # ThetaT: Transmission antenna bean angle rad 0.1656317 rad 0.1656317
794 # ThetaR: Reception antenna beam angle rad 0.36774087 rad 0.36774087
795 #
796 # """
797 #
798 # Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
799 # Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * TauW * numpy.pi * ThetaT * TheraR)
800 # RadarConstant = Numerator / Denominator
801 #
802 # return RadarConstant
1018 803
1019 804
1020 805
@@ -1037,10 +822,20 class FullSpectralAnalysis(Operation):
1037 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 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 840 nChannel = spc.shape[0]
1046 841 nProfiles = spc.shape[1]
@@ -1050,14 +845,9 class FullSpectralAnalysis(Operation):
1050 845 if dataOut.ChanDist is not None :
1051 846 ChanDist = dataOut.ChanDist
1052 847 else:
1053 ChanDist = numpy.array([[E01, N01],[E02,N02],[E12,N12]])
1054
1055 #print 'ChanDist', ChanDist
848 ChanDist = numpy.array([[Xi01, Eta01],[Xi02,Eta02],[Xi12,Eta12]])
1056 849
1057 if dataOut.VelRange is not None:
1058 VelRange= dataOut.VelRange
1059 else:
1060 VelRange= dataOut.abscissaList
850 FrecRange = dataOut.spc_range[0]
1061 851
1062 852 ySamples=numpy.ones([nChannel,nProfiles])
1063 853 phase=numpy.ones([nChannel,nProfiles])
@@ -1065,113 +855,108 class FullSpectralAnalysis(Operation):
1065 855 coherence=numpy.ones([nChannel,nProfiles])
1066 856 PhaseSlope=numpy.ones(nChannel)
1067 857 PhaseInter=numpy.ones(nChannel)
1068 dataSNR = dataOut.data_SNR
1069
1070
858 data_SNR=numpy.zeros([nProfiles])
1071 859
1072 860 data = dataOut.data_pre
1073 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)
1078 #SNRdBMean = []
1079
863 dataOut.data_SNR = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0]
1080 864
1081 #for j in range(nHeights):
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]]))
1083 # SNRdBMean = numpy.append(SNRdBMean,numpy.mean([SNRdB[0,j],SNRdB[1,j],SNRdB[2,j]]))
1084
1085 data_output=numpy.ones([3,spc.shape[2]])*numpy.NaN
865 dataOut.data_SNR[numpy.where( dataOut.data_SNR <0 )] = 1e-20
866
867
868 data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN
1086 869
1087 870 velocityX=[]
1088 871 velocityY=[]
1089 872 velocityV=[]
873 PhaseLine=[]
1090 874
1091 dbSNR = 10*numpy.log10(dataSNR)
875 dbSNR = 10*numpy.log10(dataOut.data_SNR)
1092 876 dbSNR = numpy.average(dbSNR,0)
877
1093 878 for Height in range(nHeights):
1094
1095 [Vzon,Vmer,Vver, GaussCenter]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, VelRange, dbSNR[Height], SNRlimit)
879
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 883 if abs(Vzon)<100. and abs(Vzon)> 0.:
1098 884 velocityX=numpy.append(velocityX, Vzon)#Vmag
1099 885
1100 886 else:
1101 print('Vzon',Vzon)
1102 887 velocityX=numpy.append(velocityX, numpy.NaN)
1103 888
1104 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 892 else:
1108 print('Vmer',Vmer)
1109 893 velocityY=numpy.append(velocityY, numpy.NaN)
1110 894
1111 895 if dbSNR[Height] > SNRlimit:
1112 velocityV=numpy.append(velocityV, FirstMoment[Height])
896 velocityV=numpy.append(velocityV, -Vver)#FirstMoment[Height])
1113 897 else:
1114 898 velocityV=numpy.append(velocityV, numpy.NaN)
1115 #FirstMoment[Height]= numpy.NaN
1116 # if SNRdBMean[Height] <12:
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(' ')
899
900
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 909 dataOut.data_output=data_output
910
1141 911 return
1142 912
1143 913
1144 914 def moving_average(self,x, N=2):
1145 915 return numpy.convolve(x, numpy.ones((N,))/N)[(N-1):]
1146 916
1147 def gaus(self,xSamples,a,x0,sigma):
1148 return a*numpy.exp(-(xSamples-x0)**2/(2*sigma**2))
917 def gaus(self,xSamples,Amp,Mu,Sigma):
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):
1151 for index in range(len(x)):
1152 if x[index]==value:
1153 return index
922 def Moments(self, ySamples, xSamples):
923 Pot = numpy.nansum( ySamples ) # Potencia, momento 0
924 yNorm = ySamples / Pot
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 935 ySamples=numpy.ones([spc.shape[0],spc.shape[1]])
1158 936 phase=numpy.ones([spc.shape[0],spc.shape[1]])
1159 937 CSPCSamples=numpy.ones([spc.shape[0],spc.shape[1]],dtype=numpy.complex_)
1160 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 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 950 '''Getting Eij and Nij'''
1166 951
1167 E01=ChanDist[0][0]
1168 N01=ChanDist[0][1]
952 Xi01=ChanDist[0][0]
953 Eta01=ChanDist[0][1]
1169 954
1170 E02=ChanDist[1][0]
1171 N02=ChanDist[1][1]
955 Xi02=ChanDist[1][0]
956 Eta02=ChanDist[1][1]
1172 957
1173 E12=ChanDist[2][0]
1174 N12=ChanDist[2][1]
958 Xi12=ChanDist[2][0]
959 Eta12=ChanDist[2][1]
1175 960
1176 961 z = spc.copy()
1177 962 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
@@ -1179,176 +964,197 class FullSpectralAnalysis(Operation):
1179 964 for i in range(spc.shape[0]):
1180 965
1181 966 '''****** Line of Data SPC ******'''
1182 zline=z[i,:,Height]
967 zline=z[i,:,Height].copy() - noise[i] # Se resta ruido
1183 968
1184 969 '''****** SPC is normalized ******'''
1185 FactNorm= (zline.copy()-noise[i]) / numpy.sum(zline.copy())
1186 FactNorm= FactNorm/numpy.sum(FactNorm)
970 SmoothSPC =self.moving_average(zline.copy(),N=1) # Se suaviza el ruido
971 FactNorm = SmoothSPC/numpy.nansum(SmoothSPC) # SPC Normalizado y suavizado
1187 972
1188 SmoothSPC=self.moving_average(FactNorm,N=3)
1189
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))
973 xSamples = xFrec # Se toma el rango de frecuncias
974 ySamples[i] = FactNorm # Se toman los valores de SPC normalizado
1204 975
1205 976 for i in range(spc.shape[0]):
1206 977
1207 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 984 chan_index0 = pairsList[i][0]
1212 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 990 CSPCSamples[i] = CSPCNorm
991
1218 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 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])
1225 print('CSPCFactor', CSPCFactor)#, CSPCFactor[0:20]
1226 print(numpy.sum(ySamples[chan_index0]), numpy.sum(ySamples[chan_index1]), -noise[i])
1227 print('CSPCNorm', CSPCNorm.shape, CSPCNorm[0:20])
1228 print('CSPCNorm suma', numpy.sum(CSPCNorm))
1229 print('CSPCSamples', CSPCSamples.shape, CSPCSamples[0,0:20])
998 CSPCmoments = numpy.vstack([self.Moments(numpy.abs(CSPCSamples[0]), xSamples),
999 self.Moments(numpy.abs(CSPCSamples[1]), xSamples),
1000 self.Moments(numpy.abs(CSPCSamples[2]), xSamples)])
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=[]
1234 yMean2=[]
1015 #mask = ~numpy.isnan(CSPCMask01)
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('****************************')
1244 print('len(xSamples): ',len(xSamples))
1245 print('yMean: ', yMean.shape, yMean[0:20])
1246 print('ySamples', ySamples.shape, ySamples[0,0:20])
1247 print('xSamples: ',xSamples.shape, xSamples[0:20])
1023 '''***Fit Gauss CSPC01***'''
1024 if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3 :
1025 try:
1026 popt01,pcov = curve_fit(self.gaus,xSamples[mask01],numpy.abs(CSPCMask01),p0=CSPCmoments[0])
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)
1250 print('sigma',sigma)
1048 yMoments = self.Moments(yMean, xSamples)
1251 1049
1252 #if (abs(meanGauss/sigma**2) > 0.0001) : #0.000000001):
1253 if dbSNR > SNRlimit :
1254 try:
1255 popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=[1,meanGauss,sigma])
1050 if dbSNR > SNRlimit and numpy.abs(SPCmoments[1])<3: # and abs(meanGauss/sigma2) > 0.00001:
1051 try:
1052 popt,pcov = curve_fit(self.gaus,xSamples,yMean,p0=yMoments)
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 1055 except :#RuntimeError:
1264 1056 FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
1265 1057
1266
1058
1267 1059 else:
1268 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 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]
1283 if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0):
1284 Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001
1078 if xSamples[PointFij] > xSamples[PointGauCenter]:
1079 Fij = xSamples[PointFij] - xSamples[PointGauCenter]
1080
1285 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:
1293 Range=numpy.array([Rangpos,2*GCpos-Rangpos])
1294 elif Rangpos< ( len(xFrec)- len(xFrec)*0.1):
1295 Range=numpy.array([2*GCpos-Rangpos,Rangpos])
1296 else:
1297 Range = numpy.array([0,0])
1097 PointRangeMin = numpy.where(xSamples==ClosRangeMin)[0][0]
1098 PointRangeMax = numpy.where(xSamples==ClosRangeMax)[0][0]
1099
1100 Range=numpy.array([ PointRangeMin, PointRangeMax ])
1101
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 1106 '''****** Getting SCPC Slope ******'''
1306 1107
1307 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 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 1117 if len(FrecRange) == len(PhaseRange):
1316 slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange)
1317 PhaseSlope[i]=slope
1318 PhaseInter[i]=intercept
1118 try:
1119 slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange[mask], PhaseRange[mask])
1120 PhaseSlope[i]=slope
1121 PhaseInter[i]=intercept
1122 except:
1123 PhaseSlope[i]=0
1124 PhaseInter[i]=0
1319 1125 else:
1320 1126 PhaseSlope[i]=0
1321 1127 PhaseInter[i]=0
1322 1128 else:
1323 1129 PhaseSlope[i]=0
1324 1130 PhaseInter[i]=0
1325
1131
1132
1326 1133 '''Getting constant C'''
1327 1134 cC=(Fij*numpy.pi)**2
1328 1135
1329 1136 '''****** Getting constants F and G ******'''
1330 MijEijNij=numpy.array([[E02,N02], [E12,N12]])
1137 MijEijNij=numpy.array([[Xi02,Eta02], [Xi12,Eta12]])
1331 1138 MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
1332 1139 MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
1333 1140 MijResults=numpy.array([MijResult0,MijResult1])
1334 1141 (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
1335 1142
1336 1143 '''****** Getting constants A, B and H ******'''
1337 W01=numpy.amax(coherence[0])
1338 W02=numpy.amax(coherence[1])
1339 W12=numpy.amax(coherence[2])
1144 W01=numpy.nanmax( FitGauss01 ) #numpy.abs(CSPCSamples[0]))
1145 W02=numpy.nanmax( FitGauss02 ) #numpy.abs(CSPCSamples[1]))
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))
1342 WijResult1=((cF*E02+cG*N02)**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))
1148 WijResult0=((cF*Xi01+cG*Eta01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
1149 WijResult1=((cF*Xi02+cG*Eta02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
1150 WijResult2=((cF*Xi12+cG*Eta12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
1344 1151
1345 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 1155 (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
1349 1156
1350 1157 VxVy=numpy.array([[cA,cH],[cH,cB]])
1351
1352 1158 VxVyResults=numpy.array([-cF,-cG])
1353 1159 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1354 1160
@@ -1356,10 +1162,15 class FullSpectralAnalysis(Operation):
1356 1162 Vmer = Vx
1357 1163 Vmag=numpy.sqrt(Vzon**2+Vmer**2)
1358 1164 Vang=numpy.arctan2(Vmer,Vzon)
1359 Vver=xFrec[Vpos]
1360 print('vzon y vmer', Vzon, Vmer)
1361 return Vzon, Vmer, Vver, GaussCenter
1362
1165 if numpy.abs( popt[1] ) < 3.5 and len(FrecRange)>4:
1166 Vver=popt[1]
1167 else:
1168 Vver=numpy.NaN
1169 FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12])
1170
1171
1172 return Vzon, Vmer, Vver, GaussCenter, PhaseSlope, FitGaussCSPC
1173
1363 1174 class SpectralMoments(Operation):
1364 1175
1365 1176 '''
@@ -1384,7 +1195,7 class SpectralMoments(Operation):
1384 1195 self.dataOut.noise : Noise level per channel
1385 1196
1386 1197 Affected:
1387 self.dataOut.data_param : Parameters per channel
1198 self.dataOut.moments : Parameters per channel
1388 1199 self.dataOut.data_SNR : SNR per channel
1389 1200
1390 1201 '''
@@ -1401,7 +1212,7 class SpectralMoments(Operation):
1401 1212 for ind in range(nChannel):
1402 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 1216 dataOut.data_SNR = data_param[:,0]
1406 1217 dataOut.data_DOP = data_param[:,1]
1407 1218 dataOut.data_MEAN = data_param[:,2]
@@ -1431,6 +1242,8 class SpectralMoments(Operation):
1431 1242 vec_fd = numpy.zeros(oldspec.shape[1])
1432 1243 vec_w = numpy.zeros(oldspec.shape[1])
1433 1244 vec_snr = numpy.zeros(oldspec.shape[1])
1245
1246 oldspec = numpy.ma.masked_invalid(oldspec)
1434 1247
1435 1248 for ind in range(oldspec.shape[1]):
1436 1249
@@ -1469,7 +1282,7 class SpectralMoments(Operation):
1469 1282 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
1470 1283 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
1471 1284 snr = (spec2.mean()-n0)/n0
1472
1285
1473 1286 if (snr < 1.e-20) :
1474 1287 snr = 1.e-20
1475 1288
@@ -1477,7 +1290,7 class SpectralMoments(Operation):
1477 1290 vec_fd[ind] = fd
1478 1291 vec_w[ind] = w
1479 1292 vec_snr[ind] = snr
1480
1293
1481 1294 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
1482 1295 return moments
1483 1296
@@ -1675,7 +1488,6 class SpectralFitting(Operation):
1675 1488 dataCross = dataCross**2/K
1676 1489
1677 1490 for h in range(nHeights):
1678 # print self.dataOut.heightList[h]
1679 1491
1680 1492 #Input
1681 1493 d = data[:,h]
@@ -1734,7 +1546,7 class SpectralFitting(Operation):
1734 1546
1735 1547 fm = self.dataOut.library.modelFunction(p, constants)
1736 1548 fmp=numpy.dot(LT,fm)
1737
1549
1738 1550 return dp-fmp
1739 1551
1740 1552 def __getSNR(self, z, noise):
@@ -1768,8 +1580,8 class WindProfiler(Operation):
1768 1580
1769 1581 n = None
1770 1582
1771 def __init__(self, **kwargs):
1772 Operation.__init__(self, **kwargs)
1583 def __init__(self):
1584 Operation.__init__(self)
1773 1585
1774 1586 def __calculateCosDir(self, elev, azim):
1775 1587 zen = (90 - elev)*numpy.pi/180
@@ -2071,12 +1883,9 class WindProfiler(Operation):
2071 1883
2072 1884 Parameters affected: Winds
2073 1885 '''
2074 # print arrayMeteor.shape
2075 1886 #Settings
2076 1887 nInt = (heightMax - heightMin)/2
2077 # print nInt
2078 1888 nInt = int(nInt)
2079 # print nInt
2080 1889 winds = numpy.zeros((2,nInt))*numpy.nan
2081 1890
2082 1891 #Filter errors
@@ -2475,8 +2284,8 class WindProfiler(Operation):
2475 2284
2476 2285 class EWDriftsEstimation(Operation):
2477 2286
2478 def __init__(self, **kwargs):
2479 Operation.__init__(self, **kwargs)
2287 def __init__(self):
2288 Operation.__init__(self)
2480 2289
2481 2290 def __correctValues(self, heiRang, phi, velRadial, SNR):
2482 2291 listPhi = phi.tolist()
Show file before | Show file after | Hide comments
@@ -159,9 +159,7 class SpectraProc(ProcessingUnit):
159 159 dtype='complex')
160 160
161 161 if self.dataIn.flagDataAsBlock:
162 # data dimension: [nChannels, nProfiles, nSamples]
163 162 nVoltProfiles = self.dataIn.data.shape[1]
164 # nVoltProfiles = self.dataIn.nProfiles
165 163
166 164 if nVoltProfiles == nProfiles:
167 165 self.buffer = self.dataIn.data.copy()
@@ -299,7 +297,57 class SpectraProc(ProcessingUnit):
299 297 self.__selectPairsByChannel(self.dataOut.channelList)
300 298
301 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 351 def selectHeights(self, minHei, maxHei):
304 352 """
305 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 364 1 si el metodo se ejecuto con exito caso contrario devuelve 0
317 365 """
318 366
367
319 368 if (minHei > maxHei):
320 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (
321 minHei, maxHei))
369 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei))
322 370
323 371 if (minHei < self.dataOut.heightList[0]):
324 372 minHei = self.dataOut.heightList[0]
@@ -344,6 +392,7 class SpectraProc(ProcessingUnit):
344 392 maxIndex = len(heights)
345 393
346 394 self.selectHeightsByIndex(minIndex, maxIndex)
395
347 396
348 397 return 1
349 398
@@ -389,6 +438,40 class SpectraProc(ProcessingUnit):
389 438
390 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 475 def selectHeightsByIndex(self, minIndex, maxIndex):
393 476 """
394 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 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 607 jspectra = self.dataOut.data_spc
500 608 jcspectra = self.dataOut.data_cspc
Show file before | Show file after | Hide comments
1 NO CONTENT: file was removed
General Comments 0
You need to be logged in to leave comments. Login now