##// END OF EJS Templates
schain
RSS
Update noise C extension to properly work with python 3
Juan C. Espinoza -
r1286:a59b24777aee
parent child
Show More
Side by Side Unified
Context file:
r1286:a59b24777aee Loading diff...:
Show file before | Show file after | Hide comments
@@ -0,0 +1,82
1 #include <Python.h>
2 #include <numpy/arrayobject.h>
3 #include <math.h>
4
5
6 static PyObject *hildebrand_sekhon(PyObject *self, PyObject *args) {
7 double navg;
8 PyObject *data_obj, *data_array;
9
10 if (!PyArg_ParseTuple(args, "Od", &data_obj, &navg)) {
11 return NULL;
12 }
13
14 data_array = PyArray_FROM_OTF(data_obj, NPY_FLOAT64, NPY_IN_ARRAY);
15
16 if (data_array == NULL) {
17 Py_XDECREF(data_array);
18 Py_XDECREF(data_obj);
19 return NULL;
20 }
21 double *sortdata = (double*)PyArray_DATA(data_array);
22 int lenOfData = (int)PyArray_SIZE(data_array) ;
23 double nums_min = lenOfData*0.2;
24 if (nums_min <= 5) nums_min = 5;
25 double sump = 0;
26 double sumq = 0;
27 int j = 0;
28 int cont = 1;
29 double rtest = 0;
30 while ((cont == 1) && (j < lenOfData)) {
31 sump = sump + sortdata[j];
32 sumq = sumq + pow(sortdata[j], 2);
33 if (j > nums_min) {
34 rtest = (double)j/(j-1) + 1/navg;
35 if ((sumq*j) > (rtest*pow(sump, 2))) {
36 j = j - 1;
37 sump = sump - sortdata[j];
38 sumq = sumq - pow(sortdata[j],2);
39 cont = 0;
40 }
41 }
42 j = j + 1;
43 }
44
45 double lnoise = sump / j;
46
47 Py_DECREF(data_array);
48
49 return PyLong_FromLong(lnoise);
50 //return Py_BuildValue("d", lnoise);
51 }
52
53
54 static PyMethodDef noiseMethods[] = {
55 { "hildebrand_sekhon", hildebrand_sekhon, METH_VARARGS, "Get noise with hildebrand_sekhon algorithm" },
56 { NULL, NULL, 0, NULL }
57 };
58
59 #if PY_MAJOR_VERSION >= 3
60
61 static struct PyModuleDef noisemodule = {
62 PyModuleDef_HEAD_INIT,
63 "_noise",
64 "Get noise with hildebrand_sekhon algorithm",
65 -1,
66 noiseMethods
67 };
68
69 #endif
70
71 #if PY_MAJOR_VERSION >= 3
72 PyMODINIT_FUNC PyInit__noise(void) {
73 Py_Initialize();
74 import_array();
75 return PyModule_Create(&noisemodule);
76 }
77 #else
78 PyMODINIT_FUNC init_noise() {
79 Py_InitModule("_noise", noiseMethods);
80 import_array();
81 }
82 #endif
Show file before | Show file after | Hide comments
@@ -1,1384 +1,1386
1 '''
1 '''
2
2
3 $Author: murco $
3 $Author: murco $
4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 '''
5 '''
6
6
7 import copy
7 import copy
8 import numpy
8 import numpy
9 import datetime
9 import datetime
10 import json
10 import json
11
11
12 import schainpy.admin
12 import schainpy.admin
13 from schainpy.utils import log
13 from schainpy.utils import log
14 from .jroheaderIO import SystemHeader, RadarControllerHeader
14 from .jroheaderIO import SystemHeader, RadarControllerHeader
15 from schainpy.model.data import _noise
15
16
16
17
17 def getNumpyDtype(dataTypeCode):
18 def getNumpyDtype(dataTypeCode):
18
19
19 if dataTypeCode == 0:
20 if dataTypeCode == 0:
20 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
21 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
21 elif dataTypeCode == 1:
22 elif dataTypeCode == 1:
22 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
23 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
23 elif dataTypeCode == 2:
24 elif dataTypeCode == 2:
24 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
25 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
25 elif dataTypeCode == 3:
26 elif dataTypeCode == 3:
26 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
27 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
27 elif dataTypeCode == 4:
28 elif dataTypeCode == 4:
28 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
29 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
29 elif dataTypeCode == 5:
30 elif dataTypeCode == 5:
30 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
31 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
31 else:
32 else:
32 raise ValueError('dataTypeCode was not defined')
33 raise ValueError('dataTypeCode was not defined')
33
34
34 return numpyDtype
35 return numpyDtype
35
36
36
37
37 def getDataTypeCode(numpyDtype):
38 def getDataTypeCode(numpyDtype):
38
39
39 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
40 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
40 datatype = 0
41 datatype = 0
41 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
42 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
42 datatype = 1
43 datatype = 1
43 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
44 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
44 datatype = 2
45 datatype = 2
45 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
46 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
46 datatype = 3
47 datatype = 3
47 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
48 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
48 datatype = 4
49 datatype = 4
49 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
50 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
50 datatype = 5
51 datatype = 5
51 else:
52 else:
52 datatype = None
53 datatype = None
53
54
54 return datatype
55 return datatype
55
56
56
57
57 def hildebrand_sekhon(data, navg):
58 def hildebrand_sekhon(data, navg):
58 """
59 """
59 This method is for the objective determination of the noise level in Doppler spectra. This
60 This method is for the objective determination of the noise level in Doppler spectra. This
60 implementation technique is based on the fact that the standard deviation of the spectral
61 implementation technique is based on the fact that the standard deviation of the spectral
61 densities is equal to the mean spectral density for white Gaussian noise
62 densities is equal to the mean spectral density for white Gaussian noise
62
63
63 Inputs:
64 Inputs:
64 Data : heights
65 Data : heights
65 navg : numbers of averages
66 navg : numbers of averages
66
67
67 Return:
68 Return:
68 mean : noise's level
69 mean : noise's level
69 """
70 """
70
71
71 sortdata = numpy.sort(data, axis=None)
72 sortdata = numpy.sort(data, axis=None)
73 '''
72 lenOfData = len(sortdata)
74 lenOfData = len(sortdata)
73 nums_min = lenOfData*0.2
75 nums_min = lenOfData*0.2
74
76
75 if nums_min <= 5:
77 if nums_min <= 5:
76
78
77 nums_min = 5
79 nums_min = 5
78
80
79 sump = 0.
81 sump = 0.
80 sumq = 0.
82 sumq = 0.
81
83
82 j = 0
84 j = 0
83 cont = 1
85 cont = 1
84
86
85 while((cont == 1)and(j < lenOfData)):
87 while((cont == 1)and(j < lenOfData)):
86
88
87 sump += sortdata[j]
89 sump += sortdata[j]
88 sumq += sortdata[j]**2
90 sumq += sortdata[j]**2
89
91
90 if j > nums_min:
92 if j > nums_min:
91 rtest = float(j)/(j-1) + 1.0/navg
93 rtest = float(j)/(j-1) + 1.0/navg
92 if ((sumq*j) > (rtest*sump**2)):
94 if ((sumq*j) > (rtest*sump**2)):
93 j = j - 1
95 j = j - 1
94 sump = sump - sortdata[j]
96 sump = sump - sortdata[j]
95 sumq = sumq - sortdata[j]**2
97 sumq = sumq - sortdata[j]**2
96 cont = 0
98 cont = 0
97
99
98 j += 1
100 j += 1
99
101
100 lnoise = sump / j
102 lnoise = sump / j
101
103 '''
102 return lnoise
104 return _noise.hildebrand_sekhon(sortdata, navg)
103
105
104
106
105 class Beam:
107 class Beam:
106
108
107 def __init__(self):
109 def __init__(self):
108 self.codeList = []
110 self.codeList = []
109 self.azimuthList = []
111 self.azimuthList = []
110 self.zenithList = []
112 self.zenithList = []
111
113
112
114
113 class GenericData(object):
115 class GenericData(object):
114
116
115 flagNoData = True
117 flagNoData = True
116
118
117 def copy(self, inputObj=None):
119 def copy(self, inputObj=None):
118
120
119 if inputObj == None:
121 if inputObj == None:
120 return copy.deepcopy(self)
122 return copy.deepcopy(self)
121
123
122 for key in list(inputObj.__dict__.keys()):
124 for key in list(inputObj.__dict__.keys()):
123
125
124 attribute = inputObj.__dict__[key]
126 attribute = inputObj.__dict__[key]
125
127
126 # If this attribute is a tuple or list
128 # If this attribute is a tuple or list
127 if type(inputObj.__dict__[key]) in (tuple, list):
129 if type(inputObj.__dict__[key]) in (tuple, list):
128 self.__dict__[key] = attribute[:]
130 self.__dict__[key] = attribute[:]
129 continue
131 continue
130
132
131 # If this attribute is another object or instance
133 # If this attribute is another object or instance
132 if hasattr(attribute, '__dict__'):
134 if hasattr(attribute, '__dict__'):
133 self.__dict__[key] = attribute.copy()
135 self.__dict__[key] = attribute.copy()
134 continue
136 continue
135
137
136 self.__dict__[key] = inputObj.__dict__[key]
138 self.__dict__[key] = inputObj.__dict__[key]
137
139
138 def deepcopy(self):
140 def deepcopy(self):
139
141
140 return copy.deepcopy(self)
142 return copy.deepcopy(self)
141
143
142 def isEmpty(self):
144 def isEmpty(self):
143
145
144 return self.flagNoData
146 return self.flagNoData
145
147
146 def isReady(self):
148 def isReady(self):
147
149
148 return not self.flagNoData
150 return not self.flagNoData
149
151
150
152
151 class JROData(GenericData):
153 class JROData(GenericData):
152
154
153 # m_BasicHeader = BasicHeader()
155 # m_BasicHeader = BasicHeader()
154 # m_ProcessingHeader = ProcessingHeader()
156 # m_ProcessingHeader = ProcessingHeader()
155
157
156 systemHeaderObj = SystemHeader()
158 systemHeaderObj = SystemHeader()
157 radarControllerHeaderObj = RadarControllerHeader()
159 radarControllerHeaderObj = RadarControllerHeader()
158 # data = None
160 # data = None
159 type = None
161 type = None
160 datatype = None # dtype but in string
162 datatype = None # dtype but in string
161 # dtype = None
163 # dtype = None
162 # nChannels = None
164 # nChannels = None
163 # nHeights = None
165 # nHeights = None
164 nProfiles = None
166 nProfiles = None
165 heightList = None
167 heightList = None
166 channelList = None
168 channelList = None
167 flagDiscontinuousBlock = False
169 flagDiscontinuousBlock = False
168 useLocalTime = False
170 useLocalTime = False
169 utctime = None
171 utctime = None
170 timeZone = None
172 timeZone = None
171 dstFlag = None
173 dstFlag = None
172 errorCount = None
174 errorCount = None
173 blocksize = None
175 blocksize = None
174 # nCode = None
176 # nCode = None
175 # nBaud = None
177 # nBaud = None
176 # code = None
178 # code = None
177 flagDecodeData = False # asumo q la data no esta decodificada
179 flagDecodeData = False # asumo q la data no esta decodificada
178 flagDeflipData = False # asumo q la data no esta sin flip
180 flagDeflipData = False # asumo q la data no esta sin flip
179 flagShiftFFT = False
181 flagShiftFFT = False
180 # ippSeconds = None
182 # ippSeconds = None
181 # timeInterval = None
183 # timeInterval = None
182 nCohInt = None
184 nCohInt = None
183 # noise = None
185 # noise = None
184 windowOfFilter = 1
186 windowOfFilter = 1
185 # Speed of ligth
187 # Speed of ligth
186 C = 3e8
188 C = 3e8
187 frequency = 49.92e6
189 frequency = 49.92e6
188 realtime = False
190 realtime = False
189 beacon_heiIndexList = None
191 beacon_heiIndexList = None
190 last_block = None
192 last_block = None
191 blocknow = None
193 blocknow = None
192 azimuth = None
194 azimuth = None
193 zenith = None
195 zenith = None
194 beam = Beam()
196 beam = Beam()
195 profileIndex = None
197 profileIndex = None
196 error = None
198 error = None
197 data = None
199 data = None
198 nmodes = None
200 nmodes = None
199
201
200 def __str__(self):
202 def __str__(self):
201
203
202 return '{} - {}'.format(self.type, self.getDatatime())
204 return '{} - {}'.format(self.type, self.getDatatime())
203
205
204 def getNoise(self):
206 def getNoise(self):
205
207
206 raise NotImplementedError
208 raise NotImplementedError
207
209
208 def getNChannels(self):
210 def getNChannels(self):
209
211
210 return len(self.channelList)
212 return len(self.channelList)
211
213
212 def getChannelIndexList(self):
214 def getChannelIndexList(self):
213
215
214 return list(range(self.nChannels))
216 return list(range(self.nChannels))
215
217
216 def getNHeights(self):
218 def getNHeights(self):
217
219
218 return len(self.heightList)
220 return len(self.heightList)
219
221
220 def getHeiRange(self, extrapoints=0):
222 def getHeiRange(self, extrapoints=0):
221
223
222 heis = self.heightList
224 heis = self.heightList
223 # deltah = self.heightList[1] - self.heightList[0]
225 # deltah = self.heightList[1] - self.heightList[0]
224 #
226 #
225 # heis.append(self.heightList[-1])
227 # heis.append(self.heightList[-1])
226
228
227 return heis
229 return heis
228
230
229 def getDeltaH(self):
231 def getDeltaH(self):
230
232
231 delta = self.heightList[1] - self.heightList[0]
233 delta = self.heightList[1] - self.heightList[0]
232
234
233 return delta
235 return delta
234
236
235 def getltctime(self):
237 def getltctime(self):
236
238
237 if self.useLocalTime:
239 if self.useLocalTime:
238 return self.utctime - self.timeZone * 60
240 return self.utctime - self.timeZone * 60
239
241
240 return self.utctime
242 return self.utctime
241
243
242 def getDatatime(self):
244 def getDatatime(self):
243
245
244 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
246 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
245 return datatimeValue
247 return datatimeValue
246
248
247 def getTimeRange(self):
249 def getTimeRange(self):
248
250
249 datatime = []
251 datatime = []
250
252
251 datatime.append(self.ltctime)
253 datatime.append(self.ltctime)
252 datatime.append(self.ltctime + self.timeInterval + 1)
254 datatime.append(self.ltctime + self.timeInterval + 1)
253
255
254 datatime = numpy.array(datatime)
256 datatime = numpy.array(datatime)
255
257
256 return datatime
258 return datatime
257
259
258 def getFmaxTimeResponse(self):
260 def getFmaxTimeResponse(self):
259
261
260 period = (10**-6) * self.getDeltaH() / (0.15)
262 period = (10**-6) * self.getDeltaH() / (0.15)
261
263
262 PRF = 1. / (period * self.nCohInt)
264 PRF = 1. / (period * self.nCohInt)
263
265
264 fmax = PRF
266 fmax = PRF
265
267
266 return fmax
268 return fmax
267
269
268 def getFmax(self):
270 def getFmax(self):
269 PRF = 1. / (self.ippSeconds * self.nCohInt)
271 PRF = 1. / (self.ippSeconds * self.nCohInt)
270
272
271 fmax = PRF
273 fmax = PRF
272 return fmax
274 return fmax
273
275
274 def getVmax(self):
276 def getVmax(self):
275
277
276 _lambda = self.C / self.frequency
278 _lambda = self.C / self.frequency
277
279
278 vmax = self.getFmax() * _lambda / 2
280 vmax = self.getFmax() * _lambda / 2
279
281
280 return vmax
282 return vmax
281
283
282 def get_ippSeconds(self):
284 def get_ippSeconds(self):
283 '''
285 '''
284 '''
286 '''
285 return self.radarControllerHeaderObj.ippSeconds
287 return self.radarControllerHeaderObj.ippSeconds
286
288
287 def set_ippSeconds(self, ippSeconds):
289 def set_ippSeconds(self, ippSeconds):
288 '''
290 '''
289 '''
291 '''
290
292
291 self.radarControllerHeaderObj.ippSeconds = ippSeconds
293 self.radarControllerHeaderObj.ippSeconds = ippSeconds
292
294
293 return
295 return
294
296
295 def get_dtype(self):
297 def get_dtype(self):
296 '''
298 '''
297 '''
299 '''
298 return getNumpyDtype(self.datatype)
300 return getNumpyDtype(self.datatype)
299
301
300 def set_dtype(self, numpyDtype):
302 def set_dtype(self, numpyDtype):
301 '''
303 '''
302 '''
304 '''
303
305
304 self.datatype = getDataTypeCode(numpyDtype)
306 self.datatype = getDataTypeCode(numpyDtype)
305
307
306 def get_code(self):
308 def get_code(self):
307 '''
309 '''
308 '''
310 '''
309 return self.radarControllerHeaderObj.code
311 return self.radarControllerHeaderObj.code
310
312
311 def set_code(self, code):
313 def set_code(self, code):
312 '''
314 '''
313 '''
315 '''
314 self.radarControllerHeaderObj.code = code
316 self.radarControllerHeaderObj.code = code
315
317
316 return
318 return
317
319
318 def get_ncode(self):
320 def get_ncode(self):
319 '''
321 '''
320 '''
322 '''
321 return self.radarControllerHeaderObj.nCode
323 return self.radarControllerHeaderObj.nCode
322
324
323 def set_ncode(self, nCode):
325 def set_ncode(self, nCode):
324 '''
326 '''
325 '''
327 '''
326 self.radarControllerHeaderObj.nCode = nCode
328 self.radarControllerHeaderObj.nCode = nCode
327
329
328 return
330 return
329
331
330 def get_nbaud(self):
332 def get_nbaud(self):
331 '''
333 '''
332 '''
334 '''
333 return self.radarControllerHeaderObj.nBaud
335 return self.radarControllerHeaderObj.nBaud
334
336
335 def set_nbaud(self, nBaud):
337 def set_nbaud(self, nBaud):
336 '''
338 '''
337 '''
339 '''
338 self.radarControllerHeaderObj.nBaud = nBaud
340 self.radarControllerHeaderObj.nBaud = nBaud
339
341
340 return
342 return
341
343
342 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
344 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
343 channelIndexList = property(
345 channelIndexList = property(
344 getChannelIndexList, "I'm the 'channelIndexList' property.")
346 getChannelIndexList, "I'm the 'channelIndexList' property.")
345 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
347 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
346 #noise = property(getNoise, "I'm the 'nHeights' property.")
348 #noise = property(getNoise, "I'm the 'nHeights' property.")
347 datatime = property(getDatatime, "I'm the 'datatime' property")
349 datatime = property(getDatatime, "I'm the 'datatime' property")
348 ltctime = property(getltctime, "I'm the 'ltctime' property")
350 ltctime = property(getltctime, "I'm the 'ltctime' property")
349 ippSeconds = property(get_ippSeconds, set_ippSeconds)
351 ippSeconds = property(get_ippSeconds, set_ippSeconds)
350 dtype = property(get_dtype, set_dtype)
352 dtype = property(get_dtype, set_dtype)
351 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
353 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
352 code = property(get_code, set_code)
354 code = property(get_code, set_code)
353 nCode = property(get_ncode, set_ncode)
355 nCode = property(get_ncode, set_ncode)
354 nBaud = property(get_nbaud, set_nbaud)
356 nBaud = property(get_nbaud, set_nbaud)
355
357
356
358
357 class Voltage(JROData):
359 class Voltage(JROData):
358
360
359 # data es un numpy array de 2 dmensiones (canales, alturas)
361 # data es un numpy array de 2 dmensiones (canales, alturas)
360 data = None
362 data = None
361
363
362 def __init__(self):
364 def __init__(self):
363 '''
365 '''
364 Constructor
366 Constructor
365 '''
367 '''
366
368
367 self.useLocalTime = True
369 self.useLocalTime = True
368 self.radarControllerHeaderObj = RadarControllerHeader()
370 self.radarControllerHeaderObj = RadarControllerHeader()
369 self.systemHeaderObj = SystemHeader()
371 self.systemHeaderObj = SystemHeader()
370 self.type = "Voltage"
372 self.type = "Voltage"
371 self.data = None
373 self.data = None
372 # self.dtype = None
374 # self.dtype = None
373 # self.nChannels = 0
375 # self.nChannels = 0
374 # self.nHeights = 0
376 # self.nHeights = 0
375 self.nProfiles = None
377 self.nProfiles = None
376 self.heightList = None
378 self.heightList = None
377 self.channelList = None
379 self.channelList = None
378 # self.channelIndexList = None
380 # self.channelIndexList = None
379 self.flagNoData = True
381 self.flagNoData = True
380 self.flagDiscontinuousBlock = False
382 self.flagDiscontinuousBlock = False
381 self.utctime = None
383 self.utctime = None
382 self.timeZone = None
384 self.timeZone = None
383 self.dstFlag = None
385 self.dstFlag = None
384 self.errorCount = None
386 self.errorCount = None
385 self.nCohInt = None
387 self.nCohInt = None
386 self.blocksize = None
388 self.blocksize = None
387 self.flagDecodeData = False # asumo q la data no esta decodificada
389 self.flagDecodeData = False # asumo q la data no esta decodificada
388 self.flagDeflipData = False # asumo q la data no esta sin flip
390 self.flagDeflipData = False # asumo q la data no esta sin flip
389 self.flagShiftFFT = False
391 self.flagShiftFFT = False
390 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
392 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
391 self.profileIndex = 0
393 self.profileIndex = 0
392
394
393 def getNoisebyHildebrand(self, channel=None):
395 def getNoisebyHildebrand(self, channel=None):
394 """
396 """
395 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
397 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
396
398
397 Return:
399 Return:
398 noiselevel
400 noiselevel
399 """
401 """
400
402
401 if channel != None:
403 if channel != None:
402 data = self.data[channel]
404 data = self.data[channel]
403 nChannels = 1
405 nChannels = 1
404 else:
406 else:
405 data = self.data
407 data = self.data
406 nChannels = self.nChannels
408 nChannels = self.nChannels
407
409
408 noise = numpy.zeros(nChannels)
410 noise = numpy.zeros(nChannels)
409 power = data * numpy.conjugate(data)
411 power = data * numpy.conjugate(data)
410
412
411 for thisChannel in range(nChannels):
413 for thisChannel in range(nChannels):
412 if nChannels == 1:
414 if nChannels == 1:
413 daux = power[:].real
415 daux = power[:].real
414 else:
416 else:
415 daux = power[thisChannel, :].real
417 daux = power[thisChannel, :].real
416 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
418 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
417
419
418 return noise
420 return noise
419
421
420 def getNoise(self, type=1, channel=None):
422 def getNoise(self, type=1, channel=None):
421
423
422 if type == 1:
424 if type == 1:
423 noise = self.getNoisebyHildebrand(channel)
425 noise = self.getNoisebyHildebrand(channel)
424
426
425 return noise
427 return noise
426
428
427 def getPower(self, channel=None):
429 def getPower(self, channel=None):
428
430
429 if channel != None:
431 if channel != None:
430 data = self.data[channel]
432 data = self.data[channel]
431 else:
433 else:
432 data = self.data
434 data = self.data
433
435
434 power = data * numpy.conjugate(data)
436 power = data * numpy.conjugate(data)
435 powerdB = 10 * numpy.log10(power.real)
437 powerdB = 10 * numpy.log10(power.real)
436 powerdB = numpy.squeeze(powerdB)
438 powerdB = numpy.squeeze(powerdB)
437
439
438 return powerdB
440 return powerdB
439
441
440 def getTimeInterval(self):
442 def getTimeInterval(self):
441
443
442 timeInterval = self.ippSeconds * self.nCohInt
444 timeInterval = self.ippSeconds * self.nCohInt
443
445
444 return timeInterval
446 return timeInterval
445
447
446 noise = property(getNoise, "I'm the 'nHeights' property.")
448 noise = property(getNoise, "I'm the 'nHeights' property.")
447 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
449 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
448
450
449
451
450 class Spectra(JROData):
452 class Spectra(JROData):
451
453
452 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
454 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
453 data_spc = None
455 data_spc = None
454 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
456 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
455 data_cspc = None
457 data_cspc = None
456 # data dc es un numpy array de 2 dmensiones (canales, alturas)
458 # data dc es un numpy array de 2 dmensiones (canales, alturas)
457 data_dc = None
459 data_dc = None
458 # data power
460 # data power
459 data_pwr = None
461 data_pwr = None
460 nFFTPoints = None
462 nFFTPoints = None
461 # nPairs = None
463 # nPairs = None
462 pairsList = None
464 pairsList = None
463 nIncohInt = None
465 nIncohInt = None
464 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
466 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
465 nCohInt = None # se requiere para determinar el valor de timeInterval
467 nCohInt = None # se requiere para determinar el valor de timeInterval
466 ippFactor = None
468 ippFactor = None
467 profileIndex = 0
469 profileIndex = 0
468 plotting = "spectra"
470 plotting = "spectra"
469
471
470 def __init__(self):
472 def __init__(self):
471 '''
473 '''
472 Constructor
474 Constructor
473 '''
475 '''
474
476
475 self.useLocalTime = True
477 self.useLocalTime = True
476 self.radarControllerHeaderObj = RadarControllerHeader()
478 self.radarControllerHeaderObj = RadarControllerHeader()
477 self.systemHeaderObj = SystemHeader()
479 self.systemHeaderObj = SystemHeader()
478 self.type = "Spectra"
480 self.type = "Spectra"
479 # self.data = None
481 # self.data = None
480 # self.dtype = None
482 # self.dtype = None
481 # self.nChannels = 0
483 # self.nChannels = 0
482 # self.nHeights = 0
484 # self.nHeights = 0
483 self.nProfiles = None
485 self.nProfiles = None
484 self.heightList = None
486 self.heightList = None
485 self.channelList = None
487 self.channelList = None
486 # self.channelIndexList = None
488 # self.channelIndexList = None
487 self.pairsList = None
489 self.pairsList = None
488 self.flagNoData = True
490 self.flagNoData = True
489 self.flagDiscontinuousBlock = False
491 self.flagDiscontinuousBlock = False
490 self.utctime = None
492 self.utctime = None
491 self.nCohInt = None
493 self.nCohInt = None
492 self.nIncohInt = None
494 self.nIncohInt = None
493 self.blocksize = None
495 self.blocksize = None
494 self.nFFTPoints = None
496 self.nFFTPoints = None
495 self.wavelength = None
497 self.wavelength = None
496 self.flagDecodeData = False # asumo q la data no esta decodificada
498 self.flagDecodeData = False # asumo q la data no esta decodificada
497 self.flagDeflipData = False # asumo q la data no esta sin flip
499 self.flagDeflipData = False # asumo q la data no esta sin flip
498 self.flagShiftFFT = False
500 self.flagShiftFFT = False
499 self.ippFactor = 1
501 self.ippFactor = 1
500 #self.noise = None
502 #self.noise = None
501 self.beacon_heiIndexList = []
503 self.beacon_heiIndexList = []
502 self.noise_estimation = None
504 self.noise_estimation = None
503
505
504 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
506 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
505 """
507 """
506 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
508 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
507
509
508 Return:
510 Return:
509 noiselevel
511 noiselevel
510 """
512 """
511
513
512 noise = numpy.zeros(self.nChannels)
514 noise = numpy.zeros(self.nChannels)
513
515
514 for channel in range(self.nChannels):
516 for channel in range(self.nChannels):
515 daux = self.data_spc[channel,
517 daux = self.data_spc[channel,
516 xmin_index:xmax_index, ymin_index:ymax_index]
518 xmin_index:xmax_index, ymin_index:ymax_index]
517 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
519 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
518
520
519 return noise
521 return noise
520
522
521 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
523 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
522
524
523 if self.noise_estimation is not None:
525 if self.noise_estimation is not None:
524 # this was estimated by getNoise Operation defined in jroproc_spectra.py
526 # this was estimated by getNoise Operation defined in jroproc_spectra.py
525 return self.noise_estimation
527 return self.noise_estimation
526 else:
528 else:
527 noise = self.getNoisebyHildebrand(
529 noise = self.getNoisebyHildebrand(
528 xmin_index, xmax_index, ymin_index, ymax_index)
530 xmin_index, xmax_index, ymin_index, ymax_index)
529 return noise
531 return noise
530
532
531 def getFreqRangeTimeResponse(self, extrapoints=0):
533 def getFreqRangeTimeResponse(self, extrapoints=0):
532
534
533 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
535 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
534 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
536 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
535
537
536 return freqrange
538 return freqrange
537
539
538 def getAcfRange(self, extrapoints=0):
540 def getAcfRange(self, extrapoints=0):
539
541
540 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
542 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
541 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
543 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
542
544
543 return freqrange
545 return freqrange
544
546
545 def getFreqRange(self, extrapoints=0):
547 def getFreqRange(self, extrapoints=0):
546
548
547 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
549 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
548 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
550 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
549
551
550 return freqrange
552 return freqrange
551
553
552 def getVelRange(self, extrapoints=0):
554 def getVelRange(self, extrapoints=0):
553
555
554 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
556 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
555 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
557 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
556
558
557 if self.nmodes:
559 if self.nmodes:
558 return velrange/self.nmodes
560 return velrange/self.nmodes
559 else:
561 else:
560 return velrange
562 return velrange
561
563
562 def getNPairs(self):
564 def getNPairs(self):
563
565
564 return len(self.pairsList)
566 return len(self.pairsList)
565
567
566 def getPairsIndexList(self):
568 def getPairsIndexList(self):
567
569
568 return list(range(self.nPairs))
570 return list(range(self.nPairs))
569
571
570 def getNormFactor(self):
572 def getNormFactor(self):
571
573
572 pwcode = 1
574 pwcode = 1
573
575
574 if self.flagDecodeData:
576 if self.flagDecodeData:
575 pwcode = numpy.sum(self.code[0]**2)
577 pwcode = numpy.sum(self.code[0]**2)
576 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
578 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
577 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
579 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
578
580
579 return normFactor
581 return normFactor
580
582
581 def getFlagCspc(self):
583 def getFlagCspc(self):
582
584
583 if self.data_cspc is None:
585 if self.data_cspc is None:
584 return True
586 return True
585
587
586 return False
588 return False
587
589
588 def getFlagDc(self):
590 def getFlagDc(self):
589
591
590 if self.data_dc is None:
592 if self.data_dc is None:
591 return True
593 return True
592
594
593 return False
595 return False
594
596
595 def getTimeInterval(self):
597 def getTimeInterval(self):
596
598
597 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
599 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
598 if self.nmodes:
600 if self.nmodes:
599 return self.nmodes*timeInterval
601 return self.nmodes*timeInterval
600 else:
602 else:
601 return timeInterval
603 return timeInterval
602
604
603 def getPower(self):
605 def getPower(self):
604
606
605 factor = self.normFactor
607 factor = self.normFactor
606 z = self.data_spc / factor
608 z = self.data_spc / factor
607 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
609 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
608 avg = numpy.average(z, axis=1)
610 avg = numpy.average(z, axis=1)
609
611
610 return 10 * numpy.log10(avg)
612 return 10 * numpy.log10(avg)
611
613
612 def getCoherence(self, pairsList=None, phase=False):
614 def getCoherence(self, pairsList=None, phase=False):
613
615
614 z = []
616 z = []
615 if pairsList is None:
617 if pairsList is None:
616 pairsIndexList = self.pairsIndexList
618 pairsIndexList = self.pairsIndexList
617 else:
619 else:
618 pairsIndexList = []
620 pairsIndexList = []
619 for pair in pairsList:
621 for pair in pairsList:
620 if pair not in self.pairsList:
622 if pair not in self.pairsList:
621 raise ValueError("Pair %s is not in dataOut.pairsList" % (
623 raise ValueError("Pair %s is not in dataOut.pairsList" % (
622 pair))
624 pair))
623 pairsIndexList.append(self.pairsList.index(pair))
625 pairsIndexList.append(self.pairsList.index(pair))
624 for i in range(len(pairsIndexList)):
626 for i in range(len(pairsIndexList)):
625 pair = self.pairsList[pairsIndexList[i]]
627 pair = self.pairsList[pairsIndexList[i]]
626 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
628 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
627 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
629 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
628 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
630 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
629 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
631 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
630 if phase:
632 if phase:
631 data = numpy.arctan2(avgcoherenceComplex.imag,
633 data = numpy.arctan2(avgcoherenceComplex.imag,
632 avgcoherenceComplex.real) * 180 / numpy.pi
634 avgcoherenceComplex.real) * 180 / numpy.pi
633 else:
635 else:
634 data = numpy.abs(avgcoherenceComplex)
636 data = numpy.abs(avgcoherenceComplex)
635
637
636 z.append(data)
638 z.append(data)
637
639
638 return numpy.array(z)
640 return numpy.array(z)
639
641
640 def setValue(self, value):
642 def setValue(self, value):
641
643
642 print("This property should not be initialized")
644 print("This property should not be initialized")
643
645
644 return
646 return
645
647
646 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
648 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
647 pairsIndexList = property(
649 pairsIndexList = property(
648 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
650 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
649 normFactor = property(getNormFactor, setValue,
651 normFactor = property(getNormFactor, setValue,
650 "I'm the 'getNormFactor' property.")
652 "I'm the 'getNormFactor' property.")
651 flag_cspc = property(getFlagCspc, setValue)
653 flag_cspc = property(getFlagCspc, setValue)
652 flag_dc = property(getFlagDc, setValue)
654 flag_dc = property(getFlagDc, setValue)
653 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
655 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
654 timeInterval = property(getTimeInterval, setValue,
656 timeInterval = property(getTimeInterval, setValue,
655 "I'm the 'timeInterval' property")
657 "I'm the 'timeInterval' property")
656
658
657
659
658 class SpectraHeis(Spectra):
660 class SpectraHeis(Spectra):
659
661
660 data_spc = None
662 data_spc = None
661 data_cspc = None
663 data_cspc = None
662 data_dc = None
664 data_dc = None
663 nFFTPoints = None
665 nFFTPoints = None
664 # nPairs = None
666 # nPairs = None
665 pairsList = None
667 pairsList = None
666 nCohInt = None
668 nCohInt = None
667 nIncohInt = None
669 nIncohInt = None
668
670
669 def __init__(self):
671 def __init__(self):
670
672
671 self.radarControllerHeaderObj = RadarControllerHeader()
673 self.radarControllerHeaderObj = RadarControllerHeader()
672
674
673 self.systemHeaderObj = SystemHeader()
675 self.systemHeaderObj = SystemHeader()
674
676
675 self.type = "SpectraHeis"
677 self.type = "SpectraHeis"
676
678
677 # self.dtype = None
679 # self.dtype = None
678
680
679 # self.nChannels = 0
681 # self.nChannels = 0
680
682
681 # self.nHeights = 0
683 # self.nHeights = 0
682
684
683 self.nProfiles = None
685 self.nProfiles = None
684
686
685 self.heightList = None
687 self.heightList = None
686
688
687 self.channelList = None
689 self.channelList = None
688
690
689 # self.channelIndexList = None
691 # self.channelIndexList = None
690
692
691 self.flagNoData = True
693 self.flagNoData = True
692
694
693 self.flagDiscontinuousBlock = False
695 self.flagDiscontinuousBlock = False
694
696
695 # self.nPairs = 0
697 # self.nPairs = 0
696
698
697 self.utctime = None
699 self.utctime = None
698
700
699 self.blocksize = None
701 self.blocksize = None
700
702
701 self.profileIndex = 0
703 self.profileIndex = 0
702
704
703 self.nCohInt = 1
705 self.nCohInt = 1
704
706
705 self.nIncohInt = 1
707 self.nIncohInt = 1
706
708
707 def getNormFactor(self):
709 def getNormFactor(self):
708 pwcode = 1
710 pwcode = 1
709 if self.flagDecodeData:
711 if self.flagDecodeData:
710 pwcode = numpy.sum(self.code[0]**2)
712 pwcode = numpy.sum(self.code[0]**2)
711
713
712 normFactor = self.nIncohInt * self.nCohInt * pwcode
714 normFactor = self.nIncohInt * self.nCohInt * pwcode
713
715
714 return normFactor
716 return normFactor
715
717
716 def getTimeInterval(self):
718 def getTimeInterval(self):
717
719
718 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
720 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
719
721
720 return timeInterval
722 return timeInterval
721
723
722 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
724 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
723 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
725 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
724
726
725
727
726 class Fits(JROData):
728 class Fits(JROData):
727
729
728 heightList = None
730 heightList = None
729 channelList = None
731 channelList = None
730 flagNoData = True
732 flagNoData = True
731 flagDiscontinuousBlock = False
733 flagDiscontinuousBlock = False
732 useLocalTime = False
734 useLocalTime = False
733 utctime = None
735 utctime = None
734 timeZone = None
736 timeZone = None
735 # ippSeconds = None
737 # ippSeconds = None
736 # timeInterval = None
738 # timeInterval = None
737 nCohInt = None
739 nCohInt = None
738 nIncohInt = None
740 nIncohInt = None
739 noise = None
741 noise = None
740 windowOfFilter = 1
742 windowOfFilter = 1
741 # Speed of ligth
743 # Speed of ligth
742 C = 3e8
744 C = 3e8
743 frequency = 49.92e6
745 frequency = 49.92e6
744 realtime = False
746 realtime = False
745
747
746 def __init__(self):
748 def __init__(self):
747
749
748 self.type = "Fits"
750 self.type = "Fits"
749
751
750 self.nProfiles = None
752 self.nProfiles = None
751
753
752 self.heightList = None
754 self.heightList = None
753
755
754 self.channelList = None
756 self.channelList = None
755
757
756 # self.channelIndexList = None
758 # self.channelIndexList = None
757
759
758 self.flagNoData = True
760 self.flagNoData = True
759
761
760 self.utctime = None
762 self.utctime = None
761
763
762 self.nCohInt = 1
764 self.nCohInt = 1
763
765
764 self.nIncohInt = 1
766 self.nIncohInt = 1
765
767
766 self.useLocalTime = True
768 self.useLocalTime = True
767
769
768 self.profileIndex = 0
770 self.profileIndex = 0
769
771
770 # self.utctime = None
772 # self.utctime = None
771 # self.timeZone = None
773 # self.timeZone = None
772 # self.ltctime = None
774 # self.ltctime = None
773 # self.timeInterval = None
775 # self.timeInterval = None
774 # self.header = None
776 # self.header = None
775 # self.data_header = None
777 # self.data_header = None
776 # self.data = None
778 # self.data = None
777 # self.datatime = None
779 # self.datatime = None
778 # self.flagNoData = False
780 # self.flagNoData = False
779 # self.expName = ''
781 # self.expName = ''
780 # self.nChannels = None
782 # self.nChannels = None
781 # self.nSamples = None
783 # self.nSamples = None
782 # self.dataBlocksPerFile = None
784 # self.dataBlocksPerFile = None
783 # self.comments = ''
785 # self.comments = ''
784 #
786 #
785
787
786 def getltctime(self):
788 def getltctime(self):
787
789
788 if self.useLocalTime:
790 if self.useLocalTime:
789 return self.utctime - self.timeZone * 60
791 return self.utctime - self.timeZone * 60
790
792
791 return self.utctime
793 return self.utctime
792
794
793 def getDatatime(self):
795 def getDatatime(self):
794
796
795 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
797 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
796 return datatime
798 return datatime
797
799
798 def getTimeRange(self):
800 def getTimeRange(self):
799
801
800 datatime = []
802 datatime = []
801
803
802 datatime.append(self.ltctime)
804 datatime.append(self.ltctime)
803 datatime.append(self.ltctime + self.timeInterval)
805 datatime.append(self.ltctime + self.timeInterval)
804
806
805 datatime = numpy.array(datatime)
807 datatime = numpy.array(datatime)
806
808
807 return datatime
809 return datatime
808
810
809 def getHeiRange(self):
811 def getHeiRange(self):
810
812
811 heis = self.heightList
813 heis = self.heightList
812
814
813 return heis
815 return heis
814
816
815 def getNHeights(self):
817 def getNHeights(self):
816
818
817 return len(self.heightList)
819 return len(self.heightList)
818
820
819 def getNChannels(self):
821 def getNChannels(self):
820
822
821 return len(self.channelList)
823 return len(self.channelList)
822
824
823 def getChannelIndexList(self):
825 def getChannelIndexList(self):
824
826
825 return list(range(self.nChannels))
827 return list(range(self.nChannels))
826
828
827 def getNoise(self, type=1):
829 def getNoise(self, type=1):
828
830
829 #noise = numpy.zeros(self.nChannels)
831 #noise = numpy.zeros(self.nChannels)
830
832
831 if type == 1:
833 if type == 1:
832 noise = self.getNoisebyHildebrand()
834 noise = self.getNoisebyHildebrand()
833
835
834 if type == 2:
836 if type == 2:
835 noise = self.getNoisebySort()
837 noise = self.getNoisebySort()
836
838
837 if type == 3:
839 if type == 3:
838 noise = self.getNoisebyWindow()
840 noise = self.getNoisebyWindow()
839
841
840 return noise
842 return noise
841
843
842 def getTimeInterval(self):
844 def getTimeInterval(self):
843
845
844 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
846 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
845
847
846 return timeInterval
848 return timeInterval
847
849
848 def get_ippSeconds(self):
850 def get_ippSeconds(self):
849 '''
851 '''
850 '''
852 '''
851 return self.ipp_sec
853 return self.ipp_sec
852
854
853
855
854 datatime = property(getDatatime, "I'm the 'datatime' property")
856 datatime = property(getDatatime, "I'm the 'datatime' property")
855 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
857 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
856 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
858 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
857 channelIndexList = property(
859 channelIndexList = property(
858 getChannelIndexList, "I'm the 'channelIndexList' property.")
860 getChannelIndexList, "I'm the 'channelIndexList' property.")
859 noise = property(getNoise, "I'm the 'nHeights' property.")
861 noise = property(getNoise, "I'm the 'nHeights' property.")
860
862
861 ltctime = property(getltctime, "I'm the 'ltctime' property")
863 ltctime = property(getltctime, "I'm the 'ltctime' property")
862 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
864 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
863 ippSeconds = property(get_ippSeconds, '')
865 ippSeconds = property(get_ippSeconds, '')
864
866
865 class Correlation(JROData):
867 class Correlation(JROData):
866
868
867 noise = None
869 noise = None
868 SNR = None
870 SNR = None
869 #--------------------------------------------------
871 #--------------------------------------------------
870 mode = None
872 mode = None
871 split = False
873 split = False
872 data_cf = None
874 data_cf = None
873 lags = None
875 lags = None
874 lagRange = None
876 lagRange = None
875 pairsList = None
877 pairsList = None
876 normFactor = None
878 normFactor = None
877 #--------------------------------------------------
879 #--------------------------------------------------
878 # calculateVelocity = None
880 # calculateVelocity = None
879 nLags = None
881 nLags = None
880 nPairs = None
882 nPairs = None
881 nAvg = None
883 nAvg = None
882
884
883 def __init__(self):
885 def __init__(self):
884 '''
886 '''
885 Constructor
887 Constructor
886 '''
888 '''
887 self.radarControllerHeaderObj = RadarControllerHeader()
889 self.radarControllerHeaderObj = RadarControllerHeader()
888
890
889 self.systemHeaderObj = SystemHeader()
891 self.systemHeaderObj = SystemHeader()
890
892
891 self.type = "Correlation"
893 self.type = "Correlation"
892
894
893 self.data = None
895 self.data = None
894
896
895 self.dtype = None
897 self.dtype = None
896
898
897 self.nProfiles = None
899 self.nProfiles = None
898
900
899 self.heightList = None
901 self.heightList = None
900
902
901 self.channelList = None
903 self.channelList = None
902
904
903 self.flagNoData = True
905 self.flagNoData = True
904
906
905 self.flagDiscontinuousBlock = False
907 self.flagDiscontinuousBlock = False
906
908
907 self.utctime = None
909 self.utctime = None
908
910
909 self.timeZone = None
911 self.timeZone = None
910
912
911 self.dstFlag = None
913 self.dstFlag = None
912
914
913 self.errorCount = None
915 self.errorCount = None
914
916
915 self.blocksize = None
917 self.blocksize = None
916
918
917 self.flagDecodeData = False # asumo q la data no esta decodificada
919 self.flagDecodeData = False # asumo q la data no esta decodificada
918
920
919 self.flagDeflipData = False # asumo q la data no esta sin flip
921 self.flagDeflipData = False # asumo q la data no esta sin flip
920
922
921 self.pairsList = None
923 self.pairsList = None
922
924
923 self.nPoints = None
925 self.nPoints = None
924
926
925 def getPairsList(self):
927 def getPairsList(self):
926
928
927 return self.pairsList
929 return self.pairsList
928
930
929 def getNoise(self, mode=2):
931 def getNoise(self, mode=2):
930
932
931 indR = numpy.where(self.lagR == 0)[0][0]
933 indR = numpy.where(self.lagR == 0)[0][0]
932 indT = numpy.where(self.lagT == 0)[0][0]
934 indT = numpy.where(self.lagT == 0)[0][0]
933
935
934 jspectra0 = self.data_corr[:, :, indR, :]
936 jspectra0 = self.data_corr[:, :, indR, :]
935 jspectra = copy.copy(jspectra0)
937 jspectra = copy.copy(jspectra0)
936
938
937 num_chan = jspectra.shape[0]
939 num_chan = jspectra.shape[0]
938 num_hei = jspectra.shape[2]
940 num_hei = jspectra.shape[2]
939
941
940 freq_dc = jspectra.shape[1] / 2
942 freq_dc = jspectra.shape[1] / 2
941 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
943 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
942
944
943 if ind_vel[0] < 0:
945 if ind_vel[0] < 0:
944 ind_vel[list(range(0, 1))] = ind_vel[list(
946 ind_vel[list(range(0, 1))] = ind_vel[list(
945 range(0, 1))] + self.num_prof
947 range(0, 1))] + self.num_prof
946
948
947 if mode == 1:
949 if mode == 1:
948 jspectra[:, freq_dc, :] = (
950 jspectra[:, freq_dc, :] = (
949 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
951 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
950
952
951 if mode == 2:
953 if mode == 2:
952
954
953 vel = numpy.array([-2, -1, 1, 2])
955 vel = numpy.array([-2, -1, 1, 2])
954 xx = numpy.zeros([4, 4])
956 xx = numpy.zeros([4, 4])
955
957
956 for fil in range(4):
958 for fil in range(4):
957 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
959 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
958
960
959 xx_inv = numpy.linalg.inv(xx)
961 xx_inv = numpy.linalg.inv(xx)
960 xx_aux = xx_inv[0, :]
962 xx_aux = xx_inv[0, :]
961
963
962 for ich in range(num_chan):
964 for ich in range(num_chan):
963 yy = jspectra[ich, ind_vel, :]
965 yy = jspectra[ich, ind_vel, :]
964 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
966 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
965
967
966 junkid = jspectra[ich, freq_dc, :] <= 0
968 junkid = jspectra[ich, freq_dc, :] <= 0
967 cjunkid = sum(junkid)
969 cjunkid = sum(junkid)
968
970
969 if cjunkid.any():
971 if cjunkid.any():
970 jspectra[ich, freq_dc, junkid.nonzero()] = (
972 jspectra[ich, freq_dc, junkid.nonzero()] = (
971 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
973 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
972
974
973 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
975 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
974
976
975 return noise
977 return noise
976
978
977 def getTimeInterval(self):
979 def getTimeInterval(self):
978
980
979 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
981 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
980
982
981 return timeInterval
983 return timeInterval
982
984
983 def splitFunctions(self):
985 def splitFunctions(self):
984
986
985 pairsList = self.pairsList
987 pairsList = self.pairsList
986 ccf_pairs = []
988 ccf_pairs = []
987 acf_pairs = []
989 acf_pairs = []
988 ccf_ind = []
990 ccf_ind = []
989 acf_ind = []
991 acf_ind = []
990 for l in range(len(pairsList)):
992 for l in range(len(pairsList)):
991 chan0 = pairsList[l][0]
993 chan0 = pairsList[l][0]
992 chan1 = pairsList[l][1]
994 chan1 = pairsList[l][1]
993
995
994 # Obteniendo pares de Autocorrelacion
996 # Obteniendo pares de Autocorrelacion
995 if chan0 == chan1:
997 if chan0 == chan1:
996 acf_pairs.append(chan0)
998 acf_pairs.append(chan0)
997 acf_ind.append(l)
999 acf_ind.append(l)
998 else:
1000 else:
999 ccf_pairs.append(pairsList[l])
1001 ccf_pairs.append(pairsList[l])
1000 ccf_ind.append(l)
1002 ccf_ind.append(l)
1001
1003
1002 data_acf = self.data_cf[acf_ind]
1004 data_acf = self.data_cf[acf_ind]
1003 data_ccf = self.data_cf[ccf_ind]
1005 data_ccf = self.data_cf[ccf_ind]
1004
1006
1005 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1007 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1006
1008
1007 def getNormFactor(self):
1009 def getNormFactor(self):
1008 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1010 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1009 acf_pairs = numpy.array(acf_pairs)
1011 acf_pairs = numpy.array(acf_pairs)
1010 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1012 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1011
1013
1012 for p in range(self.nPairs):
1014 for p in range(self.nPairs):
1013 pair = self.pairsList[p]
1015 pair = self.pairsList[p]
1014
1016
1015 ch0 = pair[0]
1017 ch0 = pair[0]
1016 ch1 = pair[1]
1018 ch1 = pair[1]
1017
1019
1018 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1020 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1019 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1021 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1020 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1022 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1021
1023
1022 return normFactor
1024 return normFactor
1023
1025
1024 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1026 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1025 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1027 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1026
1028
1027
1029
1028 class Parameters(Spectra):
1030 class Parameters(Spectra):
1029
1031
1030 experimentInfo = None # Information about the experiment
1032 experimentInfo = None # Information about the experiment
1031 # Information from previous data
1033 # Information from previous data
1032 inputUnit = None # Type of data to be processed
1034 inputUnit = None # Type of data to be processed
1033 operation = None # Type of operation to parametrize
1035 operation = None # Type of operation to parametrize
1034 # normFactor = None #Normalization Factor
1036 # normFactor = None #Normalization Factor
1035 groupList = None # List of Pairs, Groups, etc
1037 groupList = None # List of Pairs, Groups, etc
1036 # Parameters
1038 # Parameters
1037 data_param = None # Parameters obtained
1039 data_param = None # Parameters obtained
1038 data_pre = None # Data Pre Parametrization
1040 data_pre = None # Data Pre Parametrization
1039 data_SNR = None # Signal to Noise Ratio
1041 data_SNR = None # Signal to Noise Ratio
1040 # heightRange = None #Heights
1042 # heightRange = None #Heights
1041 abscissaList = None # Abscissa, can be velocities, lags or time
1043 abscissaList = None # Abscissa, can be velocities, lags or time
1042 # noise = None #Noise Potency
1044 # noise = None #Noise Potency
1043 utctimeInit = None # Initial UTC time
1045 utctimeInit = None # Initial UTC time
1044 paramInterval = None # Time interval to calculate Parameters in seconds
1046 paramInterval = None # Time interval to calculate Parameters in seconds
1045 useLocalTime = True
1047 useLocalTime = True
1046 # Fitting
1048 # Fitting
1047 data_error = None # Error of the estimation
1049 data_error = None # Error of the estimation
1048 constants = None
1050 constants = None
1049 library = None
1051 library = None
1050 # Output signal
1052 # Output signal
1051 outputInterval = None # Time interval to calculate output signal in seconds
1053 outputInterval = None # Time interval to calculate output signal in seconds
1052 data_output = None # Out signal
1054 data_output = None # Out signal
1053 nAvg = None
1055 nAvg = None
1054 noise_estimation = None
1056 noise_estimation = None
1055 GauSPC = None # Fit gaussian SPC
1057 GauSPC = None # Fit gaussian SPC
1056
1058
1057 def __init__(self):
1059 def __init__(self):
1058 '''
1060 '''
1059 Constructor
1061 Constructor
1060 '''
1062 '''
1061 self.radarControllerHeaderObj = RadarControllerHeader()
1063 self.radarControllerHeaderObj = RadarControllerHeader()
1062
1064
1063 self.systemHeaderObj = SystemHeader()
1065 self.systemHeaderObj = SystemHeader()
1064
1066
1065 self.type = "Parameters"
1067 self.type = "Parameters"
1066
1068
1067 def getTimeRange1(self, interval):
1069 def getTimeRange1(self, interval):
1068
1070
1069 datatime = []
1071 datatime = []
1070
1072
1071 if self.useLocalTime:
1073 if self.useLocalTime:
1072 time1 = self.utctimeInit - self.timeZone * 60
1074 time1 = self.utctimeInit - self.timeZone * 60
1073 else:
1075 else:
1074 time1 = self.utctimeInit
1076 time1 = self.utctimeInit
1075
1077
1076 datatime.append(time1)
1078 datatime.append(time1)
1077 datatime.append(time1 + interval)
1079 datatime.append(time1 + interval)
1078 datatime = numpy.array(datatime)
1080 datatime = numpy.array(datatime)
1079
1081
1080 return datatime
1082 return datatime
1081
1083
1082 def getTimeInterval(self):
1084 def getTimeInterval(self):
1083
1085
1084 if hasattr(self, 'timeInterval1'):
1086 if hasattr(self, 'timeInterval1'):
1085 return self.timeInterval1
1087 return self.timeInterval1
1086 else:
1088 else:
1087 return self.paramInterval
1089 return self.paramInterval
1088
1090
1089 def setValue(self, value):
1091 def setValue(self, value):
1090
1092
1091 print("This property should not be initialized")
1093 print("This property should not be initialized")
1092
1094
1093 return
1095 return
1094
1096
1095 def getNoise(self):
1097 def getNoise(self):
1096
1098
1097 return self.spc_noise
1099 return self.spc_noise
1098
1100
1099 timeInterval = property(getTimeInterval)
1101 timeInterval = property(getTimeInterval)
1100 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1102 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1101
1103
1102
1104
1103 class PlotterData(object):
1105 class PlotterData(object):
1104 '''
1106 '''
1105 Object to hold data to be plotted
1107 Object to hold data to be plotted
1106 '''
1108 '''
1107
1109
1108 MAXNUMX = 100
1110 MAXNUMX = 100
1109 MAXNUMY = 100
1111 MAXNUMY = 100
1110
1112
1111 def __init__(self, code, throttle_value, exp_code, buffering=True, snr=False):
1113 def __init__(self, code, throttle_value, exp_code, buffering=True, snr=False):
1112
1114
1113 self.key = code
1115 self.key = code
1114 self.throttle = throttle_value
1116 self.throttle = throttle_value
1115 self.exp_code = exp_code
1117 self.exp_code = exp_code
1116 self.buffering = buffering
1118 self.buffering = buffering
1117 self.ready = False
1119 self.ready = False
1118 self.flagNoData = False
1120 self.flagNoData = False
1119 self.localtime = False
1121 self.localtime = False
1120 self.data = {}
1122 self.data = {}
1121 self.meta = {}
1123 self.meta = {}
1122 self.__times = []
1124 self.__times = []
1123 self.__heights = []
1125 self.__heights = []
1124
1126
1125 if 'snr' in code:
1127 if 'snr' in code:
1126 self.plottypes = ['snr']
1128 self.plottypes = ['snr']
1127 elif code == 'spc':
1129 elif code == 'spc':
1128 self.plottypes = ['spc', 'noise', 'rti']
1130 self.plottypes = ['spc', 'noise', 'rti']
1129 elif code == 'rti':
1131 elif code == 'rti':
1130 self.plottypes = ['noise', 'rti']
1132 self.plottypes = ['noise', 'rti']
1131 else:
1133 else:
1132 self.plottypes = [code]
1134 self.plottypes = [code]
1133
1135
1134 if 'snr' not in self.plottypes and snr:
1136 if 'snr' not in self.plottypes and snr:
1135 self.plottypes.append('snr')
1137 self.plottypes.append('snr')
1136
1138
1137 for plot in self.plottypes:
1139 for plot in self.plottypes:
1138 self.data[plot] = {}
1140 self.data[plot] = {}
1139
1141
1140 def __str__(self):
1142 def __str__(self):
1141 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1143 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1142 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
1144 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
1143
1145
1144 def __len__(self):
1146 def __len__(self):
1145 return len(self.__times)
1147 return len(self.__times)
1146
1148
1147 def __getitem__(self, key):
1149 def __getitem__(self, key):
1148
1150
1149 if key not in self.data:
1151 if key not in self.data:
1150 raise KeyError(log.error('Missing key: {}'.format(key)))
1152 raise KeyError(log.error('Missing key: {}'.format(key)))
1151 if 'spc' in key or not self.buffering:
1153 if 'spc' in key or not self.buffering:
1152 ret = self.data[key]
1154 ret = self.data[key]
1153 elif 'scope' in key:
1155 elif 'scope' in key:
1154 ret = numpy.array(self.data[key][float(self.tm)])
1156 ret = numpy.array(self.data[key][float(self.tm)])
1155 else:
1157 else:
1156 ret = numpy.array([self.data[key][x] for x in self.times])
1158 ret = numpy.array([self.data[key][x] for x in self.times])
1157 if ret.ndim > 1:
1159 if ret.ndim > 1:
1158 ret = numpy.swapaxes(ret, 0, 1)
1160 ret = numpy.swapaxes(ret, 0, 1)
1159 return ret
1161 return ret
1160
1162
1161 def __contains__(self, key):
1163 def __contains__(self, key):
1162 return key in self.data
1164 return key in self.data
1163
1165
1164 def setup(self):
1166 def setup(self):
1165 '''
1167 '''
1166 Configure object
1168 Configure object
1167 '''
1169 '''
1168
1170
1169 self.type = ''
1171 self.type = ''
1170 self.ready = False
1172 self.ready = False
1171 self.data = {}
1173 self.data = {}
1172 self.__times = []
1174 self.__times = []
1173 self.__heights = []
1175 self.__heights = []
1174 self.__all_heights = set()
1176 self.__all_heights = set()
1175 for plot in self.plottypes:
1177 for plot in self.plottypes:
1176 if 'snr' in plot:
1178 if 'snr' in plot:
1177 plot = 'snr'
1179 plot = 'snr'
1178 elif 'spc_moments' == plot:
1180 elif 'spc_moments' == plot:
1179 plot = 'moments'
1181 plot = 'moments'
1180 self.data[plot] = {}
1182 self.data[plot] = {}
1181
1183
1182 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
1184 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
1183 self.data['noise'] = {}
1185 self.data['noise'] = {}
1184 self.data['rti'] = {}
1186 self.data['rti'] = {}
1185 if 'noise' not in self.plottypes:
1187 if 'noise' not in self.plottypes:
1186 self.plottypes.append('noise')
1188 self.plottypes.append('noise')
1187 if 'rti' not in self.plottypes:
1189 if 'rti' not in self.plottypes:
1188 self.plottypes.append('rti')
1190 self.plottypes.append('rti')
1189
1191
1190 def shape(self, key):
1192 def shape(self, key):
1191 '''
1193 '''
1192 Get the shape of the one-element data for the given key
1194 Get the shape of the one-element data for the given key
1193 '''
1195 '''
1194
1196
1195 if len(self.data[key]):
1197 if len(self.data[key]):
1196 if 'spc' in key or not self.buffering:
1198 if 'spc' in key or not self.buffering:
1197 return self.data[key].shape
1199 return self.data[key].shape
1198 return self.data[key][self.__times[0]].shape
1200 return self.data[key][self.__times[0]].shape
1199 return (0,)
1201 return (0,)
1200
1202
1201 def update(self, dataOut, tm):
1203 def update(self, dataOut, tm):
1202 '''
1204 '''
1203 Update data object with new dataOut
1205 Update data object with new dataOut
1204 '''
1206 '''
1205
1207
1206 if tm in self.__times:
1208 if tm in self.__times:
1207 return
1209 return
1208 self.profileIndex = dataOut.profileIndex
1210 self.profileIndex = dataOut.profileIndex
1209 self.tm = tm
1211 self.tm = tm
1210 self.type = dataOut.type
1212 self.type = dataOut.type
1211 self.parameters = getattr(dataOut, 'parameters', [])
1213 self.parameters = getattr(dataOut, 'parameters', [])
1212
1214
1213 if hasattr(dataOut, 'meta'):
1215 if hasattr(dataOut, 'meta'):
1214 self.meta.update(dataOut.meta)
1216 self.meta.update(dataOut.meta)
1215
1217
1216 if hasattr(dataOut, 'pairsList'):
1218 if hasattr(dataOut, 'pairsList'):
1217 self.pairs = dataOut.pairsList
1219 self.pairs = dataOut.pairsList
1218
1220
1219 self.interval = dataOut.getTimeInterval()
1221 self.interval = dataOut.getTimeInterval()
1220 self.localtime = dataOut.useLocalTime
1222 self.localtime = dataOut.useLocalTime
1221 if True in ['spc' in ptype for ptype in self.plottypes]:
1223 if True in ['spc' in ptype for ptype in self.plottypes]:
1222 self.xrange = (dataOut.getFreqRange(1)/1000.,
1224 self.xrange = (dataOut.getFreqRange(1)/1000.,
1223 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1225 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1224 self.factor = dataOut.normFactor
1226 self.factor = dataOut.normFactor
1225 self.__heights.append(dataOut.heightList)
1227 self.__heights.append(dataOut.heightList)
1226 self.__all_heights.update(dataOut.heightList)
1228 self.__all_heights.update(dataOut.heightList)
1227 self.__times.append(tm)
1229 self.__times.append(tm)
1228
1230
1229 for plot in self.plottypes:
1231 for plot in self.plottypes:
1230 if plot in ('spc', 'spc_moments', 'spc_cut'):
1232 if plot in ('spc', 'spc_moments', 'spc_cut'):
1231 z = dataOut.data_spc/dataOut.normFactor
1233 z = dataOut.data_spc/dataOut.normFactor
1232 buffer = 10*numpy.log10(z)
1234 buffer = 10*numpy.log10(z)
1233 if plot == 'cspc':
1235 if plot == 'cspc':
1234 z = dataOut.data_spc/dataOut.normFactor
1236 z = dataOut.data_spc/dataOut.normFactor
1235 buffer = (dataOut.data_spc, dataOut.data_cspc)
1237 buffer = (dataOut.data_spc, dataOut.data_cspc)
1236 if plot == 'noise':
1238 if plot == 'noise':
1237 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1239 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1238 if plot in ('rti', 'spcprofile'):
1240 if plot in ('rti', 'spcprofile'):
1239 buffer = dataOut.getPower()
1241 buffer = dataOut.getPower()
1240 if plot == 'snr_db':
1242 if plot == 'snr_db':
1241 buffer = dataOut.data_SNR
1243 buffer = dataOut.data_SNR
1242 if plot == 'snr':
1244 if plot == 'snr':
1243 buffer = 10*numpy.log10(dataOut.data_SNR)
1245 buffer = 10*numpy.log10(dataOut.data_SNR)
1244 if plot == 'dop':
1246 if plot == 'dop':
1245 buffer = dataOut.data_DOP
1247 buffer = dataOut.data_DOP
1246 if plot == 'pow':
1248 if plot == 'pow':
1247 buffer = 10*numpy.log10(dataOut.data_POW)
1249 buffer = 10*numpy.log10(dataOut.data_POW)
1248 if plot == 'width':
1250 if plot == 'width':
1249 buffer = dataOut.data_WIDTH
1251 buffer = dataOut.data_WIDTH
1250 if plot == 'coh':
1252 if plot == 'coh':
1251 buffer = dataOut.getCoherence()
1253 buffer = dataOut.getCoherence()
1252 if plot == 'phase':
1254 if plot == 'phase':
1253 buffer = dataOut.getCoherence(phase=True)
1255 buffer = dataOut.getCoherence(phase=True)
1254 if plot == 'output':
1256 if plot == 'output':
1255 buffer = dataOut.data_output
1257 buffer = dataOut.data_output
1256 if plot == 'param':
1258 if plot == 'param':
1257 buffer = dataOut.data_param
1259 buffer = dataOut.data_param
1258 if plot == 'scope':
1260 if plot == 'scope':
1259 buffer = dataOut.data
1261 buffer = dataOut.data
1260 self.flagDataAsBlock = dataOut.flagDataAsBlock
1262 self.flagDataAsBlock = dataOut.flagDataAsBlock
1261 self.nProfiles = dataOut.nProfiles
1263 self.nProfiles = dataOut.nProfiles
1262
1264
1263 if plot == 'spc':
1265 if plot == 'spc':
1264 self.data['spc'] = buffer
1266 self.data['spc'] = buffer
1265 elif plot == 'cspc':
1267 elif plot == 'cspc':
1266 self.data['spc'] = buffer[0]
1268 self.data['spc'] = buffer[0]
1267 self.data['cspc'] = buffer[1]
1269 self.data['cspc'] = buffer[1]
1268 elif plot == 'spc_moments':
1270 elif plot == 'spc_moments':
1269 self.data['spc'] = buffer
1271 self.data['spc'] = buffer
1270 self.data['moments'][tm] = dataOut.moments
1272 self.data['moments'][tm] = dataOut.moments
1271 else:
1273 else:
1272 if self.buffering:
1274 if self.buffering:
1273 self.data[plot][tm] = buffer
1275 self.data[plot][tm] = buffer
1274 else:
1276 else:
1275 self.data[plot] = buffer
1277 self.data[plot] = buffer
1276
1278
1277 if dataOut.channelList is None:
1279 if dataOut.channelList is None:
1278 self.channels = range(buffer.shape[0])
1280 self.channels = range(buffer.shape[0])
1279 else:
1281 else:
1280 self.channels = dataOut.channelList
1282 self.channels = dataOut.channelList
1281
1283
1282 if buffer is None:
1284 if buffer is None:
1283 self.flagNoData = True
1285 self.flagNoData = True
1284 raise schainpy.admin.SchainWarning('Attribute data_{} is empty'.format(self.key))
1286 raise schainpy.admin.SchainWarning('Attribute data_{} is empty'.format(self.key))
1285
1287
1286 def normalize_heights(self):
1288 def normalize_heights(self):
1287 '''
1289 '''
1288 Ensure same-dimension of the data for different heighList
1290 Ensure same-dimension of the data for different heighList
1289 '''
1291 '''
1290
1292
1291 H = numpy.array(list(self.__all_heights))
1293 H = numpy.array(list(self.__all_heights))
1292 H.sort()
1294 H.sort()
1293 for key in self.data:
1295 for key in self.data:
1294 shape = self.shape(key)[:-1] + H.shape
1296 shape = self.shape(key)[:-1] + H.shape
1295 for tm, obj in list(self.data[key].items()):
1297 for tm, obj in list(self.data[key].items()):
1296 h = self.__heights[self.__times.index(tm)]
1298 h = self.__heights[self.__times.index(tm)]
1297 if H.size == h.size:
1299 if H.size == h.size:
1298 continue
1300 continue
1299 index = numpy.where(numpy.in1d(H, h))[0]
1301 index = numpy.where(numpy.in1d(H, h))[0]
1300 dummy = numpy.zeros(shape) + numpy.nan
1302 dummy = numpy.zeros(shape) + numpy.nan
1301 if len(shape) == 2:
1303 if len(shape) == 2:
1302 dummy[:, index] = obj
1304 dummy[:, index] = obj
1303 else:
1305 else:
1304 dummy[index] = obj
1306 dummy[index] = obj
1305 self.data[key][tm] = dummy
1307 self.data[key][tm] = dummy
1306
1308
1307 self.__heights = [H for tm in self.__times]
1309 self.__heights = [H for tm in self.__times]
1308
1310
1309 def jsonify(self, plot_name, plot_type, decimate=False):
1311 def jsonify(self, plot_name, plot_type, decimate=False):
1310 '''
1312 '''
1311 Convert data to json
1313 Convert data to json
1312 '''
1314 '''
1313
1315
1314 tm = self.times[-1]
1316 tm = self.times[-1]
1315 dy = int(self.heights.size/self.MAXNUMY) + 1
1317 dy = int(self.heights.size/self.MAXNUMY) + 1
1316 if self.key in ('spc', 'cspc') or not self.buffering:
1318 if self.key in ('spc', 'cspc') or not self.buffering:
1317 dx = int(self.data[self.key].shape[1]/self.MAXNUMX) + 1
1319 dx = int(self.data[self.key].shape[1]/self.MAXNUMX) + 1
1318 data = self.roundFloats(
1320 data = self.roundFloats(
1319 self.data[self.key][::, ::dx, ::dy].tolist())
1321 self.data[self.key][::, ::dx, ::dy].tolist())
1320 else:
1322 else:
1321 data = self.roundFloats(self.data[self.key][tm].tolist())
1323 data = self.roundFloats(self.data[self.key][tm].tolist())
1322 if self.key is 'noise':
1324 if self.key is 'noise':
1323 data = [[x] for x in data]
1325 data = [[x] for x in data]
1324
1326
1325 meta = {}
1327 meta = {}
1326 ret = {
1328 ret = {
1327 'plot': plot_name,
1329 'plot': plot_name,
1328 'code': self.exp_code,
1330 'code': self.exp_code,
1329 'time': float(tm),
1331 'time': float(tm),
1330 'data': data,
1332 'data': data,
1331 }
1333 }
1332 meta['type'] = plot_type
1334 meta['type'] = plot_type
1333 meta['interval'] = float(self.interval)
1335 meta['interval'] = float(self.interval)
1334 meta['localtime'] = self.localtime
1336 meta['localtime'] = self.localtime
1335 meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1337 meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1336 if 'spc' in self.data or 'cspc' in self.data:
1338 if 'spc' in self.data or 'cspc' in self.data:
1337 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1339 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1338 else:
1340 else:
1339 meta['xrange'] = []
1341 meta['xrange'] = []
1340
1342
1341 meta.update(self.meta)
1343 meta.update(self.meta)
1342 ret['metadata'] = meta
1344 ret['metadata'] = meta
1343 return json.dumps(ret)
1345 return json.dumps(ret)
1344
1346
1345 @property
1347 @property
1346 def times(self):
1348 def times(self):
1347 '''
1349 '''
1348 Return the list of times of the current data
1350 Return the list of times of the current data
1349 '''
1351 '''
1350
1352
1351 ret = numpy.array(self.__times)
1353 ret = numpy.array(self.__times)
1352 ret.sort()
1354 ret.sort()
1353 return ret
1355 return ret
1354
1356
1355 @property
1357 @property
1356 def min_time(self):
1358 def min_time(self):
1357 '''
1359 '''
1358 Return the minimun time value
1360 Return the minimun time value
1359 '''
1361 '''
1360
1362
1361 return self.times[0]
1363 return self.times[0]
1362
1364
1363 @property
1365 @property
1364 def max_time(self):
1366 def max_time(self):
1365 '''
1367 '''
1366 Return the maximun time value
1368 Return the maximun time value
1367 '''
1369 '''
1368
1370
1369 return self.times[-1]
1371 return self.times[-1]
1370
1372
1371 @property
1373 @property
1372 def heights(self):
1374 def heights(self):
1373 '''
1375 '''
1374 Return the list of heights of the current data
1376 Return the list of heights of the current data
1375 '''
1377 '''
1376
1378
1377 return numpy.array(self.__heights[-1])
1379 return numpy.array(self.__heights[-1])
1378
1380
1379 @staticmethod
1381 @staticmethod
1380 def roundFloats(obj):
1382 def roundFloats(obj):
1381 if isinstance(obj, list):
1383 if isinstance(obj, list):
1382 return list(map(PlotterData.roundFloats, obj))
1384 return list(map(PlotterData.roundFloats, obj))
1383 elif isinstance(obj, float):
1385 elif isinstance(obj, float):
1384 return round(obj, 2)
1386 return round(obj, 2)
Show file before | Show file after | Hide comments
@@ -1,63 +1,68
1 '''
1 '''
2 Created on Jul 16, 2014
2 Created on Jul 16, 2014
3
3
4 @author: Miguel Urco
4 @author: Miguel Urco
5 @author: Juan C. Espinoza
5 '''
6 '''
6
7
7 import os
8 import os
8 from setuptools import setup, Extension
9 from setuptools import setup, Extension
9 from setuptools.command.build_ext import build_ext as _build_ext
10 from setuptools.command.build_ext import build_ext as _build_ext
10 from schainpy import __version__
11 from schainpy import __version__
11
12
12 class build_ext(_build_ext):
13 class build_ext(_build_ext):
13 def finalize_options(self):
14 def finalize_options(self):
14 _build_ext.finalize_options(self)
15 _build_ext.finalize_options(self)
15 # Prevent numpy from thinking it is still in its setup process:
16 # Prevent numpy from thinking it is still in its setup process:
16 __builtins__.__NUMPY_SETUP__ = False
17 __builtins__.__NUMPY_SETUP__ = False
17 import numpy
18 import numpy
18 self.include_dirs.append(numpy.get_include())
19 self.include_dirs.append(numpy.get_include())
19
20
20 setup(name = "schainpy",
21 setup(
21 version = __version__,
22 name = "schainpy",
22 description = "Python tools to read, write and process Jicamarca data",
23 version = __version__,
23 author = "Miguel Urco",
24 description = "Python tools to read, write and process Jicamarca data",
24 author_email = "miguel.urco@jro.igp.gob.pe",
25 author = "Miguel Urco, Juan C. Espinoza",
25 url = "http://jro.igp.gob.pe",
26 author_email = "juan.espinoza@jro.igp.gob.pe",
26 packages = {'schainpy',
27 url = "http://jro-dev.igp.gob.pe/rhodecode/schain",
27 'schainpy.model',
28 packages = {
28 'schainpy.model.data',
29 'schainpy',
29 'schainpy.model.graphics',
30 'schainpy.model',
30 'schainpy.model.io',
31 'schainpy.model.data',
31 'schainpy.model.proc',
32 'schainpy.model.graphics',
32 'schainpy.model.serializer',
33 'schainpy.model.io',
33 'schainpy.model.utils',
34 'schainpy.model.proc',
34 'schainpy.utils',
35 'schainpy.model.utils',
35 'schainpy.gui',
36 'schainpy.utils',
36 'schainpy.gui.figures',
37 'schainpy.gui',
37 'schainpy.gui.viewcontroller',
38 'schainpy.gui.figures',
38 'schainpy.gui.viewer',
39 'schainpy.gui.viewcontroller',
39 'schainpy.gui.viewer.windows',
40 'schainpy.gui.viewer',
40 'schainpy.cli'},
41 'schainpy.gui.viewer.windows',
41 ext_package = 'schainpy',
42 'schainpy.cli',
42 package_data = {'': ['schain.conf.template'],
43 },
43 'schainpy.gui.figures': ['*.png', '*.jpg'],
44 package_data = {'': ['schain.conf.template'],
44 'schainpy.files': ['*.oga']
45 'schainpy.gui.figures': ['*.png', '*.jpg'],
45 },
46 'schainpy.files': ['*.oga']
46 include_package_data = False,
47 },
47 scripts = ['schainpy/gui/schainGUI'],
48 include_package_data = False,
48 entry_points = {
49 scripts = ['schainpy/gui/schainGUI'],
49 'console_scripts': [
50 entry_points = {
50 'schain = schainpy.cli.cli:main',
51 'console_scripts': [
51 ],
52 'schain = schainpy.cli.cli:main',
52 },
53 ],
53 cmdclass = {'build_ext': build_ext},
54 },
54 setup_requires = ["numpy >= 1.11.2"],
55 cmdclass = {'build_ext': build_ext},
55 install_requires = [
56 ext_modules=[
56 "scipy",
57 Extension("schainpy.model.data._noise", ["schainc/_noise.c"]),
57 "h5py",
58 ],
58 "matplotlib",
59 setup_requires = ["numpy"],
59 "pyzmq",
60 install_requires = [
60 "fuzzywuzzy",
61 "scipy",
61 "click",
62 "h5py",
62 ],
63 "matplotlib",
64 "pyzmq",
65 "fuzzywuzzy",
66 "click",
67 ],
63 )
68 )
General Comments 0
You need to be logged in to leave comments. Login now