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