##// END OF EJS Templates
First Spectral Fitting and EW Drifts operative module inside Signal Chain TRUNK
Julio Valdez -
r513:19f921674eb5
parent child
Show More
@@ -0,0 +1,135
1 # DIAS 19 Y 20 FEB 2014
2 # Comprobacion de Resultados DBS con SA
3
4 import os, sys
5
6 path = os.path.split(os.getcwd())[0]
7 sys.path.append(path)
8
9 from controller import *
10
11 desc = "DBS Experiment Test"
12 filename = "DBStest.xml"
13
14 controllerObj = Project()
15
16 controllerObj.setup(id = '191', name='test01', description=desc)
17
18 #Experimentos
19
20 path = '/host/Jicamarca/EW_Drifts/d2012248'
21 pathFigure = '/home/propietario/workspace/Graficos/drifts'
22
23
24 path = "/home/soporte/Data/drifts"
25 pathFigure = '/home/soporte/workspace/Graficos/drifts/prueba'
26
27 xmin = 11.75
28 xmax = 14.75
29 #------------------------------------------------------------------------------------------------
30 readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
31 path=path,
32 startDate='2012/01/01',
33 endDate='2012/12/31',
34 startTime='00:00:00',
35 endTime='23:59:59',
36 online=0,
37 walk=1)
38
39 opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock')
40
41 #--------------------------------------------------------------------------------------------------
42
43 procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
44
45 opObj11 = procUnitConfObj0.addOperation(name='ProfileSelector', optype='other')
46 opObj11.addParameter(name='profileRangeList', value='0,127', format='intlist')
47
48 opObj11 = procUnitConfObj0.addOperation(name='filterByHeights')
49 opObj11.addParameter(name='window', value='3', format='int')
50
51 opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other')
52 # opObj11.addParameter(name='code', value='1,-1', format='floatlist')
53 # opObj11.addParameter(name='nCode', value='2', format='int')
54 # opObj11.addParameter(name='nBaud', value='1', format='int')
55
56 procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId())
57 procUnitConfObj1.addParameter(name='nFFTPoints', value='128', format='int')
58 procUnitConfObj1.addParameter(name='nProfiles', value='128', format='int')
59 procUnitConfObj1.addParameter(name='pairsList', value='(0,1),(2,3)', format='pairsList')#,(2,3)
60
61 opObj11 = procUnitConfObj1.addOperation(name='selectHeights')
62 # # opObj11.addParameter(name='minHei', value='320.0', format='float')
63 # # opObj11.addParameter(name='maxHei', value='350.0', format='float')
64 opObj11.addParameter(name='minHei', value='200.0', format='float')
65 opObj11.addParameter(name='maxHei', value='600.0', format='float')
66
67 opObj11 = procUnitConfObj1.addOperation(name='selectChannels')
68 opObj11.addParameter(name='channelList', value='0,1,2,3', format='intlist')
69
70 opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other')
71 opObj11.addParameter(name='timeInterval', value='300.0', format='float')
72
73 opObj13 = procUnitConfObj1.addOperation(name='removeDC')
74
75 # opObj14 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other')
76 # opObj14.addParameter(name='id', value='1', format='int')
77 # # opObj14.addParameter(name='wintitle', value='Con interf', format='str')
78 # opObj14.addParameter(name='save', value='1', format='bool')
79 # opObj14.addParameter(name='figpath', value=pathFigure, format='str')
80 # # opObj14.addParameter(name='zmin', value='5', format='int')
81 # opObj14.addParameter(name='zmax', value='30', format='int')
82 #
83 # opObj12 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other')
84 # opObj12.addParameter(name='id', value='2', format='int')
85 # opObj12.addParameter(name='wintitle', value='RTI Plot', format='str')
86 # opObj12.addParameter(name='save', value='1', format='bool')
87 # opObj12.addParameter(name='figpath', value = pathFigure, format='str')
88 # opObj12.addParameter(name='xmin', value=xmin, format='float')
89 # opObj12.addParameter(name='xmax', value=xmax, format='float')
90 # # opObj12.addParameter(name='zmin', value='5', format='int')
91 # opObj12.addParameter(name='zmax', value='30', format='int')
92
93 #--------------------------------------------------------------------------------------------------
94
95 procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId())
96 opObj20 = procUnitConfObj2.addOperation(name='SpectralFitting')
97 opObj20.addParameter(name='path', value='/home/soporte/workspace/RemoteSystemsTempFiles', format='str')
98 opObj20.addParameter(name='file', value='modelSpectralFitting', format='str')
99 opObj20.addParameter(name='groupList', value='(0,1),(2,3)',format='multiList')
100
101 opObj11 = procUnitConfObj2.addOperation(name='SpectralFittingPlot', optype='other')
102 opObj11.addParameter(name='id', value='3', format='int')
103 opObj11.addParameter(name='wintitle', value='DopplerPlot', format='str')
104 opObj11.addParameter(name='cutHeight', value='350', format='int')
105 opObj11.addParameter(name='fit', value='1', format='int')#1--True/include fit
106 opObj11.addParameter(name='save', value='1', format='bool')
107 opObj11.addParameter(name='figpath', value = pathFigure, format='str')
108
109 opObj11 = procUnitConfObj2.addOperation(name='EWDriftsEstimation', optype='other')
110 opObj11.addParameter(name='zenith', value='-3.80208,3.10658', format='floatlist')
111 opObj11.addParameter(name='zenithCorrection', value='0.183201', format='float')
112
113 opObj23 = procUnitConfObj2.addOperation(name='EWDriftsPlot', optype='other')
114 opObj23.addParameter(name='id', value='4', format='int')
115 opObj23.addParameter(name='wintitle', value='EW Drifts', format='str')
116 opObj23.addParameter(name='save', value='1', format='bool')
117 opObj23.addParameter(name='figpath', value = pathFigure, format='str')
118 opObj23.addParameter(name='zminZonal', value='-150', format='int')
119 opObj23.addParameter(name='zmaxZonal', value='150', format='int')
120 opObj23.addParameter(name='zminVertical', value='-30', format='float')
121 opObj23.addParameter(name='zmaxVertical', value='30', format='float')
122 opObj23.addParameter(name='SNR_1', value='1', format='bool')
123 opObj23.addParameter(name='SNRmax', value='5', format='int')
124 # opObj23.addParameter(name='SNRthresh', value='-50', format='float')
125 opObj23.addParameter(name='xmin', value=xmin, format='float')
126 opObj23.addParameter(name='xmax', value=xmax, format='float')
127 #--------------------------------------------------------------------------------------------------
128 print "Escribiendo el archivo XML"
129 controllerObj.writeXml(filename)
130 print "Leyendo el archivo XML"
131 controllerObj.readXml(filename)
132
133 controllerObj.createObjects()
134 controllerObj.connectObjects()
135 controllerObj.run() No newline at end of file
@@ -917,10 +917,18 class Correlation(JROData):
917 917
918 918 class Parameters(JROData):
919 919
920 #Information from previous data
921
920 922 inputUnit = None #Type of data to be processed
921 923
922 924 operation = None #Type of operation to parametrize
923 925
926 normFactor = None #Normalization Factor
927
928 groupList = None #List of Pairs, Groups, etc
929
930 #Parameters
931
924 932 data_param = None #Parameters obtained
925 933
926 934 data_pre = None #Data Pre Parametrization
@@ -933,17 +941,25 class Parameters(JROData):
933 941
934 942 SNR = None #Signal to Noise Ratio
935 943
936 pairsList = None #List of Pairs for Cross correlations or Cross spectrum
937
938 944 initUtcTime = None #Initial UTC time
939 945
940 946 paramInterval = None #Time interval to calculate Parameters in seconds
941 947
942 windsInterval = None #Time interval to calculate Winds in seconds
948 #Fitting
949
950 constants = None
951
952 error = None
953
954 library = None
955
956 #Output signal
957
958 outputInterval = None #Time interval to calculate output signal in seconds
959
960 data_output = None #Out signal
943 961
944 normFactor = None #Normalization Factor
945 962
946 winds = None #Wind estimations
947 963
948 964 def __init__(self):
949 965 '''
@@ -960,7 +976,7 class Parameters(JROData):
960 976 datatime = []
961 977
962 978 datatime.append(self.initUtcTime)
963 datatime.append(self.initUtcTime + self.windsInterval - 1)
979 datatime.append(self.initUtcTime + self.outputInterval - 1)
964 980
965 981 datatime = numpy.array(datatime)
966 982
@@ -455,7 +455,7 class WindProfilerPlot(Figure):
455 455 # y = dataOut.heightRange
456 456 y = dataOut.heightRange
457 457
458 z = dataOut.winds.copy()
458 z = dataOut.data_output.copy()
459 459 nplots = z.shape[0] #Number of wind dimensions estimated
460 460 nplotsw = nplots
461 461
@@ -771,4 +771,408 class ParametersPlot(Figure):
771 771 if x[1] >= self.axesList[0].xmax:
772 772 self.counter_imagwr = wr_period
773 773 self.__isConfig = False
774 self.figfile = None
775
776
777 class SpectralFittingPlot(Figure):
778
779 __isConfig = None
780 __nsubplots = None
781
782 WIDTHPROF = None
783 HEIGHTPROF = None
784 PREFIX = 'prm'
785
786
787 N = None
788 ippSeconds = None
789
790 def __init__(self):
791 self.__isConfig = False
792 self.__nsubplots = 1
793
794 self.WIDTH = 450
795 self.HEIGHT = 250
796 self.WIDTHPROF = 0
797 self.HEIGHTPROF = 0
798
799 def getSubplots(self):
800
801 ncol = int(numpy.sqrt(self.nplots)+0.9)
802 nrow = int(self.nplots*1./ncol + 0.9)
803
804 return nrow, ncol
805
806 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
807
808 showprofile = False
809 self.__showprofile = showprofile
810 self.nplots = nplots
811
812 ncolspan = 5
813 colspan = 4
814 if showprofile:
815 ncolspan = 5
816 colspan = 4
817 self.__nsubplots = 2
818
819 self.createFigure(id = id,
820 wintitle = wintitle,
821 widthplot = self.WIDTH + self.WIDTHPROF,
822 heightplot = self.HEIGHT + self.HEIGHTPROF,
823 show=show)
824
825 nrow, ncol = self.getSubplots()
826
827 counter = 0
828 for y in range(nrow):
829 for x in range(ncol):
830
831 if counter >= self.nplots:
832 break
833
834 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
835
836 if showprofile:
837 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
838
839 counter += 1
840
841 def run(self, dataOut, id, cutHeight=None, fit=False, wintitle="", channelList=None, showprofile=True,
842 xmin=None, xmax=None, ymin=None, ymax=None,
843 save=False, figpath='./', figfile=None, show=True):
844
845 """
846
847 Input:
848 dataOut :
849 id :
850 wintitle :
851 channelList :
852 showProfile :
853 xmin : None,
854 xmax : None,
855 zmin : None,
856 zmax : None
857 """
858
859 if cutHeight==None:
860 h=270
861 heightindex = numpy.abs(cutHeight - dataOut.heightList).argmin()
862 cutHeight = dataOut.heightList[heightindex]
863
864 factor = dataOut.normFactor
865 x = dataOut.abscissaRange[:-1]
866 #y = dataOut.getHeiRange()
867
868 z = dataOut.data_pre[:,:,heightindex]/factor
869 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
870 avg = numpy.average(z, axis=1)
871 listChannels = z.shape[0]
872
873 #Reconstruct Function
874 if fit==True:
875 groupArray = dataOut.groupList
876 listChannels = groupArray.reshape((groupArray.size))
877 listChannels.sort()
878 spcFitLine = numpy.zeros(z.shape)
879 constants = dataOut.constants
880
881 nGroups = groupArray.shape[0]
882 nChannels = groupArray.shape[1]
883 nProfiles = z.shape[1]
884
885 for f in range(nGroups):
886 groupChann = groupArray[f,:]
887 p = dataOut.data_param[f,:,heightindex]
888 # p = numpy.array([ 89.343967,0.14036615,0.17086219,18.89835291,1.58388365,1.55099167])
889 fitLineAux = dataOut.library.modelFunction(p, constants)*nProfiles
890 fitLineAux = fitLineAux.reshape((nChannels,nProfiles))
891 spcFitLine[groupChann,:] = fitLineAux
892 # spcFitLine = spcFitLine/factor
893
894 z = z[listChannels,:]
895 spcFitLine = spcFitLine[listChannels,:]
896 spcFitLinedB = 10*numpy.log10(spcFitLine)
897
898 zdB = 10*numpy.log10(z)
899 #thisDatetime = dataOut.datatime
900 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
901 title = wintitle + " Doppler Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
902 xlabel = "Velocity (m/s)"
903 ylabel = "Spectrum"
904
905 if not self.__isConfig:
906
907 nplots = listChannels.size
908
909 self.setup(id=id,
910 nplots=nplots,
911 wintitle=wintitle,
912 showprofile=showprofile,
913 show=show)
914
915 if xmin == None: xmin = numpy.nanmin(x)
916 if xmax == None: xmax = numpy.nanmax(x)
917 if ymin == None: ymin = numpy.nanmin(zdB)
918 if ymax == None: ymax = numpy.nanmax(zdB)+2
919
920 self.__isConfig = True
921
922 self.setWinTitle(title)
923 for i in range(self.nplots):
924 # title = "Channel %d: %4.2fdB" %(dataOut.channelList[i]+1, noisedB[i])
925 title = "Height %4.1f km\nChannel %d:" %(cutHeight, listChannels[i]+1)
926 axes = self.axesList[i*self.__nsubplots]
927 if fit == False:
928 axes.pline(x, zdB[i,:],
929 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
930 xlabel=xlabel, ylabel=ylabel, title=title
931 )
932 if fit == True:
933 fitline=spcFitLinedB[i,:]
934 y=numpy.vstack([zdB[i,:],fitline] )
935 legendlabels=['Data','Fitting']
936 axes.pmultilineyaxis(x, y,
937 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
938 xlabel=xlabel, ylabel=ylabel, title=title,
939 legendlabels=legendlabels, marker=None,
940 linestyle='solid', grid='both')
941
942 self.draw()
943
944 if save:
945 date = thisDatetime.strftime("%Y%m%d_%H%M%S")
946 if figfile == None:
947 figfile = self.getFilename(name = date)
948
949 self.saveFigure(figpath, figfile)
950
951
952 class EWDriftsPlot(Figure):
953
954 __isConfig = None
955 __nsubplots = None
956
957 WIDTHPROF = None
958 HEIGHTPROF = None
959 PREFIX = 'drift'
960
961 def __init__(self):
962
963 self.timerange = 2*60*60
964 self.isConfig = False
965 self.__nsubplots = 1
966
967 self.WIDTH = 800
968 self.HEIGHT = 150
969 self.WIDTHPROF = 120
970 self.HEIGHTPROF = 0
971 self.counter_imagwr = 0
972
973 self.PLOT_CODE = 0
974 self.FTP_WEI = None
975 self.EXP_CODE = None
976 self.SUB_EXP_CODE = None
977 self.PLOT_POS = None
978 self.tmin = None
979 self.tmax = None
980
981 self.xmin = None
982 self.xmax = None
983
984 self.figfile = None
985
986 def getSubplots(self):
987
988 ncol = 1
989 nrow = self.nplots
990
991 return nrow, ncol
992
993 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
994
995 self.__showprofile = showprofile
996 self.nplots = nplots
997
998 ncolspan = 1
999 colspan = 1
1000
1001 self.createFigure(id = id,
1002 wintitle = wintitle,
1003 widthplot = self.WIDTH + self.WIDTHPROF,
1004 heightplot = self.HEIGHT + self.HEIGHTPROF,
1005 show=show)
1006
1007 nrow, ncol = self.getSubplots()
1008
1009 counter = 0
1010 for y in range(nrow):
1011 if counter >= self.nplots:
1012 break
1013
1014 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
1015 counter += 1
1016
1017 def run(self, dataOut, id, wintitle="", channelList=None,
1018 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
1019 zmaxVertical = None, zminVertical = None, zmaxZonal = None, zminZonal = None,
1020 timerange=None, SNRthresh = -numpy.inf, SNRmin = None, SNRmax = None, SNR_1 = False,
1021 save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1022 server=None, folder=None, username=None, password=None,
1023 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1024 """
1025
1026 Input:
1027 dataOut :
1028 id :
1029 wintitle :
1030 channelList :
1031 showProfile :
1032 xmin : None,
1033 xmax : None,
1034 ymin : None,
1035 ymax : None,
1036 zmin : None,
1037 zmax : None
1038 """
1039
1040 if channelList == None:
1041 channelIndexList = dataOut.channelIndexList
1042 else:
1043 channelIndexList = []
1044 for channel in channelList:
1045 if channel not in dataOut.channelList:
1046 raise ValueError, "Channel %d is not in dataOut.channelList"
1047 channelIndexList.append(dataOut.channelList.index(channel))
1048
1049 if timerange != None:
1050 self.timerange = timerange
1051
1052 tmin = None
1053 tmax = None
1054
1055 x = dataOut.getTimeRange1()
1056 # y = dataOut.heightRange
1057 y = dataOut.heightList
1058
1059 z = dataOut.data_output
1060 nplots = z.shape[0] #Number of wind dimensions estimated
1061 nplotsw = nplots
1062
1063 #If there is a SNR function defined
1064 if dataOut.SNR != None:
1065 nplots += 1
1066 SNR = dataOut.SNR
1067
1068 if SNR_1:
1069 SNR += 1
1070
1071 SNRavg = numpy.average(SNR, axis=0)
1072
1073 SNRdB = 10*numpy.log10(SNR)
1074 SNRavgdB = 10*numpy.log10(SNRavg)
1075
1076 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
1077
1078 for i in range(nplotsw):
1079 z[i,ind] = numpy.nan
1080
1081
1082 showprofile = False
1083 # thisDatetime = dataOut.datatime
1084 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1])
1085 title = wintitle + " EW Drifts"
1086 xlabel = ""
1087 ylabel = "Height (Km)"
1088
1089 if not self.__isConfig:
1090
1091 self.setup(id=id,
1092 nplots=nplots,
1093 wintitle=wintitle,
1094 showprofile=showprofile,
1095 show=show)
1096
1097 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1098
1099 if ymin == None: ymin = numpy.nanmin(y)
1100 if ymax == None: ymax = numpy.nanmax(y)
1101
1102 if zmaxZonal == None: zmaxZonal = numpy.nanmax(abs(z[0,:]))
1103 if zminZonal == None: zminZonal = -zmaxZonal
1104 if zmaxVertical == None: zmaxVertical = numpy.nanmax(abs(z[1,:]))
1105 if zminVertical == None: zminVertical = -zmaxVertical
1106
1107 if dataOut.SNR != None:
1108 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
1109 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
1110
1111 self.FTP_WEI = ftp_wei
1112 self.EXP_CODE = exp_code
1113 self.SUB_EXP_CODE = sub_exp_code
1114 self.PLOT_POS = plot_pos
1115
1116 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1117 self.__isConfig = True
1118
1119
1120 self.setWinTitle(title)
1121
1122 if ((self.xmax - x[1]) < (x[1]-x[0])):
1123 x[1] = self.xmax
1124
1125 strWind = ['Zonal','Vertical']
1126 strCb = 'Velocity (m/s)'
1127 zmaxVector = [zmaxZonal, zmaxVertical]
1128 zminVector = [zminZonal, zminVertical]
1129
1130 for i in range(nplotsw):
1131
1132 title = "%s Drifts: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1133 axes = self.axesList[i*self.__nsubplots]
1134
1135 z1 = z[i,:].reshape((1,-1))
1136
1137 axes.pcolorbuffer(x, y, z1,
1138 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
1139 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1140 ticksize=9, cblabel=strCb, cbsize="1%", colormap="RdBu_r")
1141
1142 if dataOut.SNR != None:
1143 i += 1
1144 if SNR_1:
1145 title = "Signal Noise Ratio + 1 (SNR+1): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1146 else:
1147 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1148 axes = self.axesList[i*self.__nsubplots]
1149 SNRavgdB = SNRavgdB.reshape((1,-1))
1150
1151 axes.pcolorbuffer(x, y, SNRavgdB,
1152 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1153 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1154 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
1155
1156 self.draw()
1157
1158 if self.figfile == None:
1159 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
1160 self.figfile = self.getFilename(name = str_datetime)
1161
1162 if figpath != '':
1163
1164 self.counter_imagwr += 1
1165 if (self.counter_imagwr>=wr_period):
1166 # store png plot to local folder
1167 self.saveFigure(figpath, self.figfile)
1168 # store png plot to FTP server according to RT-Web format
1169 name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
1170 ftp_filename = os.path.join(figpath, name)
1171 self.saveFigure(figpath, ftp_filename)
1172
1173 self.counter_imagwr = 0
1174
1175 if x[1] >= self.axesList[0].xmax:
1176 self.counter_imagwr = wr_period
1177 self.__isConfig = False
774 1178 self.figfile = None No newline at end of file
@@ -316,7 +316,7 def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel=''
316 316 matplotlib.pyplot.ioff()
317 317
318 318 # lines = ax.plot(x, y.T, marker=marker,markersize=markersize,linestyle=linestyle)
319 lines = ax.plot(x, y.T, linestyle='None', marker='.', markersize=markersize)
319 lines = ax.plot(x, y.T, linestyle=linestyle, marker=marker, markersize=markersize)
320 320 leg = ax.legend(lines, legendlabels, loc='upper left', bbox_to_anchor=(1.01, 1.00), numpoints=1, handlelength=1.5, \
321 321 handletextpad=0.5, borderpad=0.5, labelspacing=0.5, borderaxespad=0.)
322 322
@@ -7,7 +7,9 from scipy import stats
7 7 import re
8 8 import datetime
9 9 import copy
10
10 import sys
11 import importlib
12 import itertools
11 13
12 14 from jroproc_base import ProcessingUnit, Operation
13 15 from model.data.jrodata import Parameters
@@ -60,7 +62,7 class ParametersProc(ProcessingUnit):
60 62
61 63 def run(self, nSeconds = None, nProfiles = None):
62 64
63 self.dataOut.flagNoData = True
65
64 66
65 67 if self.firstdatatime == None:
66 68 self.firstdatatime = self.dataIn.utctime
@@ -68,6 +70,7 class ParametersProc(ProcessingUnit):
68 70 #---------------------- Voltage Data ---------------------------
69 71
70 72 if self.dataIn.type == "Voltage":
73 self.dataOut.flagNoData = True
71 74 if nSeconds != None:
72 75 self.nSeconds = nSeconds
73 76 self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))
@@ -100,6 +103,7 class ParametersProc(ProcessingUnit):
100 103 self.dataOut.abscissaRange = self.dataIn.getVelRange(1)
101 104 self.dataOut.noise = self.dataIn.getNoise()
102 105 self.dataOut.normFactor = self.dataIn.normFactor
106 self.dataOut.flagNoData = False
103 107
104 108 #---------------------- Correlation Data ---------------------------
105 109
@@ -112,14 +116,14 class ParametersProc(ProcessingUnit):
112 116 self.dataOut.noise = self.dataIn.noise
113 117 self.dataOut.normFactor = self.dataIn.normFactor
114 118 self.dataOut.SNR = self.dataIn.SNR
115 self.dataOut.pairsList = self.dataIn.pairsList
119 self.dataOut.groupList = self.dataIn.pairsList
120 self.dataOut.flagNoData = False
116 121
117 122
118 123 self.__updateObjFromInput()
119 self.dataOut.flagNoData = False
120 124 self.firstdatatime = None
121 125 self.dataOut.initUtcTime = self.dataIn.ltctime
122 self.dataOut.windsInterval = self.dataIn.timeInterval
126 self.dataOut.outputInterval = self.dataIn.timeInterval
123 127
124 128 #------------------- Get Moments ----------------------------------
125 129 def GetMoments(self, channelList = None):
@@ -238,7 +242,7 class ParametersProc(ProcessingUnit):
238 242 self.dataOut.noise
239 243 self.dataOut.normFactor
240 244 self.dataOut.SNR
241 self.dataOut.pairsList
245 self.dataOut.groupList
242 246 self.dataOut.nChannels
243 247
244 248 Affected:
@@ -251,7 +255,7 class ParametersProc(ProcessingUnit):
251 255 absc = self.dataOut.abscissaRange[:-1]
252 256 noise = self.dataOut.noise
253 257 SNR = self.dataOut.SNR
254 pairsList = self.dataOut.pairsList
258 pairsList = self.dataOut.groupList
255 259 nChannels = self.dataOut.nChannels
256 260 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
257 261 self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
@@ -1109,6 +1113,137 class ParametersProc(ProcessingUnit):
1109 1113
1110 1114 return heights, error
1111 1115
1116 def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None):
1117
1118 '''
1119 Function GetMoments()
1120
1121 Input:
1122 Output:
1123 Variables modified:
1124 '''
1125 if path != None:
1126 sys.path.append(path)
1127 self.dataOut.library = importlib.import_module(file)
1128
1129 #To be inserted as a parameter
1130 groupArray = numpy.array(groupList)
1131 # groupArray = numpy.array([[0,1],[2,3]])
1132 self.dataOut.groupList = groupArray
1133
1134 nGroups = groupArray.shape[0]
1135 nChannels = self.dataIn.nChannels
1136 nHeights=self.dataIn.heightList.size
1137
1138 #Parameters Array
1139 self.dataOut.data_param = None
1140
1141 #Set constants
1142 constants = self.dataOut.library.setConstants(self.dataIn)
1143 self.dataOut.constants = constants
1144 M = self.dataIn.normFactor
1145 N = self.dataIn.nFFTPoints
1146 ippSeconds = self.dataIn.ippSeconds
1147 K = self.dataIn.nIncohInt
1148 pairsArray = numpy.array(self.dataIn.pairsList)
1149
1150 #List of possible combinations
1151 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
1152 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
1153
1154 if getSNR:
1155 listChannels = groupArray.reshape((groupArray.size))
1156 listChannels.sort()
1157 noise = self.dataIn.getNoise()
1158 self.dataOut.SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
1159
1160 for i in range(nGroups):
1161 coord = groupArray[i,:]
1162
1163 #Input data array
1164 data = self.dataIn.data_spc[coord,:,:]/(M*N)
1165 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
1166
1167 #Cross Spectra data array for Covariance Matrixes
1168 ind = 0
1169 for pairs in listComb:
1170 pairsSel = numpy.array([coord[x],coord[y]])
1171 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
1172 ind += 1
1173 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
1174 dataCross = dataCross**2/K
1175
1176 for h in range(nHeights):
1177 # print self.dataOut.heightList[h]
1178
1179 #Input
1180 d = data[:,h]
1181
1182 #Covariance Matrix
1183 D = numpy.diag(d**2/K)
1184 ind = 0
1185 for pairs in listComb:
1186 #Coordinates in Covariance Matrix
1187 x = pairs[0]
1188 y = pairs[1]
1189 #Channel Index
1190 S12 = dataCross[ind,:,h]
1191 D12 = numpy.diag(S12)
1192 #Completing Covariance Matrix with Cross Spectras
1193 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
1194 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
1195 ind += 1
1196 Dinv=numpy.linalg.inv(D)
1197 L=numpy.linalg.cholesky(Dinv)
1198 LT=L.T
1199
1200 dp = numpy.dot(LT,d)
1201
1202 #Initial values
1203 data_spc = self.dataIn.data_spc[coord,:,h]
1204 p0 = self.dataOut.library.initialValuesFunction(data_spc, constants)
1205
1206 #Least Squares
1207 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
1208 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
1209 #Chi square error
1210 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
1211 # error0 = 0
1212 #Error with Jacobian
1213 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
1214 #Save
1215 if self.dataOut.data_param == None:
1216 self.dataOut.data_param = numpy.zeros((nGroups, minp.size, nHeights))*numpy.nan
1217 self.dataOut.error = numpy.zeros((nGroups, error1.size + 1, nHeights))*numpy.nan
1218
1219 self.dataOut.error[i,:,h] = numpy.hstack((error0,error1))
1220 self.dataOut.data_param[i,:,h] = minp
1221 return
1222
1223
1224 def __residFunction(self, p, dp, LT, constants):
1225
1226 fm = self.dataOut.library.modelFunction(p, constants)
1227 fmp=numpy.dot(LT,fm)
1228
1229 return dp-fmp
1230
1231 def __getSNR(self, z, noise):
1232
1233 avg = numpy.average(z, axis=1)
1234 SNR = (avg.T-noise)/noise
1235 SNR = SNR.T
1236 return SNR
1237
1238 def __chisq(p,chindex,hindex):
1239 #similar to Resid but calculates CHI**2
1240 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
1241 dp=numpy.dot(LT,d)
1242 fmp=numpy.dot(LT,fm)
1243 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
1244 return chisq
1245
1246
1112 1247
1113 1248 class WindProfiler(Operation):
1114 1249
@@ -1359,12 +1494,12 class WindProfiler(Operation):
1359 1494 winds = correctFactor*winds
1360 1495 return winds
1361 1496
1362 def __checkTime(self, currentTime, paramInterval, windsInterval):
1497 def __checkTime(self, currentTime, paramInterval, outputInterval):
1363 1498
1364 1499 dataTime = currentTime + paramInterval
1365 1500 deltaTime = dataTime - self.__initime
1366 1501
1367 if deltaTime >= windsInterval or deltaTime < 0:
1502 if deltaTime >= outputInterval or deltaTime < 0:
1368 1503 self.__dataReady = True
1369 1504 return
1370 1505
@@ -1462,7 +1597,7 class WindProfiler(Operation):
1462 1597 theta_y = theta_y[arrayChannel]
1463 1598
1464 1599 velRadial0 = param[:,1,:] #Radial velocity
1465 dataOut.winds, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function
1600 dataOut.data_output, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function
1466 1601
1467 1602 elif technique == 'SA':
1468 1603
@@ -1483,12 +1618,12 class WindProfiler(Operation):
1483 1618
1484 1619 tau = dataOut.data_param
1485 1620 _lambda = dataOut.C/dataOut.frequency
1486 pairsList = dataOut.pairsList
1621 pairsList = dataOut.groupList
1487 1622 nChannels = dataOut.nChannels
1488 1623
1489 dataOut.winds = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1624 dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1490 1625 dataOut.initUtcTime = dataOut.ltctime
1491 dataOut.windsInterval = dataOut.timeInterval
1626 dataOut.outputInterval = dataOut.timeInterval
1492 1627
1493 1628 elif technique == 'Meteors':
1494 1629 dataOut.flagNoData = True
@@ -1511,7 +1646,7 class WindProfiler(Operation):
1511 1646 hmax = kwargs['hmax']
1512 1647 else: hmax = 110
1513 1648
1514 dataOut.windsInterval = nHours*3600
1649 dataOut.outputInterval = nHours*3600
1515 1650
1516 1651 if self.__isConfig == False:
1517 1652 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
@@ -1526,14 +1661,89 class WindProfiler(Operation):
1526 1661 else:
1527 1662 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1528 1663
1529 self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.windsInterval) #Check if the buffer is ready
1664 self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1530 1665
1531 1666 if self.__dataReady:
1532 1667 dataOut.initUtcTime = self.__initime
1533 self.__initime = self.__initime + dataOut.windsInterval #to erase time offset
1668 self.__initime = self.__initime + dataOut.outputInterval #to erase time offset
1534 1669
1535 dataOut.winds, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1670 dataOut.data_output, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1536 1671 dataOut.flagNoData = False
1537 1672 self.__buffer = None
1538 1673
1539 return No newline at end of file
1674 return
1675
1676 class EWDriftsEstimation(Operation):
1677
1678
1679 def __init__(self):
1680 Operation.__init__(self)
1681
1682 def __correctValues(self, heiRang, phi, velRadial, SNR):
1683 listPhi = phi.tolist()
1684 maxid = listPhi.index(max(listPhi))
1685 minid = listPhi.index(min(listPhi))
1686
1687 rango = range(len(phi))
1688 # rango = numpy.delete(rango,maxid)
1689
1690 heiRang1 = heiRang*math.cos(phi[maxid])
1691 heiRangAux = heiRang*math.cos(phi[minid])
1692 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1693 heiRang1 = numpy.delete(heiRang1,indOut)
1694
1695 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1696 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1697
1698 for i in rango:
1699 x = heiRang*math.cos(phi[i])
1700 y1 = velRadial[i,:]
1701 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1702
1703 x1 = heiRang1
1704 y11 = f1(x1)
1705
1706 y2 = SNR[i,:]
1707 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1708 y21 = f2(x1)
1709
1710 velRadial1[i,:] = y11
1711 SNR1[i,:] = y21
1712
1713 return heiRang1, velRadial1, SNR1
1714
1715 def run(self, dataOut, zenith, zenithCorrection):
1716 heiRang = dataOut.heightList
1717 velRadial = dataOut.data_param[:,3,:]
1718 SNR = dataOut.SNR
1719
1720 zenith = numpy.array(zenith)
1721 zenith -= zenithCorrection
1722 zenith *= numpy.pi/180
1723
1724 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
1725
1726 alp = zenith[0]
1727 bet = zenith[1]
1728
1729 w_w = velRadial1[0,:]
1730 w_e = velRadial1[1,:]
1731
1732 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
1733 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
1734
1735 winds = numpy.vstack((u,w))
1736
1737 dataOut.heightList = heiRang1
1738 dataOut.data_output = winds
1739 dataOut.SNR = SNR1
1740
1741 dataOut.initUtcTime = dataOut.ltctime
1742 dataOut.outputInterval = dataOut.timeInterval
1743 return
1744
1745
1746
1747
1748
1749 No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now