##// END OF EJS Templates
schain-jc
RSS
Last Commit foreva!
ebocanegra -
r1157:823d4012bd34
parent child
Show More
Side by Side Unified
Context file:
r1157:823d4012bd34 Loading diff...:
Show file before | Show file after | Hide comments
@@ -5,6 +5,181 import inspect
5 5 from figure import Figure, isRealtime, isTimeInHourRange
6 6 from plotting_codes import *
7 7
8 class ParamLine(Figure):
9
10 isConfig = None
11
12 def __init__(self):
13
14 self.isConfig = False
15 self.WIDTH = 300
16 self.HEIGHT = 200
17 self.counter_imagwr = 0
18
19 def getSubplots(self):
20
21 nrow = self.nplots
22 ncol = 3
23 return nrow, ncol
24
25 def setup(self, id, nplots, wintitle, show):
26
27 self.nplots = nplots
28
29 self.createFigure(id=id,
30 wintitle=wintitle,
31 show=show)
32
33 nrow,ncol = self.getSubplots()
34 colspan = 3
35 rowspan = 1
36
37 for i in range(nplots):
38 self.addAxes(nrow, ncol, i, 0, colspan, rowspan)
39
40 def plot_iq(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax):
41 yreal = y[channelIndexList,:].real
42 yimag = y[channelIndexList,:].imag
43
44 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
45 xlabel = "Range (Km)"
46 ylabel = "Intensity - IQ"
47
48 if not self.isConfig:
49 nplots = len(channelIndexList)
50
51 self.setup(id=id,
52 nplots=nplots,
53 wintitle='',
54 show=show)
55
56 if xmin == None: xmin = numpy.nanmin(x)
57 if xmax == None: xmax = numpy.nanmax(x)
58 if ymin == None: ymin = min(numpy.nanmin(yreal),numpy.nanmin(yimag))
59 if ymax == None: ymax = max(numpy.nanmax(yreal),numpy.nanmax(yimag))
60
61 self.isConfig = True
62
63 self.setWinTitle(title)
64
65 for i in range(len(self.axesList)):
66 title = "Channel %d" %(i)
67 axes = self.axesList[i]
68
69 axes.pline(x, yreal[i,:],
70 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
71 xlabel=xlabel, ylabel=ylabel, title=title)
72
73 axes.addpline(x, yimag[i,:], idline=1, color="red", linestyle="solid", lw=2)
74
75 def plot_power(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax):
76 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
77 yreal = y.real
78
79 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
80 xlabel = "Range (Km)"
81 ylabel = "Intensity"
82
83 if not self.isConfig:
84 nplots = len(channelIndexList)
85
86 self.setup(id=id,
87 nplots=nplots,
88 wintitle='',
89 show=show)
90
91 if xmin == None: xmin = numpy.nanmin(x)
92 if xmax == None: xmax = numpy.nanmax(x)
93 if ymin == None: ymin = numpy.nanmin(yreal)
94 if ymax == None: ymax = numpy.nanmax(yreal)
95
96 self.isConfig = True
97
98 self.setWinTitle(title)
99
100 for i in range(len(self.axesList)):
101 title = "Channel %d" %(i)
102 axes = self.axesList[i]
103 ychannel = yreal[i,:]
104 axes.pline(x, ychannel,
105 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
106 xlabel=xlabel, ylabel=ylabel, title=title)
107
108
109 def run(self, dataOut, id, wintitle="", channelList=None,
110 xmin=None, xmax=None, ymin=None, ymax=None, save=False,
111 figpath='./', figfile=None, show=True, wr_period=1,
112 ftp=False, server=None, folder=None, username=None, password=None):
113
114 """
115
116 Input:
117 dataOut :
118 id :
119 wintitle :
120 channelList :
121 xmin : None,
122 xmax : None,
123 ymin : None,
124 ymax : None,
125 """
126
127 if channelList == None:
128 channelIndexList = dataOut.channelIndexList
129 else:
130 channelIndexList = []
131 for channel in channelList:
132 if channel not in dataOut.channelList:
133 raise ValueError, "Channel %d is not in dataOut.channelList"
134 channelIndexList.append(dataOut.channelList.index(channel))
135
136 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
137
138 y = dataOut.RR
139
140 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
141 xlabel = "Range (Km)"
142 ylabel = "Intensity"
143
144 if not self.isConfig:
145 nplots = len(channelIndexList)
146
147 self.setup(id=id,
148 nplots=nplots,
149 wintitle='',
150 show=show)
151
152 if xmin == None: xmin = numpy.nanmin(x)
153 if xmax == None: xmax = numpy.nanmax(x)
154 if ymin == None: ymin = numpy.nanmin(y)
155 if ymax == None: ymax = numpy.nanmax(y)
156
157 self.isConfig = True
158
159 self.setWinTitle(title)
160
161 for i in range(len(self.axesList)):
162 title = "Channel %d" %(i)
163 axes = self.axesList[i]
164 ychannel = y[i,:]
165 axes.pline(x, ychannel,
166 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
167 xlabel=xlabel, ylabel=ylabel, title=title)
168
169
170 self.draw()
171
172 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + "_" + str(dataOut.profileIndex)
173 figfile = self.getFilename(name = str_datetime)
174
175 self.save(figpath=figpath,
176 figfile=figfile,
177 save=save,
178 ftp=ftp,
179 wr_period=wr_period,
180 thisDatetime=thisDatetime)
181
182
8 183
9 184 class SpcParamPlot(Figure):
10 185
@@ -800,8 +975,8 class ParametersPlot(Figure):
800 975 self.isConfig = False
801 976 self.__nsubplots = 1
802 977
803 self.WIDTH = 800
804 self.HEIGHT = 250
978 self.WIDTH = 300
979 self.HEIGHT = 550
805 980 self.WIDTHPROF = 120
806 981 self.HEIGHTPROF = 0
807 982 self.counter_imagwr = 0
@@ -910,7 +1085,7 class ParametersPlot(Figure):
910 1085 # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
911 1086 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
912 1087 xlabel = ""
913 ylabel = "Range (Km)"
1088 ylabel = "Range (km)"
914 1089
915 1090 update_figfile = False
916 1091
@@ -954,9 +1129,28 class ParametersPlot(Figure):
954 1129
955 1130 self.setWinTitle(title)
956 1131
957 for i in range(self.nchan):
1132 # for i in range(self.nchan):
1133 # index = channelIndexList[i]
1134 # title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1135 # axes = self.axesList[i*self.plotFact]
1136 # z1 = z[i,:].reshape((1,-1))
1137 # axes.pcolorbuffer(x, y, z1,
1138 # xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
1139 # xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1140 # ticksize=9, cblabel='', cbsize="1%",colormap=colormap)
1141 #
1142 # if showSNR:
1143 # title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1144 # axes = self.axesList[i*self.plotFact + 1]
1145 # SNRdB1 = SNRdB[i,:].reshape((1,-1))
1146 # axes.pcolorbuffer(x, y, SNRdB1,
1147 # xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1148 # xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1149 # ticksize=9, cblabel='', cbsize="1%",colormap='jet')
1150
1151 i=0
958 1152 index = channelIndexList[i]
959 title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1153 title = "Factor de reflectividad Z [dBZ]"
960 1154 axes = self.axesList[i*self.plotFact]
961 1155 z1 = z[i,:].reshape((1,-1))
962 1156 axes.pcolorbuffer(x, y, z1,
@@ -973,6 +1167,44 class ParametersPlot(Figure):
973 1167 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
974 1168 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
975 1169
1170 i=1
1171 index = channelIndexList[i]
1172 title = "Velocidad vertical Doppler [m/s]"
1173 axes = self.axesList[i*self.plotFact]
1174 z1 = z[i,:].reshape((1,-1))
1175 axes.pcolorbuffer(x, y, z1,
1176 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=-10, zmax=10,
1177 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1178 ticksize=9, cblabel='', cbsize="1%",colormap='seismic_r')
1179
1180 if showSNR:
1181 title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1182 axes = self.axesList[i*self.plotFact + 1]
1183 SNRdB1 = SNRdB[i,:].reshape((1,-1))
1184 axes.pcolorbuffer(x, y, SNRdB1,
1185 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1186 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1187 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
1188
1189 i=2
1190 index = channelIndexList[i]
1191 title = "Intensidad de lluvia [mm/h]"
1192 axes = self.axesList[i*self.plotFact]
1193 z1 = z[i,:].reshape((1,-1))
1194 axes.pcolorbuffer(x, y, z1,
1195 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=0, zmax=40,
1196 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1197 ticksize=9, cblabel='', cbsize="1%",colormap='ocean_r')
1198
1199 if showSNR:
1200 title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1201 axes = self.axesList[i*self.plotFact + 1]
1202 SNRdB1 = SNRdB[i,:].reshape((1,-1))
1203 axes.pcolorbuffer(x, y, SNRdB1,
1204 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1205 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1206 ticksize=9, cblabel='', cbsize="1%",colormap='jet')
1207
976 1208
977 1209 self.draw()
978 1210
Show file before | Show file after | Hide comments
@@ -568,7 +568,7 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
568 568
569 569 if self.flagNoMoreFiles:
570 570 self.dataOut.flagNoData = True
571 print 'NoData se vuelve true'
571 #print 'NoData se vuelve true'
572 572 return 0
573 573
574 574 self.fp=self.path
@@ -578,7 +578,7 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
578 578 self.dataOut.data_cspc =self.data_cspc
579 579 self.dataOut.data_output=self.data_output
580 580
581 print 'self.dataOut.data_output', shape(self.dataOut.data_output)
581 #print 'self.dataOut.data_output', shape(self.dataOut.data_output)
582 582
583 583 #self.removeDC()
584 584 return self.dataOut.data_spc
@@ -596,7 +596,7 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
596 596 '''
597 597
598 598 #The address of the folder is generated the name of the .fdt file that will be read
599 print "File: ",self.fileSelector+1
599 #print "File: ",self.fileSelector+1
600 600
601 601 if self.fileSelector < len(self.filenameList):
602 602
@@ -608,7 +608,7 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
608 608
609 609 self.readBlock() #Block reading
610 610 else:
611 print 'readFile FlagNoData becomes true'
611 #print 'readFile FlagNoData becomes true'
612 612 self.flagNoMoreFiles=True
613 613 self.dataOut.flagNoData = True
614 614 return 0
@@ -634,8 +634,8 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
634 634
635 635 '''
636 636
637 if self.BlockCounter < self.nFDTdataRecors-2:
638 print self.nFDTdataRecors, 'CONDICION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
637 if self.BlockCounter < self.nFDTdataRecors-1:
638 #print self.nFDTdataRecors, 'CONDICION'
639 639 if self.ReadMode==1:
640 640 rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter+1)
641 641 elif self.ReadMode==0:
@@ -675,7 +675,7 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
675 675 self.dataOut.outputInterval= self.dataOut.ippSeconds * self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles
676 676
677 677 self.data_output=numpy.ones([3,rheader.nHeights])*numpy.NaN
678 print 'self.data_output', shape(self.data_output)
678 #print 'self.data_output', shape(self.data_output)
679 679 self.dataOut.velocityX=[]
680 680 self.dataOut.velocityY=[]
681 681 self.dataOut.velocityV=[]
@@ -711,9 +711,17 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
711 711 if self.DualModeIndex==self.ReadMode:
712 712
713 713 self.data_fft = numpy.fromfile( startDATA, [('complex','<c8')],self.nProfiles*self.nChannels*self.nHeights )
714 #
715 # if len(self.data_fft) is not 101376:
716 #
717 # self.data_fft = numpy.empty(101376)
714 718
715 719 self.data_fft=self.data_fft.astype(numpy.dtype('complex'))
716 720
721
722
723
724
717 725 self.data_block=numpy.reshape(self.data_fft,(self.nHeights, self.nChannels, self.nProfiles ))
718 726
719 727 self.data_block = numpy.transpose(self.data_block, (1,2,0))
@@ -729,11 +737,7 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
729 737 z = self.data_spc.copy()#/factor
730 738 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
731 739 #zdB = 10*numpy.log10(z)
732 print ' '
733 print 'Z: '
734 print shape(z)
735 print ' '
736 print ' '
740
737 741
738 742 self.dataOut.data_spc=self.data_spc
739 743
@@ -789,367 +793,3 class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa
789 793 self.dataOut.ChanDist = self.ChanDist
790 794
791 795
792 # for Height in range(self.nHeights):
793 #
794 # for i in range(self.nRdPairs):
795 #
796 # '''****** Line of Data SPC ******'''
797 # zline=z[i,:,Height]
798 #
799 # '''****** DC is removed ******'''
800 # DC=Find(zline,numpy.amax(zline))
801 # zline[DC]=(zline[DC-1]+zline[DC+1])/2
802 #
803 #
804 # '''****** SPC is normalized ******'''
805 # FactNorm= zline.copy() / numpy.sum(zline.copy())
806 # FactNorm= FactNorm/numpy.sum(FactNorm)
807 #
808 # SmoothSPC=moving_average(FactNorm,N=3)
809 #
810 # xSamples = ar(range(len(SmoothSPC)))
811 # ySamples[i] = SmoothSPC-self.noise[i]
812 #
813 # for i in range(self.nRdPairs):
814 #
815 # '''****** Line of Data CSPC ******'''
816 # cspcLine=self.data_cspc[i,:,Height].copy()
817 #
818 #
819 #
820 # '''****** CSPC is normalized ******'''
821 # chan_index0 = self.dataOut.pairsList[i][0]
822 # chan_index1 = self.dataOut.pairsList[i][1]
823 # CSPCFactor= numpy.sum(ySamples[chan_index0]) * numpy.sum(ySamples[chan_index1])
824 #
825 #
826 # CSPCNorm= cspcLine.copy() / numpy.sqrt(CSPCFactor)
827 #
828 #
829 # CSPCSamples[i] = CSPCNorm-self.noise[i]
830 # coherence[i] = numpy.abs(CSPCSamples[i]) / numpy.sqrt(CSPCFactor)
831 #
832 # '''****** DC is removed ******'''
833 # DC=Find(coherence[i],numpy.amax(coherence[i]))
834 # coherence[i][DC]=(coherence[i][DC-1]+coherence[i][DC+1])/2
835 # coherence[i]= moving_average(coherence[i],N=2)
836 #
837 # phase[i] = moving_average( numpy.arctan2(CSPCSamples[i].imag, CSPCSamples[i].real),N=1)#*180/numpy.pi
838 #
839 #
840 # '''****** Getting fij width ******'''
841 #
842 # yMean=[]
843 # yMean2=[]
844 #
845 # for j in range(len(ySamples[1])):
846 # yMean=numpy.append(yMean,numpy.average([ySamples[0,j],ySamples[1,j],ySamples[2,j]]))
847 #
848 # '''******* Getting fitting Gaussian ******'''
849 # meanGauss=sum(xSamples*yMean) / len(xSamples)
850 # sigma=sum(yMean*(xSamples-meanGauss)**2) / len(xSamples)
851 # #print 'Height',Height,'SNR', meanGauss/sigma**2
852 #
853 # if (abs(meanGauss/sigma**2) > 0.0001) :
854 #
855 # try:
856 # popt,pcov = curve_fit(gaus,xSamples,yMean,p0=[1,meanGauss,sigma])
857 #
858 # if numpy.amax(popt)>numpy.amax(yMean)*0.3:
859 # FitGauss=gaus(xSamples,*popt)
860 #
861 # else:
862 # FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
863 # print 'Verificador: Dentro', Height
864 # except RuntimeError:
865 #
866 # try:
867 # for j in range(len(ySamples[1])):
868 # yMean2=numpy.append(yMean2,numpy.average([ySamples[1,j],ySamples[2,j]]))
869 # popt,pcov = curve_fit(gaus,xSamples,yMean2,p0=[1,meanGauss,sigma])
870 # FitGauss=gaus(xSamples,*popt)
871 # print 'Verificador: Exepcion1', Height
872 # except RuntimeError:
873 #
874 # try:
875 # popt,pcov = curve_fit(gaus,xSamples,ySamples[1],p0=[1,meanGauss,sigma])
876 # FitGauss=gaus(xSamples,*popt)
877 # print 'Verificador: Exepcion2', Height
878 # except RuntimeError:
879 # FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
880 # print 'Verificador: Exepcion3', Height
881 # else:
882 # FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean)
883 # #print 'Verificador: Fuera', Height
884 #
885 #
886 #
887 # Maximun=numpy.amax(yMean)
888 # eMinus1=Maximun*numpy.exp(-1)
889 #
890 # HWpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1)))
891 # HalfWidth= xFrec[HWpos]
892 # GCpos=Find(FitGauss, numpy.amax(FitGauss))
893 # Vpos=Find(FactNorm, numpy.amax(FactNorm))
894 # #Vpos=numpy.sum(FactNorm)/len(FactNorm)
895 # #Vpos=Find(FactNorm, min(FactNorm, key=lambda value:abs(value- numpy.mean(FactNorm) )))
896 # #print 'GCpos',GCpos, numpy.amax(FitGauss), 'HWpos',HWpos
897 # '''****** Getting Fij ******'''
898 #
899 # GaussCenter=xFrec[GCpos]
900 # if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0):
901 # Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001
902 # else:
903 # Fij=abs(GaussCenter-HalfWidth)+0.0000001
904 #
905 # '''****** Getting Frecuency range of significant data ******'''
906 #
907 # Rangpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10)))
908 #
909 # if Rangpos<GCpos:
910 # Range=numpy.array([Rangpos,2*GCpos-Rangpos])
911 # else:
912 # Range=numpy.array([2*GCpos-Rangpos,Rangpos])
913 #
914 # FrecRange=xFrec[Range[0]:Range[1]]
915 #
916 # #print 'FrecRange', FrecRange
917 # '''****** Getting SCPC Slope ******'''
918 #
919 # for i in range(self.nRdPairs):
920 #
921 # if len(FrecRange)>5 and len(FrecRange)<self.nProfiles*0.5:
922 # PhaseRange=moving_average(phase[i,Range[0]:Range[1]],N=3)
923 #
924 # slope, intercept, r_value, p_value, std_err = stats.linregress(FrecRange,PhaseRange)
925 # PhaseSlope[i]=slope
926 # PhaseInter[i]=intercept
927 # else:
928 # PhaseSlope[i]=0
929 # PhaseInter[i]=0
930 #
931 # # plt.figure(i+15)
932 # # plt.title('FASE ( CH%s*CH%s )' %(self.dataOut.pairsList[i][0],self.dataOut.pairsList[i][1]))
933 # # plt.xlabel('Frecuencia (KHz)')
934 # # plt.ylabel('Magnitud')
935 # # #plt.subplot(311+i)
936 # # plt.plot(FrecRange,PhaseRange,'b')
937 # # plt.plot(FrecRange,FrecRange*PhaseSlope[i]+PhaseInter[i],'r')
938 #
939 # #plt.axis([-0.6, 0.2, -3.2, 3.2])
940 #
941 #
942 # '''Getting constant C'''
943 # cC=(Fij*numpy.pi)**2
944 #
945 # # '''Getting Eij and Nij'''
946 # # (AntennaX0,AntennaY0)=pol2cart(rheader.AntennaCoord0, rheader.AntennaAngl0*numpy.pi/180)
947 # # (AntennaX1,AntennaY1)=pol2cart(rheader.AntennaCoord1, rheader.AntennaAngl1*numpy.pi/180)
948 # # (AntennaX2,AntennaY2)=pol2cart(rheader.AntennaCoord2, rheader.AntennaAngl2*numpy.pi/180)
949 # #
950 # # E01=AntennaX0-AntennaX1
951 # # N01=AntennaY0-AntennaY1
952 # #
953 # # E02=AntennaX0-AntennaX2
954 # # N02=AntennaY0-AntennaY2
955 # #
956 # # E12=AntennaX1-AntennaX2
957 # # N12=AntennaY1-AntennaY2
958 #
959 # '''****** Getting constants F and G ******'''
960 # MijEijNij=numpy.array([[E02,N02], [E12,N12]])
961 # MijResult0=(-PhaseSlope[1]*cC) / (2*numpy.pi)
962 # MijResult1=(-PhaseSlope[2]*cC) / (2*numpy.pi)
963 # MijResults=numpy.array([MijResult0,MijResult1])
964 # (cF,cG) = numpy.linalg.solve(MijEijNij, MijResults)
965 #
966 # '''****** Getting constants A, B and H ******'''
967 # W01=numpy.amax(coherence[0])
968 # W02=numpy.amax(coherence[1])
969 # W12=numpy.amax(coherence[2])
970 #
971 # WijResult0=((cF*E01+cG*N01)**2)/cC - numpy.log(W01 / numpy.sqrt(numpy.pi/cC))
972 # WijResult1=((cF*E02+cG*N02)**2)/cC - numpy.log(W02 / numpy.sqrt(numpy.pi/cC))
973 # WijResult2=((cF*E12+cG*N12)**2)/cC - numpy.log(W12 / numpy.sqrt(numpy.pi/cC))
974 #
975 # WijResults=numpy.array([WijResult0, WijResult1, WijResult2])
976 #
977 # WijEijNij=numpy.array([ [E01**2, N01**2, 2*E01*N01] , [E02**2, N02**2, 2*E02*N02] , [E12**2, N12**2, 2*E12*N12] ])
978 # (cA,cB,cH) = numpy.linalg.solve(WijEijNij, WijResults)
979 #
980 # VxVy=numpy.array([[cA,cH],[cH,cB]])
981 #
982 # VxVyResults=numpy.array([-cF,-cG])
983 # (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
984 # Vzon = Vy
985 # Vmer = Vx
986 # Vmag=numpy.sqrt(Vzon**2+Vmer**2)
987 # Vang=numpy.arctan2(Vmer,Vzon)
988 #
989 # if abs(Vy)<100 and abs(Vy)> 0.:
990 # self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, Vzon) #Vmag
991 # #print 'Vmag',Vmag
992 # else:
993 # self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, NaN)
994 #
995 # if abs(Vx)<100 and abs(Vx) > 0.:
996 # self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, Vmer) #Vang
997 # #print 'Vang',Vang
998 # else:
999 # self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, NaN)
1000 #
1001 # if abs(GaussCenter)<2:
1002 # self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, xFrec[Vpos])
1003 #
1004 # else:
1005 # self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, NaN)
1006 #
1007 #
1008 # # print '********************************************'
1009 # # print 'HalfWidth ', HalfWidth
1010 # # print 'Maximun ', Maximun
1011 # # print 'eMinus1 ', eMinus1
1012 # # print 'Rangpos ', Rangpos
1013 # # print 'GaussCenter ',GaussCenter
1014 # # print 'E01 ',E01
1015 # # print 'N01 ',N01
1016 # # print 'E02 ',E02
1017 # # print 'N02 ',N02
1018 # # print 'E12 ',E12
1019 # # print 'N12 ',N12
1020 # #print 'self.dataOut.velocityX ', self.dataOut.velocityX
1021 # # print 'Fij ', Fij
1022 # # print 'cC ', cC
1023 # # print 'cF ', cF
1024 # # print 'cG ', cG
1025 # # print 'cA ', cA
1026 # # print 'cB ', cB
1027 # # print 'cH ', cH
1028 # # print 'Vx ', Vx
1029 # # print 'Vy ', Vy
1030 # # print 'Vmag ', Vmag
1031 # # print 'Vang ', Vang*180/numpy.pi
1032 # # print 'PhaseSlope ',PhaseSlope[0]
1033 # # print 'PhaseSlope ',PhaseSlope[1]
1034 # # print 'PhaseSlope ',PhaseSlope[2]
1035 # # print '********************************************'
1036 # #print 'data_output',shape(self.dataOut.velocityX), shape(self.dataOut.velocityY)
1037 #
1038 # #print 'self.dataOut.velocityX', len(self.dataOut.velocityX)
1039 # #print 'self.dataOut.velocityY', len(self.dataOut.velocityY)
1040 # #print 'self.dataOut.velocityV', self.dataOut.velocityV
1041 #
1042 # self.data_output[0]=numpy.array(self.dataOut.velocityX)
1043 # self.data_output[1]=numpy.array(self.dataOut.velocityY)
1044 # self.data_output[2]=numpy.array(self.dataOut.velocityV)
1045 #
1046 # prin= self.data_output[0][~numpy.isnan(self.data_output[0])]
1047 # print ' '
1048 # print 'VmagAverage',numpy.mean(prin)
1049 # print ' '
1050 # # plt.figure(5)
1051 # # plt.subplot(211)
1052 # # plt.plot(self.dataOut.velocityX,'yo:')
1053 # # plt.subplot(212)
1054 # # plt.plot(self.dataOut.velocityY,'yo:')
1055 #
1056 # # plt.figure(1)
1057 # # # plt.subplot(121)
1058 # # # plt.plot(xFrec,ySamples[0],'k',label='Ch0')
1059 # # # plt.plot(xFrec,ySamples[1],'g',label='Ch1')
1060 # # # plt.plot(xFrec,ySamples[2],'r',label='Ch2')
1061 # # # plt.plot(xFrec,FitGauss,'yo:',label='fit')
1062 # # # plt.legend()
1063 # # plt.title('DATOS A ALTURA DE 2850 METROS')
1064 # #
1065 # # plt.xlabel('Frecuencia (KHz)')
1066 # # plt.ylabel('Magnitud')
1067 # # # plt.subplot(122)
1068 # # # plt.title('Fit for Time Constant')
1069 # # #plt.plot(xFrec,zline)
1070 # # #plt.plot(xFrec,SmoothSPC,'g')
1071 # # plt.plot(xFrec,FactNorm)
1072 # # plt.axis([-4, 4, 0, 0.15])
1073 # # # plt.xlabel('SelfSpectra KHz')
1074 # #
1075 # # plt.figure(10)
1076 # # # plt.subplot(121)
1077 # # plt.plot(xFrec,ySamples[0],'b',label='Ch0')
1078 # # plt.plot(xFrec,ySamples[1],'y',label='Ch1')
1079 # # plt.plot(xFrec,ySamples[2],'r',label='Ch2')
1080 # # # plt.plot(xFrec,FitGauss,'yo:',label='fit')
1081 # # plt.legend()
1082 # # plt.title('SELFSPECTRA EN CANALES')
1083 # #
1084 # # plt.xlabel('Frecuencia (KHz)')
1085 # # plt.ylabel('Magnitud')
1086 # # # plt.subplot(122)
1087 # # # plt.title('Fit for Time Constant')
1088 # # #plt.plot(xFrec,zline)
1089 # # #plt.plot(xFrec,SmoothSPC,'g')
1090 # # # plt.plot(xFrec,FactNorm)
1091 # # # plt.axis([-4, 4, 0, 0.15])
1092 # # # plt.xlabel('SelfSpectra KHz')
1093 # #
1094 # # plt.figure(9)
1095 # #
1096 # #
1097 # # plt.title('DATOS SUAVIZADOS')
1098 # # plt.xlabel('Frecuencia (KHz)')
1099 # # plt.ylabel('Magnitud')
1100 # # plt.plot(xFrec,SmoothSPC,'g')
1101 # #
1102 # # #plt.plot(xFrec,FactNorm)
1103 # # plt.axis([-4, 4, 0, 0.15])
1104 # # # plt.xlabel('SelfSpectra KHz')
1105 # # #
1106 # # plt.figure(2)
1107 # # # #plt.subplot(121)
1108 # # plt.plot(xFrec,yMean,'r',label='Mean SelfSpectra')
1109 # # plt.plot(xFrec,FitGauss,'yo:',label='Ajuste Gaussiano')
1110 # # # plt.plot(xFrec[Rangpos],FitGauss[Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.1)))],'bo')
1111 # # # #plt.plot(xFrec,phase)
1112 # # # plt.xlabel('Suavizado, promediado KHz')
1113 # # plt.title('SELFSPECTRA PROMEDIADO')
1114 # # # #plt.subplot(122)
1115 # # # #plt.plot(xSamples,zline)
1116 # # plt.xlabel('Frecuencia (KHz)')
1117 # # plt.ylabel('Magnitud')
1118 # # plt.legend()
1119 # # #
1120 # # # plt.figure(3)
1121 # # # plt.subplot(311)
1122 # # # #plt.plot(xFrec,phase[0])
1123 # # # plt.plot(xFrec,phase[0],'g')
1124 # # # plt.subplot(312)
1125 # # # plt.plot(xFrec,phase[1],'g')
1126 # # # plt.subplot(313)
1127 # # # plt.plot(xFrec,phase[2],'g')
1128 # # # #plt.plot(xFrec,phase[2])
1129 # # #
1130 # # # plt.figure(4)
1131 # # #
1132 # # # plt.plot(xSamples,coherence[0],'b')
1133 # # # plt.plot(xSamples,coherence[1],'r')
1134 # # # plt.plot(xSamples,coherence[2],'g')
1135 # # plt.show()
1136 # # #
1137 # # # plt.clf()
1138 # # # plt.cla()
1139 # # # plt.close()
1140 #
1141 # print ' '
1142
1143
1144
1145 self.BlockCounter+=2
1146
1147 else:
1148 self.fileSelector+=1
1149 self.BlockCounter=0
1150 print "Next File"
1151
1152
1153
1154
1155
Show file before | Show file after | Hide comments
@@ -645,7 +645,10 class ParamWriter(Operation):
645 645 dsDict['variable'] = self.dataList[i]
646 646 #--------------------- Conditionals ------------------------
647 647 #There is no data
648
649
648 650 if dataAux is None:
651
649 652 return 0
650 653
651 654 #Not array, just a number
@@ -1078,6 +1081,7 class ParamWriter(Operation):
1078 1081 def run(self, dataOut, **kwargs):
1079 1082
1080 1083 if not(self.isConfig):
1084
1081 1085 flagdata = self.setup(dataOut, **kwargs)
1082 1086
1083 1087 if not(flagdata):
Show file before | Show file after | Hide comments
@@ -182,7 +182,7 class ParametersProc(ProcessingUnit):
182 182 def target(tups):
183 183
184 184 obj, args = tups
185 #print 'TARGETTT', obj, args
185
186 186 return obj.FitGau(args)
187 187
188 188
@@ -238,12 +238,6 class SpectralFilters(Operation):
238 238
239 239 Delta = self.Num_Bin/2 - Breaker1R[0]
240 240
241 #Breaker1W=VelRange[numpy.abs(VelRange-(-LimitW)).argmin()]
242 #Breaker1W=numpy.where(VelRange == Breaker1W)
243
244 #Breaker2W=VelRange[numpy.abs(VelRange-(LimitW)).argmin()]
245 #Breaker2W=numpy.where(VelRange == Breaker2W)
246
247 241
248 242 '''Reacomodando SPCrange'''
249 243
@@ -279,16 +273,6 class SpectralFilters(Operation):
279 273 SPCroll[i]=SPCroll[i]-dataOut.noise[i]
280 274 SPCroll[ numpy.where( SPCroll<0 ) ] = 1e-20
281 275
282 #self.spc[i, 0:int(Breaker1W[0]) ,:] = dataOut.noise[i]
283 #self.spc[i, int(Breaker2W[0]):self.Num_Bin ,:] = dataOut.noise[i]
284
285 #self.cspc[i, 0:int(Breaker1W[0]) ,:] = dataOut.noise[i]
286 #self.cspc[i, int(Breaker2W[0]):self.Num_Bin ,:] = dataOut.noise[i]
287
288
289
290
291
292 276 SPC_ch1 = SPCroll
293 277
294 278 SPC_ch2 = SPCcut
@@ -296,9 +280,6 class SpectralFilters(Operation):
296 280 SPCparam = (SPC_ch1, SPC_ch2, self.spc)
297 281 dataOut.SPCparam = numpy.asarray(SPCparam)
298 282
299 #dataOut.data_pre= (self.spc , self.cspc)
300
301 #dataOut.data_preParam = (self.spc , self.cspc)
302 283
303 284 dataOut.spcparam_range=numpy.zeros([self.Num_Chn,self.Num_Bin+1])
304 285
@@ -307,8 +288,6 class SpectralFilters(Operation):
307 288 dataOut.spcparam_range[0]=FrecRange
308 289
309 290
310
311
312 291 class GaussianFit(Operation):
313 292
314 293 '''
@@ -338,21 +317,7 class GaussianFit(Operation):
338 317 self.spc = dataOut.data_pre[0].copy()
339 318
340 319
341 print 'SelfSpectra Shape', numpy.asarray(self.spc).shape
342
343
344 #plt.figure(50)
345 #plt.subplot(121)
346 #plt.plot(self.spc,'k',label='spc(66)')
347 #plt.plot(xFrec,ySamples[1],'g',label='Ch1')
348 #plt.plot(xFrec,ySamples[2],'r',label='Ch2')
349 #plt.plot(xFrec,FitGauss,'yo:',label='fit')
350 #plt.legend()
351 #plt.title('DATOS A ALTURA DE 7500 METROS')
352 #plt.show()
353
354 320 self.Num_Hei = self.spc.shape[2]
355 #self.Num_Bin = len(self.spc)
356 321 self.Num_Bin = self.spc.shape[1]
357 322 self.Num_Chn = self.spc.shape[0]
358 323 Vrange = dataOut.abscissaList
@@ -376,16 +341,9 class GaussianFit(Operation):
376 341 gauSPC = pool.map(target, attrs)
377 342 dataOut.SPCparam = numpy.asarray(SPCparam)
378 343
379
380
381 print '========================================================'
382 print 'total_time: ', time.time()-start_time
383
384 344 # re-normalizing spc and noise
385 345 # This part differs from gg1
386 346
387
388
389 347 ''' Parameters:
390 348 1. Amplitude
391 349 2. Shift
@@ -398,9 +356,6 class GaussianFit(Operation):
398 356 def FitGau(self, X):
399 357
400 358 Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X
401 #print 'VARSSSS', ch, pnoise, noise, num_intg
402
403 #print 'HEIGHTS', self.Num_Hei
404 359
405 360 SPCparam = []
406 361 SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei])
@@ -411,10 +366,6 class GaussianFit(Operation):
411 366
412 367
413 368 for ht in range(self.Num_Hei):
414 #print (numpy.asarray(self.spc).shape)
415
416 #print 'TTTTT', ch , ht
417 #print self.spc.shape
418 369
419 370
420 371 spc = numpy.asarray(self.spc)[ch,:,ht]
@@ -436,16 +387,13 class GaussianFit(Operation):
436 387 noisebl=wnoise*0.9;
437 388 noisebh=wnoise*1.1
438 389 spc=spc-wnoise
439 # print 'wnoise', noise_[0], spc_norm_max, wnoise
390
440 391 minx=numpy.argmin(spc)
441 392 #spcs=spc.copy()
442 393 spcs=numpy.roll(spc,-minx)
443 394 cum=numpy.cumsum(spcs)
444 395 tot_noise=wnoise * self.Num_Bin #64;
445 #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise
446 #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? '''
447 #snr=tot_signal/tot_noise
448 #snr=cum[-1]/tot_noise
396
449 397 snr = sum(spcs)/tot_noise
450 398 snrdB=10.*numpy.log10(snr)
451 399
@@ -455,8 +403,6 class GaussianFit(Operation):
455 403 SPC_ch1[:,ht] = 0#numpy.NaN
456 404 SPCparam = (SPC_ch1,SPC_ch2)
457 405 continue
458 #print 'snr',snrdB #, sum(spcs) , tot_noise
459
460 406
461 407
462 408 #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4:
@@ -602,57 +548,10 class GaussianFit(Operation):
602 548 shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN]
603 549
604 550
605 # if choice==0: # pick the single gaussian fit
606 # Amplitude0=lsq1[0][2]
607 # shift0=lsq1[0][0]
608 # width0=lsq1[0][1]
609 # p0=lsq1[0][3]
610 # Amplitude1=0.
611 # shift1=0.
612 # width1=0.
613 # p1=0.
614 # noise=lsq1[0][4]
615 # elif choice==1: # take the first one of the 2 gaussians fitted
616 # Amplitude0 = lsq2[0][2]
617 # shift0 = lsq2[0][0]
618 # width0 = lsq2[0][1]
619 # p0 = lsq2[0][3]
620 # Amplitude1 = lsq2[0][6] # This is 0 in gg1
621 # shift1 = lsq2[0][4] # This is 0 in gg1
622 # width1 = lsq2[0][5] # This is 0 in gg1
623 # p1 = lsq2[0][7] # This is 0 in gg1
624 # noise = lsq2[0][8]
625 # else: # the second one
626 # Amplitude0 = lsq2[0][6]
627 # shift0 = lsq2[0][4]
628 # width0 = lsq2[0][5]
629 # p0 = lsq2[0][7]
630 # Amplitude1 = lsq2[0][2] # This is 0 in gg1
631 # shift1 = lsq2[0][0] # This is 0 in gg1
632 # width1 = lsq2[0][1] # This is 0 in gg1
633 # p1 = lsq2[0][3] # This is 0 in gg1
634 # noise = lsq2[0][8]
635
636 #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0)
637 551 SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0
638 552 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1
639 #print 'SPC_ch1.shape',SPC_ch1.shape
640 #print 'SPC_ch2.shape',SPC_ch2.shape
641 #dataOut.data_param = SPC_ch1
642 553 SPCparam = (SPC_ch1,SPC_ch2)
643 #GauSPC[1] = SPC_ch2
644
645 # print 'shift0', shift0
646 # print 'Amplitude0', Amplitude0
647 # print 'width0', width0
648 # print 'p0', p0
649 # print '========================'
650 # print 'shift1', shift1
651 # print 'Amplitude1', Amplitude1
652 # print 'width1', width1
653 # print 'p1', p1
654 # print 'noise', noise
655 # print 's_noise', wnoise
554
656 555
657 556 return GauSPC
658 557
@@ -703,6 +602,12 class PrecipitationProc(Operation):
703 602
704 603 Parameters affected:
705 604 '''
605
606 def __init__(self, **kwargs):
607 Operation.__init__(self, **kwargs)
608 self.i=0
609
610
706 611 def gaus(self,xSamples,Amp,Mu,Sigma):
707 612 return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) ))
708 613
@@ -739,11 +644,16 class PrecipitationProc(Operation):
739 644 else:
740 645
741 646 self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() #
647
648 """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
649
650 self.spc[:,:,0:7]= numpy.NaN
651
652 """##########################################"""
653
742 654 self.Num_Hei = self.spc.shape[2]
743 655 self.Num_Bin = self.spc.shape[1]
744 656 self.Num_Chn = self.spc.shape[0]
745 print '==================== SPC SHAPE',numpy.shape(self.spc)
746
747 657
748 658 ''' Se obtiene la constante del RADAR '''
749 659
@@ -758,10 +668,8 class PrecipitationProc(Operation):
758 668
759 669 Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) )
760 670 Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR)
761 RadarConstant = 5e-27 * Numerator / Denominator #
762 print '***'
763 print '*** RadarConstant' , RadarConstant, '****'
764 print '***'
671 RadarConstant = 5e-26 * Numerator / Denominator #
672
765 673 ''' ============================= '''
766 674
767 675 self.spc[0] = (self.spc[0]-dataOut.noise[0])
@@ -781,8 +689,6 class PrecipitationProc(Operation):
781 689
782 690 VelMeteoro = numpy.mean(SPCmean,axis=0)
783 691
784 #print '==================== Vel SHAPE',VelMeteoro
785
786 692 D_range = numpy.zeros([self.Num_Bin,self.Num_Hei])
787 693 SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei])
788 694 N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei])
@@ -847,7 +753,7 class PrecipitationProc(Operation):
847 753 dBZe = 10*numpy.log10(Ze)
848 754 dBZ = 10*numpy.log10(Z)
849 755
850 dataOut.data_output = Z
756 dataOut.data_output = RR[8]
851 757 dataOut.data_param = numpy.ones([3,self.Num_Hei])
852 758 dataOut.channelList = [0,1,2]
853 759
@@ -855,31 +761,6 class PrecipitationProc(Operation):
855 761 dataOut.data_param[1]=V_mean
856 762 dataOut.data_param[2]=RR
857 763
858 #print 'VELRANGE', Velrange
859 #print 'Range', len(Range)
860 #print 'delv',del_V
861 # print 'DRANGE', D_range[:,56]
862 # #print 'NOISE', dataOut.noise[0], 10*numpy.log10(dataOut.noise[0])
863 # print 'radarconstant', RadarConstant
864 # #print 'ETAn SHAPE', ETAn.shape
865 # # print 'ETAn ', numpy.nansum(ETAn, axis=0)
866 # # print 'ETAd ', numpy.nansum(ETAd, axis=0)
867 # print 'Pr ', numpy.nansum(Pr, axis=0)
868 # print 'dataOut.SPCparam[1]', numpy.nansum(dataOut.SPCparam[1][0], axis=0)
869 # print 'Ze ', dBZe
870 print 'Z ', dBZ
871 # print 'Ndist',N_dist[:,56]
872 print 'RR2 ', RR2
873 print 'RR ', RR
874 print 'Vr', V_mean
875 #print 'RR2 ', dBRR2
876 #print 'D_mean', D_mean
877 #print 'del_V', del_V
878 #print 'D_range',D_range.shape, D_range[:,30]
879 #print 'Velrange', Velrange
880 #print 'numpy.nansum( N_dist[:,R]', numpy.nansum( N_dist, axis=0)
881 #print 'dataOut.data_param SHAPE', dataOut.data_param.shape
882
883 764
884 765 def dBZeMODE2(self, dataOut): # Processing for MIRA35C
885 766
@@ -894,7 +775,7 class PrecipitationProc(Operation):
894 775 data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN
895 776
896 777 ETA = numpy.sum(SNR,1)
897 print 'ETA' , ETA
778
898 779 ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN)
899 780
900 781 Ze = numpy.ones([self.Num_Chn, self.Num_Hei] )
@@ -956,6 +837,14 class FullSpectralAnalysis(Operation):
956 837 spc = dataOut.data_pre[0].copy()
957 838 cspc = dataOut.data_pre[1]
958 839
840 """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX"""
841
842 SNRspc = spc.copy()
843 SNRspc[:,:,0:7]= numpy.NaN
844
845 """##########################################"""
846
847
959 848 nChannel = spc.shape[0]
960 849 nProfiles = spc.shape[1]
961 850 nHeights = spc.shape[2]
@@ -979,16 +868,11 class FullSpectralAnalysis(Operation):
979 868 data = dataOut.data_pre
980 869 noise = dataOut.noise
981 870
982 dataOut.data_SNR = (numpy.mean(spc,axis=1)- noise[0]) / noise[0]
871 dataOut.data_SNR = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0]
983 872
984 873 dataOut.data_SNR[numpy.where( dataOut.data_SNR <0 )] = 1e-20
985 874
986 875
987 #FirstMoment = dataOut.moments[0,1,:]#numpy.average(dataOut.data_param[:,1,:],0)
988 #SecondMoment = numpy.average(dataOut.moments[:,2,:],0)
989
990 #SNRdBMean = []
991
992 876 data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN
993 877
994 878 velocityX=[]
@@ -1008,25 +892,19 class FullSpectralAnalysis(Operation):
1008 892 velocityX=numpy.append(velocityX, Vzon)#Vmag
1009 893
1010 894 else:
1011 #print 'Vzon',Vzon
1012 895 velocityX=numpy.append(velocityX, numpy.NaN)
1013 896
1014 897 if abs(Vmer)<100. and abs(Vmer) > 0.:
1015 898 velocityY=numpy.append(velocityY, -Vmer)#Vang
1016 899
1017 900 else:
1018 #print 'Vmer',Vmer
1019 901 velocityY=numpy.append(velocityY, numpy.NaN)
1020 902
1021 903 if dbSNR[Height] > SNRlimit:
1022 904 velocityV=numpy.append(velocityV, -Vver)#FirstMoment[Height])
1023 905 else:
1024 906 velocityV=numpy.append(velocityV, numpy.NaN)
1025 #FirstMoment[Height]= numpy.NaN
1026 # if SNRdBMean[Height] <12:
1027 # FirstMoment[Height] = numpy.NaN
1028 # velocityX[Height] = numpy.NaN
1029 # velocityY[Height] = numpy.NaN
907
1030 908
1031 909
1032 910
@@ -1034,21 +912,6 class FullSpectralAnalysis(Operation):
1034 912 data_output[1] = numpy.array(velocityY) #self.moving_average(numpy.array(velocityY) , N=1)
1035 913 data_output[2] = velocityV#FirstMoment
1036 914
1037 print 'data_output', data_output.shape
1038 #print FirstMoment
1039 # print 'velocityX',numpy.shape(data_output[0])
1040 # print 'velocityX',data_output[0]
1041 # print ' '
1042 # print 'velocityY',numpy.shape(data_output[1])
1043 # print 'velocityY',data_output[1]
1044 # print 'velocityV',data_output[2]
1045 # print 'PhaseLine',PhaseLine
1046 #print numpy.array(velocityY)
1047 #print 'SNR'
1048 #print 10*numpy.log10(dataOut.data_SNR)
1049 #print numpy.shape(10*numpy.log10(dataOut.data_SNR))
1050 print ' '
1051
1052 915 xFrec=FrecRange[0:spc.shape[1]]
1053 916
1054 917 dataOut.data_output=data_output
@@ -1144,15 +1007,6 class FullSpectralAnalysis(Operation):
1144 1007 self.Moments(numpy.abs(CSPCSamples[1]), xSamples),
1145 1008 self.Moments(numpy.abs(CSPCSamples[2]), xSamples)])
1146 1009
1147 #print '##### SUMA de SPC #####', len(ySamples)
1148 #print numpy.sum(ySamples[0])
1149 #print '##### SUMA de CSPC #####', len(coherence)
1150 #print numpy.sum(numpy.abs(CSPCNorm))
1151 #print numpy.sum(coherence[0])
1152 # print 'len',len(xSamples)
1153 # print 'CSPCmoments', numpy.shape(CSPCmoments)
1154 # print CSPCmoments
1155 # print '#######################'
1156 1010
1157 1011 popt=[1e-10,0,1e-10]
1158 1012 popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10]
@@ -1235,24 +1089,6 class FullSpectralAnalysis(Operation):
1235 1089 else:
1236 1090 Fij = xSamples[PointGauCenter] - xSamples[PointFij]
1237 1091
1238 # print 'CSPCopt'
1239 # print CSPCopt
1240 # print 'popt'
1241 # print popt
1242 # print '#######################################'
1243 #print 'dataOut.data_param', numpy.shape(data_param)
1244 #print 'dataOut.data_param0', data_param[0,0,Height]
1245 #print 'dataOut.data_param1', data_param[0,1,Height]
1246 #print 'dataOut.data_param2', data_param[0,2,Height]
1247
1248
1249 # print 'yMoments', yMoments
1250 # print 'Moments', SPCmoments
1251 # print 'Fij2 Moment', Fij
1252 # #print 'Fij', Fij, 'popt[2]/2',popt[2]/2
1253 # print 'Fijcspc',Fijcspc
1254 # print '#######################################'
1255
1256 1092
1257 1093 '''****** Taking frequency ranges from SPCs ******'''
1258 1094
@@ -1275,10 +1111,6 class FullSpectralAnalysis(Operation):
1275 1111 VelRange = xVel[ Range[0] : Range[1] ]
1276 1112
1277 1113
1278 #print 'RANGE: ', Range
1279 #print 'FrecRange', numpy.shape(FrecRange)#,FrecRange
1280 #print 'len: ', len(FrecRange)
1281
1282 1114 '''****** Getting SCPC Slope ******'''
1283 1115
1284 1116 for i in range(spc.shape[0]):
@@ -1286,13 +1118,6 class FullSpectralAnalysis(Operation):
1286 1118 if len(FrecRange)>5 and len(FrecRange)<spc.shape[1]*0.3:
1287 1119 PhaseRange=self.moving_average(phase[i,Range[0]:Range[1]],N=3)
1288 1120
1289 #print 'Ancho espectral Frecuencias', FrecRange[-1]-FrecRange[0], 'Hz'
1290 #print 'Ancho espectral Velocidades', VelRange[-1]-VelRange[0], 'm/s'
1291 #print 'FrecRange', len(FrecRange) , FrecRange
1292 #print 'VelRange', len(VelRange) , VelRange
1293 #print 'PhaseRange', numpy.shape(PhaseRange), PhaseRange
1294 #print ' '
1295
1296 1121 '''***********************VelRange******************'''
1297 1122
1298 1123 mask = ~numpy.isnan(FrecRange) & ~numpy.isnan(PhaseRange)
@@ -1341,12 +1166,6 class FullSpectralAnalysis(Operation):
1341 1166 VxVyResults=numpy.array([-cF,-cG])
1342 1167 (Vx,Vy) = numpy.linalg.solve(VxVy, VxVyResults)
1343 1168
1344 #print 'MijResults, cC, PhaseSlope', MijResults, cC, PhaseSlope
1345 #print 'W01,02,12', W01, W02, W12
1346 #print 'WijResult0,1,2',WijResult0, WijResult1, WijResult2, 'Results', WijResults
1347 #print 'cA,cB,cH, cF, cG', cA, cB, cH, cF, cG
1348 #print 'VxVy', VxVyResults
1349 #print '###########################****************************************'
1350 1169 Vzon = Vy
1351 1170 Vmer = Vx
1352 1171 Vmag=numpy.sqrt(Vzon**2+Vmer**2)
@@ -1358,42 +1177,6 class FullSpectralAnalysis(Operation):
1358 1177 FitGaussCSPC = numpy.array([FitGauss01,FitGauss02,FitGauss12])
1359 1178
1360 1179
1361 # ''' Ploteo por altura '''
1362 # if Height == 28:
1363 # for i in range(3):
1364 # #print 'FASE', numpy.shape(phase), y[25]
1365 # #print numpy.shape(coherence)
1366 # fig = plt.figure(10+self.indice)
1367 # #plt.plot( x[0:256],coherence[:,25] )
1368 # #cohAv = numpy.average(coherence[i],1)
1369 # Pendiente = FrecRange * PhaseSlope[i]
1370 # plt.plot( FrecRange, Pendiente)
1371 # plt.plot( xFrec,phase[i])
1372 #
1373 # CSPCmean = numpy.mean(numpy.abs(CSPCSamples),0)
1374 # #plt.plot(xFrec, FitGauss01)
1375 # #plt.plot(xFrec, CSPCmean)
1376 # #plt.plot(xFrec, numpy.abs(CSPCSamples[0]))
1377 # #plt.plot(xFrec, FitGauss)
1378 # #plt.plot(xFrec, yMean)
1379 # #plt.plot(xFrec, numpy.abs(coherence[0]))
1380 #
1381 # #plt.axis([-12, 12, 15, 50])
1382 # #plt.title("%s" %( '%s %s, Channel %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S") , i)))
1383 # plt.ylabel('Desfase [rad]')
1384 # #plt.ylabel('CSPC normalizado')
1385 # plt.xlabel('Frec range [Hz]')
1386
1387 #fig.savefig('/home/erick/Documents/Pics/to{}.png'.format(self.indice))
1388
1389 # plt.show()
1390 # self.indice=self.indice+1
1391
1392
1393
1394
1395
1396 # print 'vzon y vmer', Vzon, Vmer
1397 1180 return Vzon, Vmer, Vver, GaussCenter, PhaseSlope, FitGaussCSPC
1398 1181
1399 1182 class SpectralMoments(Operation):
@@ -1711,7 +1494,6 class SpectralFitting(Operation):
1711 1494 dataCross = dataCross**2/K
1712 1495
1713 1496 for h in range(nHeights):
1714 # print self.dataOut.heightList[h]
1715 1497
1716 1498 #Input
1717 1499 d = data[:,h]
Show file before | Show file after | Hide comments
@@ -319,6 +319,17 class SpectraProc(ProcessingUnit):
319 319 return 1
320 320
321 321
322 def setH0(self, h0, deltaHeight = None):
323
324 if not deltaHeight:
325 deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
326
327 nHeights = self.dataOut.nHeights
328
329 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
330
331 self.dataOut.heightList = newHeiRange
332
322 333
323 334 def selectHeights(self, minHei, maxHei):
324 335 """
Show file before | Show file after | Hide comments
@@ -1,1 +1,1
1 <Project description="Segundo Test" id="191" name="test01"><ReadUnit datatype="SpectraReader" id="1911" inputId="0" name="SpectraReader"><Operation id="19111" name="run" priority="1" type="self"><Parameter format="str" id="191111" name="datatype" value="SpectraReader" /><Parameter format="str" id="191112" name="path" value="/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/2018" /><Parameter format="date" id="191113" name="startDate" value="2018/01/26" /><Parameter format="date" id="191114" name="endDate" value="2018/01/26" /><Parameter format="time" id="191115" name="startTime" value="17:45:00" /><Parameter format="time" id="191116" name="endTime" value="23:59:00" /><Parameter format="int" id="191118" name="online" value="0" /><Parameter format="int" id="191119" name="walk" value="1" /></Operation><Operation id="19112" name="printInfo" priority="2" type="self" /><Operation id="19113" name="printNumberOfBlock" priority="3" type="self" /></ReadUnit><ProcUnit datatype="Parameters" id="1913" inputId="1912" name="ParametersProc"><Operation id="19131" name="run" priority="1" type="self" /><Operation id="19132" name="SpectralFilters" priority="2" type="other"><Parameter format="float" id="191321" name="PositiveLimit" value="1.5" /><Parameter format="float" id="191322" name="NegativeLimit" value="12.5" /></Operation><Operation id="19133" name="PrecipitationProc" priority="3" type="other" /><Operation id="19134" name="ParametersPlot" priority="4" type="other"><Parameter format="int" id="191341" name="id" value="10" /><Parameter format="str" id="191342" name="wintitle" value="First_gg" /><Parameter format="str" id="191343" name="colormap" value="ocean_r" /><Parameter format="int" id="191344" name="zmin" value="00" /><Parameter format="int" id="191345" name="zmax" value="40" /><Parameter format="int" id="191346" name="ymin" value="0" /><Parameter format="int" id="191347" name="ymax" value="11" /><Parameter format="int" id="191348" name="xmin" value="17" /><Parameter format="int" id="191349" name="xmax" value="24" /><Parameter format="int" id="191350" name="save" value="1" /><Parameter format="str" id="191351" name="figpath" value="/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/2018" /></Operation><Operation id="19135" name="SpcParamPlot" priority="5" type="other"><Parameter format="int" id="191351" name="id" value="21" /><Parameter format="str" id="191352" name="wintitle" value="Primer eco removido" /><Parameter format="str" id="191353" name="xaxis" value="velocity" /><Parameter format="int" id="191354" name="showprofile" value="1" /><Parameter format="int" id="191355" name="zmin" value="10" /><Parameter format="int" id="191356" name="zmax" value="40" /><Parameter format="int" id="191357" name="ymin" value="0" /><Parameter format="int" id="191358" name="ymax" value="10" /><Parameter format="int" id="191359" name="Selector" value="1" /></Operation></ProcUnit><ProcUnit datatype="SpectraProc" id="1912" inputId="1911" name="SpectraProc"><Operation id="19121" name="run" priority="1" type="self" /><Operation id="19122" name="setRadarFrequency" priority="2" type="self"><Parameter format="float" id="191221" name="frequency" value="445.09e6" /></Operation></ProcUnit></Project> No newline at end of file
1 <Project description="Segundo Test" id="191" name="test01"><ReadUnit datatype="SpectraReader" id="1911" inputId="0" name="SpectraReader"><Operation id="19111" name="run" priority="1" type="self"><Parameter format="str" id="191111" name="datatype" value="SpectraReader" /><Parameter format="str" id="191112" name="path" value="/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra" /><Parameter format="date" id="191113" name="startDate" value="2018/02/23" /><Parameter format="date" id="191114" name="endDate" value="2018/02/23" /><Parameter format="time" id="191115" name="startTime" value="00:00:00" /><Parameter format="time" id="191116" name="endTime" value="03:59:50" /><Parameter format="int" id="191118" name="online" value="0" /><Parameter format="int" id="191119" name="walk" value="1" /></Operation><Operation id="19112" name="printInfo" priority="2" type="self" /><Operation id="19113" name="printNumberOfBlock" priority="3" type="self" /></ReadUnit><ProcUnit datatype="Parameters" id="1913" inputId="1912" name="ParametersProc"><Operation id="19131" name="run" priority="1" type="self" /><Operation id="19132" name="SpectralFilters" priority="2" type="other"><Parameter format="float" id="191321" name="PositiveLimit" value="1.5" /><Parameter format="float" id="191322" name="NegativeLimit" value="12.5" /></Operation><Operation id="19133" name="FullSpectralAnalysis" priority="3" type="other"><Parameter format="float" id="191331" name="SNRlimit" value="-16" /><Parameter format="float" id="191332" name="Xi01" value="1.500" /><Parameter format="float" id="191333" name="Xi02" value="1.500" /><Parameter format="float" id="191334" name="Xi12" value="0" /><Parameter format="float" id="191335" name="Eta01" value="0.875" /><Parameter format="float" id="191336" name="Eta02" value="-0.875" /><Parameter format="float" id="191337" name="Eta12" value="-1.750" /></Operation><Operation id="19134" name="WindProfilerPlot" priority="4" type="other"><Parameter format="int" id="191341" name="id" value="4" /><Parameter format="str" id="191342" name="wintitle" value="Wind Profiler" /><Parameter format="float" id="191343" name="xmin" value="0" /><Parameter format="float" id="191344" name="xmax" value="4" /><Parameter format="float" id="191345" name="ymin" value="0" /><Parameter format="int" id="191346" name="ymax" value="4" /><Parameter format="float" id="191347" name="zmin" value="-20" /><Parameter format="float" id="191348" name="zmax" value="20" /><Parameter format="float" id="191349" name="SNRmin" value="-20" /><Parameter format="float" id="191350" name="SNRmax" value="20" /><Parameter format="float" id="191351" name="zmin_ver" value="-200" /><Parameter format="float" id="191352" name="zmax_ver" value="200" /><Parameter format="float" id="191353" name="SNRthresh" value="-20" /><Parameter format="int" id="191354" name="save" value="1" /><Parameter format="str" id="191355" name="figpath" value="/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/ImagenesTesis" /></Operation><Operation id="19135" name="ParamWriter" priority="5" type="other"><Parameter format="str" id="191351" name="path" value="/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/ImagenesTesis" /><Parameter format="int" id="191352" name="blocksPerFile" value="500" /><Parameter format="list" id="191353" name="metadataList" value="heightList,timeZone,paramInterval" /><Parameter format="list" id="191354" name="dataList" value="data_output,utctime,utctimeInit" /></Operation></ProcUnit><ProcUnit datatype="SpectraProc" id="1912" inputId="1911" name="SpectraProc"><Operation id="19121" name="run" priority="1" type="self" /><Operation id="19122" name="setRadarFrequency" priority="2" type="self"><Parameter format="float" id="191221" name="frequency" value="445.09e6" /></Operation><Operation id="19123" name="setH0" priority="3" type="self"><Parameter format="float" id="191231" name="h0" value="0.500" /></Operation></ProcUnit></Project> No newline at end of file
Show file before | Show file after | Hide comments
@@ -13,7 +13,7 xmax = '24'
13 13
14 14 desc = "ProcBLTR Test"
15 15 filename = "ProcBLTR.xml"
16 figpath = '/media/erick/6F60F7113095A154/BLTR'
16 figpath = '/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/2018'
17 17
18 18 controllerObj = Project()
19 19
@@ -22,10 +22,11 controllerObj.setup(id='191', name='test01', description=desc)
22 22
23 23 readUnitConfObj = controllerObj.addReadUnit(datatype='BLTRSpectraReader',
24 24 #path='/media/erick/6F60F7113095A154/BLTR/',
25 #path='/data/BLTR',
25 26 path='/home/erick/Documents/Data/BLTR_Data/fdt/',
26 endDate='2017/10/19',
27 startTime='13:00:00',
28 startDate='2016/11/8',
27 endDate='2016/11/01',
28 startTime='0:00:00',
29 startDate='2016/11/01',
29 30 endTime='23:59:59',
30 31
31 32
@@ -82,12 +83,12 opObj10 = procUnitConfObj1.addOperation(name='removeDC')
82 83 # opObj11.addParameter(name='showprofile', value='1', format='int')
83 84 # opObj11.addParameter(name='timerange', value=str(2*60*60), format='int')
84 85
85 opObj11 = procUnitConfObj1.addOperation(name='CrossSpectraPlot', optype='other')
86 procUnitConfObj1.addParameter(name='pairsList', value='(0,1),(0,2),(1,2)', format='pairsList')
87 opObj11.addParameter(name='id', value='2005', format='int')
88 opObj11.addParameter(name='wintitle', value='CrossSpectraPlot_ShortPulse', format='str')
89 # opObj11.addParameter(name='exp_code', value='13', format='int')
90 opObj11.addParameter(name='xaxis', value='Velocity', format='str')
86 # opObj11 = procUnitConfObj1.addOperation(name='CrossSpectraPlot', optype='other')
87 # procUnitConfObj1.addParameter(name='pairsList', value='(0,1),(0,2),(1,2)', format='pairsList')
88 # opObj11.addParameter(name='id', value='2005', format='int')
89 # opObj11.addParameter(name='wintitle', value='CrossSpectraPlot_ShortPulse', format='str')
90 # # opObj11.addParameter(name='exp_code', value='13', format='int')
91 # opObj11.addParameter(name='xaxis', value='Velocity', format='str')
91 92 #opObj11.addParameter(name='xmin', value='-10', format='float')
92 93 #opObj11.addParameter(name='xmax', value='10', format='float')
93 94 #opObj11.addParameter(name='ymin', value='225', format='float')
@@ -99,14 +100,14 opObj11.addParameter(name='xaxis', value='Velocity', format='str')
99 100 # procUnitConfObj2.addParameter(name='pairsList', value='(0,1),(0,2),(1,2)', format='pairsList')
100 101
101 102 procUnitConfObj2 = controllerObj.addProcUnit(datatype='Parameters', inputId=procUnitConfObj1.getId())
102 opObj11 = procUnitConfObj2.addOperation(name='SpectralMoments', optype='other')
103 #opObj11 = procUnitConfObj2.addOperation(name='SpectralMoments', optype='other')
103 104 opObj22 = procUnitConfObj2.addOperation(name='FullSpectralAnalysis', optype='other')
104 opObj22.addParameter(name='SNRlimit', value='7', format='float')
105 opObj22.addParameter(name='SNRlimit', value='2', format='float')
105 106
106 107 opObj22 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other')
107 108 opObj22.addParameter(name='id', value='4', format='int')
108 109 opObj22.addParameter(name='wintitle', value='Wind Profiler', format='str')
109 opObj22.addParameter(name='save', value='1', format='bool')
110 #opObj22.addParameter(name='save', value='1', format='bool')
110 111 # opObj22.addParameter(name='figpath', value = '/home/erick/Pictures', format='str')
111 112
112 113 opObj22.addParameter(name='zmin', value='-20', format='int')
Show file before | Show file after | Hide comments
@@ -7,20 +7,20 if __name__ == '__main__':
7 7 desc = "Segundo Test"
8 8 filename = "schain.xml"
9 9
10 pathW='/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra'
11 figpath = '/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra'
10 pathW='/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra'
11 figpath = '/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra'
12 12
13 13 controllerObj = Project()
14 14
15 15 controllerObj.setup(id='191', name='test01', description=desc)
16 16
17 17 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
18 path='/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/',
18 path='/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA',
19 19 #path='/home/erick/Documents/Data/Claire_Data/raw',
20 startDate='2018/02/01',
21 endDate='2018/02/01',
22 startTime='12:00:00',
23 endTime='20:00:00',
20 startDate='2018/02/22',
21 endDate='2018/02/24',
22 startTime='00:00:00',
23 endTime='23:59:00',
24 24 online=0,
25 25 walk=1)
26 26
General Comments 0
You need to be logged in to leave comments. Login now