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