##// END OF EJS Templates
data...
Julio Valdez -
r608:3cc55c179318
parent child
Show More
@@ -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()
@@ -1,1115 +1,1124
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10
11 11 from jroheaderIO import SystemHeader, RadarControllerHeader
12 12
13 13 def getNumpyDtype(dataTypeCode):
14 14
15 15 if dataTypeCode == 0:
16 16 numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
17 17 elif dataTypeCode == 1:
18 18 numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
19 19 elif dataTypeCode == 2:
20 20 numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
21 21 elif dataTypeCode == 3:
22 22 numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
23 23 elif dataTypeCode == 4:
24 24 numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
25 25 elif dataTypeCode == 5:
26 26 numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
27 27 else:
28 28 raise ValueError, 'dataTypeCode was not defined'
29 29
30 30 return numpyDtype
31 31
32 32 def getDataTypeCode(numpyDtype):
33 33
34 34 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
35 35 datatype = 0
36 36 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
37 37 datatype = 1
38 38 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
39 39 datatype = 2
40 40 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
41 41 datatype = 3
42 42 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
43 43 datatype = 4
44 44 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
45 45 datatype = 5
46 46 else:
47 47 datatype = None
48 48
49 49 return datatype
50 50
51 51 def hildebrand_sekhon(data, navg):
52 52 """
53 53 This method is for the objective determination of the noise level in Doppler spectra. This
54 54 implementation technique is based on the fact that the standard deviation of the spectral
55 55 densities is equal to the mean spectral density for white Gaussian noise
56 56
57 57 Inputs:
58 58 Data : heights
59 59 navg : numbers of averages
60 60
61 61 Return:
62 62 -1 : any error
63 63 anoise : noise's level
64 64 """
65 65
66 66 sortdata = numpy.sort(data,axis=None)
67 67 lenOfData = len(sortdata)
68 68 nums_min = lenOfData/10
69 69
70 70 if (lenOfData/10) > 2:
71 71 nums_min = lenOfData/10
72 72 else:
73 73 nums_min = 2
74 74
75 75 sump = 0.
76 76
77 77 sumq = 0.
78 78
79 79 j = 0
80 80
81 81 cont = 1
82 82
83 83 while((cont==1)and(j<lenOfData)):
84 84
85 85 sump += sortdata[j]
86 86
87 87 sumq += sortdata[j]**2
88 88
89 89 if j > nums_min:
90 90 rtest = float(j)/(j-1) + 1.0/navg
91 91 if ((sumq*j) > (rtest*sump**2)):
92 92 j = j - 1
93 93 sump = sump - sortdata[j]
94 94 sumq = sumq - sortdata[j]**2
95 95 cont = 0
96 96
97 97 j += 1
98 98
99 99 lnoise = sump /j
100 100 stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
101 101 return lnoise
102 102
103 103 class Beam:
104 104 def __init__(self):
105 105 self.codeList = []
106 106 self.azimuthList = []
107 107 self.zenithList = []
108 108
109 109 class GenericData(object):
110 110
111 111 flagNoData = True
112 112
113 113 def __init__(self):
114 114
115 115 raise ValueError, "This class has not been implemented"
116 116
117 117 def copy(self, inputObj=None):
118 118
119 119 if inputObj == None:
120 120 return copy.deepcopy(self)
121 121
122 122 for key in inputObj.__dict__.keys():
123 123 self.__dict__[key] = inputObj.__dict__[key]
124 124
125 125 def deepcopy(self):
126 126
127 127 return copy.deepcopy(self)
128 128
129 129 def isEmpty(self):
130 130
131 131 return self.flagNoData
132 132
133 133 class JROData(GenericData):
134 134
135 135 # m_BasicHeader = BasicHeader()
136 136 # m_ProcessingHeader = ProcessingHeader()
137 137
138 138 systemHeaderObj = SystemHeader()
139 139
140 140 radarControllerHeaderObj = RadarControllerHeader()
141 141
142 142 # data = None
143 143
144 144 type = None
145 145
146 146 datatype = None #dtype but in string
147 147
148 148 # dtype = None
149 149
150 150 # nChannels = None
151 151
152 152 # nHeights = None
153 153
154 154 nProfiles = None
155 155
156 156 heightList = None
157 157
158 158 channelList = None
159 159
160 160 flagDiscontinuousBlock = False
161 161
162 162 useLocalTime = False
163 163
164 164 utctime = None
165 165
166 166 timeZone = None
167 167
168 168 dstFlag = None
169 169
170 170 errorCount = None
171 171
172 172 blocksize = None
173 173
174 174 # nCode = None
175 175 #
176 176 # nBaud = None
177 177 #
178 178 # code = None
179 179
180 180 flagDecodeData = False #asumo q la data no esta decodificada
181 181
182 182 flagDeflipData = False #asumo q la data no esta sin flip
183 183
184 184 flagShiftFFT = False
185 185
186 186 # ippSeconds = None
187 187
188 188 # timeInterval = None
189 189
190 190 nCohInt = None
191 191
192 192 # noise = None
193 193
194 194 windowOfFilter = 1
195 195
196 196 #Speed of ligth
197 197 C = 3e8
198 198
199 199 frequency = 49.92e6
200 200
201 201 realtime = False
202 202
203 203 beacon_heiIndexList = None
204 204
205 205 last_block = None
206 206
207 207 blocknow = None
208 208
209 209 azimuth = None
210 210
211 211 zenith = None
212 212
213 213 beam = Beam()
214 214
215 215 profileIndex = None
216 216
217 217 def __init__(self):
218 218
219 219 raise ValueError, "This class has not been implemented"
220 220
221 221 def getNoise(self):
222 222
223 223 raise ValueError, "Not implemented"
224 224
225 225 def getNChannels(self):
226 226
227 227 return len(self.channelList)
228 228
229 229 def getChannelIndexList(self):
230 230
231 231 return range(self.nChannels)
232 232
233 233 def getNHeights(self):
234 234
235 235 return len(self.heightList)
236 236
237 237 def getHeiRange(self, extrapoints=0):
238 238
239 239 heis = self.heightList
240 240 # deltah = self.heightList[1] - self.heightList[0]
241 241 #
242 242 # heis.append(self.heightList[-1])
243 243
244 244 return heis
245 245
246 246 def getltctime(self):
247 247
248 248 if self.useLocalTime:
249 249 return self.utctime - self.timeZone*60
250 250
251 251 return self.utctime
252 252
253 253 def getDatatime(self):
254 254
255 255 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
256 256 return datatimeValue
257 257
258 258 def getTimeRange(self):
259 259
260 260 datatime = []
261 261
262 262 datatime.append(self.ltctime)
263 263 datatime.append(self.ltctime + self.timeInterval+60)
264 264
265 265 datatime = numpy.array(datatime)
266 266
267 267 return datatime
268 268
269 269 def getFmax(self):
270 270
271 271 PRF = 1./(self.ippSeconds * self.nCohInt)
272 272
273 273 fmax = PRF/2.
274 274
275 275 return fmax
276 276
277 277 def getVmax(self):
278 278
279 279 _lambda = self.C/self.frequency
280 280
281 281 vmax = self.getFmax() * _lambda
282 282
283 283 return vmax
284 284
285 285 def get_ippSeconds(self):
286 286 '''
287 287 '''
288 288 return self.radarControllerHeaderObj.ippSeconds
289 289
290 290 def set_ippSeconds(self, ippSeconds):
291 291 '''
292 292 '''
293 293
294 294 self.radarControllerHeaderObj.ippSeconds = ippSeconds
295 295
296 296 return
297 297
298 298 def get_dtype(self):
299 299 '''
300 300 '''
301 301 return getNumpyDtype(self.datatype)
302 302
303 303 def set_dtype(self, numpyDtype):
304 304 '''
305 305 '''
306 306
307 307 self.datatype = getDataTypeCode(numpyDtype)
308 308
309 309 def get_code(self):
310 310 '''
311 311 '''
312 312 return self.radarControllerHeaderObj.code
313 313
314 314 def set_code(self, code):
315 315 '''
316 316 '''
317 317 self.radarControllerHeaderObj.code = code
318 318
319 319 return
320 320
321 321 def get_ncode(self):
322 322 '''
323 323 '''
324 324 return self.radarControllerHeaderObj.nCode
325 325
326 326 def set_ncode(self, nCode):
327 327 '''
328 328 '''
329 329 self.radarControllerHeaderObj.nCode = nCode
330 330
331 331 return
332 332
333 333 def get_nbaud(self):
334 334 '''
335 335 '''
336 336 return self.radarControllerHeaderObj.nBaud
337 337
338 338 def set_nbaud(self, nBaud):
339 339 '''
340 340 '''
341 341 self.radarControllerHeaderObj.nBaud = nBaud
342 342
343 343 return
344 344 # def getTimeInterval(self):
345 345 #
346 346 # raise IOError, "This method should be implemented inside each Class"
347 347
348 348 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
349 349 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
350 350 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
351 351 #noise = property(getNoise, "I'm the 'nHeights' property.")
352 352 datatime = property(getDatatime, "I'm the 'datatime' property")
353 353 ltctime = property(getltctime, "I'm the 'ltctime' property")
354 354 ippSeconds = property(get_ippSeconds, set_ippSeconds)
355 355 dtype = property(get_dtype, set_dtype)
356 356 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
357 357 code = property(get_code, set_code)
358 358 nCode = property(get_ncode, set_ncode)
359 359 nBaud = property(get_nbaud, set_nbaud)
360 360
361 361 class Voltage(JROData):
362 362
363 363 #data es un numpy array de 2 dmensiones (canales, alturas)
364 364 data = None
365 365
366 366 def __init__(self):
367 367 '''
368 368 Constructor
369 369 '''
370 370
371 371 self.useLocalTime = True
372 372
373 373 self.radarControllerHeaderObj = RadarControllerHeader()
374 374
375 375 self.systemHeaderObj = SystemHeader()
376 376
377 377 self.type = "Voltage"
378 378
379 379 self.data = None
380 380
381 381 # self.dtype = None
382 382
383 383 # self.nChannels = 0
384 384
385 385 # self.nHeights = 0
386 386
387 387 self.nProfiles = None
388 388
389 389 self.heightList = None
390 390
391 391 self.channelList = None
392 392
393 393 # self.channelIndexList = None
394 394
395 395 self.flagNoData = True
396 396
397 397 self.flagDiscontinuousBlock = False
398 398
399 399 self.utctime = None
400 400
401 401 self.timeZone = None
402 402
403 403 self.dstFlag = None
404 404
405 405 self.errorCount = None
406 406
407 407 self.nCohInt = None
408 408
409 409 self.blocksize = None
410 410
411 411 self.flagDecodeData = False #asumo q la data no esta decodificada
412 412
413 413 self.flagDeflipData = False #asumo q la data no esta sin flip
414 414
415 415 self.flagShiftFFT = False
416 416
417 417 self.flagDataAsBlock = False #Asumo que la data es leida perfil a perfil
418 418
419 419 self.profileIndex = 0
420 420
421 421 def getNoisebyHildebrand(self, channel = None):
422 422 """
423 423 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
424 424
425 425 Return:
426 426 noiselevel
427 427 """
428 428
429 429 if channel != None:
430 430 data = self.data[channel]
431 431 nChannels = 1
432 432 else:
433 433 data = self.data
434 434 nChannels = self.nChannels
435 435
436 436 noise = numpy.zeros(nChannels)
437 437 power = data * numpy.conjugate(data)
438 438
439 439 for thisChannel in range(nChannels):
440 440 if nChannels == 1:
441 441 daux = power[:].real
442 442 else:
443 443 daux = power[thisChannel,:].real
444 444 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
445 445
446 446 return noise
447 447
448 448 def getNoise(self, type = 1, channel = None):
449 449
450 450 if type == 1:
451 451 noise = self.getNoisebyHildebrand(channel)
452 452
453 453 return 10*numpy.log10(noise)
454 454
455 455 def getPower(self, channel = None):
456 456
457 457 if channel != None:
458 458 data = self.data[channel]
459 459 else:
460 460 data = self.data
461 461
462 462 power = data * numpy.conjugate(data)
463 463
464 464 return 10*numpy.log10(power.real)
465 465
466 466 def getTimeInterval(self):
467 467
468 468 timeInterval = self.ippSeconds * self.nCohInt
469 469
470 470 return timeInterval
471 471
472 472 noise = property(getNoise, "I'm the 'nHeights' property.")
473 473 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
474 474
475 475 class Spectra(JROData):
476 476
477 477 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
478 478 data_spc = None
479 479
480 480 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
481 481 data_cspc = None
482 482
483 483 #data es un numpy array de 2 dmensiones (canales, alturas)
484 484 data_dc = None
485 485
486 486 nFFTPoints = None
487 487
488 488 # nPairs = None
489 489
490 490 pairsList = None
491 491
492 492 nIncohInt = None
493 493
494 494 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
495 495
496 496 nCohInt = None #se requiere para determinar el valor de timeInterval
497 497
498 498 ippFactor = None
499 499
500 500 profileIndex = 0
501 501
502 502 def __init__(self):
503 503 '''
504 504 Constructor
505 505 '''
506 506
507 507 self.useLocalTime = True
508 508
509 509 self.radarControllerHeaderObj = RadarControllerHeader()
510 510
511 511 self.systemHeaderObj = SystemHeader()
512 512
513 513 self.type = "Spectra"
514 514
515 515 # self.data = None
516 516
517 517 # self.dtype = None
518 518
519 519 # self.nChannels = 0
520 520
521 521 # self.nHeights = 0
522 522
523 523 self.nProfiles = None
524 524
525 525 self.heightList = None
526 526
527 527 self.channelList = None
528 528
529 529 # self.channelIndexList = None
530 530
531 531 self.pairsList = None
532 532
533 533 self.flagNoData = True
534 534
535 535 self.flagDiscontinuousBlock = False
536 536
537 537 self.utctime = None
538 538
539 539 self.nCohInt = None
540 540
541 541 self.nIncohInt = None
542 542
543 543 self.blocksize = None
544 544
545 545 self.nFFTPoints = None
546 546
547 547 self.wavelength = None
548 548
549 549 self.flagDecodeData = False #asumo q la data no esta decodificada
550 550
551 551 self.flagDeflipData = False #asumo q la data no esta sin flip
552 552
553 553 self.flagShiftFFT = False
554 554
555 555 self.ippFactor = 1
556 556
557 557 #self.noise = None
558 558
559 559 self.beacon_heiIndexList = []
560 560
561 561 self.noise_estimation = None
562 562
563 563
564 564 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
565 565 """
566 566 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
567 567
568 568 Return:
569 569 noiselevel
570 570 """
571 571
572 572 noise = numpy.zeros(self.nChannels)
573 573
574 574 for channel in range(self.nChannels):
575 575 daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
576 576 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
577 577
578 578 return noise
579 579
580 580 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
581 581
582 582 if self.noise_estimation != None:
583 583 return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
584 584 else:
585 585 noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
586 586 return noise
587 587
588 588
589 589 def getFreqRange(self, extrapoints=0):
590 590
591 591 deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
592 592 freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
593 593
594 594 return freqrange
595 595
596 596 def getVelRange(self, extrapoints=0):
597 597
598 598 deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
599 599 velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
600 600
601 601 return velrange
602 602
603 603 def getNPairs(self):
604 604
605 605 return len(self.pairsList)
606 606
607 607 def getPairsIndexList(self):
608 608
609 609 return range(self.nPairs)
610 610
611 611 def getNormFactor(self):
612 612 pwcode = 1
613 613 if self.flagDecodeData:
614 614 pwcode = numpy.sum(self.code[0]**2)
615 615 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
616 616 normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
617 617
618 618 return normFactor
619 619
620 620 def getFlagCspc(self):
621 621
622 622 if self.data_cspc == None:
623 623 return True
624 624
625 625 return False
626 626
627 627 def getFlagDc(self):
628 628
629 629 if self.data_dc == None:
630 630 return True
631 631
632 632 return False
633 633
634 634 def getTimeInterval(self):
635 635
636 636 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
637 637
638 638 return timeInterval
639 639
640 640 nPairs = property(getNPairs, "I'm the 'nPairs' property.")
641 641 pairsIndexList = property(getPairsIndexList, "I'm the 'pairsIndexList' property.")
642 642 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
643 643 flag_cspc = property(getFlagCspc)
644 644 flag_dc = property(getFlagDc)
645 645 noise = property(getNoise, "I'm the 'nHeights' property.")
646 646 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
647 647
648 648 class SpectraHeis(Spectra):
649 649
650 650 data_spc = None
651 651
652 652 data_cspc = None
653 653
654 654 data_dc = None
655 655
656 656 nFFTPoints = None
657 657
658 658 # nPairs = None
659 659
660 660 pairsList = None
661 661
662 662 nCohInt = None
663 663
664 664 nIncohInt = None
665 665
666 666 def __init__(self):
667 667
668 668 self.radarControllerHeaderObj = RadarControllerHeader()
669 669
670 670 self.systemHeaderObj = SystemHeader()
671 671
672 672 self.type = "SpectraHeis"
673 673
674 674 # self.dtype = None
675 675
676 676 # self.nChannels = 0
677 677
678 678 # self.nHeights = 0
679 679
680 680 self.nProfiles = None
681 681
682 682 self.heightList = None
683 683
684 684 self.channelList = None
685 685
686 686 # self.channelIndexList = None
687 687
688 688 self.flagNoData = True
689 689
690 690 self.flagDiscontinuousBlock = False
691 691
692 692 # self.nPairs = 0
693 693
694 694 self.utctime = None
695 695
696 696 self.blocksize = None
697 697
698 698 self.profileIndex = 0
699 699
700 700 self.nCohInt = 1
701 701
702 702 self.nIncohInt = 1
703 703
704 704 def getNormFactor(self):
705 705 pwcode = 1
706 706 if self.flagDecodeData:
707 707 pwcode = numpy.sum(self.code[0]**2)
708 708
709 709 normFactor = self.nIncohInt*self.nCohInt*pwcode
710 710
711 711 return normFactor
712 712
713 713 def getTimeInterval(self):
714 714
715 715 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
716 716
717 717 return timeInterval
718 718
719 719 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
720 720 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
721 721
722 722 class Fits(JROData):
723 723
724 724 heightList = None
725 725
726 726 channelList = None
727 727
728 728 flagNoData = True
729 729
730 730 flagDiscontinuousBlock = False
731 731
732 732 useLocalTime = False
733 733
734 734 utctime = None
735 735
736 736 timeZone = None
737 737
738 738 # ippSeconds = None
739 739
740 740 # timeInterval = None
741 741
742 742 nCohInt = None
743 743
744 744 nIncohInt = None
745 745
746 746 noise = None
747 747
748 748 windowOfFilter = 1
749 749
750 750 #Speed of ligth
751 751 C = 3e8
752 752
753 753 frequency = 49.92e6
754 754
755 755 realtime = False
756 756
757 757
758 758 def __init__(self):
759 759
760 760 self.type = "Fits"
761 761
762 762 self.nProfiles = None
763 763
764 764 self.heightList = None
765 765
766 766 self.channelList = None
767 767
768 768 # self.channelIndexList = None
769 769
770 770 self.flagNoData = True
771 771
772 772 self.utctime = None
773 773
774 774 self.nCohInt = 1
775 775
776 776 self.nIncohInt = 1
777 777
778 778 self.useLocalTime = True
779 779
780 780 self.profileIndex = 0
781 781
782 782 # self.utctime = None
783 783 # self.timeZone = None
784 784 # self.ltctime = None
785 785 # self.timeInterval = None
786 786 # self.header = None
787 787 # self.data_header = None
788 788 # self.data = None
789 789 # self.datatime = None
790 790 # self.flagNoData = False
791 791 # self.expName = ''
792 792 # self.nChannels = None
793 793 # self.nSamples = None
794 794 # self.dataBlocksPerFile = None
795 795 # self.comments = ''
796 796 #
797 797
798 798
799 799 def getltctime(self):
800 800
801 801 if self.useLocalTime:
802 802 return self.utctime - self.timeZone*60
803 803
804 804 return self.utctime
805 805
806 806 def getDatatime(self):
807 807
808 808 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
809 809 return datatime
810 810
811 811 def getTimeRange(self):
812 812
813 813 datatime = []
814 814
815 815 datatime.append(self.ltctime)
816 816 datatime.append(self.ltctime + self.timeInterval)
817 817
818 818 datatime = numpy.array(datatime)
819 819
820 820 return datatime
821 821
822 822 def getHeiRange(self):
823 823
824 824 heis = self.heightList
825 825
826 826 return heis
827 827
828 828 def isEmpty(self):
829 829
830 830 return self.flagNoData
831 831
832 832 def getNHeights(self):
833 833
834 834 return len(self.heightList)
835 835
836 836 def getNChannels(self):
837 837
838 838 return len(self.channelList)
839 839
840 840 def getChannelIndexList(self):
841 841
842 842 return range(self.nChannels)
843 843
844 844 def getNoise(self, type = 1):
845 845
846 846 #noise = numpy.zeros(self.nChannels)
847 847
848 848 if type == 1:
849 849 noise = self.getNoisebyHildebrand()
850 850
851 851 if type == 2:
852 852 noise = self.getNoisebySort()
853 853
854 854 if type == 3:
855 855 noise = self.getNoisebyWindow()
856 856
857 857 return noise
858 858
859 859 def getTimeInterval(self):
860 860
861 861 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
862 862
863 863 return timeInterval
864 864
865 865 datatime = property(getDatatime, "I'm the 'datatime' property")
866 866 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
867 867 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
868 868 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
869 869 noise = property(getNoise, "I'm the 'nHeights' property.")
870 870 datatime = property(getDatatime, "I'm the 'datatime' property")
871 871 ltctime = property(getltctime, "I'm the 'ltctime' property")
872 872 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
873 873
874 874 class Correlation(JROData):
875 875
876 876 noise = None
877 877
878 878 SNR = None
879 879
880 880 pairsAutoCorr = None #Pairs of Autocorrelation
881 881
882 882 #--------------------------------------------------
883 883
884 884 data_corr = None
885 885
886 886 data_volt = None
887 887
888 888 lagT = None # each element value is a profileIndex
889 889
890 890 lagR = None # each element value is in km
891 891
892 892 pairsList = None
893 893
894 894 calculateVelocity = None
895 895
896 896 nPoints = None
897 897
898 898 nAvg = None
899 899
900 900 bufferSize = None
901 901
902 902 def __init__(self):
903 903 '''
904 904 Constructor
905 905 '''
906 906 self.radarControllerHeaderObj = RadarControllerHeader()
907 907
908 908 self.systemHeaderObj = SystemHeader()
909 909
910 910 self.type = "Correlation"
911 911
912 912 self.data = None
913 913
914 914 self.dtype = None
915 915
916 916 self.nProfiles = None
917 917
918 918 self.heightList = None
919 919
920 920 self.channelList = None
921 921
922 922 self.flagNoData = True
923 923
924 924 self.flagDiscontinuousBlock = False
925 925
926 926 self.utctime = None
927 927
928 928 self.timeZone = None
929 929
930 930 self.dstFlag = None
931 931
932 932 self.errorCount = None
933 933
934 934 self.blocksize = None
935 935
936 936 self.flagDecodeData = False #asumo q la data no esta decodificada
937 937
938 938 self.flagDeflipData = False #asumo q la data no esta sin flip
939 939
940 940 self.pairsList = None
941 941
942 942 self.nPoints = None
943 943
944 944 def getLagTRange(self, extrapoints=0):
945 945
946 946 lagTRange = self.lagT
947 947 diff = lagTRange[1] - lagTRange[0]
948 948 extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1]
949 949 lagTRange = numpy.hstack((lagTRange, extra))
950 950
951 951 return lagTRange
952 952
953 953 def getLagRRange(self, extrapoints=0):
954 954
955 955 return self.lagR
956 956
957 957 def getPairsList(self):
958 958
959 959 return self.pairsList
960 960
961 961 def getCalculateVelocity(self):
962 962
963 963 return self.calculateVelocity
964 964
965 965 def getNPoints(self):
966 966
967 967 return self.nPoints
968 968
969 969 def getNAvg(self):
970 970
971 971 return self.nAvg
972 972
973 973 def getBufferSize(self):
974 974
975 975 return self.bufferSize
976 976
977 977 def getPairsAutoCorr(self):
978 978 pairsList = self.pairsList
979 979 pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan
980 980
981 981 for l in range(len(pairsList)):
982 982 firstChannel = pairsList[l][0]
983 983 secondChannel = pairsList[l][1]
984 984
985 985 #Obteniendo pares de Autocorrelacion
986 986 if firstChannel == secondChannel:
987 987 pairsAutoCorr[firstChannel] = int(l)
988 988
989 989 pairsAutoCorr = pairsAutoCorr.astype(int)
990 990
991 991 return pairsAutoCorr
992 992
993 993 def getNoise(self, mode = 2):
994 994
995 995 indR = numpy.where(self.lagR == 0)[0][0]
996 996 indT = numpy.where(self.lagT == 0)[0][0]
997 997
998 998 jspectra0 = self.data_corr[:,:,indR,:]
999 999 jspectra = copy.copy(jspectra0)
1000 1000
1001 1001 num_chan = jspectra.shape[0]
1002 1002 num_hei = jspectra.shape[2]
1003 1003
1004 1004 freq_dc = jspectra.shape[1]/2
1005 1005 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
1006 1006
1007 1007 if ind_vel[0]<0:
1008 1008 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
1009 1009
1010 1010 if mode == 1:
1011 1011 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
1012 1012
1013 1013 if mode == 2:
1014 1014
1015 1015 vel = numpy.array([-2,-1,1,2])
1016 1016 xx = numpy.zeros([4,4])
1017 1017
1018 1018 for fil in range(4):
1019 1019 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
1020 1020
1021 1021 xx_inv = numpy.linalg.inv(xx)
1022 1022 xx_aux = xx_inv[0,:]
1023 1023
1024 1024 for ich in range(num_chan):
1025 1025 yy = jspectra[ich,ind_vel,:]
1026 1026 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
1027 1027
1028 1028 junkid = jspectra[ich,freq_dc,:]<=0
1029 1029 cjunkid = sum(junkid)
1030 1030
1031 1031 if cjunkid.any():
1032 1032 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
1033 1033
1034 1034 noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
1035 1035
1036 1036 return noise
1037 1037
1038 1038 def getTimeInterval(self):
1039 1039
1040 1040 timeInterval = self.ippSeconds * self.nCohInt * self.nPoints
1041 1041
1042 1042 return timeInterval
1043 1043
1044 1044 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1045 1045 # pairsList = property(getPairsList, "I'm the 'pairsList' property.")
1046 1046 # nPoints = property(getNPoints, "I'm the 'nPoints' property.")
1047 1047 calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.")
1048 1048 nAvg = property(getNAvg, "I'm the 'nAvg' property.")
1049 1049 bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.")
1050 1050
1051 1051
1052 1052 class Parameters(JROData):
1053 1053
1054 1054 #Information from previous data
1055 1055
1056 1056 inputUnit = None #Type of data to be processed
1057 1057
1058 1058 operation = None #Type of operation to parametrize
1059 1059
1060 1060 normFactor = None #Normalization Factor
1061 1061
1062 1062 groupList = None #List of Pairs, Groups, etc
1063 1063
1064 1064 #Parameters
1065 1065
1066 1066 data_param = None #Parameters obtained
1067 1067
1068 1068 data_pre = None #Data Pre Parametrization
1069 1069
1070 1070 data_SNR = None #Signal to Noise Ratio
1071 1071
1072 heightRange = None #Heights
1072 # heightRange = None #Heights
1073 1073
1074 abscissaRange = None #Abscissa, can be velocities, lags or time
1074 abscissaList = None #Abscissa, can be velocities, lags or time
1075 1075
1076 1076 noise = None #Noise Potency
1077 1077
1078 # initUtcTime = None #Initial UTC time
1078 utctimeInit = None #Initial UTC time
1079 1079
1080 1080 paramInterval = None #Time interval to calculate Parameters in seconds
1081 1081
1082 1082 #Fitting
1083 1083
1084 1084 data_error = None #Error of the estimation
1085 1085
1086 1086 constants = None
1087 1087
1088 1088 library = None
1089 1089
1090 1090 #Output signal
1091 1091
1092 1092 outputInterval = None #Time interval to calculate output signal in seconds
1093 1093
1094 1094 data_output = None #Out signal
1095 1095
1096
1097
1096 1098 def __init__(self):
1097 1099 '''
1098 1100 Constructor
1099 1101 '''
1100 1102 self.radarControllerHeaderObj = RadarControllerHeader()
1101 1103
1102 1104 self.systemHeaderObj = SystemHeader()
1103 1105
1104 1106 self.type = "Parameters"
1105
1107
1106 1108 def getTimeRange1(self):
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
1115 return datatime
1124 return datatime
@@ -1,610 +1,610
1 1 import os
2 2 import numpy
3 3 import time, datetime
4 4 import mpldriver
5 5
6 6 from schainpy.model.proc.jroproc_base import Operation
7 7
8 8 def isRealtime(utcdatatime):
9 9 utcnow = time.mktime(time.localtime())
10 10 delta = abs(utcnow - utcdatatime) # abs
11 11 if delta >= 30.:
12 12 return False
13 13 return True
14 14
15 15 class Figure(Operation):
16 16
17 17 __driver = mpldriver
18 18 __isConfigThread = False
19 19 fig = None
20 20
21 21 id = None
22 22 wintitle = None
23 23 width = None
24 24 height = None
25 25 nplots = None
26 26 timerange = None
27 27
28 28 axesObjList = []
29 29
30 30 WIDTH = None
31 31 HEIGHT = None
32 32 PREFIX = 'fig'
33 33
34 34 xmin = None
35 35 xmax = None
36 36
37 37 counter_imagwr = 0
38 38
39 39 figfile = None
40 40
41 41 def __init__(self):
42 42
43 43 raise ValueError, "This method is not implemented"
44 44
45 45 def __del__(self):
46 46
47 47 self.__driver.closeFigure()
48 48
49 49 def getFilename(self, name, ext='.png'):
50 50
51 51 path = '%s%03d' %(self.PREFIX, self.id)
52 52 filename = '%s_%s%s' %(self.PREFIX, name, ext)
53 53 return os.path.join(path, filename)
54 54
55 55 def getAxesObjList(self):
56 56
57 57 return self.axesObjList
58 58
59 59 def getSubplots(self):
60 60
61 61 raise ValueError, "Abstract method: This method should be defined"
62 62
63 63 def getScreenDim(self, widthplot, heightplot):
64 64
65 65 nrow, ncol = self.getSubplots()
66 66
67 67 widthscreen = widthplot*ncol
68 68 heightscreen = heightplot*nrow
69 69
70 70 return widthscreen, heightscreen
71 71
72 72 def getTimeLim(self, x, xmin=None, xmax=None, timerange=None):
73 73
74 74 if self.xmin != None and self.xmax != None:
75 75 if timerange == None:
76 76 timerange = self.xmax - self.xmin
77 77 xmin = self.xmin + timerange
78 78 xmax = self.xmax + timerange
79 79
80 80 return xmin, xmax
81 81
82 82 if timerange == None and (xmin==None or xmax==None):
83 83 timerange = 14400 #seconds
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))
93 93
94 94 if timerange != None:
95 95 xmin = (thisdatetime - thisdate).seconds/(60*60.)
96 96 xmax = xmin + timerange/(60*60.)
97 97
98 98 mindt = thisdate + datetime.timedelta(hours=xmin) - datetime.timedelta(seconds=time.timezone)
99 99 xmin_sec = time.mktime(mindt.timetuple())
100 100
101 101 maxdt = thisdate + datetime.timedelta(hours=xmax) - datetime.timedelta(seconds=time.timezone)
102 102 xmax_sec = time.mktime(maxdt.timetuple())
103 103
104 104 return xmin_sec, xmax_sec
105 105
106 106 def init(self, id, nplots, wintitle):
107 107
108 108 raise ValueError, "This method has been replaced with createFigure"
109 109
110 110 def createFigure(self, id, wintitle, widthplot=None, heightplot=None, show=True):
111 111
112 112 """
113 113 Crea la figura de acuerdo al driver y parametros seleccionados seleccionados.
114 114 Las dimensiones de la pantalla es calculada a partir de los atributos self.WIDTH
115 115 y self.HEIGHT y el numero de subplots (nrow, ncol)
116 116
117 117 Input:
118 118 id : Los parametros necesarios son
119 119 wintitle :
120 120
121 121 """
122 122
123 123 if widthplot == None:
124 124 widthplot = self.WIDTH
125 125
126 126 if heightplot == None:
127 127 heightplot = self.HEIGHT
128 128
129 129 self.id = id
130 130
131 131 self.wintitle = wintitle
132 132
133 133 self.widthscreen, self.heightscreen = self.getScreenDim(widthplot, heightplot)
134 134
135 135 self.fig = self.__driver.createFigure(id=self.id,
136 136 wintitle=self.wintitle,
137 137 width=self.widthscreen,
138 138 height=self.heightscreen,
139 139 show=show)
140 140
141 141 self.axesObjList = []
142 142 self.counter_imagwr = 0
143 143
144 144
145 145 def setDriver(self, driver=mpldriver):
146 146
147 147 self.__driver = driver
148 148
149 149 def setTitle(self, title):
150 150
151 151 self.__driver.setTitle(self.fig, title)
152 152
153 153 def setWinTitle(self, title):
154 154
155 155 self.__driver.setWinTitle(self.fig, title=title)
156 156
157 157 def setTextFromAxes(self, text):
158 158
159 159 raise ValueError, "Este metodo ha sido reemplazaado con el metodo setText de la clase Axes"
160 160
161 161 def makeAxes(self, nrow, ncol, xpos, ypos, colspan, rowspan):
162 162
163 163 raise ValueError, "Este metodo ha sido reemplazaado con el metodo addAxes"
164 164
165 165 def addAxes(self, *args):
166 166 """
167 167
168 168 Input:
169 169 *args : Los parametros necesarios son
170 170 nrow, ncol, xpos, ypos, colspan, rowspan
171 171 """
172 172
173 173 axesObj = Axes(self.fig, *args)
174 174 self.axesObjList.append(axesObj)
175 175
176 176 def saveFigure(self, figpath, figfile, *args):
177 177
178 178 filename = os.path.join(figpath, figfile)
179 179
180 180 fullpath = os.path.split(filename)[0]
181 181
182 182 if not os.path.exists(fullpath):
183 183 subpath = os.path.split(fullpath)[0]
184 184
185 185 if not os.path.exists(subpath):
186 186 os.mkdir(subpath)
187 187
188 188 os.mkdir(fullpath)
189 189
190 190 self.__driver.saveFigure(self.fig, filename, *args)
191 191
192 192 def save(self, figpath, figfile=None, save=True, ftp=False, wr_period=1, thisDatetime=None, update_figfile=True):
193 193
194 194 self.counter_imagwr += 1
195 195 if self.counter_imagwr < wr_period:
196 196 return
197 197
198 198 self.counter_imagwr = 0
199 199
200 200 if save:
201 201
202 202 if not figfile:
203 203
204 204 if not thisDatetime:
205 205 raise ValueError, "Saving figure: figfile or thisDatetime should be defined"
206 206 return
207 207
208 208 str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S")
209 209 figfile = self.getFilename(name = str_datetime)
210 210
211 211 if self.figfile == None:
212 212 self.figfile = figfile
213 213
214 214 if update_figfile:
215 215 self.figfile = figfile
216 216
217 217 # store png plot to local folder
218 218 self.saveFigure(figpath, self.figfile)
219 219
220 220
221 221 if not ftp:
222 222 return
223 223
224 224 if not thisDatetime:
225 225 return
226 226
227 227 # store png plot to FTP server according to RT-Web format
228 228 ftp_filename = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS)
229 229 # ftp_filename = os.path.join(figpath, name)
230 230 self.saveFigure(figpath, ftp_filename)
231 231
232 232 def getNameToFtp(self, thisDatetime, FTP_WEI, EXP_CODE, SUB_EXP_CODE, PLOT_CODE, PLOT_POS):
233 233 YEAR_STR = '%4.4d'%thisDatetime.timetuple().tm_year
234 234 DOY_STR = '%3.3d'%thisDatetime.timetuple().tm_yday
235 235 FTP_WEI = '%2.2d'%FTP_WEI
236 236 EXP_CODE = '%3.3d'%EXP_CODE
237 237 SUB_EXP_CODE = '%2.2d'%SUB_EXP_CODE
238 238 PLOT_CODE = '%2.2d'%PLOT_CODE
239 239 PLOT_POS = '%2.2d'%PLOT_POS
240 240 name = YEAR_STR + DOY_STR + FTP_WEI + EXP_CODE + SUB_EXP_CODE + PLOT_CODE + PLOT_POS
241 241 return name
242 242
243 243 def draw(self):
244 244
245 245 self.__driver.draw(self.fig)
246 246
247 247 def run(self):
248 248
249 249 raise ValueError, "This method is not implemented"
250 250
251 251 def close(self, show=False):
252 252
253 253 self.__driver.closeFigure(show=show, fig=self.fig)
254 254
255 255 axesList = property(getAxesObjList)
256 256
257 257
258 258 class Axes:
259 259
260 260 __driver = mpldriver
261 261 fig = None
262 262 ax = None
263 263 plot = None
264 264 __missing = 1E30
265 265 __firsttime = None
266 266
267 267 __showprofile = False
268 268
269 269 xmin = None
270 270 xmax = None
271 271 ymin = None
272 272 ymax = None
273 273 zmin = None
274 274 zmax = None
275 275
276 276 x_buffer = None
277 277 z_buffer = None
278 278
279 279 decimationx = None
280 280 decimationy = None
281 281
282 282 __MAXNUMX = 300
283 283 __MAXNUMY = 150
284 284
285 285 def __init__(self, *args):
286 286
287 287 """
288 288
289 289 Input:
290 290 *args : Los parametros necesarios son
291 291 fig, nrow, ncol, xpos, ypos, colspan, rowspan
292 292 """
293 293
294 294 ax = self.__driver.createAxes(*args)
295 295 self.fig = args[0]
296 296 self.ax = ax
297 297 self.plot = None
298 298
299 299 self.__firsttime = True
300 300 self.idlineList = []
301 301
302 302 self.x_buffer = numpy.array([])
303 303 self.z_buffer = numpy.array([])
304 304
305 305 def setText(self, text):
306 306
307 307 self.__driver.setAxesText(self.ax, text)
308 308
309 309 def setXAxisAsTime(self):
310 310 pass
311 311
312 312 def pline(self, x, y,
313 313 xmin=None, xmax=None,
314 314 ymin=None, ymax=None,
315 315 xlabel='', ylabel='',
316 316 title='',
317 317 **kwargs):
318 318
319 319 """
320 320
321 321 Input:
322 322 x :
323 323 y :
324 324 xmin :
325 325 xmax :
326 326 ymin :
327 327 ymax :
328 328 xlabel :
329 329 ylabel :
330 330 title :
331 331 **kwargs : Los parametros aceptados son
332 332
333 333 ticksize
334 334 ytick_visible
335 335 """
336 336
337 337 if self.__firsttime:
338 338
339 339 if xmin == None: xmin = numpy.nanmin(x)
340 340 if xmax == None: xmax = numpy.nanmax(x)
341 341 if ymin == None: ymin = numpy.nanmin(y)
342 342 if ymax == None: ymax = numpy.nanmax(y)
343 343
344 344 self.plot = self.__driver.createPline(self.ax, x, y,
345 345 xmin, xmax,
346 346 ymin, ymax,
347 347 xlabel=xlabel,
348 348 ylabel=ylabel,
349 349 title=title,
350 350 **kwargs)
351 351
352 352 self.idlineList.append(0)
353 353 self.__firsttime = False
354 354 return
355 355
356 356 self.__driver.pline(self.plot, x, y, xlabel=xlabel,
357 357 ylabel=ylabel,
358 358 title=title)
359 359
360 360 def addpline(self, x, y, idline, **kwargs):
361 361 lines = self.ax.lines
362 362
363 363 if idline in self.idlineList:
364 364 self.__driver.set_linedata(self.ax, x, y, idline)
365 365
366 366 if idline not in(self.idlineList):
367 367 self.__driver.addpline(self.ax, x, y, **kwargs)
368 368 self.idlineList.append(idline)
369 369
370 370 return
371 371
372 372 def pmultiline(self, x, y,
373 373 xmin=None, xmax=None,
374 374 ymin=None, ymax=None,
375 375 xlabel='', ylabel='',
376 376 title='',
377 377 **kwargs):
378 378
379 379 if self.__firsttime:
380 380
381 381 if xmin == None: xmin = numpy.nanmin(x)
382 382 if xmax == None: xmax = numpy.nanmax(x)
383 383 if ymin == None: ymin = numpy.nanmin(y)
384 384 if ymax == None: ymax = numpy.nanmax(y)
385 385
386 386 self.plot = self.__driver.createPmultiline(self.ax, x, y,
387 387 xmin, xmax,
388 388 ymin, ymax,
389 389 xlabel=xlabel,
390 390 ylabel=ylabel,
391 391 title=title,
392 392 **kwargs)
393 393 self.__firsttime = False
394 394 return
395 395
396 396 self.__driver.pmultiline(self.plot, x, y, xlabel=xlabel,
397 397 ylabel=ylabel,
398 398 title=title)
399 399
400 400 def pmultilineyaxis(self, x, y,
401 401 xmin=None, xmax=None,
402 402 ymin=None, ymax=None,
403 403 xlabel='', ylabel='',
404 404 title='',
405 405 **kwargs):
406 406
407 407 if self.__firsttime:
408 408
409 409 if xmin == None: xmin = numpy.nanmin(x)
410 410 if xmax == None: xmax = numpy.nanmax(x)
411 411 if ymin == None: ymin = numpy.nanmin(y)
412 412 if ymax == None: ymax = numpy.nanmax(y)
413 413
414 414 self.plot = self.__driver.createPmultilineYAxis(self.ax, x, y,
415 415 xmin, xmax,
416 416 ymin, ymax,
417 417 xlabel=xlabel,
418 418 ylabel=ylabel,
419 419 title=title,
420 420 **kwargs)
421 421 if self.xmin == None: self.xmin = xmin
422 422 if self.xmax == None: self.xmax = xmax
423 423 if self.ymin == None: self.ymin = ymin
424 424 if self.ymax == None: self.ymax = ymax
425 425
426 426 self.__firsttime = False
427 427 return
428 428
429 429 self.__driver.pmultilineyaxis(self.plot, x, y, xlabel=xlabel,
430 430 ylabel=ylabel,
431 431 title=title)
432 432
433 433 def pcolor(self, x, y, z,
434 434 xmin=None, xmax=None,
435 435 ymin=None, ymax=None,
436 436 zmin=None, zmax=None,
437 437 xlabel='', ylabel='',
438 438 title='', rti = False, colormap='jet',
439 439 **kwargs):
440 440
441 441 """
442 442 Input:
443 443 x :
444 444 y :
445 445 x :
446 446 xmin :
447 447 xmax :
448 448 ymin :
449 449 ymax :
450 450 zmin :
451 451 zmax :
452 452 xlabel :
453 453 ylabel :
454 454 title :
455 455 **kwargs : Los parametros aceptados son
456 456 ticksize=9,
457 457 cblabel=''
458 458 rti = True or False
459 459 """
460 460
461 461 if self.__firsttime:
462 462
463 463 if xmin == None: xmin = numpy.nanmin(x)
464 464 if xmax == None: xmax = numpy.nanmax(x)
465 465 if ymin == None: ymin = numpy.nanmin(y)
466 466 if ymax == None: ymax = numpy.nanmax(y)
467 467 if zmin == None: zmin = numpy.nanmin(z)
468 468 if zmax == None: zmax = numpy.nanmax(z)
469 469
470 470
471 471 self.plot = self.__driver.createPcolor(self.ax, x, y, z,
472 472 xmin, xmax,
473 473 ymin, ymax,
474 474 zmin, zmax,
475 475 xlabel=xlabel,
476 476 ylabel=ylabel,
477 477 title=title,
478 478 colormap=colormap,
479 479 **kwargs)
480 480
481 481 if self.xmin == None: self.xmin = xmin
482 482 if self.xmax == None: self.xmax = xmax
483 483 if self.ymin == None: self.ymin = ymin
484 484 if self.ymax == None: self.ymax = ymax
485 485 if self.zmin == None: self.zmin = zmin
486 486 if self.zmax == None: self.zmax = zmax
487 487
488 488 self.__firsttime = False
489 489 return
490 490
491 491 if rti:
492 492 self.__driver.addpcolor(self.ax, x, y, z, self.zmin, self.zmax,
493 493 xlabel=xlabel,
494 494 ylabel=ylabel,
495 495 title=title,
496 496 colormap=colormap)
497 497 return
498 498
499 499 self.__driver.pcolor(self.plot, z,
500 500 xlabel=xlabel,
501 501 ylabel=ylabel,
502 502 title=title)
503 503
504 504 def pcolorbuffer(self, x, y, z,
505 505 xmin=None, xmax=None,
506 506 ymin=None, ymax=None,
507 507 zmin=None, zmax=None,
508 508 xlabel='', ylabel='',
509 509 title='', rti = True, colormap='jet',
510 510 maxNumX = None, maxNumY = None,
511 511 **kwargs):
512 512
513 513 if maxNumX == None:
514 514 maxNumX = self.__MAXNUMX
515 515
516 516 if maxNumY == None:
517 517 maxNumY = self.__MAXNUMY
518 518
519 519 if self.__firsttime:
520 520 self.z_buffer = z
521 521 self.x_buffer = numpy.hstack((self.x_buffer, x))
522 522
523 523 if xmin == None: xmin = numpy.nanmin(x)
524 524 if xmax == None: xmax = numpy.nanmax(x)
525 525 if ymin == None: ymin = numpy.nanmin(y)
526 526 if ymax == None: ymax = numpy.nanmax(y)
527 527 if zmin == None: zmin = numpy.nanmin(z)
528 528 if zmax == None: zmax = numpy.nanmax(z)
529 529
530 530
531 531 self.plot = self.__driver.createPcolor(self.ax, self.x_buffer, y, z,
532 532 xmin, xmax,
533 533 ymin, ymax,
534 534 zmin, zmax,
535 535 xlabel=xlabel,
536 536 ylabel=ylabel,
537 537 title=title,
538 538 colormap=colormap,
539 539 **kwargs)
540 540
541 541 if self.xmin == None: self.xmin = xmin
542 542 if self.xmax == None: self.xmax = xmax
543 543 if self.ymin == None: self.ymin = ymin
544 544 if self.ymax == None: self.ymax = ymax
545 545 if self.zmin == None: self.zmin = zmin
546 546 if self.zmax == None: self.zmax = zmax
547 547
548 548 self.__firsttime = False
549 549 return
550 550
551 551 self.x_buffer = numpy.hstack((self.x_buffer, x[-1]))
552 552 self.z_buffer = numpy.hstack((self.z_buffer, z))
553 553
554 554 if self.decimationx == None:
555 555 deltax = float(self.xmax - self.xmin)/maxNumX
556 556 deltay = float(self.ymax - self.ymin)/maxNumY
557 557
558 558 resolutionx = self.x_buffer[2]-self.x_buffer[0]
559 559 resolutiony = y[1]-y[0]
560 560
561 561 self.decimationx = numpy.ceil(deltax / resolutionx)
562 562 self.decimationy = numpy.ceil(deltay / resolutiony)
563 563
564 564 z_buffer = self.z_buffer.reshape(-1,len(y))
565 565
566 566 x_buffer = self.x_buffer[::self.decimationx]
567 567 y_buffer = y[::self.decimationy]
568 568 z_buffer = z_buffer[::self.decimationx, ::self.decimationy]
569 569 #===================================================
570 570
571 571 x_buffer, y_buffer, z_buffer = self.__fillGaps(x_buffer, y_buffer, z_buffer)
572 572
573 573 self.__driver.addpcolorbuffer(self.ax, x_buffer, y_buffer, z_buffer, self.zmin, self.zmax,
574 574 xlabel=xlabel,
575 575 ylabel=ylabel,
576 576 title=title,
577 577 colormap=colormap)
578 578
579 579 def polar(self, x, y,
580 580 title='', xlabel='',ylabel='',**kwargs):
581 581
582 582 if self.__firsttime:
583 583 self.plot = self.__driver.createPolar(self.ax, x, y, title = title, xlabel = xlabel, ylabel = ylabel)
584 584 self.__firsttime = False
585 585 self.x_buffer = x
586 586 self.y_buffer = y
587 587 return
588 588
589 589 self.x_buffer = numpy.hstack((self.x_buffer,x))
590 590 self.y_buffer = numpy.hstack((self.y_buffer,y))
591 591 self.__driver.polar(self.plot, self.x_buffer, self.y_buffer, xlabel=xlabel,
592 592 ylabel=ylabel,
593 593 title=title)
594 594
595 595 def __fillGaps(self, x_buffer, y_buffer, z_buffer):
596 596
597 597 deltas = x_buffer[1:] - x_buffer[0:-1]
598 598 x_median = numpy.median(deltas)
599 599
600 600 index = numpy.where(deltas >= 2*x_median)
601 601
602 602 if len(index[0]) != 0:
603 603 z_buffer[index[0],::] = self.__missing
604 604 z_buffer = numpy.ma.masked_inside(z_buffer,0.99*self.__missing,1.01*self.__missing)
605 605
606 606 return x_buffer, y_buffer, z_buffer
607 607
608 608
609 609
610 610 No newline at end of file
@@ -1,1164 +1,1362
1 1 import os
2 2 import datetime
3 3 import numpy
4 4
5 5 from figure import Figure, isRealtime
6 6 from plotting_codes import *
7 7
8 8 class MomentsPlot(Figure):
9 9
10 10 isConfig = None
11 11 __nsubplots = None
12 12
13 13 WIDTHPROF = None
14 14 HEIGHTPROF = None
15 15 PREFIX = 'prm'
16 16
17 17 def __init__(self):
18 18
19 19 self.isConfig = False
20 20 self.__nsubplots = 1
21 21
22 22 self.WIDTH = 280
23 23 self.HEIGHT = 250
24 24 self.WIDTHPROF = 120
25 25 self.HEIGHTPROF = 0
26 26 self.counter_imagwr = 0
27 27
28 28 self.PLOT_CODE = MOMENTS_CODE
29 29
30 30 self.FTP_WEI = None
31 31 self.EXP_CODE = None
32 32 self.SUB_EXP_CODE = None
33 33 self.PLOT_POS = None
34 34
35 35 def getSubplots(self):
36 36
37 37 ncol = int(numpy.sqrt(self.nplots)+0.9)
38 38 nrow = int(self.nplots*1./ncol + 0.9)
39 39
40 40 return nrow, ncol
41 41
42 42 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
43 43
44 44 self.__showprofile = showprofile
45 45 self.nplots = nplots
46 46
47 47 ncolspan = 1
48 48 colspan = 1
49 49 if showprofile:
50 50 ncolspan = 3
51 51 colspan = 2
52 52 self.__nsubplots = 2
53 53
54 54 self.createFigure(id = id,
55 55 wintitle = wintitle,
56 56 widthplot = self.WIDTH + self.WIDTHPROF,
57 57 heightplot = self.HEIGHT + self.HEIGHTPROF,
58 58 show=show)
59 59
60 60 nrow, ncol = self.getSubplots()
61 61
62 62 counter = 0
63 63 for y in range(nrow):
64 64 for x in range(ncol):
65 65
66 66 if counter >= self.nplots:
67 67 break
68 68
69 69 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
70 70
71 71 if showprofile:
72 72 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
73 73
74 74 counter += 1
75 75
76 76 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True,
77 77 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
78 78 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
79 79 server=None, folder=None, username=None, password=None,
80 80 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
81 81
82 82 """
83 83
84 84 Input:
85 85 dataOut :
86 86 id :
87 87 wintitle :
88 88 channelList :
89 89 showProfile :
90 90 xmin : None,
91 91 xmax : None,
92 92 ymin : None,
93 93 ymax : None,
94 94 zmin : None,
95 95 zmax : None
96 96 """
97 97
98 98 if dataOut.flagNoData:
99 99 return None
100 100
101 101 if realtime:
102 102 if not(isRealtime(utcdatatime = dataOut.utctime)):
103 103 print 'Skipping this plot function'
104 104 return
105 105
106 106 if channelList == None:
107 107 channelIndexList = dataOut.channelIndexList
108 108 else:
109 109 channelIndexList = []
110 110 for channel in channelList:
111 111 if channel not in dataOut.channelList:
112 112 raise ValueError, "Channel %d is not in dataOut.channelList"
113 113 channelIndexList.append(dataOut.channelList.index(channel))
114 114
115 115 factor = dataOut.normFactor
116 116 x = dataOut.abscissaList
117 117 y = dataOut.heightList
118 118
119 119 z = dataOut.data_pre[channelIndexList,:,:]/factor
120 120 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
121 121 avg = numpy.average(z, axis=1)
122 122 noise = dataOut.noise/factor
123 123
124 124 zdB = 10*numpy.log10(z)
125 125 avgdB = 10*numpy.log10(avg)
126 126 noisedB = 10*numpy.log10(noise)
127 127
128 128 #thisDatetime = dataOut.datatime
129 129 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
130 130 title = wintitle + " Parameters"
131 131 xlabel = "Velocity (m/s)"
132 132 ylabel = "Range (Km)"
133 133
134 134 if not self.isConfig:
135 135
136 136 nplots = len(channelIndexList)
137 137
138 138 self.setup(id=id,
139 139 nplots=nplots,
140 140 wintitle=wintitle,
141 141 showprofile=showprofile,
142 142 show=show)
143 143
144 144 if xmin == None: xmin = numpy.nanmin(x)
145 145 if xmax == None: xmax = numpy.nanmax(x)
146 146 if ymin == None: ymin = numpy.nanmin(y)
147 147 if ymax == None: ymax = numpy.nanmax(y)
148 148 if zmin == None: zmin = numpy.nanmin(avgdB)*0.9
149 149 if zmax == None: zmax = numpy.nanmax(avgdB)*0.9
150 150
151 151 self.FTP_WEI = ftp_wei
152 152 self.EXP_CODE = exp_code
153 153 self.SUB_EXP_CODE = sub_exp_code
154 154 self.PLOT_POS = plot_pos
155 155
156 156 self.isConfig = True
157 157
158 158 self.setWinTitle(title)
159 159
160 160 for i in range(self.nplots):
161 161 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
162 162 title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i]+1, noisedB[i], str_datetime)
163 163 axes = self.axesList[i*self.__nsubplots]
164 164 axes.pcolor(x, y, zdB[i,:,:],
165 165 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
166 166 xlabel=xlabel, ylabel=ylabel, title=title,
167 167 ticksize=9, cblabel='')
168 168 #Mean Line
169 169 mean = dataOut.data_param[i, 1, :]
170 170 axes.addpline(mean, y, idline=0, color="black", linestyle="solid", lw=1)
171 171
172 172 if self.__showprofile:
173 173 axes = self.axesList[i*self.__nsubplots +1]
174 174 axes.pline(avgdB[i], y,
175 175 xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax,
176 176 xlabel='dB', ylabel='', title='',
177 177 ytick_visible=False,
178 178 grid='x')
179 179
180 180 noiseline = numpy.repeat(noisedB[i], len(y))
181 181 axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2)
182 182
183 183 self.draw()
184 184
185 185 self.save(figpath=figpath,
186 186 figfile=figfile,
187 187 save=save,
188 188 ftp=ftp,
189 189 wr_period=wr_period,
190 190 thisDatetime=thisDatetime)
191 191
192 192
193 193
194 194 class SkyMapPlot(Figure):
195 195
196 196 __isConfig = None
197 197 __nsubplots = None
198 198
199 199 WIDTHPROF = None
200 200 HEIGHTPROF = None
201 PREFIX = 'prm'
201 PREFIX = 'mmap'
202 202
203 203 def __init__(self):
204 204
205 205 self.__isConfig = False
206 206 self.__nsubplots = 1
207 207
208 208 # self.WIDTH = 280
209 209 # self.HEIGHT = 250
210 210 self.WIDTH = 600
211 211 self.HEIGHT = 600
212 212 self.WIDTHPROF = 120
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
220 220 self.SUB_EXP_CODE = None
221 221 self.PLOT_POS = None
222 222
223 223 def getSubplots(self):
224 224
225 225 ncol = int(numpy.sqrt(self.nplots)+0.9)
226 226 nrow = int(self.nplots*1./ncol + 0.9)
227 227
228 228 return nrow, ncol
229 229
230 230 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
231 231
232 232 self.__showprofile = showprofile
233 233 self.nplots = nplots
234 234
235 235 ncolspan = 1
236 236 colspan = 1
237 237
238 238 self.createFigure(id = id,
239 239 wintitle = wintitle,
240 240 widthplot = self.WIDTH, #+ self.WIDTHPROF,
241 241 heightplot = self.HEIGHT,# + self.HEIGHTPROF,
242 242 show=show)
243 243
244 244 nrow, ncol = 1,1
245 245 counter = 0
246 246 x = 0
247 247 y = 0
248 248 self.addAxes(1, 1, 0, 0, 1, 1, True)
249 249
250 250 def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False,
251 251 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
252 252 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
253 253 server=None, folder=None, username=None, password=None,
254 254 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False):
255 255
256 256 """
257 257
258 258 Input:
259 259 dataOut :
260 260 id :
261 261 wintitle :
262 262 channelList :
263 263 showProfile :
264 264 xmin : None,
265 265 xmax : None,
266 266 ymin : None,
267 267 ymax : None,
268 268 zmin : None,
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,:]
276 276 finalAzimuth = finalMeteor[:,4]
277 277 finalZenith = finalMeteor[:,5]
278 278
279 279 x = finalAzimuth*numpy.pi/180
280 280 y = finalZenith
281 281
282 282
283 283 #thisDatetime = dataOut.datatime
284 284 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
285 285 title = wintitle + " Parameters"
286 286 xlabel = "Zonal Zenith Angle (deg) "
287 287 ylabel = "Meridional Zenith Angle (deg)"
288 288
289 289 if not self.__isConfig:
290 290
291 291 nplots = 1
292 292
293 293 self.setup(id=id,
294 294 nplots=nplots,
295 295 wintitle=wintitle,
296 296 showprofile=showprofile,
297 297 show=show)
298 298
299 299 self.FTP_WEI = ftp_wei
300 300 self.EXP_CODE = exp_code
301 301 self.SUB_EXP_CODE = sub_exp_code
302 302 self.PLOT_POS = plot_pos
303 303 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
304 304 self.firstdate = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
305 305 self.__isConfig = True
306 306
307 307 self.setWinTitle(title)
308 308
309 309 i = 0
310 310 str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S"))
311 311
312 312 axes = self.axesList[i*self.__nsubplots]
313 313 nevents = axes.x_buffer.shape[0] + x.shape[0]
314 314 title = "Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n" %(self.firstdate,str_datetime,nevents)
315 315 axes.polar(x, y,
316 316 title=title, xlabel=xlabel, ylabel=ylabel,
317 317 ticksize=9, cblabel='')
318 318
319 319 self.draw()
320 320
321 321 self.save(figpath=figpath,
322 322 figfile=figfile,
323 323 save=save,
324 324 ftp=ftp,
325 325 wr_period=wr_period,
326 326 thisDatetime=thisDatetime)
327
327 328
328 329 class WindProfilerPlot(Figure):
329 330
330 331 __isConfig = None
331 332 __nsubplots = None
332 333
333 334 WIDTHPROF = None
334 335 HEIGHTPROF = None
335 336 PREFIX = 'wind'
336 337
337 338 def __init__(self):
338 339
339 self.timerange = 2*60*60
340 self.timerange = None
340 341 self.__isConfig = False
341 342 self.__nsubplots = 1
342 343
343 344 self.WIDTH = 800
344 345 self.HEIGHT = 150
345 346 self.WIDTHPROF = 120
346 347 self.HEIGHTPROF = 0
347 348 self.counter_imagwr = 0
348 349
349 350 self.PLOT_CODE = WIND_CODE
350 351
351 352 self.FTP_WEI = None
352 353 self.EXP_CODE = None
353 354 self.SUB_EXP_CODE = None
354 355 self.PLOT_POS = None
355 356 self.tmin = None
356 357 self.tmax = None
357 358
358 359 self.xmin = None
359 360 self.xmax = None
360 361
361 362 self.figfile = None
362 363
363 364 def getSubplots(self):
364 365
365 366 ncol = 1
366 367 nrow = self.nplots
367 368
368 369 return nrow, ncol
369 370
370 371 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
371 372
372 373 self.__showprofile = showprofile
373 374 self.nplots = nplots
374 375
375 376 ncolspan = 1
376 377 colspan = 1
377 378
378 379 self.createFigure(id = id,
379 380 wintitle = wintitle,
380 381 widthplot = self.WIDTH + self.WIDTHPROF,
381 382 heightplot = self.HEIGHT + self.HEIGHTPROF,
382 383 show=show)
383 384
384 385 nrow, ncol = self.getSubplots()
385 386
386 387 counter = 0
387 388 for y in range(nrow):
388 389 if counter >= self.nplots:
389 390 break
390 391
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,
398 399 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
399 400 server=None, folder=None, username=None, password=None,
400 401 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
401 402 """
402 403
403 404 Input:
404 405 dataOut :
405 406 id :
406 407 wintitle :
407 408 channelList :
408 409 showProfile :
409 410 xmin : None,
410 411 xmax : None,
411 412 ymin : None,
412 413 ymax : None,
413 414 zmin : None,
414 415 zmax : None
415 416 """
416 417
417 418 if channelList == None:
418 419 channelIndexList = dataOut.channelIndexList
419 420 else:
420 421 channelIndexList = []
421 422 for channel in channelList:
422 423 if channel not in dataOut.channelList:
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
434 435 y = dataOut.heightList
435 436
436 437 z = dataOut.data_output.copy()
437 438 nplots = z.shape[0] #Number of wind dimensions estimated
438 439 nplotsw = nplots
439 440
440 441 #If there is a SNR function defined
441 442 if dataOut.data_SNR != None:
442 443 nplots += 1
443 444 SNR = dataOut.data_SNR
444 445 SNRavg = numpy.average(SNR, axis=0)
445 446
446 447 SNRdB = 10*numpy.log10(SNR)
447 448 SNRavgdB = 10*numpy.log10(SNRavg)
448 449
449 450 if SNRthresh == None: SNRthresh = -5.0
450 451 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
451 452
452 453 for i in range(nplotsw):
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"
460 461 xlabel = ""
461 462 ylabel = "Range (Km)"
462 463
463 464 if not self.__isConfig:
464 465
465 466 self.setup(id=id,
466 467 nplots=nplots,
467 468 wintitle=wintitle,
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)
474 478 if ymax == None: ymax = numpy.nanmax(y)
475 479
476 480 if zmax == None: zmax = numpy.nanmax(abs(z[range(2),:]))
477 481 #if numpy.isnan(zmax): zmax = 50
478 482 if zmin == None: zmin = -zmax
479 483
480 484 if nplotsw == 3:
481 485 if zmax_ver == None: zmax_ver = numpy.nanmax(abs(z[2,:]))
482 486 if zmin_ver == None: zmin_ver = -zmax_ver
483 487
484 488 if dataOut.data_SNR != None:
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
491 496 self.PLOT_POS = plot_pos
492
497
493 498 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
494 499 self.__isConfig = True
495
496
500 self.figfile = figfile
501
497 502 self.setWinTitle(title)
498
503
499 504 if ((self.xmax - x[1]) < (x[1]-x[0])):
500 505 x[1] = self.xmax
501 506
502 507 strWind = ['Zonal', 'Meridional', 'Vertical']
503 508 strCb = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)']
504 509 zmaxVector = [zmax, zmax, zmax_ver]
505 510 zminVector = [zmin, zmin, zmin_ver]
506 511 windFactor = [1,1,100]
507 512
508 513 for i in range(nplotsw):
509 514
510 515 title = "%s Wind: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
511 516 axes = self.axesList[i*self.__nsubplots]
512 517
513 518 z1 = z[i,:].reshape((1,-1))*windFactor[i]
514 519
515 520 axes.pcolorbuffer(x, y, z1,
516 521 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
517 522 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
518 523 ticksize=9, cblabel=strCb[i], cbsize="1%", colormap="RdBu_r" )
519 524
520 525 if dataOut.data_SNR != None:
521 526 i += 1
522 527 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
523 528 axes = self.axesList[i*self.__nsubplots]
524 529
525 530 SNRavgdB = SNRavgdB.reshape((1,-1))
526 531
527 532 axes.pcolorbuffer(x, y, SNRavgdB,
528 533 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
529 534 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
530 535 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
531 536
532 537 self.draw()
533 538
534 539 if x[1] >= self.axesList[0].xmax:
535 540 self.counter_imagwr = wr_period
536 541 self.__isConfig = False
537 542 self.figfile = None
538 543
539 544 self.save(figpath=figpath,
540 545 figfile=figfile,
541 546 save=save,
542 547 ftp=ftp,
543 548 wr_period=wr_period,
544 549 thisDatetime=thisDatetime,
545 550 update_figfile=False)
546 551
547 552
548 553 class ParametersPlot(Figure):
549 554
550 555 __isConfig = None
551 556 __nsubplots = None
552 557
553 558 WIDTHPROF = None
554 559 HEIGHTPROF = None
555 560 PREFIX = 'prm'
556 561
557 562 def __init__(self):
558 563
559 564 self.timerange = 2*60*60
560 565 self.__isConfig = False
561 566 self.__nsubplots = 1
562 567
563 568 self.WIDTH = 800
564 569 self.HEIGHT = 150
565 570 self.WIDTHPROF = 120
566 571 self.HEIGHTPROF = 0
567 572 self.counter_imagwr = 0
568 573
569 574 self.PLOT_CODE = PARMS_CODE
570 575
571 576 self.FTP_WEI = None
572 577 self.EXP_CODE = None
573 578 self.SUB_EXP_CODE = None
574 579 self.PLOT_POS = None
575 580 self.tmin = None
576 581 self.tmax = None
577 582
578 583 self.xmin = None
579 584 self.xmax = None
580 585
581 586 self.figfile = None
582 587
583 588 def getSubplots(self):
584 589
585 590 ncol = 1
586 591 nrow = self.nplots
587 592
588 593 return nrow, ncol
589 594
590 595 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
591 596
592 597 self.__showprofile = showprofile
593 598 self.nplots = nplots
594 599
595 600 ncolspan = 1
596 601 colspan = 1
597 602
598 603 self.createFigure(id = id,
599 604 wintitle = wintitle,
600 605 widthplot = self.WIDTH + self.WIDTHPROF,
601 606 heightplot = self.HEIGHT + self.HEIGHTPROF,
602 607 show=show)
603 608
604 609 nrow, ncol = self.getSubplots()
605 610
606 611 counter = 0
607 612 for y in range(nrow):
608 613 for x in range(ncol):
609 614
610 615 if counter >= self.nplots:
611 616 break
612 617
613 618 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
614 619
615 620 if showprofile:
616 621 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
617 622
618 623 counter += 1
619 624
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,
627 632 server=None, folder=None, username=None, password=None,
628 633 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
629 634
630 635 """
631 636
632 637 Input:
633 638 dataOut :
634 639 id :
635 640 wintitle :
636 641 channelList :
637 642 showProfile :
638 643 xmin : None,
639 644 xmax : None,
640 645 ymin : None,
641 646 ymax : None,
642 647 zmin : None,
643 648 zmax : None
644 649 """
645 650
646 651 data_param = getattr(dataOut, parameterObject)
647 652
648 653 if channelList == None:
649 654 channelIndexList = numpy.arange(data_param.shape[0])
650 655 else:
651 656 channelIndexList = numpy.array(channelList)
652 657
653 658 nchan = len(channelIndexList) #Number of channels being plotted
654 659
655 660 if nchan < 1:
656 661 return
657 662
658 663 nGraphsByChannel = 0
659 664
660 665 if SNR:
661 666 nGraphsByChannel += 1
662 667 if DOP:
663 668 nGraphsByChannel += 1
664 669
665 670 if nGraphsByChannel < 1:
666 671 return
667 672
668 673 nplots = nGraphsByChannel*nchan
669 674
670 675 if timerange != None:
671 676 self.timerange = timerange
672 677
673 678 #tmin = None
674 679 #tmax = None
675 680 if parameterIndex == None:
676 681 parameterIndex = 1
677 682
678 683 x = dataOut.getTimeRange1()
679 684 y = dataOut.heightList
680 685 z = data_param[channelIndexList,parameterIndex,:].copy()
681 686
682 687 zRange = dataOut.abscissaList
683 688 # nChannels = z.shape[0] #Number of wind dimensions estimated
684 689 # thisDatetime = dataOut.datatime
685 690
686 691 if dataOut.data_SNR != None:
687 692 SNRarray = dataOut.data_SNR[channelIndexList,:]
688 693 SNRdB = 10*numpy.log10(SNRarray)
689 694 # SNRavgdB = 10*numpy.log10(SNRavg)
690 695 ind = numpy.where(SNRdB < 10**(SNRthresh/10))
691 696 z[ind] = numpy.nan
692
697
693 698 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
694 699 title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y"))
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
701 708 else: colormap = "RdBu_r"
702 709
703 710 if not self.__isConfig:
704 711
705 712 self.setup(id=id,
706 713 nplots=nplots,
707 714 wintitle=wintitle,
708 715 showprofile=showprofile,
709 716 show=show)
710 717
711 718 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
712 719
713 720 if ymin == None: ymin = numpy.nanmin(y)
714 721 if ymax == None: ymax = numpy.nanmax(y)
715 722 if zmin == None: zmin = numpy.nanmin(zRange)
716 723 if zmax == None: zmax = numpy.nanmax(zRange)
717 724
718 725 if SNR:
719 726 if SNRmin == None: SNRmin = numpy.nanmin(SNRdB)
720 727 if SNRmax == None: SNRmax = numpy.nanmax(SNRdB)
721 728
722 729 self.FTP_WEI = ftp_wei
723 730 self.EXP_CODE = exp_code
724 731 self.SUB_EXP_CODE = sub_exp_code
725 732 self.PLOT_POS = plot_pos
726 733
727 734 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
728 735 self.__isConfig = True
729 736 self.figfile = figfile
730 737
731 738 self.setWinTitle(title)
732 739
733 740 if ((self.xmax - x[1]) < (x[1]-x[0])):
734 741 x[1] = self.xmax
735 742
736 743 for i in range(nchan):
744
745 if (SNR and not onlySNR): j = 2*i
746 else: j = i
737 747
738 748 j = nGraphsByChannel*i
749
750 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
751 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
739 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
743 764 if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)):
744 765 title = title + '_' + 'azimuth,zenith=%2.2f,%2.2f'%(dataOut.azimuth, dataOut.zenith)
745 766 axes = self.axesList[j]
746 767 z1 = z[i,:].reshape((1,-1))
747 768 axes.pcolorbuffer(x, y, z1,
748 769 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax,
749 770 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap=colormap,
750 771 ticksize=9, cblabel=zlabel, cbsize="1%")
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,
758 784 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,colormap="jet",
759 785 ticksize=9, cblabel=zlabel, cbsize="1%")
760 786
761 787
762 788
763 789 self.draw()
764 790
765 791 if x[1] >= self.axesList[0].xmax:
766 792 self.counter_imagwr = wr_period
767 793 self.__isConfig = False
768 794 self.figfile = None
769 795
770 796 self.save(figpath=figpath,
771 797 figfile=figfile,
772 798 save=save,
773 799 ftp=ftp,
774 800 wr_period=wr_period,
775 801 thisDatetime=thisDatetime,
776 802 update_figfile=False)
777 803
778 804 class SpectralFittingPlot(Figure):
779 805
780 806 __isConfig = None
781 807 __nsubplots = None
782 808
783 809 WIDTHPROF = None
784 810 HEIGHTPROF = None
785 811 PREFIX = 'prm'
786 812
787 813
788 814 N = None
789 815 ippSeconds = None
790 816
791 817 def __init__(self):
792 818 self.__isConfig = False
793 819 self.__nsubplots = 1
794 820
795 821 self.PLOT_CODE = SPECFIT_CODE
796 822
797 823 self.WIDTH = 450
798 824 self.HEIGHT = 250
799 825 self.WIDTHPROF = 0
800 826 self.HEIGHTPROF = 0
801 827
802 828 def getSubplots(self):
803 829
804 830 ncol = int(numpy.sqrt(self.nplots)+0.9)
805 831 nrow = int(self.nplots*1./ncol + 0.9)
806 832
807 833 return nrow, ncol
808 834
809 835 def setup(self, id, nplots, wintitle, showprofile=False, show=True):
810 836
811 837 showprofile = False
812 838 self.__showprofile = showprofile
813 839 self.nplots = nplots
814 840
815 841 ncolspan = 5
816 842 colspan = 4
817 843 if showprofile:
818 844 ncolspan = 5
819 845 colspan = 4
820 846 self.__nsubplots = 2
821 847
822 848 self.createFigure(id = id,
823 849 wintitle = wintitle,
824 850 widthplot = self.WIDTH + self.WIDTHPROF,
825 851 heightplot = self.HEIGHT + self.HEIGHTPROF,
826 852 show=show)
827 853
828 854 nrow, ncol = self.getSubplots()
829 855
830 856 counter = 0
831 857 for y in range(nrow):
832 858 for x in range(ncol):
833 859
834 860 if counter >= self.nplots:
835 861 break
836 862
837 863 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1)
838 864
839 865 if showprofile:
840 866 self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1)
841 867
842 868 counter += 1
843 869
844 870 def run(self, dataOut, id, cutHeight=None, fit=False, wintitle="", channelList=None, showprofile=True,
845 871 xmin=None, xmax=None, ymin=None, ymax=None,
846 872 save=False, figpath='./', figfile=None, show=True):
847 873
848 874 """
849 875
850 876 Input:
851 877 dataOut :
852 878 id :
853 879 wintitle :
854 880 channelList :
855 881 showProfile :
856 882 xmin : None,
857 883 xmax : None,
858 884 zmin : None,
859 885 zmax : None
860 886 """
861 887
862 888 if cutHeight==None:
863 889 h=270
864 890 heightindex = numpy.abs(cutHeight - dataOut.heightList).argmin()
865 891 cutHeight = dataOut.heightList[heightindex]
866 892
867 893 factor = dataOut.normFactor
868 894 x = dataOut.abscissaList[:-1]
869 895 #y = dataOut.getHeiRange()
870 896
871 897 z = dataOut.data_pre[:,:,heightindex]/factor
872 898 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
873 899 avg = numpy.average(z, axis=1)
874 900 listChannels = z.shape[0]
875 901
876 902 #Reconstruct Function
877 903 if fit==True:
878 904 groupArray = dataOut.groupList
879 905 listChannels = groupArray.reshape((groupArray.size))
880 906 listChannels.sort()
881 907 spcFitLine = numpy.zeros(z.shape)
882 908 constants = dataOut.constants
883 909
884 910 nGroups = groupArray.shape[0]
885 911 nChannels = groupArray.shape[1]
886 912 nProfiles = z.shape[1]
887 913
888 914 for f in range(nGroups):
889 915 groupChann = groupArray[f,:]
890 916 p = dataOut.data_param[f,:,heightindex]
891 917 # p = numpy.array([ 89.343967,0.14036615,0.17086219,18.89835291,1.58388365,1.55099167])
892 918 fitLineAux = dataOut.library.modelFunction(p, constants)*nProfiles
893 919 fitLineAux = fitLineAux.reshape((nChannels,nProfiles))
894 920 spcFitLine[groupChann,:] = fitLineAux
895 921 # spcFitLine = spcFitLine/factor
896 922
897 923 z = z[listChannels,:]
898 924 spcFitLine = spcFitLine[listChannels,:]
899 925 spcFitLinedB = 10*numpy.log10(spcFitLine)
900 926
901 927 zdB = 10*numpy.log10(z)
902 928 #thisDatetime = dataOut.datatime
903 929 thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0])
904 930 title = wintitle + " Doppler Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
905 931 xlabel = "Velocity (m/s)"
906 932 ylabel = "Spectrum"
907 933
908 934 if not self.__isConfig:
909 935
910 936 nplots = listChannels.size
911 937
912 938 self.setup(id=id,
913 939 nplots=nplots,
914 940 wintitle=wintitle,
915 941 showprofile=showprofile,
916 942 show=show)
917 943
918 944 if xmin == None: xmin = numpy.nanmin(x)
919 945 if xmax == None: xmax = numpy.nanmax(x)
920 946 if ymin == None: ymin = numpy.nanmin(zdB)
921 947 if ymax == None: ymax = numpy.nanmax(zdB)+2
922 948
923 949 self.__isConfig = True
924 950
925 951 self.setWinTitle(title)
926 952 for i in range(self.nplots):
927 953 # title = "Channel %d: %4.2fdB" %(dataOut.channelList[i]+1, noisedB[i])
928 954 title = "Height %4.1f km\nChannel %d:" %(cutHeight, listChannels[i]+1)
929 955 axes = self.axesList[i*self.__nsubplots]
930 956 if fit == False:
931 957 axes.pline(x, zdB[i,:],
932 958 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
933 959 xlabel=xlabel, ylabel=ylabel, title=title
934 960 )
935 961 if fit == True:
936 962 fitline=spcFitLinedB[i,:]
937 963 y=numpy.vstack([zdB[i,:],fitline] )
938 964 legendlabels=['Data','Fitting']
939 965 axes.pmultilineyaxis(x, y,
940 966 xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax,
941 967 xlabel=xlabel, ylabel=ylabel, title=title,
942 968 legendlabels=legendlabels, marker=None,
943 969 linestyle='solid', grid='both')
944 970
945 971 self.draw()
946 972
947 973 self.save(figpath=figpath,
948 974 figfile=figfile,
949 975 save=save,
950 976 ftp=ftp,
951 977 wr_period=wr_period,
952 978 thisDatetime=thisDatetime)
953 979
954 980
955 981 class EWDriftsPlot(Figure):
956 982
957 983 __isConfig = None
958 984 __nsubplots = None
959 985
960 986 WIDTHPROF = None
961 987 HEIGHTPROF = None
962 988 PREFIX = 'drift'
963 989
964 990 def __init__(self):
965 991
966 992 self.timerange = 2*60*60
967 993 self.isConfig = False
968 994 self.__nsubplots = 1
969 995
970 996 self.WIDTH = 800
971 997 self.HEIGHT = 150
972 998 self.WIDTHPROF = 120
973 999 self.HEIGHTPROF = 0
974 1000 self.counter_imagwr = 0
975 1001
976 1002 self.PLOT_CODE = EWDRIFT_CODE
977 1003
978 1004 self.FTP_WEI = None
979 1005 self.EXP_CODE = None
980 1006 self.SUB_EXP_CODE = None
981 1007 self.PLOT_POS = None
982 1008 self.tmin = None
983 1009 self.tmax = None
984 1010
985 1011 self.xmin = None
986 1012 self.xmax = None
987 1013
988 1014 self.figfile = None
989 1015
990 1016 def getSubplots(self):
991 1017
992 1018 ncol = 1
993 1019 nrow = self.nplots
994 1020
995 1021 return nrow, ncol
996 1022
997 1023 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
998 1024
999 1025 self.__showprofile = showprofile
1000 1026 self.nplots = nplots
1001 1027
1002 1028 ncolspan = 1
1003 1029 colspan = 1
1004 1030
1005 1031 self.createFigure(id = id,
1006 1032 wintitle = wintitle,
1007 1033 widthplot = self.WIDTH + self.WIDTHPROF,
1008 1034 heightplot = self.HEIGHT + self.HEIGHTPROF,
1009 1035 show=show)
1010 1036
1011 1037 nrow, ncol = self.getSubplots()
1012 1038
1013 1039 counter = 0
1014 1040 for y in range(nrow):
1015 1041 if counter >= self.nplots:
1016 1042 break
1017 1043
1018 1044 self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1)
1019 1045 counter += 1
1020 1046
1021 1047 def run(self, dataOut, id, wintitle="", channelList=None,
1022 1048 xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None,
1023 1049 zmaxVertical = None, zminVertical = None, zmaxZonal = None, zminZonal = None,
1024 1050 timerange=None, SNRthresh = -numpy.inf, SNRmin = None, SNRmax = None, SNR_1 = False,
1025 1051 save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True,
1026 1052 server=None, folder=None, username=None, password=None,
1027 1053 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1028 1054 """
1029 1055
1030 1056 Input:
1031 1057 dataOut :
1032 1058 id :
1033 1059 wintitle :
1034 1060 channelList :
1035 1061 showProfile :
1036 1062 xmin : None,
1037 1063 xmax : None,
1038 1064 ymin : None,
1039 1065 ymax : None,
1040 1066 zmin : None,
1041 1067 zmax : None
1042 1068 """
1043 1069
1044 1070 if timerange != None:
1045 1071 self.timerange = timerange
1046 1072
1047 1073 tmin = None
1048 1074 tmax = None
1049 1075
1050 1076 x = dataOut.getTimeRange1()
1051 1077 # y = dataOut.heightList
1052 1078 y = dataOut.heightList
1053 1079
1054 1080 z = dataOut.data_output
1055 1081 nplots = z.shape[0] #Number of wind dimensions estimated
1056 1082 nplotsw = nplots
1057 1083
1058 1084 #If there is a SNR function defined
1059 1085 if dataOut.data_SNR != None:
1060 1086 nplots += 1
1061 1087 SNR = dataOut.data_SNR
1062 1088
1063 1089 if SNR_1:
1064 1090 SNR += 1
1065 1091
1066 1092 SNRavg = numpy.average(SNR, axis=0)
1067 1093
1068 1094 SNRdB = 10*numpy.log10(SNR)
1069 1095 SNRavgdB = 10*numpy.log10(SNRavg)
1070 1096
1071 1097 ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0]
1072 1098
1073 1099 for i in range(nplotsw):
1074 1100 z[i,ind] = numpy.nan
1075 1101
1076 1102
1077 1103 showprofile = False
1078 1104 # thisDatetime = dataOut.datatime
1079 1105 thisDatetime = datetime.datetime.utcfromtimestamp(x[1])
1080 1106 title = wintitle + " EW Drifts"
1081 1107 xlabel = ""
1082 1108 ylabel = "Height (Km)"
1083 1109
1084 1110 if not self.__isConfig:
1085 1111
1086 1112 self.setup(id=id,
1087 1113 nplots=nplots,
1088 1114 wintitle=wintitle,
1089 1115 showprofile=showprofile,
1090 1116 show=show)
1091 1117
1092 1118 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1093 1119
1094 1120 if ymin == None: ymin = numpy.nanmin(y)
1095 1121 if ymax == None: ymax = numpy.nanmax(y)
1096 1122
1097 1123 if zmaxZonal == None: zmaxZonal = numpy.nanmax(abs(z[0,:]))
1098 1124 if zminZonal == None: zminZonal = -zmaxZonal
1099 1125 if zmaxVertical == None: zmaxVertical = numpy.nanmax(abs(z[1,:]))
1100 1126 if zminVertical == None: zminVertical = -zmaxVertical
1101 1127
1102 1128 if dataOut.data_SNR != None:
1103 1129 if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB)
1104 1130 if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB)
1105 1131
1106 1132 self.FTP_WEI = ftp_wei
1107 1133 self.EXP_CODE = exp_code
1108 1134 self.SUB_EXP_CODE = sub_exp_code
1109 1135 self.PLOT_POS = plot_pos
1110 1136
1111 1137 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1112 1138 self.__isConfig = True
1113 1139
1114 1140
1115 1141 self.setWinTitle(title)
1116 1142
1117 1143 if ((self.xmax - x[1]) < (x[1]-x[0])):
1118 1144 x[1] = self.xmax
1119 1145
1120 1146 strWind = ['Zonal','Vertical']
1121 1147 strCb = 'Velocity (m/s)'
1122 1148 zmaxVector = [zmaxZonal, zmaxVertical]
1123 1149 zminVector = [zminZonal, zminVertical]
1124 1150
1125 1151 for i in range(nplotsw):
1126 1152
1127 1153 title = "%s Drifts: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1128 1154 axes = self.axesList[i*self.__nsubplots]
1129 1155
1130 1156 z1 = z[i,:].reshape((1,-1))
1131 1157
1132 1158 axes.pcolorbuffer(x, y, z1,
1133 1159 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i],
1134 1160 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1135 1161 ticksize=9, cblabel=strCb, cbsize="1%", colormap="RdBu_r")
1136 1162
1137 1163 if dataOut.data_SNR != None:
1138 1164 i += 1
1139 1165 if SNR_1:
1140 1166 title = "Signal Noise Ratio + 1 (SNR+1): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1141 1167 else:
1142 1168 title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1143 1169 axes = self.axesList[i*self.__nsubplots]
1144 1170 SNRavgdB = SNRavgdB.reshape((1,-1))
1145 1171
1146 1172 axes.pcolorbuffer(x, y, SNRavgdB,
1147 1173 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax,
1148 1174 xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True,
1149 1175 ticksize=9, cblabel='', cbsize="1%", colormap="jet")
1150 1176
1151 1177 self.draw()
1152 1178
1153 1179 if x[1] >= self.axesList[0].xmax:
1154 1180 self.counter_imagwr = wr_period
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,
1161 1359 ftp=ftp,
1162 1360 wr_period=wr_period,
1163 1361 thisDatetime=thisDatetime,
1164 update_figfile=False) No newline at end of file
1362 update_figfile=False)
@@ -1,435 +1,437
1 1 import numpy
2 2 import datetime
3 3 import sys
4 4 import matplotlib
5 5
6 6 if 'linux' in sys.platform:
7 7 matplotlib.use("TKAgg")
8 8
9 9 if 'darwin' in sys.platform:
10 10 matplotlib.use("WXAgg")
11 11 #Qt4Agg', 'GTK', 'GTKAgg', 'ps', 'agg', 'cairo', 'MacOSX', 'GTKCairo', 'WXAgg', 'template', 'TkAgg', 'GTK3Cairo', 'GTK3Agg', 'svg', 'WebAgg', 'CocoaAgg', 'emf', 'gdk', 'WX'
12 12 import matplotlib.pyplot
13 13
14 14 from mpl_toolkits.axes_grid1 import make_axes_locatable
15 15 from matplotlib.ticker import *
16 16
17 17 ###########################################
18 18 #Actualizacion de las funciones del driver
19 19 ###########################################
20 20
21 21 def createFigure(id, wintitle, width, height, facecolor="w", show=True):
22 22
23 23 matplotlib.pyplot.ioff()
24 24 fig = matplotlib.pyplot.figure(num=id, facecolor=facecolor)
25 25 fig.canvas.manager.set_window_title(wintitle)
26 26 fig.canvas.manager.resize(width, height)
27 27 matplotlib.pyplot.ion()
28 28 if show:
29 29 matplotlib.pyplot.show()
30 30
31 31 return fig
32 32
33 33 def closeFigure(show=False, fig=None):
34 34
35 35 matplotlib.pyplot.ioff()
36 36 matplotlib.pyplot.pause(0.1)
37 37
38 38 if show:
39 39 matplotlib.pyplot.show()
40 40
41 41 if fig != None:
42 42 matplotlib.pyplot.close(fig)
43 43 matplotlib.pyplot.pause(0.1)
44 44 matplotlib.pyplot.ion()
45 45 return
46 46
47 47 matplotlib.pyplot.close("all")
48 48 matplotlib.pyplot.pause(0.1)
49 49 matplotlib.pyplot.ion()
50 50 return
51 51
52 52 def saveFigure(fig, filename):
53 53
54 54 matplotlib.pyplot.ioff()
55 55 fig.savefig(filename)
56 56 matplotlib.pyplot.ion()
57 57
58 58 def setWinTitle(fig, title):
59 59
60 60 fig.canvas.manager.set_window_title(title)
61 61
62 62 def setTitle(fig, title):
63 63
64 64 fig.suptitle(title)
65 65
66 66 def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan, polar=False):
67 67
68 68 matplotlib.pyplot.ioff()
69 69 matplotlib.pyplot.figure(fig.number)
70 70 axes = matplotlib.pyplot.subplot2grid((nrow, ncol),
71 71 (xpos, ypos),
72 72 colspan=colspan,
73 73 rowspan=rowspan,
74 74 polar=polar)
75 75
76 76 matplotlib.pyplot.ion()
77 77 return axes
78 78
79 79 def setAxesText(ax, text):
80 80
81 81 ax.annotate(text,
82 82 xy = (.1, .99),
83 83 xycoords = 'figure fraction',
84 84 horizontalalignment = 'left',
85 85 verticalalignment = 'top',
86 86 fontsize = 10)
87 87
88 88 def printLabels(ax, xlabel, ylabel, title):
89 89
90 90 ax.set_xlabel(xlabel, size=11)
91 91 ax.set_ylabel(ylabel, size=11)
92 92 ax.set_title(title, size=8)
93 93
94 94 def createPline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='',
95 95 ticksize=9, xtick_visible=True, ytick_visible=True,
96 96 nxticks=4, nyticks=10,
97 97 grid=None,color='blue'):
98 98
99 99 """
100 100
101 101 Input:
102 102 grid : None, 'both', 'x', 'y'
103 103 """
104 104
105 105 matplotlib.pyplot.ioff()
106 106
107 107 ax.set_xlim([xmin,xmax])
108 108 ax.set_ylim([ymin,ymax])
109 109
110 110 printLabels(ax, xlabel, ylabel, title)
111 111
112 112 ######################################################
113 113 if (xmax-xmin)<=1:
114 114 xtickspos = numpy.linspace(xmin,xmax,nxticks)
115 115 xtickspos = numpy.array([float("%.1f"%i) for i in xtickspos])
116 116 ax.set_xticks(xtickspos)
117 117 else:
118 118 xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
119 119 # xtickspos = numpy.arange(nxticks)*float(xmax-xmin)/float(nxticks) + int(xmin)
120 120 ax.set_xticks(xtickspos)
121 121
122 122 for tick in ax.get_xticklabels():
123 123 tick.set_visible(xtick_visible)
124 124
125 125 for tick in ax.xaxis.get_major_ticks():
126 126 tick.label.set_fontsize(ticksize)
127 127
128 128 ######################################################
129 129 for tick in ax.get_yticklabels():
130 130 tick.set_visible(ytick_visible)
131 131
132 132 for tick in ax.yaxis.get_major_ticks():
133 133 tick.label.set_fontsize(ticksize)
134 134
135 135 ax.plot(x, y, color=color)
136 136 iplot = ax.lines[-1]
137 137
138 138 ######################################################
139 139 if '0.' in matplotlib.__version__[0:2]:
140 140 print "The matplotlib version has to be updated to 1.1 or newer"
141 141 return iplot
142 142
143 143 if '1.0.' in matplotlib.__version__[0:4]:
144 144 print "The matplotlib version has to be updated to 1.1 or newer"
145 145 return iplot
146 146
147 147 if grid != None:
148 148 ax.grid(b=True, which='major', axis=grid)
149 149
150 150 matplotlib.pyplot.tight_layout()
151 151
152 152 matplotlib.pyplot.ion()
153 153
154 154 return iplot
155 155
156 156 def set_linedata(ax, x, y, idline):
157 157
158 158 ax.lines[idline].set_data(x,y)
159 159
160 160 def pline(iplot, x, y, xlabel='', ylabel='', title=''):
161 161
162 162 ax = iplot.get_axes()
163 163
164 164 printLabels(ax, xlabel, ylabel, title)
165 165
166 166 set_linedata(ax, x, y, idline=0)
167 167
168 168 def addpline(ax, x, y, color, linestyle, lw):
169 169
170 170 ax.plot(x,y,color=color,linestyle=linestyle,lw=lw)
171 171
172 172
173 173 def createPcolor(ax, x, y, z, xmin, xmax, ymin, ymax, zmin, zmax,
174 174 xlabel='', ylabel='', title='', ticksize = 9,
175 175 colormap='jet',cblabel='', cbsize="5%",
176 176 XAxisAsTime=False):
177 177
178 178 matplotlib.pyplot.ioff()
179 179
180 180 divider = make_axes_locatable(ax)
181 181 ax_cb = divider.new_horizontal(size=cbsize, pad=0.05)
182 182 fig = ax.get_figure()
183 183 fig.add_axes(ax_cb)
184 184
185 185 ax.set_xlim([xmin,xmax])
186 186 ax.set_ylim([ymin,ymax])
187 187
188 188 printLabels(ax, xlabel, ylabel, title)
189 189
190 190 imesh = ax.pcolormesh(x,y,z.T, vmin=zmin, vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
191 191 cb = matplotlib.pyplot.colorbar(imesh, cax=ax_cb)
192 192 cb.set_label(cblabel)
193 193
194 194 # for tl in ax_cb.get_yticklabels():
195 195 # tl.set_visible(True)
196 196
197 197 for tick in ax.yaxis.get_major_ticks():
198 198 tick.label.set_fontsize(ticksize)
199 199
200 200 for tick in ax.xaxis.get_major_ticks():
201 201 tick.label.set_fontsize(ticksize)
202 202
203 203 for tick in cb.ax.get_yticklabels():
204 204 tick.set_fontsize(ticksize)
205 205
206 206 ax_cb.yaxis.tick_right()
207 207
208 208 if '0.' in matplotlib.__version__[0:2]:
209 209 print "The matplotlib version has to be updated to 1.1 or newer"
210 210 return imesh
211 211
212 212 if '1.0.' in matplotlib.__version__[0:4]:
213 213 print "The matplotlib version has to be updated to 1.1 or newer"
214 214 return imesh
215 215
216 216 matplotlib.pyplot.tight_layout()
217 217
218 218 if XAxisAsTime:
219 219
220 220 func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
221 221 ax.xaxis.set_major_formatter(FuncFormatter(func))
222 222 ax.xaxis.set_major_locator(LinearLocator(7))
223 223
224 224 matplotlib.pyplot.ion()
225 225 return imesh
226 226
227 227 def pcolor(imesh, z, xlabel='', ylabel='', title=''):
228 228
229 229 z = z.T
230 230 ax = imesh.get_axes()
231 231 printLabels(ax, xlabel, ylabel, title)
232 232 imesh.set_array(z.ravel())
233 233
234 234 def addpcolor(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
235 235
236 236 printLabels(ax, xlabel, ylabel, title)
237 237
238 238 ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
239 239
240 240 def addpcolorbuffer(ax, x, y, z, zmin, zmax, xlabel='', ylabel='', title='', colormap='jet'):
241 241
242 242 printLabels(ax, xlabel, ylabel, title)
243 243
244 244 ax.collections.remove(ax.collections[0])
245 245
246 246 ax.pcolormesh(x,y,z.T,vmin=zmin,vmax=zmax, cmap=matplotlib.pyplot.get_cmap(colormap))
247 247
248 248 def createPmultiline(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
249 249 ticksize=9, xtick_visible=True, ytick_visible=True,
250 250 nxticks=4, nyticks=10,
251 251 grid=None):
252 252
253 253 """
254 254
255 255 Input:
256 256 grid : None, 'both', 'x', 'y'
257 257 """
258 258
259 259 matplotlib.pyplot.ioff()
260 260
261 261 lines = ax.plot(x.T, y)
262 262 leg = ax.legend(lines, legendlabels, loc='upper right')
263 263 leg.get_frame().set_alpha(0.5)
264 264 ax.set_xlim([xmin,xmax])
265 265 ax.set_ylim([ymin,ymax])
266 266 printLabels(ax, xlabel, ylabel, title)
267 267
268 268 xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
269 269 ax.set_xticks(xtickspos)
270 270
271 271 for tick in ax.get_xticklabels():
272 272 tick.set_visible(xtick_visible)
273 273
274 274 for tick in ax.xaxis.get_major_ticks():
275 275 tick.label.set_fontsize(ticksize)
276 276
277 277 for tick in ax.get_yticklabels():
278 278 tick.set_visible(ytick_visible)
279 279
280 280 for tick in ax.yaxis.get_major_ticks():
281 281 tick.label.set_fontsize(ticksize)
282 282
283 283 iplot = ax.lines[-1]
284 284
285 285 if '0.' in matplotlib.__version__[0:2]:
286 286 print "The matplotlib version has to be updated to 1.1 or newer"
287 287 return iplot
288 288
289 289 if '1.0.' in matplotlib.__version__[0:4]:
290 290 print "The matplotlib version has to be updated to 1.1 or newer"
291 291 return iplot
292 292
293 293 if grid != None:
294 294 ax.grid(b=True, which='major', axis=grid)
295 295
296 296 matplotlib.pyplot.tight_layout()
297 297
298 298 matplotlib.pyplot.ion()
299 299
300 300 return iplot
301 301
302 302
303 303 def pmultiline(iplot, x, y, xlabel='', ylabel='', title=''):
304 304
305 305 ax = iplot.get_axes()
306 306
307 307 printLabels(ax, xlabel, ylabel, title)
308 308
309 309 for i in range(len(ax.lines)):
310 310 line = ax.lines[i]
311 311 line.set_data(x[i,:],y)
312 312
313 313 def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='', title='', legendlabels=None,
314 314 ticksize=9, xtick_visible=True, ytick_visible=True,
315 315 nxticks=4, nyticks=10, marker='.', markersize=10, linestyle="None",
316 316 grid=None, XAxisAsTime=False):
317 317
318 318 """
319 319
320 320 Input:
321 321 grid : None, 'both', 'x', 'y'
322 322 """
323 323
324 324 matplotlib.pyplot.ioff()
325 325
326 326 # lines = ax.plot(x, y.T, marker=marker,markersize=markersize,linestyle=linestyle)
327 327 lines = ax.plot(x, y.T, linestyle=linestyle, marker=marker, markersize=markersize)
328 328 leg = ax.legend(lines, legendlabels, loc='upper left', bbox_to_anchor=(1.01, 1.00), numpoints=1, handlelength=1.5, \
329 329 handletextpad=0.5, borderpad=0.5, labelspacing=0.5, borderaxespad=0.)
330 330
331 331 for label in leg.get_texts(): label.set_fontsize(9)
332 332
333 333 ax.set_xlim([xmin,xmax])
334 334 ax.set_ylim([ymin,ymax])
335 335 printLabels(ax, xlabel, ylabel, title)
336 336
337 337 # xtickspos = numpy.arange(nxticks)*int((xmax-xmin)/(nxticks)) + int(xmin)
338 338 # ax.set_xticks(xtickspos)
339 339
340 340 for tick in ax.get_xticklabels():
341 341 tick.set_visible(xtick_visible)
342 342
343 343 for tick in ax.xaxis.get_major_ticks():
344 344 tick.label.set_fontsize(ticksize)
345 345
346 346 for tick in ax.get_yticklabels():
347 347 tick.set_visible(ytick_visible)
348 348
349 349 for tick in ax.yaxis.get_major_ticks():
350 350 tick.label.set_fontsize(ticksize)
351 351
352 352 iplot = ax.lines[-1]
353 353
354 354 if '0.' in matplotlib.__version__[0:2]:
355 355 print "The matplotlib version has to be updated to 1.1 or newer"
356 356 return iplot
357 357
358 358 if '1.0.' in matplotlib.__version__[0:4]:
359 359 print "The matplotlib version has to be updated to 1.1 or newer"
360 360 return iplot
361 361
362 362 if grid != None:
363 363 ax.grid(b=True, which='major', axis=grid)
364 364
365 365 matplotlib.pyplot.tight_layout()
366 366
367 367 if XAxisAsTime:
368 368
369 369 func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime("%H:%M:%S"))
370 370 ax.xaxis.set_major_formatter(FuncFormatter(func))
371 371 ax.xaxis.set_major_locator(LinearLocator(7))
372 372
373 373 matplotlib.pyplot.ion()
374 374
375 375 return iplot
376 376
377 377 def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''):
378 378
379 379 ax = iplot.get_axes()
380 380
381 381 printLabels(ax, xlabel, ylabel, title)
382 382
383 383 for i in range(len(ax.lines)):
384 384 line = ax.lines[i]
385 385 line.set_data(x,y[i,:])
386 386
387 387 def createPolar(ax, x, y,
388 388 xlabel='', ylabel='', title='', ticksize = 9,
389 389 colormap='jet',cblabel='', cbsize="5%",
390 390 XAxisAsTime=False):
391 391
392 392 matplotlib.pyplot.ioff()
393 393
394 394 ax.plot(x,y,'bo', markersize=5)
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]:
404 406 print "The matplotlib version has to be updated to 1.1 or newer"
405 407 return iplot
406 408
407 409 if '1.0.' in matplotlib.__version__[0:4]:
408 410 print "The matplotlib version has to be updated to 1.1 or newer"
409 411 return iplot
410 412
411 413 # if grid != None:
412 414 # ax.grid(b=True, which='major', axis=grid)
413 415
414 416 matplotlib.pyplot.tight_layout()
415 417
416 418 matplotlib.pyplot.ion()
417 419
418 420
419 421 return iplot
420 422
421 423 def polar(iplot, x, y, xlabel='', ylabel='', title=''):
422 424
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, '', title)
428 printLabels(ax, xlabel, ylabel, title)
427 429
428 430 set_linedata(ax, x, y, idline=0)
429 431
430 432 def draw(fig):
431 433
432 434 if type(fig) == 'int':
433 435 raise ValueError, "Error drawing: Fig parameter should be a matplotlib figure object figure"
434 436
435 437 fig.canvas.draw()
@@ -1,27 +1,28
1 1 '''
2 2 @author: roj-idl71
3 3 '''
4 4 #USED IN jroplot_spectra.py
5 5 RTI_CODE = 0 #Range time intensity (RTI).
6 6 SPEC_CODE = 1 #Spectra (and Cross-spectra) information.
7 7 CROSS_CODE = 2 #Cross-Correlation information.
8 8 COH_CODE = 3 #Coherence map.
9 9 BASE_CODE = 4 #Base lines graphic.
10 10 ROW_CODE = 5 #Row Spectra.
11 11 TOTAL_CODE = 6 #Total Power.
12 12 DRIFT_CODE = 7 #Drifts graphics.
13 13 HEIGHT_CODE = 8 #Height profile.
14 14 PHASE_CODE = 9 #Signal Phase.
15 15
16 16 POWER_CODE = 16
17 17 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 PARMS_CODE = 23
26 SPECFIT_CODE = 24
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
@@ -1,903 +1,1005
1 1 import numpy
2 2 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
9 10 from schainpy.model.io.jroIO_base import *
10 11
11 12
12 13 class HDF5Reader(ProcessingUnit):
13 14
14 15 ext = ".hdf5"
15 16
16 17 optchar = "D"
17 18
18 19 timezone = None
19 20
20 21 secStart = None
21 22
22 23 secEnd = None
23 24
24 25 fileIndex = None
25 26
26 27 blockIndex = None
27 28
28 29 blocksPerFile = None
29 30
30 31 path = None
31 32
32 33 #List of Files
33 34
34 35 filenameList = None
35 36
36 37 datetimeList = None
37 38
38 39 #Hdf5 File
39 40
40 41 fpMetadata = None
41 42
42 43 pathMeta = None
43 44
44 45 listMetaname = None
45 46
46 47 listMeta = None
47 48
48 49 listDataname = None
49 50
50 51 listData = None
51 52
52 53 listShapes = None
53 54
54 55 fp = None
55 56
56 57 #dataOut reconstruction
57 58
58 59 dataOut = None
59 60
60 61 nRecords = None
61 62
62 63
63 64 def __init__(self):
64 65 self.dataOut = self.__createObjByDefault()
65 66 return
66 67
67 68 def __createObjByDefault(self):
68 69
69 70 dataObj = Parameters()
70 71
71 72 return dataObj
72 73
73 74 def setup(self,path=None,
74 75 startDate=None,
75 76 endDate=None,
76 77 startTime=datetime.time(0,0,0),
77 78 endTime=datetime.time(23,59,59),
78 79 walk=True,
79 80 timezone='ut',
80 81 all=0,
81 82 online=False,
82 83 ext=None):
83 84
84 85 if ext==None:
85 86 ext = self.ext
86 87 self.timezone = timezone
87 88 # self.all = all
88 89 # self.online = online
89 90 self.path = path
90 91
91 92 startDateTime = datetime.datetime.combine(startDate,startTime)
92 93 endDateTime = datetime.datetime.combine(endDate,endTime)
93 94 secStart = (startDateTime-datetime.datetime(1970,1,1)).total_seconds()
94 95 secEnd = (endDateTime-datetime.datetime(1970,1,1)).total_seconds()
95 96
96 97 self.secStart = secStart
97 98 self.secEnd = secEnd
98 99
99 100 if not(online):
100 101 #Busqueda de archivos offline
101 102 self.__searchFilesOffline(path, startDate, endDate, ext, startTime, endTime, secStart, secEnd, walk)
102 103 else:
103 104 self.__searchFilesOnline(path, walk)
104 105
105 106 if not(self.filenameList):
106 107 print "There is no files into the folder: %s"%(path)
107 108 sys.exit(-1)
108 109
109 110 # self.__getExpParameters()
110 111
111 112 self.fileIndex = -1
112 113
113 114 self.__setNextFileOffline()
114 115
115 116 self.__readMetadata()
116 117
117 118 self.blockIndex = 0
118 119
119 120 return
120 121
121 122 def __searchFilesOffline(self,
122 123 path,
123 124 startDate,
124 125 endDate,
125 126 ext,
126 127 startTime=datetime.time(0,0,0),
127 128 endTime=datetime.time(23,59,59),
128 129 secStart = 0,
129 130 secEnd = numpy.inf,
130 131 walk=True):
131 132
132 133 # self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
133 134 #
134 135 # self.__checkPath()
135 136 #
136 137 # self.__findDataForDates()
137 138 #
138 139 # self.__selectDataForTimes()
139 140 #
140 141 # for i in range(len(self.filenameList)):
141 142 # print "%s" %(self.filenameList[i])
142 143
143 144 pathList = []
144 145
145 146 if not walk:
146 147 #pathList.append(path)
147 148 multi_path = path.split(',')
148 149 for single_path in multi_path:
149 150 pathList.append(single_path)
150 151
151 152 else:
152 153 #dirList = []
153 154 multi_path = path.split(',')
154 155 for single_path in multi_path:
155 156 dirList = []
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 isRadarFolder(thisPath):
160 if not isDoyFolder(thisPath):
160 161 continue
161 162
162 163 dirList.append(thisPath)
163 164
164 165 if not(dirList):
165 166 return None, None
166 167
167 168 thisDate = startDate
168 169
169 170 while(thisDate <= endDate):
170 171 year = thisDate.timetuple().tm_year
171 172 doy = thisDate.timetuple().tm_yday
172 173
173 174 matchlist = fnmatch.filter(dirList, '?' + '%4.4d%3.3d' % (year,doy) + '*')
174 175 if len(matchlist) == 0:
175 176 thisDate += datetime.timedelta(1)
176 177 continue
177 178 for match in matchlist:
178 179 pathList.append(os.path.join(single_path,match))
179 180
180 181 thisDate += datetime.timedelta(1)
181 182
182 183 if pathList == []:
183 184 print "Any folder was found for the date range: %s-%s" %(startDate, endDate)
184 185 return None, None
185 186
186 187 print "%d folder(s) was(were) found for the date range: %s - %s" %(len(pathList), startDate, endDate)
187 188
188 189 filenameList = []
189 190 datetimeList = []
190 191 pathDict = {}
191 192 filenameList_to_sort = []
192 193
193 194 for i in range(len(pathList)):
194 195
195 196 thisPath = pathList[i]
196 197
197 198 fileList = glob.glob1(thisPath, "*%s" %ext)
198 199 fileList.sort()
199 200 pathDict.setdefault(fileList[0])
200 201 pathDict[fileList[0]] = i
201 202 filenameList_to_sort.append(fileList[0])
202 203
203 204 filenameList_to_sort.sort()
204 205
205 206 for file in filenameList_to_sort:
206 207 thisPath = pathList[pathDict[file]]
207 208
208 209 fileList = glob.glob1(thisPath, "*%s" %ext)
209 210 fileList.sort()
210 211
211 212 for file in fileList:
212 213
213 214 filename = os.path.join(thisPath,file)
214 215 thisDatetime = self.__isFileinThisTime(filename, secStart, secEnd)
215 216
216 217 if not(thisDatetime):
217 218 continue
218 219
219 220 filenameList.append(filename)
220 221 datetimeList.append(thisDatetime)
221 222
222 223 if not(filenameList):
223 224 print "Any file was found for the time range %s - %s" %(startTime, endTime)
224 225 return None, None
225 226
226 227 print "%d file(s) was(were) found for the time range: %s - %s" %(len(filenameList), startTime, endTime)
227 228 print
228 229
229 230 for i in range(len(filenameList)):
230 231 print "%s -> [%s]" %(filenameList[i], datetimeList[i].ctime())
231 232
232 233 self.filenameList = filenameList
233 234 self.datetimeList = datetimeList
234 235
235 236 return pathList, filenameList
236 237
237 238 def __isFileinThisTime(self, filename, startSeconds, endSeconds):
238 239 """
239 240 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
240 241
241 242 Inputs:
242 243 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
243 244
244 245 startTime : tiempo inicial del rango seleccionado en formato datetime.time
245 246
246 247 endTime : tiempo final del rango seleccionado en formato datetime.time
247 248
248 249 Return:
249 250 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
250 251 fecha especificado, de lo contrario retorna False.
251 252
252 253 Excepciones:
253 254 Si el archivo no existe o no puede ser abierto
254 255 Si la cabecera no puede ser leida.
255 256
256 257 """
257 258
258 259 try:
259 260 fp = fp = h5py.File(filename,'r')
260 261 except IOError:
261 262 traceback.print_exc()
262 263 raise IOError, "The file %s can't be opened" %(filename)
263 264
264 265 grp = fp['Data']
265 266 timeAux = grp['time']
266 267 time0 = timeAux[:][0].astype(numpy.float) #Time Vector
267 268
268 269 fp.close()
269 270
270 271 if self.timezone == 'lt':
271 272 time0 -= 5*3600
272 273
273 274 boolTimer = numpy.logical_and(time0 >= startSeconds,time0 < endSeconds)
274 275
275 276 if not (numpy.any(boolTimer)):
276 277 return None
277 278
278 279 thisDatetime = datetime.datetime.utcfromtimestamp(time0[0])
279 280 return thisDatetime
280 281
281 282 def __checkPath(self):
282 283 if os.path.exists(self.path):
283 284 self.status = 1
284 285 else:
285 286 self.status = 0
286 287 print 'Path:%s does not exists'%self.path
287 288
288 289 return
289 290
290 291 def __setNextFileOffline(self):
291 292 idFile = self.fileIndex
292 293 idFile += 1
293 294
294 295 if not(idFile < len(self.filenameList)):
295 296 print "No more Files"
296 297 return 0
297 298
298 299 filename = self.filenameList[idFile]
299 300
300 301 filePointer = h5py.File(filename,'r')
301 302
302 303 self.flagIsNewFile = 1
303 304 self.fileIndex = idFile
304 305 self.filename = filename
305 306
306 307 self.fp = filePointer
307 308
308 309 print "Setting the file: %s"%self.filename
309 310
310 311 self.__readMetadata()
311 312 self.__setBlockList()
312 313 # self.nRecords = self.fp['Data'].attrs['blocksPerFile']
313 314 self.nRecords = self.fp['Data'].attrs['nRecords']
314 315 self.blockIndex = 0
315 316 return 1
316 317
317 318 def __setBlockList(self):
318 319 '''
319 320 self.fp
320 321 self.startDateTime
321 322 self.endDateTime
322 323
323 324 self.blockList
324 325 self.blocksPerFile
325 326
326 327 '''
327 328 filePointer = self.fp
328 329 secStart = self.secStart
329 330 secEnd = self.secEnd
330 331
331 332 grp = filePointer['Data']
332 333 timeVector = grp['time'].value.astype(numpy.float)[0]
333 334
334 335 if self.timezone == 'lt':
335 336 timeVector -= 5*3600
336 337
337 338 ind = numpy.where(numpy.logical_and(timeVector >= secStart , timeVector < secEnd))[0]
338 339
339 340 self.blockList = ind
340 341 self.blocksPerFile = len(ind)
341 342
342 343 return
343 344
344 345 def __readMetadata(self):
345 346 '''
346 347 self.pathMeta
347 348
348 349 self.listShapes
349 350 self.listMetaname
350 351 self.listMeta
351 352
352 353 '''
353 354
354 355 grp = self.fp['Data']
355 356 pathMeta = os.path.join(self.path, grp.attrs['metadata'])
356 357
357 358 if pathMeta == self.pathMeta:
358 359 return
359 360 else:
360 361 self.pathMeta = pathMeta
361 362
362 363 filePointer = h5py.File(self.pathMeta,'r')
363 364 groupPointer = filePointer['Metadata']
364 365
365 366 listMetaname = []
366 367 listMetadata = []
367 368 for item in groupPointer.items():
368 369 name = item[0]
369 370
370 371 if name=='array dimensions':
371 372 table = groupPointer[name][:]
372 373 listShapes = {}
373 374 for shapes in table:
374 375 listShapes[shapes[0]] = numpy.array([shapes[1],shapes[2],shapes[3],shapes[4]])
375 376 else:
376 377 data = groupPointer[name].value
377 378 listMetaname.append(name)
378 379 listMetadata.append(data)
379 380
380 381 if name=='type':
381 382 self.__initDataOut(data)
382 383
383 384 filePointer.close()
384 385
385 386 self.listShapes = listShapes
386 387 self.listMetaname = listMetaname
387 388 self.listMeta = listMetadata
388 389
389 390 return
390 391
391 392 def __readData(self):
392 393 grp = self.fp['Data']
393 394 listdataname = []
394 395 listdata = []
395 396
396 397 for item in grp.items():
397 398 name = item[0]
398 399
399 400 if name == 'time':
400 401 listdataname.append('utctime')
401 402 timeAux = grp[name].value.astype(numpy.float)[0]
402 403 listdata.append(timeAux)
403 404 continue
404 405
405 406 listdataname.append(name)
406 407 array = self.__setDataArray(self.nRecords, grp[name],self.listShapes[name])
407 408 listdata.append(array)
408 409
409 410 self.listDataname = listdataname
410 411 self.listData = listdata
411 412 return
412 413
413 414 def __setDataArray(self, nRecords, dataset, shapes):
414 415
415 416 nChannels = shapes[0] #Dimension 0
416 417
417 418 nPoints = shapes[1] #Dimension 1, number of Points or Parameters
418 419
419 420 nSamples = shapes[2] #Dimension 2, number of samples or ranges
420 421
421 422 mode = shapes[3]
422 423
423 424 # if nPoints>1:
424 425 # arrayData = numpy.zeros((nRecords,nChannels,nPoints,nSamples))
425 426 # else:
426 427 # arrayData = numpy.zeros((nRecords,nChannels,nSamples))
427 428 #
428 429 # chn = 'channel'
429 430 #
430 431 # for i in range(nChannels):
431 432 #
432 433 # data = dataset[chn + str(i)].value
433 434 #
434 435 # if nPoints>1:
435 436 # data = numpy.rollaxis(data,2)
436 437 #
437 438 # arrayData[:,i,:] = data
438 439
439 440 arrayData = numpy.zeros((nRecords,nChannels,nPoints,nSamples))
440 441 doSqueeze = False
441 442 if mode == 0:
442 443 strds = 'channel'
443 444 nDatas = nChannels
444 445 newShapes = (nRecords,nPoints,nSamples)
445 446 if nPoints == 1:
446 447 doSqueeze = True
447 448 axisSqueeze = 2
448 449 else:
449 450 strds = 'param'
450 451 nDatas = nPoints
451 452 newShapes = (nRecords,nChannels,nSamples)
452 453 if nChannels == 1:
453 454 doSqueeze = True
454 455 axisSqueeze = 1
455 456
456 457 for i in range(nDatas):
457 458
458 459 data = dataset[strds + str(i)].value
459 460 data = data.reshape(newShapes)
460 461
461 462 if mode == 0:
462 463 arrayData[:,i,:,:] = data
463 464 else:
464 465 arrayData[:,:,i,:] = data
465 466
466 467 if doSqueeze:
467 468 arrayData = numpy.squeeze(arrayData, axis=axisSqueeze)
468 469
469 470 return arrayData
470 471
471 472 def __initDataOut(self, type):
472 473
473 474 # if type =='Parameters':
474 475 # self.dataOut = Parameters()
475 476 # elif type =='Spectra':
476 477 # self.dataOut = Spectra()
477 478 # elif type =='Voltage':
478 479 # self.dataOut = Voltage()
479 480 # elif type =='Correlation':
480 481 # self.dataOut = Correlation()
481 482
482 483 return
483 484
484 485 def __setDataOut(self):
485 486 listMeta = self.listMeta
486 487 listMetaname = self.listMetaname
487 488 listDataname = self.listDataname
488 489 listData = self.listData
489 490
490 491 blockIndex = self.blockIndex
491 492 blockList = self.blockList
492 493
493 494 for i in range(len(listMeta)):
494 495 setattr(self.dataOut,listMetaname[i],listMeta[i])
495 496
496 497 for j in range(len(listData)):
497 498 if listDataname[j]=='utctime':
498 499 # setattr(self.dataOut,listDataname[j],listData[j][blockList[blockIndex]])
499 500 setattr(self.dataOut,'utctimeInit',listData[j][blockList[blockIndex]])
500 501 continue
501 502
502 503 setattr(self.dataOut,listDataname[j],listData[j][blockList[blockIndex],:])
503 504
504 505 return self.dataOut.data_param
505 506
506 507 def getData(self):
507 508
508 509 # if self.flagNoMoreFiles:
509 510 # self.dataOut.flagNoData = True
510 511 # print 'Process finished'
511 512 # return 0
512 513 #
513 514 if self.blockIndex==self.blocksPerFile:
514 515 if not( self.__setNextFileOffline() ):
515 516 self.dataOut.flagNoData = True
516 517 return 0
517 518
518 519 #
519 520 # if self.datablock == None: # setear esta condicion cuando no hayan datos por leers
520 521 # self.dataOut.flagNoData = True
521 522 # return 0
522 523
523 524 self.__readData()
524 525 self.__setDataOut()
525 526 self.dataOut.flagNoData = False
526 527
527 528 self.blockIndex += 1
528 529
529 530 return
530 531
531 532 def run(self, **kwargs):
532 533
533 534 if not(self.isConfig):
534 535 self.setup(**kwargs)
535 536 # self.setObjProperties()
536 537 self.isConfig = True
537 538
538 539 self.getData()
539 540
540 541 return
541 542
542 543 class HDF5Writer(Operation):
543 544
544 545 ext = ".hdf5"
545 546
546 547 optchar = "D"
547 548
548 549 metaoptchar = "M"
549 550
550 551 metaFile = None
551 552
553 filename = None
554
552 555 path = None
553 556
554 557 setFile = None
555 558
556 559 fp = None
557 560
558 561 grp = None
559 562
560 563 ds = None
561 564
562 565 firsttime = True
563 566
564 567 #Configurations
565 568
566 569 blocksPerFile = None
567 570
568 571 blockIndex = None
569 572
570 573 dataOut = None
571 574
572 575 #Data Arrays
573 576
574 577 dataList = None
575 578
576 579 metadataList = None
577 580
578 581 arrayDim = None
579 582
580 583 tableDim = None
581 584
582 585 # dtype = [('arrayName', 'S20'),('nChannels', 'i'), ('nPoints', 'i'), ('nSamples', 'i'),('mode', 'b')]
583 586
584 587 dtype = [('arrayName', 'S20'),('nDimensions', 'i'), ('dim2', 'i'), ('dim1', 'i'),('dim0', 'i'),('mode', 'b')]
585 588
586 589 mode = None
587 590
588 591 nDatas = None #Number of datasets to be stored per array
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)
595 600 self.isConfig = False
596 601 return
597 602
598 603
599 604 def setup(self, dataOut, **kwargs):
600 605
601 606 self.path = kwargs['path']
602 607
603 608 if kwargs.has_key('ext'):
604 609 self.ext = kwargs['ext']
605 610
606 611 if kwargs.has_key('blocksPerFile'):
607 612 self.blocksPerFile = kwargs['blocksPerFile']
608 613 else:
609 614 self.blocksPerFile = 10
610 615
611 616 self.metadataList = kwargs['metadataList']
612 617
613 618 self.dataList = kwargs['dataList']
614 619
615 620 self.dataOut = dataOut
616 621
617 622 if kwargs.has_key('mode'):
618 623 mode = kwargs['mode']
619 624
620 625 if type(mode) == int:
621 626 mode = numpy.zeros(len(self.dataList)) + mode
622 627 else:
623 628 mode = numpy.zeros(len(self.dataList))
624 629
625 630 self.mode = mode
626 631
627 632 arrayDim = numpy.zeros((len(self.dataList),5))
628 633
629 634 #Table dimensions
630 635
631 636 dtype0 = self.dtype
632 637
633 638 tableList = []
634 639
635 640 for i in range(len(self.dataList)):
636 641
637 642 dataAux = getattr(self.dataOut, self.dataList[i])
638 643
639 644 if type(dataAux)==float or type(dataAux)==int:
640 645 arrayDim[i,0] = 1
641 646 else:
642 647 arrayDim0 = dataAux.shape
643 648 arrayDim[i,0] = len(arrayDim0)
644 649 arrayDim[i,4] = mode[i]
645 650
646 651 if len(arrayDim0) == 3:
647 652 arrayDim[i,1:-1] = numpy.array(arrayDim0)
648 653 elif len(arrayDim0) == 2:
649 654 arrayDim[i,2:-1] = numpy.array(arrayDim0) #nHeights
650 655 elif len(arrayDim0) == 1:
651 656 arrayDim[i,3] = arrayDim0
652 657 elif len(arrayDim0) == 0:
653 658 arrayDim[i,0] = 1
654 659 arrayDim[i,3] = 1
655 660
656 661 table = numpy.array((self.dataList[i],) + tuple(arrayDim[i,:]),dtype = dtype0)
657 662 tableList.append(table)
658 663
659 664 self.arrayDim = arrayDim
660 665 self.tableDim = numpy.array(tableList, dtype = dtype0)
661 666 self.blockIndex = 0
662 667
663 668 return
664 669
665 670 def putMetadata(self):
666 671
667 672 fp = self.createMetadataFile()
668 673 self.writeMetadata(fp)
669 674 fp.close()
670 675 return
671 676
672 677 def createMetadataFile(self):
673 678 ext = self.ext
674 679 path = self.path
675 680 setFile = self.setFile
676 681
677 682 timeTuple = time.localtime(self.dataOut.utctime)
678 subfolder = ''
679
683
684 subfolder = ''
680 685 fullpath = os.path.join( path, subfolder )
686
681 687 if not( os.path.exists(fullpath) ):
682 688 os.mkdir(fullpath)
683 689 setFile = -1 #inicializo mi contador de seteo
690
691 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
692 fullpath = os.path.join( path, subfolder )
693
694 if not( os.path.exists(fullpath) ):
695 os.mkdir(fullpath)
696 setFile = -1 #inicializo mi contador de seteo
697
684 698 else:
685 699 filesList = os.listdir( fullpath )
700 filesList = sorted( filesList, key=str.lower )
686 701 if len( filesList ) > 0:
687 filesList = sorted( filesList, key=str.lower )
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)
691 706 # x YYYY DDD SSS .ext
692 707 if isNumber( filen[8:11] ):
693 708 setFile = int( filen[8:11] ) #inicializo mi contador de seteo al seteo del ultimo file
694 709 else:
695 710 setFile = -1
696 711 else:
697 712 setFile = -1 #inicializo mi contador de seteo
698 713
699 714 setFile += 1
700
715
701 716 file = '%s%4.4d%3.3d%3.3d%s' % (self.metaoptchar,
702 717 timeTuple.tm_year,
703 718 timeTuple.tm_yday,
704 719 setFile,
705 720 ext )
706 721
707 722 filename = os.path.join( path, subfolder, file )
708 723 self.metaFile = file
709 724 #Setting HDF5 File
710 725 fp = h5py.File(filename,'w')
711 726
712 727 return fp
713 728
714 729 def writeMetadata(self, fp):
715 730
716 731 grp = fp.create_group("Metadata")
717 732 grp.create_dataset('array dimensions', data = self.tableDim, dtype = self.dtype)
718 733
719 734 for i in range(len(self.metadataList)):
720 735 grp.create_dataset(self.metadataList[i], data=getattr(self.dataOut, self.metadataList[i]))
721 736 return
722 737
723 738 def setNextFile(self):
724 739
725 740 ext = self.ext
726 741 path = self.path
727 742 setFile = self.setFile
728 743 mode = self.mode
729
730 if self.fp != None:
731 self.fp.close()
732
744
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 os.mkdir(fullpath)
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]
745 756 # el filename debera tener el siguiente formato
746 757 # 0 1234 567 89A BCDE (hex)
747 758 # x YYYY DDD SSS .ext
748 759 if isNumber( filen[8:11] ):
749 760 setFile = int( filen[8:11] ) #inicializo mi contador de seteo al seteo del ultimo file
750 761 else:
751 762 setFile = -1
752 763 else:
753 764 setFile = -1 #inicializo mi contador de seteo
754 765
755 766 setFile += 1
756 767
757 768 file = '%s%4.4d%3.3d%3.3d%s' % (self.optchar,
758 769 timeTuple.tm_year,
759 770 timeTuple.tm_yday,
760 771 setFile,
761 772 ext )
762 773
763 774 filename = os.path.join( path, subfolder, file )
764 775
765 776 #Setting HDF5 File
766 777 fp = h5py.File(filename,'w')
767 778 grp = fp.create_group("Data")
768 779 grp.attrs['metadata'] = self.metaFile
769 780
770 781 # grp.attrs['blocksPerFile'] = 0
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,None) , chunks = True, dtype='S20')
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:
788 804 strMode = "channel"
789 805 nDatas[i] = self.arrayDim[i,1]
790 806 else:
791 807 strMode = "param"
792 808 nDatas[i] = self.arrayDim[i,2]
793 809
794 810 if nDims[i]==2:
795 811 nDatas[i] = self.arrayDim[i,2]
796 812
797 813 grp0 = grp.create_group(self.dataList[i])
798 814
799 815 for j in range(int(nDatas[i])):
800 816 tableName = strMode + str(j)
801 817
802 818 if nDims[i] == 3:
803 ds0 = grp0.create_dataset(tableName, (1,1,1) , maxshape=(None,None,None), chunks=True)
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,None), chunks=True)
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
820 838 self.setFile = setFile
821 839 self.firsttime = True
822 840 self.blockIndex = 0
823 841 return
824 842
825 843 def putData(self):
844
845 if not self.firsttime:
846 self.fp.flush()
847 self.fp.close()
848 self.readBlock()
849
850 if self.blockIndex == self.blocksPerFile:
851
852 self.setNextFile()
853
826 854 self.setBlock()
827 855 self.writeBlock()
828 856
829 if self.blockIndex == self.blocksPerFile:
830 self.setNextFile()
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
836 906
837 907
838 908 self.data
839 909 '''
840 910 #Creating Arrays
841 911 data = self.data
842 912 nDatas = self.nDatas
843 913 nDims = self.nDims
844 914 mode = self.mode
845 915 ind = 0
846 916
847 917 for i in range(len(self.dataList)):
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
859 929 data[ind] = dataAux[j,:]
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
877 947 return
878 948
879 949 def writeBlock(self):
880 950 '''
881 951 Saves the block in the HDF5 file
882 952 '''
883 953 for i in range(len(self.ds)):
884 self.ds[i].resize(self.data[i].shape)
885 self.ds[i][:] = self.data[i]
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
963 self.ds[i].resize(self.data[i].shape)
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
987
988 # if self.firsttime:
989 # self.fp.close()
990 # self.readBlock2()
886 991
887 992 self.blockIndex += 1
888
889 self.grp.attrs.modify('nRecords', self.blockIndex)
890
891 993 self.firsttime = False
892 994 return
893 995
894 996 def run(self, dataOut, **kwargs):
895 997 if not(self.isConfig):
896 998 self.setup(dataOut, **kwargs)
897 999 self.isConfig = True
898 1000 self.putMetadata()
899 1001 self.setNextFile()
900 1002
901 1003 self.putData()
902 1004 return
903 1005
@@ -1,1786 +1,2144
1 1 import numpy
2 2 import math
3 3 from scipy import optimize
4 4 from scipy import interpolate
5 5 from scipy import signal
6 6 from scipy import stats
7 7 import re
8 8 import datetime
9 9 import copy
10 10 import sys
11 11 import importlib
12 12 import itertools
13 13
14 14 from jroproc_base import ProcessingUnit, Operation
15 15 from schainpy.model.data.jrodata import Parameters
16 16
17 17
18 18 class ParametersProc(ProcessingUnit):
19 19
20 20 nSeconds = None
21 21
22 22 def __init__(self):
23 23 ProcessingUnit.__init__(self)
24 24
25 25 # self.objectDict = {}
26 26 self.buffer = None
27 27 self.firstdatatime = None
28 28 self.profIndex = 0
29 29 self.dataOut = Parameters()
30 30
31 31 def __updateObjFromInput(self):
32 32
33 33 self.dataOut.inputUnit = self.dataIn.type
34 34
35 35 self.dataOut.timeZone = self.dataIn.timeZone
36 36 self.dataOut.dstFlag = self.dataIn.dstFlag
37 37 self.dataOut.errorCount = self.dataIn.errorCount
38 38 self.dataOut.useLocalTime = self.dataIn.useLocalTime
39 39
40 40 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
41 41 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
42 42 self.dataOut.channelList = self.dataIn.channelList
43 43 self.dataOut.heightList = self.dataIn.heightList
44 44 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
45 45 # self.dataOut.nHeights = self.dataIn.nHeights
46 46 # self.dataOut.nChannels = self.dataIn.nChannels
47 47 self.dataOut.nBaud = self.dataIn.nBaud
48 48 self.dataOut.nCode = self.dataIn.nCode
49 49 self.dataOut.code = self.dataIn.code
50 50 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
51 51 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
52 52 self.dataOut.utctime = self.firstdatatime
53 53 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
54 54 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
55 55 # self.dataOut.nCohInt = self.dataIn.nCohInt
56 56 # self.dataOut.nIncohInt = 1
57 57 self.dataOut.ippSeconds = self.dataIn.ippSeconds
58 58 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
59 59 self.dataOut.timeInterval = self.dataIn.timeInterval
60 60 self.dataOut.heightList = self.dataIn.getHeiRange()
61 61 self.dataOut.frequency = self.dataIn.frequency
62 62
63 63 def run(self, nSeconds = None, nProfiles = None):
64 64
65 65
66 66
67 67 if self.firstdatatime == None:
68 68 self.firstdatatime = self.dataIn.utctime
69 69
70 70 #---------------------- Voltage Data ---------------------------
71 71
72 72 if self.dataIn.type == "Voltage":
73 73 self.dataOut.flagNoData = True
74 74 if nSeconds != None:
75 75 self.nSeconds = nSeconds
76 76 self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))
77 77
78 78 if self.buffer == None:
79 79 self.buffer = numpy.zeros((self.dataIn.nChannels,
80 80 self.nProfiles,
81 81 self.dataIn.nHeights),
82 82 dtype='complex')
83 83
84 84 self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
85 85 self.profIndex += 1
86 86
87 87 if self.profIndex == self.nProfiles:
88 88
89 89 self.__updateObjFromInput()
90 90 self.dataOut.data_pre = self.buffer.copy()
91 91 self.dataOut.paramInterval = nSeconds
92 92 self.dataOut.flagNoData = False
93 93
94 94 self.buffer = None
95 95 self.firstdatatime = None
96 96 self.profIndex = 0
97 97 return
98 98
99 99 #---------------------- Spectra Data ---------------------------
100 100
101 101 if self.dataIn.type == "Spectra":
102 102 self.dataOut.data_pre = self.dataIn.data_spc.copy()
103 103 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
104 104 self.dataOut.noise = self.dataIn.getNoise()
105 105 self.dataOut.normFactor = self.dataIn.normFactor
106 106 self.dataOut.groupList = self.dataIn.pairsList
107 107 self.dataOut.flagNoData = False
108 108
109 109 #---------------------- Correlation Data ---------------------------
110 110
111 111 if self.dataIn.type == "Correlation":
112 112 lagRRange = self.dataIn.lagR
113 113 indR = numpy.where(lagRRange == 0)[0][0]
114 114
115 115 self.dataOut.data_pre = self.dataIn.data_corr.copy()[:,:,indR,:]
116 116 self.dataOut.abscissaList = self.dataIn.getLagTRange(1)
117 117 self.dataOut.noise = self.dataIn.noise
118 118 self.dataOut.normFactor = self.dataIn.normFactor
119 119 self.dataOut.data_SNR = self.dataIn.SNR
120 120 self.dataOut.groupList = self.dataIn.pairsList
121 121 self.dataOut.flagNoData = False
122 122
123 123 #---------------------- Correlation Data ---------------------------
124 124
125 125 if self.dataIn.type == "Parameters":
126 126 self.dataOut.copy(self.dataIn)
127 127 self.dataOut.flagNoData = False
128 128
129 129 return True
130 130
131 131 self.__updateObjFromInput()
132 132 self.firstdatatime = None
133 133 self.dataOut.utctimeInit = self.dataIn.utctime
134 134 self.dataOut.outputInterval = self.dataIn.timeInterval
135 135
136 136 #------------------- Get Moments ----------------------------------
137 137 def GetMoments(self, channelList = None):
138 138 '''
139 139 Function GetMoments()
140 140
141 141 Input:
142 142 channelList : simple channel list to select e.g. [2,3,7]
143 143 self.dataOut.data_pre
144 144 self.dataOut.abscissaList
145 145 self.dataOut.noise
146 146
147 147 Affected:
148 148 self.dataOut.data_param
149 149 self.dataOut.data_SNR
150 150
151 151 '''
152 152 data = self.dataOut.data_pre
153 153 absc = self.dataOut.abscissaList[:-1]
154 154 noise = self.dataOut.noise
155 155
156 156 data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
157 157
158 158 if channelList== None:
159 159 channelList = self.dataIn.channelList
160 160 self.dataOut.channelList = channelList
161 161
162 162 for ind in channelList:
163 163 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
164 164
165 165 self.dataOut.data_param = data_param[:,1:,:]
166 166 self.dataOut.data_SNR = data_param[:,0]
167 167 return
168 168
169 169 def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
170 170
171 171 if (nicoh == None): nicoh = 1
172 172 if (graph == None): graph = 0
173 173 if (smooth == None): smooth = 0
174 174 elif (self.smooth < 3): smooth = 0
175 175
176 176 if (type1 == None): type1 = 0
177 177 if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
178 178 if (snrth == None): snrth = -3
179 179 if (dc == None): dc = 0
180 180 if (aliasing == None): aliasing = 0
181 181 if (oldfd == None): oldfd = 0
182 182 if (wwauto == None): wwauto = 0
183 183
184 184 if (n0 < 1.e-20): n0 = 1.e-20
185 185
186 186 freq = oldfreq
187 187 vec_power = numpy.zeros(oldspec.shape[1])
188 188 vec_fd = numpy.zeros(oldspec.shape[1])
189 189 vec_w = numpy.zeros(oldspec.shape[1])
190 190 vec_snr = numpy.zeros(oldspec.shape[1])
191 191
192 192 for ind in range(oldspec.shape[1]):
193 193
194 194 spec = oldspec[:,ind]
195 195 aux = spec*fwindow
196 196 max_spec = aux.max()
197 197 m = list(aux).index(max_spec)
198 198
199 199 #Smooth
200 200 if (smooth == 0): spec2 = spec
201 201 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
202 202
203 203 # Calculo de Momentos
204 204 bb = spec2[range(m,spec2.size)]
205 205 bb = (bb<n0).nonzero()
206 206 bb = bb[0]
207 207
208 208 ss = spec2[range(0,m + 1)]
209 209 ss = (ss<n0).nonzero()
210 210 ss = ss[0]
211 211
212 212 if (bb.size == 0):
213 213 bb0 = spec.size - 1 - m
214 214 else:
215 215 bb0 = bb[0] - 1
216 216 if (bb0 < 0):
217 217 bb0 = 0
218 218
219 219 if (ss.size == 0): ss1 = 1
220 220 else: ss1 = max(ss) + 1
221 221
222 222 if (ss1 > m): ss1 = m
223 223
224 224 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
225 225 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
226 226 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
227 227 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
228 228 snr = (spec2.mean()-n0)/n0
229 229
230 230 if (snr < 1.e-20) :
231 231 snr = 1.e-20
232 232
233 233 vec_power[ind] = power
234 234 vec_fd[ind] = fd
235 235 vec_w[ind] = w
236 236 vec_snr[ind] = snr
237 237
238 238 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
239 239 return moments
240 240
241 241 #------------------ Get SA Parameters --------------------------
242 def GetSAParameters(self):
243 data = self.dataOut.data_pre
244 crossdata = self.dataIn.data_cspc
245 a = 1
246
247
248 242
243 def GetSAParameters(self):
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)
258
259 #Normalizing Cross Spectra
260 norm_cspectra = numpy.zeros(cspectra.shape)
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],:])
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):
252 273 '''
253 274 Function GetMoments()
254 275
255 276 Input:
256 277 self.dataOut.data_pre
257 278 self.dataOut.abscissaList
258 279 self.dataOut.noise
259 280 self.dataOut.normFactor
260 281 self.dataOut.data_SNR
261 282 self.dataOut.groupList
262 283 self.dataOut.nChannels
263 284
264 285 Affected:
265 286 self.dataOut.data_param
266 287
267 288 '''
268 289
269 290 data = self.dataOut.data_pre
270 291 normFactor = self.dataOut.normFactor
271 292 nHeights = self.dataOut.nHeights
272 293 absc = self.dataOut.abscissaList[:-1]
273 294 noise = self.dataOut.noise
274 295 SNR = self.dataOut.data_SNR
275 296 pairsList = self.dataOut.groupList
276 297 nChannels = self.dataOut.nChannels
277 298 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
278 299 self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
279 300
280 301 dataNorm = numpy.abs(data)
281 302 for l in range(len(pairsList)):
282 303 dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]
283 304
284 305 self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)
285 306 self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, pairsAutoCorr, absc)
286 307 return
287 308
288 309 def __getPairsAutoCorr(self, pairsList, nChannels):
289 310
290 311 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
291 312
292 313 for l in range(len(pairsList)):
293 314 firstChannel = pairsList[l][0]
294 315 secondChannel = pairsList[l][1]
295 316
296 317 #Obteniendo pares de Autocorrelacion
297 318 if firstChannel == secondChannel:
298 319 pairsAutoCorr[firstChannel] = int(l)
299 320
300 321 pairsAutoCorr = pairsAutoCorr.astype(int)
301 322
302 323 pairsCrossCorr = range(len(pairsList))
303 324 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
304 325
305 326 return pairsAutoCorr, pairsCrossCorr
306 327
307 328 def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):
308 329
309 330 Pt0 = data.shape[1]/2
310 331 #Funcion de Autocorrelacion
311 332 dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0)
312 333
313 334 #Obtencion Indice de TauCross
314 335 indCross = data[pairsCrossCorr,:,:].argmax(axis = 1)
315 336 #Obtencion Indice de TauAuto
316 337 indAuto = numpy.zeros(indCross.shape,dtype = 'int')
317 338 CCValue = data[pairsCrossCorr,Pt0,:]
318 339 for i in range(pairsCrossCorr.size):
319 340 indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0)
320 341
321 342 #Obtencion de TauCross y TauAuto
322 343 tauCross = lagTRange[indCross]
323 344 tauAuto = lagTRange[indAuto]
324 345
325 346 Nan1, Nan2 = numpy.where(tauCross == lagTRange[0])
326 347
327 348 tauCross[Nan1,Nan2] = numpy.nan
328 349 tauAuto[Nan1,Nan2] = numpy.nan
329 350 tau = numpy.vstack((tauCross,tauAuto))
330 351
331 352 return tau
332 353
333 354 def __calculateLag1Phase(self, data, pairs, lagTRange):
334 355 data1 = stats.nanmean(data[pairs,:,:], axis = 0)
335 356 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
336 357
337 358 phase = numpy.angle(data1[lag1,:])
338 359
339 360 return phase
340 361 #------------------- Detect Meteors ------------------------------
341 362
342 def DetectMeteors(self, hei_ref = None, tauindex = 0,
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,
347 368 phaseThresh = 20, SNRThresh = 8,
348 369 hmin = 70, hmax=110, azimuth = 0) :
349 370
350 371 '''
351 372 Function DetectMeteors()
352 373 Project developed with paper:
353 374 HOLDSWORTH ET AL. 2004
354 375
355 376 Input:
356 377 self.dataOut.data_pre
357 378
358 379 centerReceiverIndex: From the channels, which is the center receiver
359 380
360 381 hei_ref: Height reference for the Beacon signal extraction
361 382 tauindex:
362 383 predefinedPhaseShifts: Predefined phase offset for the voltge signals
363 384
364 385 cohDetection: Whether to user Coherent detection or not
365 386 cohDet_timeStep: Coherent Detection calculation time step
366 387 cohDet_thresh: Coherent Detection phase threshold to correct phases
367 388
368 389 noise_timeStep: Noise calculation time step
369 390 noise_multiple: Noise multiple to define signal threshold
370 391
371 392 multDet_timeLimit: Multiple Detection Removal time limit in seconds
372 393 multDet_rangeLimit: Multiple Detection Removal range limit in km
373 394
374 395 phaseThresh: Maximum phase difference between receiver to be consider a meteor
375 396 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
376 397
377 398 hmin: Minimum Height of the meteor to use it in the further wind estimations
378 399 hmax: Maximum Height of the meteor to use it in the further wind estimations
379 400 azimuth: Azimuth angle correction
380 401
381 402 Affected:
382 403 self.dataOut.data_param
383 404
384 405 Rejection Criteria (Errors):
385 406 0: No error; analysis OK
386 407 1: SNR < SNR threshold
387 408 2: angle of arrival (AOA) ambiguously determined
388 409 3: AOA estimate not feasible
389 410 4: Large difference in AOAs obtained from different antenna baselines
390 411 5: echo at start or end of time series
391 412 6: echo less than 5 examples long; too short for analysis
392 413 7: echo rise exceeds 0.3s
393 414 8: echo decay time less than twice rise time
394 415 9: large power level before echo
395 416 10: large power level after echo
396 417 11: poor fit to amplitude for estimation of decay time
397 418 12: poor fit to CCF phase variation for estimation of radial drift velocity
398 419 13: height unresolvable echo: not valid height within 70 to 110 km
399 420 14: height ambiguous echo: more then one possible height within 70 to 110 km
400 421 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
401 422 16: oscilatory echo, indicating event most likely not an underdense echo
402 423
403 424 17: phase difference in meteor Reestimation
404 425
405 426 Data Storage:
406 427 Meteors for Wind Estimation (8):
407 428 Day Hour | Range Height
408 429 Azimuth Zenith errorCosDir
409 430 VelRad errorVelRad
410 431 TypeError
411 432
412 433 '''
413 434 #Get Beacon signal
414 435 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
415 436
416 437 if hei_ref != None:
417 438 newheis = numpy.where(self.dataOut.heightList>hei_ref)
418 439
419 440 heiRang = self.dataOut.getHeiRange()
420 441 #Pairs List
421 442 pairslist = []
422 443 nChannel = self.dataOut.nChannels
423 444 for i in range(nChannel):
424 445 if i != centerReceiverIndex:
425 446 pairslist.append((centerReceiverIndex,i))
426 447
427 448 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
428 449 # see if the user put in pre defined phase shifts
429 450 voltsPShift = self.dataOut.data_pre.copy()
430 451
431 452 if predefinedPhaseShifts != None:
432 453 hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
433 else:
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
444 472 voltsDC = numpy.mean(voltsPShift,1)
445 473 voltsDC = numpy.mean(voltsDC,1)
446 474 for i in range(voltsDC.shape[0]):
447 475 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
448 476
449 477 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
450 478 voltsPShift = voltsPShift[:,:,:newheis[0][0]]
451 479
452 480 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
453 481 #Coherent Detection
454 482 if cohDetection:
455 483 #use coherent detection to get the net power
456 484 cohDet_thresh = cohDet_thresh*numpy.pi/180
457 485 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist, cohDet_thresh)
458 486
459 487 #Non-coherent detection!
460 488 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
461 489 #********** END OF COH/NON-COH POWER CALCULATION**********************
462 490
463 491 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
464 492 #Get noise
465 493 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)
466 494 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
467 495 #Get signal threshold
468 496 signalThresh = noise_multiple*noise
469 497 #Meteor echoes detection
470 498 listMeteors = self.__findMeteors(powerNet, signalThresh)
471 499 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
472 500
473 501 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
474 502 #Parameters
475 503 heiRange = self.dataOut.getHeiRange()
476 504 rangeInterval = heiRange[1] - heiRange[0]
477 505 rangeLimit = multDet_rangeLimit/rangeInterval
478 506 timeLimit = multDet_timeLimit/self.dataOut.timeInterval
479 507 #Multiple detection removals
480 508 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
481 509 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
482 510
483 511 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
484 512 #Parameters
485 513 phaseThresh = phaseThresh*numpy.pi/180
486 514 thresh = [phaseThresh, noise_multiple, SNRThresh]
487 515 #Meteor reestimation (Errors N 1, 6, 12, 17)
488 516 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
489 517 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
490 518 #Estimation of decay times (Errors N 7, 8, 11)
491 519 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
492 520 #******************* END OF METEOR REESTIMATION *******************
493 521
494 522 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
495 523 #Calculating Radial Velocity (Error N 15)
496 524 radialStdThresh = 10
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 pairsList.append((0,3))
511 pairsList.append((1,2))
512 arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
513
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)
530 pairx = (0,3)
531 pairy = (1,2)
532 pairsList.append(pairx)
533 pairsList.append(pairy)
534
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 #***************************+ SAVE DATA IN HDF5 FORMAT **********************
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
526 565 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
527 566
528 567 minIndex = min(newheis[0])
529 568 maxIndex = max(newheis[0])
530 569
531 570 voltage = voltage0[:,:,minIndex:maxIndex+1]
532 571 nLength = voltage.shape[1]/n
533 572 nMin = 0
534 573 nMax = 0
535 574 phaseOffset = numpy.zeros((len(pairslist),n))
536 575
537 576 for i in range(n):
538 577 nMax += nLength
539 578 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
540 579 phaseCCF = numpy.mean(phaseCCF, axis = 2)
541 580 phaseOffset[:,i] = phaseCCF.transpose()
542 581 nMin = nMax
543 582 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
544 583
545 584 #Remove Outliers
546 585 factor = 2
547 586 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
548 587 dw = numpy.std(wt,axis = 1)
549 588 dw = dw.reshape((dw.size,1))
550 589 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
551 590 phaseOffset[ind] = numpy.nan
552 591 phaseOffset = stats.nanmean(phaseOffset, axis=1)
553 592
554 593 return phaseOffset
555 594
556 595 def __shiftPhase(self, data, phaseShift):
557 596 #this will shift the phase of a complex number
558 597 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
559 598 return dataShifted
560 599
561 600 def __estimatePhaseDifference(self, array, pairslist):
562 601 nChannel = array.shape[0]
563 602 nHeights = array.shape[2]
564 603 numPairs = len(pairslist)
565 604 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
566 605 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
567 606
568 607 #Correct phases
569 608 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
570 609 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
571 610
572 611 if indDer[0].shape[0] > 0:
573 612 for i in range(indDer[0].shape[0]):
574 613 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
575 614 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
576 615
577 616 # for j in range(numSides):
578 617 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
579 618 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
580 619 #
581 620 #Linear
582 621 phaseInt = numpy.zeros((numPairs,1))
583 622 angAllCCF = phaseCCF[:,[0,1,3,4],0]
584 623 for j in range(numPairs):
585 624 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
586 625 phaseInt[j] = fit[1]
587 626 #Phase Differences
588 627 phaseDiff = phaseInt - phaseCCF[:,2,:]
589 628 phaseArrival = phaseInt.reshape(phaseInt.size)
590 629
591 630 #Dealias
592 631 indAlias = numpy.where(phaseArrival > numpy.pi)
593 632 phaseArrival[indAlias] -= 2*numpy.pi
594 633 indAlias = numpy.where(phaseArrival < -numpy.pi)
595 634 phaseArrival[indAlias] += 2*numpy.pi
596 635
597 636 return phaseDiff, phaseArrival
598 637
599 638 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
600 639 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
601 640 #find the phase shifts of each channel over 1 second intervals
602 641 #only look at ranges below the beacon signal
603 642 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
604 643 numBlocks = int(volts.shape[1]/numProfPerBlock)
605 644 numHeights = volts.shape[2]
606 645 nChannel = volts.shape[0]
607 646 voltsCohDet = volts.copy()
608 647
609 648 pairsarray = numpy.array(pairslist)
610 649 indSides = pairsarray[:,1]
611 650 # indSides = numpy.array(range(nChannel))
612 651 # indSides = numpy.delete(indSides, indCenter)
613 652 #
614 653 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
615 654 listBlocks = numpy.array_split(volts, numBlocks, 1)
616 655
617 656 startInd = 0
618 657 endInd = 0
619 658
620 659 for i in range(numBlocks):
621 660 startInd = endInd
622 661 endInd = endInd + listBlocks[i].shape[1]
623 662
624 663 arrayBlock = listBlocks[i]
625 664 # arrayBlockCenter = listCenter[i]
626 665
627 666 #Estimate the Phase Difference
628 667 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
629 668 #Phase Difference RMS
630 669 arrayPhaseRMS = numpy.abs(phaseDiff)
631 670 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
632 671 indPhase = numpy.where(phaseRMSaux==4)
633 672 #Shifting
634 673 if indPhase[0].shape[0] > 0:
635 674 for j in range(indSides.size):
636 675 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
637 676 voltsCohDet[:,startInd:endInd,:] = arrayBlock
638 677
639 678 return voltsCohDet
640 679
641 680 def __calculateCCF(self, volts, pairslist ,laglist):
642 681
643 682 nHeights = volts.shape[2]
644 683 nPoints = volts.shape[1]
645 684 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
646 685
647 686 for i in range(len(pairslist)):
648 687 volts1 = volts[pairslist[i][0]]
649 688 volts2 = volts[pairslist[i][1]]
650 689
651 690 for t in range(len(laglist)):
652 691 idxT = laglist[t]
653 692 if idxT >= 0:
654 693 vStacked = numpy.vstack((volts2[idxT:,:],
655 694 numpy.zeros((idxT, nHeights),dtype='complex')))
656 695 else:
657 696 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
658 697 volts2[:(nPoints + idxT),:]))
659 698 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
660 699
661 700 vStacked = None
662 701 return voltsCCF
663 702
664 703 def __getNoise(self, power, timeSegment, timeInterval):
665 704 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
666 705 numBlocks = int(power.shape[0]/numProfPerBlock)
667 706 numHeights = power.shape[1]
668 707
669 708 listPower = numpy.array_split(power, numBlocks, 0)
670 709 noise = numpy.zeros((power.shape[0], power.shape[1]))
671 710 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
672 711
673 712 startInd = 0
674 713 endInd = 0
675 714
676 715 for i in range(numBlocks): #split por canal
677 716 startInd = endInd
678 717 endInd = endInd + listPower[i].shape[0]
679 718
680 719 arrayBlock = listPower[i]
681 720 noiseAux = numpy.mean(arrayBlock, 0)
682 721 # noiseAux = numpy.median(noiseAux)
683 722 # noiseAux = numpy.mean(arrayBlock)
684 723 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
685 724
686 725 noiseAux1 = numpy.mean(arrayBlock)
687 726 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
688 727
689 728 return noise, noise1
690 729
691 730 def __findMeteors(self, power, thresh):
692 731 nProf = power.shape[0]
693 732 nHeights = power.shape[1]
694 733 listMeteors = []
695 734
696 735 for i in range(nHeights):
697 736 powerAux = power[:,i]
698 737 threshAux = thresh[:,i]
699 738
700 739 indUPthresh = numpy.where(powerAux > threshAux)[0]
701 740 indDNthresh = numpy.where(powerAux <= threshAux)[0]
702 741
703 742 j = 0
704 743
705 744 while (j < indUPthresh.size - 2):
706 745 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
707 746 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
708 747 indDNthresh = indDNthresh[indDNAux]
709 748
710 749 if (indDNthresh.size > 0):
711 750 indEnd = indDNthresh[0] - 1
712 751 indInit = indUPthresh[j]
713 752
714 753 meteor = powerAux[indInit:indEnd + 1]
715 754 indPeak = meteor.argmax() + indInit
716 755 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
717 756
718 757 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
719 758 j = numpy.where(indUPthresh == indEnd)[0] + 1
720 759 else: j+=1
721 760 else: j+=1
722 761
723 762 return listMeteors
724 763
725 764 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
726 765
727 766 arrayMeteors = numpy.asarray(listMeteors)
728 767 listMeteors1 = []
729 768
730 769 while arrayMeteors.shape[0] > 0:
731 770 FLAs = arrayMeteors[:,4]
732 771 maxFLA = FLAs.argmax()
733 772 listMeteors1.append(arrayMeteors[maxFLA,:])
734 773
735 774 MeteorInitTime = arrayMeteors[maxFLA,1]
736 775 MeteorEndTime = arrayMeteors[maxFLA,3]
737 776 MeteorHeight = arrayMeteors[maxFLA,0]
738 777
739 778 #Check neighborhood
740 779 maxHeightIndex = MeteorHeight + rangeLimit
741 780 minHeightIndex = MeteorHeight - rangeLimit
742 781 minTimeIndex = MeteorInitTime - timeLimit
743 782 maxTimeIndex = MeteorEndTime + timeLimit
744 783
745 784 #Check Heights
746 785 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
747 786 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
748 787 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
749 788
750 789 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
751 790
752 791 return listMeteors1
753 792
754 793 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
755 794 numHeights = volts.shape[2]
756 795 nChannel = volts.shape[0]
757 796
758 797 thresholdPhase = thresh[0]
759 798 thresholdNoise = thresh[1]
760 799 thresholdDB = float(thresh[2])
761 800
762 801 thresholdDB1 = 10**(thresholdDB/10)
763 802 pairsarray = numpy.array(pairslist)
764 803 indSides = pairsarray[:,1]
765 804
766 805 pairslist1 = list(pairslist)
767 806 pairslist1.append((0,1))
768 807 pairslist1.append((3,4))
769 808
770 809 listMeteors1 = []
771 810 listPowerSeries = []
772 811 listVoltageSeries = []
773 812 #volts has the war data
774 813
775 814 if frequency == 30e6:
776 815 timeLag = 45*10**-3
777 816 else:
778 817 timeLag = 15*10**-3
779 818 lag = numpy.ceil(timeLag/timeInterval)
780 819
781 820 for i in range(len(listMeteors)):
782 821
783 822 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
784 823 meteorAux = numpy.zeros(16)
785 824
786 825 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
787 826 mHeight = listMeteors[i][0]
788 827 mStart = listMeteors[i][1]
789 828 mPeak = listMeteors[i][2]
790 829 mEnd = listMeteors[i][3]
791 830
792 831 #get the volt data between the start and end times of the meteor
793 832 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
794 833 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
795 834
796 835 #3.6. Phase Difference estimation
797 836 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
798 837
799 838 #3.7. Phase difference removal & meteor start, peak and end times reestimated
800 839 #meteorVolts0.- all Channels, all Profiles
801 840 meteorVolts0 = volts[:,:,mHeight]
802 841 meteorThresh = noise[:,mHeight]*thresholdNoise
803 842 meteorNoise = noise[:,mHeight]
804 843 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
805 844 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
806 845
807 846 #Times reestimation
808 847 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
809 848 if mStart1.size > 0:
810 849 mStart1 = mStart1[-1] + 1
811 850
812 851 else:
813 852 mStart1 = mPeak
814 853
815 854 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
816 855 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
817 856 if mEndDecayTime1.size == 0:
818 857 mEndDecayTime1 = powerNet0.size
819 858 else:
820 859 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
821 860 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
822 861
823 862 #meteorVolts1.- all Channels, from start to end
824 863 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
825 864 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
826 865 if meteorVolts2.shape[1] == 0:
827 866 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
828 867 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
829 868 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
830 869 ##################### END PARAMETERS REESTIMATION #########################
831 870
832 871 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
833 872 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
834 873 if meteorVolts2.shape[1] > 0:
835 874 #Phase Difference re-estimation
836 875 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
837 876 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
838 877 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
839 878 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
840 879 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
841 880
842 881 #Phase Difference RMS
843 882 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
844 883 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
845 884 #Data from Meteor
846 885 mPeak1 = powerNet1.argmax() + mStart1
847 886 mPeakPower1 = powerNet1.max()
848 887 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
849 888 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
850 889 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
851 890 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
852 891 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
853 892 #Vectorize
854 893 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
855 894 meteorAux[7:11] = phaseDiffint[0:4]
856 895
857 896 #Rejection Criterions
858 897 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
859 898 meteorAux[-1] = 17
860 899 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
861 900 meteorAux[-1] = 1
862 901
863 902
864 903 else:
865 904 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
866 905 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
867 906 PowerSeries = 0
868 907
869 908 listMeteors1.append(meteorAux)
870 909 listPowerSeries.append(PowerSeries)
871 910 listVoltageSeries.append(meteorVolts1)
872 911
873 912 return listMeteors1, listPowerSeries, listVoltageSeries
874 913
875 914 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
876 915
877 916 threshError = 10
878 917 #Depending if it is 30 or 50 MHz
879 918 if frequency == 30e6:
880 919 timeLag = 45*10**-3
881 920 else:
882 921 timeLag = 15*10**-3
883 922 lag = numpy.ceil(timeLag/timeInterval)
884 923
885 924 listMeteors1 = []
886 925
887 926 for i in range(len(listMeteors)):
888 927 meteorPower = listPower[i]
889 928 meteorAux = listMeteors[i]
890 929
891 930 if meteorAux[-1] == 0:
892 931
893 932 try:
894 933 indmax = meteorPower.argmax()
895 934 indlag = indmax + lag
896 935
897 936 y = meteorPower[indlag:]
898 937 x = numpy.arange(0, y.size)*timeLag
899 938
900 939 #first guess
901 940 a = y[0]
902 941 tau = timeLag
903 942 #exponential fit
904 943 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
905 944 y1 = self.__exponential_function(x, *popt)
906 945 #error estimation
907 946 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
908 947
909 948 decayTime = popt[1]
910 949 riseTime = indmax*timeInterval
911 950 meteorAux[11:13] = [decayTime, error]
912 951
913 952 #Table items 7, 8 and 11
914 953 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
915 954 meteorAux[-1] = 7
916 955 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
917 956 meteorAux[-1] = 8
918 957 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
919 958 meteorAux[-1] = 11
920 959
921 960
922 961 except:
923 962 meteorAux[-1] = 11
924 963
925 964
926 965 listMeteors1.append(meteorAux)
927 966
928 967 return listMeteors1
929 968
930 969 #Exponential Function
931 970
932 971 def __exponential_function(self, x, a, tau):
933 972 y = a*numpy.exp(-x/tau)
934 973 return y
935 974
936 975 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
937 976
938 977 pairslist1 = list(pairslist)
939 978 pairslist1.append((0,1))
940 979 pairslist1.append((3,4))
941 980 numPairs = len(pairslist1)
942 981 #Time Lag
943 982 timeLag = 45*10**-3
944 983 c = 3e8
945 984 lag = numpy.ceil(timeLag/timeInterval)
946 985 freq = 30e6
947 986
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
957 995
958 996 #get the volt data between the start and end times of the meteor
959 997 meteorVolts = listVolts[i]
960 998 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
961 999
962 1000 #Get CCF
963 1001 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
964 1002
965 1003 #Method 2
966 1004 slopes = numpy.zeros(numPairs)
967 1005 time = numpy.array([-2,-1,1,2])*timeInterval
968 1006 angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
969 1007
970 1008 #Correct phases
971 1009 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
972 1010 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
973 1011
974 1012 if indDer[0].shape[0] > 0:
975 1013 for i in range(indDer[0].shape[0]):
976 1014 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
977 1015 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
978 1016
979 1017 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
980 1018 for j in range(numPairs):
981 1019 fit = stats.linregress(time, angAllCCF[j,:])
982 1020 slopes[j] = fit[0]
983 1021
984 1022 #Remove Outlier
985 1023 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
986 1024 # slopes = numpy.delete(slopes,indOut)
987 1025 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
988 1026 # slopes = numpy.delete(slopes,indOut)
989 1027
990 1028 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
991 1029 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
992 1030 meteorAux[-2] = radialError
993 1031 meteorAux[-3] = radialVelocity
994 1032
995 1033 #Setting Error
996 1034 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
997 1035 if numpy.abs(radialVelocity) > 200:
998 1036 meteorAux[-1] = 15
999 1037 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
1000 1038 elif radialError > radialStdThresh:
1001 1039 meteorAux[-1] = 12
1002 1040
1003 1041 listMeteors1.append(meteorAux)
1004 1042 return listMeteors1
1005 1043
1006 1044 def __setNewArrays(self, listMeteors, date, heiRang):
1007 1045
1008 1046 #New arrays
1009 1047 arrayMeteors = numpy.array(listMeteors)
1010 arrayParameters = numpy.zeros((len(listMeteors),10))
1048 arrayParameters = numpy.zeros((len(listMeteors), 14))
1011 1049
1012 1050 #Date inclusion
1013 1051 date = re.findall(r'\((.*?)\)', date)
1014 1052 date = date[0].split(',')
1015 1053 date = map(int, date)
1016 1054 date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
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[:,0:3] = arrayMeteors[:,0:3]
1025 arrayParameters[:,-3:] = arrayMeteors[:,-3:]
1026
1027 return arrayMeteors, 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
1068
1069 return arrayParameters
1028 1070
1029 1071 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
1030 1072
1031 1073 arrayAOA = numpy.zeros((phases.shape[0],3))
1032 1074 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
1033 1075
1034 1076 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
1035 1077 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
1036 1078 arrayAOA[:,2] = cosDirError
1037 1079
1038 1080 azimuthAngle = arrayAOA[:,0]
1039 1081 zenithAngle = arrayAOA[:,1]
1040 1082
1041 1083 #Setting Error
1042 1084 #Number 3: AOA not fesible
1043 1085 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
1044 1086 error[indInvalid] = 3
1045 1087 #Number 4: Large difference in AOAs obtained from different antenna baselines
1046 1088 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
1047 1089 error[indInvalid] = 4
1048 1090 return arrayAOA, error
1049 1091
1050 1092 def __getDirectionCosines(self, arrayPhase, pairsList):
1051 1093
1052 1094 #Initializing some variables
1053 1095 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
1054 1096 ang_aux = ang_aux.reshape(1,ang_aux.size)
1055 1097
1056 1098 cosdir = numpy.zeros((arrayPhase.shape[0],2))
1057 1099 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
1058 1100
1059 1101
1060 1102 for i in range(2):
1061 1103 #First Estimation
1062 1104 phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
1063 1105 #Dealias
1064 1106 indcsi = numpy.where(phi0_aux > numpy.pi)
1065 1107 phi0_aux[indcsi] -= 2*numpy.pi
1066 1108 indcsi = numpy.where(phi0_aux < -numpy.pi)
1067 1109 phi0_aux[indcsi] += 2*numpy.pi
1068 1110 #Direction Cosine 0
1069 1111 cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
1070 1112
1071 1113 #Most-Accurate Second Estimation
1072 1114 phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
1073 1115 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
1074 1116 #Direction Cosine 1
1075 1117 cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
1076 1118
1077 1119 #Searching the correct Direction Cosine
1078 1120 cosdir0_aux = cosdir0[:,i]
1079 1121 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
1080 1122 #Minimum Distance
1081 1123 cosDiff = (cosdir1 - cosdir0_aux)**2
1082 1124 indcos = cosDiff.argmin(axis = 1)
1083 1125 #Saving Value obtained
1084 1126 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
1085 1127
1086 1128 return cosdir0, cosdir
1087 1129
1088 1130 def __calculateAOA(self, cosdir, azimuth):
1089 1131 cosdirX = cosdir[:,0]
1090 1132 cosdirY = cosdir[:,1]
1091 1133
1092 1134 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
1093 1135 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
1094 1136 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
1095 1137
1096 1138 return angles
1097 1139
1098 1140 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
1099 1141
1100 1142 Ramb = 375 #Ramb = c/(2*PRF)
1101 1143 Re = 6371 #Earth Radius
1102 1144 heights = numpy.zeros(Ranges.shape)
1103 1145
1104 1146 R_aux = numpy.array([0,1,2])*Ramb
1105 1147 R_aux = R_aux.reshape(1,R_aux.size)
1106 1148
1107 1149 Ranges = Ranges.reshape(Ranges.size,1)
1108 1150
1109 1151 Ri = Ranges + R_aux
1110 1152 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
1111 1153
1112 1154 #Check if there is a height between 70 and 110 km
1113 1155 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
1114 1156 ind_h = numpy.where(h_bool == 1)[0]
1115 1157
1116 1158 hCorr = hi[ind_h, :]
1117 1159 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
1118 1160
1119 1161 hCorr = hi[ind_hCorr]
1120 1162 heights[ind_h] = hCorr
1121 1163
1122 1164 #Setting Error
1123 1165 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
1124 1166 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
1125 1167
1126 1168 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
1127 1169 error[indInvalid2] = 14
1128 1170 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
1129 1171 error[indInvalid1] = 13
1130 1172
1131 1173 return heights, error
1132 1174
1133 1175 def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None):
1134 1176
1135 1177 '''
1136 1178 Function GetMoments()
1137 1179
1138 1180 Input:
1139 1181 Output:
1140 1182 Variables modified:
1141 1183 '''
1142 1184 if path != None:
1143 1185 sys.path.append(path)
1144 1186 self.dataOut.library = importlib.import_module(file)
1145 1187
1146 1188 #To be inserted as a parameter
1147 1189 groupArray = numpy.array(groupList)
1148 1190 # groupArray = numpy.array([[0,1],[2,3]])
1149 1191 self.dataOut.groupList = groupArray
1150 1192
1151 1193 nGroups = groupArray.shape[0]
1152 1194 nChannels = self.dataIn.nChannels
1153 1195 nHeights=self.dataIn.heightList.size
1154 1196
1155 1197 #Parameters Array
1156 1198 self.dataOut.data_param = None
1157 1199
1158 1200 #Set constants
1159 1201 constants = self.dataOut.library.setConstants(self.dataIn)
1160 1202 self.dataOut.constants = constants
1161 1203 M = self.dataIn.normFactor
1162 1204 N = self.dataIn.nFFTPoints
1163 1205 ippSeconds = self.dataIn.ippSeconds
1164 1206 K = self.dataIn.nIncohInt
1165 1207 pairsArray = numpy.array(self.dataIn.pairsList)
1166 1208
1167 1209 #List of possible combinations
1168 1210 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
1169 1211 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
1170 1212
1171 1213 if getSNR:
1172 1214 listChannels = groupArray.reshape((groupArray.size))
1173 1215 listChannels.sort()
1174 1216 noise = self.dataIn.getNoise()
1175 1217 self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
1176 1218
1177 1219 for i in range(nGroups):
1178 1220 coord = groupArray[i,:]
1179 1221
1180 1222 #Input data array
1181 1223 data = self.dataIn.data_spc[coord,:,:]/(M*N)
1182 1224 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
1183 1225
1184 1226 #Cross Spectra data array for Covariance Matrixes
1185 1227 ind = 0
1186 1228 for pairs in listComb:
1187 1229 pairsSel = numpy.array([coord[x],coord[y]])
1188 1230 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
1189 1231 ind += 1
1190 1232 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
1191 1233 dataCross = dataCross**2/K
1192 1234
1193 1235 for h in range(nHeights):
1194 1236 # print self.dataOut.heightList[h]
1195 1237
1196 1238 #Input
1197 1239 d = data[:,h]
1198 1240
1199 1241 #Covariance Matrix
1200 1242 D = numpy.diag(d**2/K)
1201 1243 ind = 0
1202 1244 for pairs in listComb:
1203 1245 #Coordinates in Covariance Matrix
1204 1246 x = pairs[0]
1205 1247 y = pairs[1]
1206 1248 #Channel Index
1207 1249 S12 = dataCross[ind,:,h]
1208 1250 D12 = numpy.diag(S12)
1209 1251 #Completing Covariance Matrix with Cross Spectras
1210 1252 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
1211 1253 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
1212 1254 ind += 1
1213 1255 Dinv=numpy.linalg.inv(D)
1214 1256 L=numpy.linalg.cholesky(Dinv)
1215 1257 LT=L.T
1216 1258
1217 1259 dp = numpy.dot(LT,d)
1218 1260
1219 1261 #Initial values
1220 1262 data_spc = self.dataIn.data_spc[coord,:,h]
1221 1263
1222 1264 if (h>0)and(error1[3]<5):
1223 1265 p0 = self.dataOut.data_param[i,:,h-1]
1224 1266 else:
1225 1267 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
1226 1268
1227 1269 try:
1228 1270 #Least Squares
1229 1271 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
1230 1272 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
1231 1273 #Chi square error
1232 1274 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
1233 1275 #Error with Jacobian
1234 1276 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
1235 1277 except:
1236 1278 minp = p0*numpy.nan
1237 1279 error0 = numpy.nan
1238 1280 error1 = p0*numpy.nan
1239 1281
1240 1282 #Save
1241 1283 if self.dataOut.data_param == None:
1242 1284 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
1243 1285 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
1244 1286
1245 1287 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
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)
1253 1294 fmp=numpy.dot(LT,fm)
1254 1295
1255 1296 return dp-fmp
1256 1297
1257 1298 def __getSNR(self, z, noise):
1258 1299
1259 1300 avg = numpy.average(z, axis=1)
1260 1301 SNR = (avg.T-noise)/noise
1261 1302 SNR = SNR.T
1262 1303 return SNR
1263 1304
1264 1305 def __chisq(p,chindex,hindex):
1265 1306 #similar to Resid but calculates CHI**2
1266 1307 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
1267 1308 dp=numpy.dot(LT,d)
1268 1309 fmp=numpy.dot(LT,fm)
1269 1310 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
1270 1311 return chisq
1271 1312
1272 1313
1273 1314
1274 1315 class WindProfiler(Operation):
1275 1316
1276 1317 __isConfig = False
1277 1318
1278 1319 __initime = None
1279 1320 __lastdatatime = None
1280 1321 __integrationtime = None
1281 1322
1282 1323 __buffer = None
1283 1324
1284 1325 __dataReady = False
1285 1326
1286 1327 __firstdata = None
1287 1328
1288 1329 n = None
1289 1330
1290 1331 def __init__(self):
1291 1332 Operation.__init__(self)
1292 1333
1293 1334 def __calculateCosDir(self, elev, azim):
1294 1335 zen = (90 - elev)*numpy.pi/180
1295 1336 azim = azim*numpy.pi/180
1296 1337 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
1297 1338 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
1298 1339
1299 1340 signX = numpy.sign(numpy.cos(azim))
1300 1341 signY = numpy.sign(numpy.sin(azim))
1301 1342
1302 1343 cosDirX = numpy.copysign(cosDirX, signX)
1303 1344 cosDirY = numpy.copysign(cosDirY, signY)
1304 1345 return cosDirX, cosDirY
1305 1346
1306 1347 def __calculateAngles(self, theta_x, theta_y, azimuth):
1307 1348
1308 1349 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
1309 1350 zenith_arr = numpy.arccos(dir_cosw)
1310 1351 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
1311 1352
1312 1353 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
1313 1354 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
1314 1355
1315 1356 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
1316 1357
1317 1358 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
1318 1359
1319 1360 #
1320 1361 if horOnly:
1321 1362 A = numpy.c_[dir_cosu,dir_cosv]
1322 1363 else:
1323 1364 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
1324 1365 A = numpy.asmatrix(A)
1325 1366 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
1326 1367
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))
1334 1374
1335 1375 rango = range(len(phi))
1336 1376 # rango = numpy.delete(rango,maxid)
1337 1377
1338 1378 heiRang1 = heiRang*math.cos(phi[maxid])
1339 1379 heiRangAux = heiRang*math.cos(phi[minid])
1340 1380 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1341 1381 heiRang1 = numpy.delete(heiRang1,indOut)
1342 1382
1343 1383 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1344 1384 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1345 1385
1346 1386 for i in rango:
1347 1387 x = heiRang*math.cos(phi[i])
1348 1388 y1 = velRadial[i,:]
1349 1389 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1350 1390
1351 1391 x1 = heiRang1
1352 1392 y11 = f1(x1)
1353 1393
1354 1394 y2 = SNR[i,:]
1355 1395 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1356 1396 y21 = f2(x1)
1357 1397
1358 1398 velRadial1[i,:] = y11
1359 1399 SNR1[i,:] = y21
1360 1400
1361 1401 return heiRang1, velRadial1, SNR1
1362 1402
1363 1403 def __calculateVelUVW(self, A, velRadial):
1364 1404
1365 1405 #Operacion Matricial
1366 1406 # velUVW = numpy.zeros((velRadial.shape[1],3))
1367 1407 # for ind in range(velRadial.shape[1]):
1368 1408 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
1369 1409 # velUVW = velUVW.transpose()
1370 1410 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
1371 1411 velUVW[:,:] = numpy.dot(A,velRadial)
1372 1412
1373 1413
1374 1414 return velUVW
1375 1415
1376 1416 def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
1377 1417 """
1378 1418 Function that implements Doppler Beam Swinging (DBS) technique.
1379 1419
1380 1420 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1381 1421 Direction correction (if necessary), Ranges and SNR
1382 1422
1383 1423 Output: Winds estimation (Zonal, Meridional and Vertical)
1384 1424
1385 1425 Parameters affected: Winds, height range, SNR
1386 1426 """
1387 1427 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth)
1388 1428 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0)
1389 1429 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
1390 1430
1391 1431 #Calculo de Componentes de la velocidad con DBS
1392 1432 winds = self.__calculateVelUVW(A,velRadial1)
1393 1433
1394 1434 return winds, heiRang1, SNR1
1395 1435
1396 1436 def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):
1397 1437
1398 1438 posx = numpy.asarray(posx)
1399 1439 posy = numpy.asarray(posy)
1400 1440
1401 1441 #Rotacion Inversa para alinear con el azimuth
1402 1442 if azimuth!= None:
1403 1443 azimuth = azimuth*math.pi/180
1404 1444 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
1405 1445 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
1406 1446 else:
1407 1447 posx1 = posx
1408 1448 posy1 = posy
1409 1449
1410 1450 #Calculo de Distancias
1411 1451 distx = numpy.zeros(pairsCrossCorr.size)
1412 1452 disty = numpy.zeros(pairsCrossCorr.size)
1413 1453 dist = numpy.zeros(pairsCrossCorr.size)
1414 1454 ang = numpy.zeros(pairsCrossCorr.size)
1415 1455
1416 1456 for i in range(pairsCrossCorr.size):
1417 1457 distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]
1418 1458 disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]
1419 1459 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
1420 1460 ang[i] = numpy.arctan2(disty[i],distx[i])
1421 1461 #Calculo de Matrices
1422 1462 nPairs = len(pairs)
1423 1463 ang1 = numpy.zeros((nPairs, 2, 1))
1424 1464 dist1 = numpy.zeros((nPairs, 2, 1))
1425 1465
1426 1466 for j in range(nPairs):
1427 1467 dist1[j,0,0] = dist[pairs[j][0]]
1428 1468 dist1[j,1,0] = dist[pairs[j][1]]
1429 1469 ang1[j,0,0] = ang[pairs[j][0]]
1430 1470 ang1[j,1,0] = ang[pairs[j][1]]
1431 1471
1432 1472 return distx,disty, dist1,ang1
1433 1473
1434 1474 def __calculateVelVer(self, phase, lagTRange, _lambda):
1435 1475
1436 1476 Ts = lagTRange[1] - lagTRange[0]
1437 1477 velW = -_lambda*phase/(4*math.pi*Ts)
1438 1478
1439 1479 return velW
1440 1480
1441 1481 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
1442 1482 nPairs = tau1.shape[0]
1443 1483 vel = numpy.zeros((nPairs,3,tau1.shape[2]))
1444 1484
1445 1485 angCos = numpy.cos(ang)
1446 1486 angSin = numpy.sin(ang)
1447 1487
1448 1488 vel0 = dist*tau1/(2*tau2**2)
1449 1489 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
1450 1490 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
1451 1491
1452 1492 ind = numpy.where(numpy.isinf(vel))
1453 1493 vel[ind] = numpy.nan
1454 1494
1455 1495 return vel
1456 1496
1457 1497 def __getPairsAutoCorr(self, pairsList, nChannels):
1458 1498
1459 1499 pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
1460 1500
1461 1501 for l in range(len(pairsList)):
1462 1502 firstChannel = pairsList[l][0]
1463 1503 secondChannel = pairsList[l][1]
1464 1504
1465 1505 #Obteniendo pares de Autocorrelacion
1466 1506 if firstChannel == secondChannel:
1467 1507 pairsAutoCorr[firstChannel] = int(l)
1468 1508
1469 1509 pairsAutoCorr = pairsAutoCorr.astype(int)
1470 1510
1471 1511 pairsCrossCorr = range(len(pairsList))
1472 1512 pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
1473 1513
1474 1514 return pairsAutoCorr, pairsCrossCorr
1475 1515
1476 1516 def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
1477 1517 """
1478 1518 Function that implements Spaced Antenna (SA) technique.
1479 1519
1480 1520 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
1481 1521 Direction correction (if necessary), Ranges and SNR
1482 1522
1483 1523 Output: Winds estimation (Zonal, Meridional and Vertical)
1484 1524
1485 1525 Parameters affected: Winds
1486 1526 """
1487 1527 #Cross Correlation pairs obtained
1488 1528 pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
1489 1529 pairsArray = numpy.array(pairsList)[pairsCrossCorr]
1490 1530 pairsSelArray = numpy.array(pairsSelected)
1491 1531 pairs = []
1492 1532
1493 1533 #Wind estimation pairs obtained
1494 1534 for i in range(pairsSelArray.shape[0]/2):
1495 1535 ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
1496 1536 ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
1497 1537 pairs.append((ind1,ind2))
1498 1538
1499 1539 indtau = tau.shape[0]/2
1500 1540 tau1 = tau[:indtau,:]
1501 1541 tau2 = tau[indtau:-1,:]
1502 1542 tau1 = tau1[pairs,:]
1503 1543 tau2 = tau2[pairs,:]
1504 1544 phase1 = tau[-1,:]
1505 1545
1506 1546 #---------------------------------------------------------------------
1507 1547 #Metodo Directo
1508 1548 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)
1509 1549 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
1510 1550 winds = stats.nanmean(winds, axis=0)
1511 1551 #---------------------------------------------------------------------
1512 1552 #Metodo General
1513 1553 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
1514 1554 # #Calculo Coeficientes de Funcion de Correlacion
1515 1555 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
1516 1556 # #Calculo de Velocidades
1517 1557 # winds = self.calculateVelUV(F,G,A,B,H)
1518 1558
1519 1559 #---------------------------------------------------------------------
1520 1560 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
1521 1561 winds = correctFactor*winds
1522 1562 return winds
1523 1563
1524 1564 def __checkTime(self, currentTime, paramInterval, outputInterval):
1525 1565
1526 1566 dataTime = currentTime + paramInterval
1527 1567 deltaTime = dataTime - self.__initime
1528 1568
1529 1569 if deltaTime >= outputInterval or deltaTime < 0:
1530 1570 self.__dataReady = True
1531 1571 return
1532 1572
1533 1573 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
1534 1574 '''
1535 1575 Function that implements winds estimation technique with detected meteors.
1536 1576
1537 1577 Input: Detected meteors, Minimum meteor quantity to wind estimation
1538 1578
1539 1579 Output: Winds estimation (Zonal and Meridional)
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]
1553 1594 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
1554 1595 nMeteorsPerI = hist[0]
1555 1596 heightPerI = hist[1]
1556 1597
1557 1598 #Sort of meteors
1558 1599 indSort = finalHeights.argsort()
1559 1600 finalMeteor2 = finalMeteor[indSort,:]
1560 1601
1561 1602 # Calculating winds
1562 1603 ind1 = 0
1563 1604 ind2 = 0
1564 1605
1565 1606 for i in range(nInt):
1566 1607 nMet = nMeteorsPerI[i]
1567 1608 ind1 = ind2
1568 1609 ind2 = ind1 + nMet
1569 1610
1570 1611 meteorAux = finalMeteor2[ind1:ind2,:]
1571 1612
1572 1613 if meteorAux.shape[0] >= meteorThresh:
1573 1614 vel = meteorAux[:, 7]
1574 1615 zen = meteorAux[:, 5]*numpy.pi/180
1575 1616 azim = meteorAux[:, 4]*numpy.pi/180
1576 1617
1577 1618 n = numpy.cos(zen)
1578 1619 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
1579 1620 # l = m*numpy.tan(azim)
1580 1621 l = numpy.sin(zen)*numpy.sin(azim)
1581 1622 m = numpy.sin(zen)*numpy.cos(azim)
1582 1623
1583 1624 A = numpy.vstack((l, m)).transpose()
1584 1625 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
1585 1626 windsAux = numpy.dot(A1, vel)
1586 1627
1587 1628 winds[0,i] = windsAux[0]
1588 1629 winds[1,i] = windsAux[1]
1589 1630
1590 1631 return winds, heightPerI[:-1]
1591 1632
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.getHeiRange()
1639 heightList = dataOut.heightList
1599 1640 SNR = dataOut.data_SNR
1600 1641
1601 1642 if technique == 'DBS':
1602 1643
1603 1644 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
1604 1645 theta_x = numpy.array(kwargs['dirCosx'])
1605 1646 theta_y = numpy.array(kwargs['dirCosy'])
1606 1647 else:
1607 1648 elev = numpy.array(kwargs['elevation'])
1608 1649 azim = numpy.array(kwargs['azimuth'])
1609 1650 theta_x, theta_y = self.__calculateCosDir(elev, azim)
1610 1651 azimuth = kwargs['correctAzimuth']
1611 1652 if kwargs.has_key('horizontalOnly'):
1612 1653 horizontalOnly = kwargs['horizontalOnly']
1613 1654 else: horizontalOnly = False
1614 1655 if kwargs.has_key('correctFactor'):
1615 1656 correctFactor = kwargs['correctFactor']
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
1635 1671
1636 1672 elif technique == 'SA':
1637 1673
1638 1674 #Parameters
1639 1675 position_x = kwargs['positionX']
1640 1676 position_y = kwargs['positionY']
1641 1677 azimuth = kwargs['azimuth']
1642 1678
1643 1679 if kwargs.has_key('crosspairsList'):
1644 1680 pairs = kwargs['crosspairsList']
1645 1681 else:
1646 1682 pairs = None
1647 1683
1648 1684 if kwargs.has_key('correctFactor'):
1649 1685 correctFactor = kwargs['correctFactor']
1650 1686 else:
1651 1687 correctFactor = 1
1652 1688
1653 1689 tau = dataOut.data_param
1654 1690 _lambda = dataOut.C/dataOut.frequency
1655 1691 pairsList = dataOut.groupList
1656 1692 nChannels = dataOut.nChannels
1657 1693
1658 1694 dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1659 1695 dataOut.utctimeInit = dataOut.utctime
1660 1696 dataOut.outputInterval = dataOut.timeInterval
1661 1697
1662 1698 elif technique == 'Meteors':
1663 1699 dataOut.flagNoData = True
1664 1700 self.__dataReady = False
1665 1701
1666 1702 if kwargs.has_key('nHours'):
1667 1703 nHours = kwargs['nHours']
1668 1704 else:
1669 1705 nHours = 1
1670 1706
1671 1707 if kwargs.has_key('meteorsPerBin'):
1672 1708 meteorThresh = kwargs['meteorsPerBin']
1673 1709 else:
1674 1710 meteorThresh = 6
1675 1711
1676 1712 if kwargs.has_key('hmin'):
1677 1713 hmin = kwargs['hmin']
1678 1714 else: hmin = 70
1679 1715 if kwargs.has_key('hmax'):
1680 1716 hmax = kwargs['hmax']
1681 1717 else: hmax = 110
1682 1718
1683 1719 dataOut.outputInterval = nHours*3600
1684 1720
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(self.dataOut.utctime)
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
1692 1728
1693 1729 if self.__buffer == None:
1694 1730 self.__buffer = dataOut.data_param
1695 1731 self.__firstdata = copy.copy(dataOut)
1696 1732
1697 1733 else:
1698 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
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
1702 1738 if self.__dataReady:
1703 1739 dataOut.utctimeInit = self.__initime
1704 1740
1705 1741 self.__initime += dataOut.outputInterval #to erase time offset
1706 1742
1707 1743 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
1708 1744 dataOut.flagNoData = False
1709 1745 self.__buffer = None
1710 1746
1711 1747 return
1712 1748
1713 1749 class EWDriftsEstimation(Operation):
1714 1750
1715 1751
1716 1752 def __init__(self):
1717 1753 Operation.__init__(self)
1718 1754
1719 1755 def __correctValues(self, heiRang, phi, velRadial, SNR):
1720 1756 listPhi = phi.tolist()
1721 1757 maxid = listPhi.index(max(listPhi))
1722 1758 minid = listPhi.index(min(listPhi))
1723 1759
1724 1760 rango = range(len(phi))
1725 1761 # rango = numpy.delete(rango,maxid)
1726 1762
1727 1763 heiRang1 = heiRang*math.cos(phi[maxid])
1728 1764 heiRangAux = heiRang*math.cos(phi[minid])
1729 1765 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1730 1766 heiRang1 = numpy.delete(heiRang1,indOut)
1731 1767
1732 1768 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1733 1769 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1734 1770
1735 1771 for i in rango:
1736 1772 x = heiRang*math.cos(phi[i])
1737 1773 y1 = velRadial[i,:]
1738 1774 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1739 1775
1740 1776 x1 = heiRang1
1741 1777 y11 = f1(x1)
1742 1778
1743 1779 y2 = SNR[i,:]
1744 1780 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1745 1781 y21 = f2(x1)
1746 1782
1747 1783 velRadial1[i,:] = y11
1748 1784 SNR1[i,:] = y21
1749 1785
1750 1786 return heiRang1, velRadial1, SNR1
1751 1787
1752 1788 def run(self, dataOut, zenith, zenithCorrection):
1753 1789 heiRang = dataOut.heightList
1754 1790 velRadial = dataOut.data_param[:,3,:]
1755 1791 SNR = dataOut.data_SNR
1756 1792
1757 1793 zenith = numpy.array(zenith)
1758 1794 zenith -= zenithCorrection
1759 1795 zenith *= numpy.pi/180
1760 1796
1761 1797 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
1762 1798
1763 1799 alp = zenith[0]
1764 1800 bet = zenith[1]
1765 1801
1766 1802 w_w = velRadial1[0,:]
1767 1803 w_e = velRadial1[1,:]
1768 1804
1769 1805 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
1770 1806 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
1771 1807
1772 1808 winds = numpy.vstack((u,w))
1773 1809
1774 1810 dataOut.heightList = heiRang1
1775 1811 dataOut.data_output = winds
1776 1812 dataOut.data_SNR = SNR1
1777 1813
1778 1814 dataOut.utctimeInit = dataOut.utctime
1779 1815 dataOut.outputInterval = dataOut.timeInterval
1780 1816 return
1781 1817
1818 class PhaseCalibration(Operation):
1819
1820 __buffer = None
1821
1822 __initime = None
1782 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):
2014
2015 arrayParameters = arrayParameters0.copy()
2016 hmin = h[0]
2017 hmax = h[1]
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
1783 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):
1784 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
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