##// END OF EJS Templates
Add ScopePlot to the new way of plotting data (jroplot_data)
George Yong -
r1202:179ac651dbc0
parent child
Show More
@@ -1,1356 +1,1363
1 '''
1 '''
2
2
3 $Author: murco $
3 $Author: murco $
4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 '''
5 '''
6
6
7 import copy
7 import copy
8 import numpy
8 import numpy
9 import datetime
9 import datetime
10 import json
10 import json
11
11
12 from schainpy.utils import log
12 from schainpy.utils import log
13 from .jroheaderIO import SystemHeader, RadarControllerHeader
13 from .jroheaderIO import SystemHeader, RadarControllerHeader
14
14
15
15
16 def getNumpyDtype(dataTypeCode):
16 def getNumpyDtype(dataTypeCode):
17
17
18 if dataTypeCode == 0:
18 if dataTypeCode == 0:
19 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
19 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
20 elif dataTypeCode == 1:
20 elif dataTypeCode == 1:
21 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
21 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
22 elif dataTypeCode == 2:
22 elif dataTypeCode == 2:
23 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
23 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
24 elif dataTypeCode == 3:
24 elif dataTypeCode == 3:
25 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
25 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
26 elif dataTypeCode == 4:
26 elif dataTypeCode == 4:
27 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
27 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
28 elif dataTypeCode == 5:
28 elif dataTypeCode == 5:
29 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
29 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
30 else:
30 else:
31 raise ValueError('dataTypeCode was not defined')
31 raise ValueError('dataTypeCode was not defined')
32
32
33 return numpyDtype
33 return numpyDtype
34
34
35
35
36 def getDataTypeCode(numpyDtype):
36 def getDataTypeCode(numpyDtype):
37
37
38 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
38 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
39 datatype = 0
39 datatype = 0
40 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
40 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
41 datatype = 1
41 datatype = 1
42 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
42 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
43 datatype = 2
43 datatype = 2
44 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
44 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
45 datatype = 3
45 datatype = 3
46 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
46 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
47 datatype = 4
47 datatype = 4
48 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
48 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
49 datatype = 5
49 datatype = 5
50 else:
50 else:
51 datatype = None
51 datatype = None
52
52
53 return datatype
53 return datatype
54
54
55
55
56 def hildebrand_sekhon(data, navg):
56 def hildebrand_sekhon(data, navg):
57 """
57 """
58 This method is for the objective determination of the noise level in Doppler spectra. This
58 This method is for the objective determination of the noise level in Doppler spectra. This
59 implementation technique is based on the fact that the standard deviation of the spectral
59 implementation technique is based on the fact that the standard deviation of the spectral
60 densities is equal to the mean spectral density for white Gaussian noise
60 densities is equal to the mean spectral density for white Gaussian noise
61
61
62 Inputs:
62 Inputs:
63 Data : heights
63 Data : heights
64 navg : numbers of averages
64 navg : numbers of averages
65
65
66 Return:
66 Return:
67 mean : noise's level
67 mean : noise's level
68 """
68 """
69
69
70 sortdata = numpy.sort(data, axis=None)
70 sortdata = numpy.sort(data, axis=None)
71 lenOfData = len(sortdata)
71 lenOfData = len(sortdata)
72 nums_min = lenOfData*0.2
72 nums_min = lenOfData*0.2
73
73
74 if nums_min <= 5:
74 if nums_min <= 5:
75
75
76 nums_min = 5
76 nums_min = 5
77
77
78 sump = 0.
78 sump = 0.
79 sumq = 0.
79 sumq = 0.
80
80
81 j = 0
81 j = 0
82 cont = 1
82 cont = 1
83
83
84 while((cont == 1)and(j < lenOfData)):
84 while((cont == 1)and(j < lenOfData)):
85
85
86 sump += sortdata[j]
86 sump += sortdata[j]
87 sumq += sortdata[j]**2
87 sumq += sortdata[j]**2
88
88
89 if j > nums_min:
89 if j > nums_min:
90 rtest = float(j)/(j-1) + 1.0/navg
90 rtest = float(j)/(j-1) + 1.0/navg
91 if ((sumq*j) > (rtest*sump**2)):
91 if ((sumq*j) > (rtest*sump**2)):
92 j = j - 1
92 j = j - 1
93 sump = sump - sortdata[j]
93 sump = sump - sortdata[j]
94 sumq = sumq - sortdata[j]**2
94 sumq = sumq - sortdata[j]**2
95 cont = 0
95 cont = 0
96
96
97 j += 1
97 j += 1
98
98
99 lnoise = sump / j
99 lnoise = sump / j
100
100
101 return lnoise
101 return lnoise
102
102
103
103
104 class Beam:
104 class Beam:
105
105
106 def __init__(self):
106 def __init__(self):
107 self.codeList = []
107 self.codeList = []
108 self.azimuthList = []
108 self.azimuthList = []
109 self.zenithList = []
109 self.zenithList = []
110
110
111
111
112 class GenericData(object):
112 class GenericData(object):
113
113
114 flagNoData = True
114 flagNoData = True
115
115
116 def copy(self, inputObj=None):
116 def copy(self, inputObj=None):
117
117
118 if inputObj == None:
118 if inputObj == None:
119 return copy.deepcopy(self)
119 return copy.deepcopy(self)
120
120
121 for key in list(inputObj.__dict__.keys()):
121 for key in list(inputObj.__dict__.keys()):
122
122
123 attribute = inputObj.__dict__[key]
123 attribute = inputObj.__dict__[key]
124
124
125 # If this attribute is a tuple or list
125 # If this attribute is a tuple or list
126 if type(inputObj.__dict__[key]) in (tuple, list):
126 if type(inputObj.__dict__[key]) in (tuple, list):
127 self.__dict__[key] = attribute[:]
127 self.__dict__[key] = attribute[:]
128 continue
128 continue
129
129
130 # If this attribute is another object or instance
130 # If this attribute is another object or instance
131 if hasattr(attribute, '__dict__'):
131 if hasattr(attribute, '__dict__'):
132 self.__dict__[key] = attribute.copy()
132 self.__dict__[key] = attribute.copy()
133 continue
133 continue
134
134
135 self.__dict__[key] = inputObj.__dict__[key]
135 self.__dict__[key] = inputObj.__dict__[key]
136
136
137 def deepcopy(self):
137 def deepcopy(self):
138
138
139 return copy.deepcopy(self)
139 return copy.deepcopy(self)
140
140
141 def isEmpty(self):
141 def isEmpty(self):
142
142
143 return self.flagNoData
143 return self.flagNoData
144
144
145
145
146 class JROData(GenericData):
146 class JROData(GenericData):
147
147
148 # m_BasicHeader = BasicHeader()
148 # m_BasicHeader = BasicHeader()
149 # m_ProcessingHeader = ProcessingHeader()
149 # m_ProcessingHeader = ProcessingHeader()
150
150
151 systemHeaderObj = SystemHeader()
151 systemHeaderObj = SystemHeader()
152 radarControllerHeaderObj = RadarControllerHeader()
152 radarControllerHeaderObj = RadarControllerHeader()
153 # data = None
153 # data = None
154 type = None
154 type = None
155 datatype = None # dtype but in string
155 datatype = None # dtype but in string
156 # dtype = None
156 # dtype = None
157 # nChannels = None
157 # nChannels = None
158 # nHeights = None
158 # nHeights = None
159 nProfiles = None
159 nProfiles = None
160 heightList = None
160 heightList = None
161 channelList = None
161 channelList = None
162 flagDiscontinuousBlock = False
162 flagDiscontinuousBlock = False
163 useLocalTime = False
163 useLocalTime = False
164 utctime = None
164 utctime = None
165 timeZone = None
165 timeZone = None
166 dstFlag = None
166 dstFlag = None
167 errorCount = None
167 errorCount = None
168 blocksize = None
168 blocksize = None
169 # nCode = None
169 # nCode = None
170 # nBaud = None
170 # nBaud = None
171 # code = None
171 # code = None
172 flagDecodeData = False # asumo q la data no esta decodificada
172 flagDecodeData = False # asumo q la data no esta decodificada
173 flagDeflipData = False # asumo q la data no esta sin flip
173 flagDeflipData = False # asumo q la data no esta sin flip
174 flagShiftFFT = False
174 flagShiftFFT = False
175 # ippSeconds = None
175 # ippSeconds = None
176 # timeInterval = None
176 # timeInterval = None
177 nCohInt = None
177 nCohInt = None
178 # noise = None
178 # noise = None
179 windowOfFilter = 1
179 windowOfFilter = 1
180 # Speed of ligth
180 # Speed of ligth
181 C = 3e8
181 C = 3e8
182 frequency = 49.92e6
182 frequency = 49.92e6
183 realtime = False
183 realtime = False
184 beacon_heiIndexList = None
184 beacon_heiIndexList = None
185 last_block = None
185 last_block = None
186 blocknow = None
186 blocknow = None
187 azimuth = None
187 azimuth = None
188 zenith = None
188 zenith = None
189 beam = Beam()
189 beam = Beam()
190 profileIndex = None
190 profileIndex = None
191 error = None
191 error = None
192 data = None
192 data = None
193 nmodes = None
193 nmodes = None
194
194
195 def __str__(self):
195 def __str__(self):
196
196
197 return '{} - {}'.format(self.type, self.getDatatime())
197 return '{} - {}'.format(self.type, self.getDatatime())
198
198
199 def getNoise(self):
199 def getNoise(self):
200
200
201 raise NotImplementedError
201 raise NotImplementedError
202
202
203 def getNChannels(self):
203 def getNChannels(self):
204
204
205 return len(self.channelList)
205 return len(self.channelList)
206
206
207 def getChannelIndexList(self):
207 def getChannelIndexList(self):
208
208
209 return list(range(self.nChannels))
209 return list(range(self.nChannels))
210
210
211 def getNHeights(self):
211 def getNHeights(self):
212
212
213 return len(self.heightList)
213 return len(self.heightList)
214
214
215 def getHeiRange(self, extrapoints=0):
215 def getHeiRange(self, extrapoints=0):
216
216
217 heis = self.heightList
217 heis = self.heightList
218 # deltah = self.heightList[1] - self.heightList[0]
218 # deltah = self.heightList[1] - self.heightList[0]
219 #
219 #
220 # heis.append(self.heightList[-1])
220 # heis.append(self.heightList[-1])
221
221
222 return heis
222 return heis
223
223
224 def getDeltaH(self):
224 def getDeltaH(self):
225
225
226 delta = self.heightList[1] - self.heightList[0]
226 delta = self.heightList[1] - self.heightList[0]
227
227
228 return delta
228 return delta
229
229
230 def getltctime(self):
230 def getltctime(self):
231
231
232 if self.useLocalTime:
232 if self.useLocalTime:
233 return self.utctime - self.timeZone * 60
233 return self.utctime - self.timeZone * 60
234
234
235 return self.utctime
235 return self.utctime
236
236
237 def getDatatime(self):
237 def getDatatime(self):
238
238
239 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
239 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
240 return datatimeValue
240 return datatimeValue
241
241
242 def getTimeRange(self):
242 def getTimeRange(self):
243
243
244 datatime = []
244 datatime = []
245
245
246 datatime.append(self.ltctime)
246 datatime.append(self.ltctime)
247 datatime.append(self.ltctime + self.timeInterval + 1)
247 datatime.append(self.ltctime + self.timeInterval + 1)
248
248
249 datatime = numpy.array(datatime)
249 datatime = numpy.array(datatime)
250
250
251 return datatime
251 return datatime
252
252
253 def getFmaxTimeResponse(self):
253 def getFmaxTimeResponse(self):
254
254
255 period = (10**-6) * self.getDeltaH() / (0.15)
255 period = (10**-6) * self.getDeltaH() / (0.15)
256
256
257 PRF = 1. / (period * self.nCohInt)
257 PRF = 1. / (period * self.nCohInt)
258
258
259 fmax = PRF
259 fmax = PRF
260
260
261 return fmax
261 return fmax
262
262
263 def getFmax(self):
263 def getFmax(self):
264 PRF = 1. / (self.ippSeconds * self.nCohInt)
264 PRF = 1. / (self.ippSeconds * self.nCohInt)
265
265
266 fmax = PRF
266 fmax = PRF
267 return fmax
267 return fmax
268
268
269 def getVmax(self):
269 def getVmax(self):
270
270
271 _lambda = self.C / self.frequency
271 _lambda = self.C / self.frequency
272
272
273 vmax = self.getFmax() * _lambda / 2
273 vmax = self.getFmax() * _lambda / 2
274
274
275 return vmax
275 return vmax
276
276
277 def get_ippSeconds(self):
277 def get_ippSeconds(self):
278 '''
278 '''
279 '''
279 '''
280 return self.radarControllerHeaderObj.ippSeconds
280 return self.radarControllerHeaderObj.ippSeconds
281
281
282 def set_ippSeconds(self, ippSeconds):
282 def set_ippSeconds(self, ippSeconds):
283 '''
283 '''
284 '''
284 '''
285
285
286 self.radarControllerHeaderObj.ippSeconds = ippSeconds
286 self.radarControllerHeaderObj.ippSeconds = ippSeconds
287
287
288 return
288 return
289
289
290 def get_dtype(self):
290 def get_dtype(self):
291 '''
291 '''
292 '''
292 '''
293 return getNumpyDtype(self.datatype)
293 return getNumpyDtype(self.datatype)
294
294
295 def set_dtype(self, numpyDtype):
295 def set_dtype(self, numpyDtype):
296 '''
296 '''
297 '''
297 '''
298
298
299 self.datatype = getDataTypeCode(numpyDtype)
299 self.datatype = getDataTypeCode(numpyDtype)
300
300
301 def get_code(self):
301 def get_code(self):
302 '''
302 '''
303 '''
303 '''
304 return self.radarControllerHeaderObj.code
304 return self.radarControllerHeaderObj.code
305
305
306 def set_code(self, code):
306 def set_code(self, code):
307 '''
307 '''
308 '''
308 '''
309 self.radarControllerHeaderObj.code = code
309 self.radarControllerHeaderObj.code = code
310
310
311 return
311 return
312
312
313 def get_ncode(self):
313 def get_ncode(self):
314 '''
314 '''
315 '''
315 '''
316 return self.radarControllerHeaderObj.nCode
316 return self.radarControllerHeaderObj.nCode
317
317
318 def set_ncode(self, nCode):
318 def set_ncode(self, nCode):
319 '''
319 '''
320 '''
320 '''
321 self.radarControllerHeaderObj.nCode = nCode
321 self.radarControllerHeaderObj.nCode = nCode
322
322
323 return
323 return
324
324
325 def get_nbaud(self):
325 def get_nbaud(self):
326 '''
326 '''
327 '''
327 '''
328 return self.radarControllerHeaderObj.nBaud
328 return self.radarControllerHeaderObj.nBaud
329
329
330 def set_nbaud(self, nBaud):
330 def set_nbaud(self, nBaud):
331 '''
331 '''
332 '''
332 '''
333 self.radarControllerHeaderObj.nBaud = nBaud
333 self.radarControllerHeaderObj.nBaud = nBaud
334
334
335 return
335 return
336
336
337 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
337 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
338 channelIndexList = property(
338 channelIndexList = property(
339 getChannelIndexList, "I'm the 'channelIndexList' property.")
339 getChannelIndexList, "I'm the 'channelIndexList' property.")
340 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
340 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
341 #noise = property(getNoise, "I'm the 'nHeights' property.")
341 #noise = property(getNoise, "I'm the 'nHeights' property.")
342 datatime = property(getDatatime, "I'm the 'datatime' property")
342 datatime = property(getDatatime, "I'm the 'datatime' property")
343 ltctime = property(getltctime, "I'm the 'ltctime' property")
343 ltctime = property(getltctime, "I'm the 'ltctime' property")
344 ippSeconds = property(get_ippSeconds, set_ippSeconds)
344 ippSeconds = property(get_ippSeconds, set_ippSeconds)
345 dtype = property(get_dtype, set_dtype)
345 dtype = property(get_dtype, set_dtype)
346 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
346 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
347 code = property(get_code, set_code)
347 code = property(get_code, set_code)
348 nCode = property(get_ncode, set_ncode)
348 nCode = property(get_ncode, set_ncode)
349 nBaud = property(get_nbaud, set_nbaud)
349 nBaud = property(get_nbaud, set_nbaud)
350
350
351
351
352 class Voltage(JROData):
352 class Voltage(JROData):
353
353
354 # data es un numpy array de 2 dmensiones (canales, alturas)
354 # data es un numpy array de 2 dmensiones (canales, alturas)
355 data = None
355 data = None
356
356
357 def __init__(self):
357 def __init__(self):
358 '''
358 '''
359 Constructor
359 Constructor
360 '''
360 '''
361
361
362 self.useLocalTime = True
362 self.useLocalTime = True
363 self.radarControllerHeaderObj = RadarControllerHeader()
363 self.radarControllerHeaderObj = RadarControllerHeader()
364 self.systemHeaderObj = SystemHeader()
364 self.systemHeaderObj = SystemHeader()
365 self.type = "Voltage"
365 self.type = "Voltage"
366 self.data = None
366 self.data = None
367 # self.dtype = None
367 # self.dtype = None
368 # self.nChannels = 0
368 # self.nChannels = 0
369 # self.nHeights = 0
369 # self.nHeights = 0
370 self.nProfiles = None
370 self.nProfiles = None
371 self.heightList = None
371 self.heightList = None
372 self.channelList = None
372 self.channelList = None
373 # self.channelIndexList = None
373 # self.channelIndexList = None
374 self.flagNoData = True
374 self.flagNoData = True
375 self.flagDiscontinuousBlock = False
375 self.flagDiscontinuousBlock = False
376 self.utctime = None
376 self.utctime = None
377 self.timeZone = None
377 self.timeZone = None
378 self.dstFlag = None
378 self.dstFlag = None
379 self.errorCount = None
379 self.errorCount = None
380 self.nCohInt = None
380 self.nCohInt = None
381 self.blocksize = None
381 self.blocksize = None
382 self.flagDecodeData = False # asumo q la data no esta decodificada
382 self.flagDecodeData = False # asumo q la data no esta decodificada
383 self.flagDeflipData = False # asumo q la data no esta sin flip
383 self.flagDeflipData = False # asumo q la data no esta sin flip
384 self.flagShiftFFT = False
384 self.flagShiftFFT = False
385 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
385 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
386 self.profileIndex = 0
386 self.profileIndex = 0
387
387
388 def getNoisebyHildebrand(self, channel=None):
388 def getNoisebyHildebrand(self, channel=None):
389 """
389 """
390 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
390 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
391
391
392 Return:
392 Return:
393 noiselevel
393 noiselevel
394 """
394 """
395
395
396 if channel != None:
396 if channel != None:
397 data = self.data[channel]
397 data = self.data[channel]
398 nChannels = 1
398 nChannels = 1
399 else:
399 else:
400 data = self.data
400 data = self.data
401 nChannels = self.nChannels
401 nChannels = self.nChannels
402
402
403 noise = numpy.zeros(nChannels)
403 noise = numpy.zeros(nChannels)
404 power = data * numpy.conjugate(data)
404 power = data * numpy.conjugate(data)
405
405
406 for thisChannel in range(nChannels):
406 for thisChannel in range(nChannels):
407 if nChannels == 1:
407 if nChannels == 1:
408 daux = power[:].real
408 daux = power[:].real
409 else:
409 else:
410 daux = power[thisChannel, :].real
410 daux = power[thisChannel, :].real
411 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
411 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
412
412
413 return noise
413 return noise
414
414
415 def getNoise(self, type=1, channel=None):
415 def getNoise(self, type=1, channel=None):
416
416
417 if type == 1:
417 if type == 1:
418 noise = self.getNoisebyHildebrand(channel)
418 noise = self.getNoisebyHildebrand(channel)
419
419
420 return noise
420 return noise
421
421
422 def getPower(self, channel=None):
422 def getPower(self, channel=None):
423
423
424 if channel != None:
424 if channel != None:
425 data = self.data[channel]
425 data = self.data[channel]
426 else:
426 else:
427 data = self.data
427 data = self.data
428
428
429 power = data * numpy.conjugate(data)
429 power = data * numpy.conjugate(data)
430 powerdB = 10 * numpy.log10(power.real)
430 powerdB = 10 * numpy.log10(power.real)
431 powerdB = numpy.squeeze(powerdB)
431 powerdB = numpy.squeeze(powerdB)
432
432
433 return powerdB
433 return powerdB
434
434
435 def getTimeInterval(self):
435 def getTimeInterval(self):
436
436
437 timeInterval = self.ippSeconds * self.nCohInt
437 timeInterval = self.ippSeconds * self.nCohInt
438
438
439 return timeInterval
439 return timeInterval
440
440
441 noise = property(getNoise, "I'm the 'nHeights' property.")
441 noise = property(getNoise, "I'm the 'nHeights' property.")
442 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
442 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
443
443
444
444
445 class Spectra(JROData):
445 class Spectra(JROData):
446
446
447 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
447 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
448 data_spc = None
448 data_spc = None
449 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
449 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
450 data_cspc = None
450 data_cspc = None
451 # data dc es un numpy array de 2 dmensiones (canales, alturas)
451 # data dc es un numpy array de 2 dmensiones (canales, alturas)
452 data_dc = None
452 data_dc = None
453 # data power
453 # data power
454 data_pwr = None
454 data_pwr = None
455 nFFTPoints = None
455 nFFTPoints = None
456 # nPairs = None
456 # nPairs = None
457 pairsList = None
457 pairsList = None
458 nIncohInt = None
458 nIncohInt = None
459 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
459 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
460 nCohInt = None # se requiere para determinar el valor de timeInterval
460 nCohInt = None # se requiere para determinar el valor de timeInterval
461 ippFactor = None
461 ippFactor = None
462 profileIndex = 0
462 profileIndex = 0
463 plotting = "spectra"
463 plotting = "spectra"
464
464
465 def __init__(self):
465 def __init__(self):
466 '''
466 '''
467 Constructor
467 Constructor
468 '''
468 '''
469
469
470 self.useLocalTime = True
470 self.useLocalTime = True
471 self.radarControllerHeaderObj = RadarControllerHeader()
471 self.radarControllerHeaderObj = RadarControllerHeader()
472 self.systemHeaderObj = SystemHeader()
472 self.systemHeaderObj = SystemHeader()
473 self.type = "Spectra"
473 self.type = "Spectra"
474 # self.data = None
474 # self.data = None
475 # self.dtype = None
475 # self.dtype = None
476 # self.nChannels = 0
476 # self.nChannels = 0
477 # self.nHeights = 0
477 # self.nHeights = 0
478 self.nProfiles = None
478 self.nProfiles = None
479 self.heightList = None
479 self.heightList = None
480 self.channelList = None
480 self.channelList = None
481 # self.channelIndexList = None
481 # self.channelIndexList = None
482 self.pairsList = None
482 self.pairsList = None
483 self.flagNoData = True
483 self.flagNoData = True
484 self.flagDiscontinuousBlock = False
484 self.flagDiscontinuousBlock = False
485 self.utctime = None
485 self.utctime = None
486 self.nCohInt = None
486 self.nCohInt = None
487 self.nIncohInt = None
487 self.nIncohInt = None
488 self.blocksize = None
488 self.blocksize = None
489 self.nFFTPoints = None
489 self.nFFTPoints = None
490 self.wavelength = None
490 self.wavelength = None
491 self.flagDecodeData = False # asumo q la data no esta decodificada
491 self.flagDecodeData = False # asumo q la data no esta decodificada
492 self.flagDeflipData = False # asumo q la data no esta sin flip
492 self.flagDeflipData = False # asumo q la data no esta sin flip
493 self.flagShiftFFT = False
493 self.flagShiftFFT = False
494 self.ippFactor = 1
494 self.ippFactor = 1
495 #self.noise = None
495 #self.noise = None
496 self.beacon_heiIndexList = []
496 self.beacon_heiIndexList = []
497 self.noise_estimation = None
497 self.noise_estimation = None
498
498
499 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
499 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
500 """
500 """
501 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
501 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
502
502
503 Return:
503 Return:
504 noiselevel
504 noiselevel
505 """
505 """
506
506
507 noise = numpy.zeros(self.nChannels)
507 noise = numpy.zeros(self.nChannels)
508
508
509 for channel in range(self.nChannels):
509 for channel in range(self.nChannels):
510 daux = self.data_spc[channel,
510 daux = self.data_spc[channel,
511 xmin_index:xmax_index, ymin_index:ymax_index]
511 xmin_index:xmax_index, ymin_index:ymax_index]
512 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
512 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
513
513
514 return noise
514 return noise
515
515
516 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
516 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
517
517
518 if self.noise_estimation is not None:
518 if self.noise_estimation is not None:
519 # this was estimated by getNoise Operation defined in jroproc_spectra.py
519 # this was estimated by getNoise Operation defined in jroproc_spectra.py
520 return self.noise_estimation
520 return self.noise_estimation
521 else:
521 else:
522 noise = self.getNoisebyHildebrand(
522 noise = self.getNoisebyHildebrand(
523 xmin_index, xmax_index, ymin_index, ymax_index)
523 xmin_index, xmax_index, ymin_index, ymax_index)
524 return noise
524 return noise
525
525
526 def getFreqRangeTimeResponse(self, extrapoints=0):
526 def getFreqRangeTimeResponse(self, extrapoints=0):
527
527
528 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
528 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
529 freqrange = deltafreq * \
529 freqrange = deltafreq * \
530 (numpy.arange(self.nFFTPoints + extrapoints) -
530 (numpy.arange(self.nFFTPoints + extrapoints) -
531 self.nFFTPoints / 2.) - deltafreq / 2
531 self.nFFTPoints / 2.) - deltafreq / 2
532
532
533 return freqrange
533 return freqrange
534
534
535 def getAcfRange(self, extrapoints=0):
535 def getAcfRange(self, extrapoints=0):
536
536
537 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
537 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
538 freqrange = deltafreq * \
538 freqrange = deltafreq * \
539 (numpy.arange(self.nFFTPoints + extrapoints) -
539 (numpy.arange(self.nFFTPoints + extrapoints) -
540 self.nFFTPoints / 2.) - deltafreq / 2
540 self.nFFTPoints / 2.) - deltafreq / 2
541
541
542 return freqrange
542 return freqrange
543
543
544 def getFreqRange(self, extrapoints=0):
544 def getFreqRange(self, extrapoints=0):
545
545
546 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
546 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
547 freqrange = deltafreq * \
547 freqrange = deltafreq * \
548 (numpy.arange(self.nFFTPoints + extrapoints) -
548 (numpy.arange(self.nFFTPoints + extrapoints) -
549 self.nFFTPoints / 2.) - deltafreq / 2
549 self.nFFTPoints / 2.) - deltafreq / 2
550
550
551 return freqrange
551 return freqrange
552
552
553 def getVelRange(self, extrapoints=0):
553 def getVelRange(self, extrapoints=0):
554
554
555 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
555 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
556 velrange = deltav * (numpy.arange(self.nFFTPoints +
556 velrange = deltav * (numpy.arange(self.nFFTPoints +
557 extrapoints) - self.nFFTPoints / 2.)
557 extrapoints) - self.nFFTPoints / 2.)
558
558
559 if self.nmodes:
559 if self.nmodes:
560 return velrange/self.nmodes
560 return velrange/self.nmodes
561 else:
561 else:
562 return velrange
562 return velrange
563
563
564 def getNPairs(self):
564 def getNPairs(self):
565
565
566 return len(self.pairsList)
566 return len(self.pairsList)
567
567
568 def getPairsIndexList(self):
568 def getPairsIndexList(self):
569
569
570 return list(range(self.nPairs))
570 return list(range(self.nPairs))
571
571
572 def getNormFactor(self):
572 def getNormFactor(self):
573
573
574 pwcode = 1
574 pwcode = 1
575
575
576 if self.flagDecodeData:
576 if self.flagDecodeData:
577 pwcode = numpy.sum(self.code[0]**2)
577 pwcode = numpy.sum(self.code[0]**2)
578 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
578 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
579 normFactor = self.nProfiles * self.nIncohInt * \
579 normFactor = self.nProfiles * self.nIncohInt * \
580 self.nCohInt * pwcode * self.windowOfFilter
580 self.nCohInt * pwcode * self.windowOfFilter
581
581
582 return normFactor
582 return normFactor
583
583
584 def getFlagCspc(self):
584 def getFlagCspc(self):
585
585
586 if self.data_cspc is None:
586 if self.data_cspc is None:
587 return True
587 return True
588
588
589 return False
589 return False
590
590
591 def getFlagDc(self):
591 def getFlagDc(self):
592
592
593 if self.data_dc is None:
593 if self.data_dc is None:
594 return True
594 return True
595
595
596 return False
596 return False
597
597
598 def getTimeInterval(self):
598 def getTimeInterval(self):
599
599
600 timeInterval = self.ippSeconds * self.nCohInt * \
600 timeInterval = self.ippSeconds * self.nCohInt * \
601 self.nIncohInt * self.nProfiles * self.ippFactor
601 self.nIncohInt * self.nProfiles * self.ippFactor
602
602
603 return timeInterval
603 return timeInterval
604
604
605 def getPower(self):
605 def getPower(self):
606
606
607 factor = self.normFactor
607 factor = self.normFactor
608 z = self.data_spc / factor
608 z = self.data_spc / factor
609 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
609 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
610 avg = numpy.average(z, axis=1)
610 avg = numpy.average(z, axis=1)
611
611
612 return 10 * numpy.log10(avg)
612 return 10 * numpy.log10(avg)
613
613
614 def getCoherence(self, pairsList=None, phase=False):
614 def getCoherence(self, pairsList=None, phase=False):
615
615
616 z = []
616 z = []
617 if pairsList is None:
617 if pairsList is None:
618 pairsIndexList = self.pairsIndexList
618 pairsIndexList = self.pairsIndexList
619 else:
619 else:
620 pairsIndexList = []
620 pairsIndexList = []
621 for pair in pairsList:
621 for pair in pairsList:
622 if pair not in self.pairsList:
622 if pair not in self.pairsList:
623 raise ValueError("Pair %s is not in dataOut.pairsList" % (
623 raise ValueError("Pair %s is not in dataOut.pairsList" % (
624 pair))
624 pair))
625 pairsIndexList.append(self.pairsList.index(pair))
625 pairsIndexList.append(self.pairsList.index(pair))
626 for i in range(len(pairsIndexList)):
626 for i in range(len(pairsIndexList)):
627 pair = self.pairsList[pairsIndexList[i]]
627 pair = self.pairsList[pairsIndexList[i]]
628 ccf = numpy.average(
628 ccf = numpy.average(
629 self.data_cspc[pairsIndexList[i], :, :], axis=0)
629 self.data_cspc[pairsIndexList[i], :, :], axis=0)
630 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
630 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
631 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
631 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
632 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
632 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
633 if phase:
633 if phase:
634 data = numpy.arctan2(avgcoherenceComplex.imag,
634 data = numpy.arctan2(avgcoherenceComplex.imag,
635 avgcoherenceComplex.real) * 180 / numpy.pi
635 avgcoherenceComplex.real) * 180 / numpy.pi
636 else:
636 else:
637 data = numpy.abs(avgcoherenceComplex)
637 data = numpy.abs(avgcoherenceComplex)
638
638
639 z.append(data)
639 z.append(data)
640
640
641 return numpy.array(z)
641 return numpy.array(z)
642
642
643 def setValue(self, value):
643 def setValue(self, value):
644
644
645 print("This property should not be initialized")
645 print("This property should not be initialized")
646
646
647 return
647 return
648
648
649 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
649 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
650 pairsIndexList = property(
650 pairsIndexList = property(
651 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
651 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
652 normFactor = property(getNormFactor, setValue,
652 normFactor = property(getNormFactor, setValue,
653 "I'm the 'getNormFactor' property.")
653 "I'm the 'getNormFactor' property.")
654 flag_cspc = property(getFlagCspc, setValue)
654 flag_cspc = property(getFlagCspc, setValue)
655 flag_dc = property(getFlagDc, setValue)
655 flag_dc = property(getFlagDc, setValue)
656 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
656 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
657 timeInterval = property(getTimeInterval, setValue,
657 timeInterval = property(getTimeInterval, setValue,
658 "I'm the 'timeInterval' property")
658 "I'm the 'timeInterval' property")
659
659
660
660
661 class SpectraHeis(Spectra):
661 class SpectraHeis(Spectra):
662
662
663 data_spc = None
663 data_spc = None
664 data_cspc = None
664 data_cspc = None
665 data_dc = None
665 data_dc = None
666 nFFTPoints = None
666 nFFTPoints = None
667 # nPairs = None
667 # nPairs = None
668 pairsList = None
668 pairsList = None
669 nCohInt = None
669 nCohInt = None
670 nIncohInt = None
670 nIncohInt = None
671
671
672 def __init__(self):
672 def __init__(self):
673
673
674 self.radarControllerHeaderObj = RadarControllerHeader()
674 self.radarControllerHeaderObj = RadarControllerHeader()
675
675
676 self.systemHeaderObj = SystemHeader()
676 self.systemHeaderObj = SystemHeader()
677
677
678 self.type = "SpectraHeis"
678 self.type = "SpectraHeis"
679
679
680 # self.dtype = None
680 # self.dtype = None
681
681
682 # self.nChannels = 0
682 # self.nChannels = 0
683
683
684 # self.nHeights = 0
684 # self.nHeights = 0
685
685
686 self.nProfiles = None
686 self.nProfiles = None
687
687
688 self.heightList = None
688 self.heightList = None
689
689
690 self.channelList = None
690 self.channelList = None
691
691
692 # self.channelIndexList = None
692 # self.channelIndexList = None
693
693
694 self.flagNoData = True
694 self.flagNoData = True
695
695
696 self.flagDiscontinuousBlock = False
696 self.flagDiscontinuousBlock = False
697
697
698 # self.nPairs = 0
698 # self.nPairs = 0
699
699
700 self.utctime = None
700 self.utctime = None
701
701
702 self.blocksize = None
702 self.blocksize = None
703
703
704 self.profileIndex = 0
704 self.profileIndex = 0
705
705
706 self.nCohInt = 1
706 self.nCohInt = 1
707
707
708 self.nIncohInt = 1
708 self.nIncohInt = 1
709
709
710 def getNormFactor(self):
710 def getNormFactor(self):
711 pwcode = 1
711 pwcode = 1
712 if self.flagDecodeData:
712 if self.flagDecodeData:
713 pwcode = numpy.sum(self.code[0]**2)
713 pwcode = numpy.sum(self.code[0]**2)
714
714
715 normFactor = self.nIncohInt * self.nCohInt * pwcode
715 normFactor = self.nIncohInt * self.nCohInt * pwcode
716
716
717 return normFactor
717 return normFactor
718
718
719 def getTimeInterval(self):
719 def getTimeInterval(self):
720
720
721 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
721 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
722
722
723 return timeInterval
723 return timeInterval
724
724
725 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
725 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
726 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
726 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
727
727
728
728
729 class Fits(JROData):
729 class Fits(JROData):
730
730
731 heightList = None
731 heightList = None
732 channelList = None
732 channelList = None
733 flagNoData = True
733 flagNoData = True
734 flagDiscontinuousBlock = False
734 flagDiscontinuousBlock = False
735 useLocalTime = False
735 useLocalTime = False
736 utctime = None
736 utctime = None
737 timeZone = None
737 timeZone = None
738 # ippSeconds = None
738 # ippSeconds = None
739 # timeInterval = None
739 # timeInterval = None
740 nCohInt = None
740 nCohInt = None
741 nIncohInt = None
741 nIncohInt = None
742 noise = None
742 noise = None
743 windowOfFilter = 1
743 windowOfFilter = 1
744 # Speed of ligth
744 # Speed of ligth
745 C = 3e8
745 C = 3e8
746 frequency = 49.92e6
746 frequency = 49.92e6
747 realtime = False
747 realtime = False
748
748
749 def __init__(self):
749 def __init__(self):
750
750
751 self.type = "Fits"
751 self.type = "Fits"
752
752
753 self.nProfiles = None
753 self.nProfiles = None
754
754
755 self.heightList = None
755 self.heightList = None
756
756
757 self.channelList = None
757 self.channelList = None
758
758
759 # self.channelIndexList = None
759 # self.channelIndexList = None
760
760
761 self.flagNoData = True
761 self.flagNoData = True
762
762
763 self.utctime = None
763 self.utctime = None
764
764
765 self.nCohInt = 1
765 self.nCohInt = 1
766
766
767 self.nIncohInt = 1
767 self.nIncohInt = 1
768
768
769 self.useLocalTime = True
769 self.useLocalTime = True
770
770
771 self.profileIndex = 0
771 self.profileIndex = 0
772
772
773 # self.utctime = None
773 # self.utctime = None
774 # self.timeZone = None
774 # self.timeZone = None
775 # self.ltctime = None
775 # self.ltctime = None
776 # self.timeInterval = None
776 # self.timeInterval = None
777 # self.header = None
777 # self.header = None
778 # self.data_header = None
778 # self.data_header = None
779 # self.data = None
779 # self.data = None
780 # self.datatime = None
780 # self.datatime = None
781 # self.flagNoData = False
781 # self.flagNoData = False
782 # self.expName = ''
782 # self.expName = ''
783 # self.nChannels = None
783 # self.nChannels = None
784 # self.nSamples = None
784 # self.nSamples = None
785 # self.dataBlocksPerFile = None
785 # self.dataBlocksPerFile = None
786 # self.comments = ''
786 # self.comments = ''
787 #
787 #
788
788
789 def getltctime(self):
789 def getltctime(self):
790
790
791 if self.useLocalTime:
791 if self.useLocalTime:
792 return self.utctime - self.timeZone * 60
792 return self.utctime - self.timeZone * 60
793
793
794 return self.utctime
794 return self.utctime
795
795
796 def getDatatime(self):
796 def getDatatime(self):
797
797
798 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
798 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
799 return datatime
799 return datatime
800
800
801 def getTimeRange(self):
801 def getTimeRange(self):
802
802
803 datatime = []
803 datatime = []
804
804
805 datatime.append(self.ltctime)
805 datatime.append(self.ltctime)
806 datatime.append(self.ltctime + self.timeInterval)
806 datatime.append(self.ltctime + self.timeInterval)
807
807
808 datatime = numpy.array(datatime)
808 datatime = numpy.array(datatime)
809
809
810 return datatime
810 return datatime
811
811
812 def getHeiRange(self):
812 def getHeiRange(self):
813
813
814 heis = self.heightList
814 heis = self.heightList
815
815
816 return heis
816 return heis
817
817
818 def getNHeights(self):
818 def getNHeights(self):
819
819
820 return len(self.heightList)
820 return len(self.heightList)
821
821
822 def getNChannels(self):
822 def getNChannels(self):
823
823
824 return len(self.channelList)
824 return len(self.channelList)
825
825
826 def getChannelIndexList(self):
826 def getChannelIndexList(self):
827
827
828 return list(range(self.nChannels))
828 return list(range(self.nChannels))
829
829
830 def getNoise(self, type=1):
830 def getNoise(self, type=1):
831
831
832 #noise = numpy.zeros(self.nChannels)
832 #noise = numpy.zeros(self.nChannels)
833
833
834 if type == 1:
834 if type == 1:
835 noise = self.getNoisebyHildebrand()
835 noise = self.getNoisebyHildebrand()
836
836
837 if type == 2:
837 if type == 2:
838 noise = self.getNoisebySort()
838 noise = self.getNoisebySort()
839
839
840 if type == 3:
840 if type == 3:
841 noise = self.getNoisebyWindow()
841 noise = self.getNoisebyWindow()
842
842
843 return noise
843 return noise
844
844
845 def getTimeInterval(self):
845 def getTimeInterval(self):
846
846
847 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
847 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
848
848
849 return timeInterval
849 return timeInterval
850
850
851 def get_ippSeconds(self):
851 def get_ippSeconds(self):
852 '''
852 '''
853 '''
853 '''
854 return self.ipp_sec
854 return self.ipp_sec
855
855
856
856
857 datatime = property(getDatatime, "I'm the 'datatime' property")
857 datatime = property(getDatatime, "I'm the 'datatime' property")
858 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
858 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
859 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
859 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
860 channelIndexList = property(
860 channelIndexList = property(
861 getChannelIndexList, "I'm the 'channelIndexList' property.")
861 getChannelIndexList, "I'm the 'channelIndexList' property.")
862 noise = property(getNoise, "I'm the 'nHeights' property.")
862 noise = property(getNoise, "I'm the 'nHeights' property.")
863
863
864 ltctime = property(getltctime, "I'm the 'ltctime' property")
864 ltctime = property(getltctime, "I'm the 'ltctime' property")
865 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
865 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
866 ippSeconds = property(get_ippSeconds, '')
866 ippSeconds = property(get_ippSeconds, '')
867
867
868 class Correlation(JROData):
868 class Correlation(JROData):
869
869
870 noise = None
870 noise = None
871 SNR = None
871 SNR = None
872 #--------------------------------------------------
872 #--------------------------------------------------
873 mode = None
873 mode = None
874 split = False
874 split = False
875 data_cf = None
875 data_cf = None
876 lags = None
876 lags = None
877 lagRange = None
877 lagRange = None
878 pairsList = None
878 pairsList = None
879 normFactor = None
879 normFactor = None
880 #--------------------------------------------------
880 #--------------------------------------------------
881 # calculateVelocity = None
881 # calculateVelocity = None
882 nLags = None
882 nLags = None
883 nPairs = None
883 nPairs = None
884 nAvg = None
884 nAvg = None
885
885
886 def __init__(self):
886 def __init__(self):
887 '''
887 '''
888 Constructor
888 Constructor
889 '''
889 '''
890 self.radarControllerHeaderObj = RadarControllerHeader()
890 self.radarControllerHeaderObj = RadarControllerHeader()
891
891
892 self.systemHeaderObj = SystemHeader()
892 self.systemHeaderObj = SystemHeader()
893
893
894 self.type = "Correlation"
894 self.type = "Correlation"
895
895
896 self.data = None
896 self.data = None
897
897
898 self.dtype = None
898 self.dtype = None
899
899
900 self.nProfiles = None
900 self.nProfiles = None
901
901
902 self.heightList = None
902 self.heightList = None
903
903
904 self.channelList = None
904 self.channelList = None
905
905
906 self.flagNoData = True
906 self.flagNoData = True
907
907
908 self.flagDiscontinuousBlock = False
908 self.flagDiscontinuousBlock = False
909
909
910 self.utctime = None
910 self.utctime = None
911
911
912 self.timeZone = None
912 self.timeZone = None
913
913
914 self.dstFlag = None
914 self.dstFlag = None
915
915
916 self.errorCount = None
916 self.errorCount = None
917
917
918 self.blocksize = None
918 self.blocksize = None
919
919
920 self.flagDecodeData = False # asumo q la data no esta decodificada
920 self.flagDecodeData = False # asumo q la data no esta decodificada
921
921
922 self.flagDeflipData = False # asumo q la data no esta sin flip
922 self.flagDeflipData = False # asumo q la data no esta sin flip
923
923
924 self.pairsList = None
924 self.pairsList = None
925
925
926 self.nPoints = None
926 self.nPoints = None
927
927
928 def getPairsList(self):
928 def getPairsList(self):
929
929
930 return self.pairsList
930 return self.pairsList
931
931
932 def getNoise(self, mode=2):
932 def getNoise(self, mode=2):
933
933
934 indR = numpy.where(self.lagR == 0)[0][0]
934 indR = numpy.where(self.lagR == 0)[0][0]
935 indT = numpy.where(self.lagT == 0)[0][0]
935 indT = numpy.where(self.lagT == 0)[0][0]
936
936
937 jspectra0 = self.data_corr[:, :, indR, :]
937 jspectra0 = self.data_corr[:, :, indR, :]
938 jspectra = copy.copy(jspectra0)
938 jspectra = copy.copy(jspectra0)
939
939
940 num_chan = jspectra.shape[0]
940 num_chan = jspectra.shape[0]
941 num_hei = jspectra.shape[2]
941 num_hei = jspectra.shape[2]
942
942
943 freq_dc = jspectra.shape[1] / 2
943 freq_dc = jspectra.shape[1] / 2
944 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
944 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
945
945
946 if ind_vel[0] < 0:
946 if ind_vel[0] < 0:
947 ind_vel[list(range(0, 1))] = ind_vel[list(
947 ind_vel[list(range(0, 1))] = ind_vel[list(
948 range(0, 1))] + self.num_prof
948 range(0, 1))] + self.num_prof
949
949
950 if mode == 1:
950 if mode == 1:
951 jspectra[:, freq_dc, :] = (
951 jspectra[:, freq_dc, :] = (
952 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
952 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
953
953
954 if mode == 2:
954 if mode == 2:
955
955
956 vel = numpy.array([-2, -1, 1, 2])
956 vel = numpy.array([-2, -1, 1, 2])
957 xx = numpy.zeros([4, 4])
957 xx = numpy.zeros([4, 4])
958
958
959 for fil in range(4):
959 for fil in range(4):
960 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
960 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
961
961
962 xx_inv = numpy.linalg.inv(xx)
962 xx_inv = numpy.linalg.inv(xx)
963 xx_aux = xx_inv[0, :]
963 xx_aux = xx_inv[0, :]
964
964
965 for ich in range(num_chan):
965 for ich in range(num_chan):
966 yy = jspectra[ich, ind_vel, :]
966 yy = jspectra[ich, ind_vel, :]
967 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
967 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
968
968
969 junkid = jspectra[ich, freq_dc, :] <= 0
969 junkid = jspectra[ich, freq_dc, :] <= 0
970 cjunkid = sum(junkid)
970 cjunkid = sum(junkid)
971
971
972 if cjunkid.any():
972 if cjunkid.any():
973 jspectra[ich, freq_dc, junkid.nonzero()] = (
973 jspectra[ich, freq_dc, junkid.nonzero()] = (
974 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
974 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
975
975
976 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
976 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
977
977
978 return noise
978 return noise
979
979
980 def getTimeInterval(self):
980 def getTimeInterval(self):
981
981
982 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
982 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
983
983
984 return timeInterval
984 return timeInterval
985
985
986 def splitFunctions(self):
986 def splitFunctions(self):
987
987
988 pairsList = self.pairsList
988 pairsList = self.pairsList
989 ccf_pairs = []
989 ccf_pairs = []
990 acf_pairs = []
990 acf_pairs = []
991 ccf_ind = []
991 ccf_ind = []
992 acf_ind = []
992 acf_ind = []
993 for l in range(len(pairsList)):
993 for l in range(len(pairsList)):
994 chan0 = pairsList[l][0]
994 chan0 = pairsList[l][0]
995 chan1 = pairsList[l][1]
995 chan1 = pairsList[l][1]
996
996
997 # Obteniendo pares de Autocorrelacion
997 # Obteniendo pares de Autocorrelacion
998 if chan0 == chan1:
998 if chan0 == chan1:
999 acf_pairs.append(chan0)
999 acf_pairs.append(chan0)
1000 acf_ind.append(l)
1000 acf_ind.append(l)
1001 else:
1001 else:
1002 ccf_pairs.append(pairsList[l])
1002 ccf_pairs.append(pairsList[l])
1003 ccf_ind.append(l)
1003 ccf_ind.append(l)
1004
1004
1005 data_acf = self.data_cf[acf_ind]
1005 data_acf = self.data_cf[acf_ind]
1006 data_ccf = self.data_cf[ccf_ind]
1006 data_ccf = self.data_cf[ccf_ind]
1007
1007
1008 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1008 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1009
1009
1010 def getNormFactor(self):
1010 def getNormFactor(self):
1011 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1011 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1012 acf_pairs = numpy.array(acf_pairs)
1012 acf_pairs = numpy.array(acf_pairs)
1013 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1013 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1014
1014
1015 for p in range(self.nPairs):
1015 for p in range(self.nPairs):
1016 pair = self.pairsList[p]
1016 pair = self.pairsList[p]
1017
1017
1018 ch0 = pair[0]
1018 ch0 = pair[0]
1019 ch1 = pair[1]
1019 ch1 = pair[1]
1020
1020
1021 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1021 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1022 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1022 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1023 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1023 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1024
1024
1025 return normFactor
1025 return normFactor
1026
1026
1027 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1027 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1028 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1028 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1029
1029
1030
1030
1031 class Parameters(Spectra):
1031 class Parameters(Spectra):
1032
1032
1033 experimentInfo = None # Information about the experiment
1033 experimentInfo = None # Information about the experiment
1034 # Information from previous data
1034 # Information from previous data
1035 inputUnit = None # Type of data to be processed
1035 inputUnit = None # Type of data to be processed
1036 operation = None # Type of operation to parametrize
1036 operation = None # Type of operation to parametrize
1037 # normFactor = None #Normalization Factor
1037 # normFactor = None #Normalization Factor
1038 groupList = None # List of Pairs, Groups, etc
1038 groupList = None # List of Pairs, Groups, etc
1039 # Parameters
1039 # Parameters
1040 data_param = None # Parameters obtained
1040 data_param = None # Parameters obtained
1041 data_pre = None # Data Pre Parametrization
1041 data_pre = None # Data Pre Parametrization
1042 data_SNR = None # Signal to Noise Ratio
1042 data_SNR = None # Signal to Noise Ratio
1043 # heightRange = None #Heights
1043 # heightRange = None #Heights
1044 abscissaList = None # Abscissa, can be velocities, lags or time
1044 abscissaList = None # Abscissa, can be velocities, lags or time
1045 # noise = None #Noise Potency
1045 # noise = None #Noise Potency
1046 utctimeInit = None # Initial UTC time
1046 utctimeInit = None # Initial UTC time
1047 paramInterval = None # Time interval to calculate Parameters in seconds
1047 paramInterval = None # Time interval to calculate Parameters in seconds
1048 useLocalTime = True
1048 useLocalTime = True
1049 # Fitting
1049 # Fitting
1050 data_error = None # Error of the estimation
1050 data_error = None # Error of the estimation
1051 constants = None
1051 constants = None
1052 library = None
1052 library = None
1053 # Output signal
1053 # Output signal
1054 outputInterval = None # Time interval to calculate output signal in seconds
1054 outputInterval = None # Time interval to calculate output signal in seconds
1055 data_output = None # Out signal
1055 data_output = None # Out signal
1056 nAvg = None
1056 nAvg = None
1057 noise_estimation = None
1057 noise_estimation = None
1058 GauSPC = None # Fit gaussian SPC
1058 GauSPC = None # Fit gaussian SPC
1059
1059
1060 def __init__(self):
1060 def __init__(self):
1061 '''
1061 '''
1062 Constructor
1062 Constructor
1063 '''
1063 '''
1064 self.radarControllerHeaderObj = RadarControllerHeader()
1064 self.radarControllerHeaderObj = RadarControllerHeader()
1065
1065
1066 self.systemHeaderObj = SystemHeader()
1066 self.systemHeaderObj = SystemHeader()
1067
1067
1068 self.type = "Parameters"
1068 self.type = "Parameters"
1069
1069
1070 def getTimeRange1(self, interval):
1070 def getTimeRange1(self, interval):
1071
1071
1072 datatime = []
1072 datatime = []
1073
1073
1074 if self.useLocalTime:
1074 if self.useLocalTime:
1075 time1 = self.utctimeInit - self.timeZone * 60
1075 time1 = self.utctimeInit - self.timeZone * 60
1076 else:
1076 else:
1077 time1 = self.utctimeInit
1077 time1 = self.utctimeInit
1078
1078
1079 datatime.append(time1)
1079 datatime.append(time1)
1080 datatime.append(time1 + interval)
1080 datatime.append(time1 + interval)
1081 datatime = numpy.array(datatime)
1081 datatime = numpy.array(datatime)
1082
1082
1083 return datatime
1083 return datatime
1084
1084
1085 def getTimeInterval(self):
1085 def getTimeInterval(self):
1086
1086
1087 if hasattr(self, 'timeInterval1'):
1087 if hasattr(self, 'timeInterval1'):
1088 return self.timeInterval1
1088 return self.timeInterval1
1089 else:
1089 else:
1090 return self.paramInterval
1090 return self.paramInterval
1091
1091
1092 def setValue(self, value):
1092 def setValue(self, value):
1093
1093
1094 print("This property should not be initialized")
1094 print("This property should not be initialized")
1095
1095
1096 return
1096 return
1097
1097
1098 def getNoise(self):
1098 def getNoise(self):
1099
1099
1100 return self.spc_noise
1100 return self.spc_noise
1101
1101
1102 timeInterval = property(getTimeInterval)
1102 timeInterval = property(getTimeInterval)
1103 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1103 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1104
1104
1105
1105
1106 class PlotterData(object):
1106 class PlotterData(object):
1107 '''
1107 '''
1108 Object to hold data to be plotted
1108 Object to hold data to be plotted
1109 '''
1109 '''
1110
1110
1111 MAXNUMX = 100
1111 MAXNUMX = 100
1112 MAXNUMY = 100
1112 MAXNUMY = 100
1113
1113
1114 def __init__(self, code, throttle_value, exp_code, buffering=True):
1114 def __init__(self, code, throttle_value, exp_code, buffering=True):
1115
1115
1116 self.throttle = throttle_value
1116 self.throttle = throttle_value
1117 self.exp_code = exp_code
1117 self.exp_code = exp_code
1118 self.buffering = buffering
1118 self.buffering = buffering
1119 self.ready = False
1119 self.ready = False
1120 self.localtime = False
1120 self.localtime = False
1121 self.data = {}
1121 self.data = {}
1122 self.meta = {}
1122 self.meta = {}
1123 self.__times = []
1123 self.__times = []
1124 self.__heights = []
1124 self.__heights = []
1125
1125
1126 if 'snr' in code:
1126 if 'snr' in code:
1127 self.plottypes = ['snr']
1127 self.plottypes = ['snr']
1128 elif code == 'spc':
1128 elif code == 'spc':
1129 self.plottypes = ['spc', 'noise', 'rti']
1129 self.plottypes = ['spc', 'noise', 'rti']
1130 elif code == 'rti':
1130 elif code == 'rti':
1131 self.plottypes = ['noise', 'rti']
1131 self.plottypes = ['noise', 'rti']
1132 else:
1132 else:
1133 self.plottypes = [code]
1133 self.plottypes = [code]
1134
1134
1135 for plot in self.plottypes:
1135 for plot in self.plottypes:
1136 self.data[plot] = {}
1136 self.data[plot] = {}
1137
1137
1138 def __str__(self):
1138 def __str__(self):
1139 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1139 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1140 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
1140 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
1141
1141
1142 def __len__(self):
1142 def __len__(self):
1143 return len(self.__times)
1143 return len(self.__times)
1144
1144
1145 def __getitem__(self, key):
1145 def __getitem__(self, key):
1146
1146
1147 if key not in self.data:
1147 if key not in self.data:
1148 raise KeyError(log.error('Missing key: {}'.format(key)))
1148 raise KeyError(log.error('Missing key: {}'.format(key)))
1149 if 'spc' in key or not self.buffering:
1149 if 'spc' in key or not self.buffering:
1150 ret = self.data[key]
1150 ret = self.data[key]
1151 elif 'scope' in key:
1152 ret = numpy.array(self.data[key][float(self.tm)])
1151 else:
1153 else:
1152 ret = numpy.array([self.data[key][x] for x in self.times])
1154 ret = numpy.array([self.data[key][x] for x in self.times])
1153 if ret.ndim > 1:
1155 if ret.ndim > 1:
1154 ret = numpy.swapaxes(ret, 0, 1)
1156 ret = numpy.swapaxes(ret, 0, 1)
1155 return ret
1157 return ret
1156
1158
1157 def __contains__(self, key):
1159 def __contains__(self, key):
1158 return key in self.data
1160 return key in self.data
1159
1161
1160 def setup(self):
1162 def setup(self):
1161 '''
1163 '''
1162 Configure object
1164 Configure object
1163 '''
1165 '''
1164
1166
1165 self.type = ''
1167 self.type = ''
1166 self.ready = False
1168 self.ready = False
1167 self.data = {}
1169 self.data = {}
1168 self.__times = []
1170 self.__times = []
1169 self.__heights = []
1171 self.__heights = []
1170 self.__all_heights = set()
1172 self.__all_heights = set()
1171 for plot in self.plottypes:
1173 for plot in self.plottypes:
1172 if 'snr' in plot:
1174 if 'snr' in plot:
1173 plot = 'snr'
1175 plot = 'snr'
1174 self.data[plot] = {}
1176 self.data[plot] = {}
1175
1177
1176 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data:
1178 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data:
1177 self.data['noise'] = {}
1179 self.data['noise'] = {}
1178 if 'noise' not in self.plottypes:
1180 if 'noise' not in self.plottypes:
1179 self.plottypes.append('noise')
1181 self.plottypes.append('noise')
1180
1182
1181 def shape(self, key):
1183 def shape(self, key):
1182 '''
1184 '''
1183 Get the shape of the one-element data for the given key
1185 Get the shape of the one-element data for the given key
1184 '''
1186 '''
1185
1187
1186 if len(self.data[key]):
1188 if len(self.data[key]):
1187 if 'spc' in key or not self.buffering:
1189 if 'spc' in key or not self.buffering:
1188 return self.data[key].shape
1190 return self.data[key].shape
1189 return self.data[key][self.__times[0]].shape
1191 return self.data[key][self.__times[0]].shape
1190 return (0,)
1192 return (0,)
1191
1193
1192 def update(self, dataOut, tm):
1194 def update(self, dataOut, tm):
1193 '''
1195 '''
1194 Update data object with new dataOut
1196 Update data object with new dataOut
1195 '''
1197 '''
1196
1198
1197 if tm in self.__times:
1199 if tm in self.__times:
1198 return
1200 return
1199
1201 self.profileIndex = dataOut.profileIndex
1202 self.tm = tm
1200 self.type = dataOut.type
1203 self.type = dataOut.type
1201 self.parameters = getattr(dataOut, 'parameters', [])
1204 self.parameters = getattr(dataOut, 'parameters', [])
1202 if hasattr(dataOut, 'pairsList'):
1205 if hasattr(dataOut, 'pairsList'):
1203 self.pairs = dataOut.pairsList
1206 self.pairs = dataOut.pairsList
1204 if hasattr(dataOut, 'meta'):
1207 if hasattr(dataOut, 'meta'):
1205 self.meta = dataOut.meta
1208 self.meta = dataOut.meta
1206 self.channels = dataOut.channelList
1209 self.channels = dataOut.channelList
1207 self.interval = dataOut.getTimeInterval()
1210 self.interval = dataOut.getTimeInterval()
1208 self.localtime = dataOut.useLocalTime
1211 self.localtime = dataOut.useLocalTime
1209 if 'spc' in self.plottypes or 'cspc' in self.plottypes:
1212 if 'spc' in self.plottypes or 'cspc' in self.plottypes:
1210 self.xrange = (dataOut.getFreqRange(1)/1000.,
1213 self.xrange = (dataOut.getFreqRange(1)/1000.,
1211 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1214 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1212 self.factor = dataOut.normFactor
1215 self.factor = dataOut.normFactor
1213 self.__heights.append(dataOut.heightList)
1216 self.__heights.append(dataOut.heightList)
1214 self.__all_heights.update(dataOut.heightList)
1217 self.__all_heights.update(dataOut.heightList)
1215 self.__times.append(tm)
1218 self.__times.append(tm)
1216
1219
1217 for plot in self.plottypes:
1220 for plot in self.plottypes:
1218 if plot == 'spc':
1221 if plot == 'spc':
1219 z = dataOut.data_spc/dataOut.normFactor
1222 z = dataOut.data_spc/dataOut.normFactor
1220 buffer = 10*numpy.log10(z)
1223 buffer = 10*numpy.log10(z)
1221 if plot == 'cspc':
1224 if plot == 'cspc':
1222 z = dataOut.data_spc/dataOut.normFactor
1225 z = dataOut.data_spc/dataOut.normFactor
1223 buffer = (dataOut.data_spc, dataOut.data_cspc)
1226 buffer = (dataOut.data_spc, dataOut.data_cspc)
1224 if plot == 'noise':
1227 if plot == 'noise':
1225 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1228 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1226 if plot == 'rti':
1229 if plot == 'rti':
1227 buffer = dataOut.getPower()
1230 buffer = dataOut.getPower()
1228 if plot == 'snr_db':
1231 if plot == 'snr_db':
1229 buffer = dataOut.data_SNR
1232 buffer = dataOut.data_SNR
1230 if plot == 'snr':
1233 if plot == 'snr':
1231 buffer = 10*numpy.log10(dataOut.data_SNR)
1234 buffer = 10*numpy.log10(dataOut.data_SNR)
1232 if plot == 'dop':
1235 if plot == 'dop':
1233 buffer = 10*numpy.log10(dataOut.data_DOP)
1236 buffer = 10*numpy.log10(dataOut.data_DOP)
1234 if plot == 'mean':
1237 if plot == 'mean':
1235 buffer = dataOut.data_MEAN
1238 buffer = dataOut.data_MEAN
1236 if plot == 'std':
1239 if plot == 'std':
1237 buffer = dataOut.data_STD
1240 buffer = dataOut.data_STD
1238 if plot == 'coh':
1241 if plot == 'coh':
1239 buffer = dataOut.getCoherence()
1242 buffer = dataOut.getCoherence()
1240 if plot == 'phase':
1243 if plot == 'phase':
1241 buffer = dataOut.getCoherence(phase=True)
1244 buffer = dataOut.getCoherence(phase=True)
1242 if plot == 'output':
1245 if plot == 'output':
1243 buffer = dataOut.data_output
1246 buffer = dataOut.data_output
1244 if plot == 'param':
1247 if plot == 'param':
1245 buffer = dataOut.data_param
1248 buffer = dataOut.data_param
1246
1249 if plot == 'scope':
1250 buffer = dataOut.data
1251 self.flagDataAsBlock = dataOut.flagDataAsBlock
1252 self.nProfiles = dataOut.nProfiles
1253
1247 if plot == 'spc':
1254 if plot == 'spc':
1248 self.data[plot] = buffer
1255 self.data[plot] = buffer
1249 elif plot == 'cspc':
1256 elif plot == 'cspc':
1250 self.data['spc'] = buffer[0]
1257 self.data['spc'] = buffer[0]
1251 self.data['cspc'] = buffer[1]
1258 self.data['cspc'] = buffer[1]
1252 else:
1259 else:
1253 if self.buffering:
1260 if self.buffering:
1254 self.data[plot][tm] = buffer
1261 self.data[plot][tm] = buffer
1255 else:
1262 else:
1256 self.data[plot] = buffer
1263 self.data[plot] = buffer
1257
1264
1258 def normalize_heights(self):
1265 def normalize_heights(self):
1259 '''
1266 '''
1260 Ensure same-dimension of the data for different heighList
1267 Ensure same-dimension of the data for different heighList
1261 '''
1268 '''
1262
1269
1263 H = numpy.array(list(self.__all_heights))
1270 H = numpy.array(list(self.__all_heights))
1264 H.sort()
1271 H.sort()
1265 for key in self.data:
1272 for key in self.data:
1266 shape = self.shape(key)[:-1] + H.shape
1273 shape = self.shape(key)[:-1] + H.shape
1267 for tm, obj in list(self.data[key].items()):
1274 for tm, obj in list(self.data[key].items()):
1268 h = self.__heights[self.__times.index(tm)]
1275 h = self.__heights[self.__times.index(tm)]
1269 if H.size == h.size:
1276 if H.size == h.size:
1270 continue
1277 continue
1271 index = numpy.where(numpy.in1d(H, h))[0]
1278 index = numpy.where(numpy.in1d(H, h))[0]
1272 dummy = numpy.zeros(shape) + numpy.nan
1279 dummy = numpy.zeros(shape) + numpy.nan
1273 if len(shape) == 2:
1280 if len(shape) == 2:
1274 dummy[:, index] = obj
1281 dummy[:, index] = obj
1275 else:
1282 else:
1276 dummy[index] = obj
1283 dummy[index] = obj
1277 self.data[key][tm] = dummy
1284 self.data[key][tm] = dummy
1278
1285
1279 self.__heights = [H for tm in self.__times]
1286 self.__heights = [H for tm in self.__times]
1280
1287
1281 def jsonify(self, decimate=False):
1288 def jsonify(self, decimate=False):
1282 '''
1289 '''
1283 Convert data to json
1290 Convert data to json
1284 '''
1291 '''
1285
1292
1286 data = {}
1293 data = {}
1287 tm = self.times[-1]
1294 tm = self.times[-1]
1288 dy = int(self.heights.size/self.MAXNUMY) + 1
1295 dy = int(self.heights.size/self.MAXNUMY) + 1
1289 for key in self.data:
1296 for key in self.data:
1290 if key in ('spc', 'cspc') or not self.buffering:
1297 if key in ('spc', 'cspc') or not self.buffering:
1291 dx = int(self.data[key].shape[1]/self.MAXNUMX) + 1
1298 dx = int(self.data[key].shape[1]/self.MAXNUMX) + 1
1292 data[key] = self.roundFloats(
1299 data[key] = self.roundFloats(
1293 self.data[key][::, ::dx, ::dy].tolist())
1300 self.data[key][::, ::dx, ::dy].tolist())
1294 else:
1301 else:
1295 data[key] = self.roundFloats(self.data[key][tm].tolist())
1302 data[key] = self.roundFloats(self.data[key][tm].tolist())
1296
1303
1297 ret = {'data': data}
1304 ret = {'data': data}
1298 ret['exp_code'] = self.exp_code
1305 ret['exp_code'] = self.exp_code
1299 ret['time'] = float(tm)
1306 ret['time'] = float(tm)
1300 ret['interval'] = float(self.interval)
1307 ret['interval'] = float(self.interval)
1301 ret['localtime'] = self.localtime
1308 ret['localtime'] = self.localtime
1302 ret['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1309 ret['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1303 if 'spc' in self.data or 'cspc' in self.data:
1310 if 'spc' in self.data or 'cspc' in self.data:
1304 ret['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1311 ret['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1305 else:
1312 else:
1306 ret['xrange'] = []
1313 ret['xrange'] = []
1307 if hasattr(self, 'pairs'):
1314 if hasattr(self, 'pairs'):
1308 ret['pairs'] = [(int(p[0]), int(p[1])) for p in self.pairs]
1315 ret['pairs'] = [(int(p[0]), int(p[1])) for p in self.pairs]
1309 else:
1316 else:
1310 ret['pairs'] = []
1317 ret['pairs'] = []
1311
1318
1312 for key, value in list(self.meta.items()):
1319 for key, value in list(self.meta.items()):
1313 ret[key] = value
1320 ret[key] = value
1314
1321
1315 return json.dumps(ret)
1322 return json.dumps(ret)
1316
1323
1317 @property
1324 @property
1318 def times(self):
1325 def times(self):
1319 '''
1326 '''
1320 Return the list of times of the current data
1327 Return the list of times of the current data
1321 '''
1328 '''
1322
1329
1323 ret = numpy.array(self.__times)
1330 ret = numpy.array(self.__times)
1324 ret.sort()
1331 ret.sort()
1325 return ret
1332 return ret
1326
1333
1327 @property
1334 @property
1328 def min_time(self):
1335 def min_time(self):
1329 '''
1336 '''
1330 Return the minimun time value
1337 Return the minimun time value
1331 '''
1338 '''
1332
1339
1333 return self.times[0]
1340 return self.times[0]
1334
1341
1335 @property
1342 @property
1336 def max_time(self):
1343 def max_time(self):
1337 '''
1344 '''
1338 Return the maximun time value
1345 Return the maximun time value
1339 '''
1346 '''
1340
1347
1341 return self.times[-1]
1348 return self.times[-1]
1342
1349
1343 @property
1350 @property
1344 def heights(self):
1351 def heights(self):
1345 '''
1352 '''
1346 Return the list of heights of the current data
1353 Return the list of heights of the current data
1347 '''
1354 '''
1348
1355
1349 return numpy.array(self.__heights[-1])
1356 return numpy.array(self.__heights[-1])
1350
1357
1351 @staticmethod
1358 @staticmethod
1352 def roundFloats(obj):
1359 def roundFloats(obj):
1353 if isinstance(obj, list):
1360 if isinstance(obj, list):
1354 return list(map(PlotterData.roundFloats, obj))
1361 return list(map(PlotterData.roundFloats, obj))
1355 elif isinstance(obj, float):
1362 elif isinstance(obj, float):
1356 return round(obj, 2)
1363 return round(obj, 2)
@@ -1,799 +1,801
1
1
2 import os
2 import os
3 import sys
3 import sys
4 import zmq
4 import zmq
5 import time
5 import time
6 import datetime
6 import datetime
7 from functools import wraps
7 from functools import wraps
8 import numpy
8 import numpy
9 import matplotlib
9 import matplotlib
10
10
11 if 'BACKEND' in os.environ:
11 if 'BACKEND' in os.environ:
12 matplotlib.use(os.environ['BACKEND'])
12 matplotlib.use(os.environ['BACKEND'])
13 elif 'linux' in sys.platform:
13 elif 'linux' in sys.platform:
14 matplotlib.use("TkAgg")
14 matplotlib.use("TkAgg")
15 elif 'darwin' in sys.platform:
15 elif 'darwin' in sys.platform:
16 matplotlib.use('TkAgg')
16 matplotlib.use('TkAgg')
17 else:
17 else:
18 from schainpy.utils import log
18 from schainpy.utils import log
19 log.warning('Using default Backend="Agg"', 'INFO')
19 log.warning('Using default Backend="Agg"', 'INFO')
20 matplotlib.use('Agg')
20 matplotlib.use('Agg')
21
21
22 import matplotlib.pyplot as plt
22 import matplotlib.pyplot as plt
23 from matplotlib.patches import Polygon
23 from matplotlib.patches import Polygon
24 from mpl_toolkits.axes_grid1 import make_axes_locatable
24 from mpl_toolkits.axes_grid1 import make_axes_locatable
25 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
25 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
26
26
27 from schainpy.model.data.jrodata import PlotterData
27 from schainpy.model.data.jrodata import PlotterData
28 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
28 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
29 from schainpy.utils import log
29 from schainpy.utils import log
30
30
31 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
31 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
32 blu_values = matplotlib.pyplot.get_cmap(
32 blu_values = matplotlib.pyplot.get_cmap(
33 'seismic_r', 20)(numpy.arange(20))[10:15]
33 'seismic_r', 20)(numpy.arange(20))[10:15]
34 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
34 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
35 'jro', numpy.vstack((blu_values, jet_values)))
35 'jro', numpy.vstack((blu_values, jet_values)))
36 matplotlib.pyplot.register_cmap(cmap=ncmap)
36 matplotlib.pyplot.register_cmap(cmap=ncmap)
37
37
38 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis',
38 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis',
39 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
39 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
40
40
41 EARTH_RADIUS = 6.3710e3
41 EARTH_RADIUS = 6.3710e3
42
42
43
43
44 def ll2xy(lat1, lon1, lat2, lon2):
44 def ll2xy(lat1, lon1, lat2, lon2):
45
45
46 p = 0.017453292519943295
46 p = 0.017453292519943295
47 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
47 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
48 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
48 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
49 r = 12742 * numpy.arcsin(numpy.sqrt(a))
49 r = 12742 * numpy.arcsin(numpy.sqrt(a))
50 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
50 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
51 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
51 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
52 theta = -theta + numpy.pi/2
52 theta = -theta + numpy.pi/2
53 return r*numpy.cos(theta), r*numpy.sin(theta)
53 return r*numpy.cos(theta), r*numpy.sin(theta)
54
54
55
55
56 def km2deg(km):
56 def km2deg(km):
57 '''
57 '''
58 Convert distance in km to degrees
58 Convert distance in km to degrees
59 '''
59 '''
60
60
61 return numpy.rad2deg(km/EARTH_RADIUS)
61 return numpy.rad2deg(km/EARTH_RADIUS)
62
62
63
63
64 def figpause(interval):
64 def figpause(interval):
65 backend = plt.rcParams['backend']
65 backend = plt.rcParams['backend']
66 if backend in matplotlib.rcsetup.interactive_bk:
66 if backend in matplotlib.rcsetup.interactive_bk:
67 figManager = matplotlib._pylab_helpers.Gcf.get_active()
67 figManager = matplotlib._pylab_helpers.Gcf.get_active()
68 if figManager is not None:
68 if figManager is not None:
69 canvas = figManager.canvas
69 canvas = figManager.canvas
70 if canvas.figure.stale:
70 if canvas.figure.stale:
71 canvas.draw()
71 canvas.draw()
72 try:
72 try:
73 canvas.start_event_loop(interval)
73 canvas.start_event_loop(interval)
74 except:
74 except:
75 pass
75 pass
76 return
76 return
77
77
78
78
79 def popup(message):
79 def popup(message):
80 '''
80 '''
81 '''
81 '''
82
82
83 fig = plt.figure(figsize=(12, 8), facecolor='r')
83 fig = plt.figure(figsize=(12, 8), facecolor='r')
84 text = '\n'.join([s.strip() for s in message.split(':')])
84 text = '\n'.join([s.strip() for s in message.split(':')])
85 fig.text(0.01, 0.5, text, ha='left', va='center',
85 fig.text(0.01, 0.5, text, ha='left', va='center',
86 size='20', weight='heavy', color='w')
86 size='20', weight='heavy', color='w')
87 fig.show()
87 fig.show()
88 figpause(1000)
88 figpause(1000)
89
89
90
90
91 class Throttle(object):
91 class Throttle(object):
92 '''
92 '''
93 Decorator that prevents a function from being called more than once every
93 Decorator that prevents a function from being called more than once every
94 time period.
94 time period.
95 To create a function that cannot be called more than once a minute, but
95 To create a function that cannot be called more than once a minute, but
96 will sleep until it can be called:
96 will sleep until it can be called:
97 @Throttle(minutes=1)
97 @Throttle(minutes=1)
98 def foo():
98 def foo():
99 pass
99 pass
100
100
101 for i in range(10):
101 for i in range(10):
102 foo()
102 foo()
103 print "This function has run %s times." % i
103 print "This function has run %s times." % i
104 '''
104 '''
105
105
106 def __init__(self, seconds=0, minutes=0, hours=0):
106 def __init__(self, seconds=0, minutes=0, hours=0):
107 self.throttle_period = datetime.timedelta(
107 self.throttle_period = datetime.timedelta(
108 seconds=seconds, minutes=minutes, hours=hours
108 seconds=seconds, minutes=minutes, hours=hours
109 )
109 )
110
110
111 self.time_of_last_call = datetime.datetime.min
111 self.time_of_last_call = datetime.datetime.min
112
112
113 def __call__(self, fn):
113 def __call__(self, fn):
114 @wraps(fn)
114 @wraps(fn)
115 def wrapper(*args, **kwargs):
115 def wrapper(*args, **kwargs):
116 coerce = kwargs.pop('coerce', None)
116 coerce = kwargs.pop('coerce', None)
117 if coerce:
117 if coerce:
118 self.time_of_last_call = datetime.datetime.now()
118 self.time_of_last_call = datetime.datetime.now()
119 return fn(*args, **kwargs)
119 return fn(*args, **kwargs)
120 else:
120 else:
121 now = datetime.datetime.now()
121 now = datetime.datetime.now()
122 time_since_last_call = now - self.time_of_last_call
122 time_since_last_call = now - self.time_of_last_call
123 time_left = self.throttle_period - time_since_last_call
123 time_left = self.throttle_period - time_since_last_call
124
124
125 if time_left > datetime.timedelta(seconds=0):
125 if time_left > datetime.timedelta(seconds=0):
126 return
126 return
127
127
128 self.time_of_last_call = datetime.datetime.now()
128 self.time_of_last_call = datetime.datetime.now()
129 return fn(*args, **kwargs)
129 return fn(*args, **kwargs)
130
130
131 return wrapper
131 return wrapper
132
132
133 def apply_throttle(value):
133 def apply_throttle(value):
134
134
135 @Throttle(seconds=value)
135 @Throttle(seconds=value)
136 def fnThrottled(fn):
136 def fnThrottled(fn):
137 fn()
137 fn()
138
138
139 return fnThrottled
139 return fnThrottled
140
140
141 @MPDecorator
141 @MPDecorator
142 class Plotter(ProcessingUnit):
142 class Plotter(ProcessingUnit):
143 '''
143 '''
144 Proccessing unit to handle plot operations
144 Proccessing unit to handle plot operations
145 '''
145 '''
146
146
147 def __init__(self):
147 def __init__(self):
148
148
149 ProcessingUnit.__init__(self)
149 ProcessingUnit.__init__(self)
150
150
151 def setup(self, **kwargs):
151 def setup(self, **kwargs):
152
152
153 self.connections = 0
153 self.connections = 0
154 self.web_address = kwargs.get('web_server', False)
154 self.web_address = kwargs.get('web_server', False)
155 self.realtime = kwargs.get('realtime', False)
155 self.realtime = kwargs.get('realtime', False)
156 self.localtime = kwargs.get('localtime', True)
156 self.localtime = kwargs.get('localtime', True)
157 self.buffering = kwargs.get('buffering', True)
157 self.buffering = kwargs.get('buffering', True)
158 self.throttle = kwargs.get('throttle', 2)
158 self.throttle = kwargs.get('throttle', 2)
159 self.exp_code = kwargs.get('exp_code', None)
159 self.exp_code = kwargs.get('exp_code', None)
160 self.set_ready = apply_throttle(self.throttle)
160 self.set_ready = apply_throttle(self.throttle)
161 self.dates = []
161 self.dates = []
162 self.data = PlotterData(
162 self.data = PlotterData(
163 self.plots, self.throttle, self.exp_code, self.buffering)
163 self.plots, self.throttle, self.exp_code, self.buffering)
164 self.isConfig = True
164 self.isConfig = True
165
165
166 def ready(self):
166 def ready(self):
167 '''
167 '''
168 Set dataOut ready
168 Set dataOut ready
169 '''
169 '''
170
170
171 self.data.ready = True
171 self.data.ready = True
172 self.dataOut.data_plt = self.data
172 self.dataOut.data_plt = self.data
173
173
174 def run(self, realtime=True, localtime=True, buffering=True,
174 def run(self, realtime=True, localtime=True, buffering=True,
175 throttle=2, exp_code=None, web_server=None):
175 throttle=2, exp_code=None, web_server=None):
176
176
177 if not self.isConfig:
177 if not self.isConfig:
178 self.setup(realtime=realtime, localtime=localtime,
178 self.setup(realtime=realtime, localtime=localtime,
179 buffering=buffering, throttle=throttle, exp_code=exp_code,
179 buffering=buffering, throttle=throttle, exp_code=exp_code,
180 web_server=web_server)
180 web_server=web_server)
181
181
182 if self.web_address:
182 if self.web_address:
183 log.success(
183 log.success(
184 'Sending to web: {}'.format(self.web_address),
184 'Sending to web: {}'.format(self.web_address),
185 self.name
185 self.name
186 )
186 )
187 self.context = zmq.Context()
187 self.context = zmq.Context()
188 self.sender_web = self.context.socket(zmq.REQ)
188 self.sender_web = self.context.socket(zmq.REQ)
189 self.sender_web.connect(self.web_address)
189 self.sender_web.connect(self.web_address)
190 self.poll = zmq.Poller()
190 self.poll = zmq.Poller()
191 self.poll.register(self.sender_web, zmq.POLLIN)
191 self.poll.register(self.sender_web, zmq.POLLIN)
192 time.sleep(1)
192 time.sleep(1)
193
193
194 # t = Thread(target=self.event_monitor, args=(monitor,))
194 # t = Thread(target=self.event_monitor, args=(monitor,))
195 # t.start()
195 # t.start()
196
196
197 self.dataOut = self.dataIn
197 self.dataOut = self.dataIn
198 self.data.ready = False
198 self.data.ready = False
199
199
200 if self.dataOut.flagNoData:
200 if self.dataOut.flagNoData:
201 coerce = True
201 coerce = True
202 else:
202 else:
203 coerce = False
203 coerce = False
204
204
205 if self.dataOut.type == 'Parameters':
205 if self.dataOut.type == 'Parameters':
206 tm = self.dataOut.utctimeInit
206 tm = self.dataOut.utctimeInit
207 else:
207 else:
208 tm = self.dataOut.utctime
208 tm = self.dataOut.utctime
209 if self.dataOut.useLocalTime:
209 if self.dataOut.useLocalTime:
210 if not self.localtime:
210 if not self.localtime:
211 tm += time.timezone
211 tm += time.timezone
212 dt = datetime.datetime.fromtimestamp(tm).date()
212 dt = datetime.datetime.fromtimestamp(tm).date()
213 else:
213 else:
214 if self.localtime:
214 if self.localtime:
215 tm -= time.timezone
215 tm -= time.timezone
216 dt = datetime.datetime.utcfromtimestamp(tm).date()
216 dt = datetime.datetime.utcfromtimestamp(tm).date()
217 if dt not in self.dates:
217 if dt not in self.dates:
218 if self.data:
218 if self.data:
219 self.ready()
219 self.ready()
220 self.data.setup()
220 self.data.setup()
221 self.dates.append(dt)
221 self.dates.append(dt)
222
222
223 self.data.update(self.dataOut, tm)
223 self.data.update(self.dataOut, tm)
224
224
225 if False: # TODO check when publishers ends
225 if False: # TODO check when publishers ends
226 self.connections -= 1
226 self.connections -= 1
227 if self.connections == 0 and dt in self.dates:
227 if self.connections == 0 and dt in self.dates:
228 self.data.ended = True
228 self.data.ended = True
229 self.ready()
229 self.ready()
230 time.sleep(1)
230 time.sleep(1)
231 else:
231 else:
232 if self.realtime:
232 if self.realtime:
233 self.ready()
233 self.ready()
234 if self.web_address:
234 if self.web_address:
235 retries = 5
235 retries = 5
236 while True:
236 while True:
237 self.sender_web.send(self.data.jsonify())
237 self.sender_web.send(self.data.jsonify())
238 socks = dict(self.poll.poll(5000))
238 socks = dict(self.poll.poll(5000))
239 if socks.get(self.sender_web) == zmq.POLLIN:
239 if socks.get(self.sender_web) == zmq.POLLIN:
240 reply = self.sender_web.recv_string()
240 reply = self.sender_web.recv_string()
241 if reply == 'ok':
241 if reply == 'ok':
242 log.log("Response from server ok", self.name)
242 log.log("Response from server ok", self.name)
243 break
243 break
244 else:
244 else:
245 log.warning(
245 log.warning(
246 "Malformed reply from server: {}".format(reply), self.name)
246 "Malformed reply from server: {}".format(reply), self.name)
247
247
248 else:
248 else:
249 log.warning(
249 log.warning(
250 "No response from server, retrying...", self.name)
250 "No response from server, retrying...", self.name)
251 self.sender_web.setsockopt(zmq.LINGER, 0)
251 self.sender_web.setsockopt(zmq.LINGER, 0)
252 self.sender_web.close()
252 self.sender_web.close()
253 self.poll.unregister(self.sender_web)
253 self.poll.unregister(self.sender_web)
254 retries -= 1
254 retries -= 1
255 if retries == 0:
255 if retries == 0:
256 log.error(
256 log.error(
257 "Server seems to be offline, abandoning", self.name)
257 "Server seems to be offline, abandoning", self.name)
258 self.sender_web = self.context.socket(zmq.REQ)
258 self.sender_web = self.context.socket(zmq.REQ)
259 self.sender_web.connect(self.web_address)
259 self.sender_web.connect(self.web_address)
260 self.poll.register(self.sender_web, zmq.POLLIN)
260 self.poll.register(self.sender_web, zmq.POLLIN)
261 time.sleep(1)
261 time.sleep(1)
262 break
262 break
263 self.sender_web = self.context.socket(zmq.REQ)
263 self.sender_web = self.context.socket(zmq.REQ)
264 self.sender_web.connect(self.web_address)
264 self.sender_web.connect(self.web_address)
265 self.poll.register(self.sender_web, zmq.POLLIN)
265 self.poll.register(self.sender_web, zmq.POLLIN)
266 time.sleep(1)
266 time.sleep(1)
267 else:
267 else:
268 self.set_ready(self.ready, coerce=coerce)
268 self.set_ready(self.ready, coerce=coerce)
269
269
270 return
270 return
271
271
272 def close(self):
272 def close(self):
273 pass
273 pass
274
274
275
275
276 @MPDecorator
276 @MPDecorator
277 class Plot(Operation):
277 class Plot(Operation):
278 '''
278 '''
279 Base class for Schain plotting operations
279 Base class for Schain plotting operations
280 '''
280 '''
281
281
282 CODE = 'Figure'
282 CODE = 'Figure'
283 colormap = 'jro'
283 colormap = 'jro'
284 bgcolor = 'white'
284 bgcolor = 'white'
285 __missing = 1E30
285 __missing = 1E30
286
286
287 __attrs__ = ['show', 'save', 'xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax',
287 __attrs__ = ['show', 'save', 'xmin', 'xmax', 'ymin', 'ymax', 'zmin', 'zmax',
288 'zlimits', 'xlabel', 'ylabel', 'xaxis', 'cb_label', 'title',
288 'zlimits', 'xlabel', 'ylabel', 'xaxis', 'cb_label', 'title',
289 'colorbar', 'bgcolor', 'width', 'height', 'localtime', 'oneFigure',
289 'colorbar', 'bgcolor', 'width', 'height', 'localtime', 'oneFigure',
290 'showprofile', 'decimation', 'pause']
290 'showprofile', 'decimation', 'pause']
291
291
292 def __init__(self):
292 def __init__(self):
293
293
294 Operation.__init__(self)
294 Operation.__init__(self)
295 self.isConfig = False
295 self.isConfig = False
296 self.isPlotConfig = False
296 self.isPlotConfig = False
297
297
298 def __fmtTime(self, x, pos):
298 def __fmtTime(self, x, pos):
299 '''
299 '''
300 '''
300 '''
301
301
302 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
302 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
303
303
304 def __setup(self, **kwargs):
304 def __setup(self, **kwargs):
305 '''
305 '''
306 Initialize variables
306 Initialize variables
307 '''
307 '''
308
308
309 self.figures = []
309 self.figures = []
310 self.axes = []
310 self.axes = []
311 self.cb_axes = []
311 self.cb_axes = []
312 self.localtime = kwargs.pop('localtime', True)
312 self.localtime = kwargs.pop('localtime', True)
313 self.show = kwargs.get('show', True)
313 self.show = kwargs.get('show', True)
314 self.save = kwargs.get('save', False)
314 self.save = kwargs.get('save', False)
315 self.ftp = kwargs.get('ftp', False)
315 self.ftp = kwargs.get('ftp', False)
316 self.colormap = kwargs.get('colormap', self.colormap)
316 self.colormap = kwargs.get('colormap', self.colormap)
317 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
317 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
318 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
318 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
319 self.colormaps = kwargs.get('colormaps', None)
319 self.colormaps = kwargs.get('colormaps', None)
320 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
320 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
321 self.showprofile = kwargs.get('showprofile', False)
321 self.showprofile = kwargs.get('showprofile', False)
322 self.title = kwargs.get('wintitle', self.CODE.upper())
322 self.title = kwargs.get('wintitle', self.CODE.upper())
323 self.cb_label = kwargs.get('cb_label', None)
323 self.cb_label = kwargs.get('cb_label', None)
324 self.cb_labels = kwargs.get('cb_labels', None)
324 self.cb_labels = kwargs.get('cb_labels', None)
325 self.labels = kwargs.get('labels', None)
325 self.labels = kwargs.get('labels', None)
326 self.xaxis = kwargs.get('xaxis', 'frequency')
326 self.xaxis = kwargs.get('xaxis', 'frequency')
327 self.zmin = kwargs.get('zmin', None)
327 self.zmin = kwargs.get('zmin', None)
328 self.zmax = kwargs.get('zmax', None)
328 self.zmax = kwargs.get('zmax', None)
329 self.zlimits = kwargs.get('zlimits', None)
329 self.zlimits = kwargs.get('zlimits', None)
330 self.xmin = kwargs.get('xmin', None)
330 self.xmin = kwargs.get('xmin', None)
331 self.xmax = kwargs.get('xmax', None)
331 self.xmax = kwargs.get('xmax', None)
332 self.xrange = kwargs.get('xrange', 12)
332 self.xrange = kwargs.get('xrange', 12)
333 self.xscale = kwargs.get('xscale', None)
333 self.xscale = kwargs.get('xscale', None)
334 self.ymin = kwargs.get('ymin', None)
334 self.ymin = kwargs.get('ymin', None)
335 self.ymax = kwargs.get('ymax', None)
335 self.ymax = kwargs.get('ymax', None)
336 self.yscale = kwargs.get('yscale', None)
336 self.yscale = kwargs.get('yscale', None)
337 self.xlabel = kwargs.get('xlabel', None)
337 self.xlabel = kwargs.get('xlabel', None)
338 self.decimation = kwargs.get('decimation', None)
338 self.decimation = kwargs.get('decimation', None)
339 self.showSNR = kwargs.get('showSNR', False)
339 self.showSNR = kwargs.get('showSNR', False)
340 self.oneFigure = kwargs.get('oneFigure', True)
340 self.oneFigure = kwargs.get('oneFigure', True)
341 self.width = kwargs.get('width', None)
341 self.width = kwargs.get('width', None)
342 self.height = kwargs.get('height', None)
342 self.height = kwargs.get('height', None)
343 self.colorbar = kwargs.get('colorbar', True)
343 self.colorbar = kwargs.get('colorbar', True)
344 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
344 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
345 self.channels = kwargs.get('channels', None)
345 self.channels = kwargs.get('channels', None)
346 self.titles = kwargs.get('titles', [])
346 self.titles = kwargs.get('titles', [])
347 self.polar = False
347 self.polar = False
348 self.type = kwargs.get('type', 'iq')
348 self.grid = kwargs.get('grid', False)
349 self.grid = kwargs.get('grid', False)
349 self.pause = kwargs.get('pause', False)
350 self.pause = kwargs.get('pause', False)
350 self.save_labels = kwargs.get('save_labels', None)
351 self.save_labels = kwargs.get('save_labels', None)
351 self.realtime = kwargs.get('realtime', True)
352 self.realtime = kwargs.get('realtime', True)
352 self.buffering = kwargs.get('buffering', True)
353 self.buffering = kwargs.get('buffering', True)
353 self.throttle = kwargs.get('throttle', 2)
354 self.throttle = kwargs.get('throttle', 2)
354 self.exp_code = kwargs.get('exp_code', None)
355 self.exp_code = kwargs.get('exp_code', None)
355 self.__throttle_plot = apply_throttle(self.throttle)
356 self.__throttle_plot = apply_throttle(self.throttle)
356 self.data = PlotterData(
357 self.data = PlotterData(
357 self.CODE, self.throttle, self.exp_code, self.buffering)
358 self.CODE, self.throttle, self.exp_code, self.buffering)
358
359
359 def __setup_plot(self):
360 def __setup_plot(self):
360 '''
361 '''
361 Common setup for all figures, here figures and axes are created
362 Common setup for all figures, here figures and axes are created
362 '''
363 '''
363
364
364 self.setup()
365 self.setup()
365
366
366 self.time_label = 'LT' if self.localtime else 'UTC'
367 self.time_label = 'LT' if self.localtime else 'UTC'
367 if self.data.localtime:
368 if self.data.localtime:
368 self.getDateTime = datetime.datetime.fromtimestamp
369 self.getDateTime = datetime.datetime.fromtimestamp
369 else:
370 else:
370 self.getDateTime = datetime.datetime.utcfromtimestamp
371 self.getDateTime = datetime.datetime.utcfromtimestamp
371
372
372 if self.width is None:
373 if self.width is None:
373 self.width = 8
374 self.width = 8
374
375
375 self.figures = []
376 self.figures = []
376 self.axes = []
377 self.axes = []
377 self.cb_axes = []
378 self.cb_axes = []
378 self.pf_axes = []
379 self.pf_axes = []
379 self.cmaps = []
380 self.cmaps = []
380
381
381 size = '15%' if self.ncols == 1 else '30%'
382 size = '15%' if self.ncols == 1 else '30%'
382 pad = '4%' if self.ncols == 1 else '8%'
383 pad = '4%' if self.ncols == 1 else '8%'
383
384
384 if self.oneFigure:
385 if self.oneFigure:
385 if self.height is None:
386 if self.height is None:
386 self.height = 1.4 * self.nrows + 1
387 self.height = 1.4 * self.nrows + 1
387 fig = plt.figure(figsize=(self.width, self.height),
388 fig = plt.figure(figsize=(self.width, self.height),
388 edgecolor='k',
389 edgecolor='k',
389 facecolor='w')
390 facecolor='w')
390 self.figures.append(fig)
391 self.figures.append(fig)
391 for n in range(self.nplots):
392 for n in range(self.nplots):
392 ax = fig.add_subplot(self.nrows, self.ncols,
393 ax = fig.add_subplot(self.nrows, self.ncols,
393 n + 1, polar=self.polar)
394 n + 1, polar=self.polar)
394 ax.tick_params(labelsize=8)
395 ax.tick_params(labelsize=8)
395 ax.firsttime = True
396 ax.firsttime = True
396 ax.index = 0
397 ax.index = 0
397 ax.press = None
398 ax.press = None
398 self.axes.append(ax)
399 self.axes.append(ax)
399 if self.showprofile:
400 if self.showprofile:
400 cax = self.__add_axes(ax, size=size, pad=pad)
401 cax = self.__add_axes(ax, size=size, pad=pad)
401 cax.tick_params(labelsize=8)
402 cax.tick_params(labelsize=8)
402 self.pf_axes.append(cax)
403 self.pf_axes.append(cax)
403 else:
404 else:
404 if self.height is None:
405 if self.height is None:
405 self.height = 3
406 self.height = 3
406 for n in range(self.nplots):
407 for n in range(self.nplots):
407 fig = plt.figure(figsize=(self.width, self.height),
408 fig = plt.figure(figsize=(self.width, self.height),
408 edgecolor='k',
409 edgecolor='k',
409 facecolor='w')
410 facecolor='w')
410 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
411 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
411 ax.tick_params(labelsize=8)
412 ax.tick_params(labelsize=8)
412 ax.firsttime = True
413 ax.firsttime = True
413 ax.index = 0
414 ax.index = 0
414 ax.press = None
415 ax.press = None
415 self.figures.append(fig)
416 self.figures.append(fig)
416 self.axes.append(ax)
417 self.axes.append(ax)
417 if self.showprofile:
418 if self.showprofile:
418 cax = self.__add_axes(ax, size=size, pad=pad)
419 cax = self.__add_axes(ax, size=size, pad=pad)
419 cax.tick_params(labelsize=8)
420 cax.tick_params(labelsize=8)
420 self.pf_axes.append(cax)
421 self.pf_axes.append(cax)
421
422
422 for n in range(self.nrows):
423 for n in range(self.nrows):
423 if self.colormaps is not None:
424 if self.colormaps is not None:
424 cmap = plt.get_cmap(self.colormaps[n])
425 cmap = plt.get_cmap(self.colormaps[n])
425 else:
426 else:
426 cmap = plt.get_cmap(self.colormap)
427 cmap = plt.get_cmap(self.colormap)
427 cmap.set_bad(self.bgcolor, 1.)
428 cmap.set_bad(self.bgcolor, 1.)
428 self.cmaps.append(cmap)
429 self.cmaps.append(cmap)
429
430
430 for fig in self.figures:
431 for fig in self.figures:
431 fig.canvas.mpl_connect('key_press_event', self.OnKeyPress)
432 fig.canvas.mpl_connect('key_press_event', self.OnKeyPress)
432 fig.canvas.mpl_connect('scroll_event', self.OnBtnScroll)
433 fig.canvas.mpl_connect('scroll_event', self.OnBtnScroll)
433 fig.canvas.mpl_connect('button_press_event', self.onBtnPress)
434 fig.canvas.mpl_connect('button_press_event', self.onBtnPress)
434 fig.canvas.mpl_connect('motion_notify_event', self.onMotion)
435 fig.canvas.mpl_connect('motion_notify_event', self.onMotion)
435 fig.canvas.mpl_connect('button_release_event', self.onBtnRelease)
436 fig.canvas.mpl_connect('button_release_event', self.onBtnRelease)
436 if self.show:
437 if self.show:
437 fig.show()
438 fig.show()
438
439
439 def OnKeyPress(self, event):
440 def OnKeyPress(self, event):
440 '''
441 '''
441 Event for pressing keys (up, down) change colormap
442 Event for pressing keys (up, down) change colormap
442 '''
443 '''
443 ax = event.inaxes
444 ax = event.inaxes
444 if ax in self.axes:
445 if ax in self.axes:
445 if event.key == 'down':
446 if event.key == 'down':
446 ax.index += 1
447 ax.index += 1
447 elif event.key == 'up':
448 elif event.key == 'up':
448 ax.index -= 1
449 ax.index -= 1
449 if ax.index < 0:
450 if ax.index < 0:
450 ax.index = len(CMAPS) - 1
451 ax.index = len(CMAPS) - 1
451 elif ax.index == len(CMAPS):
452 elif ax.index == len(CMAPS):
452 ax.index = 0
453 ax.index = 0
453 cmap = CMAPS[ax.index]
454 cmap = CMAPS[ax.index]
454 ax.cbar.set_cmap(cmap)
455 ax.cbar.set_cmap(cmap)
455 ax.cbar.draw_all()
456 ax.cbar.draw_all()
456 ax.plt.set_cmap(cmap)
457 ax.plt.set_cmap(cmap)
457 ax.cbar.patch.figure.canvas.draw()
458 ax.cbar.patch.figure.canvas.draw()
458 self.colormap = cmap.name
459 self.colormap = cmap.name
459
460
460 def OnBtnScroll(self, event):
461 def OnBtnScroll(self, event):
461 '''
462 '''
462 Event for scrolling, scale figure
463 Event for scrolling, scale figure
463 '''
464 '''
464 cb_ax = event.inaxes
465 cb_ax = event.inaxes
465 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
466 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
466 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
467 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
467 pt = ax.cbar.ax.bbox.get_points()[:, 1]
468 pt = ax.cbar.ax.bbox.get_points()[:, 1]
468 nrm = ax.cbar.norm
469 nrm = ax.cbar.norm
469 vmin, vmax, p0, p1, pS = (
470 vmin, vmax, p0, p1, pS = (
470 nrm.vmin, nrm.vmax, pt[0], pt[1], event.y)
471 nrm.vmin, nrm.vmax, pt[0], pt[1], event.y)
471 scale = 2 if event.step == 1 else 0.5
472 scale = 2 if event.step == 1 else 0.5
472 point = vmin + (vmax - vmin) / (p1 - p0) * (pS - p0)
473 point = vmin + (vmax - vmin) / (p1 - p0) * (pS - p0)
473 ax.cbar.norm.vmin = point - scale * (point - vmin)
474 ax.cbar.norm.vmin = point - scale * (point - vmin)
474 ax.cbar.norm.vmax = point - scale * (point - vmax)
475 ax.cbar.norm.vmax = point - scale * (point - vmax)
475 ax.plt.set_norm(ax.cbar.norm)
476 ax.plt.set_norm(ax.cbar.norm)
476 ax.cbar.draw_all()
477 ax.cbar.draw_all()
477 ax.cbar.patch.figure.canvas.draw()
478 ax.cbar.patch.figure.canvas.draw()
478
479
479 def onBtnPress(self, event):
480 def onBtnPress(self, event):
480 '''
481 '''
481 Event for mouse button press
482 Event for mouse button press
482 '''
483 '''
483 cb_ax = event.inaxes
484 cb_ax = event.inaxes
484 if cb_ax is None:
485 if cb_ax is None:
485 return
486 return
486
487
487 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
488 if cb_ax in [ax.cbar.ax for ax in self.axes if ax.cbar]:
488 cb_ax.press = event.x, event.y
489 cb_ax.press = event.x, event.y
489 else:
490 else:
490 cb_ax.press = None
491 cb_ax.press = None
491
492
492 def onMotion(self, event):
493 def onMotion(self, event):
493 '''
494 '''
494 Event for move inside colorbar
495 Event for move inside colorbar
495 '''
496 '''
496 cb_ax = event.inaxes
497 cb_ax = event.inaxes
497 if cb_ax is None:
498 if cb_ax is None:
498 return
499 return
499 if cb_ax not in [ax.cbar.ax for ax in self.axes if ax.cbar]:
500 if cb_ax not in [ax.cbar.ax for ax in self.axes if ax.cbar]:
500 return
501 return
501 if cb_ax.press is None:
502 if cb_ax.press is None:
502 return
503 return
503
504
504 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
505 ax = [ax for ax in self.axes if cb_ax == ax.cbar.ax][0]
505 xprev, yprev = cb_ax.press
506 xprev, yprev = cb_ax.press
506 dx = event.x - xprev
507 dx = event.x - xprev
507 dy = event.y - yprev
508 dy = event.y - yprev
508 cb_ax.press = event.x, event.y
509 cb_ax.press = event.x, event.y
509 scale = ax.cbar.norm.vmax - ax.cbar.norm.vmin
510 scale = ax.cbar.norm.vmax - ax.cbar.norm.vmin
510 perc = 0.03
511 perc = 0.03
511
512
512 if event.button == 1:
513 if event.button == 1:
513 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
514 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
514 ax.cbar.norm.vmax -= (perc * scale) * numpy.sign(dy)
515 ax.cbar.norm.vmax -= (perc * scale) * numpy.sign(dy)
515 elif event.button == 3:
516 elif event.button == 3:
516 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
517 ax.cbar.norm.vmin -= (perc * scale) * numpy.sign(dy)
517 ax.cbar.norm.vmax += (perc * scale) * numpy.sign(dy)
518 ax.cbar.norm.vmax += (perc * scale) * numpy.sign(dy)
518
519
519 ax.cbar.draw_all()
520 ax.cbar.draw_all()
520 ax.plt.set_norm(ax.cbar.norm)
521 ax.plt.set_norm(ax.cbar.norm)
521 ax.cbar.patch.figure.canvas.draw()
522 ax.cbar.patch.figure.canvas.draw()
522
523
523 def onBtnRelease(self, event):
524 def onBtnRelease(self, event):
524 '''
525 '''
525 Event for mouse button release
526 Event for mouse button release
526 '''
527 '''
527 cb_ax = event.inaxes
528 cb_ax = event.inaxes
528 if cb_ax is not None:
529 if cb_ax is not None:
529 cb_ax.press = None
530 cb_ax.press = None
530
531
531 def __add_axes(self, ax, size='30%', pad='8%'):
532 def __add_axes(self, ax, size='30%', pad='8%'):
532 '''
533 '''
533 Add new axes to the given figure
534 Add new axes to the given figure
534 '''
535 '''
535 divider = make_axes_locatable(ax)
536 divider = make_axes_locatable(ax)
536 nax = divider.new_horizontal(size=size, pad=pad)
537 nax = divider.new_horizontal(size=size, pad=pad)
537 ax.figure.add_axes(nax)
538 ax.figure.add_axes(nax)
538 return nax
539 return nax
539
540
540 def setup(self):
541 def setup(self):
541 '''
542 '''
542 This method should be implemented in the child class, the following
543 This method should be implemented in the child class, the following
543 attributes should be set:
544 attributes should be set:
544
545
545 self.nrows: number of rows
546 self.nrows: number of rows
546 self.ncols: number of cols
547 self.ncols: number of cols
547 self.nplots: number of plots (channels or pairs)
548 self.nplots: number of plots (channels or pairs)
548 self.ylabel: label for Y axes
549 self.ylabel: label for Y axes
549 self.titles: list of axes title
550 self.titles: list of axes title
550
551
551 '''
552 '''
552 raise NotImplementedError
553 raise NotImplementedError
553
554
554 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
555 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
555 '''
556 '''
556 Create a masked array for missing data
557 Create a masked array for missing data
557 '''
558 '''
558 if x_buffer.shape[0] < 2:
559 if x_buffer.shape[0] < 2:
559 return x_buffer, y_buffer, z_buffer
560 return x_buffer, y_buffer, z_buffer
560
561
561 deltas = x_buffer[1:] - x_buffer[0:-1]
562 deltas = x_buffer[1:] - x_buffer[0:-1]
562 x_median = numpy.median(deltas)
563 x_median = numpy.median(deltas)
563
564
564 index = numpy.where(deltas > 5 * x_median)
565 index = numpy.where(deltas > 5 * x_median)
565
566
566 if len(index[0]) != 0:
567 if len(index[0]) != 0:
567 z_buffer[::, index[0], ::] = self.__missing
568 z_buffer[::, index[0], ::] = self.__missing
568 z_buffer = numpy.ma.masked_inside(z_buffer,
569 z_buffer = numpy.ma.masked_inside(z_buffer,
569 0.99 * self.__missing,
570 0.99 * self.__missing,
570 1.01 * self.__missing)
571 1.01 * self.__missing)
571
572
572 return x_buffer, y_buffer, z_buffer
573 return x_buffer, y_buffer, z_buffer
573
574
574 def decimate(self):
575 def decimate(self):
575
576
576 # dx = int(len(self.x)/self.__MAXNUMX) + 1
577 # dx = int(len(self.x)/self.__MAXNUMX) + 1
577 dy = int(len(self.y) / self.decimation) + 1
578 dy = int(len(self.y) / self.decimation) + 1
578
579
579 # x = self.x[::dx]
580 # x = self.x[::dx]
580 x = self.x
581 x = self.x
581 y = self.y[::dy]
582 y = self.y[::dy]
582 z = self.z[::, ::, ::dy]
583 z = self.z[::, ::, ::dy]
583
584
584 return x, y, z
585 return x, y, z
585
586
586 def format(self):
587 def format(self):
587 '''
588 '''
588 Set min and max values, labels, ticks and titles
589 Set min and max values, labels, ticks and titles
589 '''
590 '''
590
591
591 if self.xmin is None:
592 if self.xmin is None:
592 xmin = self.data.min_time
593 xmin = self.data.min_time
593 else:
594 else:
594 if self.xaxis is 'time':
595 if self.xaxis is 'time':
595 dt = self.getDateTime(self.data.min_time)
596 dt = self.getDateTime(self.data.min_time)
596 xmin = (dt.replace(hour=int(self.xmin), minute=0, second=0) -
597 xmin = (dt.replace(hour=int(self.xmin), minute=0, second=0) -
597 datetime.datetime(1970, 1, 1)).total_seconds()
598 datetime.datetime(1970, 1, 1)).total_seconds()
598 if self.data.localtime:
599 if self.data.localtime:
599 xmin += time.timezone
600 xmin += time.timezone
600 else:
601 else:
601 xmin = self.xmin
602 xmin = self.xmin
602
603
603 if self.xmax is None:
604 if self.xmax is None:
604 xmax = xmin + self.xrange * 60 * 60
605 xmax = xmin + self.xrange * 60 * 60
605 else:
606 else:
606 if self.xaxis is 'time':
607 if self.xaxis is 'time':
607 dt = self.getDateTime(self.data.max_time)
608 dt = self.getDateTime(self.data.max_time)
608 xmax = (dt.replace(hour=int(self.xmax), minute=59, second=59) -
609 xmax = (dt.replace(hour=int(self.xmax), minute=59, second=59) -
609 datetime.datetime(1970, 1, 1) + datetime.timedelta(seconds=1)).total_seconds()
610 datetime.datetime(1970, 1, 1) + datetime.timedelta(seconds=1)).total_seconds()
610 if self.data.localtime:
611 if self.data.localtime:
611 xmax += time.timezone
612 xmax += time.timezone
612 else:
613 else:
613 xmax = self.xmax
614 xmax = self.xmax
614
615
615 ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
616 ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
616 ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
617 ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
617 Y = numpy.array([1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000])
618 Y = numpy.array([1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000, 20000, 50000])
618 i = 1 if numpy.where(
619 #i = 1 if numpy.where(
619 abs(ymax-ymin) <= Y)[0][0] < 0 else numpy.where(abs(ymax-ymin) <= Y)[0][0]
620 # abs(ymax-ymin) <= Y)[0][0] < 0 else numpy.where(abs(ymax-ymin) <= Y)[0][0]
620 ystep = Y[i] / 10.
621 #ystep = Y[i] / 10.
621
622 ystep = round(ymax,-1)//5
622 if self.xaxis is not 'time':
623 if self.xaxis is not 'time':
623 X = numpy.array([0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100,
624 X = numpy.array([0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100,
624 200, 500, 1000, 2000, 5000])/2.
625 200, 500, 1000, 2000, 5000, 10000, 20000, 50000])/2.
626
625 i = 1 if numpy.where(
627 i = 1 if numpy.where(
626 abs(xmax-xmin) <= X)[0][0] < 0 else numpy.where(abs(xmax-xmin) <= X)[0][0]
628 abs(xmax-xmin) <= X)[0][0] < 0 else numpy.where(abs(xmax-xmin) <= X)[0][0]
627 xstep = X[i] / 5.
629 xstep = X[i] / 5.
628
630
629 for n, ax in enumerate(self.axes):
631 for n, ax in enumerate(self.axes):
630 if ax.firsttime:
632 if ax.firsttime:
631 ax.set_facecolor(self.bgcolor)
633 ax.set_facecolor(self.bgcolor)
632 ax.yaxis.set_major_locator(MultipleLocator(ystep))
634 ax.yaxis.set_major_locator(MultipleLocator(ystep))
633 if self.xscale:
635 if self.xscale:
634 ax.xaxis.set_major_formatter(FuncFormatter(
636 ax.xaxis.set_major_formatter(FuncFormatter(
635 lambda x, pos: '{0:g}'.format(x*self.xscale)))
637 lambda x, pos: '{0:g}'.format(x*self.xscale)))
636 if self.xscale:
638 if self.xscale:
637 ax.yaxis.set_major_formatter(FuncFormatter(
639 ax.yaxis.set_major_formatter(FuncFormatter(
638 lambda x, pos: '{0:g}'.format(x*self.yscale)))
640 lambda x, pos: '{0:g}'.format(x*self.yscale)))
639 if self.xaxis is 'time':
641 if self.xaxis is 'time':
640 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
642 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
641 ax.xaxis.set_major_locator(LinearLocator(9))
643 ax.xaxis.set_major_locator(LinearLocator(9))
642 else:
644 else:
643 ax.xaxis.set_major_locator(MultipleLocator(xstep))
645 ax.xaxis.set_major_locator(MultipleLocator(xstep))
644 if self.xlabel is not None:
646 if self.xlabel is not None:
645 ax.set_xlabel(self.xlabel)
647 ax.set_xlabel(self.xlabel)
646 ax.set_ylabel(self.ylabel)
648 ax.set_ylabel(self.ylabel)
647 ax.firsttime = False
649 ax.firsttime = False
648 if self.showprofile:
650 if self.showprofile:
649 self.pf_axes[n].set_ylim(ymin, ymax)
651 self.pf_axes[n].set_ylim(ymin, ymax)
650 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
652 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
651 self.pf_axes[n].set_xlabel('dB')
653 self.pf_axes[n].set_xlabel('dB')
652 self.pf_axes[n].grid(b=True, axis='x')
654 self.pf_axes[n].grid(b=True, axis='x')
653 [tick.set_visible(False)
655 [tick.set_visible(False)
654 for tick in self.pf_axes[n].get_yticklabels()]
656 for tick in self.pf_axes[n].get_yticklabels()]
655 if self.colorbar:
657 if self.colorbar:
656 ax.cbar = plt.colorbar(
658 ax.cbar = plt.colorbar(
657 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
659 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
658 ax.cbar.ax.tick_params(labelsize=8)
660 ax.cbar.ax.tick_params(labelsize=8)
659 ax.cbar.ax.press = None
661 ax.cbar.ax.press = None
660 if self.cb_label:
662 if self.cb_label:
661 ax.cbar.set_label(self.cb_label, size=8)
663 ax.cbar.set_label(self.cb_label, size=8)
662 elif self.cb_labels:
664 elif self.cb_labels:
663 ax.cbar.set_label(self.cb_labels[n], size=8)
665 ax.cbar.set_label(self.cb_labels[n], size=8)
664 else:
666 else:
665 ax.cbar = None
667 ax.cbar = None
666 if self.grid:
668 if self.grid:
667 ax.grid(True)
669 ax.grid(True)
668
670
669 if not self.polar:
671 if not self.polar:
670 ax.set_xlim(xmin, xmax)
672 ax.set_xlim(xmin, xmax)
671 ax.set_ylim(ymin, ymax)
673 ax.set_ylim(ymin, ymax)
672 ax.set_title('{} {} {}'.format(
674 ax.set_title('{} {} {}'.format(
673 self.titles[n],
675 self.titles[n],
674 self.getDateTime(self.data.max_time).strftime(
676 self.getDateTime(self.data.max_time).strftime(
675 '%H:%M:%S'),
677 '%H:%M:%S'),
676 self.time_label),
678 self.time_label),
677 size=8)
679 size=8)
678 else:
680 else:
679 ax.set_title('{}'.format(self.titles[n]), size=8)
681 ax.set_title('{}'.format(self.titles[n]), size=8)
680 ax.set_ylim(0, 90)
682 ax.set_ylim(0, 90)
681 ax.set_yticks(numpy.arange(0, 90, 20))
683 ax.set_yticks(numpy.arange(0, 90, 20))
682 ax.yaxis.labelpad = 40
684 ax.yaxis.labelpad = 40
683
685
684 def clear_figures(self):
686 def clear_figures(self):
685 '''
687 '''
686 Reset axes for redraw plots
688 Reset axes for redraw plots
687 '''
689 '''
688
690
689 for ax in self.axes:
691 for ax in self.axes:
690 ax.clear()
692 ax.clear()
691 ax.firsttime = True
693 ax.firsttime = True
692 if ax.cbar:
694 if ax.cbar:
693 ax.cbar.remove()
695 ax.cbar.remove()
694
696
695 def __plot(self):
697 def __plot(self):
696 '''
698 '''
697 Main function to plot, format and save figures
699 Main function to plot, format and save figures
698 '''
700 '''
699
701
700 #try:
702 #try:
701 self.plot()
703 self.plot()
702 self.format()
704 self.format()
703 #except Exception as e:
705 #except Exception as e:
704 # log.warning('{} Plot could not be updated... check data'.format(
706 # log.warning('{} Plot could not be updated... check data'.format(
705 # self.CODE), self.name)
707 # self.CODE), self.name)
706 # log.error(str(e), '')
708 # log.error(str(e), '')
707 # return
709 # return
708
710
709 for n, fig in enumerate(self.figures):
711 for n, fig in enumerate(self.figures):
710 if self.nrows == 0 or self.nplots == 0:
712 if self.nrows == 0 or self.nplots == 0:
711 log.warning('No data', self.name)
713 log.warning('No data', self.name)
712 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
714 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
713 fig.canvas.manager.set_window_title(self.CODE)
715 fig.canvas.manager.set_window_title(self.CODE)
714 continue
716 continue
715
717
716 fig.tight_layout()
718 fig.tight_layout()
717 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
719 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
718 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
720 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
719 fig.canvas.draw()
721 fig.canvas.draw()
720
722
721 if self.save:
723 if self.save:
722
724
723 if self.save_labels:
725 if self.save_labels:
724 labels = self.save_labels
726 labels = self.save_labels
725 else:
727 else:
726 labels = list(range(self.nrows))
728 labels = list(range(self.nrows))
727
729
728 if self.oneFigure:
730 if self.oneFigure:
729 label = ''
731 label = ''
730 else:
732 else:
731 label = '-{}'.format(labels[n])
733 label = '-{}'.format(labels[n])
732 figname = os.path.join(
734 figname = os.path.join(
733 self.save,
735 self.save,
734 self.CODE,
736 self.CODE,
735 '{}{}_{}.png'.format(
737 '{}{}_{}.png'.format(
736 self.CODE,
738 self.CODE,
737 label,
739 label,
738 self.getDateTime(self.data.max_time).strftime(
740 self.getDateTime(self.data.max_time).strftime(
739 '%Y%m%d_%H%M%S'),
741 '%Y%m%d_%H%M%S'),
740 )
742 )
741 )
743 )
742 log.log('Saving figure: {}'.format(figname), self.name)
744 log.log('Saving figure: {}'.format(figname), self.name)
743 if not os.path.isdir(os.path.dirname(figname)):
745 if not os.path.isdir(os.path.dirname(figname)):
744 os.makedirs(os.path.dirname(figname))
746 os.makedirs(os.path.dirname(figname))
745 fig.savefig(figname)
747 fig.savefig(figname)
746
748
747 def plot(self):
749 def plot(self):
748 '''
750 '''
749 Must be defined in the child class
751 Must be defined in the child class
750 '''
752 '''
751 raise NotImplementedError
753 raise NotImplementedError
752
754
753 def run(self, dataOut, **kwargs):
755 def run(self, dataOut, **kwargs):
754
756
755 if dataOut.error:
757 if dataOut.error:
756 coerce = True
758 coerce = True
757 else:
759 else:
758 coerce = False
760 coerce = False
759
761
760 if self.isConfig is False:
762 if self.isConfig is False:
761 self.__setup(**kwargs)
763 self.__setup(**kwargs)
762 self.data.setup()
764 self.data.setup()
763 self.isConfig = True
765 self.isConfig = True
764
766
765 if dataOut.type == 'Parameters':
767 if dataOut.type == 'Parameters':
766 tm = dataOut.utctimeInit
768 tm = dataOut.utctimeInit
767 else:
769 else:
768 tm = dataOut.utctime
770 tm = dataOut.utctime
769
771
770 if dataOut.useLocalTime:
772 if dataOut.useLocalTime:
771 if not self.localtime:
773 if not self.localtime:
772 tm += time.timezone
774 tm += time.timezone
773 else:
775 else:
774 if self.localtime:
776 if self.localtime:
775 tm -= time.timezone
777 tm -= time.timezone
776
778
777 if self.data and (tm - self.data.min_time) >= self.xrange*60*60:
779 if self.data and (tm - self.data.min_time) >= self.xrange*60*60:
778 self.__plot()
780 self.__plot()
779 self.data.setup()
781 self.data.setup()
780 self.clear_figures()
782 self.clear_figures()
781
783
782 self.data.update(dataOut, tm)
784 self.data.update(dataOut, tm)
783
785
784 if self.isPlotConfig is False:
786 if self.isPlotConfig is False:
785 self.__setup_plot()
787 self.__setup_plot()
786 self.isPlotConfig = True
788 self.isPlotConfig = True
787
789
788 if self.realtime:
790 if self.realtime:
789 self.__plot()
791 self.__plot()
790 else:
792 else:
791 self.__throttle_plot(self.__plot, coerce=coerce)
793 self.__throttle_plot(self.__plot, coerce=coerce)
792
794
793 figpause(0.001)
795 figpause(0.001)
794
796
795 def close(self):
797 def close(self):
796
798
797 if self.data and self.pause:
799 if self.data and self.pause:
798 figpause(10)
800 figpause(10)
799
801
@@ -1,616 +1,748
1 '''
1 '''
2 New Plots Operations
2 New Plots Operations
3
3
4 @author: juan.espinoza@jro.igp.gob.pe
4 @author: juan.espinoza@jro.igp.gob.pe
5 '''
5 '''
6
6
7
7
8 import time
8 import time
9 import datetime
9 import datetime
10 import numpy
10 import numpy
11
11
12 from schainpy.model.graphics.jroplot_base import Plot, plt
12 from schainpy.model.graphics.jroplot_base import Plot, plt
13 from schainpy.utils import log
13 from schainpy.utils import log
14
14
15 EARTH_RADIUS = 6.3710e3
15 EARTH_RADIUS = 6.3710e3
16
16
17
17
18 def ll2xy(lat1, lon1, lat2, lon2):
18 def ll2xy(lat1, lon1, lat2, lon2):
19
19
20 p = 0.017453292519943295
20 p = 0.017453292519943295
21 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
21 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
22 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
22 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
23 r = 12742 * numpy.arcsin(numpy.sqrt(a))
23 r = 12742 * numpy.arcsin(numpy.sqrt(a))
24 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
24 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
25 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
25 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
26 theta = -theta + numpy.pi/2
26 theta = -theta + numpy.pi/2
27 return r*numpy.cos(theta), r*numpy.sin(theta)
27 return r*numpy.cos(theta), r*numpy.sin(theta)
28
28
29
29
30 def km2deg(km):
30 def km2deg(km):
31 '''
31 '''
32 Convert distance in km to degrees
32 Convert distance in km to degrees
33 '''
33 '''
34
34
35 return numpy.rad2deg(km/EARTH_RADIUS)
35 return numpy.rad2deg(km/EARTH_RADIUS)
36
36
37
37
38 class SpectraPlot(Plot):
38 class SpectraPlot(Plot):
39 '''
39 '''
40 Plot for Spectra data
40 Plot for Spectra data
41 '''
41 '''
42
42
43 CODE = 'spc'
43 CODE = 'spc'
44 colormap = 'jro'
44 colormap = 'jro'
45
45
46 def setup(self):
46 def setup(self):
47 self.nplots = len(self.data.channels)
47 self.nplots = len(self.data.channels)
48 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
48 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
49 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
49 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
50 self.width = 3.4 * self.ncols
50 self.width = 3.4 * self.ncols
51 self.height = 3 * self.nrows
51 self.height = 3 * self.nrows
52 self.cb_label = 'dB'
52 self.cb_label = 'dB'
53 if self.showprofile:
53 if self.showprofile:
54 self.width += 0.8 * self.ncols
54 self.width += 0.8 * self.ncols
55
55
56 self.ylabel = 'Range [km]'
56 self.ylabel = 'Range [km]'
57
57
58 def plot(self):
58 def plot(self):
59 if self.xaxis == "frequency":
59 if self.xaxis == "frequency":
60 x = self.data.xrange[0]
60 x = self.data.xrange[0]
61 self.xlabel = "Frequency (kHz)"
61 self.xlabel = "Frequency (kHz)"
62 elif self.xaxis == "time":
62 elif self.xaxis == "time":
63 x = self.data.xrange[1]
63 x = self.data.xrange[1]
64 self.xlabel = "Time (ms)"
64 self.xlabel = "Time (ms)"
65 else:
65 else:
66 x = self.data.xrange[2]
66 x = self.data.xrange[2]
67 self.xlabel = "Velocity (m/s)"
67 self.xlabel = "Velocity (m/s)"
68
68
69 if self.CODE == 'spc_mean':
69 if self.CODE == 'spc_mean':
70 x = self.data.xrange[2]
70 x = self.data.xrange[2]
71 self.xlabel = "Velocity (m/s)"
71 self.xlabel = "Velocity (m/s)"
72
72
73 self.titles = []
73 self.titles = []
74
74
75 y = self.data.heights
75 y = self.data.heights
76 self.y = y
76 self.y = y
77 z = self.data['spc']
77 z = self.data['spc']
78
78
79 for n, ax in enumerate(self.axes):
79 for n, ax in enumerate(self.axes):
80 noise = self.data['noise'][n][-1]
80 noise = self.data['noise'][n][-1]
81 if self.CODE == 'spc_mean':
81 if self.CODE == 'spc_mean':
82 mean = self.data['mean'][n][-1]
82 mean = self.data['mean'][n][-1]
83 if ax.firsttime:
83 if ax.firsttime:
84 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
84 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
85 self.xmin = self.xmin if self.xmin else -self.xmax
85 self.xmin = self.xmin if self.xmin else -self.xmax
86 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
86 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
87 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
87 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
88 ax.plt = ax.pcolormesh(x, y, z[n].T,
88 ax.plt = ax.pcolormesh(x, y, z[n].T,
89 vmin=self.zmin,
89 vmin=self.zmin,
90 vmax=self.zmax,
90 vmax=self.zmax,
91 cmap=plt.get_cmap(self.colormap)
91 cmap=plt.get_cmap(self.colormap)
92 )
92 )
93
93
94 if self.showprofile:
94 if self.showprofile:
95 ax.plt_profile = self.pf_axes[n].plot(
95 ax.plt_profile = self.pf_axes[n].plot(
96 self.data['rti'][n][-1], y)[0]
96 self.data['rti'][n][-1], y)[0]
97 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
97 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
98 color="k", linestyle="dashed", lw=1)[0]
98 color="k", linestyle="dashed", lw=1)[0]
99 if self.CODE == 'spc_mean':
99 if self.CODE == 'spc_mean':
100 ax.plt_mean = ax.plot(mean, y, color='k')[0]
100 ax.plt_mean = ax.plot(mean, y, color='k')[0]
101 else:
101 else:
102 ax.plt.set_array(z[n].T.ravel())
102 ax.plt.set_array(z[n].T.ravel())
103 if self.showprofile:
103 if self.showprofile:
104 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
104 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
105 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
105 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
106 if self.CODE == 'spc_mean':
106 if self.CODE == 'spc_mean':
107 ax.plt_mean.set_data(mean, y)
107 ax.plt_mean.set_data(mean, y)
108
108
109 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
109 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
110
110
111
111
112 class CrossSpectraPlot(Plot):
112 class CrossSpectraPlot(Plot):
113
113
114 CODE = 'cspc'
114 CODE = 'cspc'
115 colormap = 'jet'
115 colormap = 'jet'
116 zmin_coh = None
116 zmin_coh = None
117 zmax_coh = None
117 zmax_coh = None
118 zmin_phase = None
118 zmin_phase = None
119 zmax_phase = None
119 zmax_phase = None
120
120
121 def setup(self):
121 def setup(self):
122
122
123 self.ncols = 4
123 self.ncols = 4
124 self.nrows = len(self.data.pairs)
124 self.nrows = len(self.data.pairs)
125 self.nplots = self.nrows * 4
125 self.nplots = self.nrows * 4
126 self.width = 3.4 * self.ncols
126 self.width = 3.4 * self.ncols
127 self.height = 3 * self.nrows
127 self.height = 3 * self.nrows
128 self.ylabel = 'Range [km]'
128 self.ylabel = 'Range [km]'
129 self.showprofile = False
129 self.showprofile = False
130
130
131 def plot(self):
131 def plot(self):
132
132
133 if self.xaxis == "frequency":
133 if self.xaxis == "frequency":
134 x = self.data.xrange[0]
134 x = self.data.xrange[0]
135 self.xlabel = "Frequency (kHz)"
135 self.xlabel = "Frequency (kHz)"
136 elif self.xaxis == "time":
136 elif self.xaxis == "time":
137 x = self.data.xrange[1]
137 x = self.data.xrange[1]
138 self.xlabel = "Time (ms)"
138 self.xlabel = "Time (ms)"
139 else:
139 else:
140 x = self.data.xrange[2]
140 x = self.data.xrange[2]
141 self.xlabel = "Velocity (m/s)"
141 self.xlabel = "Velocity (m/s)"
142
142
143 self.titles = []
143 self.titles = []
144
144
145 y = self.data.heights
145 y = self.data.heights
146 self.y = y
146 self.y = y
147 spc = self.data['spc']
147 spc = self.data['spc']
148 cspc = self.data['cspc']
148 cspc = self.data['cspc']
149
149
150 for n in range(self.nrows):
150 for n in range(self.nrows):
151 noise = self.data['noise'][n][-1]
151 noise = self.data['noise'][n][-1]
152 pair = self.data.pairs[n]
152 pair = self.data.pairs[n]
153 ax = self.axes[4 * n]
153 ax = self.axes[4 * n]
154 spc0 = 10.*numpy.log10(spc[pair[0]]/self.data.factor)
154 spc0 = 10.*numpy.log10(spc[pair[0]]/self.data.factor)
155 if ax.firsttime:
155 if ax.firsttime:
156 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
156 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
157 self.xmin = self.xmin if self.xmin else -self.xmax
157 self.xmin = self.xmin if self.xmin else -self.xmax
158 self.zmin = self.zmin if self.zmin else numpy.nanmin(spc)
158 self.zmin = self.zmin if self.zmin else numpy.nanmin(spc)
159 self.zmax = self.zmax if self.zmax else numpy.nanmax(spc)
159 self.zmax = self.zmax if self.zmax else numpy.nanmax(spc)
160 ax.plt = ax.pcolormesh(x , y , spc0.T,
160 ax.plt = ax.pcolormesh(x , y , spc0.T,
161 vmin=self.zmin,
161 vmin=self.zmin,
162 vmax=self.zmax,
162 vmax=self.zmax,
163 cmap=plt.get_cmap(self.colormap)
163 cmap=plt.get_cmap(self.colormap)
164 )
164 )
165 else:
165 else:
166 ax.plt.set_array(spc0.T.ravel())
166 ax.plt.set_array(spc0.T.ravel())
167 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise))
167 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise))
168
168
169 ax = self.axes[4 * n + 1]
169 ax = self.axes[4 * n + 1]
170 spc1 = 10.*numpy.log10(spc[pair[1]]/self.data.factor)
170 spc1 = 10.*numpy.log10(spc[pair[1]]/self.data.factor)
171 if ax.firsttime:
171 if ax.firsttime:
172 ax.plt = ax.pcolormesh(x , y, spc1.T,
172 ax.plt = ax.pcolormesh(x , y, spc1.T,
173 vmin=self.zmin,
173 vmin=self.zmin,
174 vmax=self.zmax,
174 vmax=self.zmax,
175 cmap=plt.get_cmap(self.colormap)
175 cmap=plt.get_cmap(self.colormap)
176 )
176 )
177 else:
177 else:
178 ax.plt.set_array(spc1.T.ravel())
178 ax.plt.set_array(spc1.T.ravel())
179 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise))
179 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise))
180
180
181 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
181 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
182 coh = numpy.abs(out)
182 coh = numpy.abs(out)
183 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
183 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
184
184
185 ax = self.axes[4 * n + 2]
185 ax = self.axes[4 * n + 2]
186 if ax.firsttime:
186 if ax.firsttime:
187 ax.plt = ax.pcolormesh(x, y, coh.T,
187 ax.plt = ax.pcolormesh(x, y, coh.T,
188 vmin=0,
188 vmin=0,
189 vmax=1,
189 vmax=1,
190 cmap=plt.get_cmap(self.colormap_coh)
190 cmap=plt.get_cmap(self.colormap_coh)
191 )
191 )
192 else:
192 else:
193 ax.plt.set_array(coh.T.ravel())
193 ax.plt.set_array(coh.T.ravel())
194 self.titles.append(
194 self.titles.append(
195 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
195 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
196
196
197 ax = self.axes[4 * n + 3]
197 ax = self.axes[4 * n + 3]
198 if ax.firsttime:
198 if ax.firsttime:
199 ax.plt = ax.pcolormesh(x, y, phase.T,
199 ax.plt = ax.pcolormesh(x, y, phase.T,
200 vmin=-180,
200 vmin=-180,
201 vmax=180,
201 vmax=180,
202 cmap=plt.get_cmap(self.colormap_phase)
202 cmap=plt.get_cmap(self.colormap_phase)
203 )
203 )
204 else:
204 else:
205 ax.plt.set_array(phase.T.ravel())
205 ax.plt.set_array(phase.T.ravel())
206 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
206 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
207
207
208
208
209 class SpectraMeanPlot(SpectraPlot):
209 class SpectraMeanPlot(SpectraPlot):
210 '''
210 '''
211 Plot for Spectra and Mean
211 Plot for Spectra and Mean
212 '''
212 '''
213 CODE = 'spc_mean'
213 CODE = 'spc_mean'
214 colormap = 'jro'
214 colormap = 'jro'
215
215
216
216
217 class RTIPlot(Plot):
217 class RTIPlot(Plot):
218 '''
218 '''
219 Plot for RTI data
219 Plot for RTI data
220 '''
220 '''
221
221
222 CODE = 'rti'
222 CODE = 'rti'
223 colormap = 'jro'
223 colormap = 'jro'
224
224
225 def setup(self):
225 def setup(self):
226 self.xaxis = 'time'
226 self.xaxis = 'time'
227 self.ncols = 1
227 self.ncols = 1
228 self.nrows = len(self.data.channels)
228 self.nrows = len(self.data.channels)
229 self.nplots = len(self.data.channels)
229 self.nplots = len(self.data.channels)
230 self.ylabel = 'Range [km]'
230 self.ylabel = 'Range [km]'
231 self.cb_label = 'dB'
231 self.cb_label = 'dB'
232 self.titles = ['{} Channel {}'.format(
232 self.titles = ['{} Channel {}'.format(
233 self.CODE.upper(), x) for x in range(self.nrows)]
233 self.CODE.upper(), x) for x in range(self.nrows)]
234
234
235 def plot(self):
235 def plot(self):
236 self.x = self.data.times
236 self.x = self.data.times
237 self.y = self.data.heights
237 self.y = self.data.heights
238 self.z = self.data[self.CODE]
238 self.z = self.data[self.CODE]
239 self.z = numpy.ma.masked_invalid(self.z)
239 self.z = numpy.ma.masked_invalid(self.z)
240
240
241 if self.decimation is None:
241 if self.decimation is None:
242 x, y, z = self.fill_gaps(self.x, self.y, self.z)
242 x, y, z = self.fill_gaps(self.x, self.y, self.z)
243 else:
243 else:
244 x, y, z = self.fill_gaps(*self.decimate())
244 x, y, z = self.fill_gaps(*self.decimate())
245
245
246 for n, ax in enumerate(self.axes):
246 for n, ax in enumerate(self.axes):
247 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
247 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
248 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
248 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
249 if ax.firsttime:
249 if ax.firsttime:
250 ax.plt = ax.pcolormesh(x, y, z[n].T,
250 ax.plt = ax.pcolormesh(x, y, z[n].T,
251 vmin=self.zmin,
251 vmin=self.zmin,
252 vmax=self.zmax,
252 vmax=self.zmax,
253 cmap=plt.get_cmap(self.colormap)
253 cmap=plt.get_cmap(self.colormap)
254 )
254 )
255 if self.showprofile:
255 if self.showprofile:
256 ax.plot_profile = self.pf_axes[n].plot(
256 ax.plot_profile = self.pf_axes[n].plot(
257 self.data['rti'][n][-1], self.y)[0]
257 self.data['rti'][n][-1], self.y)[0]
258 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
258 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
259 color="k", linestyle="dashed", lw=1)[0]
259 color="k", linestyle="dashed", lw=1)[0]
260 else:
260 else:
261 ax.collections.remove(ax.collections[0])
261 ax.collections.remove(ax.collections[0])
262 ax.plt = ax.pcolormesh(x, y, z[n].T,
262 ax.plt = ax.pcolormesh(x, y, z[n].T,
263 vmin=self.zmin,
263 vmin=self.zmin,
264 vmax=self.zmax,
264 vmax=self.zmax,
265 cmap=plt.get_cmap(self.colormap)
265 cmap=plt.get_cmap(self.colormap)
266 )
266 )
267 if self.showprofile:
267 if self.showprofile:
268 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
268 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
269 ax.plot_noise.set_data(numpy.repeat(
269 ax.plot_noise.set_data(numpy.repeat(
270 self.data['noise'][n][-1], len(self.y)), self.y)
270 self.data['noise'][n][-1], len(self.y)), self.y)
271
271
272
272
273 class CoherencePlot(RTIPlot):
273 class CoherencePlot(RTIPlot):
274 '''
274 '''
275 Plot for Coherence data
275 Plot for Coherence data
276 '''
276 '''
277
277
278 CODE = 'coh'
278 CODE = 'coh'
279
279
280 def setup(self):
280 def setup(self):
281 self.xaxis = 'time'
281 self.xaxis = 'time'
282 self.ncols = 1
282 self.ncols = 1
283 self.nrows = len(self.data.pairs)
283 self.nrows = len(self.data.pairs)
284 self.nplots = len(self.data.pairs)
284 self.nplots = len(self.data.pairs)
285 self.ylabel = 'Range [km]'
285 self.ylabel = 'Range [km]'
286 if self.CODE == 'coh':
286 if self.CODE == 'coh':
287 self.cb_label = ''
287 self.cb_label = ''
288 self.titles = [
288 self.titles = [
289 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
289 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
290 else:
290 else:
291 self.cb_label = 'Degrees'
291 self.cb_label = 'Degrees'
292 self.titles = [
292 self.titles = [
293 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
293 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
294
294
295
295
296 class PhasePlot(CoherencePlot):
296 class PhasePlot(CoherencePlot):
297 '''
297 '''
298 Plot for Phase map data
298 Plot for Phase map data
299 '''
299 '''
300
300
301 CODE = 'phase'
301 CODE = 'phase'
302 colormap = 'seismic'
302 colormap = 'seismic'
303
303
304
304
305 class NoisePlot(Plot):
305 class NoisePlot(Plot):
306 '''
306 '''
307 Plot for noise
307 Plot for noise
308 '''
308 '''
309
309
310 CODE = 'noise'
310 CODE = 'noise'
311
311
312 def setup(self):
312 def setup(self):
313 self.xaxis = 'time'
313 self.xaxis = 'time'
314 self.ncols = 1
314 self.ncols = 1
315 self.nrows = 1
315 self.nrows = 1
316 self.nplots = 1
316 self.nplots = 1
317 self.ylabel = 'Intensity [dB]'
317 self.ylabel = 'Intensity [dB]'
318 self.titles = ['Noise']
318 self.titles = ['Noise']
319 self.colorbar = False
319 self.colorbar = False
320
320
321 def plot(self):
321 def plot(self):
322
322
323 x = self.data.times
323 x = self.data.times
324 xmin = self.data.min_time
324 xmin = self.data.min_time
325 xmax = xmin + self.xrange * 60 * 60
325 xmax = xmin + self.xrange * 60 * 60
326 Y = self.data[self.CODE]
326 Y = self.data[self.CODE]
327
327
328 if self.axes[0].firsttime:
328 if self.axes[0].firsttime:
329 for ch in self.data.channels:
329 for ch in self.data.channels:
330 y = Y[ch]
330 y = Y[ch]
331 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
331 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
332 plt.legend()
332 plt.legend()
333 else:
333 else:
334 for ch in self.data.channels:
334 for ch in self.data.channels:
335 y = Y[ch]
335 y = Y[ch]
336 self.axes[0].lines[ch].set_data(x, y)
336 self.axes[0].lines[ch].set_data(x, y)
337
337
338 self.ymin = numpy.nanmin(Y) - 5
338 self.ymin = numpy.nanmin(Y) - 5
339 self.ymax = numpy.nanmax(Y) + 5
339 self.ymax = numpy.nanmax(Y) + 5
340
340
341
341
342 class SnrPlot(RTIPlot):
342 class SnrPlot(RTIPlot):
343 '''
343 '''
344 Plot for SNR Data
344 Plot for SNR Data
345 '''
345 '''
346
346
347 CODE = 'snr'
347 CODE = 'snr'
348 colormap = 'jet'
348 colormap = 'jet'
349
349
350
350
351 class DopplerPlot(RTIPlot):
351 class DopplerPlot(RTIPlot):
352 '''
352 '''
353 Plot for DOPPLER Data
353 Plot for DOPPLER Data
354 '''
354 '''
355
355
356 CODE = 'dop'
356 CODE = 'dop'
357 colormap = 'jet'
357 colormap = 'jet'
358
358
359
359
360 class SkyMapPlot(Plot):
360 class SkyMapPlot(Plot):
361 '''
361 '''
362 Plot for meteors detection data
362 Plot for meteors detection data
363 '''
363 '''
364
364
365 CODE = 'param'
365 CODE = 'param'
366
366
367 def setup(self):
367 def setup(self):
368
368
369 self.ncols = 1
369 self.ncols = 1
370 self.nrows = 1
370 self.nrows = 1
371 self.width = 7.2
371 self.width = 7.2
372 self.height = 7.2
372 self.height = 7.2
373 self.nplots = 1
373 self.nplots = 1
374 self.xlabel = 'Zonal Zenith Angle (deg)'
374 self.xlabel = 'Zonal Zenith Angle (deg)'
375 self.ylabel = 'Meridional Zenith Angle (deg)'
375 self.ylabel = 'Meridional Zenith Angle (deg)'
376 self.polar = True
376 self.polar = True
377 self.ymin = -180
377 self.ymin = -180
378 self.ymax = 180
378 self.ymax = 180
379 self.colorbar = False
379 self.colorbar = False
380
380
381 def plot(self):
381 def plot(self):
382
382
383 arrayParameters = numpy.concatenate(self.data['param'])
383 arrayParameters = numpy.concatenate(self.data['param'])
384 error = arrayParameters[:, -1]
384 error = arrayParameters[:, -1]
385 indValid = numpy.where(error == 0)[0]
385 indValid = numpy.where(error == 0)[0]
386 finalMeteor = arrayParameters[indValid, :]
386 finalMeteor = arrayParameters[indValid, :]
387 finalAzimuth = finalMeteor[:, 3]
387 finalAzimuth = finalMeteor[:, 3]
388 finalZenith = finalMeteor[:, 4]
388 finalZenith = finalMeteor[:, 4]
389
389
390 x = finalAzimuth * numpy.pi / 180
390 x = finalAzimuth * numpy.pi / 180
391 y = finalZenith
391 y = finalZenith
392
392
393 ax = self.axes[0]
393 ax = self.axes[0]
394
394
395 if ax.firsttime:
395 if ax.firsttime:
396 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
396 ax.plot = ax.plot(x, y, 'bo', markersize=5)[0]
397 else:
397 else:
398 ax.plot.set_data(x, y)
398 ax.plot.set_data(x, y)
399
399
400 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
400 dt1 = self.getDateTime(self.data.min_time).strftime('%y/%m/%d %H:%M:%S')
401 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
401 dt2 = self.getDateTime(self.data.max_time).strftime('%y/%m/%d %H:%M:%S')
402 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
402 title = 'Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n' % (dt1,
403 dt2,
403 dt2,
404 len(x))
404 len(x))
405 self.titles[0] = title
405 self.titles[0] = title
406
406
407
407
408 class ParametersPlot(RTIPlot):
408 class ParametersPlot(RTIPlot):
409 '''
409 '''
410 Plot for data_param object
410 Plot for data_param object
411 '''
411 '''
412
412
413 CODE = 'param'
413 CODE = 'param'
414 colormap = 'seismic'
414 colormap = 'seismic'
415
415
416 def setup(self):
416 def setup(self):
417 self.xaxis = 'time'
417 self.xaxis = 'time'
418 self.ncols = 1
418 self.ncols = 1
419 self.nrows = self.data.shape(self.CODE)[0]
419 self.nrows = self.data.shape(self.CODE)[0]
420 self.nplots = self.nrows
420 self.nplots = self.nrows
421 if self.showSNR:
421 if self.showSNR:
422 self.nrows += 1
422 self.nrows += 1
423 self.nplots += 1
423 self.nplots += 1
424
424
425 self.ylabel = 'Height [km]'
425 self.ylabel = 'Height [km]'
426 if not self.titles:
426 if not self.titles:
427 self.titles = self.data.parameters \
427 self.titles = self.data.parameters \
428 if self.data.parameters else ['Param {}'.format(x) for x in range(self.nrows)]
428 if self.data.parameters else ['Param {}'.format(x) for x in range(self.nrows)]
429 if self.showSNR:
429 if self.showSNR:
430 self.titles.append('SNR')
430 self.titles.append('SNR')
431
431
432 def plot(self):
432 def plot(self):
433 self.data.normalize_heights()
433 self.data.normalize_heights()
434 self.x = self.data.times
434 self.x = self.data.times
435 self.y = self.data.heights
435 self.y = self.data.heights
436 if self.showSNR:
436 if self.showSNR:
437 self.z = numpy.concatenate(
437 self.z = numpy.concatenate(
438 (self.data[self.CODE], self.data['snr'])
438 (self.data[self.CODE], self.data['snr'])
439 )
439 )
440 else:
440 else:
441 self.z = self.data[self.CODE]
441 self.z = self.data[self.CODE]
442
442
443 self.z = numpy.ma.masked_invalid(self.z)
443 self.z = numpy.ma.masked_invalid(self.z)
444
444
445 if self.decimation is None:
445 if self.decimation is None:
446 x, y, z = self.fill_gaps(self.x, self.y, self.z)
446 x, y, z = self.fill_gaps(self.x, self.y, self.z)
447 else:
447 else:
448 x, y, z = self.fill_gaps(*self.decimate())
448 x, y, z = self.fill_gaps(*self.decimate())
449
449
450 for n, ax in enumerate(self.axes):
450 for n, ax in enumerate(self.axes):
451
451
452 self.zmax = self.zmax if self.zmax is not None else numpy.max(
452 self.zmax = self.zmax if self.zmax is not None else numpy.max(
453 self.z[n])
453 self.z[n])
454 self.zmin = self.zmin if self.zmin is not None else numpy.min(
454 self.zmin = self.zmin if self.zmin is not None else numpy.min(
455 self.z[n])
455 self.z[n])
456
456
457 if ax.firsttime:
457 if ax.firsttime:
458 if self.zlimits is not None:
458 if self.zlimits is not None:
459 self.zmin, self.zmax = self.zlimits[n]
459 self.zmin, self.zmax = self.zlimits[n]
460
460
461 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
461 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
462 vmin=self.zmin,
462 vmin=self.zmin,
463 vmax=self.zmax,
463 vmax=self.zmax,
464 cmap=self.cmaps[n]
464 cmap=self.cmaps[n]
465 )
465 )
466 else:
466 else:
467 if self.zlimits is not None:
467 if self.zlimits is not None:
468 self.zmin, self.zmax = self.zlimits[n]
468 self.zmin, self.zmax = self.zlimits[n]
469 ax.collections.remove(ax.collections[0])
469 ax.collections.remove(ax.collections[0])
470 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
470 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
471 vmin=self.zmin,
471 vmin=self.zmin,
472 vmax=self.zmax,
472 vmax=self.zmax,
473 cmap=self.cmaps[n]
473 cmap=self.cmaps[n]
474 )
474 )
475
475
476
476
477 class OutputPlot(ParametersPlot):
477 class OutputPlot(ParametersPlot):
478 '''
478 '''
479 Plot data_output object
479 Plot data_output object
480 '''
480 '''
481
481
482 CODE = 'output'
482 CODE = 'output'
483 colormap = 'seismic'
483 colormap = 'seismic'
484
484
485
485
486 class PolarMapPlot(Plot):
486 class PolarMapPlot(Plot):
487 '''
487 '''
488 Plot for weather radar
488 Plot for weather radar
489 '''
489 '''
490
490
491 CODE = 'param'
491 CODE = 'param'
492 colormap = 'seismic'
492 colormap = 'seismic'
493
493
494 def setup(self):
494 def setup(self):
495 self.ncols = 1
495 self.ncols = 1
496 self.nrows = 1
496 self.nrows = 1
497 self.width = 9
497 self.width = 9
498 self.height = 8
498 self.height = 8
499 self.mode = self.data.meta['mode']
499 self.mode = self.data.meta['mode']
500 if self.channels is not None:
500 if self.channels is not None:
501 self.nplots = len(self.channels)
501 self.nplots = len(self.channels)
502 self.nrows = len(self.channels)
502 self.nrows = len(self.channels)
503 else:
503 else:
504 self.nplots = self.data.shape(self.CODE)[0]
504 self.nplots = self.data.shape(self.CODE)[0]
505 self.nrows = self.nplots
505 self.nrows = self.nplots
506 self.channels = list(range(self.nplots))
506 self.channels = list(range(self.nplots))
507 if self.mode == 'E':
507 if self.mode == 'E':
508 self.xlabel = 'Longitude'
508 self.xlabel = 'Longitude'
509 self.ylabel = 'Latitude'
509 self.ylabel = 'Latitude'
510 else:
510 else:
511 self.xlabel = 'Range (km)'
511 self.xlabel = 'Range (km)'
512 self.ylabel = 'Height (km)'
512 self.ylabel = 'Height (km)'
513 self.bgcolor = 'white'
513 self.bgcolor = 'white'
514 self.cb_labels = self.data.meta['units']
514 self.cb_labels = self.data.meta['units']
515 self.lat = self.data.meta['latitude']
515 self.lat = self.data.meta['latitude']
516 self.lon = self.data.meta['longitude']
516 self.lon = self.data.meta['longitude']
517 self.xmin, self.xmax = float(
517 self.xmin, self.xmax = float(
518 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
518 km2deg(self.xmin) + self.lon), float(km2deg(self.xmax) + self.lon)
519 self.ymin, self.ymax = float(
519 self.ymin, self.ymax = float(
520 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
520 km2deg(self.ymin) + self.lat), float(km2deg(self.ymax) + self.lat)
521 # self.polar = True
521 # self.polar = True
522
522
523 def plot(self):
523 def plot(self):
524
524
525 for n, ax in enumerate(self.axes):
525 for n, ax in enumerate(self.axes):
526 data = self.data['param'][self.channels[n]]
526 data = self.data['param'][self.channels[n]]
527
527
528 zeniths = numpy.linspace(
528 zeniths = numpy.linspace(
529 0, self.data.meta['max_range'], data.shape[1])
529 0, self.data.meta['max_range'], data.shape[1])
530 if self.mode == 'E':
530 if self.mode == 'E':
531 azimuths = -numpy.radians(self.data.heights)+numpy.pi/2
531 azimuths = -numpy.radians(self.data.heights)+numpy.pi/2
532 r, theta = numpy.meshgrid(zeniths, azimuths)
532 r, theta = numpy.meshgrid(zeniths, azimuths)
533 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
533 x, y = r*numpy.cos(theta)*numpy.cos(numpy.radians(self.data.meta['elevation'])), r*numpy.sin(
534 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
534 theta)*numpy.cos(numpy.radians(self.data.meta['elevation']))
535 x = km2deg(x) + self.lon
535 x = km2deg(x) + self.lon
536 y = km2deg(y) + self.lat
536 y = km2deg(y) + self.lat
537 else:
537 else:
538 azimuths = numpy.radians(self.data.heights)
538 azimuths = numpy.radians(self.data.heights)
539 r, theta = numpy.meshgrid(zeniths, azimuths)
539 r, theta = numpy.meshgrid(zeniths, azimuths)
540 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
540 x, y = r*numpy.cos(theta), r*numpy.sin(theta)
541 self.y = zeniths
541 self.y = zeniths
542
542
543 if ax.firsttime:
543 if ax.firsttime:
544 if self.zlimits is not None:
544 if self.zlimits is not None:
545 self.zmin, self.zmax = self.zlimits[n]
545 self.zmin, self.zmax = self.zlimits[n]
546 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
546 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
547 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
547 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
548 vmin=self.zmin,
548 vmin=self.zmin,
549 vmax=self.zmax,
549 vmax=self.zmax,
550 cmap=self.cmaps[n])
550 cmap=self.cmaps[n])
551 else:
551 else:
552 if self.zlimits is not None:
552 if self.zlimits is not None:
553 self.zmin, self.zmax = self.zlimits[n]
553 self.zmin, self.zmax = self.zlimits[n]
554 ax.collections.remove(ax.collections[0])
554 ax.collections.remove(ax.collections[0])
555 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
555 ax.plt = ax.pcolormesh( # r, theta, numpy.ma.array(data, mask=numpy.isnan(data)),
556 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
556 x, y, numpy.ma.array(data, mask=numpy.isnan(data)),
557 vmin=self.zmin,
557 vmin=self.zmin,
558 vmax=self.zmax,
558 vmax=self.zmax,
559 cmap=self.cmaps[n])
559 cmap=self.cmaps[n])
560
560
561 if self.mode == 'A':
561 if self.mode == 'A':
562 continue
562 continue
563
563
564 # plot district names
564 # plot district names
565 f = open('/data/workspace/schain_scripts/distrito.csv')
565 f = open('/data/workspace/schain_scripts/distrito.csv')
566 for line in f:
566 for line in f:
567 label, lon, lat = [s.strip() for s in line.split(',') if s]
567 label, lon, lat = [s.strip() for s in line.split(',') if s]
568 lat = float(lat)
568 lat = float(lat)
569 lon = float(lon)
569 lon = float(lon)
570 # ax.plot(lon, lat, '.b', ms=2)
570 # ax.plot(lon, lat, '.b', ms=2)
571 ax.text(lon, lat, label.decode('utf8'), ha='center',
571 ax.text(lon, lat, label.decode('utf8'), ha='center',
572 va='bottom', size='8', color='black')
572 va='bottom', size='8', color='black')
573
573
574 # plot limites
574 # plot limites
575 limites = []
575 limites = []
576 tmp = []
576 tmp = []
577 for line in open('/data/workspace/schain_scripts/lima.csv'):
577 for line in open('/data/workspace/schain_scripts/lima.csv'):
578 if '#' in line:
578 if '#' in line:
579 if tmp:
579 if tmp:
580 limites.append(tmp)
580 limites.append(tmp)
581 tmp = []
581 tmp = []
582 continue
582 continue
583 values = line.strip().split(',')
583 values = line.strip().split(',')
584 tmp.append((float(values[0]), float(values[1])))
584 tmp.append((float(values[0]), float(values[1])))
585 for points in limites:
585 for points in limites:
586 ax.add_patch(
586 ax.add_patch(
587 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
587 Polygon(points, ec='k', fc='none', ls='--', lw=0.5))
588
588
589 # plot Cuencas
589 # plot Cuencas
590 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
590 for cuenca in ('rimac', 'lurin', 'mala', 'chillon', 'chilca', 'chancay-huaral'):
591 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
591 f = open('/data/workspace/schain_scripts/{}.csv'.format(cuenca))
592 values = [line.strip().split(',') for line in f]
592 values = [line.strip().split(',') for line in f]
593 points = [(float(s[0]), float(s[1])) for s in values]
593 points = [(float(s[0]), float(s[1])) for s in values]
594 ax.add_patch(Polygon(points, ec='b', fc='none'))
594 ax.add_patch(Polygon(points, ec='b', fc='none'))
595
595
596 # plot grid
596 # plot grid
597 for r in (15, 30, 45, 60):
597 for r in (15, 30, 45, 60):
598 ax.add_artist(plt.Circle((self.lon, self.lat),
598 ax.add_artist(plt.Circle((self.lon, self.lat),
599 km2deg(r), color='0.6', fill=False, lw=0.2))
599 km2deg(r), color='0.6', fill=False, lw=0.2))
600 ax.text(
600 ax.text(
601 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
601 self.lon + (km2deg(r))*numpy.cos(60*numpy.pi/180),
602 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
602 self.lat + (km2deg(r))*numpy.sin(60*numpy.pi/180),
603 '{}km'.format(r),
603 '{}km'.format(r),
604 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
604 ha='center', va='bottom', size='8', color='0.6', weight='heavy')
605
605
606 if self.mode == 'E':
606 if self.mode == 'E':
607 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
607 title = 'El={}$^\circ$'.format(self.data.meta['elevation'])
608 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
608 label = 'E{:02d}'.format(int(self.data.meta['elevation']))
609 else:
609 else:
610 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
610 title = 'Az={}$^\circ$'.format(self.data.meta['azimuth'])
611 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
611 label = 'A{:02d}'.format(int(self.data.meta['azimuth']))
612
612
613 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
613 self.save_labels = ['{}-{}'.format(lbl, label) for lbl in self.labels]
614 self.titles = ['{} {}'.format(
614 self.titles = ['{} {}'.format(
615 self.data.parameters[x], title) for x in self.channels]
615 self.data.parameters[x], title) for x in self.channels]
616 No newline at end of file
616
617 class ScopePlot(Plot):
618
619 '''
620 Plot for Scope
621 '''
622
623 CODE = 'scope'
624
625 def setup(self):
626
627 self.xaxis = 'Range (Km)'
628 self.ncols = 1
629 self.nrows = 1
630 self.nplots = 1
631 self.ylabel = 'Intensity [dB]'
632 self.titles = ['Scope']
633 self.colorbar = False
634 colspan = 3
635 rowspan = 1
636
637 def plot_iq(self, x, y, channelIndexList, thisDatetime, wintitle):
638
639 yreal = y[channelIndexList,:].real
640 yimag = y[channelIndexList,:].imag
641 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
642 self.xlabel = "Range (Km)"
643 self.ylabel = "Intensity - IQ"
644
645 self.y = yreal
646 self.x = x
647 self.xmin = min(x)
648 self.xmax = max(x)
649
650
651 self.titles[0] = title
652
653 for i,ax in enumerate(self.axes):
654 title = "Channel %d" %(i)
655 if ax.firsttime:
656 ax.plt_r = ax.plot(x, yreal[i,:], color='b')[0]
657 ax.plt_i = ax.plot(x, yimag[i,:], color='r')[0]
658 else:
659 #pass
660 ax.plt_r.set_data(x, yreal[i,:])
661 ax.plt_i.set_data(x, yimag[i,:])
662
663 def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle):
664 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
665 yreal = y.real
666 self.y = yreal
667 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
668 self.xlabel = "Range (Km)"
669 self.ylabel = "Intensity"
670 self.xmin = min(x)
671 self.xmax = max(x)
672
673
674 self.titles[0] = title
675
676 for i,ax in enumerate(self.axes):
677 title = "Channel %d" %(i)
678
679 ychannel = yreal[i,:]
680
681 if ax.firsttime:
682 ax.plt_r = ax.plot(x, ychannel)[0]
683 else:
684 #pass
685 ax.plt_r.set_data(x, ychannel)
686
687
688 def plot(self):
689
690 if self.channels:
691 channels = self.channels
692 else:
693 channels = self.data.channels
694
695
696
697 thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
698
699 scope = self.data['scope']
700
701
702 if self.data.flagDataAsBlock:
703
704 for i in range(self.data.nProfiles):
705
706 wintitle1 = " [Profile = %d] " %i
707
708 if self.type == "power":
709 self.plot_power(self.data.heights,
710 scope[:,i,:],
711 channels,
712 thisDatetime,
713 wintitle1
714 )
715
716 if self.type == "iq":
717 self.plot_iq(self.data.heights,
718 scope[:,i,:],
719 channels,
720 thisDatetime,
721 wintitle1
722 )
723
724
725
726
727
728 else:
729 wintitle = " [Profile = %d] " %self.data.profileIndex
730
731 if self.type == "power":
732 self.plot_power(self.data.heights,
733 scope,
734 channels,
735 thisDatetime,
736 wintitle
737 )
738
739 if self.type == "iq":
740 self.plot_iq(self.data.heights,
741 scope,
742 channels,
743 thisDatetime,
744 wintitle
745 )
746
747
748 No newline at end of file
General Comments 0
You need to be logged in to leave comments. Login now