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