##// END OF EJS Templates
schain
RSS
Drifts tested
joabAM -
r1740:7f5b085e2124
parent child
Show More
Side by Side Unified
Context file:
r1740:7f5b085e2124 Loading diff...:

The requested changes are too big and content was truncated. Show full diff

Show file before | Show file after | Hide comments
@@ -1,1156 +1,1159
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 # All rights reserved.
2 # All rights reserved.
3 #
3 #
4 # Distributed under the terms of the BSD 3-clause license.
4 # Distributed under the terms of the BSD 3-clause license.
5 """Definition of diferent Data objects for different types of data
5 """Definition of diferent Data objects for different types of data
6
6
7 Here you will find the diferent data objects for the different types
7 Here you will find the diferent data objects for the different types
8 of data, this data objects must be used as dataIn or dataOut objects in
8 of data, this data objects must be used as dataIn or dataOut objects in
9 processing units and operations. Currently the supported data objects are:
9 processing units and operations. Currently the supported data objects are:
10 Voltage, Spectra, SpectraHeis, Fits, Correlation and Parameters
10 Voltage, Spectra, SpectraHeis, Fits, Correlation and Parameters
11 """
11 """
12
12
13 import copy
13 import copy
14 import numpy
14 import numpy
15 import datetime
15 import datetime
16 import json
16 import json
17
17
18 import schainpy.admin
18 import schainpy.admin
19 from schainpy.utils import log
19 from schainpy.utils import log
20 from .jroheaderIO import SystemHeader, RadarControllerHeader,ProcessingHeader
20 from .jroheaderIO import SystemHeader, RadarControllerHeader,ProcessingHeader
21 from schainpy.model.data import _noise
21 from schainpy.model.data import _noise
22 SPEED_OF_LIGHT = 3e8
22 SPEED_OF_LIGHT = 3e8
23
23
24
24
25 def getNumpyDtype(dataTypeCode):
25 def getNumpyDtype(dataTypeCode):
26
26
27 if dataTypeCode == 0:
27 if dataTypeCode == 0:
28 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
28 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
29 elif dataTypeCode == 1:
29 elif dataTypeCode == 1:
30 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
30 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
31 elif dataTypeCode == 2:
31 elif dataTypeCode == 2:
32 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
32 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
33 elif dataTypeCode == 3:
33 elif dataTypeCode == 3:
34 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
34 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
35 elif dataTypeCode == 4:
35 elif dataTypeCode == 4:
36 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
36 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
37 elif dataTypeCode == 5:
37 elif dataTypeCode == 5:
38 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
38 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
39 else:
39 else:
40 raise ValueError('dataTypeCode was not defined')
40 raise ValueError('dataTypeCode was not defined')
41
41
42 return numpyDtype
42 return numpyDtype
43
43
44
44
45 def getDataTypeCode(numpyDtype):
45 def getDataTypeCode(numpyDtype):
46
46
47 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
47 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
48 datatype = 0
48 datatype = 0
49 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
49 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
50 datatype = 1
50 datatype = 1
51 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
51 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
52 datatype = 2
52 datatype = 2
53 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
53 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
54 datatype = 3
54 datatype = 3
55 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
55 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
56 datatype = 4
56 datatype = 4
57 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
57 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
58 datatype = 5
58 datatype = 5
59 else:
59 else:
60 datatype = None
60 datatype = None
61
61
62 return datatype
62 return datatype
63
63
64
64
65 def hildebrand_sekhon(data, navg):
65 def hildebrand_sekhon(data, navg):
66 """
66 """
67 This method is for the objective determination of the noise level in Doppler spectra. This
67 This method is for the objective determination of the noise level in Doppler spectra. This
68 implementation technique is based on the fact that the standard deviation of the spectral
68 implementation technique is based on the fact that the standard deviation of the spectral
69 densities is equal to the mean spectral density for white Gaussian noise
69 densities is equal to the mean spectral density for white Gaussian noise
70
70
71 Inputs:
71 Inputs:
72 Data : heights
72 Data : heights
73 navg : numbers of averages
73 navg : numbers of averages
74
74
75 Return:
75 Return:
76 mean : noise's level
76 mean : noise's level
77 """
77 """
78
78
79 sortdata = numpy.sort(data, axis=None)
79 sortdata = numpy.sort(data, axis=None)
80 '''
80 '''
81 lenOfData = len(sortdata)
81 lenOfData = len(sortdata)
82 nums_min = lenOfData*0.2
82 nums_min = lenOfData*0.2
83
83
84 if nums_min <= 5:
84 if nums_min <= 5:
85
85
86 nums_min = 5
86 nums_min = 5
87
87
88 sump = 0.
88 sump = 0.
89 sumq = 0.
89 sumq = 0.
90
90
91 j = 0
91 j = 0
92 cont = 1
92 cont = 1
93
93
94 while((cont == 1)and(j < lenOfData)):
94 while((cont == 1)and(j < lenOfData)):
95
95
96 sump += sortdata[j]
96 sump += sortdata[j]
97 sumq += sortdata[j]**2
97 sumq += sortdata[j]**2
98
98
99 if j > nums_min:
99 if j > nums_min:
100 rtest = float(j)/(j-1) + 1.0/navg
100 rtest = float(j)/(j-1) + 1.0/navg
101 if ((sumq*j) > (rtest*sump**2)):
101 if ((sumq*j) > (rtest*sump**2)):
102 j = j - 1
102 j = j - 1
103 sump = sump - sortdata[j]
103 sump = sump - sortdata[j]
104 sumq = sumq - sortdata[j]**2
104 sumq = sumq - sortdata[j]**2
105 cont = 0
105 cont = 0
106
106
107 j += 1
107 j += 1
108
108
109 lnoise = sump / j
109 lnoise = sump / j
110 '''
110 '''
111 return _noise.hildebrand_sekhon(sortdata, navg)
111 return _noise.hildebrand_sekhon(sortdata, navg)
112
112
113
113
114 class Beam:
114 class Beam:
115
115
116 def __init__(self):
116 def __init__(self):
117 self.codeList = []
117 self.codeList = []
118 self.azimuthList = []
118 self.azimuthList = []
119 self.zenithList = []
119 self.zenithList = []
120
120
121
121
122 class GenericData(object):
122 class GenericData(object):
123
123
124 flagNoData = True
124 flagNoData = True
125
125
126 def copy(self, inputObj=None):
126 def copy(self, inputObj=None):
127
127
128 if inputObj == None:
128 if inputObj == None:
129 return copy.deepcopy(self)
129 return copy.deepcopy(self)
130
130
131 for key in list(inputObj.__dict__.keys()):
131 for key in list(inputObj.__dict__.keys()):
132
132
133 attribute = inputObj.__dict__[key]
133 attribute = inputObj.__dict__[key]
134
134
135 # If this attribute is a tuple or list
135 # If this attribute is a tuple or list
136 if type(inputObj.__dict__[key]) in (tuple, list):
136 if type(inputObj.__dict__[key]) in (tuple, list):
137 self.__dict__[key] = attribute[:]
137 self.__dict__[key] = attribute[:]
138 continue
138 continue
139
139
140 # If this attribute is another object or instance
140 # If this attribute is another object or instance
141 if hasattr(attribute, '__dict__'):
141 if hasattr(attribute, '__dict__'):
142 self.__dict__[key] = attribute.copy()
142 self.__dict__[key] = attribute.copy()
143 continue
143 continue
144
144
145 self.__dict__[key] = inputObj.__dict__[key]
145 self.__dict__[key] = inputObj.__dict__[key]
146
146
147 def deepcopy(self):
147 def deepcopy(self):
148
148
149 return copy.deepcopy(self)
149 return copy.deepcopy(self)
150
150
151 def isEmpty(self):
151 def isEmpty(self):
152
152
153 return self.flagNoData
153 return self.flagNoData
154
154
155 def isReady(self):
155 def isReady(self):
156
156
157 return not self.flagNoData
157 return not self.flagNoData
158
158
159
159
160 class JROData(GenericData):
160 class JROData(GenericData):
161
161
162 systemHeaderObj = SystemHeader()
162 systemHeaderObj = SystemHeader()
163 radarControllerHeaderObj = RadarControllerHeader()
163 radarControllerHeaderObj = RadarControllerHeader()
164 type = None
164 type = None
165 datatype = None # dtype but in string
165 datatype = None # dtype but in string
166 nProfiles = None
166 nProfiles = None
167 heightList = None
167 heightList = None
168 channelList = None
168 channelList = None
169 flagDiscontinuousBlock = False
169 flagDiscontinuousBlock = False
170 useLocalTime = False
170 useLocalTime = False
171 utctime = None
171 utctime = None
172 timeZone = None
172 timeZone = None
173 dstFlag = None
173 dstFlag = None
174 errorCount = None
174 errorCount = None
175 blocksize = None
175 blocksize = None
176 flagDecodeData = False # asumo q la data no esta decodificada
176 flagDecodeData = False # asumo q la data no esta decodificada
177 flagDeflipData = False # asumo q la data no esta sin flip
177 flagDeflipData = False # asumo q la data no esta sin flip
178 flagShiftFFT = False
178 flagShiftFFT = False
179 nCohInt = None
179 nCohInt = None
180 windowOfFilter = 1
180 windowOfFilter = 1
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 metadata_list = ['heightList', 'timeZone', 'type']
194 metadata_list = ['heightList', 'timeZone', 'type']
195
195
196 ippFactor = 1 #Added to correct the freq and vel range for AMISR data
196 ippFactor = 1 #Added to correct the freq and vel range for AMISR data
197 useInputBuffer = False
197 useInputBuffer = False
198 buffer_empty = True
198 buffer_empty = True
199 codeList = []
199 codeList = []
200 azimuthList = []
200 azimuthList = []
201 elevationList = []
201 elevationList = []
202 last_noise = None
202 last_noise = None
203 __ipp = None
203 __ipp = None
204 __ippSeconds = None
204 __ippSeconds = None
205 sampled_heightsFFT = None
205 sampled_heightsFFT = None
206 pulseLength_TxA = None
206 pulseLength_TxA = None
207 deltaHeight = None
207 deltaHeight = None
208 __code = None
208 __code = None
209 __nCode = None
209 __nCode = None
210 __nBaud = None
210 __nBaud = None
211 unitsDescription = "The units of the parameters are according to the International System of units (Seconds, Meter, Hertz, ...), except \
211 unitsDescription = "The units of the parameters are according to the International System of units (Seconds, Meter, Hertz, ...), except \
212 the parameters related to distances such as heightList, or heightResolution wich are in Km"
212 the parameters related to distances such as heightList, or heightResolution wich are in Km"
213
213
214
214
215
215
216 def __str__(self):
216 def __str__(self):
217
217
218 return '{} - {}'.format(self.type, self.datatime())
218 return '{} - {}'.format(self.type, self.datatime())
219
219
220 def getNoise(self):
220 def getNoise(self):
221
221
222 raise NotImplementedError
222 raise NotImplementedError
223
223
224 @property
224 @property
225 def nChannels(self):
225 def nChannels(self):
226
226
227 return len(self.channelList)
227 return len(self.channelList)
228
228
229 @property
229 @property
230 def channelIndexList(self):
230 def channelIndexList(self):
231
231
232 return list(range(self.nChannels))
232 return list(range(self.nChannels))
233
233
234 @property
234 @property
235 def nHeights(self):
235 def nHeights(self):
236
236
237 return len(self.heightList)
237 return len(self.heightList)
238
238
239 def getDeltaH(self):
239 def getDeltaH(self):
240
240
241 return self.heightList[1] - self.heightList[0]
241 return self.heightList[1] - self.heightList[0]
242
242
243 @property
243 @property
244 def ltctime(self):
244 def ltctime(self):
245
245
246 if self.useLocalTime:
246 if self.useLocalTime:
247 return self.utctime - self.timeZone * 60
247 return self.utctime - self.timeZone * 60
248
248
249 return self.utctime
249 return self.utctime
250
250
251 @property
251 @property
252 def datatime(self):
252 def datatime(self):
253
253
254 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
254 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
255 return datatimeValue
255 return datatimeValue
256
256
257 def getTimeRange(self):
257 def getTimeRange(self):
258
258
259 datatime = []
259 datatime = []
260
260
261 datatime.append(self.ltctime)
261 datatime.append(self.ltctime)
262 datatime.append(self.ltctime + self.timeInterval + 1)
262 datatime.append(self.ltctime + self.timeInterval + 1)
263
263
264 datatime = numpy.array(datatime)
264 datatime = numpy.array(datatime)
265
265
266 return datatime
266 return datatime
267
267
268 def getFmaxTimeResponse(self):
268 def getFmaxTimeResponse(self):
269
269
270 period = (10**-6) * self.getDeltaH() / (0.15)
270 period = (10**-6) * self.getDeltaH() / (0.15)
271
271
272 PRF = 1. / (period * self.nCohInt)
272 PRF = 1. / (period * self.nCohInt)
273
273
274 fmax = PRF
274 fmax = PRF
275
275
276 return fmax
276 return fmax
277
277
278 def getFmax(self):
278 def getFmax(self):
279 PRF = 1. / (self.__ippSeconds * self.nCohInt)
279 PRF = 1. / (self.__ippSeconds * self.nCohInt)
280
280
281 fmax = PRF
281 fmax = PRF
282 return fmax
282 return fmax
283
283
284 def getVmax(self):
284 def getVmax(self):
285
285
286 _lambda = self.C / self.frequency
286 _lambda = self.C / self.frequency
287
287
288 vmax = self.getFmax() * _lambda / 2
288 vmax = self.getFmax() * _lambda / 2
289
289
290 return vmax
290 return vmax
291
291
292 ## Radar Controller Header must be immutable
292 ## Radar Controller Header must be immutable
293 @property
293 @property
294 def ippSeconds(self):
294 def ippSeconds(self):
295 '''
295 '''
296 '''
296 '''
297 #return self.radarControllerHeaderObj.ippSeconds
297 #return self.radarControllerHeaderObj.ippSeconds
298 return self.__ippSeconds
298 return self.__ippSeconds
299
299
300 @ippSeconds.setter
300 @ippSeconds.setter
301 def ippSeconds(self, ippSeconds):
301 def ippSeconds(self, ippSeconds):
302 '''
302 '''
303 '''
303 '''
304 #self.radarControllerHeaderObj.ippSeconds = ippSeconds
304 #self.radarControllerHeaderObj.ippSeconds = ippSeconds
305 self.__ippSeconds = ippSeconds
305 self.__ippSeconds = ippSeconds
306 self.__ipp = ippSeconds*SPEED_OF_LIGHT/2000.0
306 self.__ipp = ippSeconds*SPEED_OF_LIGHT/2000.0
307
307
308 @property
308 @property
309 def code(self):
309 def code(self):
310 '''
310 '''
311 '''
311 '''
312 # return self.radarControllerHeaderObj.code
312 # return self.radarControllerHeaderObj.code
313 return self.__code
313 return self.__code
314
314
315 @code.setter
315 @code.setter
316 def code(self, code):
316 def code(self, code):
317 '''
317 '''
318 '''
318 '''
319 # self.radarControllerHeaderObj.code = code
319 # self.radarControllerHeaderObj.code = code
320 self.__code = code
320 self.__code = code
321
321
322 @property
322 @property
323 def nCode(self):
323 def nCode(self):
324 '''
324 '''
325 '''
325 '''
326 # return self.radarControllerHeaderObj.nCode
326 # return self.radarControllerHeaderObj.nCode
327 return self.__nCode
327 return self.__nCode
328
328
329 @nCode.setter
329 @nCode.setter
330 def nCode(self, ncode):
330 def nCode(self, ncode):
331 '''
331 '''
332 '''
332 '''
333 # self.radarControllerHeaderObj.nCode = ncode
333 # self.radarControllerHeaderObj.nCode = ncode
334 self.__nCode = ncode
334 self.__nCode = ncode
335
335
336 @property
336 @property
337 def nBaud(self):
337 def nBaud(self):
338 '''
338 '''
339 '''
339 '''
340 # return self.radarControllerHeaderObj.nBaud
340 # return self.radarControllerHeaderObj.nBaud
341 return self.__nBaud
341 return self.__nBaud
342
342
343 @nBaud.setter
343 @nBaud.setter
344 def nBaud(self, nbaud):
344 def nBaud(self, nbaud):
345 '''
345 '''
346 '''
346 '''
347 # self.radarControllerHeaderObj.nBaud = nbaud
347 # self.radarControllerHeaderObj.nBaud = nbaud
348 self.__nBaud = nbaud
348 self.__nBaud = nbaud
349
349
350 @property
350 @property
351 def ipp(self):
351 def ipp(self):
352 '''
352 '''
353 '''
353 '''
354 # return self.radarControllerHeaderObj.ipp
354 # return self.radarControllerHeaderObj.ipp
355 return self.__ipp
355 return self.__ipp
356
356
357 @ipp.setter
357 @ipp.setter
358 def ipp(self, ipp):
358 def ipp(self, ipp):
359 '''
359 '''
360 '''
360 '''
361 # self.radarControllerHeaderObj.ipp = ipp
361 # self.radarControllerHeaderObj.ipp = ipp
362 self.__ipp = ipp
362 self.__ipp = ipp
363
363
364 @property
364 @property
365 def metadata(self):
365 def metadata(self):
366 '''
366 '''
367 '''
367 '''
368
368
369 return {attr: getattr(self, attr) for attr in self.metadata_list}
369 return {attr: getattr(self, attr) for attr in self.metadata_list}
370
370
371
371
372 class Voltage(JROData):
372 class Voltage(JROData):
373
373
374 dataPP_POW = None
374 dataPP_POW = None
375 dataPP_DOP = None
375 dataPP_DOP = None
376 dataPP_WIDTH = None
376 dataPP_WIDTH = None
377 dataPP_SNR = None
377 dataPP_SNR = None
378
378
379 # To use oper
379 # To use oper
380 flagProfilesByRange = False
380 flagProfilesByRange = False
381 nProfilesByRange = None
381 nProfilesByRange = None
382 max_nIncohInt = 1
382 max_nIncohInt = 1
383
383
384 def __init__(self):
384 def __init__(self):
385 '''
385 '''
386 Constructor
386 Constructor
387 '''
387 '''
388
388
389 self.useLocalTime = True
389 self.useLocalTime = True
390 self.radarControllerHeaderObj = RadarControllerHeader()
390 self.radarControllerHeaderObj = RadarControllerHeader()
391 self.systemHeaderObj = SystemHeader()
391 self.systemHeaderObj = SystemHeader()
392 self.processingHeaderObj = ProcessingHeader()
392 self.processingHeaderObj = ProcessingHeader()
393 self.type = "Voltage"
393 self.type = "Voltage"
394 self.data = None
394 self.data = None
395 self.nProfiles = None
395 self.nProfiles = None
396 self.heightList = None
396 self.heightList = None
397 self.channelList = None
397 self.channelList = None
398 self.flagNoData = True
398 self.flagNoData = True
399 self.flagDiscontinuousBlock = False
399 self.flagDiscontinuousBlock = False
400 self.utctime = None
400 self.utctime = None
401 self.timeZone = 0
401 self.timeZone = 0
402 self.dstFlag = None
402 self.dstFlag = None
403 self.errorCount = None
403 self.errorCount = None
404 self.nCohInt = None
404 self.nCohInt = None
405 self.blocksize = None
405 self.blocksize = None
406 self.flagCohInt = False
406 self.flagCohInt = False
407 self.flagDecodeData = False # asumo q la data no esta decodificada
407 self.flagDecodeData = False # asumo q la data no esta decodificada
408 self.flagDeflipData = False # asumo q la data no esta sin flip
408 self.flagDeflipData = False # asumo q la data no esta sin flip
409 self.flagShiftFFT = False
409 self.flagShiftFFT = False
410 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
410 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
411 self.profileIndex = 0
411 self.profileIndex = 0
412 self.ippFactor=1
412 self.ippFactor=1
413 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
413 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
414 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
414 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
415
415
416 def getNoisebyHildebrand(self, channel=None, ymin_index=None, ymax_index=None):
416 def getNoisebyHildebrand(self, channel=None, ymin_index=None, ymax_index=None):
417 """
417 """
418 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
418 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
419
419
420 Return:
420 Return:
421 noiselevel
421 noiselevel
422 """
422 """
423
423
424 if channel != None:
424 if channel != None:
425 data = self.data[channel,ymin_index:ymax_index]
425 data = self.data[channel,ymin_index:ymax_index]
426 nChannels = 1
426 nChannels = 1
427 else:
427 else:
428 data = self.data[:,ymin_index:ymax_index]
428 data = self.data[:,ymin_index:ymax_index]
429 nChannels = self.nChannels
429 nChannels = self.nChannels
430
430
431 noise = numpy.zeros(nChannels)
431 noise = numpy.zeros(nChannels)
432 power = data * numpy.conjugate(data)
432 power = data * numpy.conjugate(data)
433
433
434 for thisChannel in range(nChannels):
434 for thisChannel in range(nChannels):
435 if nChannels == 1:
435 if nChannels == 1:
436 daux = power[:].real
436 daux = power[:].real
437 else:
437 else:
438 daux = power[thisChannel, :].real
438 daux = power[thisChannel, :].real
439 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
439 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
440
440
441 return noise
441 return noise
442
442
443 def getNoise(self, type=1, channel=None,ymin_index=None, ymax_index=None):
443 def getNoise(self, type=1, channel=None,ymin_index=None, ymax_index=None):
444
444
445 if type == 1:
445 if type == 1:
446 noise = self.getNoisebyHildebrand(channel,ymin_index, ymax_index)
446 noise = self.getNoisebyHildebrand(channel,ymin_index, ymax_index)
447
447
448 return noise
448 return noise
449
449
450 def getPower(self, channel=None):
450 def getPower(self, channel=None):
451
451
452 if channel != None:
452 if channel != None:
453 data = self.data[channel]
453 data = self.data[channel]
454 else:
454 else:
455 data = self.data
455 data = self.data
456
456
457 power = data * numpy.conjugate(data)
457 power = data * numpy.conjugate(data)
458 powerdB = 10 * numpy.log10(power.real)
458 powerdB = 10 * numpy.log10(power.real)
459 powerdB = numpy.squeeze(powerdB)
459 powerdB = numpy.squeeze(powerdB)
460
460
461 return powerdB
461 return powerdB
462 @property
463 def data_pow(self):
464 return self.getPower()
462
465
463 @property
466 @property
464 def timeInterval(self):
467 def timeInterval(self):
465
468
466 return self.ippSeconds * self.nCohInt
469 return self.ippSeconds * self.nCohInt
467
470
468 noise = property(getNoise, "I'm the 'nHeights' property.")
471 noise = property(getNoise, "I'm the 'nHeights' property.")
469
472
470
473
471 class Spectra(JROData):
474 class Spectra(JROData):
472
475
473 data_outlier = None
476 data_outlier = None
474 flagProfilesByRange = False
477 flagProfilesByRange = False
475 nProfilesByRange = None
478 nProfilesByRange = None
476
479
477 def __init__(self):
480 def __init__(self):
478 '''
481 '''
479 Constructor
482 Constructor
480 '''
483 '''
481
484
482 self.data_dc = None
485 self.data_dc = None
483 self.data_spc = None
486 self.data_spc = None
484 self.data_cspc = None
487 self.data_cspc = None
485 self.useLocalTime = True
488 self.useLocalTime = True
486 self.radarControllerHeaderObj = RadarControllerHeader()
489 self.radarControllerHeaderObj = RadarControllerHeader()
487 self.systemHeaderObj = SystemHeader()
490 self.systemHeaderObj = SystemHeader()
488 self.processingHeaderObj = ProcessingHeader()
491 self.processingHeaderObj = ProcessingHeader()
489 self.type = "Spectra"
492 self.type = "Spectra"
490 self.timeZone = 0
493 self.timeZone = 0
491 self.nProfiles = None
494 self.nProfiles = None
492 self.heightList = None
495 self.heightList = None
493 self.channelList = None
496 self.channelList = None
494 self.pairsList = None
497 self.pairsList = None
495 self.flagNoData = True
498 self.flagNoData = True
496 self.flagDiscontinuousBlock = False
499 self.flagDiscontinuousBlock = False
497 self.utctime = None
500 self.utctime = None
498 self.nCohInt = None
501 self.nCohInt = None
499 self.nIncohInt = None
502 self.nIncohInt = None
500 self.blocksize = None
503 self.blocksize = None
501 self.nFFTPoints = None
504 self.nFFTPoints = None
502 self.wavelength = None
505 self.wavelength = None
503 self.flagDecodeData = False # asumo q la data no esta decodificada
506 self.flagDecodeData = False # asumo q la data no esta decodificada
504 self.flagDeflipData = False # asumo q la data no esta sin flip
507 self.flagDeflipData = False # asumo q la data no esta sin flip
505 self.flagShiftFFT = False
508 self.flagShiftFFT = False
506 self.ippFactor = 1
509 self.ippFactor = 1
507 self.beacon_heiIndexList = []
510 self.beacon_heiIndexList = []
508 self.noise_estimation = None
511 self.noise_estimation = None
509 self.codeList = []
512 self.codeList = []
510 self.azimuthList = []
513 self.azimuthList = []
511 self.elevationList = []
514 self.elevationList = []
512 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
515 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
513 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
516 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
514
517
515 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
518 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
516 """
519 """
517 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
520 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
518
521
519 Return:
522 Return:
520 noiselevel
523 noiselevel
521 """
524 """
522
525
523 noise = numpy.zeros(self.nChannels)
526 noise = numpy.zeros(self.nChannels)
524
527
525 for channel in range(self.nChannels):
528 for channel in range(self.nChannels):
526 daux = self.data_spc[channel,
529 daux = self.data_spc[channel,
527 xmin_index:xmax_index, ymin_index:ymax_index]
530 xmin_index:xmax_index, ymin_index:ymax_index]
528 # noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
531 # noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
529 noise[channel] = hildebrand_sekhon(daux, self.max_nIncohInt[channel])
532 noise[channel] = hildebrand_sekhon(daux, self.max_nIncohInt[channel])
530
533
531 return noise
534 return noise
532
535
533 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
536 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
534
537
535 if self.noise_estimation is not None:
538 if self.noise_estimation is not None:
536 # this was estimated by getNoise Operation defined in jroproc_spectra.py
539 # this was estimated by getNoise Operation defined in jroproc_spectra.py
537 return self.noise_estimation
540 return self.noise_estimation
538 else:
541 else:
539 noise = self.getNoisebyHildebrand(
542 noise = self.getNoisebyHildebrand(
540 xmin_index, xmax_index, ymin_index, ymax_index)
543 xmin_index, xmax_index, ymin_index, ymax_index)
541 return noise
544 return noise
542
545
543 def getFreqRangeTimeResponse(self, extrapoints=0):
546 def getFreqRangeTimeResponse(self, extrapoints=0):
544
547
545 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
548 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
546 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
549 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
547
550
548 return freqrange
551 return freqrange
549
552
550 def getAcfRange(self, extrapoints=0):
553 def getAcfRange(self, extrapoints=0):
551
554
552 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
555 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
553 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
556 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
554
557
555 return freqrange
558 return freqrange
556
559
557 def getFreqRange(self, extrapoints=0):
560 def getFreqRange(self, extrapoints=0):
558
561
559 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
562 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
560 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
563 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
561
564
562 return freqrange
565 return freqrange
563
566
564 def getVelRange(self, extrapoints=0):
567 def getVelRange(self, extrapoints=0):
565
568
566 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
569 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
567 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
570 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
568
571
569 if self.nmodes:
572 if self.nmodes:
570 return velrange/self.nmodes
573 return velrange/self.nmodes
571 else:
574 else:
572 return velrange
575 return velrange
573
576
574 @property
577 @property
575 def nPairs(self):
578 def nPairs(self):
576
579
577 return len(self.pairsList)
580 return len(self.pairsList)
578
581
579 @property
582 @property
580 def pairsIndexList(self):
583 def pairsIndexList(self):
581
584
582 return list(range(self.nPairs))
585 return list(range(self.nPairs))
583
586
584 @property
587 @property
585 def normFactor(self):
588 def normFactor(self):
586
589
587 pwcode = 1
590 pwcode = 1
588 if self.flagDecodeData:
591 if self.flagDecodeData:
589 try:
592 try:
590 pwcode = numpy.sum(self.code[0]**2)
593 pwcode = numpy.sum(self.code[0]**2)
591 except Exception as e:
594 except Exception as e:
592 log.warning("Failed pwcode read, setting to 1")
595 log.warning("Failed pwcode read, setting to 1")
593 pwcode = 1
596 pwcode = 1
594 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
597 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
595 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
598 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
596 if self.flagProfilesByRange:
599 if self.flagProfilesByRange:
597 normFactor *= (self.nProfilesByRange/self.nProfilesByRange.max())
600 normFactor *= (self.nProfilesByRange/self.nProfilesByRange.max())
598 return normFactor
601 return normFactor
599
602
600 @property
603 @property
601 def flag_cspc(self):
604 def flag_cspc(self):
602
605
603 if self.data_cspc is None:
606 if self.data_cspc is None:
604 return True
607 return True
605
608
606 return False
609 return False
607
610
608 @property
611 @property
609 def flag_dc(self):
612 def flag_dc(self):
610
613
611 if self.data_dc is None:
614 if self.data_dc is None:
612 return True
615 return True
613
616
614 return False
617 return False
615
618
616 @property
619 @property
617 def timeInterval(self):
620 def timeInterval(self):
618
621
619 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
622 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
620 if self.nmodes:
623 if self.nmodes:
621 return self.nmodes*timeInterval
624 return self.nmodes*timeInterval
622 else:
625 else:
623 return timeInterval
626 return timeInterval
624
627
625 def getPower(self):
628 def getPower(self):
626
629
627 factor = self.normFactor
630 factor = self.normFactor
628 power = numpy.zeros( (self.nChannels,self.nHeights) )
631 power = numpy.zeros( (self.nChannels,self.nHeights) )
629 for ch in range(self.nChannels):
632 for ch in range(self.nChannels):
630 z = None
633 z = None
631 if hasattr(factor,'shape'):
634 if hasattr(factor,'shape'):
632 if factor.ndim > 1:
635 if factor.ndim > 1:
633 z = self.data_spc[ch]/factor[ch]
636 z = self.data_spc[ch]/factor[ch]
634 else:
637 else:
635 z = self.data_spc[ch]/factor
638 z = self.data_spc[ch]/factor
636 else:
639 else:
637 z = self.data_spc[ch]/factor
640 z = self.data_spc[ch]/factor
638 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
641 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
639 avg = numpy.average(z, axis=0)
642 avg = numpy.average(z, axis=0)
640 power[ch] = 10 * numpy.log10(avg)
643 power[ch] = 10 * numpy.log10(avg)
641 return power
644 return power
642
645
643 @property
646 @property
644 def max_nIncohInt(self):
647 def max_nIncohInt(self):
645
648
646 ints = numpy.zeros(self.nChannels)
649 ints = numpy.zeros(self.nChannels)
647 for ch in range(self.nChannels):
650 for ch in range(self.nChannels):
648 if hasattr(self.nIncohInt,'shape'):
651 if hasattr(self.nIncohInt,'shape'):
649 if self.nIncohInt.ndim > 1:
652 if self.nIncohInt.ndim > 1:
650 ints[ch,] = self.nIncohInt[ch].max()
653 ints[ch,] = self.nIncohInt[ch].max()
651 else:
654 else:
652 ints[ch,] = self.nIncohInt
655 ints[ch,] = self.nIncohInt
653 self.nIncohInt = int(self.nIncohInt)
656 self.nIncohInt = int(self.nIncohInt)
654 else:
657 else:
655 ints[ch,] = self.nIncohInt
658 ints[ch,] = self.nIncohInt
656
659
657 return ints
660 return ints
658
661
659 def getCoherence(self, pairsList=None, phase=False):
662 def getCoherence(self, pairsList=None, phase=False):
660
663
661 z = []
664 z = []
662 if pairsList is None:
665 if pairsList is None:
663 pairsIndexList = self.pairsIndexList
666 pairsIndexList = self.pairsIndexList
664 else:
667 else:
665 pairsIndexList = []
668 pairsIndexList = []
666 for pair in pairsList:
669 for pair in pairsList:
667 if pair not in self.pairsList:
670 if pair not in self.pairsList:
668 raise ValueError("Pair %s is not in dataOut.pairsList" % (
671 raise ValueError("Pair %s is not in dataOut.pairsList" % (
669 pair))
672 pair))
670 pairsIndexList.append(self.pairsList.index(pair))
673 pairsIndexList.append(self.pairsList.index(pair))
671 for i in range(len(pairsIndexList)):
674 for i in range(len(pairsIndexList)):
672 pair = self.pairsList[pairsIndexList[i]]
675 pair = self.pairsList[pairsIndexList[i]]
673 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
676 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
674 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
677 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
675 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
678 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
676 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
679 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
677 if phase:
680 if phase:
678 data = numpy.arctan2(avgcoherenceComplex.imag,
681 data = numpy.arctan2(avgcoherenceComplex.imag,
679 avgcoherenceComplex.real) * 180 / numpy.pi
682 avgcoherenceComplex.real) * 180 / numpy.pi
680 else:
683 else:
681 data = numpy.abs(avgcoherenceComplex)
684 data = numpy.abs(avgcoherenceComplex)
682
685
683 z.append(data)
686 z.append(data)
684
687
685 return numpy.array(z)
688 return numpy.array(z)
686
689
687 def setValue(self, value):
690 def setValue(self, value):
688
691
689 print("This property should not be initialized", value)
692 print("This property should not be initialized", value)
690
693
691 return
694 return
692
695
693 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
696 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
694
697
695
698
696 class SpectraHeis(Spectra):
699 class SpectraHeis(Spectra):
697
700
698 def __init__(self):
701 def __init__(self):
699
702
700 self.radarControllerHeaderObj = RadarControllerHeader()
703 self.radarControllerHeaderObj = RadarControllerHeader()
701 self.systemHeaderObj = SystemHeader()
704 self.systemHeaderObj = SystemHeader()
702 self.type = "SpectraHeis"
705 self.type = "SpectraHeis"
703 self.nProfiles = None
706 self.nProfiles = None
704 self.heightList = None
707 self.heightList = None
705 self.channelList = None
708 self.channelList = None
706 self.flagNoData = True
709 self.flagNoData = True
707 self.flagDiscontinuousBlock = False
710 self.flagDiscontinuousBlock = False
708 self.utctime = None
711 self.utctime = None
709 self.blocksize = None
712 self.blocksize = None
710 self.profileIndex = 0
713 self.profileIndex = 0
711 self.nCohInt = 1
714 self.nCohInt = 1
712 self.nIncohInt = 1
715 self.nIncohInt = 1
713
716
714 @property
717 @property
715 def normFactor(self):
718 def normFactor(self):
716 pwcode = 1
719 pwcode = 1
717 if self.flagDecodeData:
720 if self.flagDecodeData:
718 pwcode = numpy.sum(self.code[0]**2)
721 pwcode = numpy.sum(self.code[0]**2)
719
722
720 normFactor = self.nIncohInt * self.nCohInt * pwcode
723 normFactor = self.nIncohInt * self.nCohInt * pwcode
721
724
722 return normFactor
725 return normFactor
723
726
724 @property
727 @property
725 def timeInterval(self):
728 def timeInterval(self):
726
729
727 return self.ippSeconds * self.nCohInt * self.nIncohInt
730 return self.ippSeconds * self.nCohInt * self.nIncohInt
728
731
729
732
730 class Fits(JROData):
733 class Fits(JROData):
731
734
732 def __init__(self):
735 def __init__(self):
733
736
734 self.type = "Fits"
737 self.type = "Fits"
735 self.nProfiles = None
738 self.nProfiles = None
736 self.heightList = None
739 self.heightList = None
737 self.channelList = None
740 self.channelList = None
738 self.flagNoData = True
741 self.flagNoData = True
739 self.utctime = None
742 self.utctime = None
740 self.nCohInt = 1
743 self.nCohInt = 1
741 self.nIncohInt = 1
744 self.nIncohInt = 1
742 self.useLocalTime = True
745 self.useLocalTime = True
743 self.profileIndex = 0
746 self.profileIndex = 0
744 self.timeZone = 0
747 self.timeZone = 0
745
748
746 def getTimeRange(self):
749 def getTimeRange(self):
747
750
748 datatime = []
751 datatime = []
749
752
750 datatime.append(self.ltctime)
753 datatime.append(self.ltctime)
751 datatime.append(self.ltctime + self.timeInterval)
754 datatime.append(self.ltctime + self.timeInterval)
752
755
753 datatime = numpy.array(datatime)
756 datatime = numpy.array(datatime)
754
757
755 return datatime
758 return datatime
756
759
757 def getChannelIndexList(self):
760 def getChannelIndexList(self):
758
761
759 return list(range(self.nChannels))
762 return list(range(self.nChannels))
760
763
761 def getNoise(self, type=1):
764 def getNoise(self, type=1):
762
765
763
766
764 if type == 1:
767 if type == 1:
765 noise = self.getNoisebyHildebrand()
768 noise = self.getNoisebyHildebrand()
766
769
767 if type == 2:
770 if type == 2:
768 noise = self.getNoisebySort()
771 noise = self.getNoisebySort()
769
772
770 if type == 3:
773 if type == 3:
771 noise = self.getNoisebyWindow()
774 noise = self.getNoisebyWindow()
772
775
773 return noise
776 return noise
774
777
775 @property
778 @property
776 def timeInterval(self):
779 def timeInterval(self):
777
780
778 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
781 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
779
782
780 return timeInterval
783 return timeInterval
781
784
782 @property
785 @property
783 def ippSeconds(self):
786 def ippSeconds(self):
784 '''
787 '''
785 '''
788 '''
786 return self.ipp_sec
789 return self.ipp_sec
787
790
788 noise = property(getNoise, "I'm the 'nHeights' property.")
791 noise = property(getNoise, "I'm the 'nHeights' property.")
789
792
790
793
791 class Correlation(JROData):
794 class Correlation(JROData):
792
795
793 def __init__(self):
796 def __init__(self):
794 '''
797 '''
795 Constructor
798 Constructor
796 '''
799 '''
797 self.radarControllerHeaderObj = RadarControllerHeader()
800 self.radarControllerHeaderObj = RadarControllerHeader()
798 self.systemHeaderObj = SystemHeader()
801 self.systemHeaderObj = SystemHeader()
799 self.type = "Correlation"
802 self.type = "Correlation"
800 self.data = None
803 self.data = None
801 self.dtype = None
804 self.dtype = None
802 self.nProfiles = None
805 self.nProfiles = None
803 self.heightList = None
806 self.heightList = None
804 self.channelList = None
807 self.channelList = None
805 self.flagNoData = True
808 self.flagNoData = True
806 self.flagDiscontinuousBlock = False
809 self.flagDiscontinuousBlock = False
807 self.utctime = None
810 self.utctime = None
808 self.timeZone = 0
811 self.timeZone = 0
809 self.dstFlag = None
812 self.dstFlag = None
810 self.errorCount = None
813 self.errorCount = None
811 self.blocksize = None
814 self.blocksize = None
812 self.flagDecodeData = False # asumo q la data no esta decodificada
815 self.flagDecodeData = False # asumo q la data no esta decodificada
813 self.flagDeflipData = False # asumo q la data no esta sin flip
816 self.flagDeflipData = False # asumo q la data no esta sin flip
814 self.pairsList = None
817 self.pairsList = None
815 self.nPoints = None
818 self.nPoints = None
816
819
817 def getPairsList(self):
820 def getPairsList(self):
818
821
819 return self.pairsList
822 return self.pairsList
820
823
821 def getNoise(self, mode=2):
824 def getNoise(self, mode=2):
822
825
823 indR = numpy.where(self.lagR == 0)[0][0]
826 indR = numpy.where(self.lagR == 0)[0][0]
824 indT = numpy.where(self.lagT == 0)[0][0]
827 indT = numpy.where(self.lagT == 0)[0][0]
825
828
826 jspectra0 = self.data_corr[:, :, indR, :]
829 jspectra0 = self.data_corr[:, :, indR, :]
827 jspectra = copy.copy(jspectra0)
830 jspectra = copy.copy(jspectra0)
828
831
829 num_chan = jspectra.shape[0]
832 num_chan = jspectra.shape[0]
830 num_hei = jspectra.shape[2]
833 num_hei = jspectra.shape[2]
831
834
832 freq_dc = jspectra.shape[1] / 2
835 freq_dc = jspectra.shape[1] / 2
833 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
836 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
834
837
835 if ind_vel[0] < 0:
838 if ind_vel[0] < 0:
836 ind_vel[list(range(0, 1))] = ind_vel[list(
839 ind_vel[list(range(0, 1))] = ind_vel[list(
837 range(0, 1))] + self.num_prof
840 range(0, 1))] + self.num_prof
838
841
839 if mode == 1:
842 if mode == 1:
840 jspectra[:, freq_dc, :] = (
843 jspectra[:, freq_dc, :] = (
841 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
844 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
842
845
843 if mode == 2:
846 if mode == 2:
844
847
845 vel = numpy.array([-2, -1, 1, 2])
848 vel = numpy.array([-2, -1, 1, 2])
846 xx = numpy.zeros([4, 4])
849 xx = numpy.zeros([4, 4])
847
850
848 for fil in range(4):
851 for fil in range(4):
849 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
852 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
850
853
851 xx_inv = numpy.linalg.inv(xx)
854 xx_inv = numpy.linalg.inv(xx)
852 xx_aux = xx_inv[0, :]
855 xx_aux = xx_inv[0, :]
853
856
854 for ich in range(num_chan):
857 for ich in range(num_chan):
855 yy = jspectra[ich, ind_vel, :]
858 yy = jspectra[ich, ind_vel, :]
856 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
859 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
857
860
858 junkid = jspectra[ich, freq_dc, :] <= 0
861 junkid = jspectra[ich, freq_dc, :] <= 0
859 cjunkid = sum(junkid)
862 cjunkid = sum(junkid)
860
863
861 if cjunkid.any():
864 if cjunkid.any():
862 jspectra[ich, freq_dc, junkid.nonzero()] = (
865 jspectra[ich, freq_dc, junkid.nonzero()] = (
863 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
866 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
864
867
865 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
868 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
866
869
867 return noise
870 return noise
868
871
869 @property
872 @property
870 def timeInterval(self):
873 def timeInterval(self):
871
874
872 return self.ippSeconds * self.nCohInt * self.nProfiles
875 return self.ippSeconds * self.nCohInt * self.nProfiles
873
876
874 def splitFunctions(self):
877 def splitFunctions(self):
875
878
876 pairsList = self.pairsList
879 pairsList = self.pairsList
877 ccf_pairs = []
880 ccf_pairs = []
878 acf_pairs = []
881 acf_pairs = []
879 ccf_ind = []
882 ccf_ind = []
880 acf_ind = []
883 acf_ind = []
881 for l in range(len(pairsList)):
884 for l in range(len(pairsList)):
882 chan0 = pairsList[l][0]
885 chan0 = pairsList[l][0]
883 chan1 = pairsList[l][1]
886 chan1 = pairsList[l][1]
884
887
885 # Obteniendo pares de Autocorrelacion
888 # Obteniendo pares de Autocorrelacion
886 if chan0 == chan1:
889 if chan0 == chan1:
887 acf_pairs.append(chan0)
890 acf_pairs.append(chan0)
888 acf_ind.append(l)
891 acf_ind.append(l)
889 else:
892 else:
890 ccf_pairs.append(pairsList[l])
893 ccf_pairs.append(pairsList[l])
891 ccf_ind.append(l)
894 ccf_ind.append(l)
892
895
893 data_acf = self.data_cf[acf_ind]
896 data_acf = self.data_cf[acf_ind]
894 data_ccf = self.data_cf[ccf_ind]
897 data_ccf = self.data_cf[ccf_ind]
895
898
896 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
899 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
897
900
898 @property
901 @property
899 def normFactor(self):
902 def normFactor(self):
900 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
903 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
901 acf_pairs = numpy.array(acf_pairs)
904 acf_pairs = numpy.array(acf_pairs)
902 normFactor = numpy.zeros((self.nPairs, self.nHeights))
905 normFactor = numpy.zeros((self.nPairs, self.nHeights))
903
906
904 for p in range(self.nPairs):
907 for p in range(self.nPairs):
905 pair = self.pairsList[p]
908 pair = self.pairsList[p]
906
909
907 ch0 = pair[0]
910 ch0 = pair[0]
908 ch1 = pair[1]
911 ch1 = pair[1]
909
912
910 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
913 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
911 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
914 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
912 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
915 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
913
916
914 return normFactor
917 return normFactor
915
918
916
919
917 class Parameters(Spectra):
920 class Parameters(Spectra):
918
921
919 groupList = None # List of Pairs, Groups, etc
922 groupList = None # List of Pairs, Groups, etc
920 data_param = None # Parameters obtained
923 data_param = None # Parameters obtained
921 data_pre = None # Data Pre Parametrization
924 data_pre = None # Data Pre Parametrization
922 data_SNR = None # Signal to Noise Ratio
925 data_SNR = None # Signal to Noise Ratio
923 abscissaList = None # Abscissa, can be velocities, lags or time
926 abscissaList = None # Abscissa, can be velocities, lags or time
924 utctimeInit = None # Initial UTC time
927 utctimeInit = None # Initial UTC time
925 paramInterval = None # Time interval to calculate Parameters in seconds
928 paramInterval = None # Time interval to calculate Parameters in seconds
926 useLocalTime = True
929 useLocalTime = True
927 # Fitting
930 # Fitting
928 data_error = None # Error of the estimation
931 data_error = None # Error of the estimation
929 constants = None
932 constants = None
930 library = None
933 library = None
931 # Output signal
934 # Output signal
932 outputInterval = None # Time interval to calculate output signal in seconds
935 outputInterval = None # Time interval to calculate output signal in seconds
933 data_output = None # Out signal
936 data_output = None # Out signal
934 nAvg = None
937 nAvg = None
935 noise_estimation = None
938 noise_estimation = None
936 GauSPC = None # Fit gaussian SPC
939 GauSPC = None # Fit gaussian SPC
937
940
938 data_outlier = None
941 data_outlier = None
939 data_vdrift = None
942 data_vdrift = None
940 radarControllerHeaderTxt=None #header Controller like text
943 radarControllerHeaderTxt=None #header Controller like text
941 txPower = None
944 txPower = None
942 flagProfilesByRange = False
945 flagProfilesByRange = False
943 nProfilesByRange = None
946 nProfilesByRange = None
944
947
945
948
946 def __init__(self):
949 def __init__(self):
947 '''
950 '''
948 Constructor
951 Constructor
949 '''
952 '''
950 self.radarControllerHeaderObj = RadarControllerHeader()
953 self.radarControllerHeaderObj = RadarControllerHeader()
951 self.systemHeaderObj = SystemHeader()
954 self.systemHeaderObj = SystemHeader()
952 self.processingHeaderObj = ProcessingHeader()
955 self.processingHeaderObj = ProcessingHeader()
953 self.type = "Parameters"
956 self.type = "Parameters"
954 self.timeZone = 0
957 self.timeZone = 0
955
958
956 def getTimeRange1(self, interval):
959 def getTimeRange1(self, interval):
957
960
958 datatime = []
961 datatime = []
959
962
960 if self.useLocalTime:
963 if self.useLocalTime:
961 time1 = self.utctimeInit - self.timeZone * 60
964 time1 = self.utctimeInit - self.timeZone * 60
962 else:
965 else:
963 time1 = self.utctimeInit
966 time1 = self.utctimeInit
964
967
965 datatime.append(time1)
968 datatime.append(time1)
966 datatime.append(time1 + interval)
969 datatime.append(time1 + interval)
967 datatime = numpy.array(datatime)
970 datatime = numpy.array(datatime)
968
971
969 return datatime
972 return datatime
970
973
971 @property
974 @property
972 def timeInterval(self):
975 def timeInterval(self):
973
976
974 if hasattr(self, 'timeInterval1'):
977 if hasattr(self, 'timeInterval1'):
975 return self.timeInterval1
978 return self.timeInterval1
976 else:
979 else:
977 return self.paramInterval
980 return self.paramInterval
978
981
979 def setValue(self, value):
982 def setValue(self, value):
980
983
981 print("This property should not be initialized")
984 print("This property should not be initialized")
982
985
983 return
986 return
984
987
985 def getNoise(self):
988 def getNoise(self):
986
989
987 return self.spc_noise
990 return self.spc_noise
988
991
989 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
992 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
990
993
991
994
992 class PlotterData(object):
995 class PlotterData(object):
993 '''
996 '''
994 Object to hold data to be plotted
997 Object to hold data to be plotted
995 '''
998 '''
996
999
997 MAXNUMX = 200
1000 MAXNUMX = 200
998 MAXNUMY = 200
1001 MAXNUMY = 200
999
1002
1000 def __init__(self, code, exp_code, localtime=True):
1003 def __init__(self, code, exp_code, localtime=True):
1001
1004
1002 self.key = code
1005 self.key = code
1003 self.exp_code = exp_code
1006 self.exp_code = exp_code
1004 self.ready = False
1007 self.ready = False
1005 self.flagNoData = False
1008 self.flagNoData = False
1006 self.localtime = localtime
1009 self.localtime = localtime
1007 self.data = {}
1010 self.data = {}
1008 self.meta = {}
1011 self.meta = {}
1009 self.__heights = []
1012 self.__heights = []
1010
1013
1011 def __str__(self):
1014 def __str__(self):
1012 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1015 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1013 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
1016 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
1014
1017
1015 def __len__(self):
1018 def __len__(self):
1016 return len(self.data)
1019 return len(self.data)
1017
1020
1018 def __getitem__(self, key):
1021 def __getitem__(self, key):
1019 if isinstance(key, int):
1022 if isinstance(key, int):
1020 return self.data[self.times[key]]
1023 return self.data[self.times[key]]
1021 elif isinstance(key, str):
1024 elif isinstance(key, str):
1022 ret = numpy.array([self.data[x][key] for x in self.times])
1025 ret = numpy.array([self.data[x][key] for x in self.times])
1023 if ret.ndim > 1:
1026 if ret.ndim > 1:
1024 ret = numpy.swapaxes(ret, 0, 1)
1027 ret = numpy.swapaxes(ret, 0, 1)
1025 return ret
1028 return ret
1026
1029
1027 def __contains__(self, key):
1030 def __contains__(self, key):
1028 return key in self.data[self.min_time]
1031 return key in self.data[self.min_time]
1029
1032
1030 def setup(self):
1033 def setup(self):
1031 '''
1034 '''
1032 Configure object
1035 Configure object
1033 '''
1036 '''
1034 self.type = ''
1037 self.type = ''
1035 self.ready = False
1038 self.ready = False
1036 del self.data
1039 del self.data
1037 self.data = {}
1040 self.data = {}
1038 self.__heights = []
1041 self.__heights = []
1039 self.__all_heights = set()
1042 self.__all_heights = set()
1040
1043
1041 def shape(self, key):
1044 def shape(self, key):
1042 '''
1045 '''
1043 Get the shape of the one-element data for the given key
1046 Get the shape of the one-element data for the given key
1044 '''
1047 '''
1045
1048
1046 if len(self.data[self.min_time][key]):
1049 if len(self.data[self.min_time][key]):
1047 return self.data[self.min_time][key].shape
1050 return self.data[self.min_time][key].shape
1048 return (0,)
1051 return (0,)
1049
1052
1050 def update(self, data, tm, meta={}):
1053 def update(self, data, tm, meta={}):
1051 '''
1054 '''
1052 Update data object with new dataOut
1055 Update data object with new dataOut
1053 '''
1056 '''
1054
1057
1055 self.data[tm] = data
1058 self.data[tm] = data
1056
1059
1057 for key, value in meta.items():
1060 for key, value in meta.items():
1058 setattr(self, key, value)
1061 setattr(self, key, value)
1059
1062
1060 def normalize_heights(self):
1063 def normalize_heights(self):
1061 '''
1064 '''
1062 Ensure same-dimension of the data for different heighList
1065 Ensure same-dimension of the data for different heighList
1063 '''
1066 '''
1064
1067
1065 H = numpy.array(list(self.__all_heights))
1068 H = numpy.array(list(self.__all_heights))
1066 H.sort()
1069 H.sort()
1067 for key in self.data:
1070 for key in self.data:
1068 shape = self.shape(key)[:-1] + H.shape
1071 shape = self.shape(key)[:-1] + H.shape
1069 for tm, obj in list(self.data[key].items()):
1072 for tm, obj in list(self.data[key].items()):
1070 h = self.__heights[self.times.tolist().index(tm)]
1073 h = self.__heights[self.times.tolist().index(tm)]
1071 if H.size == h.size:
1074 if H.size == h.size:
1072 continue
1075 continue
1073 index = numpy.where(numpy.in1d(H, h))[0]
1076 index = numpy.where(numpy.in1d(H, h))[0]
1074 dummy = numpy.zeros(shape) + numpy.nan
1077 dummy = numpy.zeros(shape) + numpy.nan
1075 if len(shape) == 2:
1078 if len(shape) == 2:
1076 dummy[:, index] = obj
1079 dummy[:, index] = obj
1077 else:
1080 else:
1078 dummy[index] = obj
1081 dummy[index] = obj
1079 self.data[key][tm] = dummy
1082 self.data[key][tm] = dummy
1080
1083
1081 self.__heights = [H for tm in self.times]
1084 self.__heights = [H for tm in self.times]
1082
1085
1083 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1086 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1084 '''
1087 '''
1085 Convert data to json
1088 Convert data to json
1086 '''
1089 '''
1087
1090
1088 meta = {}
1091 meta = {}
1089 meta['xrange'] = []
1092 meta['xrange'] = []
1090 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1093 dy = int(len(self.yrange)/self.MAXNUMY) + 1
1091 tmp = self.data[tm][self.key]
1094 tmp = self.data[tm][self.key]
1092 shape = tmp.shape
1095 shape = tmp.shape
1093 if len(shape) == 2:
1096 if len(shape) == 2:
1094 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1097 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1095 elif len(shape) == 3:
1098 elif len(shape) == 3:
1096 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1099 dx = int(self.data[tm][self.key].shape[1]/self.MAXNUMX) + 1
1097 data = self.roundFloats(
1100 data = self.roundFloats(
1098 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1101 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1099 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1102 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1100 else:
1103 else:
1101 data = self.roundFloats(self.data[tm][self.key].tolist())
1104 data = self.roundFloats(self.data[tm][self.key].tolist())
1102
1105
1103 ret = {
1106 ret = {
1104 'plot': plot_name,
1107 'plot': plot_name,
1105 'code': self.exp_code,
1108 'code': self.exp_code,
1106 'time': float(tm),
1109 'time': float(tm),
1107 'data': data,
1110 'data': data,
1108 }
1111 }
1109 meta['type'] = plot_type
1112 meta['type'] = plot_type
1110 meta['interval'] = float(self.interval)
1113 meta['interval'] = float(self.interval)
1111 meta['localtime'] = self.localtime
1114 meta['localtime'] = self.localtime
1112 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1115 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1113 meta.update(self.meta)
1116 meta.update(self.meta)
1114 ret['metadata'] = meta
1117 ret['metadata'] = meta
1115 return json.dumps(ret)
1118 return json.dumps(ret)
1116
1119
1117 @property
1120 @property
1118 def times(self):
1121 def times(self):
1119 '''
1122 '''
1120 Return the list of times of the current data
1123 Return the list of times of the current data
1121 '''
1124 '''
1122
1125
1123 ret = [t for t in self.data]
1126 ret = [t for t in self.data]
1124 ret.sort()
1127 ret.sort()
1125 return numpy.array(ret)
1128 return numpy.array(ret)
1126
1129
1127 @property
1130 @property
1128 def min_time(self):
1131 def min_time(self):
1129 '''
1132 '''
1130 Return the minimun time value
1133 Return the minimun time value
1131 '''
1134 '''
1132
1135
1133 return self.times[0]
1136 return self.times[0]
1134
1137
1135 @property
1138 @property
1136 def max_time(self):
1139 def max_time(self):
1137 '''
1140 '''
1138 Return the maximun time value
1141 Return the maximun time value
1139 '''
1142 '''
1140
1143
1141 return self.times[-1]
1144 return self.times[-1]
1142
1145
1143 # @property
1146 # @property
1144 # def heights(self):
1147 # def heights(self):
1145 # '''
1148 # '''
1146 # Return the list of heights of the current data
1149 # Return the list of heights of the current data
1147 # '''
1150 # '''
1148
1151
1149 # return numpy.array(self.__heights[-1])
1152 # return numpy.array(self.__heights[-1])
1150
1153
1151 @staticmethod
1154 @staticmethod
1152 def roundFloats(obj):
1155 def roundFloats(obj):
1153 if isinstance(obj, list):
1156 if isinstance(obj, list):
1154 return list(map(PlotterData.roundFloats, obj))
1157 return list(map(PlotterData.roundFloats, obj))
1155 elif isinstance(obj, float):
1158 elif isinstance(obj, float):
1156 return round(obj, 2)
1159 return round(obj, 2)
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
NO CONTENT: modified file
Show file before | Show file after | Hide comments
@@ -1,1928 +1,1935
1 # Copyright (c) 2012-2021 Jicamarca Radio Observatory
1 # Copyright (c) 2012-2021 Jicamarca Radio Observatory
2 # All rights reserved.
2 # All rights reserved.
3 #
3 #
4 # Distributed under the terms of the BSD 3-clause license.
4 # Distributed under the terms of the BSD 3-clause license.
5 """Classes to plot Spectra data
5 """Classes to plot Spectra data
6
6
7 """
7 """
8
8
9 import os
9 import os
10 import numpy
10 import numpy
11 import datetime
11 import datetime
12
12
13 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 from schainpy.model.graphics.jroplot_base import Plot, plt, log
14 from itertools import combinations
14 from itertools import combinations
15 from matplotlib.ticker import LinearLocator
15 from matplotlib.ticker import LinearLocator
16
16
17 from schainpy.model.utils.BField import BField
17 from schainpy.model.utils.BField import BField
18 from scipy.interpolate import splrep
18 from scipy.interpolate import splrep
19 from scipy.interpolate import splev
19 from scipy.interpolate import splev
20
20
21 from matplotlib import __version__ as plt_version
21 from matplotlib import __version__ as plt_version
22
22
23 if plt_version >='3.3.4':
23 if plt_version >='3.3.4':
24 EXTRA_POINTS = 0
24 EXTRA_POINTS = 0
25 else:
25 else:
26 EXTRA_POINTS = 1
26 EXTRA_POINTS = 1
27 class SpectraPlot(Plot):
27 class SpectraPlot(Plot):
28 '''
28 '''
29 Plot for Spectra data
29 Plot for Spectra data
30 '''
30 '''
31
31
32 CODE = 'spc'
32 CODE = 'spc'
33 colormap = 'jet'
33 colormap = 'jet'
34 plot_type = 'pcolor'
34 plot_type = 'pcolor'
35 buffering = False
35 buffering = False
36 channelList = []
36 channelList = []
37 elevationList = []
37 elevationList = []
38 azimuthList = []
38 azimuthList = []
39
39
40 def setup(self):
40 def setup(self):
41
41
42 self.nplots = len(self.data.channels)
42 self.nplots = len(self.data.channels)
43 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
43 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
44 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
44 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
45 self.height = 3.4 * self.nrows
45 self.height = 3.4 * self.nrows
46 self.cb_label = 'dB'
46 self.cb_label = 'dB'
47 if self.showprofile:
47 if self.showprofile:
48 self.width = 5.2 * self.ncols
48 self.width = 5.2 * self.ncols
49 else:
49 else:
50 self.width = 4.2* self.ncols
50 self.width = 4.2* self.ncols
51 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.12})
51 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.12})
52 self.ylabel = 'Range [km]'
52 self.ylabel = 'Range [km]'
53
53
54 def update_list(self,dataOut):
54 def update_list(self,dataOut):
55
55
56 if len(self.channelList) == 0:
56 if len(self.channelList) == 0:
57 self.channelList = dataOut.channelList
57 self.channelList = dataOut.channelList
58 if len(self.elevationList) == 0:
58 if len(self.elevationList) == 0:
59 self.elevationList = dataOut.elevationList
59 self.elevationList = dataOut.elevationList
60 if len(self.azimuthList) == 0:
60 if len(self.azimuthList) == 0:
61 self.azimuthList = dataOut.azimuthList
61 self.azimuthList = dataOut.azimuthList
62
62
63 def update(self, dataOut):
63 def update(self, dataOut):
64
64
65 self.update_list(dataOut)
65 self.update_list(dataOut)
66 data = {}
66 data = {}
67 meta = {}
67 meta = {}
68
69 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
68 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
69 if dataOut.type == "Parameters":
70 noise = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
71 spc = 10*numpy.log10(dataOut.data_spc/(dataOut.nProfiles))
72 else:
70 noise = 10*numpy.log10(dataOut.getNoise()/norm)
73 noise = 10*numpy.log10(dataOut.getNoise()/norm)
74
71 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
75 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
72 for ch in range(dataOut.nChannels):
76 for ch in range(dataOut.nChannels):
73 if hasattr(dataOut.normFactor,'ndim'):
77 if hasattr(dataOut.normFactor,'ndim'):
74 if dataOut.normFactor.ndim > 1:
78 if dataOut.normFactor.ndim > 1:
75 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
79 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
76
80
77 else:
81 else:
78 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
82 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
79 else:
83 else:
80 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
84 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
81
85
82 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
86 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
83 spc = 10*numpy.log10(z)
87 spc = 10*numpy.log10(z)
84
88
85 data['spc'] = spc
89 data['spc'] = spc
86 data['rti'] = spc.mean(axis=1)
90 data['rti'] = spc.mean(axis=1)
87 data['noise'] = noise
91 data['noise'] = noise
88 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
92 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
89 if self.CODE == 'spc_moments':
93 if self.CODE == 'spc_moments':
90 data['moments'] = dataOut.moments
94 data['moments'] = dataOut.moments
91
95
92 return data, meta
96 return data, meta
93
97
94 def plot(self):
98 def plot(self):
95
99
96 if self.xaxis == "frequency":
100 if self.xaxis == "frequency":
97 x = self.data.xrange[0]
101 x = self.data.xrange[0]
98 self.xlabel = "Frequency (kHz)"
102 self.xlabel = "Frequency (kHz)"
99 elif self.xaxis == "time":
103 elif self.xaxis == "time":
100 x = self.data.xrange[1]
104 x = self.data.xrange[1]
101 self.xlabel = "Time (ms)"
105 self.xlabel = "Time (ms)"
102 else:
106 else:
103 x = self.data.xrange[2]
107 x = self.data.xrange[2]
104 self.xlabel = "Velocity (m/s)"
108 self.xlabel = "Velocity (m/s)"
105
109
106 if (self.CODE == 'spc_moments') | (self.CODE == 'gaussian_fit'):
110 if (self.CODE == 'spc_moments') | (self.CODE == 'gaussian_fit'):
107 x = self.data.xrange[2]
111 x = self.data.xrange[2]
108 self.xlabel = "Velocity (m/s)"
112 self.xlabel = "Velocity (m/s)"
109
113
110 self.titles = []
114 self.titles = []
111
115
112 y = self.data.yrange
116 y = self.data.yrange
113 self.y = y
117 self.y = y
114
118
115 data = self.data[-1]
119 data = self.data[-1]
116 z = data['spc']
120 z = data['spc']
117
121
118 for n, ax in enumerate(self.axes):
122 for n, ax in enumerate(self.axes):
119 noise = self.data['noise'][n][0]
123 noise = self.data['noise'][n][0]
120 # noise = data['noise'][n]
124 # noise = data['noise'][n]
121
125
122 if self.CODE == 'spc_moments':
126 if self.CODE == 'spc_moments':
123 mean = data['moments'][n, 1]
127 mean = data['moments'][n, 1]
124 if self.CODE == 'gaussian_fit':
128 if self.CODE == 'gaussian_fit':
125 gau0 = data['gaussfit'][n][2,:,0]
129 gau0 = data['gaussfit'][n][2,:,0]
126 gau1 = data['gaussfit'][n][2,:,1]
130 gau1 = data['gaussfit'][n][2,:,1]
127 if ax.firsttime:
131 if ax.firsttime:
128 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
132 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
129 self.xmin = self.xmin if self.xmin else -self.xmax
133 self.xmin = self.xmin if self.xmin else -self.xmax
130 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
134 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
131 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
135 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
132 ax.plt = ax.pcolormesh(x, y, z[n].T,
136 ax.plt = ax.pcolormesh(x, y, z[n].T,
133 vmin=self.zmin,
137 vmin=self.zmin,
134 vmax=self.zmax,
138 vmax=self.zmax,
135 cmap=plt.get_cmap(self.colormap)
139 cmap=plt.get_cmap(self.colormap)
136 )
140 )
137
141
138 if self.showprofile:
142 if self.showprofile:
139 ax.plt_profile = self.pf_axes[n].plot(
143 ax.plt_profile = self.pf_axes[n].plot(
140 data['rti'][n], y)[0]
144 data['rti'][n], y)[0]
141 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
145 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
142 color="k", linestyle="dashed", lw=1)[0]
146 color="k", linestyle="dashed", lw=1)[0]
143 if self.CODE == 'spc_moments':
147 if self.CODE == 'spc_moments':
144 ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
148 ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
145 if self.CODE == 'gaussian_fit':
149 if self.CODE == 'gaussian_fit':
146 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
150 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
147 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[0]
151 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[0]
148 else:
152 else:
149 ax.plt.set_array(z[n].T.ravel())
153 ax.plt.set_array(z[n].T.ravel())
150 if self.showprofile:
154 if self.showprofile:
151 ax.plt_profile.set_data(data['rti'][n], y)
155 ax.plt_profile.set_data(data['rti'][n], y)
152 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
156 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
153 if self.CODE == 'spc_moments':
157 if self.CODE == 'spc_moments':
154 ax.plt_mean.set_data(mean, y)
158 ax.plt_mean.set_data(mean, y)
155 if self.CODE == 'gaussian_fit':
159 if self.CODE == 'gaussian_fit':
156 ax.plt_gau0.set_data(gau0, y)
160 ax.plt_gau0.set_data(gau0, y)
157 ax.plt_gau1.set_data(gau1, y)
161 ax.plt_gau1.set_data(gau1, y)
158 if len(self.azimuthList) > 0 and len(self.elevationList) > 0:
162 if len(self.azimuthList) > 0 and len(self.elevationList) > 0:
159 self.titles.append('CH {}: {:2.1f}elv {:2.1f}az {:3.2f}dB'.format(self.channelList[n], noise, self.elevationList[n], self.azimuthList[n]))
163 self.titles.append('CH {}: {:2.1f}elv {:2.1f}az {:3.2f}dB'.format(self.channelList[n], noise, self.elevationList[n], self.azimuthList[n]))
160 else:
164 else:
161 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
165 self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
162
166
163 class SpectraObliquePlot(Plot):
167 class SpectraObliquePlot(Plot):
164 '''
168 '''
165 Plot for Spectra data
169 Plot for Spectra data
166 '''
170 '''
167
171
168 CODE = 'spc_oblique'
172 CODE = 'spc_oblique'
169 colormap = 'jet'
173 colormap = 'jet'
170 plot_type = 'pcolor'
174 plot_type = 'pcolor'
171
175
172 def setup(self):
176 def setup(self):
173 self.xaxis = "oblique"
177 self.xaxis = "oblique"
174 self.nplots = len(self.data.channels)
178 self.nplots = len(self.data.channels)
175 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
179 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
176 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
180 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
177 self.height = 2.6 * self.nrows
181 self.height = 2.6 * self.nrows
178 self.cb_label = 'dB'
182 self.cb_label = 'dB'
179 if self.showprofile:
183 if self.showprofile:
180 self.width = 4 * self.ncols
184 self.width = 4 * self.ncols
181 else:
185 else:
182 self.width = 3.5 * self.ncols
186 self.width = 3.5 * self.ncols
183 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
187 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
184 self.ylabel = 'Range [km]'
188 self.ylabel = 'Range [km]'
185
189
186 def update(self, dataOut):
190 def update(self, dataOut):
187
191
188 data = {}
192 data = {}
189 meta = {}
193 meta = {}
190
194
191 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
195 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
192 data['spc'] = spc
196 data['spc'] = spc
193 data['rti'] = dataOut.getPower()
197 data['rti'] = dataOut.getPower()
194 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
198 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
195 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
199 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
196
200
197 data['shift1'] = dataOut.Dop_EEJ_T1[0]
201 data['shift1'] = dataOut.Dop_EEJ_T1[0]
198 data['shift2'] = dataOut.Dop_EEJ_T2[0]
202 data['shift2'] = dataOut.Dop_EEJ_T2[0]
199 data['max_val_2'] = dataOut.Oblique_params[0,-1,:]
203 data['max_val_2'] = dataOut.Oblique_params[0,-1,:]
200 data['shift1_error'] = dataOut.Err_Dop_EEJ_T1[0]
204 data['shift1_error'] = dataOut.Err_Dop_EEJ_T1[0]
201 data['shift2_error'] = dataOut.Err_Dop_EEJ_T2[0]
205 data['shift2_error'] = dataOut.Err_Dop_EEJ_T2[0]
202
206
203 return data, meta
207 return data, meta
204
208
205 def plot(self):
209 def plot(self):
206
210
207 if self.xaxis == "frequency":
211 if self.xaxis == "frequency":
208 x = self.data.xrange[0]
212 x = self.data.xrange[0]
209 self.xlabel = "Frequency (kHz)"
213 self.xlabel = "Frequency (kHz)"
210 elif self.xaxis == "time":
214 elif self.xaxis == "time":
211 x = self.data.xrange[1]
215 x = self.data.xrange[1]
212 self.xlabel = "Time (ms)"
216 self.xlabel = "Time (ms)"
213 else:
217 else:
214 x = self.data.xrange[2]
218 x = self.data.xrange[2]
215 self.xlabel = "Velocity (m/s)"
219 self.xlabel = "Velocity (m/s)"
216
220
217 self.titles = []
221 self.titles = []
218
222
219 y = self.data.yrange
223 y = self.data.yrange
220 self.y = y
224 self.y = y
221
225
222 data = self.data[-1]
226 data = self.data[-1]
223 z = data['spc']
227 z = data['spc']
224
228
225 for n, ax in enumerate(self.axes):
229 for n, ax in enumerate(self.axes):
226 noise = self.data['noise'][n][-1]
230 noise = self.data['noise'][n][-1]
227 shift1 = data['shift1']
231 shift1 = data['shift1']
228 shift2 = data['shift2']
232 shift2 = data['shift2']
229 max_val_2 = data['max_val_2']
233 max_val_2 = data['max_val_2']
230 err1 = data['shift1_error']
234 err1 = data['shift1_error']
231 err2 = data['shift2_error']
235 err2 = data['shift2_error']
232 if ax.firsttime:
236 if ax.firsttime:
233
237
234 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
238 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
235 self.xmin = self.xmin if self.xmin else -self.xmax
239 self.xmin = self.xmin if self.xmin else -self.xmax
236 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
240 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
237 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
241 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
238 ax.plt = ax.pcolormesh(x, y, z[n].T,
242 ax.plt = ax.pcolormesh(x, y, z[n].T,
239 vmin=self.zmin,
243 vmin=self.zmin,
240 vmax=self.zmax,
244 vmax=self.zmax,
241 cmap=plt.get_cmap(self.colormap)
245 cmap=plt.get_cmap(self.colormap)
242 )
246 )
243
247
244 if self.showprofile:
248 if self.showprofile:
245 ax.plt_profile = self.pf_axes[n].plot(
249 ax.plt_profile = self.pf_axes[n].plot(
246 self.data['rti'][n][-1], y)[0]
250 self.data['rti'][n][-1], y)[0]
247 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
251 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
248 color="k", linestyle="dashed", lw=1)[0]
252 color="k", linestyle="dashed", lw=1)[0]
249
253
250 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
254 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
251 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
255 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
252 self.ploterr3 = ax.errorbar(max_val_2, y, xerr=0, fmt='g^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
256 self.ploterr3 = ax.errorbar(max_val_2, y, xerr=0, fmt='g^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
253
257
254 else:
258 else:
255 self.ploterr1.remove()
259 self.ploterr1.remove()
256 self.ploterr2.remove()
260 self.ploterr2.remove()
257 self.ploterr3.remove()
261 self.ploterr3.remove()
258 ax.plt.set_array(z[n].T.ravel())
262 ax.plt.set_array(z[n].T.ravel())
259 if self.showprofile:
263 if self.showprofile:
260 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
264 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
261 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
265 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
262 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
266 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=2.2, marker='o', linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
263 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
267 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
264 self.ploterr3 = ax.errorbar(max_val_2, y, xerr=0, fmt='g^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
268 self.ploterr3 = ax.errorbar(max_val_2, y, xerr=0, fmt='g^',elinewidth=2.2,marker='o',linestyle='None',markersize=2.5,capsize=0.3,markeredgewidth=0.2)
265
269
266 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
270 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
267
271
268
272
269 class CrossSpectraPlot(Plot):
273 class CrossSpectraPlot(Plot):
270
274
271 CODE = 'cspc'
275 CODE = 'cspc'
272 colormap = 'jet'
276 colormap = 'jet'
273 plot_type = 'pcolor'
277 plot_type = 'pcolor'
274 zmin_coh = None
278 zmin_coh = None
275 zmax_coh = None
279 zmax_coh = None
276 zmin_phase = None
280 zmin_phase = None
277 zmax_phase = None
281 zmax_phase = None
278 realChannels = None
282 realChannels = None
279 crossPairs = None
283 crossPairs = None
280
284
281 def setup(self):
285 def setup(self):
282
286
283 self.ncols = 4
287 self.ncols = 4
284 self.nplots = len(self.data.pairs) * 2
288 self.nplots = len(self.data.pairs) * 2
285 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
289 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
286 self.width = 3.1 * self.ncols
290 self.width = 3.1 * self.ncols
287 self.height = 2.6 * self.nrows
291 self.height = 2.6 * self.nrows
288 self.ylabel = 'Range [km]'
292 self.ylabel = 'Range [km]'
289 self.showprofile = False
293 self.showprofile = False
290 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
294 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
291
295
292 def update(self, dataOut):
296 def update(self, dataOut):
293
297
294 data = {}
298 data = {}
295 meta = {}
299 meta = {}
296
300
297 spc = dataOut.data_spc
301 spc = dataOut.data_spc
298 cspc = dataOut.data_cspc
302 cspc = dataOut.data_cspc
299 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
303 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
300 rawPairs = list(combinations(list(range(dataOut.nChannels)), 2))
304 rawPairs = list(combinations(list(range(dataOut.nChannels)), 2))
301 meta['pairs'] = rawPairs
305 meta['pairs'] = rawPairs
302 if self.crossPairs == None:
306 if self.crossPairs == None:
303 self.crossPairs = dataOut.pairsList
307 self.crossPairs = dataOut.pairsList
304 tmp = []
308 tmp = []
305
309
306 for n, pair in enumerate(meta['pairs']):
310 for n, pair in enumerate(meta['pairs']):
307 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
311 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
308 coh = numpy.abs(out)
312 coh = numpy.abs(out)
309 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
313 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
310 tmp.append(coh)
314 tmp.append(coh)
311 tmp.append(phase)
315 tmp.append(phase)
312
316
313 data['cspc'] = numpy.array(tmp)
317 data['cspc'] = numpy.array(tmp)
314
318
315 return data, meta
319 return data, meta
316
320
317 def plot(self):
321 def plot(self):
318
322
319 if self.xaxis == "frequency":
323 if self.xaxis == "frequency":
320 x = self.data.xrange[0]
324 x = self.data.xrange[0]
321 self.xlabel = "Frequency (kHz)"
325 self.xlabel = "Frequency (kHz)"
322 elif self.xaxis == "time":
326 elif self.xaxis == "time":
323 x = self.data.xrange[1]
327 x = self.data.xrange[1]
324 self.xlabel = "Time (ms)"
328 self.xlabel = "Time (ms)"
325 else:
329 else:
326 x = self.data.xrange[2]
330 x = self.data.xrange[2]
327 self.xlabel = "Velocity (m/s)"
331 self.xlabel = "Velocity (m/s)"
328
332
329 self.titles = []
333 self.titles = []
330
334
331 y = self.data.yrange
335 y = self.data.yrange
332 self.y = y
336 self.y = y
333
337
334 data = self.data[-1]
338 data = self.data[-1]
335 cspc = data['cspc']
339 cspc = data['cspc']
336
340
337 for n in range(len(self.data.pairs)):
341 for n in range(len(self.data.pairs)):
338 pair = self.crossPairs[n]
342 pair = self.crossPairs[n]
339 coh = cspc[n*2]
343 coh = cspc[n*2]
340 phase = cspc[n*2+1]
344 phase = cspc[n*2+1]
341 ax = self.axes[2 * n]
345 ax = self.axes[2 * n]
342 if ax.firsttime:
346 if ax.firsttime:
343 ax.plt = ax.pcolormesh(x, y, coh.T,
347 ax.plt = ax.pcolormesh(x, y, coh.T,
344 vmin=self.zmin_coh,
348 vmin=self.zmin_coh,
345 vmax=self.zmax_coh,
349 vmax=self.zmax_coh,
346 cmap=plt.get_cmap(self.colormap_coh)
350 cmap=plt.get_cmap(self.colormap_coh)
347 )
351 )
348 else:
352 else:
349 ax.plt.set_array(coh.T.ravel())
353 ax.plt.set_array(coh.T.ravel())
350 self.titles.append(
354 self.titles.append(
351 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
355 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
352
356
353 ax = self.axes[2 * n + 1]
357 ax = self.axes[2 * n + 1]
354 if ax.firsttime:
358 if ax.firsttime:
355 ax.plt = ax.pcolormesh(x, y, phase.T,
359 ax.plt = ax.pcolormesh(x, y, phase.T,
356 vmin=-180,
360 vmin=-180,
357 vmax=180,
361 vmax=180,
358 cmap=plt.get_cmap(self.colormap_phase)
362 cmap=plt.get_cmap(self.colormap_phase)
359 )
363 )
360 else:
364 else:
361 ax.plt.set_array(phase.T.ravel())
365 ax.plt.set_array(phase.T.ravel())
362 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
366 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
363
367
364
368
365 class CrossSpectra4Plot(Plot):
369 class CrossSpectra4Plot(Plot):
366
370
367 CODE = 'cspc'
371 CODE = 'cspc'
368 colormap = 'jet'
372 colormap = 'jet'
369 plot_type = 'pcolor'
373 plot_type = 'pcolor'
370 zmin_coh = None
374 zmin_coh = None
371 zmax_coh = None
375 zmax_coh = None
372 zmin_phase = None
376 zmin_phase = None
373 zmax_phase = None
377 zmax_phase = None
374
378
375 def setup(self):
379 def setup(self):
376
380
377 self.ncols = 4
381 self.ncols = 4
378 self.nrows = len(self.data.pairs)
382 self.nrows = len(self.data.pairs)
379 self.nplots = self.nrows * 4
383 self.nplots = self.nrows * 4
380 self.width = 3.1 * self.ncols
384 self.width = 3.1 * self.ncols
381 self.height = 5 * self.nrows
385 self.height = 5 * self.nrows
382 self.ylabel = 'Range [km]'
386 self.ylabel = 'Range [km]'
383 self.showprofile = False
387 self.showprofile = False
384 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
388 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
385
389
386 def plot(self):
390 def plot(self):
387
391
388 if self.xaxis == "frequency":
392 if self.xaxis == "frequency":
389 x = self.data.xrange[0]
393 x = self.data.xrange[0]
390 self.xlabel = "Frequency (kHz)"
394 self.xlabel = "Frequency (kHz)"
391 elif self.xaxis == "time":
395 elif self.xaxis == "time":
392 x = self.data.xrange[1]
396 x = self.data.xrange[1]
393 self.xlabel = "Time (ms)"
397 self.xlabel = "Time (ms)"
394 else:
398 else:
395 x = self.data.xrange[2]
399 x = self.data.xrange[2]
396 self.xlabel = "Velocity (m/s)"
400 self.xlabel = "Velocity (m/s)"
397
401
398 self.titles = []
402 self.titles = []
399
403
400
404
401 y = self.data.heights
405 y = self.data.heights
402 self.y = y
406 self.y = y
403 nspc = self.data['spc']
407 nspc = self.data['spc']
404 spc = self.data['cspc'][0]
408 spc = self.data['cspc'][0]
405 cspc = self.data['cspc'][1]
409 cspc = self.data['cspc'][1]
406
410
407 for n in range(self.nrows):
411 for n in range(self.nrows):
408 noise = self.data['noise'][:,-1]
412 noise = self.data['noise'][:,-1]
409 pair = self.data.pairs[n]
413 pair = self.data.pairs[n]
410
414
411 ax = self.axes[4 * n]
415 ax = self.axes[4 * n]
412 if ax.firsttime:
416 if ax.firsttime:
413 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
417 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
414 self.xmin = self.xmin if self.xmin else -self.xmax
418 self.xmin = self.xmin if self.xmin else -self.xmax
415 self.zmin = self.zmin if self.zmin else numpy.nanmin(nspc)
419 self.zmin = self.zmin if self.zmin else numpy.nanmin(nspc)
416 self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
420 self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
417 ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
421 ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
418 vmin=self.zmin,
422 vmin=self.zmin,
419 vmax=self.zmax,
423 vmax=self.zmax,
420 cmap=plt.get_cmap(self.colormap)
424 cmap=plt.get_cmap(self.colormap)
421 )
425 )
422 else:
426 else:
423
427
424 ax.plt.set_array(nspc[pair[0]].T.ravel())
428 ax.plt.set_array(nspc[pair[0]].T.ravel())
425 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
429 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
426
430
427 ax = self.axes[4 * n + 1]
431 ax = self.axes[4 * n + 1]
428
432
429 if ax.firsttime:
433 if ax.firsttime:
430 ax.plt = ax.pcolormesh(x , y, numpy.flip(nspc[pair[1]],axis=0).T,
434 ax.plt = ax.pcolormesh(x , y, numpy.flip(nspc[pair[1]],axis=0).T,
431 vmin=self.zmin,
435 vmin=self.zmin,
432 vmax=self.zmax,
436 vmax=self.zmax,
433 cmap=plt.get_cmap(self.colormap)
437 cmap=plt.get_cmap(self.colormap)
434 )
438 )
435 else:
439 else:
436
440
437 ax.plt.set_array(numpy.flip(nspc[pair[1]],axis=0).T.ravel())
441 ax.plt.set_array(numpy.flip(nspc[pair[1]],axis=0).T.ravel())
438 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
442 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
439
443
440 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
444 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
441 coh = numpy.abs(out)
445 coh = numpy.abs(out)
442 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
446 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
443
447
444 ax = self.axes[4 * n + 2]
448 ax = self.axes[4 * n + 2]
445 if ax.firsttime:
449 if ax.firsttime:
446 ax.plt = ax.pcolormesh(x, y, numpy.flip(coh,axis=0).T,
450 ax.plt = ax.pcolormesh(x, y, numpy.flip(coh,axis=0).T,
447 vmin=0,
451 vmin=0,
448 vmax=1,
452 vmax=1,
449 cmap=plt.get_cmap(self.colormap_coh)
453 cmap=plt.get_cmap(self.colormap_coh)
450 )
454 )
451 else:
455 else:
452 ax.plt.set_array(numpy.flip(coh,axis=0).T.ravel())
456 ax.plt.set_array(numpy.flip(coh,axis=0).T.ravel())
453 self.titles.append(
457 self.titles.append(
454 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
458 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
455
459
456 ax = self.axes[4 * n + 3]
460 ax = self.axes[4 * n + 3]
457 if ax.firsttime:
461 if ax.firsttime:
458 ax.plt = ax.pcolormesh(x, y, numpy.flip(phase,axis=0).T,
462 ax.plt = ax.pcolormesh(x, y, numpy.flip(phase,axis=0).T,
459 vmin=-180,
463 vmin=-180,
460 vmax=180,
464 vmax=180,
461 cmap=plt.get_cmap(self.colormap_phase)
465 cmap=plt.get_cmap(self.colormap_phase)
462 )
466 )
463 else:
467 else:
464 ax.plt.set_array(numpy.flip(phase,axis=0).T.ravel())
468 ax.plt.set_array(numpy.flip(phase,axis=0).T.ravel())
465 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
469 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
466
470
467
471
468 class CrossSpectra2Plot(Plot):
472 class CrossSpectra2Plot(Plot):
469
473
470 CODE = 'cspc'
474 CODE = 'cspc'
471 colormap = 'jet'
475 colormap = 'jet'
472 plot_type = 'pcolor'
476 plot_type = 'pcolor'
473 zmin_coh = None
477 zmin_coh = None
474 zmax_coh = None
478 zmax_coh = None
475 zmin_phase = None
479 zmin_phase = None
476 zmax_phase = None
480 zmax_phase = None
477
481
478 def setup(self):
482 def setup(self):
479
483
480 self.ncols = 1
484 self.ncols = 1
481 self.nrows = len(self.data.pairs)
485 self.nrows = len(self.data.pairs)
482 self.nplots = self.nrows * 1
486 self.nplots = self.nrows * 1
483 self.width = 3.1 * self.ncols
487 self.width = 3.1 * self.ncols
484 self.height = 5 * self.nrows
488 self.height = 5 * self.nrows
485 self.ylabel = 'Range [km]'
489 self.ylabel = 'Range [km]'
486 self.showprofile = False
490 self.showprofile = False
487 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
491 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
488
492
489 def plot(self):
493 def plot(self):
490
494
491 if self.xaxis == "frequency":
495 if self.xaxis == "frequency":
492 x = self.data.xrange[0]
496 x = self.data.xrange[0]
493 self.xlabel = "Frequency (kHz)"
497 self.xlabel = "Frequency (kHz)"
494 elif self.xaxis == "time":
498 elif self.xaxis == "time":
495 x = self.data.xrange[1]
499 x = self.data.xrange[1]
496 self.xlabel = "Time (ms)"
500 self.xlabel = "Time (ms)"
497 else:
501 else:
498 x = self.data.xrange[2]
502 x = self.data.xrange[2]
499 self.xlabel = "Velocity (m/s)"
503 self.xlabel = "Velocity (m/s)"
500
504
501 self.titles = []
505 self.titles = []
502
506
503
507
504 y = self.data.heights
508 y = self.data.heights
505 self.y = y
509 self.y = y
506 cspc = self.data['cspc'][1]
510 cspc = self.data['cspc'][1]
507
511
508 for n in range(self.nrows):
512 for n in range(self.nrows):
509 noise = self.data['noise'][:,-1]
513 noise = self.data['noise'][:,-1]
510 pair = self.data.pairs[n]
514 pair = self.data.pairs[n]
511 out = cspc[n]
515 out = cspc[n]
512 cross = numpy.abs(out)
516 cross = numpy.abs(out)
513 z = cross/self.data.nFactor
517 z = cross/self.data.nFactor
514 cross = 10*numpy.log10(z)
518 cross = 10*numpy.log10(z)
515
519
516 ax = self.axes[1 * n]
520 ax = self.axes[1 * n]
517 if ax.firsttime:
521 if ax.firsttime:
518 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
522 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
519 self.xmin = self.xmin if self.xmin else -self.xmax
523 self.xmin = self.xmin if self.xmin else -self.xmax
520 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
524 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
521 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
525 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
522 ax.plt = ax.pcolormesh(x, y, cross.T,
526 ax.plt = ax.pcolormesh(x, y, cross.T,
523 vmin=self.zmin,
527 vmin=self.zmin,
524 vmax=self.zmax,
528 vmax=self.zmax,
525 cmap=plt.get_cmap(self.colormap)
529 cmap=plt.get_cmap(self.colormap)
526 )
530 )
527 else:
531 else:
528 ax.plt.set_array(cross.T.ravel())
532 ax.plt.set_array(cross.T.ravel())
529 self.titles.append(
533 self.titles.append(
530 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
534 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
531
535
532
536
533 class CrossSpectra3Plot(Plot):
537 class CrossSpectra3Plot(Plot):
534
538
535 CODE = 'cspc'
539 CODE = 'cspc'
536 colormap = 'jet'
540 colormap = 'jet'
537 plot_type = 'pcolor'
541 plot_type = 'pcolor'
538 zmin_coh = None
542 zmin_coh = None
539 zmax_coh = None
543 zmax_coh = None
540 zmin_phase = None
544 zmin_phase = None
541 zmax_phase = None
545 zmax_phase = None
542
546
543 def setup(self):
547 def setup(self):
544
548
545 self.ncols = 3
549 self.ncols = 3
546 self.nrows = len(self.data.pairs)
550 self.nrows = len(self.data.pairs)
547 self.nplots = self.nrows * 3
551 self.nplots = self.nrows * 3
548 self.width = 3.1 * self.ncols
552 self.width = 3.1 * self.ncols
549 self.height = 5 * self.nrows
553 self.height = 5 * self.nrows
550 self.ylabel = 'Range [km]'
554 self.ylabel = 'Range [km]'
551 self.showprofile = False
555 self.showprofile = False
552 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
556 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
553
557
554 def plot(self):
558 def plot(self):
555
559
556 if self.xaxis == "frequency":
560 if self.xaxis == "frequency":
557 x = self.data.xrange[0]
561 x = self.data.xrange[0]
558 self.xlabel = "Frequency (kHz)"
562 self.xlabel = "Frequency (kHz)"
559 elif self.xaxis == "time":
563 elif self.xaxis == "time":
560 x = self.data.xrange[1]
564 x = self.data.xrange[1]
561 self.xlabel = "Time (ms)"
565 self.xlabel = "Time (ms)"
562 else:
566 else:
563 x = self.data.xrange[2]
567 x = self.data.xrange[2]
564 self.xlabel = "Velocity (m/s)"
568 self.xlabel = "Velocity (m/s)"
565
569
566 self.titles = []
570 self.titles = []
567
571
568
572
569 y = self.data.heights
573 y = self.data.heights
570 self.y = y
574 self.y = y
571
575
572 cspc = self.data['cspc'][1]
576 cspc = self.data['cspc'][1]
573
577
574 for n in range(self.nrows):
578 for n in range(self.nrows):
575 noise = self.data['noise'][:,-1]
579 noise = self.data['noise'][:,-1]
576 pair = self.data.pairs[n]
580 pair = self.data.pairs[n]
577 out = cspc[n]
581 out = cspc[n]
578
582
579 cross = numpy.abs(out)
583 cross = numpy.abs(out)
580 z = cross/self.data.nFactor
584 z = cross/self.data.nFactor
581 cross = 10*numpy.log10(z)
585 cross = 10*numpy.log10(z)
582
586
583 out_r= out.real/self.data.nFactor
587 out_r= out.real/self.data.nFactor
584
588
585 out_i= out.imag/self.data.nFactor
589 out_i= out.imag/self.data.nFactor
586
590
587 ax = self.axes[3 * n]
591 ax = self.axes[3 * n]
588 if ax.firsttime:
592 if ax.firsttime:
589 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
593 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
590 self.xmin = self.xmin if self.xmin else -self.xmax
594 self.xmin = self.xmin if self.xmin else -self.xmax
591 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
595 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
592 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
596 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
593 ax.plt = ax.pcolormesh(x, y, cross.T,
597 ax.plt = ax.pcolormesh(x, y, cross.T,
594 vmin=self.zmin,
598 vmin=self.zmin,
595 vmax=self.zmax,
599 vmax=self.zmax,
596 cmap=plt.get_cmap(self.colormap)
600 cmap=plt.get_cmap(self.colormap)
597 )
601 )
598 else:
602 else:
599 ax.plt.set_array(cross.T.ravel())
603 ax.plt.set_array(cross.T.ravel())
600 self.titles.append(
604 self.titles.append(
601 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
605 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
602
606
603 ax = self.axes[3 * n + 1]
607 ax = self.axes[3 * n + 1]
604 if ax.firsttime:
608 if ax.firsttime:
605 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
609 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
606 self.xmin = self.xmin if self.xmin else -self.xmax
610 self.xmin = self.xmin if self.xmin else -self.xmax
607 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
611 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
608 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
612 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
609 ax.plt = ax.pcolormesh(x, y, out_r.T,
613 ax.plt = ax.pcolormesh(x, y, out_r.T,
610 vmin=-1.e6,
614 vmin=-1.e6,
611 vmax=0,
615 vmax=0,
612 cmap=plt.get_cmap(self.colormap)
616 cmap=plt.get_cmap(self.colormap)
613 )
617 )
614 else:
618 else:
615 ax.plt.set_array(out_r.T.ravel())
619 ax.plt.set_array(out_r.T.ravel())
616 self.titles.append(
620 self.titles.append(
617 'Cross Spectra Real Ch{} * Ch{}'.format(pair[0], pair[1]))
621 'Cross Spectra Real Ch{} * Ch{}'.format(pair[0], pair[1]))
618
622
619 ax = self.axes[3 * n + 2]
623 ax = self.axes[3 * n + 2]
620
624
621
625
622 if ax.firsttime:
626 if ax.firsttime:
623 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
627 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
624 self.xmin = self.xmin if self.xmin else -self.xmax
628 self.xmin = self.xmin if self.xmin else -self.xmax
625 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
629 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
626 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
630 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
627 ax.plt = ax.pcolormesh(x, y, out_i.T,
631 ax.plt = ax.pcolormesh(x, y, out_i.T,
628 vmin=-1.e6,
632 vmin=-1.e6,
629 vmax=1.e6,
633 vmax=1.e6,
630 cmap=plt.get_cmap(self.colormap)
634 cmap=plt.get_cmap(self.colormap)
631 )
635 )
632 else:
636 else:
633 ax.plt.set_array(out_i.T.ravel())
637 ax.plt.set_array(out_i.T.ravel())
634 self.titles.append(
638 self.titles.append(
635 'Cross Spectra Imag Ch{} * Ch{}'.format(pair[0], pair[1]))
639 'Cross Spectra Imag Ch{} * Ch{}'.format(pair[0], pair[1]))
636
640
637 class RTIPlot(Plot):
641 class RTIPlot(Plot):
638 '''
642 '''
639 Plot for RTI data
643 Plot for RTI data
640 '''
644 '''
641
645
642 CODE = 'rti'
646 CODE = 'rti'
643 colormap = 'jet'
647 colormap = 'jet'
644 plot_type = 'pcolorbuffer'
648 plot_type = 'pcolorbuffer'
645 titles = None
649 titles = None
646 channelList = []
650 channelList = []
647 elevationList = []
651 elevationList = []
648 azimuthList = []
652 azimuthList = []
649
653
650 def setup(self):
654 def setup(self):
651 self.xaxis = 'time'
655 self.xaxis = 'time'
652 self.ncols = 1
656 self.ncols = 1
653 self.nrows = len(self.data.channels)
657 self.nrows = len(self.data.channels)
654 self.nplots = len(self.data.channels)
658 self.nplots = len(self.data.channels)
655 self.ylabel = 'Range [km]'
659 self.ylabel = 'Range [km]'
656 #self.xlabel = 'Time'
660 #self.xlabel = 'Time'
657 self.cb_label = 'dB'
661 self.cb_label = 'dB'
658 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
662 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
659 self.titles = ['{} Channel {}'.format(
663 self.titles = ['{} Channel {}'.format(
660 self.CODE.upper(), x) for x in range(self.nplots)]
664 self.CODE.upper(), x) for x in range(self.nplots)]
661
665
662 def update_list(self,dataOut):
666 def update_list(self,dataOut):
663
667
664 if len(self.channelList) == 0:
668 if len(self.channelList) == 0:
665 self.channelList = dataOut.channelList
669 self.channelList = dataOut.channelList
666 if len(self.elevationList) == 0:
670 if len(self.elevationList) == 0:
667 self.elevationList = dataOut.elevationList
671 self.elevationList = dataOut.elevationList
668 if len(self.azimuthList) == 0:
672 if len(self.azimuthList) == 0:
669 self.azimuthList = dataOut.azimuthList
673 self.azimuthList = dataOut.azimuthList
670
674
671
675
672 def update(self, dataOut):
676 def update(self, dataOut):
673
677
674 if len(self.channelList) == 0:
678 if len(self.channelList) == 0:
675 self.update_list(dataOut)
679 self.update_list(dataOut)
676 data = {}
680 data = {}
677 meta = {}
681 meta = {}
678 data['rti'] = dataOut.getPower()
682 data['rti'] = dataOut.getPower()
679 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
683 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
680 noise = 10*numpy.log10(dataOut.getNoise()/norm)
684 noise = 10*numpy.log10(dataOut.getNoise()/norm)
681 data['noise'] = noise
685 data['noise'] = noise
682
686
683 return data, meta
687 return data, meta
684
688
685 def plot(self):
689 def plot(self):
686
690
687 self.x = self.data.times
691 self.x = self.data.times
688 self.y = self.data.yrange
692 self.y = self.data.yrange
689 self.z = self.data[self.CODE]
693 self.z = self.data[self.CODE]
690 self.z = numpy.array(self.z, dtype=float)
694 self.z = numpy.array(self.z, dtype=float)
691 self.z = numpy.ma.masked_invalid(self.z)
695 self.z = numpy.ma.masked_invalid(self.z)
692
696
693 try:
697 try:
694 if self.channelList != None:
698 if self.channelList != None:
695 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
699 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
696 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
700 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
697 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
701 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
698 else:
702 else:
699 self.titles = ['{} Channel {}'.format(
703 self.titles = ['{} Channel {}'.format(
700 self.CODE.upper(), x) for x in self.channelList]
704 self.CODE.upper(), x) for x in self.channelList]
701 except:
705 except:
702 if self.channelList.any() != None:
706 if self.channelList.any() != None:
703 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
707 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
704 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
708 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
705 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
709 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
706 else:
710 else:
707 self.titles = ['{} Channel {}'.format(
711 self.titles = ['{} Channel {}'.format(
708 self.CODE.upper(), x) for x in self.channelList]
712 self.CODE.upper(), x) for x in self.channelList]
709
713
710 if self.decimation is None:
714 if self.decimation is None:
711 x, y, z = self.fill_gaps(self.x, self.y, self.z)
715 x, y, z = self.fill_gaps(self.x, self.y, self.z)
712 else:
716 else:
713 x, y, z = self.fill_gaps(*self.decimate())
717 x, y, z = self.fill_gaps(*self.decimate())
714
718
715 for n, ax in enumerate(self.axes):
719 for n, ax in enumerate(self.axes):
716
720
717 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
721 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
718 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
722 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
719 data = self.data[-1]
723 data = self.data[-1]
720 if ax.firsttime:
724 if ax.firsttime:
721 ax.plt = ax.pcolormesh(x, y, z[n].T,
725 ax.plt = ax.pcolormesh(x, y, z[n].T,
722 vmin=self.zmin,
726 vmin=self.zmin,
723 vmax=self.zmax,
727 vmax=self.zmax,
724 cmap=plt.get_cmap(self.colormap)
728 cmap=plt.get_cmap(self.colormap)
725 )
729 )
726 if self.showprofile:
730 if self.showprofile:
727 ax.plot_profile = self.pf_axes[n].plot(
731 ax.plot_profile = self.pf_axes[n].plot(
728 data['rti'][n], self.y)[0]
732 data[self.CODE][n], self.y)[0]
729 if "noise" in self.data:
733 if "noise" in self.data:
730 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
734 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
731 color="k", linestyle="dashed", lw=1)[0]
735 color="k", linestyle="dashed", lw=1)[0]
732 else:
736 else:
733 # ax.collections.remove(ax.collections[0]) # error while running
737 # ax.collections.remove(ax.collections[0]) # error while running
734 ax.plt = ax.pcolormesh(x, y, z[n].T,
738 ax.plt = ax.pcolormesh(x, y, z[n].T,
735 vmin=self.zmin,
739 vmin=self.zmin,
736 vmax=self.zmax,
740 vmax=self.zmax,
737 cmap=plt.get_cmap(self.colormap)
741 cmap=plt.get_cmap(self.colormap)
738 )
742 )
739 if self.showprofile:
743 if self.showprofile:
740 ax.plot_profile.set_data(data['rti'][n], self.y)
744 ax.plot_profile.set_data(data[self.CODE][n], self.y)
741 if "noise" in self.data:
745 if "noise" in self.data:
742 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
746 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
743 color="k", linestyle="dashed", lw=1)[0]
747 color="k", linestyle="dashed", lw=1)[0]
744
748
745 class SpectrogramPlot(Plot):
749 class SpectrogramPlot(Plot):
746 '''
750 '''
747 Plot for Spectrogram data
751 Plot for Spectrogram data
748 '''
752 '''
749
753
750 CODE = 'Spectrogram_Profile'
754 CODE = 'Spectrogram_Profile'
751 colormap = 'binary'
755 colormap = 'binary'
752 plot_type = 'pcolorbuffer'
756 plot_type = 'pcolorbuffer'
753
757
754 def setup(self):
758 def setup(self):
755 self.xaxis = 'time'
759 self.xaxis = 'time'
756 self.ncols = 1
760 self.ncols = 1
757 self.nrows = len(self.data.channels)
761 self.nrows = len(self.data.channels)
758 self.nplots = len(self.data.channels)
762 self.nplots = len(self.data.channels)
759 self.xlabel = 'Time'
763 self.xlabel = 'Time'
760 self.plots_adjust.update({'hspace':1.2, 'left': 0.1, 'bottom': 0.12, 'right':0.95})
764 self.plots_adjust.update({'hspace':1.2, 'left': 0.1, 'bottom': 0.12, 'right':0.95})
761 self.titles = []
765 self.titles = []
762
766
763 self.titles = ['{} Channel {}'.format(
767 self.titles = ['{} Channel {}'.format(
764 self.CODE.upper(), x) for x in range(self.nrows)]
768 self.CODE.upper(), x) for x in range(self.nrows)]
765
769
766
770
767 def update(self, dataOut):
771 def update(self, dataOut):
768 data = {}
772 data = {}
769 meta = {}
773 meta = {}
770
774
771 maxHei = 1620#+12000
775 maxHei = 1620#+12000
772 indb = numpy.where(dataOut.heightList <= maxHei)
776 indb = numpy.where(dataOut.heightList <= maxHei)
773 hei = indb[0][-1]
777 hei = indb[0][-1]
774
778
775 factor = dataOut.nIncohInt
779 factor = dataOut.nIncohInt
776 z = dataOut.data_spc[:,:,hei] / factor
780 z = dataOut.data_spc[:,:,hei] / factor
777 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
781 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
778
782
779 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
783 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
780 data['Spectrogram_Profile'] = 10 * numpy.log10(z)
784 data['Spectrogram_Profile'] = 10 * numpy.log10(z)
781
785
782 data['hei'] = hei
786 data['hei'] = hei
783 data['DH'] = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
787 data['DH'] = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
784 data['nProfiles'] = dataOut.nProfiles
788 data['nProfiles'] = dataOut.nProfiles
785
789
786 return data, meta
790 return data, meta
787
791
788 def plot(self):
792 def plot(self):
789
793
790 self.x = self.data.times
794 self.x = self.data.times
791 self.z = self.data[self.CODE]
795 self.z = self.data[self.CODE]
792 self.y = self.data.xrange[0]
796 self.y = self.data.xrange[0]
793
797
794 hei = self.data['hei'][-1]
798 hei = self.data['hei'][-1]
795 DH = self.data['DH'][-1]
799 DH = self.data['DH'][-1]
796 nProfiles = self.data['nProfiles'][-1]
800 nProfiles = self.data['nProfiles'][-1]
797
801
798 self.ylabel = "Frequency (kHz)"
802 self.ylabel = "Frequency (kHz)"
799
803
800 self.z = numpy.ma.masked_invalid(self.z)
804 self.z = numpy.ma.masked_invalid(self.z)
801
805
802 if self.decimation is None:
806 if self.decimation is None:
803 x, y, z = self.fill_gaps(self.x, self.y, self.z)
807 x, y, z = self.fill_gaps(self.x, self.y, self.z)
804 else:
808 else:
805 x, y, z = self.fill_gaps(*self.decimate())
809 x, y, z = self.fill_gaps(*self.decimate())
806
810
807 for n, ax in enumerate(self.axes):
811 for n, ax in enumerate(self.axes):
808 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
812 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
809 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
813 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
810 data = self.data[-1]
814 data = self.data[-1]
811 if ax.firsttime:
815 if ax.firsttime:
812 ax.plt = ax.pcolormesh(x, y, z[n].T,
816 ax.plt = ax.pcolormesh(x, y, z[n].T,
813 vmin=self.zmin,
817 vmin=self.zmin,
814 vmax=self.zmax,
818 vmax=self.zmax,
815 cmap=plt.get_cmap(self.colormap)
819 cmap=plt.get_cmap(self.colormap)
816 )
820 )
817 else:
821 else:
818 # ax.collections.remove(ax.collections[0]) # error while running
822 # ax.collections.remove(ax.collections[0]) # error while running
819 ax.plt = ax.pcolormesh(x, y, z[n].T,
823 ax.plt = ax.pcolormesh(x, y, z[n].T,
820 vmin=self.zmin,
824 vmin=self.zmin,
821 vmax=self.zmax,
825 vmax=self.zmax,
822 cmap=plt.get_cmap(self.colormap)
826 cmap=plt.get_cmap(self.colormap)
823 )
827 )
824
828
825
829
826
830
827 class CoherencePlot(RTIPlot):
831 class CoherencePlot(RTIPlot):
828 '''
832 '''
829 Plot for Coherence data
833 Plot for Coherence data
830 '''
834 '''
831
835
832 CODE = 'coh'
836 CODE = 'coh'
833 titles = None
837 titles = None
834
838
835 def setup(self):
839 def setup(self):
836 self.xaxis = 'time'
840 self.xaxis = 'time'
837 self.ncols = 1
841 self.ncols = 1
838 self.nrows = len(self.data.pairs)
842 self.nrows = len(self.data.pairs)
839 self.nplots = len(self.data.pairs)
843 self.nplots = len(self.data.pairs)
840 self.ylabel = 'Range [km]'
844 self.ylabel = 'Range [km]'
841 self.xlabel = 'Time'
845 self.xlabel = 'Time'
842 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
846 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
843 if self.CODE == 'coh':
847 if self.CODE == 'coh':
844 self.cb_label = ''
848 self.cb_label = ''
845 self.titles = [
849 self.titles = [
846 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
850 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
847 else:
851 else:
848 self.cb_label = 'Degrees'
852 self.cb_label = 'Degrees'
849 self.titles = [
853 self.titles = [
850 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
854 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
851
855
852 def update(self, dataOut):
856 def update(self, dataOut):
853
857
854 data = {}
858 data = {}
855 meta = {}
859 meta = {}
856 data['coh'] = dataOut.getCoherence()
860 data['coh'] = dataOut.getCoherence()
857 meta['pairs'] = dataOut.pairsList
861 meta['pairs'] = dataOut.pairsList
858
862
859 return data, meta
863 return data, meta
860
864
861 class PhasePlot(CoherencePlot):
865 class PhasePlot(CoherencePlot):
862 '''
866 '''
863 Plot for Phase map data
867 Plot for Phase map data
864 '''
868 '''
865
869
866 CODE = 'phase'
870 CODE = 'phase'
867 colormap = 'seismic'
871 colormap = 'seismic'
868
872
869 def update(self, dataOut):
873 def update(self, dataOut):
870
874
871 data = {}
875 data = {}
872 meta = {}
876 meta = {}
873 data['phase'] = dataOut.getCoherence(phase=True)
877 data['phase'] = dataOut.getCoherence(phase=True)
874 meta['pairs'] = dataOut.pairsList
878 meta['pairs'] = dataOut.pairsList
875
879
876 return data, meta
880 return data, meta
877
881
878 class NoisePlot(Plot):
882 class NoisePlot(Plot):
879 '''
883 '''
880 Plot for noise
884 Plot for noise
881 '''
885 '''
882
886
883 CODE = 'noise'
887 CODE = 'noise'
884 plot_type = 'scatterbuffer'
888 plot_type = 'scatterbuffer'
885
889
886 def setup(self):
890 def setup(self):
887 self.xaxis = 'time'
891 self.xaxis = 'time'
888 self.ncols = 1
892 self.ncols = 1
889 self.nrows = 1
893 self.nrows = 1
890 self.nplots = 1
894 self.nplots = 1
891 self.ylabel = 'Intensity [dB]'
895 self.ylabel = 'Intensity [dB]'
892 self.xlabel = 'Time'
896 self.xlabel = 'Time'
893 self.titles = ['Noise']
897 self.titles = ['Noise']
894 self.colorbar = False
898 self.colorbar = False
895 self.plots_adjust.update({'right': 0.85 })
899 self.plots_adjust.update({'right': 0.85 })
900 self.titles = ['Noise Plot']
896
901
897 def update(self, dataOut):
902 def update(self, dataOut):
898
903
899 data = {}
904 data = {}
900 meta = {}
905 meta = {}
901 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
906 noise = 10*numpy.log10(dataOut.getNoise())
907 noise = noise.reshape(dataOut.nChannels, 1)
908 data['noise'] = noise
902 meta['yrange'] = numpy.array([])
909 meta['yrange'] = numpy.array([])
903
910
904 return data, meta
911 return data, meta
905
912
906 def plot(self):
913 def plot(self):
907
914
908 x = self.data.times
915 x = self.data.times
909 xmin = self.data.min_time
916 xmin = self.data.min_time
910 xmax = xmin + self.xrange * 60 * 60
917 xmax = xmin + self.xrange * 60 * 60
911 Y = self.data['noise']
918 Y = self.data['noise']
912
919
913 if self.axes[0].firsttime:
920 if self.axes[0].firsttime:
914 self.ymin = numpy.nanmin(Y) - 5
921 self.ymin = numpy.nanmin(Y) - 5
915 self.ymax = numpy.nanmax(Y) + 5
922 self.ymax = numpy.nanmax(Y) + 5
916 for ch in self.data.channels:
923 for ch in self.data.channels:
917 y = Y[ch]
924 y = Y[ch]
918 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
925 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
919 plt.legend(bbox_to_anchor=(1.18, 1.0))
926 plt.legend(bbox_to_anchor=(1.18, 1.0))
920 else:
927 else:
921 for ch in self.data.channels:
928 for ch in self.data.channels:
922 y = Y[ch]
929 y = Y[ch]
923 self.axes[0].lines[ch].set_data(x, y)
930 self.axes[0].lines[ch].set_data(x, y)
924
931
925 class PowerProfilePlot(Plot):
932 class PowerProfilePlot(Plot):
926
933
927 CODE = 'pow_profile'
934 CODE = 'pow_profile'
928 plot_type = 'scatter'
935 plot_type = 'scatter'
929
936
930 def setup(self):
937 def setup(self):
931
938
932 self.ncols = 1
939 self.ncols = 1
933 self.nrows = 1
940 self.nrows = 1
934 self.nplots = 1
941 self.nplots = 1
935 self.height = 4
942 self.height = 4
936 self.width = 3
943 self.width = 3
937 self.ylabel = 'Range [km]'
944 self.ylabel = 'Range [km]'
938 self.xlabel = 'Intensity [dB]'
945 self.xlabel = 'Intensity [dB]'
939 self.titles = ['Power Profile']
946 self.titles = ['Power Profile']
940 self.colorbar = False
947 self.colorbar = False
941
948
942 def update(self, dataOut):
949 def update(self, dataOut):
943
950
944 data = {}
951 data = {}
945 meta = {}
952 meta = {}
946 data[self.CODE] = dataOut.getPower()
953 data[self.CODE] = dataOut.getPower()
947
954
948 return data, meta
955 return data, meta
949
956
950 def plot(self):
957 def plot(self):
951
958
952 y = self.data.yrange
959 y = self.data.yrange
953 self.y = y
960 self.y = y
954
961
955 x = self.data[-1][self.CODE]
962 x = self.data[-1][self.CODE]
956
963
957 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
964 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
958 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
965 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
959
966
960 if self.axes[0].firsttime:
967 if self.axes[0].firsttime:
961 for ch in self.data.channels:
968 for ch in self.data.channels:
962 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
969 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
963 plt.legend()
970 plt.legend()
964 else:
971 else:
965 for ch in self.data.channels:
972 for ch in self.data.channels:
966 self.axes[0].lines[ch].set_data(x[ch], y)
973 self.axes[0].lines[ch].set_data(x[ch], y)
967
974
968
975
969 class SpectraCutPlot(Plot):
976 class SpectraCutPlot(Plot):
970
977
971 CODE = 'spc_cut'
978 CODE = 'spc_cut'
972 plot_type = 'scatter'
979 plot_type = 'scatter'
973 buffering = False
980 buffering = False
974 heights = []
981 heights = []
975 channelList = []
982 channelList = []
976 maintitle = "Spectra Cuts"
983 maintitle = "Spectra Cuts"
977 flag_setIndex = False
984 flag_setIndex = False
978
985
979 def setup(self):
986 def setup(self):
980
987
981 self.nplots = len(self.data.channels)
988 self.nplots = len(self.data.channels)
982 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
989 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
983 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
990 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
984 self.width = 4.5 * self.ncols + 2.5
991 self.width = 4.5 * self.ncols + 2.5
985 self.height = 4.8 * self.nrows
992 self.height = 4.8 * self.nrows
986 self.ylabel = 'Power [dB]'
993 self.ylabel = 'Power [dB]'
987 self.colorbar = False
994 self.colorbar = False
988 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.9, 'bottom':0.08})
995 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.9, 'bottom':0.08})
989
996
990 if len(self.selectedHeightsList) > 0:
997 if len(self.selectedHeightsList) > 0:
991 self.maintitle = "Spectra Cut"# for %d km " %(int(self.selectedHeight))
998 self.maintitle = "Spectra Cut"# for %d km " %(int(self.selectedHeight))
992
999
993
1000
994
1001
995 def update(self, dataOut):
1002 def update(self, dataOut):
996 if len(self.channelList) == 0:
1003 if len(self.channelList) == 0:
997 self.channelList = dataOut.channelList
1004 self.channelList = dataOut.channelList
998
1005
999 self.heights = dataOut.heightList
1006 self.heights = dataOut.heightList
1000 #print("sels: ",self.selectedHeightsList)
1007 #print("sels: ",self.selectedHeightsList)
1001 if len(self.selectedHeightsList)>0 and not self.flag_setIndex:
1008 if len(self.selectedHeightsList)>0 and not self.flag_setIndex:
1002
1009
1003 for sel_height in self.selectedHeightsList:
1010 for sel_height in self.selectedHeightsList:
1004 index_list = numpy.where(self.heights >= sel_height)
1011 index_list = numpy.where(self.heights >= sel_height)
1005 index_list = index_list[0]
1012 index_list = index_list[0]
1006 self.height_index.append(index_list[0])
1013 self.height_index.append(index_list[0])
1007 #print("sels i:"", self.height_index)
1014 #print("sels i:"", self.height_index)
1008 self.flag_setIndex = True
1015 self.flag_setIndex = True
1009 #print(self.height_index)
1016 #print(self.height_index)
1010 data = {}
1017 data = {}
1011 meta = {}
1018 meta = {}
1012
1019
1013 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter#*dataOut.nFFTPoints
1020 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter#*dataOut.nFFTPoints
1014 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1021 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1015 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
1022 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
1016
1023
1017
1024
1018 z = []
1025 z = []
1019 for ch in range(dataOut.nChannels):
1026 for ch in range(dataOut.nChannels):
1020 if hasattr(dataOut.normFactor,'shape'):
1027 if hasattr(dataOut.normFactor,'shape'):
1021 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
1028 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
1022 else:
1029 else:
1023 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1030 z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1024
1031
1025 z = numpy.asarray(z)
1032 z = numpy.asarray(z)
1026 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1033 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1027 spc = 10*numpy.log10(z)
1034 spc = 10*numpy.log10(z)
1028
1035
1029
1036
1030 data['spc'] = spc - noise
1037 data['spc'] = spc - noise
1031 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
1038 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
1032
1039
1033 return data, meta
1040 return data, meta
1034
1041
1035 def plot(self):
1042 def plot(self):
1036 if self.xaxis == "frequency":
1043 if self.xaxis == "frequency":
1037 x = self.data.xrange[0][0:]
1044 x = self.data.xrange[0][0:]
1038 self.xlabel = "Frequency (kHz)"
1045 self.xlabel = "Frequency (kHz)"
1039 elif self.xaxis == "time":
1046 elif self.xaxis == "time":
1040 x = self.data.xrange[1]
1047 x = self.data.xrange[1]
1041 self.xlabel = "Time (ms)"
1048 self.xlabel = "Time (ms)"
1042 else:
1049 else:
1043 x = self.data.xrange[2]
1050 x = self.data.xrange[2]
1044 self.xlabel = "Velocity (m/s)"
1051 self.xlabel = "Velocity (m/s)"
1045
1052
1046 self.titles = []
1053 self.titles = []
1047
1054
1048 y = self.data.yrange
1055 y = self.data.yrange
1049 z = self.data[-1]['spc']
1056 z = self.data[-1]['spc']
1050 #print(z.shape)
1057 #print(z.shape)
1051 if len(self.height_index) > 0:
1058 if len(self.height_index) > 0:
1052 index = self.height_index
1059 index = self.height_index
1053 else:
1060 else:
1054 index = numpy.arange(0, len(y), int((len(y))/9))
1061 index = numpy.arange(0, len(y), int((len(y))/9))
1055 #print("inde x ", index, self.axes)
1062 #print("inde x ", index, self.axes)
1056
1063
1057 for n, ax in enumerate(self.axes):
1064 for n, ax in enumerate(self.axes):
1058
1065
1059 if ax.firsttime:
1066 if ax.firsttime:
1060
1067
1061
1068
1062 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1069 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1063 self.xmin = self.xmin if self.xmin else -self.xmax
1070 self.xmin = self.xmin if self.xmin else -self.xmax
1064 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
1071 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
1065 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
1072 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
1066
1073
1067
1074
1068 ax.plt = ax.plot(x, z[n, :, index].T)
1075 ax.plt = ax.plot(x, z[n, :, index].T)
1069 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
1076 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
1070 self.figures[0].legend(ax.plt, labels, loc='center right', prop={'size': 8})
1077 self.figures[0].legend(ax.plt, labels, loc='center right', prop={'size': 8})
1071 ax.minorticks_on()
1078 ax.minorticks_on()
1072 ax.grid(which='major', axis='both')
1079 ax.grid(which='major', axis='both')
1073 ax.grid(which='minor', axis='x')
1080 ax.grid(which='minor', axis='x')
1074 else:
1081 else:
1075 for i, line in enumerate(ax.plt):
1082 for i, line in enumerate(ax.plt):
1076 line.set_data(x, z[n, :, index[i]])
1083 line.set_data(x, z[n, :, index[i]])
1077
1084
1078
1085
1079 self.titles.append('CH {}'.format(self.channelList[n]))
1086 self.titles.append('CH {}'.format(self.channelList[n]))
1080 plt.suptitle(self.maintitle, fontsize=10)
1087 plt.suptitle(self.maintitle, fontsize=10)
1081
1088
1082
1089
1083 class BeaconPhase(Plot):
1090 class BeaconPhase(Plot):
1084
1091
1085 __isConfig = None
1092 __isConfig = None
1086 __nsubplots = None
1093 __nsubplots = None
1087
1094
1088 PREFIX = 'beacon_phase'
1095 PREFIX = 'beacon_phase'
1089
1096
1090 def __init__(self):
1097 def __init__(self):
1091 Plot.__init__(self)
1098 Plot.__init__(self)
1092 self.timerange = 24*60*60
1099 self.timerange = 24*60*60
1093 self.isConfig = False
1100 self.isConfig = False
1094 self.__nsubplots = 1
1101 self.__nsubplots = 1
1095 self.counter_imagwr = 0
1102 self.counter_imagwr = 0
1096 self.WIDTH = 800
1103 self.WIDTH = 800
1097 self.HEIGHT = 400
1104 self.HEIGHT = 400
1098 self.WIDTHPROF = 120
1105 self.WIDTHPROF = 120
1099 self.HEIGHTPROF = 0
1106 self.HEIGHTPROF = 0
1100 self.xdata = None
1107 self.xdata = None
1101 self.ydata = None
1108 self.ydata = None
1102
1109
1103 self.PLOT_CODE = BEACON_CODE
1110 self.PLOT_CODE = BEACON_CODE
1104
1111
1105 self.FTP_WEI = None
1112 self.FTP_WEI = None
1106 self.EXP_CODE = None
1113 self.EXP_CODE = None
1107 self.SUB_EXP_CODE = None
1114 self.SUB_EXP_CODE = None
1108 self.PLOT_POS = None
1115 self.PLOT_POS = None
1109
1116
1110 self.filename_phase = None
1117 self.filename_phase = None
1111
1118
1112 self.figfile = None
1119 self.figfile = None
1113
1120
1114 self.xmin = None
1121 self.xmin = None
1115 self.xmax = None
1122 self.xmax = None
1116
1123
1117 def getSubplots(self):
1124 def getSubplots(self):
1118
1125
1119 ncol = 1
1126 ncol = 1
1120 nrow = 1
1127 nrow = 1
1121
1128
1122 return nrow, ncol
1129 return nrow, ncol
1123
1130
1124 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1131 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1125
1132
1126 self.__showprofile = showprofile
1133 self.__showprofile = showprofile
1127 self.nplots = nplots
1134 self.nplots = nplots
1128
1135
1129 ncolspan = 7
1136 ncolspan = 7
1130 colspan = 6
1137 colspan = 6
1131 self.__nsubplots = 2
1138 self.__nsubplots = 2
1132
1139
1133 self.createFigure(id = id,
1140 self.createFigure(id = id,
1134 wintitle = wintitle,
1141 wintitle = wintitle,
1135 widthplot = self.WIDTH+self.WIDTHPROF,
1142 widthplot = self.WIDTH+self.WIDTHPROF,
1136 heightplot = self.HEIGHT+self.HEIGHTPROF,
1143 heightplot = self.HEIGHT+self.HEIGHTPROF,
1137 show=show)
1144 show=show)
1138
1145
1139 nrow, ncol = self.getSubplots()
1146 nrow, ncol = self.getSubplots()
1140
1147
1141 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1148 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1142
1149
1143 def save_phase(self, filename_phase):
1150 def save_phase(self, filename_phase):
1144 f = open(filename_phase,'w+')
1151 f = open(filename_phase,'w+')
1145 f.write('\n\n')
1152 f.write('\n\n')
1146 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1153 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1147 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1154 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1148 f.close()
1155 f.close()
1149
1156
1150 def save_data(self, filename_phase, data, data_datetime):
1157 def save_data(self, filename_phase, data, data_datetime):
1151 f=open(filename_phase,'a')
1158 f=open(filename_phase,'a')
1152 timetuple_data = data_datetime.timetuple()
1159 timetuple_data = data_datetime.timetuple()
1153 day = str(timetuple_data.tm_mday)
1160 day = str(timetuple_data.tm_mday)
1154 month = str(timetuple_data.tm_mon)
1161 month = str(timetuple_data.tm_mon)
1155 year = str(timetuple_data.tm_year)
1162 year = str(timetuple_data.tm_year)
1156 hour = str(timetuple_data.tm_hour)
1163 hour = str(timetuple_data.tm_hour)
1157 minute = str(timetuple_data.tm_min)
1164 minute = str(timetuple_data.tm_min)
1158 second = str(timetuple_data.tm_sec)
1165 second = str(timetuple_data.tm_sec)
1159 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1166 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1160 f.close()
1167 f.close()
1161
1168
1162 def plot(self):
1169 def plot(self):
1163 log.warning('TODO: Not yet implemented...')
1170 log.warning('TODO: Not yet implemented...')
1164
1171
1165 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1172 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1166 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1173 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1167 timerange=None,
1174 timerange=None,
1168 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1175 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1169 server=None, folder=None, username=None, password=None,
1176 server=None, folder=None, username=None, password=None,
1170 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1177 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1171
1178
1172 if dataOut.flagNoData:
1179 if dataOut.flagNoData:
1173 return dataOut
1180 return dataOut
1174
1181
1175 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1182 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1176 return
1183 return
1177
1184
1178 if pairsList == None:
1185 if pairsList == None:
1179 pairsIndexList = dataOut.pairsIndexList[:10]
1186 pairsIndexList = dataOut.pairsIndexList[:10]
1180 else:
1187 else:
1181 pairsIndexList = []
1188 pairsIndexList = []
1182 for pair in pairsList:
1189 for pair in pairsList:
1183 if pair not in dataOut.pairsList:
1190 if pair not in dataOut.pairsList:
1184 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
1191 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
1185 pairsIndexList.append(dataOut.pairsList.index(pair))
1192 pairsIndexList.append(dataOut.pairsList.index(pair))
1186
1193
1187 if pairsIndexList == []:
1194 if pairsIndexList == []:
1188 return
1195 return
1189
1196
1190 # if len(pairsIndexList) > 4:
1197 # if len(pairsIndexList) > 4:
1191 # pairsIndexList = pairsIndexList[0:4]
1198 # pairsIndexList = pairsIndexList[0:4]
1192
1199
1193 hmin_index = None
1200 hmin_index = None
1194 hmax_index = None
1201 hmax_index = None
1195
1202
1196 if hmin != None and hmax != None:
1203 if hmin != None and hmax != None:
1197 indexes = numpy.arange(dataOut.nHeights)
1204 indexes = numpy.arange(dataOut.nHeights)
1198 hmin_list = indexes[dataOut.heightList >= hmin]
1205 hmin_list = indexes[dataOut.heightList >= hmin]
1199 hmax_list = indexes[dataOut.heightList <= hmax]
1206 hmax_list = indexes[dataOut.heightList <= hmax]
1200
1207
1201 if hmin_list.any():
1208 if hmin_list.any():
1202 hmin_index = hmin_list[0]
1209 hmin_index = hmin_list[0]
1203
1210
1204 if hmax_list.any():
1211 if hmax_list.any():
1205 hmax_index = hmax_list[-1]+1
1212 hmax_index = hmax_list[-1]+1
1206
1213
1207 x = dataOut.getTimeRange()
1214 x = dataOut.getTimeRange()
1208
1215
1209 thisDatetime = dataOut.datatime
1216 thisDatetime = dataOut.datatime
1210
1217
1211 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1218 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1212 xlabel = "Local Time"
1219 xlabel = "Local Time"
1213 ylabel = "Phase (degrees)"
1220 ylabel = "Phase (degrees)"
1214
1221
1215 update_figfile = False
1222 update_figfile = False
1216
1223
1217 nplots = len(pairsIndexList)
1224 nplots = len(pairsIndexList)
1218 phase_beacon = numpy.zeros(len(pairsIndexList))
1225 phase_beacon = numpy.zeros(len(pairsIndexList))
1219 for i in range(nplots):
1226 for i in range(nplots):
1220 pair = dataOut.pairsList[pairsIndexList[i]]
1227 pair = dataOut.pairsList[pairsIndexList[i]]
1221 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1228 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1222 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1229 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1223 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1230 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1224 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1231 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1225 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1232 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1226
1233
1227 if dataOut.beacon_heiIndexList:
1234 if dataOut.beacon_heiIndexList:
1228 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1235 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1229 else:
1236 else:
1230 phase_beacon[i] = numpy.average(phase)
1237 phase_beacon[i] = numpy.average(phase)
1231
1238
1232 if not self.isConfig:
1239 if not self.isConfig:
1233
1240
1234 nplots = len(pairsIndexList)
1241 nplots = len(pairsIndexList)
1235
1242
1236 self.setup(id=id,
1243 self.setup(id=id,
1237 nplots=nplots,
1244 nplots=nplots,
1238 wintitle=wintitle,
1245 wintitle=wintitle,
1239 showprofile=showprofile,
1246 showprofile=showprofile,
1240 show=show)
1247 show=show)
1241
1248
1242 if timerange != None:
1249 if timerange != None:
1243 self.timerange = timerange
1250 self.timerange = timerange
1244
1251
1245 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1252 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1246
1253
1247 if ymin == None: ymin = 0
1254 if ymin == None: ymin = 0
1248 if ymax == None: ymax = 360
1255 if ymax == None: ymax = 360
1249
1256
1250 self.FTP_WEI = ftp_wei
1257 self.FTP_WEI = ftp_wei
1251 self.EXP_CODE = exp_code
1258 self.EXP_CODE = exp_code
1252 self.SUB_EXP_CODE = sub_exp_code
1259 self.SUB_EXP_CODE = sub_exp_code
1253 self.PLOT_POS = plot_pos
1260 self.PLOT_POS = plot_pos
1254
1261
1255 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1262 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1256 self.isConfig = True
1263 self.isConfig = True
1257 self.figfile = figfile
1264 self.figfile = figfile
1258 self.xdata = numpy.array([])
1265 self.xdata = numpy.array([])
1259 self.ydata = numpy.array([])
1266 self.ydata = numpy.array([])
1260
1267
1261 update_figfile = True
1268 update_figfile = True
1262
1269
1263 #open file beacon phase
1270 #open file beacon phase
1264 path = '%s%03d' %(self.PREFIX, self.id)
1271 path = '%s%03d' %(self.PREFIX, self.id)
1265 beacon_file = os.path.join(path,'%s.txt'%self.name)
1272 beacon_file = os.path.join(path,'%s.txt'%self.name)
1266 self.filename_phase = os.path.join(figpath,beacon_file)
1273 self.filename_phase = os.path.join(figpath,beacon_file)
1267
1274
1268 self.setWinTitle(title)
1275 self.setWinTitle(title)
1269
1276
1270
1277
1271 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1278 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1272
1279
1273 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1280 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1274
1281
1275 axes = self.axesList[0]
1282 axes = self.axesList[0]
1276
1283
1277 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1284 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1278
1285
1279 if len(self.ydata)==0:
1286 if len(self.ydata)==0:
1280 self.ydata = phase_beacon.reshape(-1,1)
1287 self.ydata = phase_beacon.reshape(-1,1)
1281 else:
1288 else:
1282 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1289 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1283
1290
1284
1291
1285 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1292 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1286 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1293 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1287 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1294 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1288 XAxisAsTime=True, grid='both'
1295 XAxisAsTime=True, grid='both'
1289 )
1296 )
1290
1297
1291 self.draw()
1298 self.draw()
1292
1299
1293 if dataOut.ltctime >= self.xmax:
1300 if dataOut.ltctime >= self.xmax:
1294 self.counter_imagwr = wr_period
1301 self.counter_imagwr = wr_period
1295 self.isConfig = False
1302 self.isConfig = False
1296 update_figfile = True
1303 update_figfile = True
1297
1304
1298 self.save(figpath=figpath,
1305 self.save(figpath=figpath,
1299 figfile=figfile,
1306 figfile=figfile,
1300 save=save,
1307 save=save,
1301 ftp=ftp,
1308 ftp=ftp,
1302 wr_period=wr_period,
1309 wr_period=wr_period,
1303 thisDatetime=thisDatetime,
1310 thisDatetime=thisDatetime,
1304 update_figfile=update_figfile)
1311 update_figfile=update_figfile)
1305
1312
1306 return dataOut
1313 return dataOut
1307
1314
1308 #####################################
1315 #####################################
1309 class NoiselessSpectraPlot(Plot):
1316 class NoiselessSpectraPlot(Plot):
1310 '''
1317 '''
1311 Plot for Spectra data, subtracting
1318 Plot for Spectra data, subtracting
1312 the noise in all channels, using for
1319 the noise in all channels, using for
1313 amisr-14 data
1320 amisr-14 data
1314 '''
1321 '''
1315
1322
1316 CODE = 'noiseless_spc'
1323 CODE = 'noiseless_spc'
1317 colormap = 'jet'
1324 colormap = 'jet'
1318 plot_type = 'pcolor'
1325 plot_type = 'pcolor'
1319 buffering = False
1326 buffering = False
1320 channelList = []
1327 channelList = []
1321 last_noise = None
1328 last_noise = None
1322
1329
1323 def setup(self):
1330 def setup(self):
1324
1331
1325 self.nplots = len(self.data.channels)
1332 self.nplots = len(self.data.channels)
1326 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
1333 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
1327 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
1334 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
1328 self.height = 3.5 * self.nrows
1335 self.height = 3.5 * self.nrows
1329
1336
1330 self.cb_label = 'dB'
1337 self.cb_label = 'dB'
1331 if self.showprofile:
1338 if self.showprofile:
1332 self.width = 5.8 * self.ncols
1339 self.width = 5.8 * self.ncols
1333 else:
1340 else:
1334 self.width = 4.8* self.ncols
1341 self.width = 4.8* self.ncols
1335 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.92, 'bottom': 0.12})
1342 self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.92, 'bottom': 0.12})
1336
1343
1337 self.ylabel = 'Range [km]'
1344 self.ylabel = 'Range [km]'
1338
1345
1339
1346
1340 def update_list(self,dataOut):
1347 def update_list(self,dataOut):
1341 if len(self.channelList) == 0:
1348 if len(self.channelList) == 0:
1342 self.channelList = dataOut.channelList
1349 self.channelList = dataOut.channelList
1343
1350
1344 def update(self, dataOut):
1351 def update(self, dataOut):
1345
1352
1346 self.update_list(dataOut)
1353 self.update_list(dataOut)
1347 data = {}
1354 data = {}
1348 meta = {}
1355 meta = {}
1349
1356
1350 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1357 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1351 n0 = (dataOut.getNoise()/norm)
1358 n0 = (dataOut.getNoise()/norm)
1352 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
1359 noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights)
1353 noise = 10*numpy.log10(noise)
1360 noise = 10*numpy.log10(noise)
1354
1361
1355 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
1362 z = numpy.zeros((dataOut.nChannels, dataOut.nFFTPoints, dataOut.nHeights))
1356 for ch in range(dataOut.nChannels):
1363 for ch in range(dataOut.nChannels):
1357 if hasattr(dataOut.normFactor,'ndim'):
1364 if hasattr(dataOut.normFactor,'ndim'):
1358 if dataOut.normFactor.ndim > 1:
1365 if dataOut.normFactor.ndim > 1:
1359 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
1366 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch]))
1360 else:
1367 else:
1361 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1368 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1362 else:
1369 else:
1363 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1370 z[ch] = (numpy.divide(dataOut.data_spc[ch],dataOut.normFactor))
1364
1371
1365 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1372 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
1366 spc = 10*numpy.log10(z)
1373 spc = 10*numpy.log10(z)
1367
1374
1368
1375
1369 data['spc'] = spc - noise
1376 data['spc'] = spc - noise
1370 #print(spc.shape)
1377 #print(spc.shape)
1371 data['rti'] = spc.mean(axis=1)
1378 data['rti'] = spc.mean(axis=1)
1372 data['noise'] = noise
1379 data['noise'] = noise
1373
1380
1374
1381
1375
1382
1376 # data['noise'] = noise
1383 # data['noise'] = noise
1377 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
1384 meta['xrange'] = (dataOut.getFreqRange(EXTRA_POINTS)/1000., dataOut.getAcfRange(EXTRA_POINTS), dataOut.getVelRange(EXTRA_POINTS))
1378
1385
1379 return data, meta
1386 return data, meta
1380
1387
1381 def plot(self):
1388 def plot(self):
1382 if self.xaxis == "frequency":
1389 if self.xaxis == "frequency":
1383 x = self.data.xrange[0]
1390 x = self.data.xrange[0]
1384 self.xlabel = "Frequency (kHz)"
1391 self.xlabel = "Frequency (kHz)"
1385 elif self.xaxis == "time":
1392 elif self.xaxis == "time":
1386 x = self.data.xrange[1]
1393 x = self.data.xrange[1]
1387 self.xlabel = "Time (ms)"
1394 self.xlabel = "Time (ms)"
1388 else:
1395 else:
1389 x = self.data.xrange[2]
1396 x = self.data.xrange[2]
1390 self.xlabel = "Velocity (m/s)"
1397 self.xlabel = "Velocity (m/s)"
1391
1398
1392 self.titles = []
1399 self.titles = []
1393 y = self.data.yrange
1400 y = self.data.yrange
1394 self.y = y
1401 self.y = y
1395
1402
1396 data = self.data[-1]
1403 data = self.data[-1]
1397 z = data['spc']
1404 z = data['spc']
1398
1405
1399 for n, ax in enumerate(self.axes):
1406 for n, ax in enumerate(self.axes):
1400 #noise = data['noise'][n]
1407 #noise = data['noise'][n]
1401
1408
1402 if ax.firsttime:
1409 if ax.firsttime:
1403 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1410 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1404 self.xmin = self.xmin if self.xmin else -self.xmax
1411 self.xmin = self.xmin if self.xmin else -self.xmax
1405 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
1412 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
1406 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
1413 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
1407 ax.plt = ax.pcolormesh(x, y, z[n].T,
1414 ax.plt = ax.pcolormesh(x, y, z[n].T,
1408 vmin=self.zmin,
1415 vmin=self.zmin,
1409 vmax=self.zmax,
1416 vmax=self.zmax,
1410 cmap=plt.get_cmap(self.colormap)
1417 cmap=plt.get_cmap(self.colormap)
1411 )
1418 )
1412
1419
1413 if self.showprofile:
1420 if self.showprofile:
1414 ax.plt_profile = self.pf_axes[n].plot(
1421 ax.plt_profile = self.pf_axes[n].plot(
1415 data['rti'][n], y)[0]
1422 data['rti'][n], y)[0]
1416
1423
1417
1424
1418 else:
1425 else:
1419 ax.plt.set_array(z[n].T.ravel())
1426 ax.plt.set_array(z[n].T.ravel())
1420 if self.showprofile:
1427 if self.showprofile:
1421 ax.plt_profile.set_data(data['rti'][n], y)
1428 ax.plt_profile.set_data(data['rti'][n], y)
1422
1429
1423
1430
1424 self.titles.append('CH {}'.format(self.channelList[n]))
1431 self.titles.append('CH {}'.format(self.channelList[n]))
1425
1432
1426
1433
1427 class NoiselessRTIPlot(RTIPlot):
1434 class NoiselessRTIPlot(RTIPlot):
1428 '''
1435 '''
1429 Plot for RTI data
1436 Plot for RTI data
1430 '''
1437 '''
1431
1438
1432 CODE = 'noiseless_rti'
1439 CODE = 'noiseless_rti'
1433 colormap = 'jet'
1440 colormap = 'jet'
1434 plot_type = 'pcolorbuffer'
1441 plot_type = 'pcolorbuffer'
1435 titles = None
1442 titles = None
1436 channelList = []
1443 channelList = []
1437 elevationList = []
1444 elevationList = []
1438 azimuthList = []
1445 azimuthList = []
1439 last_noise = None
1446 last_noise = None
1440
1447
1441 def setup(self):
1448 def setup(self):
1442 self.xaxis = 'time'
1449 self.xaxis = 'time'
1443 self.ncols = 1
1450 self.ncols = 1
1444 #print("dataChannels ",self.data.channels)
1451 #print("dataChannels ",self.data.channels)
1445 self.nrows = len(self.data.channels)
1452 self.nrows = len(self.data.channels)
1446 self.nplots = len(self.data.channels)
1453 self.nplots = len(self.data.channels)
1447 self.ylabel = 'Range [km]'
1454 self.ylabel = 'Range [km]'
1448 #self.xlabel = 'Time'
1455 #self.xlabel = 'Time'
1449 self.cb_label = 'dB'
1456 self.cb_label = 'dB'
1450 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1457 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1451 self.titles = ['{} Channel {}'.format(
1458 self.titles = ['{} Channel {}'.format(
1452 self.CODE.upper(), x) for x in range(self.nplots)]
1459 self.CODE.upper(), x) for x in range(self.nplots)]
1453
1460
1454 def update_list(self,dataOut):
1461 def update_list(self,dataOut):
1455 if len(self.channelList) == 0:
1462 if len(self.channelList) == 0:
1456 self.channelList = dataOut.channelList
1463 self.channelList = dataOut.channelList
1457 if len(self.elevationList) == 0:
1464 if len(self.elevationList) == 0:
1458 self.elevationList = dataOut.elevationList
1465 self.elevationList = dataOut.elevationList
1459 if len(self.azimuthList) == 0:
1466 if len(self.azimuthList) == 0:
1460 self.azimuthList = dataOut.azimuthList
1467 self.azimuthList = dataOut.azimuthList
1461
1468
1462 def update(self, dataOut):
1469 def update(self, dataOut):
1463 if len(self.channelList) == 0:
1470 if len(self.channelList) == 0:
1464 self.update_list(dataOut)
1471 self.update_list(dataOut)
1465
1472
1466 data = {}
1473 data = {}
1467 meta = {}
1474 meta = {}
1468 #print(dataOut.max_nIncohInt, dataOut.nIncohInt)
1475 #print(dataOut.max_nIncohInt, dataOut.nIncohInt)
1469 #print(dataOut.windowOfFilter,dataOut.nCohInt,dataOut.nProfiles,dataOut.max_nIncohInt,dataOut.nIncohInt
1476 #print(dataOut.windowOfFilter,dataOut.nCohInt,dataOut.nProfiles,dataOut.max_nIncohInt,dataOut.nIncohInt
1470 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1477 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1471 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1478 n0 = 10*numpy.log10(dataOut.getNoise()/norm)
1472 data['noise'] = n0
1479 data['noise'] = n0
1473 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
1480 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights)
1474 noiseless_data = dataOut.getPower() - noise
1481 noiseless_data = dataOut.getPower() - noise
1475
1482
1476 #print("power, noise:", dataOut.getPower(), n0)
1483 #print("power, noise:", dataOut.getPower(), n0)
1477 #print(noise)
1484 #print(noise)
1478 #print(noiseless_data)
1485 #print(noiseless_data)
1479
1486
1480 data['noiseless_rti'] = noiseless_data
1487 data['noiseless_rti'] = noiseless_data
1481
1488
1482 return data, meta
1489 return data, meta
1483
1490
1484 def plot(self):
1491 def plot(self):
1485 from matplotlib import pyplot as plt
1492 from matplotlib import pyplot as plt
1486 self.x = self.data.times
1493 self.x = self.data.times
1487 self.y = self.data.yrange
1494 self.y = self.data.yrange
1488 self.z = self.data['noiseless_rti']
1495 self.z = self.data['noiseless_rti']
1489 self.z = numpy.array(self.z, dtype=float)
1496 self.z = numpy.array(self.z, dtype=float)
1490 self.z = numpy.ma.masked_invalid(self.z)
1497 self.z = numpy.ma.masked_invalid(self.z)
1491
1498
1492
1499
1493 try:
1500 try:
1494 if self.channelList != None:
1501 if self.channelList != None:
1495 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1502 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1496 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1503 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1497 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1504 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1498 else:
1505 else:
1499 self.titles = ['{} Channel {}'.format(
1506 self.titles = ['{} Channel {}'.format(
1500 self.CODE.upper(), x) for x in self.channelList]
1507 self.CODE.upper(), x) for x in self.channelList]
1501 except:
1508 except:
1502 if self.channelList.any() != None:
1509 if self.channelList.any() != None:
1503 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1510 if len(self.elevationList) > 0 and len(self.azimuthList) > 0:
1504 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1511 self.titles = ['{} Channel {} ({:2.1f} Elev, {:2.1f} Azth)'.format(
1505 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1512 self.CODE.upper(), x, self.elevationList[x], self.azimuthList[x]) for x in self.channelList]
1506 else:
1513 else:
1507 self.titles = ['{} Channel {}'.format(
1514 self.titles = ['{} Channel {}'.format(
1508 self.CODE.upper(), x) for x in self.channelList]
1515 self.CODE.upper(), x) for x in self.channelList]
1509
1516
1510
1517
1511 if self.decimation is None:
1518 if self.decimation is None:
1512 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1519 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1513 else:
1520 else:
1514 x, y, z = self.fill_gaps(*self.decimate())
1521 x, y, z = self.fill_gaps(*self.decimate())
1515
1522
1516 dummy_var = self.axes #Extrañamente esto actualiza el valor axes
1523 dummy_var = self.axes #Extrañamente esto actualiza el valor axes
1517 #print("plot shapes ", z.shape, x.shape, y.shape)
1524 #print("plot shapes ", z.shape, x.shape, y.shape)
1518 #print(self.axes)
1525 #print(self.axes)
1519 for n, ax in enumerate(self.axes):
1526 for n, ax in enumerate(self.axes):
1520
1527
1521
1528
1522 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
1529 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
1523 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
1530 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
1524 data = self.data[-1]
1531 data = self.data[-1]
1525 if ax.firsttime:
1532 if ax.firsttime:
1526 if (n+1) == len(self.channelList):
1533 if (n+1) == len(self.channelList):
1527 ax.set_xlabel('Time')
1534 ax.set_xlabel('Time')
1528 ax.plt = ax.pcolormesh(x, y, z[n].T,
1535 ax.plt = ax.pcolormesh(x, y, z[n].T,
1529 vmin=self.zmin,
1536 vmin=self.zmin,
1530 vmax=self.zmax,
1537 vmax=self.zmax,
1531 cmap=plt.get_cmap(self.colormap)
1538 cmap=plt.get_cmap(self.colormap)
1532 )
1539 )
1533 if self.showprofile:
1540 if self.showprofile:
1534 ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0]
1541 ax.plot_profile = self.pf_axes[n].plot(data['noiseless_rti'][n], self.y)[0]
1535
1542
1536 else:
1543 else:
1537 # ax.collections.remove(ax.collections[0]) # error while running
1544 # ax.collections.remove(ax.collections[0]) # error while running
1538 ax.plt = ax.pcolormesh(x, y, z[n].T,
1545 ax.plt = ax.pcolormesh(x, y, z[n].T,
1539 vmin=self.zmin,
1546 vmin=self.zmin,
1540 vmax=self.zmax,
1547 vmax=self.zmax,
1541 cmap=plt.get_cmap(self.colormap)
1548 cmap=plt.get_cmap(self.colormap)
1542 )
1549 )
1543 if self.showprofile:
1550 if self.showprofile:
1544 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
1551 ax.plot_profile.set_data(data['noiseless_rti'][n], self.y)
1545 # if "noise" in self.data:
1552 # if "noise" in self.data:
1546 # #ax.plot_noise.set_data(numpy.repeat(data['noise'][n], len(self.y)), self.y)
1553 # #ax.plot_noise.set_data(numpy.repeat(data['noise'][n], len(self.y)), self.y)
1547 # ax.plot_noise.set_data(data['noise'][n], self.y)
1554 # ax.plot_noise.set_data(data['noise'][n], self.y)
1548
1555
1549
1556
1550 class OutliersRTIPlot(Plot):
1557 class OutliersRTIPlot(Plot):
1551 '''
1558 '''
1552 Plot for data_xxxx object
1559 Plot for data_xxxx object
1553 '''
1560 '''
1554
1561
1555 CODE = 'outlier_rtc' # Range Time Counts
1562 CODE = 'outlier_rtc' # Range Time Counts
1556 colormap = 'cool'
1563 colormap = 'cool'
1557 plot_type = 'pcolorbuffer'
1564 plot_type = 'pcolorbuffer'
1558
1565
1559 def setup(self):
1566 def setup(self):
1560 self.xaxis = 'time'
1567 self.xaxis = 'time'
1561 self.ncols = 1
1568 self.ncols = 1
1562 self.nrows = self.data.shape('outlier_rtc')[0]
1569 self.nrows = self.data.shape('outlier_rtc')[0]
1563 self.nplots = self.nrows
1570 self.nplots = self.nrows
1564 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1571 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1565
1572
1566
1573
1567 if not self.xlabel:
1574 if not self.xlabel:
1568 self.xlabel = 'Time'
1575 self.xlabel = 'Time'
1569
1576
1570 self.ylabel = 'Height [km]'
1577 self.ylabel = 'Height [km]'
1571 if not self.titles:
1578 if not self.titles:
1572 self.titles = ['Outliers Ch:{}'.format(x) for x in range(self.nrows)]
1579 self.titles = ['Outliers Ch:{}'.format(x) for x in range(self.nrows)]
1573
1580
1574 def update(self, dataOut):
1581 def update(self, dataOut):
1575
1582
1576 data = {}
1583 data = {}
1577 data['outlier_rtc'] = dataOut.data_outlier
1584 data['outlier_rtc'] = dataOut.data_outlier
1578
1585
1579 meta = {}
1586 meta = {}
1580
1587
1581 return data, meta
1588 return data, meta
1582
1589
1583 def plot(self):
1590 def plot(self):
1584 # self.data.normalize_heights()
1591 # self.data.normalize_heights()
1585 self.x = self.data.times
1592 self.x = self.data.times
1586 self.y = self.data.yrange
1593 self.y = self.data.yrange
1587 self.z = self.data['outlier_rtc']
1594 self.z = self.data['outlier_rtc']
1588
1595
1589 #self.z = numpy.ma.masked_invalid(self.z)
1596 #self.z = numpy.ma.masked_invalid(self.z)
1590
1597
1591 if self.decimation is None:
1598 if self.decimation is None:
1592 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1599 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1593 else:
1600 else:
1594 x, y, z = self.fill_gaps(*self.decimate())
1601 x, y, z = self.fill_gaps(*self.decimate())
1595
1602
1596 for n, ax in enumerate(self.axes):
1603 for n, ax in enumerate(self.axes):
1597
1604
1598 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1605 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1599 self.z[n])
1606 self.z[n])
1600 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1607 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1601 self.z[n])
1608 self.z[n])
1602 data = self.data[-1]
1609 data = self.data[-1]
1603 if ax.firsttime:
1610 if ax.firsttime:
1604 if self.zlimits is not None:
1611 if self.zlimits is not None:
1605 self.zmin, self.zmax = self.zlimits[n]
1612 self.zmin, self.zmax = self.zlimits[n]
1606
1613
1607 ax.plt = ax.pcolormesh(x, y, z[n].T,
1614 ax.plt = ax.pcolormesh(x, y, z[n].T,
1608 vmin=self.zmin,
1615 vmin=self.zmin,
1609 vmax=self.zmax,
1616 vmax=self.zmax,
1610 cmap=self.cmaps[n]
1617 cmap=self.cmaps[n]
1611 )
1618 )
1612 if self.showprofile:
1619 if self.showprofile:
1613 ax.plot_profile = self.pf_axes[n].plot(data['outlier_rtc'][n], self.y)[0]
1620 ax.plot_profile = self.pf_axes[n].plot(data['outlier_rtc'][n], self.y)[0]
1614 self.pf_axes[n].set_xlabel('')
1621 self.pf_axes[n].set_xlabel('')
1615 else:
1622 else:
1616 if self.zlimits is not None:
1623 if self.zlimits is not None:
1617 self.zmin, self.zmax = self.zlimits[n]
1624 self.zmin, self.zmax = self.zlimits[n]
1618 # ax.collections.remove(ax.collections[0]) # error while running
1625 # ax.collections.remove(ax.collections[0]) # error while running
1619 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1626 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1620 vmin=self.zmin,
1627 vmin=self.zmin,
1621 vmax=self.zmax,
1628 vmax=self.zmax,
1622 cmap=self.cmaps[n]
1629 cmap=self.cmaps[n]
1623 )
1630 )
1624 if self.showprofile:
1631 if self.showprofile:
1625 ax.plot_profile.set_data(data['outlier_rtc'][n], self.y)
1632 ax.plot_profile.set_data(data['outlier_rtc'][n], self.y)
1626 self.pf_axes[n].set_xlabel('')
1633 self.pf_axes[n].set_xlabel('')
1627
1634
1628 class NIncohIntRTIPlot(Plot):
1635 class NIncohIntRTIPlot(Plot):
1629 '''
1636 '''
1630 Plot for data_xxxx object
1637 Plot for data_xxxx object
1631 '''
1638 '''
1632
1639
1633 CODE = 'integrations_rtc' # Range Time Counts
1640 CODE = 'integrations_rtc' # Range Time Counts
1634 colormap = 'BuGn'
1641 colormap = 'BuGn'
1635 plot_type = 'pcolorbuffer'
1642 plot_type = 'pcolorbuffer'
1636
1643
1637 def setup(self):
1644 def setup(self):
1638 self.xaxis = 'time'
1645 self.xaxis = 'time'
1639 self.ncols = 1
1646 self.ncols = 1
1640 self.nrows = self.data.shape('integrations_rtc')[0]
1647 self.nrows = self.data.shape('integrations_rtc')[0]
1641 self.nplots = self.nrows
1648 self.nplots = self.nrows
1642 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1649 self.plots_adjust.update({'hspace':0.8, 'left': 0.08, 'bottom': 0.2, 'right':0.94})
1643
1650
1644
1651
1645 if not self.xlabel:
1652 if not self.xlabel:
1646 self.xlabel = 'Time'
1653 self.xlabel = 'Time'
1647
1654
1648 self.ylabel = 'Height [km]'
1655 self.ylabel = 'Height [km]'
1649 if not self.titles:
1656 if not self.titles:
1650 self.titles = ['Integration Ch:{}'.format(x) for x in range(self.nrows)]
1657 self.titles = ['Integration Ch:{}'.format(x) for x in range(self.nrows)]
1651
1658
1652 def update(self, dataOut):
1659 def update(self, dataOut):
1653
1660
1654 data = {}
1661 data = {}
1655 data['integrations_rtc'] = dataOut.nIncohInt
1662 data['integrations_rtc'] = dataOut.nIncohInt
1656
1663
1657 meta = {}
1664 meta = {}
1658
1665
1659 return data, meta
1666 return data, meta
1660
1667
1661 def plot(self):
1668 def plot(self):
1662 # self.data.normalize_heights()
1669 # self.data.normalize_heights()
1663 self.x = self.data.times
1670 self.x = self.data.times
1664 self.y = self.data.yrange
1671 self.y = self.data.yrange
1665 self.z = self.data['integrations_rtc']
1672 self.z = self.data['integrations_rtc']
1666
1673
1667 #self.z = numpy.ma.masked_invalid(self.z)
1674 #self.z = numpy.ma.masked_invalid(self.z)
1668
1675
1669 if self.decimation is None:
1676 if self.decimation is None:
1670 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1677 x, y, z = self.fill_gaps(self.x, self.y, self.z)
1671 else:
1678 else:
1672 x, y, z = self.fill_gaps(*self.decimate())
1679 x, y, z = self.fill_gaps(*self.decimate())
1673
1680
1674 for n, ax in enumerate(self.axes):
1681 for n, ax in enumerate(self.axes):
1675
1682
1676 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1683 self.zmax = self.zmax if self.zmax is not None else numpy.max(
1677 self.z[n])
1684 self.z[n])
1678 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1685 self.zmin = self.zmin if self.zmin is not None else numpy.min(
1679 self.z[n])
1686 self.z[n])
1680 data = self.data[-1]
1687 data = self.data[-1]
1681 if ax.firsttime:
1688 if ax.firsttime:
1682 if self.zlimits is not None:
1689 if self.zlimits is not None:
1683 self.zmin, self.zmax = self.zlimits[n]
1690 self.zmin, self.zmax = self.zlimits[n]
1684
1691
1685 ax.plt = ax.pcolormesh(x, y, z[n].T,
1692 ax.plt = ax.pcolormesh(x, y, z[n].T,
1686 vmin=self.zmin,
1693 vmin=self.zmin,
1687 vmax=self.zmax,
1694 vmax=self.zmax,
1688 cmap=self.cmaps[n]
1695 cmap=self.cmaps[n]
1689 )
1696 )
1690 if self.showprofile:
1697 if self.showprofile:
1691 ax.plot_profile = self.pf_axes[n].plot(data['integrations_rtc'][n], self.y)[0]
1698 ax.plot_profile = self.pf_axes[n].plot(data['integrations_rtc'][n], self.y)[0]
1692 self.pf_axes[n].set_xlabel('')
1699 self.pf_axes[n].set_xlabel('')
1693 else:
1700 else:
1694 if self.zlimits is not None:
1701 if self.zlimits is not None:
1695 self.zmin, self.zmax = self.zlimits[n]
1702 self.zmin, self.zmax = self.zlimits[n]
1696 # ax.collections.remove(ax.collections[0]) # error while running
1703 # ax.collections.remove(ax.collections[0]) # error while running
1697 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1704 ax.plt = ax.pcolormesh(x, y, z[n].T ,
1698 vmin=self.zmin,
1705 vmin=self.zmin,
1699 vmax=self.zmax,
1706 vmax=self.zmax,
1700 cmap=self.cmaps[n]
1707 cmap=self.cmaps[n]
1701 )
1708 )
1702 if self.showprofile:
1709 if self.showprofile:
1703 ax.plot_profile.set_data(data['integrations_rtc'][n], self.y)
1710 ax.plot_profile.set_data(data['integrations_rtc'][n], self.y)
1704 self.pf_axes[n].set_xlabel('')
1711 self.pf_axes[n].set_xlabel('')
1705
1712
1706
1713
1707
1714
1708 class RTIMapPlot(Plot):
1715 class RTIMapPlot(Plot):
1709 '''
1716 '''
1710 Plot for RTI data
1717 Plot for RTI data
1711
1718
1712 Example:
1719 Example:
1713
1720
1714 controllerObj = Project()
1721 controllerObj = Project()
1715 controllerObj.setup(id = '11', name='eej_proc', description=desc)
1722 controllerObj.setup(id = '11', name='eej_proc', description=desc)
1716 ##.......................................................................................
1723 ##.......................................................................................
1717 ##.......................................................................................
1724 ##.......................................................................................
1718 readUnitConfObj = controllerObj.addReadUnit(datatype='AMISRReader', path=inPath, startDate='2023/05/24',endDate='2023/05/24',
1725 readUnitConfObj = controllerObj.addReadUnit(datatype='AMISRReader', path=inPath, startDate='2023/05/24',endDate='2023/05/24',
1719 startTime='12:00:00',endTime='12:45:59',walk=1,timezone='lt',margin_days=1,code = code,nCode = nCode,
1726 startTime='12:00:00',endTime='12:45:59',walk=1,timezone='lt',margin_days=1,code = code,nCode = nCode,
1720 nBaud = nBaud,nOsamp = nosamp,nChannels=nChannels,nFFT=NFFT,
1727 nBaud = nBaud,nOsamp = nosamp,nChannels=nChannels,nFFT=NFFT,
1721 syncronization=False,shiftChannels=0)
1728 syncronization=False,shiftChannels=0)
1722
1729
1723 volts_proc = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
1730 volts_proc = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
1724
1731
1725 opObj01 = volts_proc.addOperation(name='Decoder', optype='other')
1732 opObj01 = volts_proc.addOperation(name='Decoder', optype='other')
1726 opObj01.addParameter(name='code', value=code, format='floatlist')
1733 opObj01.addParameter(name='code', value=code, format='floatlist')
1727 opObj01.addParameter(name='nCode', value=1, format='int')
1734 opObj01.addParameter(name='nCode', value=1, format='int')
1728 opObj01.addParameter(name='nBaud', value=nBaud, format='int')
1735 opObj01.addParameter(name='nBaud', value=nBaud, format='int')
1729 opObj01.addParameter(name='osamp', value=nosamp, format='int')
1736 opObj01.addParameter(name='osamp', value=nosamp, format='int')
1730
1737
1731 opObj12 = volts_proc.addOperation(name='selectHeights', optype='self')
1738 opObj12 = volts_proc.addOperation(name='selectHeights', optype='self')
1732 opObj12.addParameter(name='minHei', value='90', format='float')
1739 opObj12.addParameter(name='minHei', value='90', format='float')
1733 opObj12.addParameter(name='maxHei', value='150', format='float')
1740 opObj12.addParameter(name='maxHei', value='150', format='float')
1734
1741
1735 proc_spc = controllerObj.addProcUnit(datatype='SpectraProc', inputId=volts_proc.getId())
1742 proc_spc = controllerObj.addProcUnit(datatype='SpectraProc', inputId=volts_proc.getId())
1736 proc_spc.addParameter(name='nFFTPoints', value='8', format='int')
1743 proc_spc.addParameter(name='nFFTPoints', value='8', format='int')
1737
1744
1738 opObj11 = proc_spc.addOperation(name='IncohInt', optype='other')
1745 opObj11 = proc_spc.addOperation(name='IncohInt', optype='other')
1739 opObj11.addParameter(name='n', value='1', format='int')
1746 opObj11.addParameter(name='n', value='1', format='int')
1740
1747
1741 beamMapFile = "/home/japaza/Documents/AMISR_sky_mapper/UMET_beamcodes.csv"
1748 beamMapFile = "/home/japaza/Documents/AMISR_sky_mapper/UMET_beamcodes.csv"
1742
1749
1743 opObj12 = proc_spc.addOperation(name='RTIMapPlot', optype='external')
1750 opObj12 = proc_spc.addOperation(name='RTIMapPlot', optype='external')
1744 opObj12.addParameter(name='selectedHeightsList', value='95, 100, 105, 110 ', format='int')
1751 opObj12.addParameter(name='selectedHeightsList', value='95, 100, 105, 110 ', format='int')
1745 opObj12.addParameter(name='bField', value='100', format='int')
1752 opObj12.addParameter(name='bField', value='100', format='int')
1746 opObj12.addParameter(name='filename', value=beamMapFile, format='str')
1753 opObj12.addParameter(name='filename', value=beamMapFile, format='str')
1747
1754
1748 '''
1755 '''
1749
1756
1750 CODE = 'rti_skymap'
1757 CODE = 'rti_skymap'
1751
1758
1752 plot_type = 'scatter'
1759 plot_type = 'scatter'
1753 titles = None
1760 titles = None
1754 colormap = 'jet'
1761 colormap = 'jet'
1755 channelList = []
1762 channelList = []
1756 elevationList = []
1763 elevationList = []
1757 azimuthList = []
1764 azimuthList = []
1758 last_noise = None
1765 last_noise = None
1759 flag_setIndex = False
1766 flag_setIndex = False
1760 heights = []
1767 heights = []
1761 dcosx = []
1768 dcosx = []
1762 dcosy = []
1769 dcosy = []
1763 fullDcosy = None
1770 fullDcosy = None
1764 fullDcosy = None
1771 fullDcosy = None
1765 hindex = []
1772 hindex = []
1766 mapFile = False
1773 mapFile = False
1767 ##### BField ####
1774 ##### BField ####
1768 flagBField = False
1775 flagBField = False
1769 dcosxB = []
1776 dcosxB = []
1770 dcosyB = []
1777 dcosyB = []
1771 Bmarker = ['+','*','D','x','s','>','o','^']
1778 Bmarker = ['+','*','D','x','s','>','o','^']
1772
1779
1773
1780
1774 def setup(self):
1781 def setup(self):
1775
1782
1776 self.xaxis = 'Range (Km)'
1783 self.xaxis = 'Range (Km)'
1777 if len(self.selectedHeightsList) > 0:
1784 if len(self.selectedHeightsList) > 0:
1778 self.nplots = len(self.selectedHeightsList)
1785 self.nplots = len(self.selectedHeightsList)
1779 else:
1786 else:
1780 self.nplots = 4
1787 self.nplots = 4
1781 self.ncols = int(numpy.ceil(self.nplots/2))
1788 self.ncols = int(numpy.ceil(self.nplots/2))
1782 self.nrows = int(numpy.ceil(self.nplots/self.ncols))
1789 self.nrows = int(numpy.ceil(self.nplots/self.ncols))
1783 self.ylabel = 'dcosy'
1790 self.ylabel = 'dcosy'
1784 self.xlabel = 'dcosx'
1791 self.xlabel = 'dcosx'
1785 self.colorbar = True
1792 self.colorbar = True
1786 self.width = 6 + 4.1*self.nrows
1793 self.width = 6 + 4.1*self.nrows
1787 self.height = 3 + 3.5*self.ncols
1794 self.height = 3 + 3.5*self.ncols
1788
1795
1789
1796
1790 if self.extFile!=None:
1797 if self.extFile!=None:
1791 try:
1798 try:
1792 pointings = numpy.genfromtxt(self.extFile, delimiter=',')
1799 pointings = numpy.genfromtxt(self.extFile, delimiter=',')
1793 full_azi = pointings[:,1]
1800 full_azi = pointings[:,1]
1794 full_elev = pointings[:,2]
1801 full_elev = pointings[:,2]
1795 self.fullDcosx = numpy.cos(numpy.radians(full_elev))*numpy.sin(numpy.radians(full_azi))
1802 self.fullDcosx = numpy.cos(numpy.radians(full_elev))*numpy.sin(numpy.radians(full_azi))
1796 self.fullDcosy = numpy.cos(numpy.radians(full_elev))*numpy.cos(numpy.radians(full_azi))
1803 self.fullDcosy = numpy.cos(numpy.radians(full_elev))*numpy.cos(numpy.radians(full_azi))
1797 mapFile = True
1804 mapFile = True
1798 except Exception as e:
1805 except Exception as e:
1799 self.extFile = None
1806 self.extFile = None
1800 print(e)
1807 print(e)
1801
1808
1802
1809
1803 def update_list(self,dataOut):
1810 def update_list(self,dataOut):
1804 if len(self.channelList) == 0:
1811 if len(self.channelList) == 0:
1805 self.channelList = dataOut.channelList
1812 self.channelList = dataOut.channelList
1806 if len(self.elevationList) == 0:
1813 if len(self.elevationList) == 0:
1807 self.elevationList = dataOut.elevationList
1814 self.elevationList = dataOut.elevationList
1808 if len(self.azimuthList) == 0:
1815 if len(self.azimuthList) == 0:
1809 self.azimuthList = dataOut.azimuthList
1816 self.azimuthList = dataOut.azimuthList
1810 a = numpy.radians(numpy.asarray(self.azimuthList))
1817 a = numpy.radians(numpy.asarray(self.azimuthList))
1811 e = numpy.radians(numpy.asarray(self.elevationList))
1818 e = numpy.radians(numpy.asarray(self.elevationList))
1812 self.heights = dataOut.heightList
1819 self.heights = dataOut.heightList
1813 self.dcosx = numpy.cos(e)*numpy.sin(a)
1820 self.dcosx = numpy.cos(e)*numpy.sin(a)
1814 self.dcosy = numpy.cos(e)*numpy.cos(a)
1821 self.dcosy = numpy.cos(e)*numpy.cos(a)
1815
1822
1816 if len(self.bFieldList)>0:
1823 if len(self.bFieldList)>0:
1817 datetObj = datetime.datetime.fromtimestamp(dataOut.utctime)
1824 datetObj = datetime.datetime.fromtimestamp(dataOut.utctime)
1818 doy = datetObj.timetuple().tm_yday
1825 doy = datetObj.timetuple().tm_yday
1819 year = datetObj.year
1826 year = datetObj.year
1820 # self.dcosxB, self.dcosyB
1827 # self.dcosxB, self.dcosyB
1821 ObjB = BField(year=year,doy=doy,site=2,heights=self.bFieldList)
1828 ObjB = BField(year=year,doy=doy,site=2,heights=self.bFieldList)
1822 [dcos, alpha, nlon, nlat] = ObjB.getBField()
1829 [dcos, alpha, nlon, nlat] = ObjB.getBField()
1823
1830
1824 alpha_location = numpy.zeros((nlon,2,len(self.bFieldList)))
1831 alpha_location = numpy.zeros((nlon,2,len(self.bFieldList)))
1825 for ih in range(len(self.bFieldList)):
1832 for ih in range(len(self.bFieldList)):
1826 alpha_location[:,0,ih] = dcos[:,0,ih,0]
1833 alpha_location[:,0,ih] = dcos[:,0,ih,0]
1827 for ilon in numpy.arange(nlon):
1834 for ilon in numpy.arange(nlon):
1828 myx = (alpha[ilon,:,ih])[::-1]
1835 myx = (alpha[ilon,:,ih])[::-1]
1829 myy = (dcos[ilon,:,ih,0])[::-1]
1836 myy = (dcos[ilon,:,ih,0])[::-1]
1830 tck = splrep(myx,myy,s=0)
1837 tck = splrep(myx,myy,s=0)
1831 mydcosx = splev(ObjB.alpha_i,tck,der=0)
1838 mydcosx = splev(ObjB.alpha_i,tck,der=0)
1832
1839
1833 myx = (alpha[ilon,:,ih])[::-1]
1840 myx = (alpha[ilon,:,ih])[::-1]
1834 myy = (dcos[ilon,:,ih,1])[::-1]
1841 myy = (dcos[ilon,:,ih,1])[::-1]
1835 tck = splrep(myx,myy,s=0)
1842 tck = splrep(myx,myy,s=0)
1836 mydcosy = splev(ObjB.alpha_i,tck,der=0)
1843 mydcosy = splev(ObjB.alpha_i,tck,der=0)
1837 alpha_location[ilon,:,ih] = numpy.array([mydcosx, mydcosy])
1844 alpha_location[ilon,:,ih] = numpy.array([mydcosx, mydcosy])
1838 self.dcosxB.append(alpha_location[:,0,ih])
1845 self.dcosxB.append(alpha_location[:,0,ih])
1839 self.dcosyB.append(alpha_location[:,1,ih])
1846 self.dcosyB.append(alpha_location[:,1,ih])
1840 self.flagBField = True
1847 self.flagBField = True
1841
1848
1842 if len(self.celestialList)>0:
1849 if len(self.celestialList)>0:
1843 #getBField(self.bFieldList, date)
1850 #getBField(self.bFieldList, date)
1844 #pass = kwargs.get('celestial', [])
1851 #pass = kwargs.get('celestial', [])
1845 pass
1852 pass
1846
1853
1847
1854
1848 def update(self, dataOut):
1855 def update(self, dataOut):
1849
1856
1850 if len(self.channelList) == 0:
1857 if len(self.channelList) == 0:
1851 self.update_list(dataOut)
1858 self.update_list(dataOut)
1852
1859
1853 if not self.flag_setIndex:
1860 if not self.flag_setIndex:
1854 if len(self.selectedHeightsList)>0:
1861 if len(self.selectedHeightsList)>0:
1855 for sel_height in self.selectedHeightsList:
1862 for sel_height in self.selectedHeightsList:
1856 index_list = numpy.where(self.heights >= sel_height)
1863 index_list = numpy.where(self.heights >= sel_height)
1857 index_list = index_list[0]
1864 index_list = index_list[0]
1858 self.hindex.append(index_list[0])
1865 self.hindex.append(index_list[0])
1859 self.flag_setIndex = True
1866 self.flag_setIndex = True
1860
1867
1861 data = {}
1868 data = {}
1862 meta = {}
1869 meta = {}
1863
1870
1864 data['rti_skymap'] = dataOut.getPower()
1871 data['rti_skymap'] = dataOut.getPower()
1865 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1872 norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter
1866 noise = 10*numpy.log10(dataOut.getNoise()/norm)
1873 noise = 10*numpy.log10(dataOut.getNoise()/norm)
1867 data['noise'] = noise
1874 data['noise'] = noise
1868
1875
1869 return data, meta
1876 return data, meta
1870
1877
1871 def plot(self):
1878 def plot(self):
1872
1879
1873 self.x = self.dcosx
1880 self.x = self.dcosx
1874 self.y = self.dcosy
1881 self.y = self.dcosy
1875 self.z = self.data[-1]['rti_skymap']
1882 self.z = self.data[-1]['rti_skymap']
1876 self.z = numpy.array(self.z, dtype=float)
1883 self.z = numpy.array(self.z, dtype=float)
1877
1884
1878 if len(self.hindex) > 0:
1885 if len(self.hindex) > 0:
1879 index = self.hindex
1886 index = self.hindex
1880 else:
1887 else:
1881 index = numpy.arange(0, len(self.heights), int((len(self.heights))/4.2))
1888 index = numpy.arange(0, len(self.heights), int((len(self.heights))/4.2))
1882
1889
1883 self.titles = ['Height {:.2f} km '.format(self.heights[i])+" " for i in index]
1890 self.titles = ['Height {:.2f} km '.format(self.heights[i])+" " for i in index]
1884 for n, ax in enumerate(self.axes):
1891 for n, ax in enumerate(self.axes):
1885
1892
1886 if ax.firsttime:
1893 if ax.firsttime:
1887
1894
1888 self.xmax = self.xmax if self.xmax else numpy.nanmax(self.x)
1895 self.xmax = self.xmax if self.xmax else numpy.nanmax(self.x)
1889 self.xmin = self.xmin if self.xmin else numpy.nanmin(self.x)
1896 self.xmin = self.xmin if self.xmin else numpy.nanmin(self.x)
1890 self.ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
1897 self.ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
1891 self.ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
1898 self.ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
1892 self.zmax = self.zmax if self.zmax else numpy.nanmax(self.z)
1899 self.zmax = self.zmax if self.zmax else numpy.nanmax(self.z)
1893 self.zmin = self.zmin if self.zmin else numpy.nanmin(self.z)
1900 self.zmin = self.zmin if self.zmin else numpy.nanmin(self.z)
1894
1901
1895 if self.extFile!=None:
1902 if self.extFile!=None:
1896 ax.scatter(self.fullDcosx, self.fullDcosy, marker="+", s=20)
1903 ax.scatter(self.fullDcosx, self.fullDcosy, marker="+", s=20)
1897
1904
1898 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1905 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1899 s=60, marker="s", vmax = self.zmax)
1906 s=60, marker="s", vmax = self.zmax)
1900
1907
1901
1908
1902 ax.minorticks_on()
1909 ax.minorticks_on()
1903 ax.grid(which='major', axis='both')
1910 ax.grid(which='major', axis='both')
1904 ax.grid(which='minor', axis='x')
1911 ax.grid(which='minor', axis='x')
1905
1912
1906 if self.flagBField :
1913 if self.flagBField :
1907
1914
1908 for ih in range(len(self.bFieldList)):
1915 for ih in range(len(self.bFieldList)):
1909 label = str(self.bFieldList[ih]) + ' km'
1916 label = str(self.bFieldList[ih]) + ' km'
1910 ax.plot(self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1917 ax.plot(self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1911 label=label, linestyle='--', ms=4.0,lw=0.5)
1918 label=label, linestyle='--', ms=4.0,lw=0.5)
1912 handles, labels = ax.get_legend_handles_labels()
1919 handles, labels = ax.get_legend_handles_labels()
1913 a = -0.05
1920 a = -0.05
1914 b = 1.15 - 1.19*(self.nrows)
1921 b = 1.15 - 1.19*(self.nrows)
1915 self.axes[0].legend(handles,labels, bbox_to_anchor=(a,b), prop={'size': (5.8+ 1.1*self.nplots)}, title='B Field ⊥')
1922 self.axes[0].legend(handles,labels, bbox_to_anchor=(a,b), prop={'size': (5.8+ 1.1*self.nplots)}, title='B Field ⊥')
1916
1923
1917 else:
1924 else:
1918
1925
1919 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1926 ax.plt = ax.scatter(self.x, self.y, c=self.z[:,index[n]], cmap = 'jet',vmin = self.zmin,
1920 s=80, marker="s", vmax = self.zmax)
1927 s=80, marker="s", vmax = self.zmax)
1921
1928
1922 if self.flagBField :
1929 if self.flagBField :
1923 for ih in range(len(self.bFieldList)):
1930 for ih in range(len(self.bFieldList)):
1924 ax.plot (self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1931 ax.plot (self.dcosxB[ih], self.dcosyB[ih], color='k', marker=self.Bmarker[ih % 8],
1925 linestyle='--', ms=4.0,lw=0.5)
1932 linestyle='--', ms=4.0,lw=0.5)
1926
1933
1927
1934
1928
1935
Show file before | Show file after | Hide comments
@@ -1,815 +1,819
1 import os
1 import os
2 import time
2 import time
3 import datetime
3 import datetime
4
4
5 import numpy
5 import numpy
6 import h5py
6 import h5py
7
7
8 import schainpy.admin
8 import schainpy.admin
9 from schainpy.model.data.jrodata import *
9 from schainpy.model.data.jrodata import *
10 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
10 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
11 from schainpy.model.io.jroIO_base import *
11 from schainpy.model.io.jroIO_base import *
12 from schainpy.utils import log
12 from schainpy.utils import log
13
13
14
14
15 class HDFReader(Reader, ProcessingUnit):
15 class HDFReader(Reader, ProcessingUnit):
16 """Processing unit to read HDF5 format files
16 """Processing unit to read HDF5 format files
17
17
18 This unit reads HDF5 files created with `HDFWriter` operation contains
18 This unit reads HDF5 files created with `HDFWriter` operation contains
19 by default two groups Data and Metadata all variables would be saved as `dataOut`
19 by default two groups Data and Metadata all variables would be saved as `dataOut`
20 attributes.
20 attributes.
21 It is possible to read any HDF5 file by given the structure in the `description`
21 It is possible to read any HDF5 file by given the structure in the `description`
22 parameter, also you can add extra values to metadata with the parameter `extras`.
22 parameter, also you can add extra values to metadata with the parameter `extras`.
23
23
24 Parameters:
24 Parameters:
25 -----------
25 -----------
26 path : str
26 path : str
27 Path where files are located.
27 Path where files are located.
28 startDate : date
28 startDate : date
29 Start date of the files
29 Start date of the files
30 endDate : list
30 endDate : list
31 End date of the files
31 End date of the files
32 startTime : time
32 startTime : time
33 Start time of the files
33 Start time of the files
34 endTime : time
34 endTime : time
35 End time of the files
35 End time of the files
36 description : dict, optional
36 description : dict, optional
37 Dictionary with the description of the HDF5 file
37 Dictionary with the description of the HDF5 file
38 extras : dict, optional
38 extras : dict, optional
39 Dictionary with extra metadata to be be added to `dataOut`
39 Dictionary with extra metadata to be be added to `dataOut`
40
40
41 Attention: Be carefull, add attribute utcoffset, in the last part of reader in order to work in Local Time without time problems.
41 Attention: Be carefull, add attribute utcoffset, in the last part of reader in order to work in Local Time without time problems.
42
42
43 -----------
43 -----------
44 utcoffset='-18000'
44 utcoffset='-18000'
45
45
46
46
47 Examples
47 Examples
48 --------
48 --------
49
49
50 desc = {
50 desc = {
51 'Data': {
51 'Data': {
52 'data_output': ['u', 'v', 'w'],
52 'data_output': ['u', 'v', 'w'],
53 'utctime': 'timestamps',
53 'utctime': 'timestamps',
54 } ,
54 } ,
55 'Metadata': {
55 'Metadata': {
56 'heightList': 'heights'
56 'heightList': 'heights'
57 }
57 }
58 }
58 }
59
59
60 desc = {
60 desc = {
61 'Data': {
61 'Data': {
62 'data_output': 'winds',
62 'data_output': 'winds',
63 'utctime': 'timestamps'
63 'utctime': 'timestamps'
64 },
64 },
65 'Metadata': {
65 'Metadata': {
66 'heightList': 'heights'
66 'heightList': 'heights'
67 }
67 }
68 }
68 }
69
69
70 extras = {
70 extras = {
71 'timeZone': 300
71 'timeZone': 300
72 }
72 }
73
73
74 reader = project.addReadUnit(
74 reader = project.addReadUnit(
75 name='HDFReader',
75 name='HDFReader',
76 path='/path/to/files',
76 path='/path/to/files',
77 startDate='2019/01/01',
77 startDate='2019/01/01',
78 endDate='2019/01/31',
78 endDate='2019/01/31',
79 startTime='00:00:00',
79 startTime='00:00:00',
80 endTime='23:59:59',
80 endTime='23:59:59',
81 utcoffset='-18000'
81 utcoffset='-18000'
82 # description=json.dumps(desc),
82 # description=json.dumps(desc),
83 # extras=json.dumps(extras),
83 # extras=json.dumps(extras),
84 )
84 )
85
85
86 """
86 """
87
87
88 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'description', 'extras']
88 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'description', 'extras']
89
89
90 def __init__(self):
90 def __init__(self):
91
91
92 ProcessingUnit.__init__(self)
92 ProcessingUnit.__init__(self)
93 self.ext = ".hdf5"
93 self.ext = ".hdf5"
94 self.optchar = "D"
94 self.optchar = "D"
95 self.meta = {}
95 self.meta = {}
96 self.data = {}
96 self.data = {}
97 self.open_file = h5py.File
97 self.open_file = h5py.File
98 self.open_mode = 'r'
98 self.open_mode = 'r'
99 self.description = {}
99 self.description = {}
100 self.extras = {}
100 self.extras = {}
101 self.filefmt = "*%Y%j***"
101 self.filefmt = "*%Y%j***"
102 self.folderfmt = "*%Y%j"
102 self.folderfmt = "*%Y%j"
103 self.utcoffset = 0
103 self.utcoffset = 0
104 self.flagUpdateDataOut = False
104 self.flagUpdateDataOut = False
105 self.dataOut = Parameters()
105 self.dataOut = Parameters()
106 self.dataOut.error=False ## NOTE: Importante definir esto antes inicio
106 self.dataOut.error=False ## NOTE: Importante definir esto antes inicio
107 self.dataOut.flagNoData = True
107 self.dataOut.flagNoData = True
108
108
109 def setup(self, **kwargs):
109 def setup(self, **kwargs):
110
110
111 self.set_kwargs(**kwargs)
111 self.set_kwargs(**kwargs)
112 if not self.ext.startswith('.'):
112 if not self.ext.startswith('.'):
113 self.ext = '.{}'.format(self.ext)
113 self.ext = '.{}'.format(self.ext)
114
114
115 if self.online:
115 if self.online:
116 log.log("Searching files in online mode...", self.name)
116 log.log("Searching files in online mode...", self.name)
117
117
118 for nTries in range(self.nTries):
118 for nTries in range(self.nTries):
119 fullpath = self.searchFilesOnLine(self.path, self.startDate,
119 fullpath = self.searchFilesOnLine(self.path, self.startDate,
120 self.endDate, self.expLabel, self.ext, self.walk,
120 self.endDate, self.expLabel, self.ext, self.walk,
121 self.filefmt, self.folderfmt)
121 self.filefmt, self.folderfmt)
122 pathname, filename = os.path.split(fullpath)
122 pathname, filename = os.path.split(fullpath)
123 try:
123 try:
124 fullpath = next(fullpath)
124 fullpath = next(fullpath)
125 except:
125 except:
126 fullpath = None
126 fullpath = None
127
127
128 if fullpath:
128 if fullpath:
129 break
129 break
130
130
131 log.warning(
131 log.warning(
132 'Waiting {} sec for a valid file in {}: try {} ...'.format(
132 'Waiting {} sec for a valid file in {}: try {} ...'.format(
133 self.delay, self.path, nTries + 1),
133 self.delay, self.path, nTries + 1),
134 self.name)
134 self.name)
135 time.sleep(self.delay)
135 time.sleep(self.delay)
136
136
137 if not(fullpath):
137 if not(fullpath):
138 raise schainpy.admin.SchainError(
138 raise schainpy.admin.SchainError(
139 'There isn\'t any valid file in {}'.format(self.path))
139 'There isn\'t any valid file in {}'.format(self.path))
140
140
141 pathname, filename = os.path.split(fullpath)
141 pathname, filename = os.path.split(fullpath)
142 self.year = int(filename[1:5])
142 self.year = int(filename[1:5])
143 self.doy = int(filename[5:8])
143 self.doy = int(filename[5:8])
144 self.set = int(filename[8:11]) - 1
144 self.set = int(filename[8:11]) - 1
145 else:
145 else:
146 log.log("Searching files in {}".format(self.path), self.name)
146 log.log("Searching files in {}".format(self.path), self.name)
147 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
147 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
148 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
148 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
149
149
150 self.setNextFile()
150 self.setNextFile()
151
151
152 return
152 return
153
153
154 # def readFirstHeader(self):
154 # def readFirstHeader(self):
155 # '''Read metadata and data'''
155 # '''Read metadata and data'''
156
156
157 # self.__readMetadata()
157 # self.__readMetadata()
158 # self.__readData()
158 # self.__readData()
159 # self.__setBlockList()
159 # self.__setBlockList()
160
160
161 # if 'type' in self.meta:
161 # if 'type' in self.meta:
162 # self.dataOut = eval(self.meta['type'])()
162 # self.dataOut = eval(self.meta['type'])()
163
163
164 # for attr in self.meta:
164 # for attr in self.meta:
165 # setattr(self.dataOut, attr, self.meta[attr])
165 # setattr(self.dataOut, attr, self.meta[attr])
166
166
167 # self.blockIndex = 0
167 # self.blockIndex = 0
168
168
169 # return
169 # return
170
170
171 def readFirstHeader(self):
171 def readFirstHeader(self):
172 '''Read metadata and data'''
172 '''Read metadata and data'''
173
173
174 self.__readMetadata2()
174 self.__readMetadata2()
175 self.__readData()
175 self.__readData()
176 self.__setBlockList()
176 self.__setBlockList()
177 if 'type' in self.meta:
177 if 'type' in self.meta:
178 self.dataOut = eval(self.meta['type'])()
178 self.dataOut = eval(self.meta['type'])()
179
179
180 for attr in self.meta:
180 for attr in self.meta:
181 if "processingHeaderObj" in attr:
181 if "processingHeaderObj" in attr:
182 self.flagUpdateDataOut=True
182 self.flagUpdateDataOut=True
183 at = attr.split('.')
183 at = attr.split('.')
184 if len(at) > 1:
184 if len(at) > 1:
185 setattr(eval("self.dataOut."+at[0]),at[1], self.meta[attr])
185 setattr(eval("self.dataOut."+at[0]),at[1], self.meta[attr])
186 else:
186 else:
187 setattr(self.dataOut, attr, self.meta[attr])
187 setattr(self.dataOut, attr, self.meta[attr])
188 self.blockIndex = 0
188 self.blockIndex = 0
189
189
190 if self.flagUpdateDataOut:
190 if self.flagUpdateDataOut:
191 self.updateDataOut()
191 self.updateDataOut()
192
192
193 return
193 return
194
194
195 def updateDataOut(self):
195 def updateDataOut(self):
196
196
197 self.dataOut.azimuthList = self.dataOut.processingHeaderObj.azimuthList
197 self.dataOut.azimuthList = self.dataOut.processingHeaderObj.azimuthList
198 self.dataOut.elevationList = self.dataOut.processingHeaderObj.elevationList
198 self.dataOut.elevationList = self.dataOut.processingHeaderObj.elevationList
199 self.dataOut.heightList = self.dataOut.processingHeaderObj.heightList
199 self.dataOut.heightList = self.dataOut.processingHeaderObj.heightList
200 self.dataOut.ippSeconds = self.dataOut.processingHeaderObj.ipp
200 self.dataOut.ippSeconds = self.dataOut.processingHeaderObj.ipp
201 self.dataOut.elevationList = self.dataOut.processingHeaderObj.elevationList
201 self.dataOut.elevationList = self.dataOut.processingHeaderObj.elevationList
202 self.dataOut.channelList = self.dataOut.processingHeaderObj.channelList
202 self.dataOut.channelList = self.dataOut.processingHeaderObj.channelList
203 self.dataOut.nCohInt = self.dataOut.processingHeaderObj.nCohInt
203 self.dataOut.nCohInt = self.dataOut.processingHeaderObj.nCohInt
204 self.dataOut.nFFTPoints = self.dataOut.processingHeaderObj.nFFTPoints
204 self.dataOut.nFFTPoints = self.dataOut.processingHeaderObj.nFFTPoints
205 self.flagUpdateDataOut = False
205 self.flagUpdateDataOut = False
206 self.dataOut.frequency = self.dataOut.radarControllerHeaderObj.frequency
206 self.dataOut.frequency = self.dataOut.radarControllerHeaderObj.frequency
207 #self.dataOut.heightList = self.dataOut.processingHeaderObj.heightList
207 #self.dataOut.heightList = self.dataOut.processingHeaderObj.heightList
208
208
209 def __setBlockList(self):
209 def __setBlockList(self):
210 '''
210 '''
211 Selects the data within the times defined
211 Selects the data within the times defined
212
212
213 self.fp
213 self.fp
214 self.startTime
214 self.startTime
215 self.endTime
215 self.endTime
216 self.blockList
216 self.blockList
217 self.blocksPerFile
217 self.blocksPerFile
218
218
219 '''
219 '''
220
220
221 startTime = self.startTime
221 startTime = self.startTime
222 endTime = self.endTime
222 endTime = self.endTime
223 thisUtcTime = self.data['utctime'] + self.utcoffset
223 thisUtcTime = self.data['utctime'] + self.utcoffset
224 # self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1])
224 # self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1])
225 thisDatetime = datetime.datetime.utcfromtimestamp(thisUtcTime[0])
225 thisDatetime = datetime.datetime.utcfromtimestamp(thisUtcTime[0])
226 self.startFileDatetime = thisDatetime
226 self.startFileDatetime = thisDatetime
227 thisDate = thisDatetime.date()
227 thisDate = thisDatetime.date()
228 thisTime = thisDatetime.time()
228 thisTime = thisDatetime.time()
229 startUtcTime = (datetime.datetime.combine(thisDate, startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
229 startUtcTime = (datetime.datetime.combine(thisDate, startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
230 endUtcTime = (datetime.datetime.combine(thisDate, endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
230 endUtcTime = (datetime.datetime.combine(thisDate, endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
231 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
231 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
232
232
233 self.blockList = ind
233 self.blockList = ind
234 self.blocksPerFile = len(ind)
234 self.blocksPerFile = len(ind)
235 # self.blocksPerFile = len(thisUtcTime)
235 # self.blocksPerFile = len(thisUtcTime)
236 if len(ind)==0:
236 if len(ind)==0:
237 print("[Reading] Block No. %d/%d -> %s [Skipping]" % (self.blockIndex,
237 print("[Reading] Block No. %d/%d -> %s [Skipping]" % (self.blockIndex,
238 self.blocksPerFile,
238 self.blocksPerFile,
239 thisDatetime))
239 thisDatetime))
240 self.setNextFile()
240 self.setNextFile()
241
241
242 return
242 return
243
243
244 def __readMetadata(self):
244 def __readMetadata(self):
245 '''
245 '''
246 Reads Metadata
246 Reads Metadata
247 '''
247 '''
248
248
249 meta = {}
249 meta = {}
250
250
251 if self.description:
251 if self.description:
252 for key, value in self.description['Metadata'].items():
252 for key, value in self.description['Metadata'].items():
253 meta[key] = self.fp[value][()]
253 meta[key] = self.fp[value][()]
254 else:
254 else:
255 grp = self.fp['Metadata']
255 grp = self.fp['Metadata']
256 for name in grp:
256 for name in grp:
257 meta[name] = grp[name][()]
257 meta[name] = grp[name][()]
258
258
259 if self.extras:
259 if self.extras:
260 for key, value in self.extras.items():
260 for key, value in self.extras.items():
261 meta[key] = value
261 meta[key] = value
262 self.meta = meta
262 self.meta = meta
263
263
264 return
264 return
265
265
266 def __readMetadata2(self):
266 def __readMetadata2(self):
267 '''
267 '''
268 Reads Metadata
268 Reads Metadata
269 '''
269 '''
270 meta = {}
270 meta = {}
271 if self.description:
271 if self.description:
272 for key, value in self.description['Metadata'].items():
272 for key, value in self.description['Metadata'].items():
273 meta[key] = self.fp[value][()]
273 meta[key] = self.fp[value][()]
274 else:
274 else:
275 grp = self.fp['Metadata']
275 grp = self.fp['Metadata']
276 for item in grp.values():
276 for item in grp.values():
277 name = item.name
277 name = item.name
278 if isinstance(item, h5py.Dataset):
278 if isinstance(item, h5py.Dataset):
279 name = name.split("/")[-1]
279 name = name.split("/")[-1]
280 meta[name] = item[()]
280 meta[name] = item[()]
281 else:
281 else:
282 grp2 = self.fp[name]
282 grp2 = self.fp[name]
283 Obj = name.split("/")[-1]
283 Obj = name.split("/")[-1]
284
284
285 for item2 in grp2.values():
285 for item2 in grp2.values():
286 name2 = Obj+"."+item2.name.split("/")[-1]
286 name2 = Obj+"."+item2.name.split("/")[-1]
287 meta[name2] = item2[()]
287 meta[name2] = item2[()]
288
288
289 if self.extras:
289 if self.extras:
290 for key, value in self.extras.items():
290 for key, value in self.extras.items():
291 meta[key] = value
291 meta[key] = value
292 self.meta = meta
292 self.meta = meta
293
293
294 return
294 return
295
295
296 def __readData(self):
296 def __readData(self):
297
297
298 data = {}
298 data = {}
299
299
300 if self.description:
300 if self.description:
301 for key, value in self.description['Data'].items():
301 for key, value in self.description['Data'].items():
302 if isinstance(value, str):
302 if isinstance(value, str):
303 if isinstance(self.fp[value], h5py.Dataset):
303 if isinstance(self.fp[value], h5py.Dataset):
304 data[key] = self.fp[value][()]
304 data[key] = self.fp[value][()]
305 elif isinstance(self.fp[value], h5py.Group):
305 elif isinstance(self.fp[value], h5py.Group):
306 array = []
306 array = []
307 for ch in self.fp[value]:
307 for ch in self.fp[value]:
308 array.append(self.fp[value][ch][()])
308 array.append(self.fp[value][ch][()])
309 data[key] = numpy.array(array)
309 data[key] = numpy.array(array)
310 elif isinstance(value, list):
310 elif isinstance(value, list):
311 array = []
311 array = []
312 for ch in value:
312 for ch in value:
313 array.append(self.fp[ch][()])
313 array.append(self.fp[ch][()])
314 data[key] = numpy.array(array)
314 data[key] = numpy.array(array)
315 else:
315 else:
316 grp = self.fp['Data']
316 grp = self.fp['Data']
317 for name in grp:
317 for name in grp:
318 if isinstance(grp[name], h5py.Dataset):
318 if isinstance(grp[name], h5py.Dataset):
319 array = grp[name][()]
319 array = grp[name][()]
320 elif isinstance(grp[name], h5py.Group):
320 elif isinstance(grp[name], h5py.Group):
321 array = []
321 array = []
322 for ch in grp[name]:
322 for ch in grp[name]:
323 array.append(grp[name][ch][()])
323 array.append(grp[name][ch][()])
324 array = numpy.array(array)
324 array = numpy.array(array)
325 else:
325 else:
326 log.warning('Unknown type: {}'.format(name))
326 log.warning('Unknown type: {}'.format(name))
327
327
328 if name in self.description:
328 if name in self.description:
329 key = self.description[name]
329 key = self.description[name]
330 else:
330 else:
331 key = name
331 key = name
332 data[key] = array
332 data[key] = array
333
333
334 self.data = data
334 self.data = data
335 return
335 return
336
336
337 def getData(self):
337 def getData(self):
338
338
339 if not self.isDateTimeInRange(self.startFileDatetime, self.startDate, self.endDate, self.startTime, self.endTime):
339 if not self.isDateTimeInRange(self.startFileDatetime, self.startDate, self.endDate, self.startTime, self.endTime):
340 self.dataOut.flagNoData = True
340 self.dataOut.flagNoData = True
341 self.blockIndex = self.blocksPerFile
341 self.blockIndex = self.blocksPerFile
342 self.dataOut.error = True # TERMINA EL PROGRAMA
342 self.dataOut.error = True # TERMINA EL PROGRAMA
343 return
343 return
344 for attr in self.data:
344 for attr in self.data:
345
345
346 if self.data[attr].ndim == 1:
346 if self.data[attr].ndim == 1:
347 setattr(self.dataOut, attr, self.data[attr][self.blockIndex])
347 setattr(self.dataOut, attr, self.data[attr][self.blockIndex])
348 else:
348 else:
349 setattr(self.dataOut, attr, self.data[attr][:, self.blockIndex])
349 setattr(self.dataOut, attr, self.data[attr][:, self.blockIndex])
350
350
351
351
352 self.blockIndex += 1
352 self.blockIndex += 1
353
353
354 if self.blockIndex == 1:
354 if self.blockIndex == 1:
355 log.log("Block No. {}/{} -> {}".format(
355 log.log("Block No. {}/{} -> {}".format(
356 self.blockIndex,
356 self.blockIndex,
357 self.blocksPerFile,
357 self.blocksPerFile,
358 self.dataOut.datatime.ctime()), self.name)
358 self.dataOut.datatime.ctime()), self.name)
359 else:
359 else:
360 log.log("Block No. {}/{} ".format(
360 log.log("Block No. {}/{} ".format(
361 self.blockIndex,
361 self.blockIndex,
362 self.blocksPerFile),self.name)
362 self.blocksPerFile),self.name)
363
363
364 if self.blockIndex == self.blocksPerFile:
364 if self.blockIndex == self.blocksPerFile:
365 self.setNextFile()
365 self.setNextFile()
366
366
367 self.dataOut.flagNoData = False
367 self.dataOut.flagNoData = False
368
368
369 return
369 return
370
370
371 def run(self, **kwargs):
371 def run(self, **kwargs):
372
372
373 if not(self.isConfig):
373 if not(self.isConfig):
374 self.setup(**kwargs)
374 self.setup(**kwargs)
375 self.isConfig = True
375 self.isConfig = True
376
376
377 if self.blockIndex == self.blocksPerFile:
377 if self.blockIndex == self.blocksPerFile:
378 self.setNextFile()
378 self.setNextFile()
379
379
380 self.getData()
380 self.getData()
381
381
382 return
382 return
383
383
384 @MPDecorator
384 @MPDecorator
385 class HDFWriter(Operation):
385 class HDFWriter(Operation):
386 """Operation to write HDF5 files.
386 """Operation to write HDF5 files.
387
387
388 The HDF5 file contains by default two groups Data and Metadata where
388 The HDF5 file contains by default two groups Data and Metadata where
389 you can save any `dataOut` attribute specified by `dataList` and `metadataList`
389 you can save any `dataOut` attribute specified by `dataList` and `metadataList`
390 parameters, data attributes are normaly time dependent where the metadata
390 parameters, data attributes are normaly time dependent where the metadata
391 are not.
391 are not.
392 It is possible to customize the structure of the HDF5 file with the
392 It is possible to customize the structure of the HDF5 file with the
393 optional description parameter see the examples.
393 optional description parameter see the examples.
394
394
395 Parameters:
395 Parameters:
396 -----------
396 -----------
397 path : str
397 path : str
398 Path where files will be saved.
398 Path where files will be saved.
399 blocksPerFile : int
399 blocksPerFile : int
400 Number of blocks per file
400 Number of blocks per file
401 metadataList : list
401 metadataList : list
402 List of the dataOut attributes that will be saved as metadata
402 List of the dataOut attributes that will be saved as metadata
403 dataList : int
403 dataList : int
404 List of the dataOut attributes that will be saved as data
404 List of the dataOut attributes that will be saved as data
405 setType : bool
405 setType : bool
406 If True the name of the files corresponds to the timestamp of the data
406 If True the name of the files corresponds to the timestamp of the data
407 description : dict, optional
407 description : dict, optional
408 Dictionary with the desired description of the HDF5 file
408 Dictionary with the desired description of the HDF5 file
409
409
410 Examples
410 Examples
411 --------
411 --------
412
412
413 desc = {
413 desc = {
414 'data_output': {'winds': ['z', 'w', 'v']},
414 'data_output': {'winds': ['z', 'w', 'v']},
415 'utctime': 'timestamps',
415 'utctime': 'timestamps',
416 'heightList': 'heights'
416 'heightList': 'heights'
417 }
417 }
418 desc = {
418 desc = {
419 'data_output': ['z', 'w', 'v'],
419 'data_output': ['z', 'w', 'v'],
420 'utctime': 'timestamps',
420 'utctime': 'timestamps',
421 'heightList': 'heights'
421 'heightList': 'heights'
422 }
422 }
423 desc = {
423 desc = {
424 'Data': {
424 'Data': {
425 'data_output': 'winds',
425 'data_output': 'winds',
426 'utctime': 'timestamps'
426 'utctime': 'timestamps'
427 },
427 },
428 'Metadata': {
428 'Metadata': {
429 'heightList': 'heights'
429 'heightList': 'heights'
430 }
430 }
431 }
431 }
432
432
433 writer = proc_unit.addOperation(name='HDFWriter')
433 writer = proc_unit.addOperation(name='HDFWriter')
434 writer.addParameter(name='path', value='/path/to/file')
434 writer.addParameter(name='path', value='/path/to/file')
435 writer.addParameter(name='blocksPerFile', value='32')
435 writer.addParameter(name='blocksPerFile', value='32')
436 writer.addParameter(name='metadataList', value='heightList,timeZone')
436 writer.addParameter(name='metadataList', value='heightList,timeZone')
437 writer.addParameter(name='dataList',value='data_output,utctime')
437 writer.addParameter(name='dataList',value='data_output,utctime')
438 # writer.addParameter(name='description',value=json.dumps(desc))
438 # writer.addParameter(name='description',value=json.dumps(desc))
439
439
440 """
440 """
441
441
442 ext = ".hdf5"
442 ext = ".hdf5"
443 optchar = "D"
443 optchar = "D"
444 filename = None
444 filename = None
445 path = None
445 path = None
446 setFile = None
446 setFile = None
447 fp = None
447 fp = None
448 ds = None
448 firsttime = True
449 firsttime = True
449 #Configurations
450 #Configurations
450 blocksPerFile = None
451 blocksPerFile = None
451 blockIndex = None
452 blockIndex = None
452 dataOut = None #eval ??????
453 dataOut = None #eval ??????
453 #Data Arrays
454 #Data Arrays
454 dataList = None
455 dataList = None
455 metadataList = None
456 metadataList = None
456 currentDay = None
457 currentDay = None
457 lastTime = None
458 lastTime = None
458 timeZone = "ut"
459 timeZone = "ut"
459 hourLimit = 3
460 hourLimit = 3
460 breakDays = True
461 breakDays = True
461
462
462 def __init__(self):
463 def __init__(self):
463
464
464 Operation.__init__(self)
465 Operation.__init__(self)
465 return
466 return
466
467
467 def set_kwargs(self, **kwargs):
468 def set_kwargs(self, **kwargs):
468
469
469 for key, value in kwargs.items():
470 for key, value in kwargs.items():
470 setattr(self, key, value)
471 setattr(self, key, value)
471
472
472 def set_kwargs_obj(self, obj, **kwargs):
473 def set_kwargs_obj(self, obj, **kwargs):
473
474
474 for key, value in kwargs.items():
475 for key, value in kwargs.items():
475 setattr(obj, key, value)
476 setattr(obj, key, value)
476
477
477 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None,
478 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None,
478 description={},timeZone = "ut",hourLimit = 3, breakDays=True, **kwargs):
479 description={},timeZone = "ut",hourLimit = 3, breakDays=True, **kwargs):
479 self.path = path
480 self.path = path
480 self.blocksPerFile = blocksPerFile
481 self.blocksPerFile = blocksPerFile
481 self.metadataList = metadataList
482 self.metadataList = metadataList
482 self.dataList = [s.strip() for s in dataList]
483 self.dataList = [s.strip() for s in dataList]
483 self.setType = setType
484 self.setType = setType
484 self.description = description
485 self.description = description
485 self.timeZone = timeZone
486 self.timeZone = timeZone
486 self.hourLimit = hourLimit
487 self.hourLimit = hourLimit
487 self.breakDays = breakDays
488 self.breakDays = breakDays
488 self.set_kwargs(**kwargs)
489 self.set_kwargs(**kwargs)
489
490
490 if self.metadataList is None:
491 if self.metadataList is None:
491 self.metadataList = self.dataOut.metadata_list
492 self.metadataList = self.dataOut.metadata_list
492
493
494 self.metadataList = list(set(self.metadataList))
495
493 tableList = []
496 tableList = []
494 dsList = []
497 dsList = []
495
498
496 for i in range(len(self.dataList)):
499 for i in range(len(self.dataList)):
497 dsDict = {}
500 dsDict = {}
498 if hasattr(self.dataOut, self.dataList[i]):
501 if hasattr(self.dataOut, self.dataList[i]):
499 dataAux = getattr(self.dataOut, self.dataList[i])
502 dataAux = getattr(self.dataOut, self.dataList[i])
500 dsDict['variable'] = self.dataList[i]
503 dsDict['variable'] = self.dataList[i]
501 else:
504 else:
502 log.warning('Attribute {} not found in dataOut'.format(self.dataList[i]),self.name)
505 log.warning('Attribute {} not found in dataOut'.format(self.dataList[i]),self.name)
503 continue
506 continue
504
507
505 if dataAux is None:
508 if dataAux is None:
506 continue
509 continue
507 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
510 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float_)):
508 dsDict['nDim'] = 0
511 dsDict['nDim'] = 0
509 else:
512 else:
510 dsDict['nDim'] = len(dataAux.shape)
513 dsDict['nDim'] = len(dataAux.shape)
511 dsDict['shape'] = dataAux.shape
514 dsDict['shape'] = dataAux.shape
512 dsDict['dsNumber'] = dataAux.shape[0]
515 dsDict['dsNumber'] = dataAux.shape[0]
513 dsDict['dtype'] = dataAux.dtype
516 dsDict['dtype'] = dataAux.dtype
514
517
515 dsList.append(dsDict)
518 dsList.append(dsDict)
516
519
520 self.blockIndex = 0
517 self.dsList = dsList
521 self.dsList = dsList
518 self.currentDay = self.dataOut.datatime.date()
522 self.currentDay = self.dataOut.datatime.date()
519
523
520 def timeFlag(self):
524 def timeFlag(self):
521 currentTime = self.dataOut.utctime
525 currentTime = self.dataOut.utctime
522 timeTuple = None
526 timeTuple = None
523 if self.timeZone == "lt":
527 if self.timeZone == "lt":
524 timeTuple = time.localtime(currentTime)
528 timeTuple = time.localtime(currentTime)
525 else :
529 else :
526 timeTuple = time.gmtime(currentTime)
530 timeTuple = time.gmtime(currentTime)
527 dataDay = timeTuple.tm_yday
531 dataDay = timeTuple.tm_yday
528
532
529 if self.lastTime is None:
533 if self.lastTime is None:
530 self.lastTime = currentTime
534 self.lastTime = currentTime
531 self.currentDay = dataDay
535 self.currentDay = dataDay
532 return False
536 return False
533
537
534 timeDiff = currentTime - self.lastTime
538 timeDiff = currentTime - self.lastTime
535
539
536 # Si el dia es diferente o si la diferencia entre un
540 # Si el dia es diferente o si la diferencia entre un
537 # dato y otro supera self.hourLimit
541 # dato y otro supera self.hourLimit
538 if (dataDay != self.currentDay) and self.breakDays:
542 if (dataDay != self.currentDay) and self.breakDays:
539 self.currentDay = dataDay
543 self.currentDay = dataDay
540 return True
544 return True
541 elif timeDiff > self.hourLimit*60*60:
545 elif timeDiff > self.hourLimit*60*60:
542 self.lastTime = currentTime
546 self.lastTime = currentTime
543 return True
547 return True
544 else:
548 else:
545 self.lastTime = currentTime
549 self.lastTime = currentTime
546 return False
550 return False
547
551
548 def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
552 def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
549 dataList=[], setType=None, description={}, **kwargs):
553 dataList=[], setType=None, description={}, **kwargs):
550
554
551 self.dataOut = dataOut
555 self.dataOut = dataOut
552 self.set_kwargs_obj(self.dataOut, **kwargs)
556 self.set_kwargs_obj(self.dataOut, **kwargs)
553 if not(self.isConfig):
557 if not(self.isConfig):
554 self.setup(path=path, blocksPerFile=blocksPerFile,
558 self.setup(path=path, blocksPerFile=blocksPerFile,
555 metadataList=metadataList, dataList=dataList,
559 metadataList=metadataList, dataList=dataList,
556 setType=setType, description=description, **kwargs)
560 setType=setType, description=description, **kwargs)
557
561
558 self.isConfig = True
562 self.isConfig = True
559 self.setNextFile()
563 self.setNextFile()
560
564
561 self.putData()
565 self.putData()
562 return
566 return
563
567
564 def setNextFile(self):
568 def setNextFile(self):
565
569
566 ext = self.ext
570 ext = self.ext
567 path = self.path
571 path = self.path
568 setFile = self.setFile
572 setFile = self.setFile
569 timeTuple = None
573 timeTuple = None
570 if self.timeZone == "lt":
574 if self.timeZone == "lt":
571 timeTuple = time.localtime(self.dataOut.utctime)
575 timeTuple = time.localtime(self.dataOut.utctime)
572 elif self.timeZone == "ut":
576 elif self.timeZone == "ut":
573 timeTuple = time.gmtime(self.dataOut.utctime)
577 timeTuple = time.gmtime(self.dataOut.utctime)
574 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
578 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
575 fullpath = os.path.join(path, subfolder)
579 fullpath = os.path.join(path, subfolder)
576
580
577 if os.path.exists(fullpath):
581 if os.path.exists(fullpath):
578 filesList = os.listdir(fullpath)
582 filesList = os.listdir(fullpath)
579 filesList = [k for k in filesList if k.startswith(self.optchar)]
583 filesList = [k for k in filesList if k.startswith(self.optchar)]
580 if len(filesList) > 0:
584 if len(filesList) > 0:
581 filesList = sorted(filesList, key=str.lower)
585 filesList = sorted(filesList, key=str.lower)
582 filen = filesList[-1]
586 filen = filesList[-1]
583 # el filename debera tener el siguiente formato
587 # el filename debera tener el siguiente formato
584 # 0 1234 567 89A BCDE (hex)
588 # 0 1234 567 89A BCDE (hex)
585 # x YYYY DDD SSS .ext
589 # x YYYY DDD SSS .ext
586 if isNumber(filen[8:11]):
590 if isNumber(filen[8:11]):
587 setFile = int(filen[8:11]) #inicializo mi contador de seteo al seteo del ultimo file
591 setFile = int(filen[8:11]) #inicializo mi contador de seteo al seteo del ultimo file
588 else:
592 else:
589 setFile = -1
593 setFile = -1
590 else:
594 else:
591 setFile = -1 #inicializo mi contador de seteo
595 setFile = -1 #inicializo mi contador de seteo
592 else:
596 else:
593 os.makedirs(fullpath)
597 os.makedirs(fullpath)
594 setFile = -1 #inicializo mi contador de seteo
598 setFile = -1 #inicializo mi contador de seteo
595
599
596 if self.setType is None:
600 if self.setType is None:
597 setFile += 1
601 setFile += 1
598 file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
602 file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
599 timeTuple.tm_year,
603 timeTuple.tm_year,
600 timeTuple.tm_yday,
604 timeTuple.tm_yday,
601 setFile,
605 setFile,
602 ext)
606 ext)
603 else:
607 else:
604 setFile = timeTuple.tm_hour*60+timeTuple.tm_min
608 setFile = timeTuple.tm_hour*60+timeTuple.tm_min
605 file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
609 file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
606 timeTuple.tm_year,
610 timeTuple.tm_year,
607 timeTuple.tm_yday,
611 timeTuple.tm_yday,
608 setFile,
612 setFile,
609 ext)
613 ext)
610
614
611 self.filename = os.path.join(path, subfolder, file)
615 self.filename = os.path.join(path, subfolder, file)
612
616
613
617
614
618
615 def getLabel(self, name, x=None):
619 def getLabel(self, name, x=None):
616
620
617 if x is None:
621 if x is None:
618 if 'Data' in self.description:
622 if 'Data' in self.description:
619 data = self.description['Data']
623 data = self.description['Data']
620 if 'Metadata' in self.description:
624 if 'Metadata' in self.description:
621 data.update(self.description['Metadata'])
625 data.update(self.description['Metadata'])
622 else:
626 else:
623 data = self.description
627 data = self.description
624 if name in data:
628 if name in data:
625 if isinstance(data[name], str):
629 if isinstance(data[name], str):
626 return data[name]
630 return data[name]
627 elif isinstance(data[name], list):
631 elif isinstance(data[name], list):
628 return None
632 return None
629 elif isinstance(data[name], dict):
633 elif isinstance(data[name], dict):
630 for key, value in data[name].items():
634 for key, value in data[name].items():
631 return key
635 return key
632 return name
636 return name
633 else:
637 else:
634 if 'Metadata' in self.description:
638 if 'Metadata' in self.description:
635 meta = self.description['Metadata']
639 meta = self.description['Metadata']
636 else:
640 else:
637 meta = self.description
641 meta = self.description
638 if name in meta:
642 if name in meta:
639 if isinstance(meta[name], list):
643 if isinstance(meta[name], list):
640 return meta[name][x]
644 return meta[name][x]
641 elif isinstance(meta[name], dict):
645 elif isinstance(meta[name], dict):
642 for key, value in meta[name].items():
646 for key, value in meta[name].items():
643 return value[x]
647 return value[x]
644 if 'cspc' in name:
648 if 'cspc' in name:
645 return 'pair{:02d}'.format(x)
649 return 'pair{:02d}'.format(x)
646 else:
650 else:
647 return 'channel{:02d}'.format(x)
651 return 'channel{:02d}'.format(x)
648
652
649 def writeMetadata(self, fp):
653 def writeMetadata(self, fp):
650
654
651 if self.description:
655 if self.description:
652 if 'Metadata' in self.description:
656 if 'Metadata' in self.description:
653 grp = fp.create_group('Metadata')
657 grp = fp.create_group('Metadata')
654 else:
658 else:
655 grp = fp
659 grp = fp
656 else:
660 else:
657 grp = fp.create_group('Metadata')
661 grp = fp.create_group('Metadata')
658
662
659 for i in range(len(self.metadataList)):
663 for i in range(len(self.metadataList)):
660 if not hasattr(self.dataOut, self.metadataList[i]):
664 if not hasattr(self.dataOut, self.metadataList[i]):
661 log.warning('Metadata: `{}` not found'.format(self.metadataList[i]), self.name)
665 log.warning('Metadata: `{}` not found'.format(self.metadataList[i]), self.name)
662 continue
666 continue
663 value = getattr(self.dataOut, self.metadataList[i])
667 value = getattr(self.dataOut, self.metadataList[i])
664 if isinstance(value, bool):
668 if isinstance(value, bool):
665 if value is True:
669 if value is True:
666 value = 1
670 value = 1
667 else:
671 else:
668 value = 0
672 value = 0
669 grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
673 grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
670 return
674 return
671
675
672 def writeMetadata2(self, fp):
676 def writeMetadata2(self, fp):
673
677
674 if self.description:
678 if self.description:
675 if 'Metadata' in self.description:
679 if 'Metadata' in self.description:
676 grp = fp.create_group('Metadata')
680 grp = fp.create_group('Metadata')
677 else:
681 else:
678 grp = fp
682 grp = fp
679 else:
683 else:
680 grp = fp.create_group('Metadata')
684 grp = fp.create_group('Metadata')
681
685
682 for i in range(len(self.metadataList)):
686 for i in range(len(self.metadataList)):
683
687
684 attribute = self.metadataList[i]
688 attribute = self.metadataList[i]
685 attr = attribute.split('.')
689 attr = attribute.split('.')
686 if len(attr) > 1:
690 if len(attr) > 1:
687 if not hasattr(eval("self.dataOut."+attr[0]),attr[1]):
691 if not hasattr(eval("self.dataOut."+attr[0]),attr[1]):
688 log.warning('Metadata: {}.{} not found'.format(attr[0],attr[1]), self.name)
692 log.warning('Metadata: {}.{} not found'.format(attr[0],attr[1]), self.name)
689 continue
693 continue
690 value = getattr(eval("self.dataOut."+attr[0]),attr[1])
694 value = getattr(eval("self.dataOut."+attr[0]),attr[1])
691 if isinstance(value, bool):
695 if isinstance(value, bool):
692 if value is True:
696 if value is True:
693 value = 1
697 value = 1
694 else:
698 else:
695 value = 0
699 value = 0
696 if isinstance(value,type(None)):
700 if isinstance(value,type(None)):
697 log.warning("Invalid value detected, {} is None".format(attribute), self.name)
701 log.warning("Invalid value detected, {} is None".format(attribute), self.name)
698 value = 0
702 value = 0
699 grp2 = None
703 grp2 = None
700 if not 'Metadata/'+attr[0] in fp:
704 if not 'Metadata/'+attr[0] in fp:
701 grp2 = fp.create_group('Metadata/'+attr[0])
705 grp2 = fp.create_group('Metadata/'+attr[0])
702 else:
706 else:
703 grp2 = fp['Metadata/'+attr[0]]
707 grp2 = fp['Metadata/'+attr[0]]
704 grp2.create_dataset(attr[1], data=value)
708 grp2.create_dataset(attr[1], data=value)
705
709
706 else:
710 else:
707 if not hasattr(self.dataOut, attr[0] ):
711 if not hasattr(self.dataOut, attr[0] ):
708 log.warning('Metadata: `{}` not found'.format(attribute), self.name)
712 log.warning('Metadata: `{}` not found'.format(attribute), self.name)
709 continue
713 continue
710 value = getattr(self.dataOut, attr[0])
714 value = getattr(self.dataOut, attr[0])
711 if isinstance(value, bool):
715 if isinstance(value, bool):
712 if value is True:
716 if value is True:
713 value = 1
717 value = 1
714 else:
718 else:
715 value = 0
719 value = 0
716 if isinstance(value, type(None)):
720 if isinstance(value, type(None)):
717 log.error("Value {} is None".format(attribute),self.name)
721 log.error("Value {} is None".format(attribute),self.name)
718
722
719 grp.create_dataset(self.getLabel(attribute), data=value)
723 grp.create_dataset(self.getLabel(attribute), data=value)
720
724
721 return
725 return
722
726
723 def writeData(self, fp):
727 def writeData(self, fp):
724
728
725 if self.description:
729 if self.description:
726 if 'Data' in self.description:
730 if 'Data' in self.description:
727 grp = fp.create_group('Data')
731 grp = fp.create_group('Data')
728 else:
732 else:
729 grp = fp
733 grp = fp
730 else:
734 else:
731 grp = fp.create_group('Data')
735 grp = fp.create_group('Data')
732
736
733 dtsets = []
737 dtsets = []
734 data = []
738 data = []
735
739
736 for dsInfo in self.dsList:
740 for dsInfo in self.dsList:
737 if dsInfo['nDim'] == 0:
741 if dsInfo['nDim'] == 0:
738 ds = grp.create_dataset(
742 ds = grp.create_dataset(
739 self.getLabel(dsInfo['variable']),
743 self.getLabel(dsInfo['variable']),
740 (self.blocksPerFile,),
744 (self.blocksPerFile,),
741 chunks=True,
745 chunks=True,
742 dtype=numpy.float64)
746 dtype=numpy.float64)
743 dtsets.append(ds)
747 dtsets.append(ds)
744 data.append((dsInfo['variable'], -1))
748 data.append((dsInfo['variable'], -1))
745 else:
749 else:
746 label = self.getLabel(dsInfo['variable'])
750 label = self.getLabel(dsInfo['variable'])
747 if label is not None:
751 if label is not None:
748 sgrp = grp.create_group(label)
752 sgrp = grp.create_group(label)
749 else:
753 else:
750 sgrp = grp
754 sgrp = grp
751 for i in range(dsInfo['dsNumber']):
755 for i in range(dsInfo['dsNumber']):
752 ds = sgrp.create_dataset(
756 ds = sgrp.create_dataset(
753 self.getLabel(dsInfo['variable'], i),
757 self.getLabel(dsInfo['variable'], i),
754 (self.blocksPerFile,) + dsInfo['shape'][1:],
758 (self.blocksPerFile,) + dsInfo['shape'][1:],
755 chunks=True,
759 chunks=True,
756 dtype=dsInfo['dtype'])
760 dtype=dsInfo['dtype'])
757 dtsets.append(ds)
761 dtsets.append(ds)
758 data.append((dsInfo['variable'], i))
762 data.append((dsInfo['variable'], i))
759 fp.flush()
763 fp.flush()
760
764
761 log.log('Creating file: {}'.format(fp.filename), self.name)
765 log.log('Creating file: {}'.format(fp.filename), self.name)
762
766
763 self.ds = dtsets
767 self.ds = dtsets
764 self.data = data
768 self.data = data
765 self.firsttime = True
769 self.firsttime = True
766 self.blockIndex = 0
770
767 return
771 return
768
772
769 def putData(self):
773 def putData(self):
770
774
771 if (self.blockIndex == self.blocksPerFile) or self.timeFlag():
775 if (self.blockIndex == self.blocksPerFile) or self.timeFlag():
772 self.closeFile()
776 self.closeFile()
773 self.setNextFile()
777 self.setNextFile()
774 self.dataOut.flagNoData = False
778 self.dataOut.flagNoData = False
775 self.blockIndex = 0
779 self.blockIndex = 0
776
780
777 if self.blockIndex == 0:
781 if self.blockIndex == 0:
778 #Setting HDF5 File
782 #Setting HDF5 File
779 self.fp = h5py.File(self.filename, 'w')
783 self.fp = h5py.File(self.filename, 'w')
780 #write metadata
784 #write metadata
781 self.writeMetadata2(self.fp)
785 self.writeMetadata2(self.fp)
782 #Write data
786 #Write data
783 self.writeData(self.fp)
787 self.writeData(self.fp)
784 log.log('Block No. {}/{} --> {}'.format(self.blockIndex+1, self.blocksPerFile,self.dataOut.datatime.ctime()), self.name)
788 log.log('Block No. {}/{} --> {}'.format(self.blockIndex+1, self.blocksPerFile,self.dataOut.datatime.ctime()), self.name)
785 elif (self.blockIndex % 10 ==0):
789 elif (self.blockIndex % 10 ==0):
786 log.log('Block No. {}/{} --> {}'.format(self.blockIndex+1, self.blocksPerFile,self.dataOut.datatime.ctime()), self.name)
790 log.log('Block No. {}/{} --> {}'.format(self.blockIndex+1, self.blocksPerFile,self.dataOut.datatime.ctime()), self.name)
787 else:
791 else:
788
792
789 log.log('Block No. {}/{}'.format(self.blockIndex+1, self.blocksPerFile), self.name)
793 log.log('Block No. {}/{}'.format(self.blockIndex+1, self.blocksPerFile), self.name)
790
794
791 for i, ds in enumerate(self.ds):
795 for i, ds in enumerate(self.ds):
792 attr, ch = self.data[i]
796 attr, ch = self.data[i]
793 if ch == -1:
797 if ch == -1:
794 ds[self.blockIndex] = getattr(self.dataOut, attr)
798 ds[self.blockIndex] = getattr(self.dataOut, attr)
795 else:
799 else:
796 ds[self.blockIndex] = getattr(self.dataOut, attr)[ch]
800 ds[self.blockIndex] = getattr(self.dataOut, attr)[ch]
797
801
798 self.blockIndex += 1
802 self.blockIndex += 1
799
803
800 self.fp.flush()
804 self.fp.flush()
801 self.dataOut.flagNoData = True
805 self.dataOut.flagNoData = True
802
806
803 def closeFile(self):
807 def closeFile(self):
804
808
805 if self.blockIndex != self.blocksPerFile:
809 if self.blockIndex != self.blocksPerFile:
806 for ds in self.ds:
810 for ds in self.ds:
807 ds.resize(self.blockIndex, axis=0)
811 ds.resize(self.blockIndex, axis=0)
808
812
809 if self.fp:
813 if self.fp:
810 self.fp.flush()
814 self.fp.flush()
811 self.fp.close()
815 self.fp.close()
812
816
813 def close(self):
817 def close(self):
814
818
815 self.closeFile()
819 self.closeFile()
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
General Comments 0
You need to be logged in to leave comments. Login now