##// END OF EJS Templates
schain
RSS
Bug in jrodata ncode nbaud
jespinoza -
r1340:4f5903c4ea97
parent child
Show More
Side by Side Unified
Context file:
r1340:4f5903c4ea97 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1193 +1,1193
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
20 from .jroheaderIO import SystemHeader, RadarControllerHeader
21 from schainpy.model.data import _noise
21 from schainpy.model.data import _noise
22
22
23
23
24 def getNumpyDtype(dataTypeCode):
24 def getNumpyDtype(dataTypeCode):
25
25
26 if dataTypeCode == 0:
26 if dataTypeCode == 0:
27 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
27 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
28 elif dataTypeCode == 1:
28 elif dataTypeCode == 1:
29 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
29 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
30 elif dataTypeCode == 2:
30 elif dataTypeCode == 2:
31 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
31 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
32 elif dataTypeCode == 3:
32 elif dataTypeCode == 3:
33 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
33 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
34 elif dataTypeCode == 4:
34 elif dataTypeCode == 4:
35 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
35 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
36 elif dataTypeCode == 5:
36 elif dataTypeCode == 5:
37 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
37 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
38 else:
38 else:
39 raise ValueError('dataTypeCode was not defined')
39 raise ValueError('dataTypeCode was not defined')
40
40
41 return numpyDtype
41 return numpyDtype
42
42
43
43
44 def getDataTypeCode(numpyDtype):
44 def getDataTypeCode(numpyDtype):
45
45
46 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
46 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
47 datatype = 0
47 datatype = 0
48 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
48 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
49 datatype = 1
49 datatype = 1
50 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
50 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
51 datatype = 2
51 datatype = 2
52 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
52 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
53 datatype = 3
53 datatype = 3
54 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
54 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
55 datatype = 4
55 datatype = 4
56 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
56 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
57 datatype = 5
57 datatype = 5
58 else:
58 else:
59 datatype = None
59 datatype = None
60
60
61 return datatype
61 return datatype
62
62
63
63
64 def hildebrand_sekhon(data, navg):
64 def hildebrand_sekhon(data, navg):
65 """
65 """
66 This method is for the objective determination of the noise level in Doppler spectra. This
66 This method is for the objective determination of the noise level in Doppler spectra. This
67 implementation technique is based on the fact that the standard deviation of the spectral
67 implementation technique is based on the fact that the standard deviation of the spectral
68 densities is equal to the mean spectral density for white Gaussian noise
68 densities is equal to the mean spectral density for white Gaussian noise
69
69
70 Inputs:
70 Inputs:
71 Data : heights
71 Data : heights
72 navg : numbers of averages
72 navg : numbers of averages
73
73
74 Return:
74 Return:
75 mean : noise's level
75 mean : noise's level
76 """
76 """
77
77
78 sortdata = numpy.sort(data, axis=None)
78 sortdata = numpy.sort(data, axis=None)
79 '''
79 '''
80 lenOfData = len(sortdata)
80 lenOfData = len(sortdata)
81 nums_min = lenOfData*0.2
81 nums_min = lenOfData*0.2
82
82
83 if nums_min <= 5:
83 if nums_min <= 5:
84
84
85 nums_min = 5
85 nums_min = 5
86
86
87 sump = 0.
87 sump = 0.
88 sumq = 0.
88 sumq = 0.
89
89
90 j = 0
90 j = 0
91 cont = 1
91 cont = 1
92
92
93 while((cont == 1)and(j < lenOfData)):
93 while((cont == 1)and(j < lenOfData)):
94
94
95 sump += sortdata[j]
95 sump += sortdata[j]
96 sumq += sortdata[j]**2
96 sumq += sortdata[j]**2
97
97
98 if j > nums_min:
98 if j > nums_min:
99 rtest = float(j)/(j-1) + 1.0/navg
99 rtest = float(j)/(j-1) + 1.0/navg
100 if ((sumq*j) > (rtest*sump**2)):
100 if ((sumq*j) > (rtest*sump**2)):
101 j = j - 1
101 j = j - 1
102 sump = sump - sortdata[j]
102 sump = sump - sortdata[j]
103 sumq = sumq - sortdata[j]**2
103 sumq = sumq - sortdata[j]**2
104 cont = 0
104 cont = 0
105
105
106 j += 1
106 j += 1
107
107
108 lnoise = sump / j
108 lnoise = sump / j
109 '''
109 '''
110 return _noise.hildebrand_sekhon(sortdata, navg)
110 return _noise.hildebrand_sekhon(sortdata, navg)
111
111
112
112
113 class Beam:
113 class Beam:
114
114
115 def __init__(self):
115 def __init__(self):
116 self.codeList = []
116 self.codeList = []
117 self.azimuthList = []
117 self.azimuthList = []
118 self.zenithList = []
118 self.zenithList = []
119
119
120
120
121 class GenericData(object):
121 class GenericData(object):
122
122
123 flagNoData = True
123 flagNoData = True
124
124
125 def copy(self, inputObj=None):
125 def copy(self, inputObj=None):
126
126
127 if inputObj == None:
127 if inputObj == None:
128 return copy.deepcopy(self)
128 return copy.deepcopy(self)
129
129
130 for key in list(inputObj.__dict__.keys()):
130 for key in list(inputObj.__dict__.keys()):
131
131
132 attribute = inputObj.__dict__[key]
132 attribute = inputObj.__dict__[key]
133
133
134 # If this attribute is a tuple or list
134 # If this attribute is a tuple or list
135 if type(inputObj.__dict__[key]) in (tuple, list):
135 if type(inputObj.__dict__[key]) in (tuple, list):
136 self.__dict__[key] = attribute[:]
136 self.__dict__[key] = attribute[:]
137 continue
137 continue
138
138
139 # If this attribute is another object or instance
139 # If this attribute is another object or instance
140 if hasattr(attribute, '__dict__'):
140 if hasattr(attribute, '__dict__'):
141 self.__dict__[key] = attribute.copy()
141 self.__dict__[key] = attribute.copy()
142 continue
142 continue
143
143
144 self.__dict__[key] = inputObj.__dict__[key]
144 self.__dict__[key] = inputObj.__dict__[key]
145
145
146 def deepcopy(self):
146 def deepcopy(self):
147
147
148 return copy.deepcopy(self)
148 return copy.deepcopy(self)
149
149
150 def isEmpty(self):
150 def isEmpty(self):
151
151
152 return self.flagNoData
152 return self.flagNoData
153
153
154 def isReady(self):
154 def isReady(self):
155
155
156 return not self.flagNoData
156 return not self.flagNoData
157
157
158
158
159 class JROData(GenericData):
159 class JROData(GenericData):
160
160
161 systemHeaderObj = SystemHeader()
161 systemHeaderObj = SystemHeader()
162 radarControllerHeaderObj = RadarControllerHeader()
162 radarControllerHeaderObj = RadarControllerHeader()
163 type = None
163 type = None
164 datatype = None # dtype but in string
164 datatype = None # dtype but in string
165 nProfiles = None
165 nProfiles = None
166 heightList = None
166 heightList = None
167 channelList = None
167 channelList = None
168 flagDiscontinuousBlock = False
168 flagDiscontinuousBlock = False
169 useLocalTime = False
169 useLocalTime = False
170 utctime = None
170 utctime = None
171 timeZone = None
171 timeZone = None
172 dstFlag = None
172 dstFlag = None
173 errorCount = None
173 errorCount = None
174 blocksize = None
174 blocksize = None
175 flagDecodeData = False # asumo q la data no esta decodificada
175 flagDecodeData = False # asumo q la data no esta decodificada
176 flagDeflipData = False # asumo q la data no esta sin flip
176 flagDeflipData = False # asumo q la data no esta sin flip
177 flagShiftFFT = False
177 flagShiftFFT = False
178 nCohInt = None
178 nCohInt = None
179 windowOfFilter = 1
179 windowOfFilter = 1
180 C = 3e8
180 C = 3e8
181 frequency = 49.92e6
181 frequency = 49.92e6
182 realtime = False
182 realtime = False
183 beacon_heiIndexList = None
183 beacon_heiIndexList = None
184 last_block = None
184 last_block = None
185 blocknow = None
185 blocknow = None
186 azimuth = None
186 azimuth = None
187 zenith = None
187 zenith = None
188 beam = Beam()
188 beam = Beam()
189 profileIndex = None
189 profileIndex = None
190 error = None
190 error = None
191 data = None
191 data = None
192 nmodes = None
192 nmodes = None
193 metadata_list = ['heightList', 'timeZone', 'type']
193 metadata_list = ['heightList', 'timeZone', 'type']
194
194
195 def __str__(self):
195 def __str__(self):
196
196
197 return '{} - {}'.format(self.type, self.datatime())
197 return '{} - {}'.format(self.type, self.datatime())
198
198
199 def getNoise(self):
199 def getNoise(self):
200
200
201 raise NotImplementedError
201 raise NotImplementedError
202
202
203 @property
203 @property
204 def nChannels(self):
204 def nChannels(self):
205
205
206 return len(self.channelList)
206 return len(self.channelList)
207
207
208 @property
208 @property
209 def channelIndexList(self):
209 def channelIndexList(self):
210
210
211 return list(range(self.nChannels))
211 return list(range(self.nChannels))
212
212
213 @property
213 @property
214 def nHeights(self):
214 def nHeights(self):
215
215
216 return len(self.heightList)
216 return len(self.heightList)
217
217
218 def getDeltaH(self):
218 def getDeltaH(self):
219
219
220 return self.heightList[1] - self.heightList[0]
220 return self.heightList[1] - self.heightList[0]
221
221
222 @property
222 @property
223 def ltctime(self):
223 def ltctime(self):
224
224
225 if self.useLocalTime:
225 if self.useLocalTime:
226 return self.utctime - self.timeZone * 60
226 return self.utctime - self.timeZone * 60
227
227
228 return self.utctime
228 return self.utctime
229
229
230 @property
230 @property
231 def datatime(self):
231 def datatime(self):
232
232
233 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
233 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
234 return datatimeValue
234 return datatimeValue
235
235
236 def getTimeRange(self):
236 def getTimeRange(self):
237
237
238 datatime = []
238 datatime = []
239
239
240 datatime.append(self.ltctime)
240 datatime.append(self.ltctime)
241 datatime.append(self.ltctime + self.timeInterval + 1)
241 datatime.append(self.ltctime + self.timeInterval + 1)
242
242
243 datatime = numpy.array(datatime)
243 datatime = numpy.array(datatime)
244
244
245 return datatime
245 return datatime
246
246
247 def getFmaxTimeResponse(self):
247 def getFmaxTimeResponse(self):
248
248
249 period = (10**-6) * self.getDeltaH() / (0.15)
249 period = (10**-6) * self.getDeltaH() / (0.15)
250
250
251 PRF = 1. / (period * self.nCohInt)
251 PRF = 1. / (period * self.nCohInt)
252
252
253 fmax = PRF
253 fmax = PRF
254
254
255 return fmax
255 return fmax
256
256
257 def getFmax(self):
257 def getFmax(self):
258 PRF = 1. / (self.ippSeconds * self.nCohInt)
258 PRF = 1. / (self.ippSeconds * self.nCohInt)
259
259
260 fmax = PRF
260 fmax = PRF
261 return fmax
261 return fmax
262
262
263 def getVmax(self):
263 def getVmax(self):
264
264
265 _lambda = self.C / self.frequency
265 _lambda = self.C / self.frequency
266
266
267 vmax = self.getFmax() * _lambda / 2
267 vmax = self.getFmax() * _lambda / 2
268
268
269 return vmax
269 return vmax
270
270
271 @property
271 @property
272 def ippSeconds(self):
272 def ippSeconds(self):
273 '''
273 '''
274 '''
274 '''
275 return self.radarControllerHeaderObj.ippSeconds
275 return self.radarControllerHeaderObj.ippSeconds
276
276
277 @ippSeconds.setter
277 @ippSeconds.setter
278 def ippSeconds(self, ippSeconds):
278 def ippSeconds(self, ippSeconds):
279 '''
279 '''
280 '''
280 '''
281 self.radarControllerHeaderObj.ippSeconds = ippSeconds
281 self.radarControllerHeaderObj.ippSeconds = ippSeconds
282
282
283 @property
283 @property
284 def code(self):
284 def code(self):
285 '''
285 '''
286 '''
286 '''
287 return self.radarControllerHeaderObj.code
287 return self.radarControllerHeaderObj.code
288
288
289 @code.setter
289 @code.setter
290 def code(self, code):
290 def code(self, code):
291 '''
291 '''
292 '''
292 '''
293 self.radarControllerHeaderObj.code = code
293 self.radarControllerHeaderObj.code = code
294
294
295 @property
295 @property
296 def ncode(self):
296 def nCode(self):
297 '''
297 '''
298 '''
298 '''
299 return self.radarControllerHeaderObj.nCode
299 return self.radarControllerHeaderObj.nCode
300
300
301 @ncode.setter
301 @nCode.setter
302 def ncode(self, ncode):
302 def nCode(self, ncode):
303 '''
303 '''
304 '''
304 '''
305 self.radarControllerHeaderObj.nCode = ncode
305 self.radarControllerHeaderObj.nCode = ncode
306
306
307 @property
307 @property
308 def nbaud(self):
308 def nBaud(self):
309 '''
309 '''
310 '''
310 '''
311 return self.radarControllerHeaderObj.nBaud
311 return self.radarControllerHeaderObj.nBaud
312
312
313 @nbaud.setter
313 @nBaud.setter
314 def nbaud(self, nbaud):
314 def nBaud(self, nbaud):
315 '''
315 '''
316 '''
316 '''
317 self.radarControllerHeaderObj.nBaud = nbaud
317 self.radarControllerHeaderObj.nBaud = nbaud
318
318
319 @property
319 @property
320 def ipp(self):
320 def ipp(self):
321 '''
321 '''
322 '''
322 '''
323 return self.radarControllerHeaderObj.ipp
323 return self.radarControllerHeaderObj.ipp
324
324
325 @ipp.setter
325 @ipp.setter
326 def ipp(self, ipp):
326 def ipp(self, ipp):
327 '''
327 '''
328 '''
328 '''
329 self.radarControllerHeaderObj.ipp = ipp
329 self.radarControllerHeaderObj.ipp = ipp
330
330
331 @property
331 @property
332 def metadata(self):
332 def metadata(self):
333 '''
333 '''
334 '''
334 '''
335
335
336 return {attr: getattr(self, attr) for attr in self.metadata_list}
336 return {attr: getattr(self, attr) for attr in self.metadata_list}
337
337
338
338
339 class Voltage(JROData):
339 class Voltage(JROData):
340
340
341 dataPP_POW = None
341 dataPP_POW = None
342 dataPP_DOP = None
342 dataPP_DOP = None
343 dataPP_WIDTH = None
343 dataPP_WIDTH = None
344 dataPP_SNR = None
344 dataPP_SNR = None
345
345
346 def __init__(self):
346 def __init__(self):
347 '''
347 '''
348 Constructor
348 Constructor
349 '''
349 '''
350
350
351 self.useLocalTime = True
351 self.useLocalTime = True
352 self.radarControllerHeaderObj = RadarControllerHeader()
352 self.radarControllerHeaderObj = RadarControllerHeader()
353 self.systemHeaderObj = SystemHeader()
353 self.systemHeaderObj = SystemHeader()
354 self.type = "Voltage"
354 self.type = "Voltage"
355 self.data = None
355 self.data = None
356 self.nProfiles = None
356 self.nProfiles = None
357 self.heightList = None
357 self.heightList = None
358 self.channelList = None
358 self.channelList = None
359 self.flagNoData = True
359 self.flagNoData = True
360 self.flagDiscontinuousBlock = False
360 self.flagDiscontinuousBlock = False
361 self.utctime = None
361 self.utctime = None
362 self.timeZone = 0
362 self.timeZone = 0
363 self.dstFlag = None
363 self.dstFlag = None
364 self.errorCount = None
364 self.errorCount = None
365 self.nCohInt = None
365 self.nCohInt = None
366 self.blocksize = None
366 self.blocksize = None
367 self.flagCohInt = False
367 self.flagCohInt = False
368 self.flagDecodeData = False # asumo q la data no esta decodificada
368 self.flagDecodeData = False # asumo q la data no esta decodificada
369 self.flagDeflipData = False # asumo q la data no esta sin flip
369 self.flagDeflipData = False # asumo q la data no esta sin flip
370 self.flagShiftFFT = False
370 self.flagShiftFFT = False
371 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
371 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
372 self.profileIndex = 0
372 self.profileIndex = 0
373 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
373 self.metadata_list = ['type', 'heightList', 'timeZone', 'nProfiles', 'channelList', 'nCohInt',
374 'code', 'ncode', 'nbaud', 'ippSeconds', 'ipp']
374 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp']
375
375
376 def getNoisebyHildebrand(self, channel=None):
376 def getNoisebyHildebrand(self, channel=None):
377 """
377 """
378 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
378 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
379
379
380 Return:
380 Return:
381 noiselevel
381 noiselevel
382 """
382 """
383
383
384 if channel != None:
384 if channel != None:
385 data = self.data[channel]
385 data = self.data[channel]
386 nChannels = 1
386 nChannels = 1
387 else:
387 else:
388 data = self.data
388 data = self.data
389 nChannels = self.nChannels
389 nChannels = self.nChannels
390
390
391 noise = numpy.zeros(nChannels)
391 noise = numpy.zeros(nChannels)
392 power = data * numpy.conjugate(data)
392 power = data * numpy.conjugate(data)
393
393
394 for thisChannel in range(nChannels):
394 for thisChannel in range(nChannels):
395 if nChannels == 1:
395 if nChannels == 1:
396 daux = power[:].real
396 daux = power[:].real
397 else:
397 else:
398 daux = power[thisChannel, :].real
398 daux = power[thisChannel, :].real
399 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
399 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
400
400
401 return noise
401 return noise
402
402
403 def getNoise(self, type=1, channel=None):
403 def getNoise(self, type=1, channel=None):
404
404
405 if type == 1:
405 if type == 1:
406 noise = self.getNoisebyHildebrand(channel)
406 noise = self.getNoisebyHildebrand(channel)
407
407
408 return noise
408 return noise
409
409
410 def getPower(self, channel=None):
410 def getPower(self, channel=None):
411
411
412 if channel != None:
412 if channel != None:
413 data = self.data[channel]
413 data = self.data[channel]
414 else:
414 else:
415 data = self.data
415 data = self.data
416
416
417 power = data * numpy.conjugate(data)
417 power = data * numpy.conjugate(data)
418 powerdB = 10 * numpy.log10(power.real)
418 powerdB = 10 * numpy.log10(power.real)
419 powerdB = numpy.squeeze(powerdB)
419 powerdB = numpy.squeeze(powerdB)
420
420
421 return powerdB
421 return powerdB
422
422
423 @property
423 @property
424 def timeInterval(self):
424 def timeInterval(self):
425
425
426 return self.ippSeconds * self.nCohInt
426 return self.ippSeconds * self.nCohInt
427
427
428 noise = property(getNoise, "I'm the 'nHeights' property.")
428 noise = property(getNoise, "I'm the 'nHeights' property.")
429
429
430
430
431 class Spectra(JROData):
431 class Spectra(JROData):
432
432
433 def __init__(self):
433 def __init__(self):
434 '''
434 '''
435 Constructor
435 Constructor
436 '''
436 '''
437
437
438 self.useLocalTime = True
438 self.useLocalTime = True
439 self.radarControllerHeaderObj = RadarControllerHeader()
439 self.radarControllerHeaderObj = RadarControllerHeader()
440 self.systemHeaderObj = SystemHeader()
440 self.systemHeaderObj = SystemHeader()
441 self.type = "Spectra"
441 self.type = "Spectra"
442 self.timeZone = 0
442 self.timeZone = 0
443 self.nProfiles = None
443 self.nProfiles = None
444 self.heightList = None
444 self.heightList = None
445 self.channelList = None
445 self.channelList = None
446 self.pairsList = None
446 self.pairsList = None
447 self.flagNoData = True
447 self.flagNoData = True
448 self.flagDiscontinuousBlock = False
448 self.flagDiscontinuousBlock = False
449 self.utctime = None
449 self.utctime = None
450 self.nCohInt = None
450 self.nCohInt = None
451 self.nIncohInt = None
451 self.nIncohInt = None
452 self.blocksize = None
452 self.blocksize = None
453 self.nFFTPoints = None
453 self.nFFTPoints = None
454 self.wavelength = None
454 self.wavelength = None
455 self.flagDecodeData = False # asumo q la data no esta decodificada
455 self.flagDecodeData = False # asumo q la data no esta decodificada
456 self.flagDeflipData = False # asumo q la data no esta sin flip
456 self.flagDeflipData = False # asumo q la data no esta sin flip
457 self.flagShiftFFT = False
457 self.flagShiftFFT = False
458 self.ippFactor = 1
458 self.ippFactor = 1
459 self.beacon_heiIndexList = []
459 self.beacon_heiIndexList = []
460 self.noise_estimation = None
460 self.noise_estimation = None
461 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
461 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
462 'code', 'ncode', 'nbaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
462 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles']
463
463
464 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
464 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
465 """
465 """
466 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
466 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
467
467
468 Return:
468 Return:
469 noiselevel
469 noiselevel
470 """
470 """
471
471
472 noise = numpy.zeros(self.nChannels)
472 noise = numpy.zeros(self.nChannels)
473
473
474 for channel in range(self.nChannels):
474 for channel in range(self.nChannels):
475 daux = self.data_spc[channel,
475 daux = self.data_spc[channel,
476 xmin_index:xmax_index, ymin_index:ymax_index]
476 xmin_index:xmax_index, ymin_index:ymax_index]
477 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
477 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
478
478
479 return noise
479 return noise
480
480
481 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
481 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
482
482
483 if self.noise_estimation is not None:
483 if self.noise_estimation is not None:
484 # this was estimated by getNoise Operation defined in jroproc_spectra.py
484 # this was estimated by getNoise Operation defined in jroproc_spectra.py
485 return self.noise_estimation
485 return self.noise_estimation
486 else:
486 else:
487 noise = self.getNoisebyHildebrand(
487 noise = self.getNoisebyHildebrand(
488 xmin_index, xmax_index, ymin_index, ymax_index)
488 xmin_index, xmax_index, ymin_index, ymax_index)
489 return noise
489 return noise
490
490
491 def getFreqRangeTimeResponse(self, extrapoints=0):
491 def getFreqRangeTimeResponse(self, extrapoints=0):
492
492
493 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
493 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
494 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
494 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
495
495
496 return freqrange
496 return freqrange
497
497
498 def getAcfRange(self, extrapoints=0):
498 def getAcfRange(self, extrapoints=0):
499
499
500 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
500 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
501 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
501 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
502
502
503 return freqrange
503 return freqrange
504
504
505 def getFreqRange(self, extrapoints=0):
505 def getFreqRange(self, extrapoints=0):
506
506
507 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
507 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
508 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
508 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
509
509
510 return freqrange
510 return freqrange
511
511
512 def getVelRange(self, extrapoints=0):
512 def getVelRange(self, extrapoints=0):
513
513
514 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
514 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
515 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
515 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
516
516
517 if self.nmodes:
517 if self.nmodes:
518 return velrange/self.nmodes
518 return velrange/self.nmodes
519 else:
519 else:
520 return velrange
520 return velrange
521
521
522 @property
522 @property
523 def nPairs(self):
523 def nPairs(self):
524
524
525 return len(self.pairsList)
525 return len(self.pairsList)
526
526
527 @property
527 @property
528 def pairsIndexList(self):
528 def pairsIndexList(self):
529
529
530 return list(range(self.nPairs))
530 return list(range(self.nPairs))
531
531
532 @property
532 @property
533 def normFactor(self):
533 def normFactor(self):
534
534
535 pwcode = 1
535 pwcode = 1
536
536
537 if self.flagDecodeData:
537 if self.flagDecodeData:
538 pwcode = numpy.sum(self.code[0]**2)
538 pwcode = numpy.sum(self.code[0]**2)
539 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
539 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
540 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
540 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
541
541
542 return normFactor
542 return normFactor
543
543
544 @property
544 @property
545 def flag_cspc(self):
545 def flag_cspc(self):
546
546
547 if self.data_cspc is None:
547 if self.data_cspc is None:
548 return True
548 return True
549
549
550 return False
550 return False
551
551
552 @property
552 @property
553 def flag_dc(self):
553 def flag_dc(self):
554
554
555 if self.data_dc is None:
555 if self.data_dc is None:
556 return True
556 return True
557
557
558 return False
558 return False
559
559
560 @property
560 @property
561 def timeInterval(self):
561 def timeInterval(self):
562
562
563 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
563 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
564 if self.nmodes:
564 if self.nmodes:
565 return self.nmodes*timeInterval
565 return self.nmodes*timeInterval
566 else:
566 else:
567 return timeInterval
567 return timeInterval
568
568
569 def getPower(self):
569 def getPower(self):
570
570
571 factor = self.normFactor
571 factor = self.normFactor
572 z = self.data_spc / factor
572 z = self.data_spc / factor
573 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
573 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
574 avg = numpy.average(z, axis=1)
574 avg = numpy.average(z, axis=1)
575
575
576 return 10 * numpy.log10(avg)
576 return 10 * numpy.log10(avg)
577
577
578 def getCoherence(self, pairsList=None, phase=False):
578 def getCoherence(self, pairsList=None, phase=False):
579
579
580 z = []
580 z = []
581 if pairsList is None:
581 if pairsList is None:
582 pairsIndexList = self.pairsIndexList
582 pairsIndexList = self.pairsIndexList
583 else:
583 else:
584 pairsIndexList = []
584 pairsIndexList = []
585 for pair in pairsList:
585 for pair in pairsList:
586 if pair not in self.pairsList:
586 if pair not in self.pairsList:
587 raise ValueError("Pair %s is not in dataOut.pairsList" % (
587 raise ValueError("Pair %s is not in dataOut.pairsList" % (
588 pair))
588 pair))
589 pairsIndexList.append(self.pairsList.index(pair))
589 pairsIndexList.append(self.pairsList.index(pair))
590 for i in range(len(pairsIndexList)):
590 for i in range(len(pairsIndexList)):
591 pair = self.pairsList[pairsIndexList[i]]
591 pair = self.pairsList[pairsIndexList[i]]
592 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
592 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
593 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
593 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
594 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
594 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
595 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
595 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
596 if phase:
596 if phase:
597 data = numpy.arctan2(avgcoherenceComplex.imag,
597 data = numpy.arctan2(avgcoherenceComplex.imag,
598 avgcoherenceComplex.real) * 180 / numpy.pi
598 avgcoherenceComplex.real) * 180 / numpy.pi
599 else:
599 else:
600 data = numpy.abs(avgcoherenceComplex)
600 data = numpy.abs(avgcoherenceComplex)
601
601
602 z.append(data)
602 z.append(data)
603
603
604 return numpy.array(z)
604 return numpy.array(z)
605
605
606 def setValue(self, value):
606 def setValue(self, value):
607
607
608 print("This property should not be initialized")
608 print("This property should not be initialized")
609
609
610 return
610 return
611
611
612 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
612 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
613
613
614
614
615 class SpectraHeis(Spectra):
615 class SpectraHeis(Spectra):
616
616
617 def __init__(self):
617 def __init__(self):
618
618
619 self.radarControllerHeaderObj = RadarControllerHeader()
619 self.radarControllerHeaderObj = RadarControllerHeader()
620 self.systemHeaderObj = SystemHeader()
620 self.systemHeaderObj = SystemHeader()
621 self.type = "SpectraHeis"
621 self.type = "SpectraHeis"
622 self.nProfiles = None
622 self.nProfiles = None
623 self.heightList = None
623 self.heightList = None
624 self.channelList = None
624 self.channelList = None
625 self.flagNoData = True
625 self.flagNoData = True
626 self.flagDiscontinuousBlock = False
626 self.flagDiscontinuousBlock = False
627 self.utctime = None
627 self.utctime = None
628 self.blocksize = None
628 self.blocksize = None
629 self.profileIndex = 0
629 self.profileIndex = 0
630 self.nCohInt = 1
630 self.nCohInt = 1
631 self.nIncohInt = 1
631 self.nIncohInt = 1
632
632
633 @property
633 @property
634 def normFactor(self):
634 def normFactor(self):
635 pwcode = 1
635 pwcode = 1
636 if self.flagDecodeData:
636 if self.flagDecodeData:
637 pwcode = numpy.sum(self.code[0]**2)
637 pwcode = numpy.sum(self.code[0]**2)
638
638
639 normFactor = self.nIncohInt * self.nCohInt * pwcode
639 normFactor = self.nIncohInt * self.nCohInt * pwcode
640
640
641 return normFactor
641 return normFactor
642
642
643 @property
643 @property
644 def timeInterval(self):
644 def timeInterval(self):
645
645
646 return self.ippSeconds * self.nCohInt * self.nIncohInt
646 return self.ippSeconds * self.nCohInt * self.nIncohInt
647
647
648
648
649 class Fits(JROData):
649 class Fits(JROData):
650
650
651 def __init__(self):
651 def __init__(self):
652
652
653 self.type = "Fits"
653 self.type = "Fits"
654 self.nProfiles = None
654 self.nProfiles = None
655 self.heightList = None
655 self.heightList = None
656 self.channelList = None
656 self.channelList = None
657 self.flagNoData = True
657 self.flagNoData = True
658 self.utctime = None
658 self.utctime = None
659 self.nCohInt = 1
659 self.nCohInt = 1
660 self.nIncohInt = 1
660 self.nIncohInt = 1
661 self.useLocalTime = True
661 self.useLocalTime = True
662 self.profileIndex = 0
662 self.profileIndex = 0
663 self.timeZone = 0
663 self.timeZone = 0
664
664
665 def getTimeRange(self):
665 def getTimeRange(self):
666
666
667 datatime = []
667 datatime = []
668
668
669 datatime.append(self.ltctime)
669 datatime.append(self.ltctime)
670 datatime.append(self.ltctime + self.timeInterval)
670 datatime.append(self.ltctime + self.timeInterval)
671
671
672 datatime = numpy.array(datatime)
672 datatime = numpy.array(datatime)
673
673
674 return datatime
674 return datatime
675
675
676 def getChannelIndexList(self):
676 def getChannelIndexList(self):
677
677
678 return list(range(self.nChannels))
678 return list(range(self.nChannels))
679
679
680 def getNoise(self, type=1):
680 def getNoise(self, type=1):
681
681
682
682
683 if type == 1:
683 if type == 1:
684 noise = self.getNoisebyHildebrand()
684 noise = self.getNoisebyHildebrand()
685
685
686 if type == 2:
686 if type == 2:
687 noise = self.getNoisebySort()
687 noise = self.getNoisebySort()
688
688
689 if type == 3:
689 if type == 3:
690 noise = self.getNoisebyWindow()
690 noise = self.getNoisebyWindow()
691
691
692 return noise
692 return noise
693
693
694 @property
694 @property
695 def timeInterval(self):
695 def timeInterval(self):
696
696
697 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
697 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
698
698
699 return timeInterval
699 return timeInterval
700
700
701 @property
701 @property
702 def ippSeconds(self):
702 def ippSeconds(self):
703 '''
703 '''
704 '''
704 '''
705 return self.ipp_sec
705 return self.ipp_sec
706
706
707 noise = property(getNoise, "I'm the 'nHeights' property.")
707 noise = property(getNoise, "I'm the 'nHeights' property.")
708
708
709
709
710 class Correlation(JROData):
710 class Correlation(JROData):
711
711
712 def __init__(self):
712 def __init__(self):
713 '''
713 '''
714 Constructor
714 Constructor
715 '''
715 '''
716 self.radarControllerHeaderObj = RadarControllerHeader()
716 self.radarControllerHeaderObj = RadarControllerHeader()
717 self.systemHeaderObj = SystemHeader()
717 self.systemHeaderObj = SystemHeader()
718 self.type = "Correlation"
718 self.type = "Correlation"
719 self.data = None
719 self.data = None
720 self.dtype = None
720 self.dtype = None
721 self.nProfiles = None
721 self.nProfiles = None
722 self.heightList = None
722 self.heightList = None
723 self.channelList = None
723 self.channelList = None
724 self.flagNoData = True
724 self.flagNoData = True
725 self.flagDiscontinuousBlock = False
725 self.flagDiscontinuousBlock = False
726 self.utctime = None
726 self.utctime = None
727 self.timeZone = 0
727 self.timeZone = 0
728 self.dstFlag = None
728 self.dstFlag = None
729 self.errorCount = None
729 self.errorCount = None
730 self.blocksize = None
730 self.blocksize = None
731 self.flagDecodeData = False # asumo q la data no esta decodificada
731 self.flagDecodeData = False # asumo q la data no esta decodificada
732 self.flagDeflipData = False # asumo q la data no esta sin flip
732 self.flagDeflipData = False # asumo q la data no esta sin flip
733 self.pairsList = None
733 self.pairsList = None
734 self.nPoints = None
734 self.nPoints = None
735
735
736 def getPairsList(self):
736 def getPairsList(self):
737
737
738 return self.pairsList
738 return self.pairsList
739
739
740 def getNoise(self, mode=2):
740 def getNoise(self, mode=2):
741
741
742 indR = numpy.where(self.lagR == 0)[0][0]
742 indR = numpy.where(self.lagR == 0)[0][0]
743 indT = numpy.where(self.lagT == 0)[0][0]
743 indT = numpy.where(self.lagT == 0)[0][0]
744
744
745 jspectra0 = self.data_corr[:, :, indR, :]
745 jspectra0 = self.data_corr[:, :, indR, :]
746 jspectra = copy.copy(jspectra0)
746 jspectra = copy.copy(jspectra0)
747
747
748 num_chan = jspectra.shape[0]
748 num_chan = jspectra.shape[0]
749 num_hei = jspectra.shape[2]
749 num_hei = jspectra.shape[2]
750
750
751 freq_dc = jspectra.shape[1] / 2
751 freq_dc = jspectra.shape[1] / 2
752 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
752 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
753
753
754 if ind_vel[0] < 0:
754 if ind_vel[0] < 0:
755 ind_vel[list(range(0, 1))] = ind_vel[list(
755 ind_vel[list(range(0, 1))] = ind_vel[list(
756 range(0, 1))] + self.num_prof
756 range(0, 1))] + self.num_prof
757
757
758 if mode == 1:
758 if mode == 1:
759 jspectra[:, freq_dc, :] = (
759 jspectra[:, freq_dc, :] = (
760 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
760 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
761
761
762 if mode == 2:
762 if mode == 2:
763
763
764 vel = numpy.array([-2, -1, 1, 2])
764 vel = numpy.array([-2, -1, 1, 2])
765 xx = numpy.zeros([4, 4])
765 xx = numpy.zeros([4, 4])
766
766
767 for fil in range(4):
767 for fil in range(4):
768 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
768 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
769
769
770 xx_inv = numpy.linalg.inv(xx)
770 xx_inv = numpy.linalg.inv(xx)
771 xx_aux = xx_inv[0, :]
771 xx_aux = xx_inv[0, :]
772
772
773 for ich in range(num_chan):
773 for ich in range(num_chan):
774 yy = jspectra[ich, ind_vel, :]
774 yy = jspectra[ich, ind_vel, :]
775 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
775 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
776
776
777 junkid = jspectra[ich, freq_dc, :] <= 0
777 junkid = jspectra[ich, freq_dc, :] <= 0
778 cjunkid = sum(junkid)
778 cjunkid = sum(junkid)
779
779
780 if cjunkid.any():
780 if cjunkid.any():
781 jspectra[ich, freq_dc, junkid.nonzero()] = (
781 jspectra[ich, freq_dc, junkid.nonzero()] = (
782 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
782 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
783
783
784 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
784 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
785
785
786 return noise
786 return noise
787
787
788 @property
788 @property
789 def timeInterval(self):
789 def timeInterval(self):
790
790
791 return self.ippSeconds * self.nCohInt * self.nProfiles
791 return self.ippSeconds * self.nCohInt * self.nProfiles
792
792
793 def splitFunctions(self):
793 def splitFunctions(self):
794
794
795 pairsList = self.pairsList
795 pairsList = self.pairsList
796 ccf_pairs = []
796 ccf_pairs = []
797 acf_pairs = []
797 acf_pairs = []
798 ccf_ind = []
798 ccf_ind = []
799 acf_ind = []
799 acf_ind = []
800 for l in range(len(pairsList)):
800 for l in range(len(pairsList)):
801 chan0 = pairsList[l][0]
801 chan0 = pairsList[l][0]
802 chan1 = pairsList[l][1]
802 chan1 = pairsList[l][1]
803
803
804 # Obteniendo pares de Autocorrelacion
804 # Obteniendo pares de Autocorrelacion
805 if chan0 == chan1:
805 if chan0 == chan1:
806 acf_pairs.append(chan0)
806 acf_pairs.append(chan0)
807 acf_ind.append(l)
807 acf_ind.append(l)
808 else:
808 else:
809 ccf_pairs.append(pairsList[l])
809 ccf_pairs.append(pairsList[l])
810 ccf_ind.append(l)
810 ccf_ind.append(l)
811
811
812 data_acf = self.data_cf[acf_ind]
812 data_acf = self.data_cf[acf_ind]
813 data_ccf = self.data_cf[ccf_ind]
813 data_ccf = self.data_cf[ccf_ind]
814
814
815 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
815 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
816
816
817 @property
817 @property
818 def normFactor(self):
818 def normFactor(self):
819 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
819 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
820 acf_pairs = numpy.array(acf_pairs)
820 acf_pairs = numpy.array(acf_pairs)
821 normFactor = numpy.zeros((self.nPairs, self.nHeights))
821 normFactor = numpy.zeros((self.nPairs, self.nHeights))
822
822
823 for p in range(self.nPairs):
823 for p in range(self.nPairs):
824 pair = self.pairsList[p]
824 pair = self.pairsList[p]
825
825
826 ch0 = pair[0]
826 ch0 = pair[0]
827 ch1 = pair[1]
827 ch1 = pair[1]
828
828
829 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
829 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
830 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
830 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
831 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
831 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
832
832
833 return normFactor
833 return normFactor
834
834
835
835
836 class Parameters(Spectra):
836 class Parameters(Spectra):
837
837
838 groupList = None # List of Pairs, Groups, etc
838 groupList = None # List of Pairs, Groups, etc
839 data_param = None # Parameters obtained
839 data_param = None # Parameters obtained
840 data_pre = None # Data Pre Parametrization
840 data_pre = None # Data Pre Parametrization
841 data_SNR = None # Signal to Noise Ratio
841 data_SNR = None # Signal to Noise Ratio
842 abscissaList = None # Abscissa, can be velocities, lags or time
842 abscissaList = None # Abscissa, can be velocities, lags or time
843 utctimeInit = None # Initial UTC time
843 utctimeInit = None # Initial UTC time
844 paramInterval = None # Time interval to calculate Parameters in seconds
844 paramInterval = None # Time interval to calculate Parameters in seconds
845 useLocalTime = True
845 useLocalTime = True
846 # Fitting
846 # Fitting
847 data_error = None # Error of the estimation
847 data_error = None # Error of the estimation
848 constants = None
848 constants = None
849 library = None
849 library = None
850 # Output signal
850 # Output signal
851 outputInterval = None # Time interval to calculate output signal in seconds
851 outputInterval = None # Time interval to calculate output signal in seconds
852 data_output = None # Out signal
852 data_output = None # Out signal
853 nAvg = None
853 nAvg = None
854 noise_estimation = None
854 noise_estimation = None
855 GauSPC = None # Fit gaussian SPC
855 GauSPC = None # Fit gaussian SPC
856
856
857 def __init__(self):
857 def __init__(self):
858 '''
858 '''
859 Constructor
859 Constructor
860 '''
860 '''
861 self.radarControllerHeaderObj = RadarControllerHeader()
861 self.radarControllerHeaderObj = RadarControllerHeader()
862 self.systemHeaderObj = SystemHeader()
862 self.systemHeaderObj = SystemHeader()
863 self.type = "Parameters"
863 self.type = "Parameters"
864 self.timeZone = 0
864 self.timeZone = 0
865
865
866 def getTimeRange1(self, interval):
866 def getTimeRange1(self, interval):
867
867
868 datatime = []
868 datatime = []
869
869
870 if self.useLocalTime:
870 if self.useLocalTime:
871 time1 = self.utctimeInit - self.timeZone * 60
871 time1 = self.utctimeInit - self.timeZone * 60
872 else:
872 else:
873 time1 = self.utctimeInit
873 time1 = self.utctimeInit
874
874
875 datatime.append(time1)
875 datatime.append(time1)
876 datatime.append(time1 + interval)
876 datatime.append(time1 + interval)
877 datatime = numpy.array(datatime)
877 datatime = numpy.array(datatime)
878
878
879 return datatime
879 return datatime
880
880
881 @property
881 @property
882 def timeInterval(self):
882 def timeInterval(self):
883
883
884 if hasattr(self, 'timeInterval1'):
884 if hasattr(self, 'timeInterval1'):
885 return self.timeInterval1
885 return self.timeInterval1
886 else:
886 else:
887 return self.paramInterval
887 return self.paramInterval
888
888
889 def setValue(self, value):
889 def setValue(self, value):
890
890
891 print("This property should not be initialized")
891 print("This property should not be initialized")
892
892
893 return
893 return
894
894
895 def getNoise(self):
895 def getNoise(self):
896
896
897 return self.spc_noise
897 return self.spc_noise
898
898
899 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
899 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
900
900
901
901
902 class PlotterData(object):
902 class PlotterData(object):
903 '''
903 '''
904 Object to hold data to be plotted
904 Object to hold data to be plotted
905 '''
905 '''
906
906
907 MAXNUMX = 200
907 MAXNUMX = 200
908 MAXNUMY = 200
908 MAXNUMY = 200
909
909
910 def __init__(self, code, throttle_value, exp_code, localtime=True, buffering=True, snr=False):
910 def __init__(self, code, throttle_value, exp_code, localtime=True, buffering=True, snr=False):
911
911
912 self.key = code
912 self.key = code
913 self.throttle = throttle_value
913 self.throttle = throttle_value
914 self.exp_code = exp_code
914 self.exp_code = exp_code
915 self.buffering = buffering
915 self.buffering = buffering
916 self.ready = False
916 self.ready = False
917 self.flagNoData = False
917 self.flagNoData = False
918 self.localtime = localtime
918 self.localtime = localtime
919 self.data = {}
919 self.data = {}
920 self.meta = {}
920 self.meta = {}
921 self.__heights = []
921 self.__heights = []
922
922
923 if 'snr' in code:
923 if 'snr' in code:
924 self.plottypes = ['snr']
924 self.plottypes = ['snr']
925 elif code == 'spc':
925 elif code == 'spc':
926 self.plottypes = ['spc', 'noise', 'rti']
926 self.plottypes = ['spc', 'noise', 'rti']
927 elif code == 'cspc':
927 elif code == 'cspc':
928 self.plottypes = ['cspc', 'spc', 'noise', 'rti']
928 self.plottypes = ['cspc', 'spc', 'noise', 'rti']
929 elif code == 'rti':
929 elif code == 'rti':
930 self.plottypes = ['noise', 'rti']
930 self.plottypes = ['noise', 'rti']
931 else:
931 else:
932 self.plottypes = [code]
932 self.plottypes = [code]
933
933
934 if 'snr' not in self.plottypes and snr:
934 if 'snr' not in self.plottypes and snr:
935 self.plottypes.append('snr')
935 self.plottypes.append('snr')
936
936
937 for plot in self.plottypes:
937 for plot in self.plottypes:
938 self.data[plot] = {}
938 self.data[plot] = {}
939
939
940 def __str__(self):
940 def __str__(self):
941 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
941 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
942 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
942 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
943
943
944 def __len__(self):
944 def __len__(self):
945 return len(self.data[self.key])
945 return len(self.data[self.key])
946
946
947 def __getitem__(self, key):
947 def __getitem__(self, key):
948
948
949 if key not in self.data:
949 if key not in self.data:
950 raise KeyError(log.error('Missing key: {}'.format(key)))
950 raise KeyError(log.error('Missing key: {}'.format(key)))
951 if 'spc' in key or not self.buffering:
951 if 'spc' in key or not self.buffering:
952 ret = self.data[key][self.tm]
952 ret = self.data[key][self.tm]
953 elif 'scope' in key:
953 elif 'scope' in key:
954 ret = numpy.array(self.data[key][float(self.tm)])
954 ret = numpy.array(self.data[key][float(self.tm)])
955 else:
955 else:
956 ret = numpy.array([self.data[key][x] for x in self.times])
956 ret = numpy.array([self.data[key][x] for x in self.times])
957 if ret.ndim > 1:
957 if ret.ndim > 1:
958 ret = numpy.swapaxes(ret, 0, 1)
958 ret = numpy.swapaxes(ret, 0, 1)
959 return ret
959 return ret
960
960
961 def __contains__(self, key):
961 def __contains__(self, key):
962 return key in self.data
962 return key in self.data
963
963
964 def setup(self):
964 def setup(self):
965 '''
965 '''
966 Configure object
966 Configure object
967 '''
967 '''
968 self.type = ''
968 self.type = ''
969 self.ready = False
969 self.ready = False
970 del self.data
970 del self.data
971 self.data = {}
971 self.data = {}
972 self.__heights = []
972 self.__heights = []
973 self.__all_heights = set()
973 self.__all_heights = set()
974 for plot in self.plottypes:
974 for plot in self.plottypes:
975 if 'snr' in plot:
975 if 'snr' in plot:
976 plot = 'snr'
976 plot = 'snr'
977 elif 'spc_moments' == plot:
977 elif 'spc_moments' == plot:
978 plot = 'moments'
978 plot = 'moments'
979 self.data[plot] = {}
979 self.data[plot] = {}
980
980
981 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
981 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
982 self.data['noise'] = {}
982 self.data['noise'] = {}
983 self.data['rti'] = {}
983 self.data['rti'] = {}
984 if 'noise' not in self.plottypes:
984 if 'noise' not in self.plottypes:
985 self.plottypes.append('noise')
985 self.plottypes.append('noise')
986 if 'rti' not in self.plottypes:
986 if 'rti' not in self.plottypes:
987 self.plottypes.append('rti')
987 self.plottypes.append('rti')
988
988
989 def shape(self, key):
989 def shape(self, key):
990 '''
990 '''
991 Get the shape of the one-element data for the given key
991 Get the shape of the one-element data for the given key
992 '''
992 '''
993
993
994 if len(self.data[key]):
994 if len(self.data[key]):
995 if 'spc' in key or not self.buffering:
995 if 'spc' in key or not self.buffering:
996 return self.data[key].shape
996 return self.data[key].shape
997 return self.data[key][self.times[0]].shape
997 return self.data[key][self.times[0]].shape
998 return (0,)
998 return (0,)
999
999
1000 def update(self, dataOut, tm):
1000 def update(self, dataOut, tm):
1001 '''
1001 '''
1002 Update data object with new dataOut
1002 Update data object with new dataOut
1003 '''
1003 '''
1004
1004
1005 self.profileIndex = dataOut.profileIndex
1005 self.profileIndex = dataOut.profileIndex
1006 self.tm = tm
1006 self.tm = tm
1007 self.type = dataOut.type
1007 self.type = dataOut.type
1008 self.parameters = getattr(dataOut, 'parameters', [])
1008 self.parameters = getattr(dataOut, 'parameters', [])
1009
1009
1010 if hasattr(dataOut, 'meta'):
1010 if hasattr(dataOut, 'meta'):
1011 self.meta.update(dataOut.meta)
1011 self.meta.update(dataOut.meta)
1012
1012
1013 if hasattr(dataOut, 'pairsList'):
1013 if hasattr(dataOut, 'pairsList'):
1014 self.pairs = dataOut.pairsList
1014 self.pairs = dataOut.pairsList
1015
1015
1016 self.interval = dataOut.timeInterval
1016 self.interval = dataOut.timeInterval
1017 if True in ['spc' in ptype for ptype in self.plottypes]:
1017 if True in ['spc' in ptype for ptype in self.plottypes]:
1018 self.xrange = (dataOut.getFreqRange(1)/1000.,
1018 self.xrange = (dataOut.getFreqRange(1)/1000.,
1019 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1019 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1020 self.__heights.append(dataOut.heightList)
1020 self.__heights.append(dataOut.heightList)
1021 self.__all_heights.update(dataOut.heightList)
1021 self.__all_heights.update(dataOut.heightList)
1022
1022
1023 for plot in self.plottypes:
1023 for plot in self.plottypes:
1024 if plot in ('spc', 'spc_moments', 'spc_cut'):
1024 if plot in ('spc', 'spc_moments', 'spc_cut'):
1025 z = dataOut.data_spc/dataOut.normFactor
1025 z = dataOut.data_spc/dataOut.normFactor
1026 buffer = 10*numpy.log10(z)
1026 buffer = 10*numpy.log10(z)
1027 if plot == 'cspc':
1027 if plot == 'cspc':
1028 buffer = (dataOut.data_spc, dataOut.data_cspc)
1028 buffer = (dataOut.data_spc, dataOut.data_cspc)
1029 if plot == 'noise':
1029 if plot == 'noise':
1030 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1030 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1031 if plot in ('rti', 'spcprofile'):
1031 if plot in ('rti', 'spcprofile'):
1032 buffer = dataOut.getPower()
1032 buffer = dataOut.getPower()
1033 if plot == 'snr_db':
1033 if plot == 'snr_db':
1034 buffer = dataOut.data_SNR
1034 buffer = dataOut.data_SNR
1035 if plot == 'snr':
1035 if plot == 'snr':
1036 buffer = 10*numpy.log10(dataOut.data_SNR)
1036 buffer = 10*numpy.log10(dataOut.data_SNR)
1037 if plot == 'dop':
1037 if plot == 'dop':
1038 buffer = dataOut.data_DOP
1038 buffer = dataOut.data_DOP
1039 if plot == 'pow':
1039 if plot == 'pow':
1040 buffer = 10*numpy.log10(dataOut.data_POW)
1040 buffer = 10*numpy.log10(dataOut.data_POW)
1041 if plot == 'width':
1041 if plot == 'width':
1042 buffer = dataOut.data_WIDTH
1042 buffer = dataOut.data_WIDTH
1043 if plot == 'coh':
1043 if plot == 'coh':
1044 buffer = dataOut.getCoherence()
1044 buffer = dataOut.getCoherence()
1045 if plot == 'phase':
1045 if plot == 'phase':
1046 buffer = dataOut.getCoherence(phase=True)
1046 buffer = dataOut.getCoherence(phase=True)
1047 if plot == 'output':
1047 if plot == 'output':
1048 buffer = dataOut.data_output
1048 buffer = dataOut.data_output
1049 if plot == 'param':
1049 if plot == 'param':
1050 buffer = dataOut.data_param
1050 buffer = dataOut.data_param
1051 if plot == 'scope':
1051 if plot == 'scope':
1052 buffer = dataOut.data
1052 buffer = dataOut.data
1053 self.flagDataAsBlock = dataOut.flagDataAsBlock
1053 self.flagDataAsBlock = dataOut.flagDataAsBlock
1054 self.nProfiles = dataOut.nProfiles
1054 self.nProfiles = dataOut.nProfiles
1055 if plot == 'pp_power':
1055 if plot == 'pp_power':
1056 buffer = dataOut.dataPP_POWER
1056 buffer = dataOut.dataPP_POWER
1057 self.flagDataAsBlock = dataOut.flagDataAsBlock
1057 self.flagDataAsBlock = dataOut.flagDataAsBlock
1058 self.nProfiles = dataOut.nProfiles
1058 self.nProfiles = dataOut.nProfiles
1059 if plot == 'pp_signal':
1059 if plot == 'pp_signal':
1060 buffer = dataOut.dataPP_POW
1060 buffer = dataOut.dataPP_POW
1061 self.flagDataAsBlock = dataOut.flagDataAsBlock
1061 self.flagDataAsBlock = dataOut.flagDataAsBlock
1062 self.nProfiles = dataOut.nProfiles
1062 self.nProfiles = dataOut.nProfiles
1063 if plot == 'pp_velocity':
1063 if plot == 'pp_velocity':
1064 buffer = dataOut.dataPP_DOP
1064 buffer = dataOut.dataPP_DOP
1065 self.flagDataAsBlock = dataOut.flagDataAsBlock
1065 self.flagDataAsBlock = dataOut.flagDataAsBlock
1066 self.nProfiles = dataOut.nProfiles
1066 self.nProfiles = dataOut.nProfiles
1067 if plot == 'pp_specwidth':
1067 if plot == 'pp_specwidth':
1068 buffer = dataOut.dataPP_WIDTH
1068 buffer = dataOut.dataPP_WIDTH
1069 self.flagDataAsBlock = dataOut.flagDataAsBlock
1069 self.flagDataAsBlock = dataOut.flagDataAsBlock
1070 self.nProfiles = dataOut.nProfiles
1070 self.nProfiles = dataOut.nProfiles
1071
1071
1072 if plot == 'spc':
1072 if plot == 'spc':
1073 self.data['spc'][tm] = buffer
1073 self.data['spc'][tm] = buffer
1074 elif plot == 'cspc':
1074 elif plot == 'cspc':
1075 self.data['cspc'][tm] = buffer
1075 self.data['cspc'][tm] = buffer
1076 elif plot == 'spc_moments':
1076 elif plot == 'spc_moments':
1077 self.data['spc'][tm] = buffer
1077 self.data['spc'][tm] = buffer
1078 self.data['moments'][tm] = dataOut.moments
1078 self.data['moments'][tm] = dataOut.moments
1079 else:
1079 else:
1080 if self.buffering:
1080 if self.buffering:
1081 self.data[plot][tm] = buffer
1081 self.data[plot][tm] = buffer
1082 else:
1082 else:
1083 self.data[plot][tm] = buffer
1083 self.data[plot][tm] = buffer
1084
1084
1085 if dataOut.channelList is None:
1085 if dataOut.channelList is None:
1086 self.channels = range(buffer.shape[0])
1086 self.channels = range(buffer.shape[0])
1087 else:
1087 else:
1088 self.channels = dataOut.channelList
1088 self.channels = dataOut.channelList
1089
1089
1090 if buffer is None:
1090 if buffer is None:
1091 self.flagNoData = True
1091 self.flagNoData = True
1092 raise schainpy.admin.SchainWarning('Attribute data_{} is empty'.format(self.key))
1092 raise schainpy.admin.SchainWarning('Attribute data_{} is empty'.format(self.key))
1093
1093
1094 def normalize_heights(self):
1094 def normalize_heights(self):
1095 '''
1095 '''
1096 Ensure same-dimension of the data for different heighList
1096 Ensure same-dimension of the data for different heighList
1097 '''
1097 '''
1098
1098
1099 H = numpy.array(list(self.__all_heights))
1099 H = numpy.array(list(self.__all_heights))
1100 H.sort()
1100 H.sort()
1101 for key in self.data:
1101 for key in self.data:
1102 shape = self.shape(key)[:-1] + H.shape
1102 shape = self.shape(key)[:-1] + H.shape
1103 for tm, obj in list(self.data[key].items()):
1103 for tm, obj in list(self.data[key].items()):
1104 h = self.__heights[self.times.tolist().index(tm)]
1104 h = self.__heights[self.times.tolist().index(tm)]
1105 if H.size == h.size:
1105 if H.size == h.size:
1106 continue
1106 continue
1107 index = numpy.where(numpy.in1d(H, h))[0]
1107 index = numpy.where(numpy.in1d(H, h))[0]
1108 dummy = numpy.zeros(shape) + numpy.nan
1108 dummy = numpy.zeros(shape) + numpy.nan
1109 if len(shape) == 2:
1109 if len(shape) == 2:
1110 dummy[:, index] = obj
1110 dummy[:, index] = obj
1111 else:
1111 else:
1112 dummy[index] = obj
1112 dummy[index] = obj
1113 self.data[key][tm] = dummy
1113 self.data[key][tm] = dummy
1114
1114
1115 self.__heights = [H for tm in self.times]
1115 self.__heights = [H for tm in self.times]
1116
1116
1117 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1117 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1118 '''
1118 '''
1119 Convert data to json
1119 Convert data to json
1120 '''
1120 '''
1121
1121
1122 dy = int(self.heights.size/self.MAXNUMY) + 1
1122 dy = int(self.heights.size/self.MAXNUMY) + 1
1123 if self.key in ('spc', 'cspc'):
1123 if self.key in ('spc', 'cspc'):
1124 dx = int(self.data[self.key][tm].shape[1]/self.MAXNUMX) + 1
1124 dx = int(self.data[self.key][tm].shape[1]/self.MAXNUMX) + 1
1125 data = self.roundFloats(
1125 data = self.roundFloats(
1126 self.data[self.key][tm][::, ::dx, ::dy].tolist())
1126 self.data[self.key][tm][::, ::dx, ::dy].tolist())
1127 else:
1127 else:
1128 if self.key is 'noise':
1128 if self.key is 'noise':
1129 data = [[x] for x in self.roundFloats(self.data[self.key][tm].tolist())]
1129 data = [[x] for x in self.roundFloats(self.data[self.key][tm].tolist())]
1130 else:
1130 else:
1131 data = self.roundFloats(self.data[self.key][tm][::, ::dy].tolist())
1131 data = self.roundFloats(self.data[self.key][tm][::, ::dy].tolist())
1132
1132
1133 meta = {}
1133 meta = {}
1134 ret = {
1134 ret = {
1135 'plot': plot_name,
1135 'plot': plot_name,
1136 'code': self.exp_code,
1136 'code': self.exp_code,
1137 'time': float(tm),
1137 'time': float(tm),
1138 'data': data,
1138 'data': data,
1139 }
1139 }
1140 meta['type'] = plot_type
1140 meta['type'] = plot_type
1141 meta['interval'] = float(self.interval)
1141 meta['interval'] = float(self.interval)
1142 meta['localtime'] = self.localtime
1142 meta['localtime'] = self.localtime
1143 meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1143 meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1144 if 'spc' in self.data or 'cspc' in self.data:
1144 if 'spc' in self.data or 'cspc' in self.data:
1145 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1145 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1146 else:
1146 else:
1147 meta['xrange'] = []
1147 meta['xrange'] = []
1148
1148
1149 meta.update(self.meta)
1149 meta.update(self.meta)
1150 ret['metadata'] = meta
1150 ret['metadata'] = meta
1151 return json.dumps(ret)
1151 return json.dumps(ret)
1152
1152
1153 @property
1153 @property
1154 def times(self):
1154 def times(self):
1155 '''
1155 '''
1156 Return the list of times of the current data
1156 Return the list of times of the current data
1157 '''
1157 '''
1158
1158
1159 ret = numpy.array([t for t in self.data[self.key]])
1159 ret = numpy.array([t for t in self.data[self.key]])
1160 if self:
1160 if self:
1161 ret.sort()
1161 ret.sort()
1162 return ret
1162 return ret
1163
1163
1164 @property
1164 @property
1165 def min_time(self):
1165 def min_time(self):
1166 '''
1166 '''
1167 Return the minimun time value
1167 Return the minimun time value
1168 '''
1168 '''
1169
1169
1170 return self.times[0]
1170 return self.times[0]
1171
1171
1172 @property
1172 @property
1173 def max_time(self):
1173 def max_time(self):
1174 '''
1174 '''
1175 Return the maximun time value
1175 Return the maximun time value
1176 '''
1176 '''
1177
1177
1178 return self.times[-1]
1178 return self.times[-1]
1179
1179
1180 @property
1180 @property
1181 def heights(self):
1181 def heights(self):
1182 '''
1182 '''
1183 Return the list of heights of the current data
1183 Return the list of heights of the current data
1184 '''
1184 '''
1185
1185
1186 return numpy.array(self.__heights[-1])
1186 return numpy.array(self.__heights[-1])
1187
1187
1188 @staticmethod
1188 @staticmethod
1189 def roundFloats(obj):
1189 def roundFloats(obj):
1190 if isinstance(obj, list):
1190 if isinstance(obj, list):
1191 return list(map(PlotterData.roundFloats, obj))
1191 return list(map(PlotterData.roundFloats, obj))
1192 elif isinstance(obj, float):
1192 elif isinstance(obj, float):
1193 return round(obj, 2)
1193 return round(obj, 2)
General Comments 0
You need to be logged in to leave comments. Login now