##// END OF EJS Templates
schain
RSS
Modifications for JULIA Processing
imanay -
r1764:66b4c9f36252
parent child
Show More
Side by Side Unified
Context file:
r1764:66b4c9f36252 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1182 +1,1187
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 # JULIA processing
546 self.data_diffcspc = 0
547 self.nDiffIncohInt = 0
548 # JULIA processing
546 self.useLocalTime = True
549 self.useLocalTime = True
547 self.radarControllerHeaderObj = RadarControllerHeader()
550 self.radarControllerHeaderObj = RadarControllerHeader()
548 self.systemHeaderObj = SystemHeader()
551 self.systemHeaderObj = SystemHeader()
549 self.type = "Spectra"
552 self.type = "Spectra"
550 self.timeZone = 0
553 self.timeZone = 0
551 self.nProfiles = None
554 self.nProfiles = None
552 self.heightList = None
555 self.heightList = None
553 self.channelList = None
556 self.channelList = None
554 self.pairsList = None
557 self.pairsList = None
555 self.flagNoData = True
558 self.flagNoData = True
556 self.flagDiscontinuousBlock = False
559 self.flagDiscontinuousBlock = False
557 self.utctime = None
560 self.utctime = None
558 self.nCohInt = None
561 self.nCohInt = None
559 self.nIncohInt = None
562 self.nIncohInt = None
560 self.blocksize = None
563 self.blocksize = None
561 self.nFFTPoints = None
564 self.nFFTPoints = None
562 self.wavelength = None
565 self.wavelength = None
563 self.flagDecodeData = False # asumo q la data no esta decodificada
566 self.flagDecodeData = False # asumo q la data no esta decodificada
564 self.flagDeflipData = False # asumo q la data no esta sin flip
567 self.flagDeflipData = False # asumo q la data no esta sin flip
565 self.flagShiftFFT = False
568 self.flagShiftFFT = False
566 self.ippFactor = 1
569 self.ippFactor = 1
567 self.beacon_heiIndexList = []
570 self.beacon_heiIndexList = []
568 self.noise_estimation = None
571 self.noise_estimation = None
569 self.spc_noise = None
572 self.spc_noise = None
570 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
573 self.metadata_list = ['type', 'heightList', 'timeZone', 'pairsList', 'channelList', 'nCohInt',
571 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp', 'nIncohInt', 'nFFTPoints', 'nProfiles', 'flagDecodeData']
574 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp', 'nIncohInt', 'nFFTPoints', 'nProfiles', 'flagDecodeData']
572
575
573
576
574 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
577 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
575 """
578 """
576 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
579 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
577
580
578 Return:
581 Return:
579 noiselevel
582 noiselevel
580 """
583 """
581
584
582 noise = numpy.zeros(self.nChannels)
585 noise = numpy.zeros(self.nChannels)
583
586
584 for channel in range(self.nChannels):
587 for channel in range(self.nChannels):
585 daux = self.data_spc[channel,
588 daux = self.data_spc[channel,
586 xmin_index:xmax_index, ymin_index:ymax_index]
589 xmin_index:xmax_index, ymin_index:ymax_index]
587 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
590 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
588
591
589 return noise
592 return noise
590
593
591 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
594 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
592
595
593 if self.spc_noise is not None:
596 if self.spc_noise is not None:
594 # this was estimated by getNoise Operation defined in jroproc_parameters.py
597 # this was estimated by getNoise Operation defined in jroproc_parameters.py
595 return self.spc_noise
598 return self.spc_noise
596 elif self.noise_estimation is not None:
599 elif self.noise_estimation is not None:
597 # this was estimated by getNoise Operation defined in jroproc_spectra.py
600 # this was estimated by getNoise Operation defined in jroproc_spectra.py
598 return self.noise_estimation
601 return self.noise_estimation
599 else:
602 else:
600 noise = self.getNoisebyHildebrand(
603 noise = self.getNoisebyHildebrand(
601 xmin_index, xmax_index, ymin_index, ymax_index)
604 xmin_index, xmax_index, ymin_index, ymax_index)
602 return noise
605 return noise
603
606
604 def getFreqRangeTimeResponse(self, extrapoints=0):
607 def getFreqRangeTimeResponse(self, extrapoints=0):
605
608
606 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
609 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
607 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
610 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
608
611
609 return freqrange
612 return freqrange
610
613
611 def getAcfRange(self, extrapoints=0):
614 def getAcfRange(self, extrapoints=0):
612
615
613 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
616 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
614 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
617 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
615
618
616 return freqrange
619 return freqrange
617
620
618 def getFreqRange(self, extrapoints=0):
621 def getFreqRange(self, extrapoints=0):
619
622
620 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
623 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
621 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
624 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
622
625
623 return freqrange
626 return freqrange
624
627
625 def getVelRange(self, extrapoints=0):
628 def getVelRange(self, extrapoints=0):
626
629
627 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
630 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
628 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
631 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
629
632
630 if self.nmodes:
633 if self.nmodes:
631 return velrange / self.nmodes
634 return velrange / self.nmodes
632 else:
635 else:
633 return velrange
636 return velrange
634
637
635 @property
638 @property
636 def nPairs(self):
639 def nPairs(self):
637
640
638 return len(self.pairsList)
641 return len(self.pairsList)
639
642
640 @property
643 @property
641 def pairsIndexList(self):
644 def pairsIndexList(self):
642
645
643 return list(range(self.nPairs))
646 return list(range(self.nPairs))
644
647
645 @property
648 @property
646 def normFactor(self):
649 def normFactor(self):
647
650
648 pwcode = 1
651 pwcode = 1
649
652
650 if self.flagDecodeData:
653 if self.flagDecodeData:
651 pwcode = numpy.sum(self.code[0] ** 2)
654 pwcode = numpy.sum(self.code[0] ** 2)
652 # normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
655 # normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
653 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
656 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
654
657
655 return normFactor
658 return normFactor
656
659
657 @property
660 @property
658 def flag_cspc(self):
661 def flag_cspc(self):
659
662
660 if self.data_cspc is None:
663 if self.data_cspc is None:
661 return True
664 return True
662
665
663 return False
666 return False
664
667
665 @property
668 @property
666 def flag_dc(self):
669 def flag_dc(self):
667
670
668 if self.data_dc is None:
671 if self.data_dc is None:
669 return True
672 return True
670
673
671 return False
674 return False
672
675
673 @property
676 @property
674 def timeInterval(self):
677 def timeInterval(self):
675
678
676 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
679 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
677 if self.nmodes:
680 if self.nmodes:
678 return self.nmodes * timeInterval
681 return self.nmodes * timeInterval
679 else:
682 else:
680 return timeInterval
683 return timeInterval
681
684
682 def getPower(self):
685 def getPower(self):
683
686
684 factor = self.normFactor
687 factor = self.normFactor
685 z = self.data_spc / factor
688 z = self.data_spc / factor
686 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
689 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
687 avg = numpy.average(z, axis=1)
690 avg = numpy.average(z, axis=1)
688
691
689 return 10 * numpy.log10(avg)
692 return 10 * numpy.log10(avg)
690
693
691 def getCoherence(self, pairsList=None, phase=False):
694 def getCoherence(self, pairsList=None, phase=False):
692
695
693 z = []
696 z = []
694 if pairsList is None:
697 if pairsList is None:
695 pairsIndexList = self.pairsIndexList
698 pairsIndexList = self.pairsIndexList
696 else:
699 else:
697 pairsIndexList = []
700 pairsIndexList = []
698 for pair in pairsList:
701 for pair in pairsList:
699 if pair not in self.pairsList:
702 if pair not in self.pairsList:
700 raise ValueError("Pair %s is not in dataOut.pairsList" % (
703 raise ValueError("Pair %s is not in dataOut.pairsList" % (
701 pair))
704 pair))
702 pairsIndexList.append(self.pairsList.index(pair))
705 pairsIndexList.append(self.pairsList.index(pair))
703 for i in range(len(pairsIndexList)):
706 for i in range(len(pairsIndexList)):
704 pair = self.pairsList[pairsIndexList[i]]
707 pair = self.pairsList[pairsIndexList[i]]
705 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
708 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
706 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
709 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
707 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
710 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
708 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
711 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
709 if phase:
712 if phase:
710 data = numpy.arctan2(avgcoherenceComplex.imag,
713 data = numpy.arctan2(avgcoherenceComplex.imag,
711 avgcoherenceComplex.real) * 180 / numpy.pi
714 avgcoherenceComplex.real) * 180 / numpy.pi
712 else:
715 else:
713 data = numpy.abs(avgcoherenceComplex)
716 data = numpy.abs(avgcoherenceComplex)
714
717
715 z.append(data)
718 z.append(data)
716
719
717 return numpy.array(z)
720 return numpy.array(z)
718
721
719 def setValue(self, value):
722 def setValue(self, value):
720
723
721 print("This property should not be initialized")
724 print("This property should not be initialized")
722
725
723 return
726 return
724
727
725 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
728 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
726
729
727
730
728 class SpectraHeis(Spectra):
731 class SpectraHeis(Spectra):
729
732
730 def __init__(self):
733 def __init__(self):
731
734
732 self.radarControllerHeaderObj = RadarControllerHeader()
735 self.radarControllerHeaderObj = RadarControllerHeader()
733 self.systemHeaderObj = SystemHeader()
736 self.systemHeaderObj = SystemHeader()
734 self.type = "SpectraHeis"
737 self.type = "SpectraHeis"
735 self.nProfiles = None
738 self.nProfiles = None
736 self.heightList = None
739 self.heightList = None
737 self.channelList = None
740 self.channelList = None
738 self.flagNoData = True
741 self.flagNoData = True
739 self.flagDiscontinuousBlock = False
742 self.flagDiscontinuousBlock = False
740 self.utctime = None
743 self.utctime = None
741 self.blocksize = None
744 self.blocksize = None
742 self.profileIndex = 0
745 self.profileIndex = 0
743 self.nCohInt = 1
746 self.nCohInt = 1
744 self.nIncohInt = 1
747 self.nIncohInt = 1
745
748
746 @property
749 @property
747 def normFactor(self):
750 def normFactor(self):
748 pwcode = 1
751 pwcode = 1
749 if self.flagDecodeData:
752 if self.flagDecodeData:
750 pwcode = numpy.sum(self.code[0] ** 2)
753 pwcode = numpy.sum(self.code[0] ** 2)
751
754
752 normFactor = self.nIncohInt * self.nCohInt * pwcode
755 normFactor = self.nIncohInt * self.nCohInt * pwcode
753
756
754 return normFactor
757 return normFactor
755
758
756 @property
759 @property
757 def timeInterval(self):
760 def timeInterval(self):
758
761
759 return self.ippSeconds * self.nCohInt * self.nIncohInt
762 return self.ippSeconds * self.nCohInt * self.nIncohInt
760
763
761
764
762 class Fits(JROData):
765 class Fits(JROData):
763
766
764 def __init__(self):
767 def __init__(self):
765
768
766 self.type = "Fits"
769 self.type = "Fits"
767 self.nProfiles = None
770 self.nProfiles = None
768 self.heightList = None
771 self.heightList = None
769 self.channelList = None
772 self.channelList = None
770 self.flagNoData = True
773 self.flagNoData = True
771 self.utctime = None
774 self.utctime = None
772 self.nCohInt = 1
775 self.nCohInt = 1
773 self.nIncohInt = 1
776 self.nIncohInt = 1
774 self.useLocalTime = True
777 self.useLocalTime = True
775 self.profileIndex = 0
778 self.profileIndex = 0
776 self.timeZone = 0
779 self.timeZone = 0
777
780
778 def getTimeRange(self):
781 def getTimeRange(self):
779
782
780 datatime = []
783 datatime = []
781
784
782 datatime.append(self.ltctime)
785 datatime.append(self.ltctime)
783 datatime.append(self.ltctime + self.timeInterval)
786 datatime.append(self.ltctime + self.timeInterval)
784
787
785 datatime = numpy.array(datatime)
788 datatime = numpy.array(datatime)
786
789
787 return datatime
790 return datatime
788
791
789 def getChannelIndexList(self):
792 def getChannelIndexList(self):
790
793
791 return list(range(self.nChannels))
794 return list(range(self.nChannels))
792
795
793 def getNoise(self, type=1):
796 def getNoise(self, type=1):
794
797
795
798
796 if type == 1:
799 if type == 1:
797 noise = self.getNoisebyHildebrand()
800 noise = self.getNoisebyHildebrand()
798
801
799 if type == 2:
802 if type == 2:
800 noise = self.getNoisebySort()
803 noise = self.getNoisebySort()
801
804
802 if type == 3:
805 if type == 3:
803 noise = self.getNoisebyWindow()
806 noise = self.getNoisebyWindow()
804
807
805 return noise
808 return noise
806
809
807 @property
810 @property
808 def timeInterval(self):
811 def timeInterval(self):
809
812
810 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
813 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
811
814
812 return timeInterval
815 return timeInterval
813
816
814 @property
817 @property
815 def ippSeconds(self):
818 def ippSeconds(self):
816 '''
819 '''
817 '''
820 '''
818 return self.ipp_sec
821 return self.ipp_sec
819
822
820 noise = property(getNoise, "I'm the 'nHeights' property.")
823 noise = property(getNoise, "I'm the 'nHeights' property.")
821
824
822
825
823 class Correlation(JROData):
826 class Correlation(JROData):
824
827
825 def __init__(self):
828 def __init__(self):
826 '''
829 '''
827 Constructor
830 Constructor
828 '''
831 '''
829 self.radarControllerHeaderObj = RadarControllerHeader()
832 self.radarControllerHeaderObj = RadarControllerHeader()
830 self.systemHeaderObj = SystemHeader()
833 self.systemHeaderObj = SystemHeader()
831 self.type = "Correlation"
834 self.type = "Correlation"
832 self.data = None
835 self.data = None
833 self.dtype = None
836 self.dtype = None
834 self.nProfiles = None
837 self.nProfiles = None
835 self.heightList = None
838 self.heightList = None
836 self.channelList = None
839 self.channelList = None
837 self.flagNoData = True
840 self.flagNoData = True
838 self.flagDiscontinuousBlock = False
841 self.flagDiscontinuousBlock = False
839 self.utctime = None
842 self.utctime = None
840 self.timeZone = 0
843 self.timeZone = 0
841 self.dstFlag = None
844 self.dstFlag = None
842 self.errorCount = None
845 self.errorCount = None
843 self.blocksize = None
846 self.blocksize = None
844 self.flagDecodeData = False # asumo q la data no esta decodificada
847 self.flagDecodeData = False # asumo q la data no esta decodificada
845 self.flagDeflipData = False # asumo q la data no esta sin flip
848 self.flagDeflipData = False # asumo q la data no esta sin flip
846 self.pairsList = None
849 self.pairsList = None
847 self.nPoints = None
850 self.nPoints = None
848
851
849 def getPairsList(self):
852 def getPairsList(self):
850
853
851 return self.pairsList
854 return self.pairsList
852
855
853 def getNoise(self, mode=2):
856 def getNoise(self, mode=2):
854
857
855 indR = numpy.where(self.lagR == 0)[0][0]
858 indR = numpy.where(self.lagR == 0)[0][0]
856 indT = numpy.where(self.lagT == 0)[0][0]
859 indT = numpy.where(self.lagT == 0)[0][0]
857
860
858 jspectra0 = self.data_corr[:, :, indR, :]
861 jspectra0 = self.data_corr[:, :, indR, :]
859 jspectra = copy.copy(jspectra0)
862 jspectra = copy.copy(jspectra0)
860
863
861 num_chan = jspectra.shape[0]
864 num_chan = jspectra.shape[0]
862 num_hei = jspectra.shape[2]
865 num_hei = jspectra.shape[2]
863
866
864 freq_dc = jspectra.shape[1] / 2
867 freq_dc = jspectra.shape[1] / 2
865 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
868 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
866
869
867 if ind_vel[0] < 0:
870 if ind_vel[0] < 0:
868 ind_vel[list(range(0, 1))] = ind_vel[list(
871 ind_vel[list(range(0, 1))] = ind_vel[list(
869 range(0, 1))] + self.num_prof
872 range(0, 1))] + self.num_prof
870
873
871 if mode == 1:
874 if mode == 1:
872 jspectra[:, freq_dc, :] = (
875 jspectra[:, freq_dc, :] = (
873 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
876 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
874
877
875 if mode == 2:
878 if mode == 2:
876
879
877 vel = numpy.array([-2, -1, 1, 2])
880 vel = numpy.array([-2, -1, 1, 2])
878 xx = numpy.zeros([4, 4])
881 xx = numpy.zeros([4, 4])
879
882
880 for fil in range(4):
883 for fil in range(4):
881 xx[fil, :] = vel[fil] ** numpy.asarray(list(range(4)))
884 xx[fil, :] = vel[fil] ** numpy.asarray(list(range(4)))
882
885
883 xx_inv = numpy.linalg.inv(xx)
886 xx_inv = numpy.linalg.inv(xx)
884 xx_aux = xx_inv[0, :]
887 xx_aux = xx_inv[0, :]
885
888
886 for ich in range(num_chan):
889 for ich in range(num_chan):
887 yy = jspectra[ich, ind_vel, :]
890 yy = jspectra[ich, ind_vel, :]
888 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
891 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
889
892
890 junkid = jspectra[ich, freq_dc, :] <= 0
893 junkid = jspectra[ich, freq_dc, :] <= 0
891 cjunkid = sum(junkid)
894 cjunkid = sum(junkid)
892
895
893 if cjunkid.any():
896 if cjunkid.any():
894 jspectra[ich, freq_dc, junkid.nonzero()] = (
897 jspectra[ich, freq_dc, junkid.nonzero()] = (
895 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
898 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
896
899
897 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
900 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
898
901
899 return noise
902 return noise
900
903
901 @property
904 @property
902 def timeInterval(self):
905 def timeInterval(self):
903
906
904 return self.ippSeconds * self.nCohInt * self.nProfiles
907 return self.ippSeconds * self.nCohInt * self.nProfiles
905
908
906 def splitFunctions(self):
909 def splitFunctions(self):
907
910
908 pairsList = self.pairsList
911 pairsList = self.pairsList
909 ccf_pairs = []
912 ccf_pairs = []
910 acf_pairs = []
913 acf_pairs = []
911 ccf_ind = []
914 ccf_ind = []
912 acf_ind = []
915 acf_ind = []
913 for l in range(len(pairsList)):
916 for l in range(len(pairsList)):
914 chan0 = pairsList[l][0]
917 chan0 = pairsList[l][0]
915 chan1 = pairsList[l][1]
918 chan1 = pairsList[l][1]
916
919
917 # Obteniendo pares de Autocorrelacion
920 # Obteniendo pares de Autocorrelacion
918 if chan0 == chan1:
921 if chan0 == chan1:
919 acf_pairs.append(chan0)
922 acf_pairs.append(chan0)
920 acf_ind.append(l)
923 acf_ind.append(l)
921 else:
924 else:
922 ccf_pairs.append(pairsList[l])
925 ccf_pairs.append(pairsList[l])
923 ccf_ind.append(l)
926 ccf_ind.append(l)
924
927
925 data_acf = self.data_cf[acf_ind]
928 data_acf = self.data_cf[acf_ind]
926 data_ccf = self.data_cf[ccf_ind]
929 data_ccf = self.data_cf[ccf_ind]
927
930
928 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
931 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
929
932
930 @property
933 @property
931 def normFactor(self):
934 def normFactor(self):
932 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
935 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
933 acf_pairs = numpy.array(acf_pairs)
936 acf_pairs = numpy.array(acf_pairs)
934 normFactor = numpy.zeros((self.nPairs, self.nHeights))
937 normFactor = numpy.zeros((self.nPairs, self.nHeights))
935
938
936 for p in range(self.nPairs):
939 for p in range(self.nPairs):
937 pair = self.pairsList[p]
940 pair = self.pairsList[p]
938
941
939 ch0 = pair[0]
942 ch0 = pair[0]
940 ch1 = pair[1]
943 ch1 = pair[1]
941
944
942 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
945 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
943 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
946 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
944 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
947 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
945
948
946 return normFactor
949 return normFactor
947
950
948
951
949 class Parameters(Spectra):
952 class Parameters(Spectra):
950
953
951 groupList = None # List of Pairs, Groups, etc
954 groupList = None # List of Pairs, Groups, etc
952 data_param = None # Parameters obtained
955 data_param = None # Parameters obtained
953 data_pre = None # Data Pre Parametrization
956 data_pre = None # Data Pre Parametrization
954 data_SNR = None # Signal to Noise Ratio
957 data_SNR = None # Signal to Noise Ratio
955 abscissaList = None # Abscissa, can be velocities, lags or time
958 abscissaList = None # Abscissa, can be velocities, lags or time
956 utctimeInit = None # Initial UTC time
959 utctimeInit = None # Initial UTC time
957 paramInterval = None # Time interval to calculate Parameters in seconds
960 paramInterval = None # Time interval to calculate Parameters in seconds
958 useLocalTime = True
961 useLocalTime = True
959 # Fitting
962 # Fitting
960 data_error = None # Error of the estimation
963 data_error = None # Error of the estimation
961 constants = None
964 constants = None
962 library = None
965 library = None
963 # Output signal
966 # Output signal
964 outputInterval = None # Time interval to calculate output signal in seconds
967 outputInterval = None # Time interval to calculate output signal in seconds
965 data_output = None # Out signal
968 data_output = None # Out signal
966 nAvg = None
969 nAvg = None
967 noise_estimation = None
970 noise_estimation = None
968 GauSPC = None # Fit gaussian SPC
971 GauSPC = None # Fit gaussian SPC
969 spc_noise = None
972 spc_noise = None
973 avg_output = None # for 150Km processing
970 #
974 #
971
975
972
976
973 def __init__(self):
977 def __init__(self):
974 '''
978 '''
975 Constructor
979 Constructor
976 '''
980 '''
977 self.radarControllerHeaderObj = RadarControllerHeader()
981 self.radarControllerHeaderObj = RadarControllerHeader()
978 self.systemHeaderObj = SystemHeader()
982 self.systemHeaderObj = SystemHeader()
979 self.type = "Parameters"
983 self.type = "Parameters"
980 self.timeZone = 0
984 self.timeZone = 0
981 self.ippFactor = 1
985 self.ippFactor = 1
982 # JULIA processing
986 # JULIA processing
983 self.diffcspectra = None
987 self.diffcspectra = None
988 self.nDiffIncohInt = None
984 self.ccfpar = None
989 self.ccfpar = None
985 # JULIA processing
990 # JULIA processing
986
991
987 def getTimeRange1(self, interval):
992 def getTimeRange1(self, interval):
988
993
989 datatime = []
994 datatime = []
990
995
991 if self.useLocalTime:
996 if self.useLocalTime:
992 time1 = self.utctimeInit - self.timeZone * 60
997 time1 = self.utctimeInit - self.timeZone * 60
993 else:
998 else:
994 time1 = self.utctimeInit
999 time1 = self.utctimeInit
995
1000
996 datatime.append(time1)
1001 datatime.append(time1)
997 datatime.append(time1 + interval)
1002 datatime.append(time1 + interval)
998 datatime = numpy.array(datatime)
1003 datatime = numpy.array(datatime)
999
1004
1000 return datatime
1005 return datatime
1001
1006
1002 @property
1007 @property
1003 def timeInterval(self):
1008 def timeInterval(self):
1004
1009
1005 if hasattr(self, 'timeInterval1'):
1010 if hasattr(self, 'timeInterval1'):
1006 return self.timeInterval1
1011 return self.timeInterval1
1007 else:
1012 else:
1008 return self.paramInterval
1013 return self.paramInterval
1009
1014
1010
1015
1011 def setValue(self, value):
1016 def setValue(self, value):
1012
1017
1013 print("This property should not be initialized")
1018 print("This property should not be initialized")
1014
1019
1015 return
1020 return
1016
1021
1017
1022
1018 class PlotterData(object):
1023 class PlotterData(object):
1019 '''
1024 '''
1020 Object to hold data to be plotted
1025 Object to hold data to be plotted
1021 '''
1026 '''
1022
1027
1023 MAXNUMX = 200
1028 MAXNUMX = 200
1024 MAXNUMY = 200
1029 MAXNUMY = 200
1025
1030
1026 def __init__(self, code, exp_code, localtime=True):
1031 def __init__(self, code, exp_code, localtime=True):
1027
1032
1028 self.key = code
1033 self.key = code
1029 self.exp_code = exp_code
1034 self.exp_code = exp_code
1030 self.ready = False
1035 self.ready = False
1031 self.flagNoData = False
1036 self.flagNoData = False
1032 self.localtime = localtime
1037 self.localtime = localtime
1033 self.data = {}
1038 self.data = {}
1034 self.meta = {}
1039 self.meta = {}
1035 self.__heights = []
1040 self.__heights = []
1036
1041
1037 def __str__(self):
1042 def __str__(self):
1038 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1043 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1039 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
1044 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
1040
1045
1041 def __len__(self):
1046 def __len__(self):
1042 return len(self.data)
1047 return len(self.data)
1043
1048
1044 def __getitem__(self, key):
1049 def __getitem__(self, key):
1045 if isinstance(key, int):
1050 if isinstance(key, int):
1046 return self.data[self.times[key]]
1051 return self.data[self.times[key]]
1047 elif isinstance(key, str):
1052 elif isinstance(key, str):
1048 ret = numpy.array([self.data[x][key] for x in self.times])
1053 ret = numpy.array([self.data[x][key] for x in self.times])
1049 if ret.ndim > 1:
1054 if ret.ndim > 1:
1050 ret = numpy.swapaxes(ret, 0, 1)
1055 ret = numpy.swapaxes(ret, 0, 1)
1051 return ret
1056 return ret
1052
1057
1053 def __contains__(self, key):
1058 def __contains__(self, key):
1054 return key in self.data[self.min_time]
1059 return key in self.data[self.min_time]
1055
1060
1056 def setup(self):
1061 def setup(self):
1057 '''
1062 '''
1058 Configure object
1063 Configure object
1059 '''
1064 '''
1060 self.type = ''
1065 self.type = ''
1061 self.ready = False
1066 self.ready = False
1062 del self.data
1067 del self.data
1063 self.data = {}
1068 self.data = {}
1064 self.__heights = []
1069 self.__heights = []
1065 self.__all_heights = set()
1070 self.__all_heights = set()
1066
1071
1067 def shape(self, key):
1072 def shape(self, key):
1068 '''
1073 '''
1069 Get the shape of the one-element data for the given key
1074 Get the shape of the one-element data for the given key
1070 '''
1075 '''
1071
1076
1072 if len(self.data[self.min_time][key]):
1077 if len(self.data[self.min_time][key]):
1073 return self.data[self.min_time][key].shape
1078 return self.data[self.min_time][key].shape
1074 return (0,)
1079 return (0,)
1075
1080
1076 def update(self, data, tm, meta={}):
1081 def update(self, data, tm, meta={}):
1077 '''
1082 '''
1078 Update data object with new dataOut
1083 Update data object with new dataOut
1079 '''
1084 '''
1080
1085
1081 self.data[tm] = data
1086 self.data[tm] = data
1082
1087
1083 for key, value in meta.items():
1088 for key, value in meta.items():
1084 setattr(self, key, value)
1089 setattr(self, key, value)
1085
1090
1086 def normalize_heights(self):
1091 def normalize_heights(self):
1087 '''
1092 '''
1088 Ensure same-dimension of the data for different heighList
1093 Ensure same-dimension of the data for different heighList
1089 '''
1094 '''
1090
1095
1091 H = numpy.array(list(self.__all_heights))
1096 H = numpy.array(list(self.__all_heights))
1092 H.sort()
1097 H.sort()
1093 for key in self.data:
1098 for key in self.data:
1094 shape = self.shape(key)[:-1] + H.shape
1099 shape = self.shape(key)[:-1] + H.shape
1095 for tm, obj in list(self.data[key].items()):
1100 for tm, obj in list(self.data[key].items()):
1096 h = self.__heights[self.times.tolist().index(tm)]
1101 h = self.__heights[self.times.tolist().index(tm)]
1097 if H.size == h.size:
1102 if H.size == h.size:
1098 continue
1103 continue
1099 index = numpy.where(numpy.in1d(H, h))[0]
1104 index = numpy.where(numpy.in1d(H, h))[0]
1100 dummy = numpy.zeros(shape) + numpy.nan
1105 dummy = numpy.zeros(shape) + numpy.nan
1101 if len(shape) == 2:
1106 if len(shape) == 2:
1102 dummy[:, index] = obj
1107 dummy[:, index] = obj
1103 else:
1108 else:
1104 dummy[index] = obj
1109 dummy[index] = obj
1105 self.data[key][tm] = dummy
1110 self.data[key][tm] = dummy
1106
1111
1107 self.__heights = [H for tm in self.times]
1112 self.__heights = [H for tm in self.times]
1108
1113
1109 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1114 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1110 '''
1115 '''
1111 Convert data to json
1116 Convert data to json
1112 '''
1117 '''
1113
1118
1114 meta = {}
1119 meta = {}
1115 meta['xrange'] = []
1120 meta['xrange'] = []
1116 dy = int(len(self.yrange) / self.MAXNUMY) + 1
1121 dy = int(len(self.yrange) / self.MAXNUMY) + 1
1117 tmp = self.data[tm][self.key]
1122 tmp = self.data[tm][self.key]
1118 shape = tmp.shape
1123 shape = tmp.shape
1119 if len(shape) == 2:
1124 if len(shape) == 2:
1120 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1125 data = self.roundFloats(self.data[tm][self.key][::, ::dy].tolist())
1121 elif len(shape) == 3:
1126 elif len(shape) == 3:
1122 dx = int(self.data[tm][self.key].shape[1] / self.MAXNUMX) + 1
1127 dx = int(self.data[tm][self.key].shape[1] / self.MAXNUMX) + 1
1123 data = self.roundFloats(
1128 data = self.roundFloats(
1124 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1129 self.data[tm][self.key][::, ::dx, ::dy].tolist())
1125 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1130 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1126 else:
1131 else:
1127 data = self.roundFloats(self.data[tm][self.key].tolist())
1132 data = self.roundFloats(self.data[tm][self.key].tolist())
1128
1133
1129 ret = {
1134 ret = {
1130 'plot': plot_name,
1135 'plot': plot_name,
1131 'code': self.exp_code,
1136 'code': self.exp_code,
1132 'time': float(tm),
1137 'time': float(tm),
1133 'data': data,
1138 'data': data,
1134 }
1139 }
1135 meta['type'] = plot_type
1140 meta['type'] = plot_type
1136 meta['interval'] = float(self.interval)
1141 meta['interval'] = float(self.interval)
1137 meta['localtime'] = self.localtime
1142 meta['localtime'] = self.localtime
1138 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1143 meta['yrange'] = self.roundFloats(self.yrange[::dy].tolist())
1139 meta.update(self.meta)
1144 meta.update(self.meta)
1140 ret['metadata'] = meta
1145 ret['metadata'] = meta
1141 return json.dumps(ret)
1146 return json.dumps(ret)
1142
1147
1143 @property
1148 @property
1144 def times(self):
1149 def times(self):
1145 '''
1150 '''
1146 Return the list of times of the current data
1151 Return the list of times of the current data
1147 '''
1152 '''
1148
1153
1149 ret = [t for t in self.data]
1154 ret = [t for t in self.data]
1150 ret.sort()
1155 ret.sort()
1151 return numpy.array(ret)
1156 return numpy.array(ret)
1152
1157
1153 @property
1158 @property
1154 def min_time(self):
1159 def min_time(self):
1155 '''
1160 '''
1156 Return the minimun time value
1161 Return the minimun time value
1157 '''
1162 '''
1158
1163
1159 return self.times[0]
1164 return self.times[0]
1160
1165
1161 @property
1166 @property
1162 def max_time(self):
1167 def max_time(self):
1163 '''
1168 '''
1164 Return the maximun time value
1169 Return the maximun time value
1165 '''
1170 '''
1166
1171
1167 return self.times[-1]
1172 return self.times[-1]
1168
1173
1169 # @property
1174 # @property
1170 # def heights(self):
1175 # def heights(self):
1171 # '''
1176 # '''
1172 # Return the list of heights of the current data
1177 # Return the list of heights of the current data
1173 # '''
1178 # '''
1174
1179
1175 # return numpy.array(self.__heights[-1])
1180 # return numpy.array(self.__heights[-1])
1176
1181
1177 @staticmethod
1182 @staticmethod
1178 def roundFloats(obj):
1183 def roundFloats(obj):
1179 if isinstance(obj, list):
1184 if isinstance(obj, list):
1180 return list(map(PlotterData.roundFloats, obj))
1185 return list(map(PlotterData.roundFloats, obj))
1181 elif isinstance(obj, float):
1186 elif isinstance(obj, float):
1182 return round(obj, 2)
1187 return round(obj, 2)
General Comments 0
You need to be logged in to leave comments. Login now