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