@@ -0,0 +1,188 | |||
|
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 | path = os.path.split(path)[0] | |
|
8 | ||
|
9 | sys.path.insert(0, path) | |
|
10 | ||
|
11 | from schainpy.controller import Project | |
|
12 | ||
|
13 | desc = "JASMET Experiment Test" | |
|
14 | filename = "JASMETtest.xml" | |
|
15 | ||
|
16 | controllerObj = Project() | |
|
17 | ||
|
18 | controllerObj.setup(id = '191', name='test01', description=desc) | |
|
19 | ||
|
20 | #Experimentos | |
|
21 | ||
|
22 | #2014051 20 Feb 2014 | |
|
23 | # path = '/home/soporte/Data/JASMET/JASMET_30/2014106' | |
|
24 | # pathFigure = '/home/soporte/workspace/Graficos/JASMET/prueba1' | |
|
25 | ||
|
26 | remotefolder = "/home/wmaster/graficos" | |
|
27 | path = '/media/joscanoa/84A65E64A65E5730/soporte/Data/JASMET/JASMET_30/' | |
|
28 | ||
|
29 | pathfile1 = os.path.join(os.environ['HOME'],'Pictures/JASMET30/meteor') | |
|
30 | pathfile2 = os.path.join(os.environ['HOME'],'Pictures/JASMET30/wind') | |
|
31 | pathfile3 = os.path.join(os.environ['HOME'],'Pictures/JASMET30/phase') | |
|
32 | ||
|
33 | pathfig = os.path.join(os.environ['HOME'],'Pictures/JASMET30/graph') | |
|
34 | ||
|
35 | startTime = '00:00:00' | |
|
36 | endTime = '23:59:59' | |
|
37 | xmin ='17.0' | |
|
38 | xmax = '34.0' | |
|
39 | ||
|
40 | #------------------------------------------------------------------------------------------------ | |
|
41 | readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader', | |
|
42 | path=path, | |
|
43 | startDate='2014/04/15', | |
|
44 | endDate='2014/04/16', | |
|
45 | startTime=startTime, | |
|
46 | endTime=endTime, | |
|
47 | online=0, | |
|
48 | delay=5, | |
|
49 | walk=1) | |
|
50 | ||
|
51 | opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock') | |
|
52 | ||
|
53 | ||
|
54 | #-------------------------------------------------------------------------------------------------- | |
|
55 | ||
|
56 | procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId()) | |
|
57 | ||
|
58 | opObj00 = procUnitConfObj0.addOperation(name='selectChannels') | |
|
59 | opObj00.addParameter(name='channelList', value='0, 1, 2, 3, 4', format='intlist') | |
|
60 | ||
|
61 | opObj01 = procUnitConfObj0.addOperation(name='setRadarFrequency') | |
|
62 | # opObj01.addParameter(name='frequency', value='30.e6', format='float') | |
|
63 | opObj01.addParameter(name='frequency', value='50.e6', format='float') | |
|
64 | ||
|
65 | #opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other') | |
|
66 | #-------------------------------------------------------------------------------------------------- | |
|
67 | ||
|
68 | procUnitConfObj1 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj0.getId()) | |
|
69 | procUnitConfObj1.addParameter(name='nSeconds', value='100', format='int') | |
|
70 | ||
|
71 | opObj10 = procUnitConfObj1.addOperation(name='MeteorDetection') | |
|
72 | # opObj10.addParameter(name='predefinedPhaseShifts', value='-89.5, 41.5, 0.0, -138.0, -85.5', format='floatlist') | |
|
73 | opObj10.addParameter(name='predefinedPhaseShifts', value='0, 0, 0, 0, 0', format='floatlist') | |
|
74 | opObj10.addParameter(name='cohDetection', value='0', format='bool') | |
|
75 | opObj10.addParameter(name='noise_multiple', value='4', format='int') | |
|
76 | opObj10.addParameter(name='SNRThresh', value='5', format='float') | |
|
77 | opObj10.addParameter(name='phaseThresh', value='20', format='float') | |
|
78 | opObj10.addParameter(name='azimuth', value='45', format='float') | |
|
79 | opObj10.addParameter(name='hmin', value='68', format='float') | |
|
80 | opObj10.addParameter(name='hmax', value='112', format='float') | |
|
81 | opObj10.addParameter(name='saveAll', value='1', format='bool') | |
|
82 | ||
|
83 | opObj12 = procUnitConfObj1.addOperation(name='HDF5Writer', optype='other') | |
|
84 | opObj12.addParameter(name='path', value=pathfile1) | |
|
85 | opObj12.addParameter(name='blocksPerFile', value='1000', format='int') | |
|
86 | opObj12.addParameter(name='metadataList',value='type,inputUnit,heightList,paramInterval',format='list') | |
|
87 | opObj12.addParameter(name='dataList',value='data_param',format='list') | |
|
88 | opObj12.addParameter(name='mode',value='0',format='int') | |
|
89 | #Tiene que ser de 3 dimensiones, append en lugar de aumentar una dimension | |
|
90 | ||
|
91 | opObj13 = procUnitConfObj1.addOperation(name='SkyMapPlot', optype='other') | |
|
92 | opObj13.addParameter(name='id', value='1', format='int') | |
|
93 | opObj13.addParameter(name='wintitle', value='Sky Map', format='str') | |
|
94 | opObj13.addParameter(name='save', value='1', format='bool') | |
|
95 | opObj13.addParameter(name='figpath', value=pathfig, format='str') | |
|
96 | opObj13.addParameter(name='ftp', value='1', format='int') | |
|
97 | opObj13.addParameter(name='exp_code', value='15', format='int') | |
|
98 | opObj13.addParameter(name='sub_exp_code', value='1', format='int') | |
|
99 | ||
|
100 | ||
|
101 | #-------------------------------------------------------------------------------------------------- | |
|
102 | procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId()) | |
|
103 | ||
|
104 | opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other') | |
|
105 | opObj22.addParameter(name='technique', value='Meteors', format='str') | |
|
106 | opObj22.addParameter(name='nHours', value='0.25', format='float') | |
|
107 | opObj22.addParameter(name='SNRThresh', value='12.0', format='float') | |
|
108 | ||
|
109 | opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other') | |
|
110 | opObj23.addParameter(name='id', value='2', format='int') | |
|
111 | opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str') | |
|
112 | opObj23.addParameter(name='save', value='1', format='bool') | |
|
113 | opObj23.addParameter(name='figpath', value = pathfig, format='str') | |
|
114 | opObj23.addParameter(name='zmin', value='-120', format='int') | |
|
115 | opObj23.addParameter(name='zmax', value='120', format='int') | |
|
116 | # opObj12.addParameter(name='zmin_ver', value='-0.8', format='float') | |
|
117 | # opObj12.addParameter(name='zmax_ver', value='0.8', format='float') | |
|
118 | # opObj23.addParameter(name='SNRmin', value='-10', format='int') | |
|
119 | # opObj23.addParameter(name='SNRmax', value='60', format='int') | |
|
120 | # opObj23.addParameter(name='SNRthresh', value='0', format='float') | |
|
121 | opObj23.addParameter(name='xmin', value=xmin, format='float') | |
|
122 | opObj23.addParameter(name='xmax', value=xmax, format='float') | |
|
123 | opObj23.addParameter(name='ftp', value='1', format='int') | |
|
124 | opObj23.addParameter(name='exp_code', value='15', format='int') | |
|
125 | opObj23.addParameter(name='sub_exp_code', value='1', format='int') | |
|
126 | ||
|
127 | opObj24 = procUnitConfObj2.addOperation(name='HDF5Writer', optype='other') | |
|
128 | opObj24.addParameter(name='path', value=pathfile2) | |
|
129 | opObj24.addParameter(name='blocksPerFile', value='1000', format='int') | |
|
130 | opObj24.addParameter(name='metadataList',value='type,inputUnit,outputInterval',format='list') | |
|
131 | opObj24.addParameter(name='dataList',value='data_output,utctime',format='list') | |
|
132 | ||
|
133 | #-------------------------------------------------------------------------------------------------- | |
|
134 | procUnitConfObj3 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId()) | |
|
135 | ||
|
136 | ||
|
137 | opObj31 = procUnitConfObj3.addOperation(name='PhaseCalibration', optype='other') | |
|
138 | opObj31.addParameter(name='nHours', value='0.25', format='float') | |
|
139 | opObj31.addParameter(name='distances', value='-15, -15, 12, 12', format='intlist') | |
|
140 | # opObj31.addParameter(name='distances', value='-25, -25, 20, 20', format='intlist') | |
|
141 | opObj31.addParameter(name='pairs', value='(0,3),(1,2)', format='pairslist') | |
|
142 | opObj31.addParameter(name='hmin', value='68', format='float') | |
|
143 | opObj31.addParameter(name='hmax', value='112', format='float') | |
|
144 | ||
|
145 | opObj32 = procUnitConfObj3.addOperation(name='PhasePlot', optype='other') | |
|
146 | opObj32.addParameter(name='id', value='201', format='int') | |
|
147 | opObj32.addParameter(name='wintitle', value='PhaseCalibration', format='str') | |
|
148 | # opObj32.addParameter(name='timerange', value='300', format='int') | |
|
149 | opObj32.addParameter(name='xmin', value=xmin, format='float') | |
|
150 | opObj32.addParameter(name='xmax', value=xmax, format='float') | |
|
151 | # opObj32.addParameter(name='xmin', value='19', format='float') | |
|
152 | # opObj32.addParameter(name='xmax', value='20', format='float') | |
|
153 | opObj32.addParameter(name='ymin', value='-180', format='float') | |
|
154 | opObj32.addParameter(name='ymax', value='180', format='float') | |
|
155 | opObj32.addParameter(name='figpath', value=pathfig, format='str') | |
|
156 | opObj32.addParameter(name='ftp', value='1', format='int') | |
|
157 | opObj32.addParameter(name='exp_code', value='15', format='int') | |
|
158 | opObj32.addParameter(name='sub_exp_code', value='1', format='int') | |
|
159 | ||
|
160 | opObj33 = procUnitConfObj3.addOperation(name='HDF5Writer', optype='other') | |
|
161 | opObj33.addParameter(name='path', value=pathfile3) | |
|
162 | opObj33.addParameter(name='blocksPerFile', value='1000', format='int') | |
|
163 | opObj33.addParameter(name='metadataList',value='type,inputUnit,outputInterval',format='list') | |
|
164 | opObj33.addParameter(name='dataList',value='data_output,utctime',format='list') | |
|
165 | # opObj25.addParameter(name='mode',value='1,0,0',format='intlist') | |
|
166 | ||
|
167 | #-------------------------------------------------------------------------------------------------- | |
|
168 | ||
|
169 | procUnitConfObj4 = controllerObj.addProcUnit(name='SendToServer') | |
|
170 | procUnitConfObj4.addParameter(name='server', value='jro-app.igp.gob.pe', format='str') | |
|
171 | procUnitConfObj4.addParameter(name='username', value='wmaster', format='str') | |
|
172 | procUnitConfObj4.addParameter(name='password', value='mst2010vhf', format='str') | |
|
173 | procUnitConfObj4.addParameter(name='localfolder', value=pathfig, format='str') | |
|
174 | procUnitConfObj4.addParameter(name='remotefolder', value=remotefolder, format='str') | |
|
175 | procUnitConfObj4.addParameter(name='ext', value='.png', format='str') | |
|
176 | procUnitConfObj4.addParameter(name='period', value=20, format='int') | |
|
177 | procUnitConfObj4.addParameter(name='protocol', value='ftp', format='str') | |
|
178 | ||
|
179 | #-------------------------------------------------------------------------------------------------- | |
|
180 | ||
|
181 | print "Escribiendo el archivo XML" | |
|
182 | controllerObj.writeXml(filename) | |
|
183 | print "Leyendo el archivo XML" | |
|
184 | controllerObj.readXml(filename) | |
|
185 | ||
|
186 | controllerObj.createObjects() | |
|
187 | controllerObj.connectObjects() | |
|
188 | controllerObj.run() No newline at end of file |
@@ -0,0 +1,56 | |||
|
1 | """ | |
|
2 | Se debe verficar que el disco de datos se encuentra montado en el sistema | |
|
3 | """ | |
|
4 | import os, sys | |
|
5 | ||
|
6 | path = os.path.split(os.getcwd())[0] | |
|
7 | path = os.path.split(path)[0] | |
|
8 | ||
|
9 | sys.path.insert(0, path) | |
|
10 | ||
|
11 | from schainpy.controller import Project | |
|
12 | ||
|
13 | desc = "Meteor Experiment Test" | |
|
14 | filename = "meteor20130812.xml" | |
|
15 | ||
|
16 | controllerObj = Project() | |
|
17 | controllerObj.setup(id = '191', name='meteor_test01', description=desc) | |
|
18 | ||
|
19 | path = '/home/dsuarez/.gvfs/data on 10.10.20.13/Jasmet50' | |
|
20 | ||
|
21 | ||
|
22 | readUnitConfObj = controllerObj.addReadUnit(datatype='Voltage', | |
|
23 | path=path, | |
|
24 | startDate='2014/04/16', | |
|
25 | endDate='2014/04/16', | |
|
26 | startTime='00:00:00', | |
|
27 | endTime='23:59:59', | |
|
28 | online=0, | |
|
29 | walk=1) | |
|
30 | ||
|
31 | opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock') | |
|
32 | ||
|
33 | procUnitConfObj0 = controllerObj.addProcUnit(datatype='Voltage', inputId=readUnitConfObj.getId()) | |
|
34 | ||
|
35 | ||
|
36 | opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other') | |
|
37 | ||
|
38 | ||
|
39 | opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other') | |
|
40 | opObj11.addParameter(name='n', value='2', format='int') | |
|
41 | ||
|
42 | opObj11 = procUnitConfObj0.addOperation(name='VoltageWriter', optype='other') | |
|
43 | opObj11.addParameter(name='path', value='/home/jasmet/jasmet30_abril') | |
|
44 | opObj11.addParameter(name='blocksPerFile', value='100', format='int') | |
|
45 | opObj11.addParameter(name='profilesPerBlock', value='200', format='int') | |
|
46 | ||
|
47 | ||
|
48 | ||
|
49 | print "Escribiendo el archivo XML" | |
|
50 | controllerObj.writeXml(filename) | |
|
51 | print "Leyendo el archivo XML" | |
|
52 | controllerObj.readXml(filename) | |
|
53 | ||
|
54 | controllerObj.createObjects() | |
|
55 | controllerObj.connectObjects() | |
|
56 | controllerObj.run() |
@@ -1069,13 +1069,13 class Parameters(JROData): | |||
|
1069 | 1069 | |
|
1070 | 1070 | data_SNR = None #Signal to Noise Ratio |
|
1071 | 1071 | |
|
1072 | heightRange = None #Heights | |
|
1072 | # heightRange = None #Heights | |
|
1073 | 1073 | |
|
1074 |
abscissa |
|
|
1074 | abscissaList = None #Abscissa, can be velocities, lags or time | |
|
1075 | 1075 | |
|
1076 | 1076 | noise = None #Noise Potency |
|
1077 | 1077 | |
|
1078 |
|
|
|
1078 | utctimeInit = None #Initial UTC time | |
|
1079 | 1079 | |
|
1080 | 1080 | paramInterval = None #Time interval to calculate Parameters in seconds |
|
1081 | 1081 | |
@@ -1093,6 +1093,8 class Parameters(JROData): | |||
|
1093 | 1093 | |
|
1094 | 1094 | data_output = None #Out signal |
|
1095 | 1095 | |
|
1096 | ||
|
1097 | ||
|
1096 | 1098 | def __init__(self): |
|
1097 | 1099 | ''' |
|
1098 | 1100 | Constructor |
@@ -1107,8 +1109,15 class Parameters(JROData): | |||
|
1107 | 1109 | |
|
1108 | 1110 | datatime = [] |
|
1109 | 1111 | |
|
1110 | datatime.append(self.ltctime) | |
|
1111 | datatime.append(self.ltctime + self.outputInterval - 1) | |
|
1112 | if self.useLocalTime: | |
|
1113 | time1 = self.utctimeInit - self.timeZone*60 | |
|
1114 | else: | |
|
1115 | time1 = utctimeInit | |
|
1116 | ||
|
1117 | # datatime.append(self.utctimeInit) | |
|
1118 | # datatime.append(self.utctimeInit + self.outputInterval) | |
|
1119 | datatime.append(time1) | |
|
1120 | datatime.append(time1 + self.outputInterval) | |
|
1112 | 1121 | |
|
1113 | 1122 | datatime = numpy.array(datatime) |
|
1114 | 1123 |
@@ -84,9 +84,9 class Figure(Operation): | |||
|
84 | 84 | #raise ValueError, "(timerange) or (xmin & xmax) should be defined" |
|
85 | 85 | |
|
86 | 86 | if timerange != None: |
|
87 | txmin = x[0] - x[0] % min(timerange/10, 10*60) | |
|
87 | txmin = x[0] #- x[0] % min(timerange/10, 10*60) | |
|
88 | 88 | else: |
|
89 | txmin = x[0] - x[0] % 10*60 | |
|
89 | txmin = x[0] #- x[0] % 10*60 | |
|
90 | 90 | |
|
91 | 91 | thisdatetime = datetime.datetime.utcfromtimestamp(txmin) |
|
92 | 92 | thisdate = datetime.datetime.combine(thisdatetime.date(), datetime.time(0,0,0)) |
@@ -198,7 +198,7 class SkyMapPlot(Figure): | |||
|
198 | 198 | |
|
199 | 199 | WIDTHPROF = None |
|
200 | 200 | HEIGHTPROF = None |
|
201 |
PREFIX = 'p |
|
|
201 | PREFIX = 'mmap' | |
|
202 | 202 | |
|
203 | 203 | def __init__(self): |
|
204 | 204 | |
@@ -213,7 +213,7 class SkyMapPlot(Figure): | |||
|
213 | 213 | self.HEIGHTPROF = 0 |
|
214 | 214 | self.counter_imagwr = 0 |
|
215 | 215 | |
|
216 | self.PLOT_CODE = SKYMAP_CODE | |
|
216 | self.PLOT_CODE = MSKYMAP_CODE | |
|
217 | 217 | |
|
218 | 218 | self.FTP_WEI = None |
|
219 | 219 | self.EXP_CODE = None |
@@ -269,7 +269,7 class SkyMapPlot(Figure): | |||
|
269 | 269 | zmax : None |
|
270 | 270 | """ |
|
271 | 271 | |
|
272 | arrayParameters = dataOut.data_param | |
|
272 | arrayParameters = dataOut.data_param[0,:] | |
|
273 | 273 | error = arrayParameters[:,-1] |
|
274 | 274 | indValid = numpy.where(error == 0)[0] |
|
275 | 275 | finalMeteor = arrayParameters[indValid,:] |
@@ -325,6 +325,7 class SkyMapPlot(Figure): | |||
|
325 | 325 | wr_period=wr_period, |
|
326 | 326 | thisDatetime=thisDatetime) |
|
327 | 327 | |
|
328 | ||
|
328 | 329 | class WindProfilerPlot(Figure): |
|
329 | 330 | |
|
330 | 331 | __isConfig = None |
@@ -336,7 +337,7 class WindProfilerPlot(Figure): | |||
|
336 | 337 | |
|
337 | 338 | def __init__(self): |
|
338 | 339 | |
|
339 |
self.timerange = |
|
|
340 | self.timerange = None | |
|
340 | 341 | self.__isConfig = False |
|
341 | 342 | self.__nsubplots = 1 |
|
342 | 343 | |
@@ -391,7 +392,7 class WindProfilerPlot(Figure): | |||
|
391 | 392 | self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1) |
|
392 | 393 | counter += 1 |
|
393 | 394 | |
|
394 | def run(self, dataOut, id, wintitle="", channelList=None, | |
|
395 | def run(self, dataOut, id, wintitle="", channelList=None, showprofile='False', | |
|
395 | 396 | xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, |
|
396 | 397 | zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None, |
|
397 | 398 | timerange=None, SNRthresh = None, |
@@ -423,11 +424,11 class WindProfilerPlot(Figure): | |||
|
423 | 424 | raise ValueError, "Channel %d is not in dataOut.channelList" |
|
424 | 425 | channelIndexList.append(dataOut.channelList.index(channel)) |
|
425 | 426 | |
|
426 | if timerange != None: | |
|
427 | self.timerange = timerange | |
|
428 | ||
|
429 | tmin = None | |
|
430 | tmax = None | |
|
427 | # if timerange != None: | |
|
428 | # self.timerange = timerange | |
|
429 | # | |
|
430 | # tmin = None | |
|
431 | # tmax = None | |
|
431 | 432 | |
|
432 | 433 | x = dataOut.getTimeRange1() |
|
433 | 434 | # y = dataOut.heightList |
@@ -453,7 +454,7 class WindProfilerPlot(Figure): | |||
|
453 | 454 | z[i,ind] = numpy.nan |
|
454 | 455 | |
|
455 | 456 | |
|
456 | showprofile = False | |
|
457 | # showprofile = False | |
|
457 | 458 | # thisDatetime = dataOut.datatime |
|
458 | 459 | thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) |
|
459 | 460 | title = wintitle + "Wind" |
@@ -468,6 +469,9 class WindProfilerPlot(Figure): | |||
|
468 | 469 | showprofile=showprofile, |
|
469 | 470 | show=show) |
|
470 | 471 | |
|
472 | if timerange != None: | |
|
473 | self.timerange = timerange | |
|
474 | ||
|
471 | 475 | self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange) |
|
472 | 476 | |
|
473 | 477 | if ymin == None: ymin = numpy.nanmin(y) |
@@ -485,6 +489,7 class WindProfilerPlot(Figure): | |||
|
485 | 489 | if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB) |
|
486 | 490 | if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB) |
|
487 | 491 | |
|
492 | ||
|
488 | 493 | self.FTP_WEI = ftp_wei |
|
489 | 494 | self.EXP_CODE = exp_code |
|
490 | 495 | self.SUB_EXP_CODE = sub_exp_code |
@@ -492,7 +497,7 class WindProfilerPlot(Figure): | |||
|
492 | 497 |
|
|
493 | 498 | self.name = thisDatetime.strftime("%Y%m%d_%H%M%S") |
|
494 | 499 | self.__isConfig = True |
|
495 | ||
|
500 | self.figfile = figfile | |
|
496 | 501 |
|
|
497 | 502 | self.setWinTitle(title) |
|
498 | 503 |
|
@@ -620,7 +625,7 class ParametersPlot(Figure): | |||
|
620 | 625 | def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False, |
|
621 | 626 | xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,timerange=None, |
|
622 | 627 | parameterIndex = None, onlyPositive = False, |
|
623 | SNRthresh = -numpy.inf, SNR = True, SNRmin = None, SNRmax = None, | |
|
628 | SNRthresh = -numpy.inf, SNR = True, SNRmin = None, SNRmax = None, onlySNR = False, | |
|
624 | 629 | DOP = True, |
|
625 | 630 | zlabel = "", parameterName = "", parameterObject = "data_param", |
|
626 | 631 | save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True, |
@@ -695,6 +700,8 class ParametersPlot(Figure): | |||
|
695 | 700 | xlabel = "" |
|
696 | 701 | ylabel = "Range (Km)" |
|
697 | 702 | |
|
703 | if (SNR and not onlySNR): nplots = 2*nplots | |
|
704 | ||
|
698 | 705 | if onlyPositive: |
|
699 | 706 | colormap = "jet" |
|
700 | 707 | zmin = 0 |
@@ -735,8 +742,22 class ParametersPlot(Figure): | |||
|
735 | 742 | |
|
736 | 743 | for i in range(nchan): |
|
737 | 744 | |
|
745 | if (SNR and not onlySNR): j = 2*i | |
|
746 | else: j = i | |
|
747 | ||
|
738 | 748 | j = nGraphsByChannel*i |
|
739 | 749 | |
|
750 | if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)): | |
|
751 | title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith) | |
|
752 | ||
|
753 | if not onlySNR: | |
|
754 | axes = self.axesList[j*self.__nsubplots] | |
|
755 | z1 = z[i,:].reshape((1,-1)) | |
|
756 | axes.pcolorbuffer(x, y, z1, | |
|
757 | xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, | |
|
758 | xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap, | |
|
759 | ticksize=9, cblabel=zlabel, cbsize="1%") | |
|
760 | ||
|
740 | 761 | if DOP: |
|
741 | 762 | title = "%s Channel %d: %s" %(parameterName, channelIndexList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) |
|
742 | 763 | |
@@ -751,7 +772,12 class ParametersPlot(Figure): | |||
|
751 | 772 | |
|
752 | 773 | if SNR: |
|
753 | 774 | title = "Channel %d Signal Noise Ratio (SNR): %s" %(channelIndexList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) |
|
775 | axes = self.axesList[(j)*self.__nsubplots] | |
|
776 | if not onlySNR: | |
|
777 | axes = self.axesList[(j + 1)*self.__nsubplots] | |
|
778 | ||
|
754 | 779 | axes = self.axesList[(j + nGraphsByChannel-1)] |
|
780 | ||
|
755 | 781 | z1 = SNRdB[i,:].reshape((1,-1)) |
|
756 | 782 | axes.pcolorbuffer(x, y, z1, |
|
757 | 783 | xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, |
@@ -1155,6 +1181,178 class EWDriftsPlot(Figure): | |||
|
1155 | 1181 | self.__isConfig = False |
|
1156 | 1182 | self.figfile = None |
|
1157 | 1183 | |
|
1184 | ||
|
1185 | ||
|
1186 | ||
|
1187 | class PhasePlot(Figure): | |
|
1188 | ||
|
1189 | __isConfig = None | |
|
1190 | __nsubplots = None | |
|
1191 | ||
|
1192 | PREFIX = 'mphase' | |
|
1193 | ||
|
1194 | def __init__(self): | |
|
1195 | ||
|
1196 | self.timerange = 24*60*60 | |
|
1197 | self.__isConfig = False | |
|
1198 | self.__nsubplots = 1 | |
|
1199 | self.counter_imagwr = 0 | |
|
1200 | self.WIDTH = 600 | |
|
1201 | self.HEIGHT = 300 | |
|
1202 | self.WIDTHPROF = 120 | |
|
1203 | self.HEIGHTPROF = 0 | |
|
1204 | self.xdata = None | |
|
1205 | self.ydata = None | |
|
1206 | ||
|
1207 | self.PLOT_CODE = MPHASE_CODE | |
|
1208 | ||
|
1209 | self.FTP_WEI = None | |
|
1210 | self.EXP_CODE = None | |
|
1211 | self.SUB_EXP_CODE = None | |
|
1212 | self.PLOT_POS = None | |
|
1213 | ||
|
1214 | ||
|
1215 | self.filename_phase = None | |
|
1216 | ||
|
1217 | self.figfile = None | |
|
1218 | ||
|
1219 | def getSubplots(self): | |
|
1220 | ||
|
1221 | ncol = 1 | |
|
1222 | nrow = 1 | |
|
1223 | ||
|
1224 | return nrow, ncol | |
|
1225 | ||
|
1226 | def setup(self, id, nplots, wintitle, showprofile=True, show=True): | |
|
1227 | ||
|
1228 | self.__showprofile = showprofile | |
|
1229 | self.nplots = nplots | |
|
1230 | ||
|
1231 | ncolspan = 7 | |
|
1232 | colspan = 6 | |
|
1233 | self.__nsubplots = 2 | |
|
1234 | ||
|
1235 | self.createFigure(id = id, | |
|
1236 | wintitle = wintitle, | |
|
1237 | widthplot = self.WIDTH+self.WIDTHPROF, | |
|
1238 | heightplot = self.HEIGHT+self.HEIGHTPROF, | |
|
1239 | show=show) | |
|
1240 | ||
|
1241 | nrow, ncol = self.getSubplots() | |
|
1242 | ||
|
1243 | self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1) | |
|
1244 | ||
|
1245 | ||
|
1246 | def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True', | |
|
1247 | xmin=None, xmax=None, ymin=None, ymax=None, | |
|
1248 | timerange=None, | |
|
1249 | save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1, | |
|
1250 | server=None, folder=None, username=None, password=None, | |
|
1251 | ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): | |
|
1252 | ||
|
1253 | ||
|
1254 | if timerange != None: | |
|
1255 | self.timerange = timerange | |
|
1256 | ||
|
1257 | tmin = None | |
|
1258 | tmax = None | |
|
1259 | x = dataOut.getTimeRange1() | |
|
1260 | y = dataOut.getHeiRange() | |
|
1261 | ||
|
1262 | ||
|
1263 | #thisDatetime = dataOut.datatime | |
|
1264 | thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) | |
|
1265 | title = wintitle + " Phase of Beacon Signal" # : %s" %(thisDatetime.strftime("%d-%b-%Y")) | |
|
1266 | xlabel = "Local Time" | |
|
1267 | ylabel = "Phase" | |
|
1268 | ||
|
1269 | ||
|
1270 | #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList))) | |
|
1271 | phase_beacon = dataOut.data_output | |
|
1272 | ||
|
1273 | ||
|
1274 | if not self.__isConfig: | |
|
1275 | ||
|
1276 | self.nplots = phase_beacon.size | |
|
1277 | ||
|
1278 | self.setup(id=id, | |
|
1279 | nplots=self.nplots, | |
|
1280 | wintitle=wintitle, | |
|
1281 | showprofile=showprofile, | |
|
1282 | show=show) | |
|
1283 | ||
|
1284 | tmin, tmax = self.getTimeLim(x, xmin, xmax) | |
|
1285 | if ymin == None: ymin = numpy.nanmin(phase_beacon) - 10.0 | |
|
1286 | if ymax == None: ymax = numpy.nanmax(phase_beacon) + 10.0 | |
|
1287 | ||
|
1288 | self.FTP_WEI = ftp_wei | |
|
1289 | self.EXP_CODE = exp_code | |
|
1290 | self.SUB_EXP_CODE = sub_exp_code | |
|
1291 | self.PLOT_POS = plot_pos | |
|
1292 | ||
|
1293 | self.name = thisDatetime.strftime("%Y%m%d_%H%M%S") | |
|
1294 | self.__isConfig = True | |
|
1295 | self.figfile = figfile | |
|
1296 | self.xdata = numpy.array([]) | |
|
1297 | self.ydata = numpy.array([]) | |
|
1298 | ||
|
1299 | #open file beacon phase | |
|
1300 | path = '%s%03d' %(self.PREFIX, self.id) | |
|
1301 | beacon_file = os.path.join(path,'%s.txt'%self.name) | |
|
1302 | self.filename_phase = os.path.join(figpath,beacon_file) | |
|
1303 | #self.save_phase(self.filename_phase) | |
|
1304 | ||
|
1305 | ||
|
1306 | #store data beacon phase | |
|
1307 | #self.save_data(self.filename_phase, phase_beacon, thisDatetime) | |
|
1308 | ||
|
1309 | self.setWinTitle(title) | |
|
1310 | ||
|
1311 | ||
|
1312 | title = "Phase Offset %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) | |
|
1313 | ||
|
1314 | legendlabels = ["phase %d"%(chan) for chan in numpy.arange(self.nplots)] | |
|
1315 | ||
|
1316 | axes = self.axesList[0] | |
|
1317 | ||
|
1318 | self.xdata = numpy.hstack((self.xdata, x[0:1])) | |
|
1319 | ||
|
1320 | if len(self.ydata)==0: | |
|
1321 | self.ydata = phase_beacon.reshape(-1,1) | |
|
1322 | else: | |
|
1323 | self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1))) | |
|
1324 | ||
|
1325 | ||
|
1326 | axes.pmultilineyaxis(x=self.xdata, y=self.ydata, | |
|
1327 | xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax, | |
|
1328 | xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid", | |
|
1329 | XAxisAsTime=True, grid='both' | |
|
1330 | ) | |
|
1331 | ||
|
1332 | self.draw() | |
|
1333 | ||
|
1334 | if x[1] >= self.axesList[0].xmax: | |
|
1335 | self.counter_imagwr = wr_period | |
|
1336 | del self.xdata | |
|
1337 | del self.ydata | |
|
1338 | self.__isConfig = False | |
|
1339 | ||
|
1340 | if self.figfile == None: | |
|
1341 | str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") | |
|
1342 | self.figfile = self.getFilename(name = str_datetime) | |
|
1343 | ||
|
1344 | if figpath != '': | |
|
1345 | self.counter_imagwr += 1 | |
|
1346 | if (self.counter_imagwr>=wr_period): | |
|
1347 | # store png plot to local folder | |
|
1348 | self.saveFigure(figpath, self.figfile) | |
|
1349 | # store png plot to FTP server according to RT-Web format | |
|
1350 | name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS) | |
|
1351 | ftp_filename = os.path.join(figpath, name) | |
|
1352 | self.saveFigure(figpath, ftp_filename) | |
|
1353 | self.counter_imagwr = 0 | |
|
1354 | self.figfile = None | |
|
1355 | ||
|
1158 | 1356 | self.save(figpath=figpath, |
|
1159 | 1357 | figfile=figfile, |
|
1160 | 1358 | save=save, |
@@ -395,9 +395,11 def createPolar(ax, x, y, | |||
|
395 | 395 | # ax.set_rmax(90) |
|
396 | 396 | ax.set_ylim(0,90) |
|
397 | 397 | ax.set_yticks(numpy.arange(0,90,20)) |
|
398 | ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11') | |
|
398 | # ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11') | |
|
399 | # ax.text(0, 50, ylabel, rotation='vertical', va ='center', ha = 'left' ,size='11') | |
|
399 | 400 | # ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical') |
|
400 | printLabels(ax, xlabel, '', title) | |
|
401 | ax.yaxis.labelpad = 230 | |
|
402 | printLabels(ax, xlabel, ylabel, title) | |
|
401 | 403 | iplot = ax.lines[-1] |
|
402 | 404 | |
|
403 | 405 | if '0.' in matplotlib.__version__[0:2]: |
@@ -423,7 +425,7 def polar(iplot, x, y, xlabel='', ylabel='', title=''): | |||
|
423 | 425 | ax = iplot.get_axes() |
|
424 | 426 | |
|
425 | 427 | # ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11') |
|
426 |
printLabels(ax, xlabel, |
|
|
428 | printLabels(ax, xlabel, ylabel, title) | |
|
427 | 429 | |
|
428 | 430 | set_linedata(ax, x, y, idline=0) |
|
429 | 431 |
@@ -18,10 +18,11 NOISE_CODE = 17 | |||
|
18 | 18 | BEACON_CODE = 18 |
|
19 | 19 | |
|
20 | 20 | #USED IN jroplot_parameters.py |
|
21 | ||
|
22 | MOMENTS_CODE = 20 | |
|
23 | SKYMAP_CODE = 21 | |
|
24 | 21 | WIND_CODE = 22 |
|
25 |
|
|
|
26 |
S |
|
|
27 | EWDRIFT_CODE = 25 | |
|
22 | MSKYMAP_CODE = 23 | |
|
23 | MPHASE_CODE = 24 | |
|
24 | ||
|
25 | MOMENTS_CODE = 25 | |
|
26 | PARMS_CODE = 26 | |
|
27 | SPECFIT_CODE = 27 | |
|
28 | EWDRIFT_CODE = 28 |
@@ -3,6 +3,7 import time | |||
|
3 | 3 | import os |
|
4 | 4 | import h5py |
|
5 | 5 | import re |
|
6 | import tables | |
|
6 | 7 | |
|
7 | 8 | from schainpy.model.data.jrodata import * |
|
8 | 9 | from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation |
@@ -156,7 +157,7 class HDF5Reader(ProcessingUnit): | |||
|
156 | 157 | for thisPath in os.listdir(single_path): |
|
157 | 158 | if not os.path.isdir(os.path.join(single_path,thisPath)): |
|
158 | 159 | continue |
|
159 |
if not is |
|
|
160 | if not isDoyFolder(thisPath): | |
|
160 | 161 | continue |
|
161 | 162 | |
|
162 | 163 | dirList.append(thisPath) |
@@ -549,6 +550,8 class HDF5Writer(Operation): | |||
|
549 | 550 | |
|
550 | 551 | metaFile = None |
|
551 | 552 | |
|
553 | filename = None | |
|
554 | ||
|
552 | 555 | path = None |
|
553 | 556 | |
|
554 | 557 | setFile = None |
@@ -589,6 +592,8 class HDF5Writer(Operation): | |||
|
589 | 592 | |
|
590 | 593 | nDims = None #Number Dimensions in each dataset |
|
591 | 594 | |
|
595 | nDimsForDs = None | |
|
596 | ||
|
592 | 597 | def __init__(self): |
|
593 | 598 | |
|
594 | 599 | Operation.__init__(self) |
@@ -675,16 +680,26 class HDF5Writer(Operation): | |||
|
675 | 680 | setFile = self.setFile |
|
676 | 681 | |
|
677 | 682 | timeTuple = time.localtime(self.dataOut.utctime) |
|
683 | ||
|
678 | 684 | subfolder = '' |
|
685 | fullpath = os.path.join( path, subfolder ) | |
|
679 | 686 | |
|
687 | if not( os.path.exists(fullpath) ): | |
|
688 | os.mkdir(fullpath) | |
|
689 | setFile = -1 #inicializo mi contador de seteo | |
|
690 | ||
|
691 | subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday) | |
|
680 | 692 | fullpath = os.path.join( path, subfolder ) |
|
693 | ||
|
681 | 694 | if not( os.path.exists(fullpath) ): |
|
682 | 695 | os.mkdir(fullpath) |
|
683 | 696 | setFile = -1 #inicializo mi contador de seteo |
|
697 | ||
|
684 | 698 | else: |
|
685 | 699 | filesList = os.listdir( fullpath ) |
|
686 | if len( filesList ) > 0: | |
|
687 | 700 |
|
|
701 | if len( filesList ) > 0: | |
|
702 | filesList = [k for k in filesList if 'M' in k] | |
|
688 | 703 | filen = filesList[-1] |
|
689 | 704 | # el filename debera tener el siguiente formato |
|
690 | 705 | # 0 1234 567 89A BCDE (hex) |
@@ -727,18 +742,14 class HDF5Writer(Operation): | |||
|
727 | 742 | setFile = self.setFile |
|
728 | 743 | mode = self.mode |
|
729 | 744 | |
|
730 | if self.fp != None: | |
|
731 | self.fp.close() | |
|
732 | ||
|
733 | 745 | timeTuple = time.localtime(self.dataOut.utctime) |
|
734 | 746 | subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday) |
|
735 | 747 | |
|
736 | 748 | fullpath = os.path.join( path, subfolder ) |
|
737 | if not( os.path.exists(fullpath) ): | |
|
738 |
|
|
|
739 | setFile = -1 #inicializo mi contador de seteo | |
|
740 | else: | |
|
749 | ||
|
750 | if os.path.exists(fullpath): | |
|
741 | 751 | filesList = os.listdir( fullpath ) |
|
752 | filesList = [k for k in filesList if 'D' in k] | |
|
742 | 753 | if len( filesList ) > 0: |
|
743 | 754 | filesList = sorted( filesList, key=str.lower ) |
|
744 | 755 | filen = filesList[-1] |
@@ -771,17 +782,22 class HDF5Writer(Operation): | |||
|
771 | 782 | |
|
772 | 783 | ds = [] |
|
773 | 784 | data = [] |
|
785 | nDimsForDs = [] | |
|
774 | 786 | |
|
775 | 787 | nDatas = numpy.zeros(len(self.dataList)) |
|
776 | 788 | nDims = self.arrayDim[:,0] |
|
777 | 789 | |
|
790 | nDim1 = self.arrayDim[:,2] | |
|
791 | nDim0 = self.arrayDim[:,3] | |
|
792 | ||
|
778 | 793 | for i in range(len(self.dataList)): |
|
779 | 794 | |
|
780 | 795 | if nDims[i]==1: |
|
781 |
ds0 = grp.create_dataset(self.dataList[i], (1,1), maxshape=(1, |
|
|
796 | # ds0 = grp.create_dataset(self.dataList[i], (1,1), maxshape=(1,self.blocksPerFile) , chunks = True, dtype='S20') | |
|
797 | ds0 = grp.create_dataset(self.dataList[i], (1,1), maxshape=(1,self.blocksPerFile) , chunks = True, dtype=numpy.float64) | |
|
782 | 798 | ds.append(ds0) |
|
783 | 799 | data.append([]) |
|
784 | ||
|
800 | nDimsForDs.append(nDims[i]) | |
|
785 | 801 | else: |
|
786 | 802 | |
|
787 | 803 | if mode[i]==0: |
@@ -800,20 +816,22 class HDF5Writer(Operation): | |||
|
800 | 816 | tableName = strMode + str(j) |
|
801 | 817 | |
|
802 | 818 | if nDims[i] == 3: |
|
803 |
ds0 = grp0.create_dataset(tableName, (1, |
|
|
819 | ds0 = grp0.create_dataset(tableName, (nDim1[i],nDim0[i],1) , data = numpy.zeros((nDim1[i],nDim0[i],1)) ,maxshape=(None,nDim0[i],None), chunks=True) | |
|
804 | 820 | else: |
|
805 |
ds0 = grp0.create_dataset(tableName, (1,1) , maxshape=(None, |
|
|
821 | ds0 = grp0.create_dataset(tableName, (1,nDim0[i]), data = numpy.zeros((1,nDim0[i])) , maxshape=(None,nDim0[i]), chunks=True) | |
|
806 | 822 | |
|
807 | 823 | ds.append(ds0) |
|
808 | 824 | data.append([]) |
|
809 | ||
|
825 | nDimsForDs.append(nDims[i]) | |
|
810 | 826 | self.nDatas = nDatas |
|
811 | 827 | self.nDims = nDims |
|
812 | ||
|
828 | self.nDimsForDs = nDimsForDs | |
|
813 | 829 | #Saving variables |
|
814 | 830 | print 'Writing the file: %s'%filename |
|
831 | self.filename = filename | |
|
815 | 832 | self.fp = fp |
|
816 | 833 | self.grp = grp |
|
834 | self.grp.attrs.modify('nRecords', 1) | |
|
817 | 835 | self.ds = ds |
|
818 | 836 | self.data = data |
|
819 | 837 | |
@@ -823,13 +841,65 class HDF5Writer(Operation): | |||
|
823 | 841 | return |
|
824 | 842 | |
|
825 | 843 | def putData(self): |
|
826 | self.setBlock() | |
|
827 | self.writeBlock() | |
|
844 | ||
|
845 | if not self.firsttime: | |
|
846 | self.fp.flush() | |
|
847 | self.fp.close() | |
|
848 | self.readBlock() | |
|
828 | 849 | |
|
829 | 850 | if self.blockIndex == self.blocksPerFile: |
|
851 | ||
|
830 | 852 |
self.setNextFile() |
|
853 | ||
|
854 | self.setBlock() | |
|
855 | self.writeBlock() | |
|
856 | ||
|
857 | ||
|
831 | 858 | return |
|
832 | 859 | |
|
860 | def readBlock(self): | |
|
861 | ||
|
862 | ''' | |
|
863 | data Array configured | |
|
864 | ||
|
865 | ||
|
866 | self.data | |
|
867 | ''' | |
|
868 | ds = self.ds | |
|
869 | #Setting HDF5 File | |
|
870 | fp = h5py.File(self.filename,'r+') | |
|
871 | grp = fp["Data"] | |
|
872 | ind = 0 | |
|
873 | ||
|
874 | # grp.attrs['blocksPerFile'] = 0 | |
|
875 | for i in range(len(self.dataList)): | |
|
876 | ||
|
877 | if self.nDims[i]==1: | |
|
878 | ds0 = grp[self.dataList[i]] | |
|
879 | ds[ind] = ds0 | |
|
880 | ind += 1 | |
|
881 | else: | |
|
882 | if self.mode[i]==0: | |
|
883 | strMode = "channel" | |
|
884 | else: | |
|
885 | strMode = "param" | |
|
886 | ||
|
887 | grp0 = grp[self.dataList[i]] | |
|
888 | ||
|
889 | for j in range(int(self.nDatas[i])): | |
|
890 | tableName = strMode + str(j) | |
|
891 | ds0 = grp0[tableName] | |
|
892 | ds[ind] = ds0 | |
|
893 | ind += 1 | |
|
894 | ||
|
895 | ||
|
896 | self.fp = fp | |
|
897 | self.grp = grp | |
|
898 | self.ds = ds | |
|
899 | ||
|
900 | return | |
|
901 | ||
|
902 | ||
|
833 | 903 | def setBlock(self): |
|
834 | 904 | ''' |
|
835 | 905 | data Array configured |
@@ -848,11 +918,11 class HDF5Writer(Operation): | |||
|
848 | 918 | dataAux = getattr(self.dataOut,self.dataList[i]) |
|
849 | 919 | |
|
850 | 920 | if nDims[i] == 1: |
|
851 | data[ind] = numpy.array([str(dataAux)]).reshape((1,1)) | |
|
852 | if not self.firsttime: | |
|
853 | data[ind] = numpy.hstack((self.ds[ind][:], self.data[ind])) | |
|
921 | # data[ind] = numpy.array([str(dataAux)]).reshape((1,1)) | |
|
922 | data[ind] = dataAux | |
|
923 | # if not self.firsttime: | |
|
924 | # data[ind] = numpy.hstack((self.ds[ind][:], self.data[ind])) | |
|
854 | 925 | ind += 1 |
|
855 | ||
|
856 | 926 | else: |
|
857 | 927 | for j in range(int(nDatas[i])): |
|
858 | 928 | if (mode[i] == 0) or (nDims[i] == 2): #In case division per channel or Dimensions is only 1 |
@@ -860,17 +930,17 class HDF5Writer(Operation): | |||
|
860 | 930 | else: |
|
861 | 931 | data[ind] = dataAux[:,j,:] |
|
862 | 932 | |
|
863 | if nDims[i] == 3: | |
|
864 | data[ind] = data[ind].reshape((data[ind].shape[0],data[ind].shape[1],1)) | |
|
933 | # if nDims[i] == 3: | |
|
934 | # data[ind] = data[ind].reshape((data[ind].shape[0],data[ind].shape[1],1)) | |
|
865 | 935 | |
|
866 | if not self.firsttime: | |
|
867 | data[ind] = numpy.dstack((self.ds[ind][:], data[ind])) | |
|
936 | # if not self.firsttime: | |
|
937 | # data[ind] = numpy.dstack((self.ds[ind][:], data[ind])) | |
|
868 | 938 | |
|
869 | else: | |
|
870 | data[ind] = data[ind].reshape((1,data[ind].shape[0])) | |
|
939 | # else: | |
|
940 | # data[ind] = data[ind].reshape((1,data[ind].shape[0])) | |
|
871 | 941 | |
|
872 | if not self.firsttime: | |
|
873 | data[ind] = numpy.vstack((self.ds[ind][:], data[ind])) | |
|
942 | # if not self.firsttime: | |
|
943 | # data[ind] = numpy.vstack((self.ds[ind][:], data[ind])) | |
|
874 | 944 | ind += 1 |
|
875 | 945 |
|
|
876 | 946 | self.data = data |
@@ -881,13 +951,45 class HDF5Writer(Operation): | |||
|
881 | 951 | Saves the block in the HDF5 file |
|
882 | 952 | ''' |
|
883 | 953 | for i in range(len(self.ds)): |
|
954 | if self.firsttime: | |
|
955 | # self.ds[i].resize(self.data[i].shape) | |
|
956 | # self.ds[i][self.blockIndex,:] = self.data[i] | |
|
957 | if type(self.data[i]) == numpy.ndarray: | |
|
958 | nDims1 = len(self.ds[i].shape) | |
|
959 | ||
|
960 | if nDims1 == 3: | |
|
961 | self.data[i] = self.data[i].reshape((self.data[i].shape[0],self.data[i].shape[1],1)) | |
|
962 | ||
|
884 | 963 |
self.ds[i].resize(self.data[i].shape) |
|
885 | 964 | self.ds[i][:] = self.data[i] |
|
965 | else: | |
|
966 | if self.nDimsForDs[i] == 1: | |
|
967 | self.ds[i].resize((self.ds[i].shape[0], self.ds[i].shape[1] + 1)) | |
|
968 | self.ds[i][0,-1] = self.data[i] | |
|
969 | elif self.nDimsForDs[i] == 2: | |
|
970 | self.ds[i].resize((self.ds[i].shape[0] + 1,self.ds[i].shape[1])) | |
|
971 | self.ds[i][self.blockIndex,:] = self.data[i] | |
|
972 | elif self.nDimsForDs[i] == 3: | |
|
973 | ||
|
974 | dataShape = self.data[i].shape | |
|
975 | dsShape = self.ds[i].shape | |
|
976 | ||
|
977 | if dataShape[0]==dsShape[0]: | |
|
978 | self.ds[i].resize((self.ds[i].shape[0],self.ds[i].shape[1],self.ds[i].shape[2]+1)) | |
|
979 | self.ds[i][:,:,-1] = self.data[i] | |
|
980 | else: | |
|
981 | self.ds[i].resize((self.ds[i].shape[0] + dataShape[0],self.ds[i].shape[1],self.ds[i].shape[2])) | |
|
982 | self.ds[i][dsShape[0]:,:,0] = self.data[i] | |
|
983 | # self.ds[i].append(self.data[i]) | |
|
984 | # self.fp.flush() | |
|
985 | # if not self.firsttime: | |
|
986 | # self.fp.root.Data._v_attrs.nRecords = self.blockIndex | |
|
886 | 987 | |
|
887 | self.blockIndex += 1 | |
|
888 | ||
|
889 | self.grp.attrs.modify('nRecords', self.blockIndex) | |
|
988 | # if self.firsttime: | |
|
989 | # self.fp.close() | |
|
990 | # self.readBlock2() | |
|
890 | 991 | |
|
992 | self.blockIndex += 1 | |
|
891 | 993 | self.firsttime = False |
|
892 | 994 | return |
|
893 | 995 |
@@ -239,13 +239,34 class ParametersProc(ProcessingUnit): | |||
|
239 | 239 | return moments |
|
240 | 240 | |
|
241 | 241 | #------------------ Get SA Parameters -------------------------- |
|
242 | ||
|
242 | 243 | def GetSAParameters(self): |
|
243 |
|
|
|
244 | crossdata = self.dataIn.data_cspc | |
|
245 |
|
|
|
244 | pairslist = self.dataOut.groupList | |
|
245 | num_pairs = len(pairslist) | |
|
246 | ||
|
247 | vel = self.dataOut.abscissaList | |
|
248 | spectra = self.dataOut.data_pre | |
|
249 | cspectra = self.dataIn.data_cspc | |
|
250 | delta_v = vel[1] - vel[0] | |
|
251 | ||
|
252 | #Calculating the power spectrum | |
|
253 | spc_pow = numpy.sum(spectra, 3)*delta_v | |
|
254 | #Normalizing Spectra | |
|
255 | norm_spectra = spectra/spc_pow | |
|
256 | #Calculating the norm_spectra at peak | |
|
257 | max_spectra = numpy.max(norm_spectra, 3) | |
|
246 | 258 | |
|
259 | #Normalizing Cross Spectra | |
|
260 | norm_cspectra = numpy.zeros(cspectra.shape) | |
|
247 | 261 | |
|
262 | for i in range(num_chan): | |
|
263 | norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:]) | |
|
248 | 264 | |
|
265 | max_cspectra = numpy.max(norm_cspectra,2) | |
|
266 | max_cspectra_index = numpy.argmax(norm_cspectra, 2) | |
|
267 | ||
|
268 | for i in range(num_pairs): | |
|
269 | cspc_par[i,:,:] = __calculateMoments(norm_cspectra) | |
|
249 | 270 | #------------------- Get Lags ---------------------------------- |
|
250 | 271 | |
|
251 | 272 | def GetLags(self): |
@@ -339,8 +360,8 class ParametersProc(ProcessingUnit): | |||
|
339 | 360 | return phase |
|
340 | 361 | #------------------- Detect Meteors ------------------------------ |
|
341 | 362 | |
|
342 |
def |
|
|
343 | predefinedPhaseShifts = None, centerReceiverIndex = 2, | |
|
363 | def MeteorDetection(self, hei_ref = None, tauindex = 0, | |
|
364 | predefinedPhaseShifts = None, centerReceiverIndex = 2, saveAll = False, | |
|
344 | 365 | cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25, |
|
345 | 366 | noise_timeStep = 4, noise_multiple = 4, |
|
346 | 367 | multDet_timeLimit = 1, multDet_rangeLimit = 3, |
@@ -430,14 +451,21 class ParametersProc(ProcessingUnit): | |||
|
430 | 451 | |
|
431 | 452 | if predefinedPhaseShifts != None: |
|
432 | 453 | hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180 |
|
433 |
|
|
|
454 | ||
|
455 | elif beaconPhaseShifts: | |
|
434 | 456 | #get hardware phase shifts using beacon signal |
|
435 | 457 | hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10) |
|
436 | 458 | hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0) |
|
437 | 459 | |
|
460 | else: | |
|
461 | hardwarePhaseShifts = numpy.zeros(5) | |
|
462 | ||
|
463 | ||
|
438 | 464 | voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex') |
|
439 | 465 | for i in range(self.dataOut.data_pre.shape[0]): |
|
440 | 466 | voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i]) |
|
467 | ||
|
468 | ||
|
441 | 469 | #******************END OF REMOVING HARDWARE PHASE DIFFERENCES********* |
|
442 | 470 | |
|
443 | 471 | #Remove DC |
@@ -497,29 +525,40 class ParametersProc(ProcessingUnit): | |||
|
497 | 525 | listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval) |
|
498 | 526 | |
|
499 | 527 | if len(listMeteors4) > 0: |
|
500 | #Setting New Array | |
|
501 | date = repr(self.dataOut.datatime) | |
|
502 | arrayMeteors4, arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang) | |
|
503 | 528 | |
|
504 | #Calculate AOA (Error N 3, 4) | |
|
505 | #JONES ET AL. 1998 | |
|
506 | AOAthresh = numpy.pi/8 | |
|
507 | error = arrayParameters[:,-1] | |
|
508 | phases = -arrayMeteors4[:,9:13] | |
|
509 | 529 | pairsList = [] |
|
510 |
pair |
|
|
511 |
pair |
|
|
512 | arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth) | |
|
530 | pairx = (0,3) | |
|
531 | pairy = (1,2) | |
|
532 | pairsList.append(pairx) | |
|
533 | pairsList.append(pairy) | |
|
513 | 534 | |
|
514 | #Calculate Heights (Error N 13 and 14) | |
|
515 | error = arrayParameters[:,-1] | |
|
516 | Ranges = arrayParameters[:,2] | |
|
517 | zenith = arrayParameters[:,5] | |
|
518 | arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax) | |
|
535 | #Setting New Array | |
|
536 | date = repr(self.dataOut.datatime) | |
|
537 | arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang) | |
|
538 | ||
|
539 | meteorOps = MeteorOperations() | |
|
540 | jph = numpy.array([0,0,0,0]) | |
|
541 | h = (hmin,hmax) | |
|
542 | arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, jph) | |
|
543 | ||
|
544 | # #Calculate AOA (Error N 3, 4) | |
|
545 | # #JONES ET AL. 1998 | |
|
546 | # error = arrayParameters[:,-1] | |
|
547 | # AOAthresh = numpy.pi/8 | |
|
548 | # phases = -arrayParameters[:,9:13] | |
|
549 | # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth) | |
|
550 | # | |
|
551 | # #Calculate Heights (Error N 13 and 14) | |
|
552 | # error = arrayParameters[:,-1] | |
|
553 | # Ranges = arrayParameters[:,2] | |
|
554 | # zenith = arrayParameters[:,5] | |
|
555 | # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax) | |
|
556 | # error = arrayParameters[:,-1] | |
|
519 | 557 | #********************* END OF PARAMETERS CALCULATION ************************** |
|
520 | 558 | |
|
521 |
#***************************+ |
|
|
522 | self.dataOut.data_param = arrayParameters | |
|
559 | #***************************+ PASS DATA TO NEXT STEP ********************** | |
|
560 | arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1])) | |
|
561 | self.dataOut.data_param = arrayFinal | |
|
523 | 562 | |
|
524 | 563 | return |
|
525 | 564 | |
@@ -948,9 +987,8 class ParametersProc(ProcessingUnit): | |||
|
948 | 987 | listMeteors1 = [] |
|
949 | 988 | |
|
950 | 989 | for i in range(len(listMeteors)): |
|
951 | meteor = listMeteors[i] | |
|
952 | meteorAux = numpy.hstack((meteor[:-1], 0, 0, meteor[-1])) | |
|
953 | if meteor[-1] == 0: | |
|
990 | meteorAux = listMeteors[i] | |
|
991 | if meteorAux[-1] == 0: | |
|
954 | 992 | mStart = listMeteors[i][1] |
|
955 | 993 | mPeak = listMeteors[i][2] |
|
956 | 994 | mLag = mPeak - mStart + lag |
@@ -1007,7 +1045,7 class ParametersProc(ProcessingUnit): | |||
|
1007 | 1045 | |
|
1008 | 1046 | #New arrays |
|
1009 | 1047 | arrayMeteors = numpy.array(listMeteors) |
|
1010 |
arrayParameters = numpy.zeros((len(listMeteors), |
|
|
1048 | arrayParameters = numpy.zeros((len(listMeteors), 14)) | |
|
1011 | 1049 | |
|
1012 | 1050 | #Date inclusion |
|
1013 | 1051 | date = re.findall(r'\((.*?)\)', date) |
@@ -1017,14 +1055,18 class ParametersProc(ProcessingUnit): | |||
|
1017 | 1055 | arrayDate = numpy.tile(date, (len(listMeteors), 1)) |
|
1018 | 1056 | |
|
1019 | 1057 | #Meteor array |
|
1020 | arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)] | |
|
1021 | arrayMeteors = numpy.hstack((arrayDate, arrayMeteors)) | |
|
1058 | # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)] | |
|
1059 | # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors)) | |
|
1022 | 1060 | |
|
1023 | 1061 | #Parameters Array |
|
1024 |
arrayParameters[:, |
|
|
1025 |
arrayParameters[:, |
|
|
1062 | arrayParameters[:,:2] = arrayDate #Date | |
|
1063 | arrayParameters[:,2] = heiRang[arrayMeteors[:,0].astype(int)] #Range | |
|
1064 | arrayParameters[:,7:9] = arrayMeteors[:,-3:-1] #Radial velocity and its error | |
|
1065 | arrayParameters[:,9:13] = arrayMeteors[:,7:11] #Phases | |
|
1066 | arrayParameters[:,-1] = arrayMeteors[:,-1] #Error | |
|
1067 | ||
|
1026 | 1068 | |
|
1027 |
return |
|
|
1069 | return arrayParameters | |
|
1028 | 1070 | |
|
1029 | 1071 | def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth): |
|
1030 | 1072 | |
@@ -1246,7 +1288,6 class ParametersProc(ProcessingUnit): | |||
|
1246 | 1288 | self.dataOut.data_param[i,:,h] = minp |
|
1247 | 1289 | return |
|
1248 | 1290 | |
|
1249 | ||
|
1250 | 1291 | def __residFunction(self, p, dp, LT, constants): |
|
1251 | 1292 | |
|
1252 | 1293 | fm = self.dataOut.library.modelFunction(p, constants) |
@@ -1327,7 +1368,6 class WindProfiler(Operation): | |||
|
1327 | 1368 | return A1 |
|
1328 | 1369 | |
|
1329 | 1370 | def __correctValues(self, heiRang, phi, velRadial, SNR): |
|
1330 | ||
|
1331 | 1371 | listPhi = phi.tolist() |
|
1332 | 1372 | maxid = listPhi.index(max(listPhi)) |
|
1333 | 1373 | minid = listPhi.index(min(listPhi)) |
@@ -1540,13 +1580,14 class WindProfiler(Operation): | |||
|
1540 | 1580 | |
|
1541 | 1581 | Parameters affected: Winds |
|
1542 | 1582 |
''' |
|
1583 | print arrayMeteor.shape | |
|
1543 | 1584 | #Settings |
|
1544 | 1585 | nInt = (heightMax - heightMin)/2 |
|
1545 | 1586 | winds = numpy.zeros((2,nInt))*numpy.nan |
|
1546 | 1587 | |
|
1547 | 1588 | #Filter errors |
|
1548 | error = numpy.where(arrayMeteor[:,-1] == 0)[0] | |
|
1549 | finalMeteor = arrayMeteor[error,:] | |
|
1589 | error = numpy.where(arrayMeteor[0,:,-1] == 0)[0] | |
|
1590 | finalMeteor = arrayMeteor[0,error,:] | |
|
1550 | 1591 | |
|
1551 | 1592 | #Meteor Histogram |
|
1552 | 1593 | finalHeights = finalMeteor[:,3] |
@@ -1592,10 +1633,10 class WindProfiler(Operation): | |||
|
1592 | 1633 | def run(self, dataOut, technique, **kwargs): |
|
1593 | 1634 | |
|
1594 | 1635 | param = dataOut.data_param |
|
1595 | if dataOut.abscissaList != None: | |
|
1596 | absc = dataOut.abscissaList[:-1] | |
|
1636 | # if dataOut.abscissaList != None: | |
|
1637 | # absc = dataOut.abscissaList[:-1] | |
|
1597 | 1638 | noise = dataOut.noise |
|
1598 |
heightList = dataOut. |
|
|
1639 | heightList = dataOut.heightList | |
|
1599 | 1640 | SNR = dataOut.data_SNR |
|
1600 | 1641 | |
|
1601 | 1642 | if technique == 'DBS': |
@@ -1616,19 +1657,14 class WindProfiler(Operation): | |||
|
1616 | 1657 | else: correctFactor = 1 |
|
1617 | 1658 | if kwargs.has_key('channelList'): |
|
1618 | 1659 | channelList = kwargs['channelList'] |
|
1660 | if len(channelList) == 2: | |
|
1661 | horizontalOnly = True | |
|
1619 | 1662 | arrayChannel = numpy.array(channelList) |
|
1620 | 1663 | param = param[arrayChannel,:,:] |
|
1621 | 1664 | theta_x = theta_x[arrayChannel] |
|
1622 | 1665 | theta_y = theta_y[arrayChannel] |
|
1623 | 1666 | |
|
1624 | 1667 | velRadial0 = param[:,1,:] #Radial velocity |
|
1625 | ||
|
1626 | if velRadial0.shape[0] != theta_x.shape[0] or velRadial0.shape[0] != theta_y.shape[0]: | |
|
1627 | raise ValueError, "The max number of channels is %d, and the length of cosine director is %d. Please check: dirCosX, dirCosY, elevation or azimuth arguments" %(velRadial0.shape[0], theta_x.shape[0]) | |
|
1628 | ||
|
1629 | if theta_x.shape[0] == 2: | |
|
1630 | horizontalOnly = True | |
|
1631 | ||
|
1632 | 1668 | dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function |
|
1633 | 1669 | dataOut.utctimeInit = dataOut.utctime |
|
1634 | 1670 | dataOut.outputInterval = dataOut.timeInterval |
@@ -1685,7 +1721,7 class WindProfiler(Operation): | |||
|
1685 | 1721 | if self.__isConfig == False: |
|
1686 | 1722 | # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) |
|
1687 | 1723 | #Get Initial LTC time |
|
1688 |
self.__initime = datetime.datetime.utcfromtimestamp( |
|
|
1724 | self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime) | |
|
1689 | 1725 | self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() |
|
1690 | 1726 | |
|
1691 | 1727 | self.__isConfig = True |
@@ -1695,7 +1731,7 class WindProfiler(Operation): | |||
|
1695 | 1731 | self.__firstdata = copy.copy(dataOut) |
|
1696 | 1732 | |
|
1697 | 1733 | else: |
|
1698 |
self.__buffer = numpy. |
|
|
1734 | self.__buffer = numpy.hstack((self.__buffer, dataOut.data_param)) | |
|
1699 | 1735 | |
|
1700 | 1736 | self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready |
|
1701 | 1737 | |
@@ -1779,8 +1815,330 class EWDriftsEstimation(Operation): | |||
|
1779 | 1815 | dataOut.outputInterval = dataOut.timeInterval |
|
1780 | 1816 | return |
|
1781 | 1817 | |
|
1818 | class PhaseCalibration(Operation): | |
|
1819 | ||
|
1820 | __buffer = None | |
|
1821 | ||
|
1822 | __initime = None | |
|
1823 | ||
|
1824 | __dataReady = False | |
|
1825 | ||
|
1826 | __isConfig = False | |
|
1827 | ||
|
1828 | def __checkTime(self, currentTime, initTime, paramInterval, outputInterval): | |
|
1829 | ||
|
1830 | dataTime = currentTime + paramInterval | |
|
1831 | deltaTime = dataTime - initTime | |
|
1832 | ||
|
1833 | if deltaTime >= outputInterval or deltaTime < 0: | |
|
1834 | return True | |
|
1835 | ||
|
1836 | return False | |
|
1837 | ||
|
1838 | def __getGammas(self, pairs, k, d, phases): | |
|
1839 | gammas = numpy.zeros(2) | |
|
1840 | ||
|
1841 | for i in range(len(pairs)): | |
|
1842 | ||
|
1843 | pairi = pairs[i] | |
|
1844 | ||
|
1845 | #Calculating gamma | |
|
1846 | jdcos = phases[:,pairi[1]]/(k*d[pairi[1]]) | |
|
1847 | jgamma = numpy.angle(numpy.exp(1j*(k*d[pairi[0]]*jdcos - phases[:,pairi[0]]))) | |
|
1848 | ||
|
1849 | #Revised distribution | |
|
1850 | jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi)) | |
|
1851 | ||
|
1852 | #Histogram | |
|
1853 | nBins = 64.0 | |
|
1854 | rmin = -0.5*numpy.pi | |
|
1855 | rmax = 0.5*numpy.pi | |
|
1856 | phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax)) | |
|
1857 | ||
|
1858 | meteorsY = phaseHisto[0] | |
|
1859 | phasesX = phaseHisto[1][:-1] | |
|
1860 | width = phasesX[1] - phasesX[0] | |
|
1861 | phasesX += width/2 | |
|
1862 | ||
|
1863 | #Gaussian aproximation | |
|
1864 | bpeak = meteorsY.argmax() | |
|
1865 | peak = meteorsY.max() | |
|
1866 | jmin = bpeak - 5 | |
|
1867 | jmax = bpeak + 5 + 1 | |
|
1868 | ||
|
1869 | if jmin<0: | |
|
1870 | jmin = 0 | |
|
1871 | jmax = 6 | |
|
1872 | elif jmax > meteorsY.size: | |
|
1873 | jmin = meteorsY.size - 6 | |
|
1874 | jmax = meteorsY.size | |
|
1875 | ||
|
1876 | x0 = numpy.array([peak,bpeak,50]) | |
|
1877 | coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax])) | |
|
1878 | ||
|
1879 | #Gammas | |
|
1880 | gammas[i] = coeff[0][1] | |
|
1881 | # gammas[i] = bpeak | |
|
1882 | ||
|
1883 | return gammas | |
|
1884 | ||
|
1885 | def __residualFunction(self, coeffs, y, t): | |
|
1886 | ||
|
1887 | return y - self.__gauss_function(t, coeffs) | |
|
1888 | ||
|
1889 | def __gauss_function(self, t, coeffs): | |
|
1890 | ||
|
1891 | return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2) | |
|
1892 | ||
|
1893 | def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray): | |
|
1894 | meteorOps = MeteorOperations() | |
|
1895 | nchan = 4 | |
|
1896 | pairx = pairsList[0] | |
|
1897 | pairy = pairsList[1] | |
|
1898 | center_xangle = 0 | |
|
1899 | center_yangle = 0 | |
|
1900 | range_angle = numpy.array([8*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4]) | |
|
1901 | ntimes = len(range_angle) | |
|
1902 | ||
|
1903 | nstepsx = 20.0 | |
|
1904 | nstepsy = 20.0 | |
|
1905 | ||
|
1906 | for iz in range(ntimes): | |
|
1907 | min_xangle = -range_angle[iz]/2 + center_xangle | |
|
1908 | max_xangle = range_angle[iz]/2 + center_xangle | |
|
1909 | min_yangle = -range_angle[iz]/2 + center_yangle | |
|
1910 | max_yangle = range_angle[iz]/2 + center_yangle | |
|
1911 | ||
|
1912 | inc_x = (max_xangle-min_xangle)/nstepsx | |
|
1913 | inc_y = (max_yangle-min_yangle)/nstepsy | |
|
1914 | ||
|
1915 | alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle | |
|
1916 | alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle | |
|
1917 | penalty = numpy.zeros((nstepsx,nstepsy)) | |
|
1918 | jph_array = numpy.zeros((nchan,nstepsx,nstepsy)) | |
|
1919 | jph = numpy.zeros(nchan) | |
|
1920 | ||
|
1921 | # Iterations looking for the offset | |
|
1922 | for iy in range(int(nstepsy)): | |
|
1923 | for ix in range(int(nstepsx)): | |
|
1924 | jph[pairy[1]] = alpha_y[iy] | |
|
1925 | jph[pairy[0]] = -gammas[1] + alpha_y[iy]*d[pairy[0]]/d[pairy[1]] | |
|
1926 | ||
|
1927 | jph[pairx[1]] = alpha_x[ix] | |
|
1928 | jph[pairx[0]] = -gammas[0] + alpha_x[ix]*d[pairx[0]]/d[pairx[1]] | |
|
1929 | ||
|
1930 | jph_array[:,ix,iy] = jph | |
|
1931 | ||
|
1932 | meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, jph) | |
|
1933 | error = meteorsArray1[:,-1] | |
|
1934 | ind1 = numpy.where(error==0)[0] | |
|
1935 | penalty[ix,iy] = ind1.size | |
|
1936 | ||
|
1937 | i,j = numpy.unravel_index(penalty.argmax(), penalty.shape) | |
|
1938 | phOffset = jph_array[:,i,j] | |
|
1939 | ||
|
1940 | center_xangle = phOffset[pairx[1]] | |
|
1941 | center_yangle = phOffset[pairy[1]] | |
|
1942 | ||
|
1943 | phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j])) | |
|
1944 | phOffset = phOffset*180/numpy.pi | |
|
1945 | return phOffset | |
|
1946 | ||
|
1947 | ||
|
1948 | def run(self, dataOut, pairs, distances, hmin, hmax, nHours = None): | |
|
1949 | ||
|
1950 | dataOut.flagNoData = True | |
|
1951 | self.__dataReady = False | |
|
1952 | ||
|
1953 | if nHours == None: | |
|
1954 | nHours = 1 | |
|
1955 | ||
|
1956 | dataOut.outputInterval = nHours*3600 | |
|
1957 | ||
|
1958 | if self.__isConfig == False: | |
|
1959 | # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) | |
|
1960 | #Get Initial LTC time | |
|
1961 | self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime) | |
|
1962 | self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() | |
|
1963 | ||
|
1964 | self.__isConfig = True | |
|
1965 | ||
|
1966 | if self.__buffer == None: | |
|
1967 | self.__buffer = dataOut.data_param.copy() | |
|
1968 | ||
|
1969 | else: | |
|
1970 | self.__buffer = numpy.hstack((self.__buffer, dataOut.data_param)) | |
|
1971 | ||
|
1972 | self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready | |
|
1973 | ||
|
1974 | if self.__dataReady: | |
|
1975 | dataOut.utctimeInit = self.__initime | |
|
1976 | self.__initime += dataOut.outputInterval #to erase time offset | |
|
1977 | ||
|
1978 | freq = dataOut.frequency | |
|
1979 | c = dataOut.C #m/s | |
|
1980 | lamb = c/freq | |
|
1981 | k = 2*numpy.pi/lamb | |
|
1982 | azimuth = 0 | |
|
1983 | h = (hmin, hmax) | |
|
1984 | pairsList = ((0,3),(1,2)) | |
|
1985 | ||
|
1986 | meteorsArray = self.__buffer[0,:,:] | |
|
1987 | error = meteorsArray[:,-1] | |
|
1988 | boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14) | |
|
1989 | ind1 = numpy.where(boolError)[0] | |
|
1990 | meteorsArray = meteorsArray[ind1,:] | |
|
1991 | meteorsArray[:,-1] = 0 | |
|
1992 | phases = meteorsArray[:,9:13] | |
|
1993 | ||
|
1994 | #Calculate Gammas | |
|
1995 | gammas = self.__getGammas(pairs, k, distances, phases) | |
|
1996 | # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180 | |
|
1997 | #Calculate Phases | |
|
1998 | phasesOff = self.__getPhases(azimuth, h, pairsList, distances, gammas, meteorsArray) | |
|
1999 | phasesOff = phasesOff.reshape((1,phasesOff.size)) | |
|
2000 | dataOut.data_output = -phasesOff | |
|
2001 | dataOut.flagNoData = False | |
|
2002 | self.__buffer = None | |
|
2003 | ||
|
2004 | ||
|
2005 | return | |
|
2006 | ||
|
2007 | class MeteorOperations(): | |
|
2008 | ||
|
2009 | def __init__(self): | |
|
2010 | ||
|
2011 | return | |
|
2012 | ||
|
2013 | def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, jph): | |
|
1782 | 2014 | |
|
2015 | arrayParameters = arrayParameters0.copy() | |
|
2016 | hmin = h[0] | |
|
2017 | hmax = h[1] | |
|
1783 | 2018 | |
|
2019 | #Calculate AOA (Error N 3, 4) | |
|
2020 | #JONES ET AL. 1998 | |
|
2021 | AOAthresh = numpy.pi/8 | |
|
2022 | error = arrayParameters[:,-1] | |
|
2023 | phases = -arrayParameters[:,9:13] + jph | |
|
2024 | arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth) | |
|
2025 | ||
|
2026 | #Calculate Heights (Error N 13 and 14) | |
|
2027 | error = arrayParameters[:,-1] | |
|
2028 | Ranges = arrayParameters[:,2] | |
|
2029 | zenith = arrayParameters[:,5] | |
|
2030 | arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax) | |
|
2031 | ||
|
2032 | #----------------------- Get Final data ------------------------------------ | |
|
2033 | # error = arrayParameters[:,-1] | |
|
2034 | # ind1 = numpy.where(error==0)[0] | |
|
2035 | # arrayParameters = arrayParameters[ind1,:] | |
|
2036 | ||
|
2037 | return arrayParameters | |
|
2038 | ||
|
2039 | def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth): | |
|
2040 | ||
|
2041 | arrayAOA = numpy.zeros((phases.shape[0],3)) | |
|
2042 | cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList) | |
|
2043 | ||
|
2044 | arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth) | |
|
2045 | cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1) | |
|
2046 | arrayAOA[:,2] = cosDirError | |
|
2047 | ||
|
2048 | azimuthAngle = arrayAOA[:,0] | |
|
2049 | zenithAngle = arrayAOA[:,1] | |
|
2050 | ||
|
2051 | #Setting Error | |
|
2052 | #Number 3: AOA not fesible | |
|
2053 | indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0] | |
|
2054 | error[indInvalid] = 3 | |
|
2055 | #Number 4: Large difference in AOAs obtained from different antenna baselines | |
|
2056 | indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0] | |
|
2057 | error[indInvalid] = 4 | |
|
2058 | return arrayAOA, error | |
|
2059 | ||
|
2060 | def __getDirectionCosines(self, arrayPhase, pairsList): | |
|
2061 | ||
|
2062 | #Initializing some variables | |
|
2063 | ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi | |
|
2064 | ang_aux = ang_aux.reshape(1,ang_aux.size) | |
|
2065 | ||
|
2066 | cosdir = numpy.zeros((arrayPhase.shape[0],2)) | |
|
2067 | cosdir0 = numpy.zeros((arrayPhase.shape[0],2)) | |
|
2068 | ||
|
2069 | ||
|
2070 | for i in range(2): | |
|
2071 | #First Estimation | |
|
2072 | phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]] | |
|
2073 | #Dealias | |
|
2074 | indcsi = numpy.where(phi0_aux > numpy.pi) | |
|
2075 | phi0_aux[indcsi] -= 2*numpy.pi | |
|
2076 | indcsi = numpy.where(phi0_aux < -numpy.pi) | |
|
2077 | phi0_aux[indcsi] += 2*numpy.pi | |
|
2078 | #Direction Cosine 0 | |
|
2079 | cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5) | |
|
2080 | ||
|
2081 | #Most-Accurate Second Estimation | |
|
2082 | phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]] | |
|
2083 | phi1_aux = phi1_aux.reshape(phi1_aux.size,1) | |
|
2084 | #Direction Cosine 1 | |
|
2085 | cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5) | |
|
2086 | ||
|
2087 | #Searching the correct Direction Cosine | |
|
2088 | cosdir0_aux = cosdir0[:,i] | |
|
2089 | cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1) | |
|
2090 | #Minimum Distance | |
|
2091 | cosDiff = (cosdir1 - cosdir0_aux)**2 | |
|
2092 | indcos = cosDiff.argmin(axis = 1) | |
|
2093 | #Saving Value obtained | |
|
2094 | cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos] | |
|
2095 | ||
|
2096 | return cosdir0, cosdir | |
|
2097 | ||
|
2098 | def __calculateAOA(self, cosdir, azimuth): | |
|
2099 | cosdirX = cosdir[:,0] | |
|
2100 | cosdirY = cosdir[:,1] | |
|
2101 | ||
|
2102 | zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi | |
|
2103 | azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east | |
|
2104 | angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose() | |
|
2105 | ||
|
2106 | return angles | |
|
2107 | ||
|
2108 | def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight): | |
|
2109 | ||
|
2110 | Ramb = 375 #Ramb = c/(2*PRF) | |
|
2111 | Re = 6371 #Earth Radius | |
|
2112 | heights = numpy.zeros(Ranges.shape) | |
|
2113 | ||
|
2114 | R_aux = numpy.array([0,1,2])*Ramb | |
|
2115 | R_aux = R_aux.reshape(1,R_aux.size) | |
|
2116 | ||
|
2117 | Ranges = Ranges.reshape(Ranges.size,1) | |
|
2118 | ||
|
2119 | Ri = Ranges + R_aux | |
|
2120 | hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re | |
|
2121 | ||
|
2122 | #Check if there is a height between 70 and 110 km | |
|
2123 | h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1) | |
|
2124 | ind_h = numpy.where(h_bool == 1)[0] | |
|
2125 | ||
|
2126 | hCorr = hi[ind_h, :] | |
|
2127 | ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight)) | |
|
2128 | ||
|
2129 | hCorr = hi[ind_hCorr] | |
|
2130 | heights[ind_h] = hCorr | |
|
2131 | ||
|
2132 | #Setting Error | |
|
2133 | #Number 13: Height unresolvable echo: not valid height within 70 to 110 km | |
|
2134 | #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km | |
|
2135 | ||
|
2136 | indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0] | |
|
2137 | error[indInvalid2] = 14 | |
|
2138 | indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0] | |
|
2139 | error[indInvalid1] = 13 | |
|
2140 | ||
|
2141 | return heights, error | |
|
1784 | 2142 | |
|
1785 | 2143 | |
|
1786 | 2144 | No newline at end of file |
General Comments 0
You need to be logged in to leave comments.
Login now