##// END OF EJS Templates
AMISR correcciones, procesamientos
joabAM -
r1417:1ff56cfe5bef
parent child
Show More
@@ -1,657 +1,659
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 # All rights reserved.
2 # All rights reserved.
3 #
3 #
4 # Distributed under the terms of the BSD 3-clause license.
4 # Distributed under the terms of the BSD 3-clause license.
5 """API to create signal chain projects
5 """API to create signal chain projects
6
6
7 The API is provide through class: Project
7 The API is provide through class: Project
8 """
8 """
9
9
10 import re
10 import re
11 import sys
11 import sys
12 import ast
12 import ast
13 import datetime
13 import datetime
14 import traceback
14 import traceback
15 import time
15 import time
16 import multiprocessing
16 import multiprocessing
17 from multiprocessing import Process, Queue
17 from multiprocessing import Process, Queue
18 from threading import Thread
18 from threading import Thread
19 from xml.etree.ElementTree import ElementTree, Element, SubElement
19 from xml.etree.ElementTree import ElementTree, Element, SubElement
20
20
21 from schainpy.admin import Alarm, SchainWarning
21 from schainpy.admin import Alarm, SchainWarning
22 from schainpy.model import *
22 from schainpy.model import *
23 from schainpy.utils import log
23 from schainpy.utils import log
24
24
25 if 'darwin' in sys.platform and sys.version_info[0] == 3 and sys.version_info[1] > 7:
25 if 'darwin' in sys.platform and sys.version_info[0] == 3 and sys.version_info[1] > 7:
26 multiprocessing.set_start_method('fork')
26 multiprocessing.set_start_method('fork')
27
27
28 class ConfBase():
28 class ConfBase():
29
29
30 def __init__(self):
30 def __init__(self):
31
31
32 self.id = '0'
32 self.id = '0'
33 self.name = None
33 self.name = None
34 self.priority = None
34 self.priority = None
35 self.parameters = {}
35 self.parameters = {}
36 self.object = None
36 self.object = None
37 self.operations = []
37 self.operations = []
38
38
39 def getId(self):
39 def getId(self):
40
40
41 return self.id
41 return self.id
42
42
43 def getNewId(self):
43 def getNewId(self):
44
44
45 return int(self.id) * 10 + len(self.operations) + 1
45 return int(self.id) * 10 + len(self.operations) + 1
46
46
47 def updateId(self, new_id):
47 def updateId(self, new_id):
48
48
49 self.id = str(new_id)
49 self.id = str(new_id)
50
50
51 n = 1
51 n = 1
52 for conf in self.operations:
52 for conf in self.operations:
53 conf_id = str(int(new_id) * 10 + n)
53 conf_id = str(int(new_id) * 10 + n)
54 conf.updateId(conf_id)
54 conf.updateId(conf_id)
55 n += 1
55 n += 1
56
56
57 def getKwargs(self):
57 def getKwargs(self):
58
58
59 params = {}
59 params = {}
60
60
61 for key, value in self.parameters.items():
61 for key, value in self.parameters.items():
62 if value not in (None, '', ' '):
62 if value not in (None, '', ' '):
63 params[key] = value
63 params[key] = value
64
64
65 return params
65 return params
66
66
67 def update(self, **kwargs):
67 def update(self, **kwargs):
68
68
69 for key, value in kwargs.items():
69 for key, value in kwargs.items():
70 self.addParameter(name=key, value=value)
70 self.addParameter(name=key, value=value)
71
71
72 def addParameter(self, name, value, format=None):
72 def addParameter(self, name, value, format=None):
73 '''
73 '''
74 '''
74 '''
75
75
76 if isinstance(value, str) and re.search(r'(\d+/\d+/\d+)', value):
76 if isinstance(value, str) and re.search(r'(\d+/\d+/\d+)', value):
77 self.parameters[name] = datetime.date(*[int(x) for x in value.split('/')])
77 self.parameters[name] = datetime.date(*[int(x) for x in value.split('/')])
78 elif isinstance(value, str) and re.search(r'(\d+:\d+:\d+)', value):
78 elif isinstance(value, str) and re.search(r'(\d+:\d+:\d+)', value):
79 self.parameters[name] = datetime.time(*[int(x) for x in value.split(':')])
79 self.parameters[name] = datetime.time(*[int(x) for x in value.split(':')])
80 else:
80 else:
81 try:
81 try:
82 self.parameters[name] = ast.literal_eval(value)
82 self.parameters[name] = ast.literal_eval(value)
83 except:
83 except:
84 if isinstance(value, str) and ',' in value:
84 if isinstance(value, str) and ',' in value:
85 self.parameters[name] = value.split(',')
85 self.parameters[name] = value.split(',')
86 else:
86 else:
87 self.parameters[name] = value
87 self.parameters[name] = value
88
88
89 def getParameters(self):
89 def getParameters(self):
90
90
91 params = {}
91 params = {}
92 for key, value in self.parameters.items():
92 for key, value in self.parameters.items():
93 s = type(value).__name__
93 s = type(value).__name__
94 if s == 'date':
94 if s == 'date':
95 params[key] = value.strftime('%Y/%m/%d')
95 params[key] = value.strftime('%Y/%m/%d')
96 elif s == 'time':
96 elif s == 'time':
97 params[key] = value.strftime('%H:%M:%S')
97 params[key] = value.strftime('%H:%M:%S')
98 else:
98 else:
99 params[key] = str(value)
99 params[key] = str(value)
100
100
101 return params
101 return params
102
102
103 def makeXml(self, element):
103 def makeXml(self, element):
104
104
105 xml = SubElement(element, self.ELEMENTNAME)
105 xml = SubElement(element, self.ELEMENTNAME)
106 for label in self.xml_labels:
106 for label in self.xml_labels:
107 xml.set(label, str(getattr(self, label)))
107 xml.set(label, str(getattr(self, label)))
108
108
109 for key, value in self.getParameters().items():
109 for key, value in self.getParameters().items():
110 xml_param = SubElement(xml, 'Parameter')
110 xml_param = SubElement(xml, 'Parameter')
111 xml_param.set('name', key)
111 xml_param.set('name', key)
112 xml_param.set('value', value)
112 xml_param.set('value', value)
113
113
114 for conf in self.operations:
114 for conf in self.operations:
115 conf.makeXml(xml)
115 conf.makeXml(xml)
116
116
117 def __str__(self):
117 def __str__(self):
118
118
119 if self.ELEMENTNAME == 'Operation':
119 if self.ELEMENTNAME == 'Operation':
120 s = ' {}[id={}]\n'.format(self.name, self.id)
120 s = ' {}[id={}]\n'.format(self.name, self.id)
121 else:
121 else:
122 s = '{}[id={}, inputId={}]\n'.format(self.name, self.id, self.inputId)
122 s = '{}[id={}, inputId={}]\n'.format(self.name, self.id, self.inputId)
123
123
124 for key, value in self.parameters.items():
124 for key, value in self.parameters.items():
125 if self.ELEMENTNAME == 'Operation':
125 if self.ELEMENTNAME == 'Operation':
126 s += ' {}: {}\n'.format(key, value)
126 s += ' {}: {}\n'.format(key, value)
127 else:
127 else:
128 s += ' {}: {}\n'.format(key, value)
128 s += ' {}: {}\n'.format(key, value)
129
129
130 for conf in self.operations:
130 for conf in self.operations:
131 s += str(conf)
131 s += str(conf)
132
132
133 return s
133 return s
134
134
135 class OperationConf(ConfBase):
135 class OperationConf(ConfBase):
136
136
137 ELEMENTNAME = 'Operation'
137 ELEMENTNAME = 'Operation'
138 xml_labels = ['id', 'name']
138 xml_labels = ['id', 'name']
139
139
140 def setup(self, id, name, priority, project_id, err_queue):
140 def setup(self, id, name, priority, project_id, err_queue):
141
141
142 self.id = str(id)
142 self.id = str(id)
143 self.project_id = project_id
143 self.project_id = project_id
144 self.name = name
144 self.name = name
145 self.type = 'other'
145 self.type = 'other'
146 self.err_queue = err_queue
146 self.err_queue = err_queue
147
147
148 def readXml(self, element, project_id, err_queue):
148 def readXml(self, element, project_id, err_queue):
149
149
150 self.id = element.get('id')
150 self.id = element.get('id')
151 self.name = element.get('name')
151 self.name = element.get('name')
152 self.type = 'other'
152 self.type = 'other'
153 self.project_id = str(project_id)
153 self.project_id = str(project_id)
154 self.err_queue = err_queue
154 self.err_queue = err_queue
155
155
156 for elm in element.iter('Parameter'):
156 for elm in element.iter('Parameter'):
157 self.addParameter(elm.get('name'), elm.get('value'))
157 self.addParameter(elm.get('name'), elm.get('value'))
158
158
159 def createObject(self):
159 def createObject(self):
160
160
161 className = eval(self.name)
161 className = eval(self.name)
162
162
163 if 'Plot' in self.name or 'Writer' in self.name or 'Send' in self.name or 'print' in self.name:
163 if 'Plot' in self.name or 'Writer' in self.name or 'Send' in self.name or 'print' in self.name:
164 kwargs = self.getKwargs()
164 kwargs = self.getKwargs()
165 opObj = className(self.name, self.id, self.project_id, self.err_queue, **kwargs)
165 opObj = className(self.name, self.id, self.project_id, self.err_queue, **kwargs)
166 opObj.start()
166 opObj.start()
167 self.type = 'external'
167 self.type = 'external'
168 else:
168 else:
169 opObj = className()
169 opObj = className()
170
170
171 self.object = opObj
171 self.object = opObj
172 return opObj
172 return opObj
173
173
174 class ProcUnitConf(ConfBase):
174 class ProcUnitConf(ConfBase):
175
175
176 ELEMENTNAME = 'ProcUnit'
176 ELEMENTNAME = 'ProcUnit'
177 xml_labels = ['id', 'inputId', 'name']
177 xml_labels = ['id', 'inputId', 'name']
178
178
179 def setup(self, project_id, id, name, datatype, inputId, err_queue):
179 def setup(self, project_id, id, name, datatype, inputId, err_queue):
180 '''
180 '''
181 '''
181 '''
182
182
183 if datatype == None and name == None:
183 if datatype == None and name == None:
184 raise ValueError('datatype or name should be defined')
184 raise ValueError('datatype or name should be defined')
185
185
186 if name == None:
186 if name == None:
187 if 'Proc' in datatype:
187 if 'Proc' in datatype:
188 name = datatype
188 name = datatype
189 else:
189 else:
190 name = '%sProc' % (datatype)
190 name = '%sProc' % (datatype)
191
191
192 if datatype == None:
192 if datatype == None:
193 datatype = name.replace('Proc', '')
193 datatype = name.replace('Proc', '')
194
194
195 self.id = str(id)
195 self.id = str(id)
196 self.project_id = project_id
196 self.project_id = project_id
197 self.name = name
197 self.name = name
198 self.datatype = datatype
198 self.datatype = datatype
199 self.inputId = inputId
199 self.inputId = inputId
200 self.err_queue = err_queue
200 self.err_queue = err_queue
201 self.operations = []
201 self.operations = []
202 self.parameters = {}
202 self.parameters = {}
203
203
204 def removeOperation(self, id):
204 def removeOperation(self, id):
205
205
206 i = [1 if x.id==id else 0 for x in self.operations]
206 i = [1 if x.id==id else 0 for x in self.operations]
207 self.operations.pop(i.index(1))
207 self.operations.pop(i.index(1))
208
208
209 def getOperation(self, id):
209 def getOperation(self, id):
210
210
211 for conf in self.operations:
211 for conf in self.operations:
212 if conf.id == id:
212 if conf.id == id:
213 return conf
213 return conf
214
214
215 def addOperation(self, name, optype='self'):
215 def addOperation(self, name, optype='self'):
216 '''
216 '''
217 '''
217 '''
218
218
219 id = self.getNewId()
219 id = self.getNewId()
220 conf = OperationConf()
220 conf = OperationConf()
221 conf.setup(id, name=name, priority='0', project_id=self.project_id, err_queue=self.err_queue)
221 conf.setup(id, name=name, priority='0', project_id=self.project_id, err_queue=self.err_queue)
222 self.operations.append(conf)
222 self.operations.append(conf)
223
223
224 return conf
224 return conf
225
225
226 def readXml(self, element, project_id, err_queue):
226 def readXml(self, element, project_id, err_queue):
227
227
228 self.id = element.get('id')
228 self.id = element.get('id')
229 self.name = element.get('name')
229 self.name = element.get('name')
230 self.inputId = None if element.get('inputId') == 'None' else element.get('inputId')
230 self.inputId = None if element.get('inputId') == 'None' else element.get('inputId')
231 self.datatype = element.get('datatype', self.name.replace(self.ELEMENTNAME.replace('Unit', ''), ''))
231 self.datatype = element.get('datatype', self.name.replace(self.ELEMENTNAME.replace('Unit', ''), ''))
232 self.project_id = str(project_id)
232 self.project_id = str(project_id)
233 self.err_queue = err_queue
233 self.err_queue = err_queue
234 self.operations = []
234 self.operations = []
235 self.parameters = {}
235 self.parameters = {}
236
236
237 for elm in element:
237 for elm in element:
238 if elm.tag == 'Parameter':
238 if elm.tag == 'Parameter':
239 self.addParameter(elm.get('name'), elm.get('value'))
239 self.addParameter(elm.get('name'), elm.get('value'))
240 elif elm.tag == 'Operation':
240 elif elm.tag == 'Operation':
241 conf = OperationConf()
241 conf = OperationConf()
242 conf.readXml(elm, project_id, err_queue)
242 conf.readXml(elm, project_id, err_queue)
243 self.operations.append(conf)
243 self.operations.append(conf)
244
244
245 def createObjects(self):
245 def createObjects(self):
246 '''
246 '''
247 Instancia de unidades de procesamiento.
247 Instancia de unidades de procesamiento.
248 '''
248 '''
249
249
250 className = eval(self.name)
250 className = eval(self.name)
251 kwargs = self.getKwargs()
251 kwargs = self.getKwargs()
252 procUnitObj = className()
252 procUnitObj = className()
253 procUnitObj.name = self.name
253 procUnitObj.name = self.name
254 log.success('creating process... id: {}, inputId: {}'.format(self.id,self.inputId), self.name)
254 log.success('creating process... id: {}, inputId: {}'.format(self.id,self.inputId), self.name)
255
255
256 for conf in self.operations:
256 for conf in self.operations:
257
257
258 opObj = conf.createObject()
258 opObj = conf.createObject()
259
259
260 log.success('adding operation: {}, type:{}'.format(conf.name,conf.type), self.name)
260 log.success('adding operation: {}, type:{}'.format(conf.name,conf.type), self.name)
261
261
262 procUnitObj.addOperation(conf, opObj)
262 procUnitObj.addOperation(conf, opObj)
263
263
264 self.object = procUnitObj
264 self.object = procUnitObj
265
265
266 def run(self):
266 def run(self):
267 '''
267 '''
268 '''
268 '''
269
269
270 return self.object.call(**self.getKwargs())
270 return self.object.call(**self.getKwargs())
271
271
272
272
273 class ReadUnitConf(ProcUnitConf):
273 class ReadUnitConf(ProcUnitConf):
274
274
275 ELEMENTNAME = 'ReadUnit'
275 ELEMENTNAME = 'ReadUnit'
276
276
277 def __init__(self):
277 def __init__(self):
278
278
279 self.id = None
279 self.id = None
280 self.datatype = None
280 self.datatype = None
281 self.name = None
281 self.name = None
282 self.inputId = None
282 self.inputId = None
283 self.operations = []
283 self.operations = []
284 self.parameters = {}
284 self.parameters = {}
285
285
286 def setup(self, project_id, id, name, datatype, err_queue, path='', startDate='', endDate='',
286 def setup(self, project_id, id, name, datatype, err_queue, path='', startDate='', endDate='',
287 startTime='', endTime='', server=None, **kwargs):
287 startTime='', endTime='', server=None, **kwargs):
288
288
289 if datatype == None and name == None:
289 if datatype == None and name == None:
290 raise ValueError('datatype or name should be defined')
290 raise ValueError('datatype or name should be defined')
291 if name == None:
291 if name == None:
292 if 'Reader' in datatype:
292 if 'Reader' in datatype:
293 name = datatype
293 name = datatype
294 datatype = name.replace('Reader','')
294 datatype = name.replace('Reader','')
295 else:
295 else:
296 name = '{}Reader'.format(datatype)
296 name = '{}Reader'.format(datatype)
297 if datatype == None:
297 if datatype == None:
298 if 'Reader' in name:
298 if 'Reader' in name:
299 datatype = name.replace('Reader','')
299 datatype = name.replace('Reader','')
300 else:
300 else:
301 datatype = name
301 datatype = name
302 name = '{}Reader'.format(name)
302 name = '{}Reader'.format(name)
303
303
304 self.id = id
304 self.id = id
305 self.project_id = project_id
305 self.project_id = project_id
306 self.name = name
306 self.name = name
307 self.datatype = datatype
307 self.datatype = datatype
308 self.err_queue = err_queue
308 self.err_queue = err_queue
309
309
310 self.addParameter(name='path', value=path)
310 self.addParameter(name='path', value=path)
311 self.addParameter(name='startDate', value=startDate)
311 self.addParameter(name='startDate', value=startDate)
312 self.addParameter(name='endDate', value=endDate)
312 self.addParameter(name='endDate', value=endDate)
313 self.addParameter(name='startTime', value=startTime)
313 self.addParameter(name='startTime', value=startTime)
314 self.addParameter(name='endTime', value=endTime)
314 self.addParameter(name='endTime', value=endTime)
315
315
316 for key, value in kwargs.items():
316 for key, value in kwargs.items():
317 self.addParameter(name=key, value=value)
317 self.addParameter(name=key, value=value)
318
318
319
319
320 class Project(Process):
320 class Project(Process):
321 """API to create signal chain projects"""
321 """API to create signal chain projects"""
322
322
323 ELEMENTNAME = 'Project'
323 ELEMENTNAME = 'Project'
324
324
325 def __init__(self, name=''):
325 def __init__(self, name=''):
326
326
327 Process.__init__(self)
327 Process.__init__(self)
328 self.id = '1'
328 self.id = '1'
329 if name:
329 if name:
330 self.name = '{} ({})'.format(Process.__name__, name)
330 self.name = '{} ({})'.format(Process.__name__, name)
331 self.filename = None
331 self.filename = None
332 self.description = None
332 self.description = None
333 self.email = None
333 self.email = None
334 self.alarm = []
334 self.alarm = []
335 self.configurations = {}
335 self.configurations = {}
336 # self.err_queue = Queue()
336 # self.err_queue = Queue()
337 self.err_queue = None
337 self.err_queue = None
338 self.started = False
338 self.started = False
339
339
340 def getNewId(self):
340 def getNewId(self):
341
341
342 idList = list(self.configurations.keys())
342 idList = list(self.configurations.keys())
343 id = int(self.id) * 10
343 id = int(self.id) * 10
344
344
345 while True:
345 while True:
346 id += 1
346 id += 1
347
347
348 if str(id) in idList:
348 if str(id) in idList:
349 continue
349 continue
350
350
351 break
351 break
352
352
353 return str(id)
353 return str(id)
354
354
355 def updateId(self, new_id):
355 def updateId(self, new_id):
356
356
357 self.id = str(new_id)
357 self.id = str(new_id)
358
358
359 keyList = list(self.configurations.keys())
359 keyList = list(self.configurations.keys())
360 keyList.sort()
360 keyList.sort()
361
361
362 n = 1
362 n = 1
363 new_confs = {}
363 new_confs = {}
364
364
365 for procKey in keyList:
365 for procKey in keyList:
366
366
367 conf = self.configurations[procKey]
367 conf = self.configurations[procKey]
368 idProcUnit = str(int(self.id) * 10 + n)
368 idProcUnit = str(int(self.id) * 10 + n)
369 conf.updateId(idProcUnit)
369 conf.updateId(idProcUnit)
370 new_confs[idProcUnit] = conf
370 new_confs[idProcUnit] = conf
371 n += 1
371 n += 1
372
372
373 self.configurations = new_confs
373 self.configurations = new_confs
374
374
375 def setup(self, id=1, name='', description='', email=None, alarm=[]):
375 def setup(self, id=1, name='', description='', email=None, alarm=[]):
376
376
377 self.id = str(id)
377 self.id = str(id)
378 self.description = description
378 self.description = description
379 self.email = email
379 self.email = email
380 self.alarm = alarm
380 self.alarm = alarm
381 if name:
381 if name:
382 self.name = '{} ({})'.format(Process.__name__, name)
382 self.name = '{} ({})'.format(Process.__name__, name)
383
383
384 def update(self, **kwargs):
384 def update(self, **kwargs):
385
385
386 for key, value in kwargs.items():
386 for key, value in kwargs.items():
387 setattr(self, key, value)
387 setattr(self, key, value)
388
388
389 def clone(self):
389 def clone(self):
390
390
391 p = Project()
391 p = Project()
392 p.id = self.id
392 p.id = self.id
393 p.name = self.name
393 p.name = self.name
394 p.description = self.description
394 p.description = self.description
395 p.configurations = self.configurations.copy()
395 p.configurations = self.configurations.copy()
396
396
397 return p
397 return p
398
398
399 def addReadUnit(self, id=None, datatype=None, name=None, **kwargs):
399 def addReadUnit(self, id=None, datatype=None, name=None, **kwargs):
400
400
401 '''
401 '''
402 '''
402 '''
403
403
404 if id is None:
404 if id is None:
405 idReadUnit = self.getNewId()
405 idReadUnit = self.getNewId()
406 else:
406 else:
407 idReadUnit = str(id)
407 idReadUnit = str(id)
408
408
409 conf = ReadUnitConf()
409 conf = ReadUnitConf()
410 conf.setup(self.id, idReadUnit, name, datatype, self.err_queue, **kwargs)
410 conf.setup(self.id, idReadUnit, name, datatype, self.err_queue, **kwargs)
411 self.configurations[conf.id] = conf
411 self.configurations[conf.id] = conf
412
412
413 return conf
413 return conf
414
414
415 def addProcUnit(self, id=None, inputId='0', datatype=None, name=None):
415 def addProcUnit(self, id=None, inputId='0', datatype=None, name=None):
416
416
417 '''
417 '''
418 '''
418 '''
419
419
420 if id is None:
420 if id is None:
421 idProcUnit = self.getNewId()
421 idProcUnit = self.getNewId()
422 else:
422 else:
423 idProcUnit = id
423 idProcUnit = id
424
424
425 conf = ProcUnitConf()
425 conf = ProcUnitConf()
426 conf.setup(self.id, idProcUnit, name, datatype, inputId, self.err_queue)
426 conf.setup(self.id, idProcUnit, name, datatype, inputId, self.err_queue)
427 self.configurations[conf.id] = conf
427 self.configurations[conf.id] = conf
428
428
429 return conf
429 return conf
430
430
431 def removeProcUnit(self, id):
431 def removeProcUnit(self, id):
432
432
433 if id in self.configurations:
433 if id in self.configurations:
434 self.configurations.pop(id)
434 self.configurations.pop(id)
435
435
436 def getReadUnit(self):
436 def getReadUnit(self):
437
437
438 for obj in list(self.configurations.values()):
438 for obj in list(self.configurations.values()):
439 if obj.ELEMENTNAME == 'ReadUnit':
439 if obj.ELEMENTNAME == 'ReadUnit':
440 return obj
440 return obj
441
441
442 return None
442 return None
443
443
444 def getProcUnit(self, id):
444 def getProcUnit(self, id):
445
445
446 return self.configurations[id]
446 return self.configurations[id]
447
447
448 def getUnits(self):
448 def getUnits(self):
449
449
450 keys = list(self.configurations)
450 keys = list(self.configurations)
451 keys.sort()
451 keys.sort()
452
452
453 for key in keys:
453 for key in keys:
454 yield self.configurations[key]
454 yield self.configurations[key]
455
455
456 def updateUnit(self, id, **kwargs):
456 def updateUnit(self, id, **kwargs):
457
457
458 conf = self.configurations[id].update(**kwargs)
458 conf = self.configurations[id].update(**kwargs)
459
459
460 def makeXml(self):
460 def makeXml(self):
461
461
462 xml = Element('Project')
462 xml = Element('Project')
463 xml.set('id', str(self.id))
463 xml.set('id', str(self.id))
464 xml.set('name', self.name)
464 xml.set('name', self.name)
465 xml.set('description', self.description)
465 xml.set('description', self.description)
466
466
467 for conf in self.configurations.values():
467 for conf in self.configurations.values():
468 conf.makeXml(xml)
468 conf.makeXml(xml)
469
469
470 self.xml = xml
470 self.xml = xml
471
471
472 def writeXml(self, filename=None):
472 def writeXml(self, filename=None):
473
473
474 if filename == None:
474 if filename == None:
475 if self.filename:
475 if self.filename:
476 filename = self.filename
476 filename = self.filename
477 else:
477 else:
478 filename = 'schain.xml'
478 filename = 'schain.xml'
479
479
480 if not filename:
480 if not filename:
481 print('filename has not been defined. Use setFilename(filename) for do it.')
481 print('filename has not been defined. Use setFilename(filename) for do it.')
482 return 0
482 return 0
483
483
484 abs_file = os.path.abspath(filename)
484 abs_file = os.path.abspath(filename)
485
485
486 if not os.access(os.path.dirname(abs_file), os.W_OK):
486 if not os.access(os.path.dirname(abs_file), os.W_OK):
487 print('No write permission on %s' % os.path.dirname(abs_file))
487 print('No write permission on %s' % os.path.dirname(abs_file))
488 return 0
488 return 0
489
489
490 if os.path.isfile(abs_file) and not(os.access(abs_file, os.W_OK)):
490 if os.path.isfile(abs_file) and not(os.access(abs_file, os.W_OK)):
491 print('File %s already exists and it could not be overwriten' % abs_file)
491 print('File %s already exists and it could not be overwriten' % abs_file)
492 return 0
492 return 0
493
493
494 self.makeXml()
494 self.makeXml()
495
495
496 ElementTree(self.xml).write(abs_file, method='xml')
496 ElementTree(self.xml).write(abs_file, method='xml')
497
497
498 self.filename = abs_file
498 self.filename = abs_file
499
499
500 return 1
500 return 1
501
501
502 def readXml(self, filename):
502 def readXml(self, filename):
503
503
504 abs_file = os.path.abspath(filename)
504 abs_file = os.path.abspath(filename)
505
505
506 self.configurations = {}
506 self.configurations = {}
507
507
508 try:
508 try:
509 self.xml = ElementTree().parse(abs_file)
509 self.xml = ElementTree().parse(abs_file)
510 except:
510 except:
511 log.error('Error reading %s, verify file format' % filename)
511 log.error('Error reading %s, verify file format' % filename)
512 return 0
512 return 0
513
513
514 self.id = self.xml.get('id')
514 self.id = self.xml.get('id')
515 self.name = self.xml.get('name')
515 self.name = self.xml.get('name')
516 self.description = self.xml.get('description')
516 self.description = self.xml.get('description')
517
517
518 for element in self.xml:
518 for element in self.xml:
519 if element.tag == 'ReadUnit':
519 if element.tag == 'ReadUnit':
520 conf = ReadUnitConf()
520 conf = ReadUnitConf()
521 conf.readXml(element, self.id, self.err_queue)
521 conf.readXml(element, self.id, self.err_queue)
522 self.configurations[conf.id] = conf
522 self.configurations[conf.id] = conf
523 elif element.tag == 'ProcUnit':
523 elif element.tag == 'ProcUnit':
524 conf = ProcUnitConf()
524 conf = ProcUnitConf()
525 input_proc = self.configurations[element.get('inputId')]
525 input_proc = self.configurations[element.get('inputId')]
526 conf.readXml(element, self.id, self.err_queue)
526 conf.readXml(element, self.id, self.err_queue)
527 self.configurations[conf.id] = conf
527 self.configurations[conf.id] = conf
528
528
529 self.filename = abs_file
529 self.filename = abs_file
530
530
531 return 1
531 return 1
532
532
533 def __str__(self):
533 def __str__(self):
534
534
535 text = '\nProject[id=%s, name=%s, description=%s]\n\n' % (
535 text = '\nProject[id=%s, name=%s, description=%s]\n\n' % (
536 self.id,
536 self.id,
537 self.name,
537 self.name,
538 self.description,
538 self.description,
539 )
539 )
540
540
541 for conf in self.configurations.values():
541 for conf in self.configurations.values():
542 text += '{}'.format(conf)
542 text += '{}'.format(conf)
543
543
544 return text
544 return text
545
545
546 def createObjects(self):
546 def createObjects(self):
547
547
548 keys = list(self.configurations.keys())
548 keys = list(self.configurations.keys())
549 keys.sort()
549 keys.sort()
550 for key in keys:
550 for key in keys:
551 conf = self.configurations[key]
551 conf = self.configurations[key]
552 conf.createObjects()
552 conf.createObjects()
553 if conf.inputId is not None:
553 if conf.inputId is not None:
554 conf.object.setInput(self.configurations[conf.inputId].object)
554 conf.object.setInput(self.configurations[conf.inputId].object)
555
555
556 def monitor(self):
556 def monitor(self):
557
557
558 t = Thread(target=self._monitor, args=(self.err_queue, self.ctx))
558 t = Thread(target=self._monitor, args=(self.err_queue, self.ctx))
559 t.start()
559 t.start()
560
560
561 def _monitor(self, queue, ctx):
561 def _monitor(self, queue, ctx):
562
562
563 import socket
563 import socket
564
564
565 procs = 0
565 procs = 0
566 err_msg = ''
566 err_msg = ''
567
567
568 while True:
568 while True:
569 msg = queue.get()
569 msg = queue.get()
570 if '#_start_#' in msg:
570 if '#_start_#' in msg:
571 procs += 1
571 procs += 1
572 elif '#_end_#' in msg:
572 elif '#_end_#' in msg:
573 procs -=1
573 procs -=1
574 else:
574 else:
575 err_msg = msg
575 err_msg = msg
576
576
577 if procs == 0 or 'Traceback' in err_msg:
577 if procs == 0 or 'Traceback' in err_msg:
578 break
578 break
579 time.sleep(0.1)
579 time.sleep(0.1)
580
580
581 if '|' in err_msg:
581 if '|' in err_msg:
582 name, err = err_msg.split('|')
582 name, err = err_msg.split('|')
583 if 'SchainWarning' in err:
583 if 'SchainWarning' in err:
584 log.warning(err.split('SchainWarning:')[-1].split('\n')[0].strip(), name)
584 log.warning(err.split('SchainWarning:')[-1].split('\n')[0].strip(), name)
585 elif 'SchainError' in err:
585 elif 'SchainError' in err:
586 log.error(err.split('SchainError:')[-1].split('\n')[0].strip(), name)
586 log.error(err.split('SchainError:')[-1].split('\n')[0].strip(), name)
587 else:
587 else:
588 log.error(err, name)
588 log.error(err, name)
589 else:
589 else:
590 name, err = self.name, err_msg
590 name, err = self.name, err_msg
591
591
592 time.sleep(1)
592 time.sleep(1)
593
593
594 ctx.term()
594 ctx.term()
595
595
596 message = ''.join(err)
596 message = ''.join(err)
597
597
598 if err_msg:
598 if err_msg:
599 subject = 'SChain v%s: Error running %s\n' % (
599 subject = 'SChain v%s: Error running %s\n' % (
600 schainpy.__version__, self.name)
600 schainpy.__version__, self.name)
601
601
602 subtitle = 'Hostname: %s\n' % socket.gethostbyname(
602 subtitle = 'Hostname: %s\n' % socket.gethostbyname(
603 socket.gethostname())
603 socket.gethostname())
604 subtitle += 'Working directory: %s\n' % os.path.abspath('./')
604 subtitle += 'Working directory: %s\n' % os.path.abspath('./')
605 subtitle += 'Configuration file: %s\n' % self.filename
605 subtitle += 'Configuration file: %s\n' % self.filename
606 subtitle += 'Time: %s\n' % str(datetime.datetime.now())
606 subtitle += 'Time: %s\n' % str(datetime.datetime.now())
607
607
608 readUnitConfObj = self.getReadUnit()
608 readUnitConfObj = self.getReadUnit()
609 if readUnitConfObj:
609 if readUnitConfObj:
610 subtitle += '\nInput parameters:\n'
610 subtitle += '\nInput parameters:\n'
611 subtitle += '[Data path = %s]\n' % readUnitConfObj.parameters['path']
611 subtitle += '[Data path = %s]\n' % readUnitConfObj.parameters['path']
612 subtitle += '[Start date = %s]\n' % readUnitConfObj.parameters['startDate']
612 subtitle += '[Start date = %s]\n' % readUnitConfObj.parameters['startDate']
613 subtitle += '[End date = %s]\n' % readUnitConfObj.parameters['endDate']
613 subtitle += '[End date = %s]\n' % readUnitConfObj.parameters['endDate']
614 subtitle += '[Start time = %s]\n' % readUnitConfObj.parameters['startTime']
614 subtitle += '[Start time = %s]\n' % readUnitConfObj.parameters['startTime']
615 subtitle += '[End time = %s]\n' % readUnitConfObj.parameters['endTime']
615 subtitle += '[End time = %s]\n' % readUnitConfObj.parameters['endTime']
616
616
617 a = Alarm(
617 a = Alarm(
618 modes=self.alarm,
618 modes=self.alarm,
619 email=self.email,
619 email=self.email,
620 message=message,
620 message=message,
621 subject=subject,
621 subject=subject,
622 subtitle=subtitle,
622 subtitle=subtitle,
623 filename=self.filename
623 filename=self.filename
624 )
624 )
625
625
626 a.start()
626 a.start()
627
627
628 def setFilename(self, filename):
628 def setFilename(self, filename):
629
629
630 self.filename = filename
630 self.filename = filename
631
631
632 def runProcs(self):
632 def runProcs(self):
633
633
634 err = False
634 err = False
635 n = len(self.configurations)
635 n = len(self.configurations)
636
636
637 while not err:
637 while not err:
638 #print("STAR")
638 for conf in self.getUnits():
639 for conf in self.getUnits():
640 #print("CONF: ",conf)
639 ok = conf.run()
641 ok = conf.run()
640 if ok == 'Error':
642 if ok == 'Error':
641 n -= 1
643 n -= 1
642 continue
644 continue
643 elif not ok:
645 elif not ok:
644 break
646 break
645 if n == 0:
647 if n == 0:
646 err = True
648 err = True
647
649
648 def run(self):
650 def run(self):
649
651
650 log.success('\nStarting Project {} [id={}]'.format(self.name, self.id), tag='')
652 log.success('\nStarting Project {} [id={}]'.format(self.name, self.id), tag='')
651 self.started = True
653 self.started = True
652 self.start_time = time.time()
654 self.start_time = time.time()
653 self.createObjects()
655 self.createObjects()
654 self.runProcs()
656 self.runProcs()
655 log.success('{} Done (Time: {:4.2f}s)'.format(
657 log.success('{} Done (Time: {:4.2f}s)'.format(
656 self.name,
658 self.name,
657 time.time()-self.start_time), '')
659 time.time()-self.start_time), '')
@@ -1,209 +1,212
1 '''
1 '''
2 Base clases to create Processing units and operations, the MPDecorator
2 Base clases to create Processing units and operations, the MPDecorator
3 must be used in plotting and writing operations to allow to run as an
3 must be used in plotting and writing operations to allow to run as an
4 external process.
4 external process.
5 '''
5 '''
6 import os
6 import os
7 import inspect
7 import inspect
8 import zmq
8 import zmq
9 import time
9 import time
10 import pickle
10 import pickle
11 import traceback
11 import traceback
12 from threading import Thread
12 from threading import Thread
13 from multiprocessing import Process, Queue
13 from multiprocessing import Process, Queue
14 from schainpy.utils import log
14 from schainpy.utils import log
15 import copy
15 import copy
16 QUEUE_SIZE = int(os.environ.get('QUEUE_MAX_SIZE', '10'))
16 QUEUE_SIZE = int(os.environ.get('QUEUE_MAX_SIZE', '10'))
17 class ProcessingUnit(object):
17 class ProcessingUnit(object):
18 '''
18 '''
19 Base class to create Signal Chain Units
19 Base class to create Signal Chain Units
20 '''
20 '''
21
21
22 proc_type = 'processing'
22 proc_type = 'processing'
23
23
24 def __init__(self):
24 def __init__(self):
25
25
26 self.dataIn = None
26 self.dataIn = None
27 self.dataOut = None
27 self.dataOut = None
28 self.isConfig = False
28 self.isConfig = False
29 self.operations = []
29 self.operations = []
30
30
31 def setInput(self, unit):
31 def setInput(self, unit):
32
32
33 self.dataIn = unit.dataOut
33 self.dataIn = unit.dataOut
34
34
35
35
36 def getAllowedArgs(self):
36 def getAllowedArgs(self):
37 if hasattr(self, '__attrs__'):
37 if hasattr(self, '__attrs__'):
38 return self.__attrs__
38 return self.__attrs__
39 else:
39 else:
40 return inspect.getargspec(self.run).args
40 return inspect.getargspec(self.run).args
41
41
42 def addOperation(self, conf, operation):
42 def addOperation(self, conf, operation):
43 '''
43 '''
44 '''
44 '''
45
45
46 self.operations.append((operation, conf.type, conf.getKwargs()))
46 self.operations.append((operation, conf.type, conf.getKwargs()))
47
47
48 def getOperationObj(self, objId):
48 def getOperationObj(self, objId):
49
49
50 if objId not in list(self.operations.keys()):
50 if objId not in list(self.operations.keys()):
51 return None
51 return None
52
52
53 return self.operations[objId]
53 return self.operations[objId]
54
54
55 def call(self, **kwargs):
55 def call(self, **kwargs):
56 '''
56 '''
57 '''
57 '''
58
58
59 try:
59 try:
60 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
60 # if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
61 return self.dataIn.isReady()
61 # return self.dataIn.isReady()
62 elif self.dataIn is None or not self.dataIn.error:
62 #dataIn=None es unidades de Lectura, segunda parte unidades de procesamiento
63 if self.dataIn is None or (not self.dataIn.error and not self.dataIn.flagNoData):
63 self.run(**kwargs)
64 self.run(**kwargs)
64 elif self.dataIn.error:
65 elif self.dataIn.error:
65 self.dataOut.error = self.dataIn.error
66 self.dataOut.error = self.dataIn.error
66 self.dataOut.flagNoData = True
67 self.dataOut.flagNoData = True
68
67 except:
69 except:
68
70
69 err = traceback.format_exc()
71 err = traceback.format_exc()
70 if 'SchainWarning' in err:
72 if 'SchainWarning' in err:
71 log.warning(err.split('SchainWarning:')[-1].split('\n')[0].strip(), self.name)
73 log.warning(err.split('SchainWarning:')[-1].split('\n')[0].strip(), self.name)
72 elif 'SchainError' in err:
74 elif 'SchainError' in err:
73 log.error(err.split('SchainError:')[-1].split('\n')[0].strip(), self.name)
75 log.error(err.split('SchainError:')[-1].split('\n')[0].strip(), self.name)
74 else:
76 else:
75 log.error(err, self.name)
77 log.error(err, self.name)
76 self.dataOut.error = True
78 self.dataOut.error = True
77
79
78 for op, optype, opkwargs in self.operations:
80 for op, optype, opkwargs in self.operations:
79 if optype == 'other' and not self.dataOut.flagNoData:
81 if optype == 'other' and self.dataOut.isReady():
80 try:
82 try:
81 self.dataOut = op.run(self.dataOut, **opkwargs)
83 self.dataOut = op.run(self.dataOut, **opkwargs)
82 except Exception as e:
84 except Exception as e:
85 print(e)
83 self.dataOut.error = True
86 self.dataOut.error = True
84 return 'Error'
87 return 'Error'
85 elif optype == 'external' and self.dataOut.isReady() :
88 elif optype == 'external' and self.dataOut.isReady() :
86 op.queue.put(copy.deepcopy(self.dataOut))
89 op.queue.put(copy.deepcopy(self.dataOut))
87 elif optype == 'external' and self.dataOut.error:
90 elif optype == 'external' and self.dataOut.error:
88 op.queue.put(self.dataOut)
91 op.queue.put(copy.deepcopy(self.dataOut))
89
92
90 return 'Error' if self.dataOut.error else self.dataOut.isReady()
93 return 'Error' if self.dataOut.error else True#self.dataOut.isReady()
91
94
92 def setup(self):
95 def setup(self):
93
96
94 raise NotImplementedError
97 raise NotImplementedError
95
98
96 def run(self):
99 def run(self):
97
100
98 raise NotImplementedError
101 raise NotImplementedError
99
102
100 def close(self):
103 def close(self):
101
104
102 return
105 return
103
106
104
107
105 class Operation(object):
108 class Operation(object):
106
109
107 '''
110 '''
108 '''
111 '''
109
112
110 proc_type = 'operation'
113 proc_type = 'operation'
111
114
112 def __init__(self):
115 def __init__(self):
113
116
114 self.id = None
117 self.id = None
115 self.isConfig = False
118 self.isConfig = False
116
119
117 if not hasattr(self, 'name'):
120 if not hasattr(self, 'name'):
118 self.name = self.__class__.__name__
121 self.name = self.__class__.__name__
119
122
120 def getAllowedArgs(self):
123 def getAllowedArgs(self):
121 if hasattr(self, '__attrs__'):
124 if hasattr(self, '__attrs__'):
122 return self.__attrs__
125 return self.__attrs__
123 else:
126 else:
124 return inspect.getargspec(self.run).args
127 return inspect.getargspec(self.run).args
125
128
126 def setup(self):
129 def setup(self):
127
130
128 self.isConfig = True
131 self.isConfig = True
129
132
130 raise NotImplementedError
133 raise NotImplementedError
131
134
132 def run(self, dataIn, **kwargs):
135 def run(self, dataIn, **kwargs):
133 """
136 """
134 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los
137 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los
135 atributos del objeto dataIn.
138 atributos del objeto dataIn.
136
139
137 Input:
140 Input:
138
141
139 dataIn : objeto del tipo JROData
142 dataIn : objeto del tipo JROData
140
143
141 Return:
144 Return:
142
145
143 None
146 None
144
147
145 Affected:
148 Affected:
146 __buffer : buffer de recepcion de datos.
149 __buffer : buffer de recepcion de datos.
147
150
148 """
151 """
149 if not self.isConfig:
152 if not self.isConfig:
150 self.setup(**kwargs)
153 self.setup(**kwargs)
151
154
152 raise NotImplementedError
155 raise NotImplementedError
153
156
154 def close(self):
157 def close(self):
155
158
156 return
159 return
157
160
158
161
159 def MPDecorator(BaseClass):
162 def MPDecorator(BaseClass):
160 """
163 """
161 Multiprocessing class decorator
164 Multiprocessing class decorator
162
165
163 This function add multiprocessing features to a BaseClass.
166 This function add multiprocessing features to a BaseClass.
164 """
167 """
165
168
166 class MPClass(BaseClass, Process):
169 class MPClass(BaseClass, Process):
167
170
168 def __init__(self, *args, **kwargs):
171 def __init__(self, *args, **kwargs):
169 super(MPClass, self).__init__()
172 super(MPClass, self).__init__()
170 Process.__init__(self)
173 Process.__init__(self)
171
174
172 self.args = args
175 self.args = args
173 self.kwargs = kwargs
176 self.kwargs = kwargs
174 self.t = time.time()
177 self.t = time.time()
175 self.op_type = 'external'
178 self.op_type = 'external'
176 self.name = BaseClass.__name__
179 self.name = BaseClass.__name__
177 self.__doc__ = BaseClass.__doc__
180 self.__doc__ = BaseClass.__doc__
178
181
179 if 'plot' in self.name.lower() and not self.name.endswith('_'):
182 if 'plot' in self.name.lower() and not self.name.endswith('_'):
180 self.name = '{}{}'.format(self.CODE.upper(), 'Plot')
183 self.name = '{}{}'.format(self.CODE.upper(), 'Plot')
181
184
182 self.start_time = time.time()
185 self.start_time = time.time()
183 self.err_queue = args[3]
186 self.err_queue = args[3]
184 self.queue = Queue(maxsize=QUEUE_SIZE)
187 self.queue = Queue(maxsize=QUEUE_SIZE)
185 self.myrun = BaseClass.run
188 self.myrun = BaseClass.run
186
189
187 def run(self):
190 def run(self):
188
191
189 while True:
192 while True:
190
193
191 dataOut = self.queue.get()
194 dataOut = self.queue.get()
192
195
193 if not dataOut.error:
196 if not dataOut.error:
194 try:
197 try:
195 BaseClass.run(self, dataOut, **self.kwargs)
198 BaseClass.run(self, dataOut, **self.kwargs)
196 except:
199 except:
197 err = traceback.format_exc()
200 err = traceback.format_exc()
198 log.error(err, self.name)
201 log.error(err, self.name)
199 else:
202 else:
200 break
203 break
201
204
202 self.close()
205 self.close()
203
206
204 def close(self):
207 def close(self):
205
208
206 BaseClass.close(self)
209 BaseClass.close(self)
207 log.success('Done...(Time:{:4.2f} secs)'.format(time.time()-self.start_time), self.name)
210 log.success('Done...(Time:{:4.2f} secs)'.format(time.time()-self.start_time), self.name)
208
211
209 return MPClass
212 return MPClass
@@ -1,1689 +1,1689
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 # All rights reserved.
2 # All rights reserved.
3 #
3 #
4 # Distributed under the terms of the BSD 3-clause license.
4 # Distributed under the terms of the BSD 3-clause license.
5 """Spectra processing Unit and operations
5 """Spectra processing Unit and operations
6
6
7 Here you will find the processing unit `SpectraProc` and several operations
7 Here you will find the processing unit `SpectraProc` and several operations
8 to work with Spectra data type
8 to work with Spectra data type
9 """
9 """
10
10
11 import time
11 import time
12 import itertools
12 import itertools
13
13
14 import numpy
14 import numpy
15 import math
15 import math
16
16
17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
18 from schainpy.model.data.jrodata import Spectra
18 from schainpy.model.data.jrodata import Spectra
19 from schainpy.model.data.jrodata import hildebrand_sekhon
19 from schainpy.model.data.jrodata import hildebrand_sekhon
20 from schainpy.utils import log
20 from schainpy.utils import log
21
21
22 from scipy.optimize import curve_fit
22 from scipy.optimize import curve_fit
23
23
24 class SpectraProc(ProcessingUnit):
24 class SpectraProc(ProcessingUnit):
25
25
26 def __init__(self):
26 def __init__(self):
27
27
28 ProcessingUnit.__init__(self)
28 ProcessingUnit.__init__(self)
29
29
30 self.buffer = None
30 self.buffer = None
31 self.firstdatatime = None
31 self.firstdatatime = None
32 self.profIndex = 0
32 self.profIndex = 0
33 self.dataOut = Spectra()
33 self.dataOut = Spectra()
34 self.id_min = None
34 self.id_min = None
35 self.id_max = None
35 self.id_max = None
36 self.setupReq = False #Agregar a todas las unidades de proc
36 self.setupReq = False #Agregar a todas las unidades de proc
37
37
38 def __updateSpecFromVoltage(self):
38 def __updateSpecFromVoltage(self):
39
39
40 self.dataOut.timeZone = self.dataIn.timeZone
40 self.dataOut.timeZone = self.dataIn.timeZone
41 self.dataOut.dstFlag = self.dataIn.dstFlag
41 self.dataOut.dstFlag = self.dataIn.dstFlag
42 self.dataOut.errorCount = self.dataIn.errorCount
42 self.dataOut.errorCount = self.dataIn.errorCount
43 self.dataOut.useLocalTime = self.dataIn.useLocalTime
43 self.dataOut.useLocalTime = self.dataIn.useLocalTime
44 try:
44 try:
45 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
45 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
46 except:
46 except:
47 pass
47 pass
48 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
48 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
49 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
49 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
50 self.dataOut.channelList = self.dataIn.channelList
50 self.dataOut.channelList = self.dataIn.channelList
51 self.dataOut.heightList = self.dataIn.heightList
51 self.dataOut.heightList = self.dataIn.heightList
52 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
52 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
53 self.dataOut.nProfiles = self.dataOut.nFFTPoints
53 self.dataOut.nProfiles = self.dataOut.nFFTPoints
54 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
54 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
55 self.dataOut.utctime = self.firstdatatime
55 self.dataOut.utctime = self.firstdatatime
56 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
56 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
57 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
57 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
58 self.dataOut.flagShiftFFT = False
58 self.dataOut.flagShiftFFT = False
59 self.dataOut.nCohInt = self.dataIn.nCohInt
59 self.dataOut.nCohInt = self.dataIn.nCohInt
60 self.dataOut.nIncohInt = 1
60 self.dataOut.nIncohInt = 1
61 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
61 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
62 self.dataOut.frequency = self.dataIn.frequency
62 self.dataOut.frequency = self.dataIn.frequency
63 self.dataOut.realtime = self.dataIn.realtime
63 self.dataOut.realtime = self.dataIn.realtime
64 self.dataOut.azimuth = self.dataIn.azimuth
64 self.dataOut.azimuth = self.dataIn.azimuth
65 self.dataOut.zenith = self.dataIn.zenith
65 self.dataOut.zenith = self.dataIn.zenith
66 self.dataOut.codeList = self.dataIn.codeList
66 self.dataOut.codeList = self.dataIn.codeList
67 self.dataOut.azimuthList = self.dataIn.azimuthList
67 self.dataOut.azimuthList = self.dataIn.azimuthList
68 self.dataOut.elevationList = self.dataIn.elevationList
68 self.dataOut.elevationList = self.dataIn.elevationList
69
69
70
70
71
71
72 def __getFft(self):
72 def __getFft(self):
73 """
73 """
74 Convierte valores de Voltaje a Spectra
74 Convierte valores de Voltaje a Spectra
75
75
76 Affected:
76 Affected:
77 self.dataOut.data_spc
77 self.dataOut.data_spc
78 self.dataOut.data_cspc
78 self.dataOut.data_cspc
79 self.dataOut.data_dc
79 self.dataOut.data_dc
80 self.dataOut.heightList
80 self.dataOut.heightList
81 self.profIndex
81 self.profIndex
82 self.buffer
82 self.buffer
83 self.dataOut.flagNoData
83 self.dataOut.flagNoData
84 """
84 """
85 fft_volt = numpy.fft.fft(
85 fft_volt = numpy.fft.fft(
86 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
86 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
87 fft_volt = fft_volt.astype(numpy.dtype('complex'))
87 fft_volt = fft_volt.astype(numpy.dtype('complex'))
88 dc = fft_volt[:, 0, :]
88 dc = fft_volt[:, 0, :]
89
89
90 # calculo de self-spectra
90 # calculo de self-spectra
91 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
91 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
92 spc = fft_volt * numpy.conjugate(fft_volt)
92 spc = fft_volt * numpy.conjugate(fft_volt)
93 spc = spc.real
93 spc = spc.real
94
94
95 blocksize = 0
95 blocksize = 0
96 blocksize += dc.size
96 blocksize += dc.size
97 blocksize += spc.size
97 blocksize += spc.size
98
98
99 cspc = None
99 cspc = None
100 pairIndex = 0
100 pairIndex = 0
101 if self.dataOut.pairsList != None:
101 if self.dataOut.pairsList != None:
102 # calculo de cross-spectra
102 # calculo de cross-spectra
103 cspc = numpy.zeros(
103 cspc = numpy.zeros(
104 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
104 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
105 for pair in self.dataOut.pairsList:
105 for pair in self.dataOut.pairsList:
106 if pair[0] not in self.dataOut.channelList:
106 if pair[0] not in self.dataOut.channelList:
107 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
107 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
108 str(pair), str(self.dataOut.channelList)))
108 str(pair), str(self.dataOut.channelList)))
109 if pair[1] not in self.dataOut.channelList:
109 if pair[1] not in self.dataOut.channelList:
110 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
110 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
111 str(pair), str(self.dataOut.channelList)))
111 str(pair), str(self.dataOut.channelList)))
112
112
113 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
113 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
114 numpy.conjugate(fft_volt[pair[1], :, :])
114 numpy.conjugate(fft_volt[pair[1], :, :])
115 pairIndex += 1
115 pairIndex += 1
116 blocksize += cspc.size
116 blocksize += cspc.size
117
117
118 self.dataOut.data_spc = spc
118 self.dataOut.data_spc = spc
119 self.dataOut.data_cspc = cspc
119 self.dataOut.data_cspc = cspc
120 self.dataOut.data_dc = dc
120 self.dataOut.data_dc = dc
121 self.dataOut.blockSize = blocksize
121 self.dataOut.blockSize = blocksize
122 self.dataOut.flagShiftFFT = False
122 self.dataOut.flagShiftFFT = False
123
123
124 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
124 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
125
125
126 if self.dataIn.type == "Spectra":
126 if self.dataIn.type == "Spectra":
127
127
128 try:
128 try:
129 self.dataOut.copy(self.dataIn)
129 self.dataOut.copy(self.dataIn)
130
130
131 except Exception as e:
131 except Exception as e:
132 print(e)
132 print(e)
133
133
134 if shift_fft:
134 if shift_fft:
135 #desplaza a la derecha en el eje 2 determinadas posiciones
135 #desplaza a la derecha en el eje 2 determinadas posiciones
136 shift = int(self.dataOut.nFFTPoints/2)
136 shift = int(self.dataOut.nFFTPoints/2)
137 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
137 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
138
138
139 if self.dataOut.data_cspc is not None:
139 if self.dataOut.data_cspc is not None:
140 #desplaza a la derecha en el eje 2 determinadas posiciones
140 #desplaza a la derecha en el eje 2 determinadas posiciones
141 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
141 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
142 if pairsList:
142 if pairsList:
143 self.__selectPairs(pairsList)
143 self.__selectPairs(pairsList)
144
144
145
145
146 elif self.dataIn.type == "Voltage":
146 elif self.dataIn.type == "Voltage":
147
147
148 self.dataOut.flagNoData = True
148 self.dataOut.flagNoData = True
149
149
150 if nFFTPoints == None:
150 if nFFTPoints == None:
151 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
151 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
152
152
153 if nProfiles == None:
153 if nProfiles == None:
154 nProfiles = nFFTPoints
154 nProfiles = nFFTPoints
155
155
156 if ippFactor == None:
156 if ippFactor == None:
157 self.dataOut.ippFactor = 1
157 self.dataOut.ippFactor = 1
158
158
159 self.dataOut.nFFTPoints = nFFTPoints
159 self.dataOut.nFFTPoints = nFFTPoints
160
160
161 if self.buffer is None:
161 if self.buffer is None:
162 self.buffer = numpy.zeros((self.dataIn.nChannels,
162 self.buffer = numpy.zeros((self.dataIn.nChannels,
163 nProfiles,
163 nProfiles,
164 self.dataIn.nHeights),
164 self.dataIn.nHeights),
165 dtype='complex')
165 dtype='complex')
166
166
167 if self.dataIn.flagDataAsBlock:
167 if self.dataIn.flagDataAsBlock:
168 nVoltProfiles = self.dataIn.data.shape[1]
168 nVoltProfiles = self.dataIn.data.shape[1]
169
169
170 if nVoltProfiles == nProfiles:
170 if nVoltProfiles == nProfiles:
171 self.buffer = self.dataIn.data.copy()
171 self.buffer = self.dataIn.data.copy()
172 self.profIndex = nVoltProfiles
172 self.profIndex = nVoltProfiles
173
173
174 elif nVoltProfiles < nProfiles:
174 elif nVoltProfiles < nProfiles:
175
175
176 if self.profIndex == 0:
176 if self.profIndex == 0:
177 self.id_min = 0
177 self.id_min = 0
178 self.id_max = nVoltProfiles
178 self.id_max = nVoltProfiles
179
179
180 self.buffer[:, self.id_min:self.id_max,
180 self.buffer[:, self.id_min:self.id_max,
181 :] = self.dataIn.data
181 :] = self.dataIn.data
182 self.profIndex += nVoltProfiles
182 self.profIndex += nVoltProfiles
183 self.id_min += nVoltProfiles
183 self.id_min += nVoltProfiles
184 self.id_max += nVoltProfiles
184 self.id_max += nVoltProfiles
185 else:
185 else:
186 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
186 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
187 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
187 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
188 self.dataOut.flagNoData = True
188 self.dataOut.flagNoData = True
189 else:
189 else:
190 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
190 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
191 self.profIndex += 1
191 self.profIndex += 1
192
192
193 if self.firstdatatime == None:
193 if self.firstdatatime == None:
194 self.firstdatatime = self.dataIn.utctime
194 self.firstdatatime = self.dataIn.utctime
195
195
196 if self.profIndex == nProfiles:
196 if self.profIndex == nProfiles:
197 self.__updateSpecFromVoltage()
197 self.__updateSpecFromVoltage()
198 if pairsList == None:
198 if pairsList == None:
199 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
199 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
200 else:
200 else:
201 self.dataOut.pairsList = pairsList
201 self.dataOut.pairsList = pairsList
202 self.__getFft()
202 self.__getFft()
203 self.dataOut.flagNoData = False
203 self.dataOut.flagNoData = False
204 self.firstdatatime = None
204 self.firstdatatime = None
205 self.profIndex = 0
205 self.profIndex = 0
206 else:
206 else:
207 raise ValueError("The type of input object '%s' is not valid".format(
207 raise ValueError("The type of input object '%s' is not valid".format(
208 self.dataIn.type))
208 self.dataIn.type))
209
209
210 def __selectPairs(self, pairsList):
210 def __selectPairs(self, pairsList):
211
211
212 if not pairsList:
212 if not pairsList:
213 return
213 return
214
214
215 pairs = []
215 pairs = []
216 pairsIndex = []
216 pairsIndex = []
217
217
218 for pair in pairsList:
218 for pair in pairsList:
219 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
219 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
220 continue
220 continue
221 pairs.append(pair)
221 pairs.append(pair)
222 pairsIndex.append(pairs.index(pair))
222 pairsIndex.append(pairs.index(pair))
223
223
224 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
224 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
225 self.dataOut.pairsList = pairs
225 self.dataOut.pairsList = pairs
226
226
227 return
227 return
228
228
229 def selectFFTs(self, minFFT, maxFFT ):
229 def selectFFTs(self, minFFT, maxFFT ):
230 """
230 """
231 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
231 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
232 minFFT<= FFT <= maxFFT
232 minFFT<= FFT <= maxFFT
233 """
233 """
234
234
235 if (minFFT > maxFFT):
235 if (minFFT > maxFFT):
236 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
236 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
237
237
238 if (minFFT < self.dataOut.getFreqRange()[0]):
238 if (minFFT < self.dataOut.getFreqRange()[0]):
239 minFFT = self.dataOut.getFreqRange()[0]
239 minFFT = self.dataOut.getFreqRange()[0]
240
240
241 if (maxFFT > self.dataOut.getFreqRange()[-1]):
241 if (maxFFT > self.dataOut.getFreqRange()[-1]):
242 maxFFT = self.dataOut.getFreqRange()[-1]
242 maxFFT = self.dataOut.getFreqRange()[-1]
243
243
244 minIndex = 0
244 minIndex = 0
245 maxIndex = 0
245 maxIndex = 0
246 FFTs = self.dataOut.getFreqRange()
246 FFTs = self.dataOut.getFreqRange()
247
247
248 inda = numpy.where(FFTs >= minFFT)
248 inda = numpy.where(FFTs >= minFFT)
249 indb = numpy.where(FFTs <= maxFFT)
249 indb = numpy.where(FFTs <= maxFFT)
250
250
251 try:
251 try:
252 minIndex = inda[0][0]
252 minIndex = inda[0][0]
253 except:
253 except:
254 minIndex = 0
254 minIndex = 0
255
255
256 try:
256 try:
257 maxIndex = indb[0][-1]
257 maxIndex = indb[0][-1]
258 except:
258 except:
259 maxIndex = len(FFTs)
259 maxIndex = len(FFTs)
260
260
261 self.selectFFTsByIndex(minIndex, maxIndex)
261 self.selectFFTsByIndex(minIndex, maxIndex)
262
262
263 return 1
263 return 1
264
264
265 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
265 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
266 newheis = numpy.where(
266 newheis = numpy.where(
267 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
267 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
268
268
269 if hei_ref != None:
269 if hei_ref != None:
270 newheis = numpy.where(self.dataOut.heightList > hei_ref)
270 newheis = numpy.where(self.dataOut.heightList > hei_ref)
271
271
272 minIndex = min(newheis[0])
272 minIndex = min(newheis[0])
273 maxIndex = max(newheis[0])
273 maxIndex = max(newheis[0])
274 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
274 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
275 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
275 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
276
276
277 # determina indices
277 # determina indices
278 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
278 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
279 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
279 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
280 avg_dB = 10 * \
280 avg_dB = 10 * \
281 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
281 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
282 beacon_dB = numpy.sort(avg_dB)[-nheis:]
282 beacon_dB = numpy.sort(avg_dB)[-nheis:]
283 beacon_heiIndexList = []
283 beacon_heiIndexList = []
284 for val in avg_dB.tolist():
284 for val in avg_dB.tolist():
285 if val >= beacon_dB[0]:
285 if val >= beacon_dB[0]:
286 beacon_heiIndexList.append(avg_dB.tolist().index(val))
286 beacon_heiIndexList.append(avg_dB.tolist().index(val))
287
287
288 #data_spc = data_spc[:,:,beacon_heiIndexList]
288 #data_spc = data_spc[:,:,beacon_heiIndexList]
289 data_cspc = None
289 data_cspc = None
290 if self.dataOut.data_cspc is not None:
290 if self.dataOut.data_cspc is not None:
291 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
291 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
292 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
292 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
293
293
294 data_dc = None
294 data_dc = None
295 if self.dataOut.data_dc is not None:
295 if self.dataOut.data_dc is not None:
296 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
296 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
297 #data_dc = data_dc[:,beacon_heiIndexList]
297 #data_dc = data_dc[:,beacon_heiIndexList]
298
298
299 self.dataOut.data_spc = data_spc
299 self.dataOut.data_spc = data_spc
300 self.dataOut.data_cspc = data_cspc
300 self.dataOut.data_cspc = data_cspc
301 self.dataOut.data_dc = data_dc
301 self.dataOut.data_dc = data_dc
302 self.dataOut.heightList = heightList
302 self.dataOut.heightList = heightList
303 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
303 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
304
304
305 return 1
305 return 1
306
306
307 def selectFFTsByIndex(self, minIndex, maxIndex):
307 def selectFFTsByIndex(self, minIndex, maxIndex):
308 """
308 """
309
309
310 """
310 """
311
311
312 if (minIndex < 0) or (minIndex > maxIndex):
312 if (minIndex < 0) or (minIndex > maxIndex):
313 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
313 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
314
314
315 if (maxIndex >= self.dataOut.nProfiles):
315 if (maxIndex >= self.dataOut.nProfiles):
316 maxIndex = self.dataOut.nProfiles-1
316 maxIndex = self.dataOut.nProfiles-1
317
317
318 #Spectra
318 #Spectra
319 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
319 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
320
320
321 data_cspc = None
321 data_cspc = None
322 if self.dataOut.data_cspc is not None:
322 if self.dataOut.data_cspc is not None:
323 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
323 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
324
324
325 data_dc = None
325 data_dc = None
326 if self.dataOut.data_dc is not None:
326 if self.dataOut.data_dc is not None:
327 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
327 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
328
328
329 self.dataOut.data_spc = data_spc
329 self.dataOut.data_spc = data_spc
330 self.dataOut.data_cspc = data_cspc
330 self.dataOut.data_cspc = data_cspc
331 self.dataOut.data_dc = data_dc
331 self.dataOut.data_dc = data_dc
332
332
333 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
333 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
334 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
334 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
335 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
335 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
336
336
337 return 1
337 return 1
338
338
339 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
339 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
340 # validacion de rango
340 # validacion de rango
341 if minHei == None:
341 if minHei == None:
342 minHei = self.dataOut.heightList[0]
342 minHei = self.dataOut.heightList[0]
343
343
344 if maxHei == None:
344 if maxHei == None:
345 maxHei = self.dataOut.heightList[-1]
345 maxHei = self.dataOut.heightList[-1]
346
346
347 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
347 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
348 print('minHei: %.2f is out of the heights range' % (minHei))
348 print('minHei: %.2f is out of the heights range' % (minHei))
349 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
349 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
350 minHei = self.dataOut.heightList[0]
350 minHei = self.dataOut.heightList[0]
351
351
352 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
352 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
353 print('maxHei: %.2f is out of the heights range' % (maxHei))
353 print('maxHei: %.2f is out of the heights range' % (maxHei))
354 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
354 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
355 maxHei = self.dataOut.heightList[-1]
355 maxHei = self.dataOut.heightList[-1]
356
356
357 # validacion de velocidades
357 # validacion de velocidades
358 velrange = self.dataOut.getVelRange(1)
358 velrange = self.dataOut.getVelRange(1)
359
359
360 if minVel == None:
360 if minVel == None:
361 minVel = velrange[0]
361 minVel = velrange[0]
362
362
363 if maxVel == None:
363 if maxVel == None:
364 maxVel = velrange[-1]
364 maxVel = velrange[-1]
365
365
366 if (minVel < velrange[0]) or (minVel > maxVel):
366 if (minVel < velrange[0]) or (minVel > maxVel):
367 print('minVel: %.2f is out of the velocity range' % (minVel))
367 print('minVel: %.2f is out of the velocity range' % (minVel))
368 print('minVel is setting to %.2f' % (velrange[0]))
368 print('minVel is setting to %.2f' % (velrange[0]))
369 minVel = velrange[0]
369 minVel = velrange[0]
370
370
371 if (maxVel > velrange[-1]) or (maxVel < minVel):
371 if (maxVel > velrange[-1]) or (maxVel < minVel):
372 print('maxVel: %.2f is out of the velocity range' % (maxVel))
372 print('maxVel: %.2f is out of the velocity range' % (maxVel))
373 print('maxVel is setting to %.2f' % (velrange[-1]))
373 print('maxVel is setting to %.2f' % (velrange[-1]))
374 maxVel = velrange[-1]
374 maxVel = velrange[-1]
375
375
376 # seleccion de indices para rango
376 # seleccion de indices para rango
377 minIndex = 0
377 minIndex = 0
378 maxIndex = 0
378 maxIndex = 0
379 heights = self.dataOut.heightList
379 heights = self.dataOut.heightList
380
380
381 inda = numpy.where(heights >= minHei)
381 inda = numpy.where(heights >= minHei)
382 indb = numpy.where(heights <= maxHei)
382 indb = numpy.where(heights <= maxHei)
383
383
384 try:
384 try:
385 minIndex = inda[0][0]
385 minIndex = inda[0][0]
386 except:
386 except:
387 minIndex = 0
387 minIndex = 0
388
388
389 try:
389 try:
390 maxIndex = indb[0][-1]
390 maxIndex = indb[0][-1]
391 except:
391 except:
392 maxIndex = len(heights)
392 maxIndex = len(heights)
393
393
394 if (minIndex < 0) or (minIndex > maxIndex):
394 if (minIndex < 0) or (minIndex > maxIndex):
395 raise ValueError("some value in (%d,%d) is not valid" % (
395 raise ValueError("some value in (%d,%d) is not valid" % (
396 minIndex, maxIndex))
396 minIndex, maxIndex))
397
397
398 if (maxIndex >= self.dataOut.nHeights):
398 if (maxIndex >= self.dataOut.nHeights):
399 maxIndex = self.dataOut.nHeights - 1
399 maxIndex = self.dataOut.nHeights - 1
400
400
401 # seleccion de indices para velocidades
401 # seleccion de indices para velocidades
402 indminvel = numpy.where(velrange >= minVel)
402 indminvel = numpy.where(velrange >= minVel)
403 indmaxvel = numpy.where(velrange <= maxVel)
403 indmaxvel = numpy.where(velrange <= maxVel)
404 try:
404 try:
405 minIndexVel = indminvel[0][0]
405 minIndexVel = indminvel[0][0]
406 except:
406 except:
407 minIndexVel = 0
407 minIndexVel = 0
408
408
409 try:
409 try:
410 maxIndexVel = indmaxvel[0][-1]
410 maxIndexVel = indmaxvel[0][-1]
411 except:
411 except:
412 maxIndexVel = len(velrange)
412 maxIndexVel = len(velrange)
413
413
414 # seleccion del espectro
414 # seleccion del espectro
415 data_spc = self.dataOut.data_spc[:,
415 data_spc = self.dataOut.data_spc[:,
416 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
416 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
417 # estimacion de ruido
417 # estimacion de ruido
418 noise = numpy.zeros(self.dataOut.nChannels)
418 noise = numpy.zeros(self.dataOut.nChannels)
419
419
420 for channel in range(self.dataOut.nChannels):
420 for channel in range(self.dataOut.nChannels):
421 daux = data_spc[channel, :, :]
421 daux = data_spc[channel, :, :]
422 sortdata = numpy.sort(daux, axis=None)
422 sortdata = numpy.sort(daux, axis=None)
423 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
423 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
424
424
425 self.dataOut.noise_estimation = noise.copy()
425 self.dataOut.noise_estimation = noise.copy()
426
426
427 return 1
427 return 1
428
428
429 class removeDC(Operation):
429 class removeDC(Operation):
430
430
431 def run(self, dataOut, mode=2):
431 def run(self, dataOut, mode=2):
432 self.dataOut = dataOut
432 self.dataOut = dataOut
433 jspectra = self.dataOut.data_spc
433 jspectra = self.dataOut.data_spc
434 jcspectra = self.dataOut.data_cspc
434 jcspectra = self.dataOut.data_cspc
435
435
436 num_chan = jspectra.shape[0]
436 num_chan = jspectra.shape[0]
437 num_hei = jspectra.shape[2]
437 num_hei = jspectra.shape[2]
438
438
439 if jcspectra is not None:
439 if jcspectra is not None:
440 jcspectraExist = True
440 jcspectraExist = True
441 num_pairs = jcspectra.shape[0]
441 num_pairs = jcspectra.shape[0]
442 else:
442 else:
443 jcspectraExist = False
443 jcspectraExist = False
444
444
445 freq_dc = int(jspectra.shape[1] / 2)
445 freq_dc = int(jspectra.shape[1] / 2)
446 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
446 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
447 ind_vel = ind_vel.astype(int)
447 ind_vel = ind_vel.astype(int)
448
448
449 if ind_vel[0] < 0:
449 if ind_vel[0] < 0:
450 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
450 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
451
451
452 if mode == 1:
452 if mode == 1:
453 jspectra[:, freq_dc, :] = (
453 jspectra[:, freq_dc, :] = (
454 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
454 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
455
455
456 if jcspectraExist:
456 if jcspectraExist:
457 jcspectra[:, freq_dc, :] = (
457 jcspectra[:, freq_dc, :] = (
458 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
458 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
459
459
460 if mode == 2:
460 if mode == 2:
461
461
462 vel = numpy.array([-2, -1, 1, 2])
462 vel = numpy.array([-2, -1, 1, 2])
463 xx = numpy.zeros([4, 4])
463 xx = numpy.zeros([4, 4])
464
464
465 for fil in range(4):
465 for fil in range(4):
466 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
466 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
467
467
468 xx_inv = numpy.linalg.inv(xx)
468 xx_inv = numpy.linalg.inv(xx)
469 xx_aux = xx_inv[0, :]
469 xx_aux = xx_inv[0, :]
470
470
471 for ich in range(num_chan):
471 for ich in range(num_chan):
472 yy = jspectra[ich, ind_vel, :]
472 yy = jspectra[ich, ind_vel, :]
473 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
473 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
474
474
475 junkid = jspectra[ich, freq_dc, :] <= 0
475 junkid = jspectra[ich, freq_dc, :] <= 0
476 cjunkid = sum(junkid)
476 cjunkid = sum(junkid)
477
477
478 if cjunkid.any():
478 if cjunkid.any():
479 jspectra[ich, freq_dc, junkid.nonzero()] = (
479 jspectra[ich, freq_dc, junkid.nonzero()] = (
480 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
480 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
481
481
482 if jcspectraExist:
482 if jcspectraExist:
483 for ip in range(num_pairs):
483 for ip in range(num_pairs):
484 yy = jcspectra[ip, ind_vel, :]
484 yy = jcspectra[ip, ind_vel, :]
485 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
485 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
486
486
487 self.dataOut.data_spc = jspectra
487 self.dataOut.data_spc = jspectra
488 self.dataOut.data_cspc = jcspectra
488 self.dataOut.data_cspc = jcspectra
489
489
490 return self.dataOut
490 return self.dataOut
491
491
492 # import matplotlib.pyplot as plt
492 # import matplotlib.pyplot as plt
493
493
494 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
494 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
495 z = (x - a1) / a2
495 z = (x - a1) / a2
496 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
496 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
497 return y
497 return y
498
498
499
499
500 class CleanRayleigh(Operation):
500 class CleanRayleigh(Operation):
501
501
502 def __init__(self):
502 def __init__(self):
503
503
504 Operation.__init__(self)
504 Operation.__init__(self)
505 self.i=0
505 self.i=0
506 self.isConfig = False
506 self.isConfig = False
507 self.__dataReady = False
507 self.__dataReady = False
508 self.__profIndex = 0
508 self.__profIndex = 0
509 self.byTime = False
509 self.byTime = False
510 self.byProfiles = False
510 self.byProfiles = False
511
511
512 self.bloques = None
512 self.bloques = None
513 self.bloque0 = None
513 self.bloque0 = None
514
514
515 self.index = 0
515 self.index = 0
516
516
517 self.buffer = 0
517 self.buffer = 0
518 self.buffer2 = 0
518 self.buffer2 = 0
519 self.buffer3 = 0
519 self.buffer3 = 0
520
520
521
521
522 def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
522 def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
523
523
524 self.nChannels = dataOut.nChannels
524 self.nChannels = dataOut.nChannels
525 self.nProf = dataOut.nProfiles
525 self.nProf = dataOut.nProfiles
526 self.nPairs = dataOut.data_cspc.shape[0]
526 self.nPairs = dataOut.data_cspc.shape[0]
527 self.pairsArray = numpy.array(dataOut.pairsList)
527 self.pairsArray = numpy.array(dataOut.pairsList)
528 self.spectra = dataOut.data_spc
528 self.spectra = dataOut.data_spc
529 self.cspectra = dataOut.data_cspc
529 self.cspectra = dataOut.data_cspc
530 self.heights = dataOut.heightList #alturas totales
530 self.heights = dataOut.heightList #alturas totales
531 self.nHeights = len(self.heights)
531 self.nHeights = len(self.heights)
532 self.min_hei = min_hei
532 self.min_hei = min_hei
533 self.max_hei = max_hei
533 self.max_hei = max_hei
534 if (self.min_hei == None):
534 if (self.min_hei == None):
535 self.min_hei = 0
535 self.min_hei = 0
536 if (self.max_hei == None):
536 if (self.max_hei == None):
537 self.max_hei = dataOut.heightList[-1]
537 self.max_hei = dataOut.heightList[-1]
538 self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
538 self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
539 self.heightsClean = self.heights[self.hval] #alturas filtradas
539 self.heightsClean = self.heights[self.hval] #alturas filtradas
540 self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
540 self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
541 self.nHeightsClean = len(self.heightsClean)
541 self.nHeightsClean = len(self.heightsClean)
542 self.channels = dataOut.channelList
542 self.channels = dataOut.channelList
543 self.nChan = len(self.channels)
543 self.nChan = len(self.channels)
544 self.nIncohInt = dataOut.nIncohInt
544 self.nIncohInt = dataOut.nIncohInt
545 self.__initime = dataOut.utctime
545 self.__initime = dataOut.utctime
546 self.maxAltInd = self.hval[-1]+1
546 self.maxAltInd = self.hval[-1]+1
547 self.minAltInd = self.hval[0]
547 self.minAltInd = self.hval[0]
548
548
549 self.crosspairs = dataOut.pairsList
549 self.crosspairs = dataOut.pairsList
550 self.nPairs = len(self.crosspairs)
550 self.nPairs = len(self.crosspairs)
551 self.normFactor = dataOut.normFactor
551 self.normFactor = dataOut.normFactor
552 self.nFFTPoints = dataOut.nFFTPoints
552 self.nFFTPoints = dataOut.nFFTPoints
553 self.ippSeconds = dataOut.ippSeconds
553 self.ippSeconds = dataOut.ippSeconds
554 self.currentTime = self.__initime
554 self.currentTime = self.__initime
555 self.pairsArray = numpy.array(dataOut.pairsList)
555 self.pairsArray = numpy.array(dataOut.pairsList)
556 self.factor_stdv = factor_stdv
556 self.factor_stdv = factor_stdv
557
557
558 if n != None :
558 if n != None :
559 self.byProfiles = True
559 self.byProfiles = True
560 self.nIntProfiles = n
560 self.nIntProfiles = n
561 else:
561 else:
562 self.__integrationtime = timeInterval
562 self.__integrationtime = timeInterval
563
563
564 self.__dataReady = False
564 self.__dataReady = False
565 self.isConfig = True
565 self.isConfig = True
566
566
567
567
568
568
569 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
569 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
570
570
571 if not self.isConfig :
571 if not self.isConfig :
572
572
573 self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
573 self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
574
574
575 tini=dataOut.utctime
575 tini=dataOut.utctime
576
576
577 if self.byProfiles:
577 if self.byProfiles:
578 if self.__profIndex == self.nIntProfiles:
578 if self.__profIndex == self.nIntProfiles:
579 self.__dataReady = True
579 self.__dataReady = True
580 else:
580 else:
581 if (tini - self.__initime) >= self.__integrationtime:
581 if (tini - self.__initime) >= self.__integrationtime:
582
582
583 self.__dataReady = True
583 self.__dataReady = True
584 self.__initime = tini
584 self.__initime = tini
585
585
586 #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
586 #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
587
587
588 if self.__dataReady:
588 if self.__dataReady:
589
589
590 self.__profIndex = 0
590 self.__profIndex = 0
591 jspc = self.buffer
591 jspc = self.buffer
592 jcspc = self.buffer2
592 jcspc = self.buffer2
593 #jnoise = self.buffer3
593 #jnoise = self.buffer3
594 self.buffer = dataOut.data_spc
594 self.buffer = dataOut.data_spc
595 self.buffer2 = dataOut.data_cspc
595 self.buffer2 = dataOut.data_cspc
596 #self.buffer3 = dataOut.noise
596 #self.buffer3 = dataOut.noise
597 self.currentTime = dataOut.utctime
597 self.currentTime = dataOut.utctime
598 if numpy.any(jspc) :
598 if numpy.any(jspc) :
599 #print( jspc.shape, jcspc.shape)
599 #print( jspc.shape, jcspc.shape)
600 jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
600 jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
601 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
601 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
602 self.__dataReady = False
602 self.__dataReady = False
603 #print( jspc.shape, jcspc.shape)
603 #print( jspc.shape, jcspc.shape)
604 dataOut.flagNoData = False
604 dataOut.flagNoData = False
605 else:
605 else:
606 dataOut.flagNoData = True
606 dataOut.flagNoData = True
607 self.__dataReady = False
607 self.__dataReady = False
608 return dataOut
608 return dataOut
609 else:
609 else:
610 #print( len(self.buffer))
610 #print( len(self.buffer))
611 if numpy.any(self.buffer):
611 if numpy.any(self.buffer):
612 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
612 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
613 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
613 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
614 self.buffer3 += dataOut.data_dc
614 self.buffer3 += dataOut.data_dc
615 else:
615 else:
616 self.buffer = dataOut.data_spc
616 self.buffer = dataOut.data_spc
617 self.buffer2 = dataOut.data_cspc
617 self.buffer2 = dataOut.data_cspc
618 self.buffer3 = dataOut.data_dc
618 self.buffer3 = dataOut.data_dc
619 #print self.index, self.fint
619 #print self.index, self.fint
620 #print self.buffer2.shape
620 #print self.buffer2.shape
621 dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
621 dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
622 self.__profIndex += 1
622 self.__profIndex += 1
623 return dataOut ## NOTE: REV
623 return dataOut ## NOTE: REV
624
624
625
625
626 #index = tini.tm_hour*12+tini.tm_min/5
626 #index = tini.tm_hour*12+tini.tm_min/5
627 '''REVISAR'''
627 '''REVISAR'''
628 # jspc = jspc/self.nFFTPoints/self.normFactor
628 # jspc = jspc/self.nFFTPoints/self.normFactor
629 # jcspc = jcspc/self.nFFTPoints/self.normFactor
629 # jcspc = jcspc/self.nFFTPoints/self.normFactor
630
630
631
631
632
632
633 tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
633 tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
634 dataOut.data_spc = tmp_spectra
634 dataOut.data_spc = tmp_spectra
635 dataOut.data_cspc = tmp_cspectra
635 dataOut.data_cspc = tmp_cspectra
636
636
637 #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
637 #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
638
638
639 dataOut.data_dc = self.buffer3
639 dataOut.data_dc = self.buffer3
640 dataOut.nIncohInt *= self.nIntProfiles
640 dataOut.nIncohInt *= self.nIntProfiles
641 dataOut.utctime = self.currentTime #tiempo promediado
641 dataOut.utctime = self.currentTime #tiempo promediado
642 #print("Time: ",time.localtime(dataOut.utctime))
642 #print("Time: ",time.localtime(dataOut.utctime))
643 # dataOut.data_spc = sat_spectra
643 # dataOut.data_spc = sat_spectra
644 # dataOut.data_cspc = sat_cspectra
644 # dataOut.data_cspc = sat_cspectra
645 self.buffer = 0
645 self.buffer = 0
646 self.buffer2 = 0
646 self.buffer2 = 0
647 self.buffer3 = 0
647 self.buffer3 = 0
648
648
649 return dataOut
649 return dataOut
650
650
651 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
651 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
652 #print("OP cleanRayleigh")
652 #print("OP cleanRayleigh")
653 #import matplotlib.pyplot as plt
653 #import matplotlib.pyplot as plt
654 #for k in range(149):
654 #for k in range(149):
655 #channelsProcssd = []
655 #channelsProcssd = []
656 #channelA_ok = False
656 #channelA_ok = False
657 #rfunc = cspectra.copy() #self.bloques
657 #rfunc = cspectra.copy() #self.bloques
658 rfunc = spectra.copy()
658 rfunc = spectra.copy()
659 #rfunc = cspectra
659 #rfunc = cspectra
660 #val_spc = spectra*0.0 #self.bloque0*0.0
660 #val_spc = spectra*0.0 #self.bloque0*0.0
661 #val_cspc = cspectra*0.0 #self.bloques*0.0
661 #val_cspc = cspectra*0.0 #self.bloques*0.0
662 #in_sat_spectra = spectra.copy() #self.bloque0
662 #in_sat_spectra = spectra.copy() #self.bloque0
663 #in_sat_cspectra = cspectra.copy() #self.bloques
663 #in_sat_cspectra = cspectra.copy() #self.bloques
664
664
665
665
666 ###ONLY FOR TEST:
666 ###ONLY FOR TEST:
667 raxs = math.ceil(math.sqrt(self.nPairs))
667 raxs = math.ceil(math.sqrt(self.nPairs))
668 caxs = math.ceil(self.nPairs/raxs)
668 caxs = math.ceil(self.nPairs/raxs)
669 if self.nPairs <4:
669 if self.nPairs <4:
670 raxs = 2
670 raxs = 2
671 caxs = 2
671 caxs = 2
672 #print(raxs, caxs)
672 #print(raxs, caxs)
673 fft_rev = 14 #nFFT to plot
673 fft_rev = 14 #nFFT to plot
674 hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot
674 hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot
675 hei_rev = hei_rev[0]
675 hei_rev = hei_rev[0]
676 #print(hei_rev)
676 #print(hei_rev)
677
677
678 #print numpy.absolute(rfunc[:,0,0,14])
678 #print numpy.absolute(rfunc[:,0,0,14])
679
679
680 gauss_fit, covariance = None, None
680 gauss_fit, covariance = None, None
681 for ih in range(self.minAltInd,self.maxAltInd):
681 for ih in range(self.minAltInd,self.maxAltInd):
682 for ifreq in range(self.nFFTPoints):
682 for ifreq in range(self.nFFTPoints):
683 '''
683 '''
684 ###ONLY FOR TEST:
684 ###ONLY FOR TEST:
685 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
685 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
686 fig, axs = plt.subplots(raxs, caxs)
686 fig, axs = plt.subplots(raxs, caxs)
687 fig2, axs2 = plt.subplots(raxs, caxs)
687 fig2, axs2 = plt.subplots(raxs, caxs)
688 col_ax = 0
688 col_ax = 0
689 row_ax = 0
689 row_ax = 0
690 '''
690 '''
691 #print(self.nPairs)
691 #print(self.nPairs)
692 for ii in range(self.nChan): #PARES DE CANALES SELF y CROSS
692 for ii in range(self.nChan): #PARES DE CANALES SELF y CROSS
693 # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
693 # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
694 # continue
694 # continue
695 # if not self.crosspairs[ii][0] in channelsProcssd:
695 # if not self.crosspairs[ii][0] in channelsProcssd:
696 # channelA_ok = True
696 # channelA_ok = True
697 #print("pair: ",self.crosspairs[ii])
697 #print("pair: ",self.crosspairs[ii])
698 '''
698 '''
699 ###ONLY FOR TEST:
699 ###ONLY FOR TEST:
700 if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1):
700 if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1):
701 col_ax = 0
701 col_ax = 0
702 row_ax += 1
702 row_ax += 1
703 '''
703 '''
704 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
704 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
705 #print(func2clean.shape)
705 #print(func2clean.shape)
706 val = (numpy.isfinite(func2clean)==True).nonzero()
706 val = (numpy.isfinite(func2clean)==True).nonzero()
707
707
708 if len(val)>0: #limitador
708 if len(val)>0: #limitador
709 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
709 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
710 if min_val <= -40 :
710 if min_val <= -40 :
711 min_val = -40
711 min_val = -40
712 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
712 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
713 if max_val >= 200 :
713 if max_val >= 200 :
714 max_val = 200
714 max_val = 200
715 #print min_val, max_val
715 #print min_val, max_val
716 step = 1
716 step = 1
717 #print("Getting bins and the histogram")
717 #print("Getting bins and the histogram")
718 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
718 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
719 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
719 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
720 #print(len(y_dist),len(binstep[:-1]))
720 #print(len(y_dist),len(binstep[:-1]))
721 #print(row_ax,col_ax, " ..")
721 #print(row_ax,col_ax, " ..")
722 #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
722 #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
723 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
723 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
724 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
724 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
725 parg = [numpy.amax(y_dist),mean,sigma]
725 parg = [numpy.amax(y_dist),mean,sigma]
726
726
727 newY = None
727 newY = None
728
728
729 try :
729 try :
730 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
730 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
731 mode = gauss_fit[1]
731 mode = gauss_fit[1]
732 stdv = gauss_fit[2]
732 stdv = gauss_fit[2]
733 #print(" FIT OK",gauss_fit)
733 #print(" FIT OK",gauss_fit)
734 '''
734 '''
735 ###ONLY FOR TEST:
735 ###ONLY FOR TEST:
736 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
736 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
737 newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
737 newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
738 axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
738 axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
739 axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
739 axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
740 axs[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
740 axs[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
741 '''
741 '''
742 except:
742 except:
743 mode = mean
743 mode = mean
744 stdv = sigma
744 stdv = sigma
745 #print("FIT FAIL")
745 #print("FIT FAIL")
746 #continue
746 #continue
747
747
748
748
749 #print(mode,stdv)
749 #print(mode,stdv)
750 #Removing echoes greater than mode + std_factor*stdv
750 #Removing echoes greater than mode + std_factor*stdv
751 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
751 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
752 #noval tiene los indices que se van a remover
752 #noval tiene los indices que se van a remover
753 #print("Chan ",ii," novals: ",len(noval[0]))
753 #print("Chan ",ii," novals: ",len(noval[0]))
754 if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
754 if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
755 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
755 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
756 #print(novall)
756 #print(novall)
757 #print(" ",self.pairsArray[ii])
757 #print(" ",self.pairsArray[ii])
758 #cross_pairs = self.pairsArray[ii]
758 #cross_pairs = self.pairsArray[ii]
759 #Getting coherent echoes which are removed.
759 #Getting coherent echoes which are removed.
760 # if len(novall[0]) > 0:
760 # if len(novall[0]) > 0:
761 #
761 #
762 # val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
762 # val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
763 # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
763 # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
764 # val_cspc[novall[0],ii,ifreq,ih] = 1
764 # val_cspc[novall[0],ii,ifreq,ih] = 1
765 #print("OUT NOVALL 1")
765 #print("OUT NOVALL 1")
766 try:
766 try:
767 pair = (self.channels[ii],self.channels[ii + 1])
767 pair = (self.channels[ii],self.channels[ii + 1])
768 except:
768 except:
769 pair = (99,99)
769 pair = (99,99)
770 #print("par ", pair)
770 #print("par ", pair)
771 if ( pair in self.crosspairs):
771 if ( pair in self.crosspairs):
772 q = self.crosspairs.index(pair)
772 q = self.crosspairs.index(pair)
773 #print("estΓ‘ aqui: ", q, (ii,ii + 1))
773 #print("estΓ‘ aqui: ", q, (ii,ii + 1))
774 new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0])
774 new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0])
775 cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
775 cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
776
776
777 #if channelA_ok:
777 #if channelA_ok:
778 #chA = self.channels.index(cross_pairs[0])
778 #chA = self.channels.index(cross_pairs[0])
779 new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0])
779 new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0])
780 spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
780 spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
781 #channelA_ok = False
781 #channelA_ok = False
782
782
783 # chB = self.channels.index(cross_pairs[1])
783 # chB = self.channels.index(cross_pairs[1])
784 # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
784 # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
785 # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
785 # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
786 #
786 #
787 # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A
787 # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A
788 # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B
788 # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B
789 '''
789 '''
790 ###ONLY FOR TEST:
790 ###ONLY FOR TEST:
791 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
791 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
792 func2clean = 10*numpy.log10(numpy.absolute(spectra[:,ii,ifreq,ih]))
792 func2clean = 10*numpy.log10(numpy.absolute(spectra[:,ii,ifreq,ih]))
793 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
793 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
794 axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
794 axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
795 axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
795 axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
796 axs2[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
796 axs2[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
797 '''
797 '''
798 '''
798 '''
799 ###ONLY FOR TEST:
799 ###ONLY FOR TEST:
800 col_ax += 1 #contador de ploteo columnas
800 col_ax += 1 #contador de ploteo columnas
801 ##print(col_ax)
801 ##print(col_ax)
802 ###ONLY FOR TEST:
802 ###ONLY FOR TEST:
803 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
803 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
804 title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
804 title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
805 title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
805 title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
806 fig.suptitle(title)
806 fig.suptitle(title)
807 fig2.suptitle(title2)
807 fig2.suptitle(title2)
808 plt.show()
808 plt.show()
809 '''
809 '''
810 ##################################################################################################
810 ##################################################################################################
811
811
812 #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
812 #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
813 out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
813 out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
814 out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
814 out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
815 for ih in range(self.nHeights):
815 for ih in range(self.nHeights):
816 for ifreq in range(self.nFFTPoints):
816 for ifreq in range(self.nFFTPoints):
817 for ich in range(self.nChan):
817 for ich in range(self.nChan):
818 tmp = spectra[:,ich,ifreq,ih]
818 tmp = spectra[:,ich,ifreq,ih]
819 valid = (numpy.isfinite(tmp[:])==True).nonzero()
819 valid = (numpy.isfinite(tmp[:])==True).nonzero()
820
820
821 if len(valid[0]) >0 :
821 if len(valid[0]) >0 :
822 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
822 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
823
823
824 for icr in range(self.nPairs):
824 for icr in range(self.nPairs):
825 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
825 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
826 valid = (numpy.isfinite(tmp)==True).nonzero()
826 valid = (numpy.isfinite(tmp)==True).nonzero()
827 if len(valid[0]) > 0:
827 if len(valid[0]) > 0:
828 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
828 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
829
829
830 return out_spectra, out_cspectra
830 return out_spectra, out_cspectra
831
831
832 def REM_ISOLATED_POINTS(self,array,rth):
832 def REM_ISOLATED_POINTS(self,array,rth):
833 # import matplotlib.pyplot as plt
833 # import matplotlib.pyplot as plt
834 if rth == None :
834 if rth == None :
835 rth = 4
835 rth = 4
836 #print("REM ISO")
836 #print("REM ISO")
837 num_prof = len(array[0,:,0])
837 num_prof = len(array[0,:,0])
838 num_hei = len(array[0,0,:])
838 num_hei = len(array[0,0,:])
839 n2d = len(array[:,0,0])
839 n2d = len(array[:,0,0])
840
840
841 for ii in range(n2d) :
841 for ii in range(n2d) :
842 #print ii,n2d
842 #print ii,n2d
843 tmp = array[ii,:,:]
843 tmp = array[ii,:,:]
844 #print tmp.shape, array[ii,101,:],array[ii,102,:]
844 #print tmp.shape, array[ii,101,:],array[ii,102,:]
845
845
846 # fig = plt.figure(figsize=(6,5))
846 # fig = plt.figure(figsize=(6,5))
847 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
847 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
848 # ax = fig.add_axes([left, bottom, width, height])
848 # ax = fig.add_axes([left, bottom, width, height])
849 # x = range(num_prof)
849 # x = range(num_prof)
850 # y = range(num_hei)
850 # y = range(num_hei)
851 # cp = ax.contour(y,x,tmp)
851 # cp = ax.contour(y,x,tmp)
852 # ax.clabel(cp, inline=True,fontsize=10)
852 # ax.clabel(cp, inline=True,fontsize=10)
853 # plt.show()
853 # plt.show()
854
854
855 #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
855 #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
856 tmp = numpy.reshape(tmp,num_prof*num_hei)
856 tmp = numpy.reshape(tmp,num_prof*num_hei)
857 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
857 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
858 indxs2 = (tmp > 0).nonzero()
858 indxs2 = (tmp > 0).nonzero()
859
859
860 indxs1 = (indxs1[0])
860 indxs1 = (indxs1[0])
861 indxs2 = indxs2[0]
861 indxs2 = indxs2[0]
862 #indxs1 = numpy.array(indxs1[0])
862 #indxs1 = numpy.array(indxs1[0])
863 #indxs2 = numpy.array(indxs2[0])
863 #indxs2 = numpy.array(indxs2[0])
864 indxs = None
864 indxs = None
865 #print indxs1 , indxs2
865 #print indxs1 , indxs2
866 for iv in range(len(indxs2)):
866 for iv in range(len(indxs2)):
867 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
867 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
868 #print len(indxs2), indv
868 #print len(indxs2), indv
869 if len(indv[0]) > 0 :
869 if len(indv[0]) > 0 :
870 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
870 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
871 # print indxs
871 # print indxs
872 indxs = indxs[1:]
872 indxs = indxs[1:]
873 #print(indxs, len(indxs))
873 #print(indxs, len(indxs))
874 if len(indxs) < 4 :
874 if len(indxs) < 4 :
875 array[ii,:,:] = 0.
875 array[ii,:,:] = 0.
876 return
876 return
877
877
878 xpos = numpy.mod(indxs ,num_hei)
878 xpos = numpy.mod(indxs ,num_hei)
879 ypos = (indxs / num_hei)
879 ypos = (indxs / num_hei)
880 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
880 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
881 #print sx
881 #print sx
882 xpos = xpos[sx]
882 xpos = xpos[sx]
883 ypos = ypos[sx]
883 ypos = ypos[sx]
884
884
885 # *********************************** Cleaning isolated points **********************************
885 # *********************************** Cleaning isolated points **********************************
886 ic = 0
886 ic = 0
887 while True :
887 while True :
888 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
888 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
889 #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
889 #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
890 #plt.plot(r)
890 #plt.plot(r)
891 #plt.show()
891 #plt.show()
892 no_coh1 = (numpy.isfinite(r)==True).nonzero()
892 no_coh1 = (numpy.isfinite(r)==True).nonzero()
893 no_coh2 = (r <= rth).nonzero()
893 no_coh2 = (r <= rth).nonzero()
894 #print r, no_coh1, no_coh2
894 #print r, no_coh1, no_coh2
895 no_coh1 = numpy.array(no_coh1[0])
895 no_coh1 = numpy.array(no_coh1[0])
896 no_coh2 = numpy.array(no_coh2[0])
896 no_coh2 = numpy.array(no_coh2[0])
897 no_coh = None
897 no_coh = None
898 #print valid1 , valid2
898 #print valid1 , valid2
899 for iv in range(len(no_coh2)):
899 for iv in range(len(no_coh2)):
900 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
900 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
901 if len(indv[0]) > 0 :
901 if len(indv[0]) > 0 :
902 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
902 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
903 no_coh = no_coh[1:]
903 no_coh = no_coh[1:]
904 #print len(no_coh), no_coh
904 #print len(no_coh), no_coh
905 if len(no_coh) < 4 :
905 if len(no_coh) < 4 :
906 #print xpos[ic], ypos[ic], ic
906 #print xpos[ic], ypos[ic], ic
907 # plt.plot(r)
907 # plt.plot(r)
908 # plt.show()
908 # plt.show()
909 xpos[ic] = numpy.nan
909 xpos[ic] = numpy.nan
910 ypos[ic] = numpy.nan
910 ypos[ic] = numpy.nan
911
911
912 ic = ic + 1
912 ic = ic + 1
913 if (ic == len(indxs)) :
913 if (ic == len(indxs)) :
914 break
914 break
915 #print( xpos, ypos)
915 #print( xpos, ypos)
916
916
917 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
917 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
918 #print indxs[0]
918 #print indxs[0]
919 if len(indxs[0]) < 4 :
919 if len(indxs[0]) < 4 :
920 array[ii,:,:] = 0.
920 array[ii,:,:] = 0.
921 return
921 return
922
922
923 xpos = xpos[indxs[0]]
923 xpos = xpos[indxs[0]]
924 ypos = ypos[indxs[0]]
924 ypos = ypos[indxs[0]]
925 for i in range(0,len(ypos)):
925 for i in range(0,len(ypos)):
926 ypos[i]=int(ypos[i])
926 ypos[i]=int(ypos[i])
927 junk = tmp
927 junk = tmp
928 tmp = junk*0.0
928 tmp = junk*0.0
929
929
930 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
930 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
931 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
931 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
932
932
933 #print array.shape
933 #print array.shape
934 #tmp = numpy.reshape(tmp,(num_prof,num_hei))
934 #tmp = numpy.reshape(tmp,(num_prof,num_hei))
935 #print tmp.shape
935 #print tmp.shape
936
936
937 # fig = plt.figure(figsize=(6,5))
937 # fig = plt.figure(figsize=(6,5))
938 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
938 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
939 # ax = fig.add_axes([left, bottom, width, height])
939 # ax = fig.add_axes([left, bottom, width, height])
940 # x = range(num_prof)
940 # x = range(num_prof)
941 # y = range(num_hei)
941 # y = range(num_hei)
942 # cp = ax.contour(y,x,array[ii,:,:])
942 # cp = ax.contour(y,x,array[ii,:,:])
943 # ax.clabel(cp, inline=True,fontsize=10)
943 # ax.clabel(cp, inline=True,fontsize=10)
944 # plt.show()
944 # plt.show()
945 return array
945 return array
946
946
947
947
948 class IntegrationFaradaySpectra(Operation):
948 class IntegrationFaradaySpectra(Operation):
949
949
950 __profIndex = 0
950 __profIndex = 0
951 __withOverapping = False
951 __withOverapping = False
952
952
953 __byTime = False
953 __byTime = False
954 __initime = None
954 __initime = None
955 __lastdatatime = None
955 __lastdatatime = None
956 __integrationtime = None
956 __integrationtime = None
957
957
958 __buffer_spc = None
958 __buffer_spc = None
959 __buffer_cspc = None
959 __buffer_cspc = None
960 __buffer_dc = None
960 __buffer_dc = None
961
961
962 __dataReady = False
962 __dataReady = False
963
963
964 __timeInterval = None
964 __timeInterval = None
965
965
966 n = None
966 n = None
967
967
968 def __init__(self):
968 def __init__(self):
969
969
970 Operation.__init__(self)
970 Operation.__init__(self)
971
971
972 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None):
972 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None):
973 """
973 """
974 Set the parameters of the integration class.
974 Set the parameters of the integration class.
975
975
976 Inputs:
976 Inputs:
977
977
978 n : Number of coherent integrations
978 n : Number of coherent integrations
979 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
979 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
980 overlapping :
980 overlapping :
981
981
982 """
982 """
983
983
984 self.__initime = None
984 self.__initime = None
985 self.__lastdatatime = 0
985 self.__lastdatatime = 0
986
986
987 self.__buffer_spc = []
987 self.__buffer_spc = []
988 self.__buffer_cspc = []
988 self.__buffer_cspc = []
989 self.__buffer_dc = 0
989 self.__buffer_dc = 0
990
990
991 self.__profIndex = 0
991 self.__profIndex = 0
992 self.__dataReady = False
992 self.__dataReady = False
993 self.__byTime = False
993 self.__byTime = False
994
994
995 #self.ByLags = dataOut.ByLags ###REDEFINIR
995 #self.ByLags = dataOut.ByLags ###REDEFINIR
996 self.ByLags = False
996 self.ByLags = False
997
997
998 if DPL != None:
998 if DPL != None:
999 self.DPL=DPL
999 self.DPL=DPL
1000 else:
1000 else:
1001 #self.DPL=dataOut.DPL ###REDEFINIR
1001 #self.DPL=dataOut.DPL ###REDEFINIR
1002 self.DPL=0
1002 self.DPL=0
1003
1003
1004 if n is None and timeInterval is None:
1004 if n is None and timeInterval is None:
1005 raise ValueError("n or timeInterval should be specified ...")
1005 raise ValueError("n or timeInterval should be specified ...")
1006
1006
1007 if n is not None:
1007 if n is not None:
1008 self.n = int(n)
1008 self.n = int(n)
1009 else:
1009 else:
1010
1010
1011 self.__integrationtime = int(timeInterval)
1011 self.__integrationtime = int(timeInterval)
1012 self.n = None
1012 self.n = None
1013 self.__byTime = True
1013 self.__byTime = True
1014
1014
1015 def putData(self, data_spc, data_cspc, data_dc):
1015 def putData(self, data_spc, data_cspc, data_dc):
1016 """
1016 """
1017 Add a profile to the __buffer_spc and increase in one the __profileIndex
1017 Add a profile to the __buffer_spc and increase in one the __profileIndex
1018
1018
1019 """
1019 """
1020
1020
1021 self.__buffer_spc.append(data_spc)
1021 self.__buffer_spc.append(data_spc)
1022
1022
1023 if data_cspc is None:
1023 if data_cspc is None:
1024 self.__buffer_cspc = None
1024 self.__buffer_cspc = None
1025 else:
1025 else:
1026 self.__buffer_cspc.append(data_cspc)
1026 self.__buffer_cspc.append(data_cspc)
1027
1027
1028 if data_dc is None:
1028 if data_dc is None:
1029 self.__buffer_dc = None
1029 self.__buffer_dc = None
1030 else:
1030 else:
1031 self.__buffer_dc += data_dc
1031 self.__buffer_dc += data_dc
1032
1032
1033 self.__profIndex += 1
1033 self.__profIndex += 1
1034
1034
1035 return
1035 return
1036
1036
1037 def hildebrand_sekhon_Integration(self,data,navg):
1037 def hildebrand_sekhon_Integration(self,data,navg):
1038
1038
1039 sortdata = numpy.sort(data, axis=None)
1039 sortdata = numpy.sort(data, axis=None)
1040 sortID=data.argsort()
1040 sortID=data.argsort()
1041 lenOfData = len(sortdata)
1041 lenOfData = len(sortdata)
1042 nums_min = lenOfData*0.75
1042 nums_min = lenOfData*0.75
1043 if nums_min <= 5:
1043 if nums_min <= 5:
1044 nums_min = 5
1044 nums_min = 5
1045 sump = 0.
1045 sump = 0.
1046 sumq = 0.
1046 sumq = 0.
1047 j = 0
1047 j = 0
1048 cont = 1
1048 cont = 1
1049 while((cont == 1)and(j < lenOfData)):
1049 while((cont == 1)and(j < lenOfData)):
1050 sump += sortdata[j]
1050 sump += sortdata[j]
1051 sumq += sortdata[j]**2
1051 sumq += sortdata[j]**2
1052 if j > nums_min:
1052 if j > nums_min:
1053 rtest = float(j)/(j-1) + 1.0/navg
1053 rtest = float(j)/(j-1) + 1.0/navg
1054 if ((sumq*j) > (rtest*sump**2)):
1054 if ((sumq*j) > (rtest*sump**2)):
1055 j = j - 1
1055 j = j - 1
1056 sump = sump - sortdata[j]
1056 sump = sump - sortdata[j]
1057 sumq = sumq - sortdata[j]**2
1057 sumq = sumq - sortdata[j]**2
1058 cont = 0
1058 cont = 0
1059 j += 1
1059 j += 1
1060 #lnoise = sump / j
1060 #lnoise = sump / j
1061
1061
1062 return j,sortID
1062 return j,sortID
1063
1063
1064 def pushData(self):
1064 def pushData(self):
1065 """
1065 """
1066 Return the sum of the last profiles and the profiles used in the sum.
1066 Return the sum of the last profiles and the profiles used in the sum.
1067
1067
1068 Affected:
1068 Affected:
1069
1069
1070 self.__profileIndex
1070 self.__profileIndex
1071
1071
1072 """
1072 """
1073 bufferH=None
1073 bufferH=None
1074 buffer=None
1074 buffer=None
1075 buffer1=None
1075 buffer1=None
1076 buffer_cspc=None
1076 buffer_cspc=None
1077 self.__buffer_spc=numpy.array(self.__buffer_spc)
1077 self.__buffer_spc=numpy.array(self.__buffer_spc)
1078 self.__buffer_cspc=numpy.array(self.__buffer_cspc)
1078 self.__buffer_cspc=numpy.array(self.__buffer_cspc)
1079 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1079 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1080 #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
1080 #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
1081 for k in range(7,self.nHeights):
1081 for k in range(7,self.nHeights):
1082 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
1082 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
1083 outliers_IDs_cspc=[]
1083 outliers_IDs_cspc=[]
1084 cspc_outliers_exist=False
1084 cspc_outliers_exist=False
1085 for i in range(self.nChannels):#dataOut.nChannels):
1085 for i in range(self.nChannels):#dataOut.nChannels):
1086
1086
1087 buffer1=numpy.copy(self.__buffer_spc[:,i,:,k])
1087 buffer1=numpy.copy(self.__buffer_spc[:,i,:,k])
1088 indexes=[]
1088 indexes=[]
1089 #sortIDs=[]
1089 #sortIDs=[]
1090 outliers_IDs=[]
1090 outliers_IDs=[]
1091
1091
1092 for j in range(self.nProfiles):
1092 for j in range(self.nProfiles):
1093 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1093 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1094 # continue
1094 # continue
1095 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1095 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1096 # continue
1096 # continue
1097 buffer=buffer1[:,j]
1097 buffer=buffer1[:,j]
1098 index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
1098 index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
1099
1099
1100 indexes.append(index)
1100 indexes.append(index)
1101 #sortIDs.append(sortID)
1101 #sortIDs.append(sortID)
1102 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1102 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1103
1103
1104 outliers_IDs=numpy.array(outliers_IDs)
1104 outliers_IDs=numpy.array(outliers_IDs)
1105 outliers_IDs=outliers_IDs.ravel()
1105 outliers_IDs=outliers_IDs.ravel()
1106 outliers_IDs=numpy.unique(outliers_IDs)
1106 outliers_IDs=numpy.unique(outliers_IDs)
1107 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
1107 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
1108 indexes=numpy.array(indexes)
1108 indexes=numpy.array(indexes)
1109 indexmin=numpy.min(indexes)
1109 indexmin=numpy.min(indexes)
1110
1110
1111 if indexmin != buffer1.shape[0]:
1111 if indexmin != buffer1.shape[0]:
1112 cspc_outliers_exist=True
1112 cspc_outliers_exist=True
1113 ###sortdata=numpy.sort(buffer1,axis=0)
1113 ###sortdata=numpy.sort(buffer1,axis=0)
1114 ###avg2=numpy.mean(sortdata[:indexmin,:],axis=0)
1114 ###avg2=numpy.mean(sortdata[:indexmin,:],axis=0)
1115 lt=outliers_IDs
1115 lt=outliers_IDs
1116 avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1116 avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1117
1117
1118 for p in list(outliers_IDs):
1118 for p in list(outliers_IDs):
1119 buffer1[p,:]=avg
1119 buffer1[p,:]=avg
1120
1120
1121 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1121 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1122 ###cspc IDs
1122 ###cspc IDs
1123 #indexmin_cspc+=indexmin_cspc
1123 #indexmin_cspc+=indexmin_cspc
1124 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1124 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1125
1125
1126 #if not breakFlag:
1126 #if not breakFlag:
1127 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1127 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1128 if cspc_outliers_exist:
1128 if cspc_outliers_exist:
1129 #sortdata=numpy.sort(buffer_cspc,axis=0)
1129 #sortdata=numpy.sort(buffer_cspc,axis=0)
1130 #avg=numpy.mean(sortdata[:indexmin_cpsc,:],axis=0)
1130 #avg=numpy.mean(sortdata[:indexmin_cpsc,:],axis=0)
1131 lt=outliers_IDs_cspc
1131 lt=outliers_IDs_cspc
1132
1132
1133 avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
1133 avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
1134 for p in list(outliers_IDs_cspc):
1134 for p in list(outliers_IDs_cspc):
1135 buffer_cspc[p,:]=avg
1135 buffer_cspc[p,:]=avg
1136
1136
1137 self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
1137 self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
1138 #else:
1138 #else:
1139 #break
1139 #break
1140
1140
1141
1141
1142
1142
1143
1143
1144 buffer=None
1144 buffer=None
1145 bufferH=None
1145 bufferH=None
1146 buffer1=None
1146 buffer1=None
1147 buffer_cspc=None
1147 buffer_cspc=None
1148
1148
1149 #print("cpsc",self.__buffer_cspc[:,0,0,0,0])
1149 #print("cpsc",self.__buffer_cspc[:,0,0,0,0])
1150 #print(self.__profIndex)
1150 #print(self.__profIndex)
1151 #exit()
1151 #exit()
1152
1152
1153 buffer=None
1153 buffer=None
1154 #print(self.__buffer_spc[:,1,3,20,0])
1154 #print(self.__buffer_spc[:,1,3,20,0])
1155 #print(self.__buffer_spc[:,1,5,37,0])
1155 #print(self.__buffer_spc[:,1,5,37,0])
1156 data_spc = numpy.sum(self.__buffer_spc,axis=0)
1156 data_spc = numpy.sum(self.__buffer_spc,axis=0)
1157 data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
1157 data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
1158
1158
1159 #print(numpy.shape(data_spc))
1159 #print(numpy.shape(data_spc))
1160 #data_spc[1,4,20,0]=numpy.nan
1160 #data_spc[1,4,20,0]=numpy.nan
1161
1161
1162 #data_cspc = self.__buffer_cspc
1162 #data_cspc = self.__buffer_cspc
1163 data_dc = self.__buffer_dc
1163 data_dc = self.__buffer_dc
1164 n = self.__profIndex
1164 n = self.__profIndex
1165
1165
1166 self.__buffer_spc = []
1166 self.__buffer_spc = []
1167 self.__buffer_cspc = []
1167 self.__buffer_cspc = []
1168 self.__buffer_dc = 0
1168 self.__buffer_dc = 0
1169 self.__profIndex = 0
1169 self.__profIndex = 0
1170
1170
1171 return data_spc, data_cspc, data_dc, n
1171 return data_spc, data_cspc, data_dc, n
1172
1172
1173 def byProfiles(self, *args):
1173 def byProfiles(self, *args):
1174
1174
1175 self.__dataReady = False
1175 self.__dataReady = False
1176 avgdata_spc = None
1176 avgdata_spc = None
1177 avgdata_cspc = None
1177 avgdata_cspc = None
1178 avgdata_dc = None
1178 avgdata_dc = None
1179
1179
1180 self.putData(*args)
1180 self.putData(*args)
1181
1181
1182 if self.__profIndex == self.n:
1182 if self.__profIndex == self.n:
1183
1183
1184 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1184 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1185 self.n = n
1185 self.n = n
1186 self.__dataReady = True
1186 self.__dataReady = True
1187
1187
1188 return avgdata_spc, avgdata_cspc, avgdata_dc
1188 return avgdata_spc, avgdata_cspc, avgdata_dc
1189
1189
1190 def byTime(self, datatime, *args):
1190 def byTime(self, datatime, *args):
1191
1191
1192 self.__dataReady = False
1192 self.__dataReady = False
1193 avgdata_spc = None
1193 avgdata_spc = None
1194 avgdata_cspc = None
1194 avgdata_cspc = None
1195 avgdata_dc = None
1195 avgdata_dc = None
1196
1196
1197 self.putData(*args)
1197 self.putData(*args)
1198
1198
1199 if (datatime - self.__initime) >= self.__integrationtime:
1199 if (datatime - self.__initime) >= self.__integrationtime:
1200 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1200 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1201 self.n = n
1201 self.n = n
1202 self.__dataReady = True
1202 self.__dataReady = True
1203
1203
1204 return avgdata_spc, avgdata_cspc, avgdata_dc
1204 return avgdata_spc, avgdata_cspc, avgdata_dc
1205
1205
1206 def integrate(self, datatime, *args):
1206 def integrate(self, datatime, *args):
1207
1207
1208 if self.__profIndex == 0:
1208 if self.__profIndex == 0:
1209 self.__initime = datatime
1209 self.__initime = datatime
1210
1210
1211 if self.__byTime:
1211 if self.__byTime:
1212 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1212 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1213 datatime, *args)
1213 datatime, *args)
1214 else:
1214 else:
1215 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1215 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1216
1216
1217 if not self.__dataReady:
1217 if not self.__dataReady:
1218 return None, None, None, None
1218 return None, None, None, None
1219
1219
1220 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1220 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1221
1221
1222 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False):
1222 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False):
1223 if n == 1:
1223 if n == 1:
1224 return dataOut
1224 return dataOut
1225
1225
1226 dataOut.flagNoData = True
1226 dataOut.flagNoData = True
1227
1227
1228 if not self.isConfig:
1228 if not self.isConfig:
1229 self.setup(dataOut, n, timeInterval, overlapping,DPL )
1229 self.setup(dataOut, n, timeInterval, overlapping,DPL )
1230 self.isConfig = True
1230 self.isConfig = True
1231
1231
1232 if not self.ByLags:
1232 if not self.ByLags:
1233 self.nProfiles=dataOut.nProfiles
1233 self.nProfiles=dataOut.nProfiles
1234 self.nChannels=dataOut.nChannels
1234 self.nChannels=dataOut.nChannels
1235 self.nHeights=dataOut.nHeights
1235 self.nHeights=dataOut.nHeights
1236 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1236 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1237 dataOut.data_spc,
1237 dataOut.data_spc,
1238 dataOut.data_cspc,
1238 dataOut.data_cspc,
1239 dataOut.data_dc)
1239 dataOut.data_dc)
1240 else:
1240 else:
1241 self.nProfiles=dataOut.nProfiles
1241 self.nProfiles=dataOut.nProfiles
1242 self.nChannels=dataOut.nChannels
1242 self.nChannels=dataOut.nChannels
1243 self.nHeights=dataOut.nHeights
1243 self.nHeights=dataOut.nHeights
1244 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1244 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1245 dataOut.dataLag_spc,
1245 dataOut.dataLag_spc,
1246 dataOut.dataLag_cspc,
1246 dataOut.dataLag_cspc,
1247 dataOut.dataLag_dc)
1247 dataOut.dataLag_dc)
1248
1248
1249 if self.__dataReady:
1249 if self.__dataReady:
1250
1250
1251 if not self.ByLags:
1251 if not self.ByLags:
1252
1252
1253 dataOut.data_spc = numpy.squeeze(avgdata_spc)
1253 dataOut.data_spc = numpy.squeeze(avgdata_spc)
1254 dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1254 dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1255 dataOut.data_dc = avgdata_dc
1255 dataOut.data_dc = avgdata_dc
1256 else:
1256 else:
1257 dataOut.dataLag_spc = avgdata_spc
1257 dataOut.dataLag_spc = avgdata_spc
1258 dataOut.dataLag_cspc = avgdata_cspc
1258 dataOut.dataLag_cspc = avgdata_cspc
1259 dataOut.dataLag_dc = avgdata_dc
1259 dataOut.dataLag_dc = avgdata_dc
1260
1260
1261 dataOut.data_spc=dataOut.dataLag_spc[:,:,:,dataOut.LagPlot]
1261 dataOut.data_spc=dataOut.dataLag_spc[:,:,:,dataOut.LagPlot]
1262 dataOut.data_cspc=dataOut.dataLag_cspc[:,:,:,dataOut.LagPlot]
1262 dataOut.data_cspc=dataOut.dataLag_cspc[:,:,:,dataOut.LagPlot]
1263 dataOut.data_dc=dataOut.dataLag_dc[:,:,dataOut.LagPlot]
1263 dataOut.data_dc=dataOut.dataLag_dc[:,:,dataOut.LagPlot]
1264
1264
1265
1265
1266 dataOut.nIncohInt *= self.n
1266 dataOut.nIncohInt *= self.n
1267 dataOut.utctime = avgdatatime
1267 dataOut.utctime = avgdatatime
1268 dataOut.flagNoData = False
1268 dataOut.flagNoData = False
1269
1269
1270 return dataOut
1270 return dataOut
1271
1271
1272 class removeInterference(Operation):
1272 class removeInterference(Operation):
1273
1273
1274 def removeInterference2(self):
1274 def removeInterference2(self):
1275
1275
1276 cspc = self.dataOut.data_cspc
1276 cspc = self.dataOut.data_cspc
1277 spc = self.dataOut.data_spc
1277 spc = self.dataOut.data_spc
1278 Heights = numpy.arange(cspc.shape[2])
1278 Heights = numpy.arange(cspc.shape[2])
1279 realCspc = numpy.abs(cspc)
1279 realCspc = numpy.abs(cspc)
1280
1280
1281 for i in range(cspc.shape[0]):
1281 for i in range(cspc.shape[0]):
1282 LinePower= numpy.sum(realCspc[i], axis=0)
1282 LinePower= numpy.sum(realCspc[i], axis=0)
1283 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
1283 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
1284 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
1284 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
1285 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
1285 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
1286 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
1286 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
1287 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
1287 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
1288
1288
1289
1289
1290 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
1290 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
1291 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
1291 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
1292 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
1292 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
1293 cspc[i,InterferenceRange,:] = numpy.NaN
1293 cspc[i,InterferenceRange,:] = numpy.NaN
1294
1294
1295 self.dataOut.data_cspc = cspc
1295 self.dataOut.data_cspc = cspc
1296
1296
1297 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
1297 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
1298
1298
1299 jspectra = self.dataOut.data_spc
1299 jspectra = self.dataOut.data_spc
1300 jcspectra = self.dataOut.data_cspc
1300 jcspectra = self.dataOut.data_cspc
1301 jnoise = self.dataOut.getNoise()
1301 jnoise = self.dataOut.getNoise()
1302 num_incoh = self.dataOut.nIncohInt
1302 num_incoh = self.dataOut.nIncohInt
1303
1303
1304 num_channel = jspectra.shape[0]
1304 num_channel = jspectra.shape[0]
1305 num_prof = jspectra.shape[1]
1305 num_prof = jspectra.shape[1]
1306 num_hei = jspectra.shape[2]
1306 num_hei = jspectra.shape[2]
1307
1307
1308 # hei_interf
1308 # hei_interf
1309 if hei_interf is None:
1309 if hei_interf is None:
1310 count_hei = int(num_hei / 2)
1310 count_hei = int(num_hei / 2)
1311 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
1311 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
1312 hei_interf = numpy.asarray(hei_interf)[0]
1312 hei_interf = numpy.asarray(hei_interf)[0]
1313 # nhei_interf
1313 # nhei_interf
1314 if (nhei_interf == None):
1314 if (nhei_interf == None):
1315 nhei_interf = 5
1315 nhei_interf = 5
1316 if (nhei_interf < 1):
1316 if (nhei_interf < 1):
1317 nhei_interf = 1
1317 nhei_interf = 1
1318 if (nhei_interf > count_hei):
1318 if (nhei_interf > count_hei):
1319 nhei_interf = count_hei
1319 nhei_interf = count_hei
1320 if (offhei_interf == None):
1320 if (offhei_interf == None):
1321 offhei_interf = 0
1321 offhei_interf = 0
1322
1322
1323 ind_hei = list(range(num_hei))
1323 ind_hei = list(range(num_hei))
1324 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
1324 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
1325 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
1325 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
1326 mask_prof = numpy.asarray(list(range(num_prof)))
1326 mask_prof = numpy.asarray(list(range(num_prof)))
1327 num_mask_prof = mask_prof.size
1327 num_mask_prof = mask_prof.size
1328 comp_mask_prof = [0, num_prof / 2]
1328 comp_mask_prof = [0, num_prof / 2]
1329
1329
1330 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
1330 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
1331 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
1331 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
1332 jnoise = numpy.nan
1332 jnoise = numpy.nan
1333 noise_exist = jnoise[0] < numpy.Inf
1333 noise_exist = jnoise[0] < numpy.Inf
1334
1334
1335 # Subrutina de Remocion de la Interferencia
1335 # Subrutina de Remocion de la Interferencia
1336 for ich in range(num_channel):
1336 for ich in range(num_channel):
1337 # Se ordena los espectros segun su potencia (menor a mayor)
1337 # Se ordena los espectros segun su potencia (menor a mayor)
1338 power = jspectra[ich, mask_prof, :]
1338 power = jspectra[ich, mask_prof, :]
1339 power = power[:, hei_interf]
1339 power = power[:, hei_interf]
1340 power = power.sum(axis=0)
1340 power = power.sum(axis=0)
1341 psort = power.ravel().argsort()
1341 psort = power.ravel().argsort()
1342
1342
1343 # Se estima la interferencia promedio en los Espectros de Potencia empleando
1343 # Se estima la interferencia promedio en los Espectros de Potencia empleando
1344 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
1344 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
1345 offhei_interf, nhei_interf + offhei_interf))]]]
1345 offhei_interf, nhei_interf + offhei_interf))]]]
1346
1346
1347 if noise_exist:
1347 if noise_exist:
1348 # tmp_noise = jnoise[ich] / num_prof
1348 # tmp_noise = jnoise[ich] / num_prof
1349 tmp_noise = jnoise[ich]
1349 tmp_noise = jnoise[ich]
1350 junkspc_interf = junkspc_interf - tmp_noise
1350 junkspc_interf = junkspc_interf - tmp_noise
1351 #junkspc_interf[:,comp_mask_prof] = 0
1351 #junkspc_interf[:,comp_mask_prof] = 0
1352
1352
1353 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
1353 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
1354 jspc_interf = jspc_interf.transpose()
1354 jspc_interf = jspc_interf.transpose()
1355 # Calculando el espectro de interferencia promedio
1355 # Calculando el espectro de interferencia promedio
1356 noiseid = numpy.where(
1356 noiseid = numpy.where(
1357 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
1357 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
1358 noiseid = noiseid[0]
1358 noiseid = noiseid[0]
1359 cnoiseid = noiseid.size
1359 cnoiseid = noiseid.size
1360 interfid = numpy.where(
1360 interfid = numpy.where(
1361 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
1361 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
1362 interfid = interfid[0]
1362 interfid = interfid[0]
1363 cinterfid = interfid.size
1363 cinterfid = interfid.size
1364
1364
1365 if (cnoiseid > 0):
1365 if (cnoiseid > 0):
1366 jspc_interf[noiseid] = 0
1366 jspc_interf[noiseid] = 0
1367
1367
1368 # Expandiendo los perfiles a limpiar
1368 # Expandiendo los perfiles a limpiar
1369 if (cinterfid > 0):
1369 if (cinterfid > 0):
1370 new_interfid = (
1370 new_interfid = (
1371 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
1371 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
1372 new_interfid = numpy.asarray(new_interfid)
1372 new_interfid = numpy.asarray(new_interfid)
1373 new_interfid = {x for x in new_interfid}
1373 new_interfid = {x for x in new_interfid}
1374 new_interfid = numpy.array(list(new_interfid))
1374 new_interfid = numpy.array(list(new_interfid))
1375 new_cinterfid = new_interfid.size
1375 new_cinterfid = new_interfid.size
1376 else:
1376 else:
1377 new_cinterfid = 0
1377 new_cinterfid = 0
1378
1378
1379 for ip in range(new_cinterfid):
1379 for ip in range(new_cinterfid):
1380 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
1380 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
1381 jspc_interf[new_interfid[ip]
1381 jspc_interf[new_interfid[ip]
1382 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
1382 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
1383
1383
1384 jspectra[ich, :, ind_hei] = jspectra[ich, :,
1384 jspectra[ich, :, ind_hei] = jspectra[ich, :,
1385 ind_hei] - jspc_interf # Corregir indices
1385 ind_hei] - jspc_interf # Corregir indices
1386
1386
1387 # Removiendo la interferencia del punto de mayor interferencia
1387 # Removiendo la interferencia del punto de mayor interferencia
1388 ListAux = jspc_interf[mask_prof].tolist()
1388 ListAux = jspc_interf[mask_prof].tolist()
1389 maxid = ListAux.index(max(ListAux))
1389 maxid = ListAux.index(max(ListAux))
1390
1390
1391 if cinterfid > 0:
1391 if cinterfid > 0:
1392 for ip in range(cinterfid * (interf == 2) - 1):
1392 for ip in range(cinterfid * (interf == 2) - 1):
1393 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
1393 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
1394 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
1394 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
1395 cind = len(ind)
1395 cind = len(ind)
1396
1396
1397 if (cind > 0):
1397 if (cind > 0):
1398 jspectra[ich, interfid[ip], ind] = tmp_noise * \
1398 jspectra[ich, interfid[ip], ind] = tmp_noise * \
1399 (1 + (numpy.random.uniform(cind) - 0.5) /
1399 (1 + (numpy.random.uniform(cind) - 0.5) /
1400 numpy.sqrt(num_incoh))
1400 numpy.sqrt(num_incoh))
1401
1401
1402 ind = numpy.array([-2, -1, 1, 2])
1402 ind = numpy.array([-2, -1, 1, 2])
1403 xx = numpy.zeros([4, 4])
1403 xx = numpy.zeros([4, 4])
1404
1404
1405 for id1 in range(4):
1405 for id1 in range(4):
1406 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1406 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1407
1407
1408 xx_inv = numpy.linalg.inv(xx)
1408 xx_inv = numpy.linalg.inv(xx)
1409 xx = xx_inv[:, 0]
1409 xx = xx_inv[:, 0]
1410 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1410 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1411 yy = jspectra[ich, mask_prof[ind], :]
1411 yy = jspectra[ich, mask_prof[ind], :]
1412 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
1412 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
1413 yy.transpose(), xx)
1413 yy.transpose(), xx)
1414
1414
1415 indAux = (jspectra[ich, :, :] < tmp_noise *
1415 indAux = (jspectra[ich, :, :] < tmp_noise *
1416 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
1416 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
1417 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
1417 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
1418 (1 - 1 / numpy.sqrt(num_incoh))
1418 (1 - 1 / numpy.sqrt(num_incoh))
1419
1419
1420 # Remocion de Interferencia en el Cross Spectra
1420 # Remocion de Interferencia en el Cross Spectra
1421 if jcspectra is None:
1421 if jcspectra is None:
1422 return jspectra, jcspectra
1422 return jspectra, jcspectra
1423 num_pairs = int(jcspectra.size / (num_prof * num_hei))
1423 num_pairs = int(jcspectra.size / (num_prof * num_hei))
1424 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
1424 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
1425
1425
1426 for ip in range(num_pairs):
1426 for ip in range(num_pairs):
1427
1427
1428 #-------------------------------------------
1428 #-------------------------------------------
1429
1429
1430 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
1430 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
1431 cspower = cspower[:, hei_interf]
1431 cspower = cspower[:, hei_interf]
1432 cspower = cspower.sum(axis=0)
1432 cspower = cspower.sum(axis=0)
1433
1433
1434 cspsort = cspower.ravel().argsort()
1434 cspsort = cspower.ravel().argsort()
1435 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
1435 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
1436 offhei_interf, nhei_interf + offhei_interf))]]]
1436 offhei_interf, nhei_interf + offhei_interf))]]]
1437 junkcspc_interf = junkcspc_interf.transpose()
1437 junkcspc_interf = junkcspc_interf.transpose()
1438 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
1438 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
1439
1439
1440 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1440 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1441
1441
1442 median_real = int(numpy.median(numpy.real(
1442 median_real = int(numpy.median(numpy.real(
1443 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1443 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1444 median_imag = int(numpy.median(numpy.imag(
1444 median_imag = int(numpy.median(numpy.imag(
1445 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1445 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1446 comp_mask_prof = [int(e) for e in comp_mask_prof]
1446 comp_mask_prof = [int(e) for e in comp_mask_prof]
1447 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
1447 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
1448 median_real, median_imag)
1448 median_real, median_imag)
1449
1449
1450 for iprof in range(num_prof):
1450 for iprof in range(num_prof):
1451 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
1451 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
1452 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
1452 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
1453
1453
1454 # Removiendo la Interferencia
1454 # Removiendo la Interferencia
1455 jcspectra[ip, :, ind_hei] = jcspectra[ip,
1455 jcspectra[ip, :, ind_hei] = jcspectra[ip,
1456 :, ind_hei] - jcspc_interf
1456 :, ind_hei] - jcspc_interf
1457
1457
1458 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1458 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1459 maxid = ListAux.index(max(ListAux))
1459 maxid = ListAux.index(max(ListAux))
1460
1460
1461 ind = numpy.array([-2, -1, 1, 2])
1461 ind = numpy.array([-2, -1, 1, 2])
1462 xx = numpy.zeros([4, 4])
1462 xx = numpy.zeros([4, 4])
1463
1463
1464 for id1 in range(4):
1464 for id1 in range(4):
1465 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1465 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1466
1466
1467 xx_inv = numpy.linalg.inv(xx)
1467 xx_inv = numpy.linalg.inv(xx)
1468 xx = xx_inv[:, 0]
1468 xx = xx_inv[:, 0]
1469
1469
1470 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1470 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1471 yy = jcspectra[ip, mask_prof[ind], :]
1471 yy = jcspectra[ip, mask_prof[ind], :]
1472 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1472 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1473
1473
1474 # Guardar Resultados
1474 # Guardar Resultados
1475 self.dataOut.data_spc = jspectra
1475 self.dataOut.data_spc = jspectra
1476 self.dataOut.data_cspc = jcspectra
1476 self.dataOut.data_cspc = jcspectra
1477
1477
1478 return 1
1478 return 1
1479
1479
1480 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
1480 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
1481
1481
1482 self.dataOut = dataOut
1482 self.dataOut = dataOut
1483
1483
1484 if mode == 1:
1484 if mode == 1:
1485 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
1485 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
1486 elif mode == 2:
1486 elif mode == 2:
1487 self.removeInterference2()
1487 self.removeInterference2()
1488
1488
1489 return self.dataOut
1489 return self.dataOut
1490
1490
1491
1491
1492 class IncohInt(Operation):
1492 class IncohInt(Operation):
1493
1493
1494 __profIndex = 0
1494 __profIndex = 0
1495 __withOverapping = False
1495 __withOverapping = False
1496
1496
1497 __byTime = False
1497 __byTime = False
1498 __initime = None
1498 __initime = None
1499 __lastdatatime = None
1499 __lastdatatime = None
1500 __integrationtime = None
1500 __integrationtime = None
1501
1501
1502 __buffer_spc = None
1502 __buffer_spc = None
1503 __buffer_cspc = None
1503 __buffer_cspc = None
1504 __buffer_dc = None
1504 __buffer_dc = None
1505
1505
1506 __dataReady = False
1506 __dataReady = False
1507
1507
1508 __timeInterval = None
1508 __timeInterval = None
1509
1509
1510 n = None
1510 n = None
1511
1511
1512 def __init__(self):
1512 def __init__(self):
1513
1513
1514 Operation.__init__(self)
1514 Operation.__init__(self)
1515
1515
1516 def setup(self, n=None, timeInterval=None, overlapping=False):
1516 def setup(self, n=None, timeInterval=None, overlapping=False):
1517 """
1517 """
1518 Set the parameters of the integration class.
1518 Set the parameters of the integration class.
1519
1519
1520 Inputs:
1520 Inputs:
1521
1521
1522 n : Number of coherent integrations
1522 n : Number of coherent integrations
1523 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1523 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1524 overlapping :
1524 overlapping :
1525
1525
1526 """
1526 """
1527
1527
1528 self.__initime = None
1528 self.__initime = None
1529 self.__lastdatatime = 0
1529 self.__lastdatatime = 0
1530
1530
1531 self.__buffer_spc = 0
1531 self.__buffer_spc = 0
1532 self.__buffer_cspc = 0
1532 self.__buffer_cspc = 0
1533 self.__buffer_dc = 0
1533 self.__buffer_dc = 0
1534
1534
1535 self.__profIndex = 0
1535 self.__profIndex = 0
1536 self.__dataReady = False
1536 self.__dataReady = False
1537 self.__byTime = False
1537 self.__byTime = False
1538
1538
1539 if n is None and timeInterval is None:
1539 if n is None and timeInterval is None:
1540 raise ValueError("n or timeInterval should be specified ...")
1540 raise ValueError("n or timeInterval should be specified ...")
1541
1541
1542 if n is not None:
1542 if n is not None:
1543 self.n = int(n)
1543 self.n = int(n)
1544 else:
1544 else:
1545
1545
1546 self.__integrationtime = int(timeInterval)
1546 self.__integrationtime = int(timeInterval)
1547 self.n = None
1547 self.n = None
1548 self.__byTime = True
1548 self.__byTime = True
1549
1549
1550 def putData(self, data_spc, data_cspc, data_dc):
1550 def putData(self, data_spc, data_cspc, data_dc):
1551 """
1551 """
1552 Add a profile to the __buffer_spc and increase in one the __profileIndex
1552 Add a profile to the __buffer_spc and increase in one the __profileIndex
1553
1553
1554 """
1554 """
1555
1555
1556 self.__buffer_spc += data_spc
1556 self.__buffer_spc += data_spc
1557
1557
1558 if data_cspc is None:
1558 if data_cspc is None:
1559 self.__buffer_cspc = None
1559 self.__buffer_cspc = None
1560 else:
1560 else:
1561 self.__buffer_cspc += data_cspc
1561 self.__buffer_cspc += data_cspc
1562
1562
1563 if data_dc is None:
1563 if data_dc is None:
1564 self.__buffer_dc = None
1564 self.__buffer_dc = None
1565 else:
1565 else:
1566 self.__buffer_dc += data_dc
1566 self.__buffer_dc += data_dc
1567
1567
1568 self.__profIndex += 1
1568 self.__profIndex += 1
1569
1569
1570 return
1570 return
1571
1571
1572 def pushData(self):
1572 def pushData(self):
1573 """
1573 """
1574 Return the sum of the last profiles and the profiles used in the sum.
1574 Return the sum of the last profiles and the profiles used in the sum.
1575
1575
1576 Affected:
1576 Affected:
1577
1577
1578 self.__profileIndex
1578 self.__profileIndex
1579
1579
1580 """
1580 """
1581
1581
1582 data_spc = self.__buffer_spc
1582 data_spc = self.__buffer_spc
1583 data_cspc = self.__buffer_cspc
1583 data_cspc = self.__buffer_cspc
1584 data_dc = self.__buffer_dc
1584 data_dc = self.__buffer_dc
1585 n = self.__profIndex
1585 n = self.__profIndex
1586
1586
1587 self.__buffer_spc = 0
1587 self.__buffer_spc = 0
1588 self.__buffer_cspc = 0
1588 self.__buffer_cspc = 0
1589 self.__buffer_dc = 0
1589 self.__buffer_dc = 0
1590 self.__profIndex = 0
1590 self.__profIndex = 0
1591
1591
1592 return data_spc, data_cspc, data_dc, n
1592 return data_spc, data_cspc, data_dc, n
1593
1593
1594 def byProfiles(self, *args):
1594 def byProfiles(self, *args):
1595
1595
1596 self.__dataReady = False
1596 self.__dataReady = False
1597 avgdata_spc = None
1597 avgdata_spc = None
1598 avgdata_cspc = None
1598 avgdata_cspc = None
1599 avgdata_dc = None
1599 avgdata_dc = None
1600
1600
1601 self.putData(*args)
1601 self.putData(*args)
1602
1602
1603 if self.__profIndex == self.n:
1603 if self.__profIndex == self.n:
1604
1604
1605 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1605 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1606 self.n = n
1606 self.n = n
1607 self.__dataReady = True
1607 self.__dataReady = True
1608
1608
1609 return avgdata_spc, avgdata_cspc, avgdata_dc
1609 return avgdata_spc, avgdata_cspc, avgdata_dc
1610
1610
1611 def byTime(self, datatime, *args):
1611 def byTime(self, datatime, *args):
1612
1612
1613 self.__dataReady = False
1613 self.__dataReady = False
1614 avgdata_spc = None
1614 avgdata_spc = None
1615 avgdata_cspc = None
1615 avgdata_cspc = None
1616 avgdata_dc = None
1616 avgdata_dc = None
1617
1617
1618 self.putData(*args)
1618 self.putData(*args)
1619
1619
1620 if (datatime - self.__initime) >= self.__integrationtime:
1620 if (datatime - self.__initime) >= self.__integrationtime:
1621 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1621 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1622 self.n = n
1622 self.n = n
1623 self.__dataReady = True
1623 self.__dataReady = True
1624
1624
1625 return avgdata_spc, avgdata_cspc, avgdata_dc
1625 return avgdata_spc, avgdata_cspc, avgdata_dc
1626
1626
1627 def integrate(self, datatime, *args):
1627 def integrate(self, datatime, *args):
1628
1628
1629 if self.__profIndex == 0:
1629 if self.__profIndex == 0:
1630 self.__initime = datatime
1630 self.__initime = datatime
1631
1631
1632 if self.__byTime:
1632 if self.__byTime:
1633 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1633 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1634 datatime, *args)
1634 datatime, *args)
1635 else:
1635 else:
1636 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1636 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1637
1637
1638 if not self.__dataReady:
1638 if not self.__dataReady:
1639 return None, None, None, None
1639 return None, None, None, None
1640
1640
1641 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1641 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1642
1642
1643 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1643 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1644 if n == 1:
1644 if n == 1:
1645 return dataOut
1645 return dataOut
1646
1646
1647 dataOut.flagNoData = True
1647 dataOut.flagNoData = True
1648
1648
1649 if not self.isConfig:
1649 if not self.isConfig:
1650 self.setup(n, timeInterval, overlapping)
1650 self.setup(n, timeInterval, overlapping)
1651 self.isConfig = True
1651 self.isConfig = True
1652
1652
1653 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1653 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1654 dataOut.data_spc,
1654 dataOut.data_spc,
1655 dataOut.data_cspc,
1655 dataOut.data_cspc,
1656 dataOut.data_dc)
1656 dataOut.data_dc)
1657
1657
1658 if self.__dataReady:
1658 if self.__dataReady:
1659
1659
1660 dataOut.data_spc = avgdata_spc
1660 dataOut.data_spc = avgdata_spc
1661 dataOut.data_cspc = avgdata_cspc
1661 dataOut.data_cspc = avgdata_cspc
1662 dataOut.data_dc = avgdata_dc
1662 dataOut.data_dc = avgdata_dc
1663 dataOut.nIncohInt *= self.n
1663 dataOut.nIncohInt *= self.n
1664 dataOut.utctime = avgdatatime
1664 dataOut.utctime = avgdatatime
1665 dataOut.flagNoData = False
1665 dataOut.flagNoData = False
1666
1666
1667 return dataOut
1667 return dataOut
1668
1668
1669 class dopplerFlip(Operation):
1669 class dopplerFlip(Operation):
1670
1670
1671 def run(self, dataOut):
1671 def run(self, dataOut):
1672 # arreglo 1: (num_chan, num_profiles, num_heights)
1672 # arreglo 1: (num_chan, num_profiles, num_heights)
1673 self.dataOut = dataOut
1673 self.dataOut = dataOut
1674 # JULIA-oblicua, indice 2
1674 # JULIA-oblicua, indice 2
1675 # arreglo 2: (num_profiles, num_heights)
1675 # arreglo 2: (num_profiles, num_heights)
1676 jspectra = self.dataOut.data_spc[2]
1676 jspectra = self.dataOut.data_spc[2]
1677 jspectra_tmp = numpy.zeros(jspectra.shape)
1677 jspectra_tmp = numpy.zeros(jspectra.shape)
1678 num_profiles = jspectra.shape[0]
1678 num_profiles = jspectra.shape[0]
1679 freq_dc = int(num_profiles / 2)
1679 freq_dc = int(num_profiles / 2)
1680 # Flip con for
1680 # Flip con for
1681 for j in range(num_profiles):
1681 for j in range(num_profiles):
1682 jspectra_tmp[num_profiles-j-1]= jspectra[j]
1682 jspectra_tmp[num_profiles-j-1]= jspectra[j]
1683 # Intercambio perfil de DC con perfil inmediato anterior
1683 # Intercambio perfil de DC con perfil inmediato anterior
1684 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
1684 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
1685 jspectra_tmp[freq_dc]= jspectra[freq_dc]
1685 jspectra_tmp[freq_dc]= jspectra[freq_dc]
1686 # canal modificado es re-escrito en el arreglo de canales
1686 # canal modificado es re-escrito en el arreglo de canales
1687 self.dataOut.data_spc[2] = jspectra_tmp
1687 self.dataOut.data_spc[2] = jspectra_tmp
1688
1688
1689 return self.dataOut
1689 return self.dataOut
@@ -1,1633 +1,1638
1 import sys
1 import sys
2 import numpy,math
2 import numpy,math
3 from scipy import interpolate
3 from scipy import interpolate
4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
5 from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
6 from schainpy.utils import log
6 from schainpy.utils import log
7 from time import time
7 from time import time
8
8 import numpy
9
9
10
10
11 class VoltageProc(ProcessingUnit):
11 class VoltageProc(ProcessingUnit):
12
12
13 def __init__(self):
13 def __init__(self):
14
14
15 ProcessingUnit.__init__(self)
15 ProcessingUnit.__init__(self)
16
16
17 self.dataOut = Voltage()
17 self.dataOut = Voltage()
18 self.flip = 1
18 self.flip = 1
19 self.setupReq = False
19 self.setupReq = False
20
20
21 def run(self):
21 def run(self):
22
22
23 if self.dataIn.type == 'AMISR':
23 if self.dataIn.type == 'AMISR':
24 self.__updateObjFromAmisrInput()
24 self.__updateObjFromAmisrInput()
25
25
26 if self.dataIn.type == 'Voltage':
26 if self.dataIn.type == 'Voltage':
27 self.dataOut.copy(self.dataIn)
27 self.dataOut.copy(self.dataIn)
28
28
29 def __updateObjFromAmisrInput(self):
29 def __updateObjFromAmisrInput(self):
30
30
31 self.dataOut.timeZone = self.dataIn.timeZone
31 self.dataOut.timeZone = self.dataIn.timeZone
32 self.dataOut.dstFlag = self.dataIn.dstFlag
32 self.dataOut.dstFlag = self.dataIn.dstFlag
33 self.dataOut.errorCount = self.dataIn.errorCount
33 self.dataOut.errorCount = self.dataIn.errorCount
34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
35
35
36 self.dataOut.flagNoData = self.dataIn.flagNoData
36 self.dataOut.flagNoData = self.dataIn.flagNoData
37 self.dataOut.data = self.dataIn.data
37 self.dataOut.data = self.dataIn.data
38 self.dataOut.utctime = self.dataIn.utctime
38 self.dataOut.utctime = self.dataIn.utctime
39 self.dataOut.channelList = self.dataIn.channelList
39 self.dataOut.channelList = self.dataIn.channelList
40 #self.dataOut.timeInterval = self.dataIn.timeInterval
40 #self.dataOut.timeInterval = self.dataIn.timeInterval
41 self.dataOut.heightList = self.dataIn.heightList
41 self.dataOut.heightList = self.dataIn.heightList
42 self.dataOut.nProfiles = self.dataIn.nProfiles
42 self.dataOut.nProfiles = self.dataIn.nProfiles
43
43
44 self.dataOut.nCohInt = self.dataIn.nCohInt
44 self.dataOut.nCohInt = self.dataIn.nCohInt
45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
46 self.dataOut.frequency = self.dataIn.frequency
46 self.dataOut.frequency = self.dataIn.frequency
47
47
48 self.dataOut.azimuth = self.dataIn.azimuth
48 self.dataOut.azimuth = self.dataIn.azimuth
49 self.dataOut.zenith = self.dataIn.zenith
49 self.dataOut.zenith = self.dataIn.zenith
50
50
51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
54
54
55
55
56 class selectChannels(Operation):
56 class selectChannels(Operation):
57
57
58 def run(self, dataOut, channelList=None):
58 def run(self, dataOut, channelList=None):
59 self.channelList = channelList
59 self.channelList = channelList
60 if self.channelList == None:
60 if self.channelList == None:
61 print("Missing channelList")
61 print("Missing channelList")
62 return dataOut
62 return dataOut
63 channelIndexList = []
63 channelIndexList = []
64 if type(self.channelList) is list:
65 pass
66 elif type(self.channelList) is tuple:
67 self.channelList = list(self.channelList)
68 else:
69 self.channelList = self.channelList.tolist()
70
71 if type(dataOut.channelList) is not list:
72 dataOut.channelList = dataOut.channelList.tolist()
73
64
65 if type(dataOut.channelList) is not list: #leer array desde HDF5
66 try:
67 dataOut.channelList = dataOut.channelList.tolist()
68 except Exception as e:
69 print("Select Channels: ",e)
74 for channel in self.channelList:
70 for channel in self.channelList:
75 if channel not in dataOut.channelList:
71 if channel not in dataOut.channelList:
76 raise ValueError("Channel %d is not in %s" %(channel, str(dataOut.channelList)))
72 raise ValueError("Channel %d is not in %s" %(channel, str(dataOut.channelList)))
77
73
78 index = dataOut.channelList.index(channel)
74 index = dataOut.channelList.index(channel)
79 channelIndexList.append(index)
75 channelIndexList.append(index)
80 dataOut = self.selectChannelsByIndex(dataOut,channelIndexList)
76 dataOut = self.selectChannelsByIndex(dataOut,channelIndexList)
81 return dataOut
77 return dataOut
82
78
83 def selectChannelsByIndex(self, dataOut, channelIndexList):
79 def selectChannelsByIndex(self, dataOut, channelIndexList):
84 """
80 """
85 Selecciona un bloque de datos en base a canales segun el channelIndexList
81 Selecciona un bloque de datos en base a canales segun el channelIndexList
86
82
87 Input:
83 Input:
88 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
84 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
89
85
90 Affected:
86 Affected:
91 dataOut.data
87 dataOut.data
92 dataOut.channelIndexList
88 dataOut.channelIndexList
93 dataOut.nChannels
89 dataOut.nChannels
94 dataOut.m_ProcessingHeader.totalSpectra
90 dataOut.m_ProcessingHeader.totalSpectra
95 dataOut.systemHeaderObj.numChannels
91 dataOut.systemHeaderObj.numChannels
96 dataOut.m_ProcessingHeader.blockSize
92 dataOut.m_ProcessingHeader.blockSize
97
93
98 Return:
94 Return:
99 None
95 None
100 """
96 """
101
97 #print("selectChannelsByIndex")
102 # for channelIndex in channelIndexList:
98 # for channelIndex in channelIndexList:
103 # if channelIndex not in dataOut.channelIndexList:
99 # if channelIndex not in dataOut.channelIndexList:
104 # raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
100 # raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
105
101
106 if dataOut.type == 'Voltage':
102 if dataOut.type == 'Voltage':
107 if dataOut.flagDataAsBlock:
103 if dataOut.flagDataAsBlock:
108 """
104 """
109 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
105 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
110 """
106 """
111 data = dataOut.data[channelIndexList,:,:]
107 data = dataOut.data[channelIndexList,:,:]
112 else:
108 else:
113 data = dataOut.data[channelIndexList,:]
109 data = dataOut.data[channelIndexList,:]
114
110
115 dataOut.data = data
111 dataOut.data = data
116 # dataOut.channelList = [dataOut.channelList[i] for i in channelIndexList]
112 # dataOut.channelList = [dataOut.channelList[i] for i in channelIndexList]
117 dataOut.channelList = range(len(channelIndexList))
113 dataOut.channelList = range(len(channelIndexList))
118
114
119 elif dataOut.type == 'Spectra':
115 elif dataOut.type == 'Spectra':
120 data_spc = dataOut.data_spc[channelIndexList, :]
116 if hasattr(dataOut, 'data_spc'):
121 dataOut.data_spc = data_spc
117 if dataOut.data_spc is None:
122 if hasattr(dataOut, 'data_dc') and dataOut.data_dc != None:
118 raise ValueError("data_spc is None")
123 data_dc = dataOut.data_dc[channelIndexList, :]
119 return dataOut
124 dataOut.data_dc = data_dc
120 else:
121 data_spc = dataOut.data_spc[channelIndexList, :]
122 dataOut.data_spc = data_spc
123
124 # if hasattr(dataOut, 'data_dc') :# and
125 # if dataOut.data_dc is None:
126 # raise ValueError("data_dc is None")
127 # return dataOut
128 # else:
129 # data_dc = dataOut.data_dc[channelIndexList, :]
130 # dataOut.data_dc = data_dc
125 # dataOut.channelList = [dataOut.channelList[i] for i in channelIndexList]
131 # dataOut.channelList = [dataOut.channelList[i] for i in channelIndexList]
126 dataOut.channelList = channelIndexList
132 dataOut.channelList = channelIndexList
127 dataOut = self.__selectPairsByChannel(dataOut,channelIndexList)
133 dataOut = self.__selectPairsByChannel(dataOut,channelIndexList)
128
134
129 return dataOut
135 return dataOut
130
136
131 def __selectPairsByChannel(self, dataOut, channelList=None):
137 def __selectPairsByChannel(self, dataOut, channelList=None):
132
138 #print("__selectPairsByChannel")
133 if channelList == None:
139 if channelList == None:
134 return
140 return
135
141
136 pairsIndexListSelected = []
142 pairsIndexListSelected = []
137 for pairIndex in dataOut.pairsIndexList:
143 for pairIndex in dataOut.pairsIndexList:
138 # First pair
144 # First pair
139 if dataOut.pairsList[pairIndex][0] not in channelList:
145 if dataOut.pairsList[pairIndex][0] not in channelList:
140 continue
146 continue
141 # Second pair
147 # Second pair
142 if dataOut.pairsList[pairIndex][1] not in channelList:
148 if dataOut.pairsList[pairIndex][1] not in channelList:
143 continue
149 continue
144
150
145 pairsIndexListSelected.append(pairIndex)
151 pairsIndexListSelected.append(pairIndex)
146
147 if not pairsIndexListSelected:
152 if not pairsIndexListSelected:
148 dataOut.data_cspc = None
153 dataOut.data_cspc = None
149 dataOut.pairsList = []
154 dataOut.pairsList = []
150 return
155 return
151
156
152 dataOut.data_cspc = dataOut.data_cspc[pairsIndexListSelected]
157 dataOut.data_cspc = dataOut.data_cspc[pairsIndexListSelected]
153 dataOut.pairsList = [dataOut.pairsList[i]
158 dataOut.pairsList = [dataOut.pairsList[i]
154 for i in pairsIndexListSelected]
159 for i in pairsIndexListSelected]
155
160
156 return dataOut
161 return dataOut
157
162
158 class selectHeights(Operation):
163 class selectHeights(Operation):
159
164
160 def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
165 def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
161 """
166 """
162 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
167 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
163 minHei <= height <= maxHei
168 minHei <= height <= maxHei
164
169
165 Input:
170 Input:
166 minHei : valor minimo de altura a considerar
171 minHei : valor minimo de altura a considerar
167 maxHei : valor maximo de altura a considerar
172 maxHei : valor maximo de altura a considerar
168
173
169 Affected:
174 Affected:
170 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
175 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
171
176
172 Return:
177 Return:
173 1 si el metodo se ejecuto con exito caso contrario devuelve 0
178 1 si el metodo se ejecuto con exito caso contrario devuelve 0
174 """
179 """
175
180
176 dataOut = dataOut
181 dataOut = dataOut
177
182
178 if minHei and maxHei:
183 if minHei and maxHei:
179
184
180 if (minHei < dataOut.heightList[0]):
185 if (minHei < dataOut.heightList[0]):
181 minHei = dataOut.heightList[0]
186 minHei = dataOut.heightList[0]
182
187
183 if (maxHei > dataOut.heightList[-1]):
188 if (maxHei > dataOut.heightList[-1]):
184 maxHei = dataOut.heightList[-1]
189 maxHei = dataOut.heightList[-1]
185
190
186 minIndex = 0
191 minIndex = 0
187 maxIndex = 0
192 maxIndex = 0
188 heights = dataOut.heightList
193 heights = dataOut.heightList
189
194
190 inda = numpy.where(heights >= minHei)
195 inda = numpy.where(heights >= minHei)
191 indb = numpy.where(heights <= maxHei)
196 indb = numpy.where(heights <= maxHei)
192
197
193 try:
198 try:
194 minIndex = inda[0][0]
199 minIndex = inda[0][0]
195 except:
200 except:
196 minIndex = 0
201 minIndex = 0
197
202
198 try:
203 try:
199 maxIndex = indb[0][-1]
204 maxIndex = indb[0][-1]
200 except:
205 except:
201 maxIndex = len(heights)
206 maxIndex = len(heights)
202
207
203 self.selectHeightsByIndex(minIndex, maxIndex)
208 self.selectHeightsByIndex(minIndex, maxIndex)
204
209
205 return self.dataOut
210 return self.dataOut
206
211
207 def selectHeightsByIndex(self, minIndex, maxIndex):
212 def selectHeightsByIndex(self, minIndex, maxIndex):
208 """
213 """
209 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
214 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
210 minIndex <= index <= maxIndex
215 minIndex <= index <= maxIndex
211
216
212 Input:
217 Input:
213 minIndex : valor de indice minimo de altura a considerar
218 minIndex : valor de indice minimo de altura a considerar
214 maxIndex : valor de indice maximo de altura a considerar
219 maxIndex : valor de indice maximo de altura a considerar
215
220
216 Affected:
221 Affected:
217 self.dataOut.data
222 self.dataOut.data
218 self.dataOut.heightList
223 self.dataOut.heightList
219
224
220 Return:
225 Return:
221 1 si el metodo se ejecuto con exito caso contrario devuelve 0
226 1 si el metodo se ejecuto con exito caso contrario devuelve 0
222 """
227 """
223
228
224 if self.dataOut.type == 'Voltage':
229 if self.dataOut.type == 'Voltage':
225 if (minIndex < 0) or (minIndex > maxIndex):
230 if (minIndex < 0) or (minIndex > maxIndex):
226 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
231 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
227
232
228 if (maxIndex >= self.dataOut.nHeights):
233 if (maxIndex >= self.dataOut.nHeights):
229 maxIndex = self.dataOut.nHeights
234 maxIndex = self.dataOut.nHeights
230
235
231 #voltage
236 #voltage
232 if self.dataOut.flagDataAsBlock:
237 if self.dataOut.flagDataAsBlock:
233 """
238 """
234 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
239 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
235 """
240 """
236 data = self.dataOut.data[:,:, minIndex:maxIndex]
241 data = self.dataOut.data[:,:, minIndex:maxIndex]
237 else:
242 else:
238 data = self.dataOut.data[:, minIndex:maxIndex]
243 data = self.dataOut.data[:, minIndex:maxIndex]
239
244
240 # firstHeight = self.dataOut.heightList[minIndex]
245 # firstHeight = self.dataOut.heightList[minIndex]
241
246
242 self.dataOut.data = data
247 self.dataOut.data = data
243 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
248 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
244
249
245 if self.dataOut.nHeights <= 1:
250 if self.dataOut.nHeights <= 1:
246 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
251 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
247 elif self.dataOut.type == 'Spectra':
252 elif self.dataOut.type == 'Spectra':
248 if (minIndex < 0) or (minIndex > maxIndex):
253 if (minIndex < 0) or (minIndex > maxIndex):
249 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
254 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
250 minIndex, maxIndex))
255 minIndex, maxIndex))
251
256
252 if (maxIndex >= self.dataOut.nHeights):
257 if (maxIndex >= self.dataOut.nHeights):
253 maxIndex = self.dataOut.nHeights - 1
258 maxIndex = self.dataOut.nHeights - 1
254
259
255 # Spectra
260 # Spectra
256 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
261 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
257
262
258 data_cspc = None
263 data_cspc = None
259 if self.dataOut.data_cspc is not None:
264 if self.dataOut.data_cspc is not None:
260 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
265 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
261
266
262 data_dc = None
267 data_dc = None
263 if self.dataOut.data_dc is not None:
268 if self.dataOut.data_dc is not None:
264 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
269 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
265
270
266 self.dataOut.data_spc = data_spc
271 self.dataOut.data_spc = data_spc
267 self.dataOut.data_cspc = data_cspc
272 self.dataOut.data_cspc = data_cspc
268 self.dataOut.data_dc = data_dc
273 self.dataOut.data_dc = data_dc
269
274
270 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
275 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
271
276
272 return 1
277 return 1
273
278
274
279
275 class filterByHeights(Operation):
280 class filterByHeights(Operation):
276
281
277 def run(self, dataOut, window):
282 def run(self, dataOut, window):
278
283
279 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
284 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
280
285
281 if window == None:
286 if window == None:
282 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
287 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
283
288
284 newdelta = deltaHeight * window
289 newdelta = deltaHeight * window
285 r = dataOut.nHeights % window
290 r = dataOut.nHeights % window
286 newheights = (dataOut.nHeights-r)/window
291 newheights = (dataOut.nHeights-r)/window
287
292
288 if newheights <= 1:
293 if newheights <= 1:
289 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
294 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
290
295
291 if dataOut.flagDataAsBlock:
296 if dataOut.flagDataAsBlock:
292 """
297 """
293 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
298 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
294 """
299 """
295 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
300 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
296 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
301 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
297 buffer = numpy.sum(buffer,3)
302 buffer = numpy.sum(buffer,3)
298
303
299 else:
304 else:
300 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
305 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
301 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
306 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
302 buffer = numpy.sum(buffer,2)
307 buffer = numpy.sum(buffer,2)
303
308
304 dataOut.data = buffer
309 dataOut.data = buffer
305 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
310 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
306 dataOut.windowOfFilter = window
311 dataOut.windowOfFilter = window
307
312
308 return dataOut
313 return dataOut
309
314
310
315
311 class setH0(Operation):
316 class setH0(Operation):
312
317
313 def run(self, dataOut, h0, deltaHeight = None):
318 def run(self, dataOut, h0, deltaHeight = None):
314
319
315 if not deltaHeight:
320 if not deltaHeight:
316 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
321 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
317
322
318 nHeights = dataOut.nHeights
323 nHeights = dataOut.nHeights
319
324
320 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
325 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
321
326
322 dataOut.heightList = newHeiRange
327 dataOut.heightList = newHeiRange
323
328
324 return dataOut
329 return dataOut
325
330
326
331
327 class deFlip(Operation):
332 class deFlip(Operation):
328
333
329 def run(self, dataOut, channelList = []):
334 def run(self, dataOut, channelList = []):
330
335
331 data = dataOut.data.copy()
336 data = dataOut.data.copy()
332
337
333 if dataOut.flagDataAsBlock:
338 if dataOut.flagDataAsBlock:
334 flip = self.flip
339 flip = self.flip
335 profileList = list(range(dataOut.nProfiles))
340 profileList = list(range(dataOut.nProfiles))
336
341
337 if not channelList:
342 if not channelList:
338 for thisProfile in profileList:
343 for thisProfile in profileList:
339 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
344 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
340 flip *= -1.0
345 flip *= -1.0
341 else:
346 else:
342 for thisChannel in channelList:
347 for thisChannel in channelList:
343 if thisChannel not in dataOut.channelList:
348 if thisChannel not in dataOut.channelList:
344 continue
349 continue
345
350
346 for thisProfile in profileList:
351 for thisProfile in profileList:
347 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
352 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
348 flip *= -1.0
353 flip *= -1.0
349
354
350 self.flip = flip
355 self.flip = flip
351
356
352 else:
357 else:
353 if not channelList:
358 if not channelList:
354 data[:,:] = data[:,:]*self.flip
359 data[:,:] = data[:,:]*self.flip
355 else:
360 else:
356 for thisChannel in channelList:
361 for thisChannel in channelList:
357 if thisChannel not in dataOut.channelList:
362 if thisChannel not in dataOut.channelList:
358 continue
363 continue
359
364
360 data[thisChannel,:] = data[thisChannel,:]*self.flip
365 data[thisChannel,:] = data[thisChannel,:]*self.flip
361
366
362 self.flip *= -1.
367 self.flip *= -1.
363
368
364 dataOut.data = data
369 dataOut.data = data
365
370
366 return dataOut
371 return dataOut
367
372
368
373
369 class setAttribute(Operation):
374 class setAttribute(Operation):
370 '''
375 '''
371 Set an arbitrary attribute(s) to dataOut
376 Set an arbitrary attribute(s) to dataOut
372 '''
377 '''
373
378
374 def __init__(self):
379 def __init__(self):
375
380
376 Operation.__init__(self)
381 Operation.__init__(self)
377 self._ready = False
382 self._ready = False
378
383
379 def run(self, dataOut, **kwargs):
384 def run(self, dataOut, **kwargs):
380
385
381 for key, value in kwargs.items():
386 for key, value in kwargs.items():
382 setattr(dataOut, key, value)
387 setattr(dataOut, key, value)
383
388
384 return dataOut
389 return dataOut
385
390
386
391
387 @MPDecorator
392 @MPDecorator
388 class printAttribute(Operation):
393 class printAttribute(Operation):
389 '''
394 '''
390 Print an arbitrary attribute of dataOut
395 Print an arbitrary attribute of dataOut
391 '''
396 '''
392
397
393 def __init__(self):
398 def __init__(self):
394
399
395 Operation.__init__(self)
400 Operation.__init__(self)
396
401
397 def run(self, dataOut, attributes):
402 def run(self, dataOut, attributes):
398
403
399 if isinstance(attributes, str):
404 if isinstance(attributes, str):
400 attributes = [attributes]
405 attributes = [attributes]
401 for attr in attributes:
406 for attr in attributes:
402 if hasattr(dataOut, attr):
407 if hasattr(dataOut, attr):
403 log.log(getattr(dataOut, attr), attr)
408 log.log(getattr(dataOut, attr), attr)
404
409
405
410
406 class interpolateHeights(Operation):
411 class interpolateHeights(Operation):
407
412
408 def run(self, dataOut, topLim, botLim):
413 def run(self, dataOut, topLim, botLim):
409 #69 al 72 para julia
414 #69 al 72 para julia
410 #82-84 para meteoros
415 #82-84 para meteoros
411 if len(numpy.shape(dataOut.data))==2:
416 if len(numpy.shape(dataOut.data))==2:
412 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
417 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
413 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
418 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
414 #dataOut.data[:,botLim:limSup+1] = sampInterp
419 #dataOut.data[:,botLim:limSup+1] = sampInterp
415 dataOut.data[:,botLim:topLim+1] = sampInterp
420 dataOut.data[:,botLim:topLim+1] = sampInterp
416 else:
421 else:
417 nHeights = dataOut.data.shape[2]
422 nHeights = dataOut.data.shape[2]
418 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
423 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
419 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
424 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
420 f = interpolate.interp1d(x, y, axis = 2)
425 f = interpolate.interp1d(x, y, axis = 2)
421 xnew = numpy.arange(botLim,topLim+1)
426 xnew = numpy.arange(botLim,topLim+1)
422 ynew = f(xnew)
427 ynew = f(xnew)
423 dataOut.data[:,:,botLim:topLim+1] = ynew
428 dataOut.data[:,:,botLim:topLim+1] = ynew
424
429
425 return dataOut
430 return dataOut
426
431
427
432
428 class CohInt(Operation):
433 class CohInt(Operation):
429
434
430 isConfig = False
435 isConfig = False
431 __profIndex = 0
436 __profIndex = 0
432 __byTime = False
437 __byTime = False
433 __initime = None
438 __initime = None
434 __lastdatatime = None
439 __lastdatatime = None
435 __integrationtime = None
440 __integrationtime = None
436 __buffer = None
441 __buffer = None
437 __bufferStride = []
442 __bufferStride = []
438 __dataReady = False
443 __dataReady = False
439 __profIndexStride = 0
444 __profIndexStride = 0
440 __dataToPutStride = False
445 __dataToPutStride = False
441 n = None
446 n = None
442
447
443 def __init__(self, **kwargs):
448 def __init__(self, **kwargs):
444
449
445 Operation.__init__(self, **kwargs)
450 Operation.__init__(self, **kwargs)
446
451
447 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
452 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
448 """
453 """
449 Set the parameters of the integration class.
454 Set the parameters of the integration class.
450
455
451 Inputs:
456 Inputs:
452
457
453 n : Number of coherent integrations
458 n : Number of coherent integrations
454 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
459 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
455 overlapping :
460 overlapping :
456 """
461 """
457
462
458 self.__initime = None
463 self.__initime = None
459 self.__lastdatatime = 0
464 self.__lastdatatime = 0
460 self.__buffer = None
465 self.__buffer = None
461 self.__dataReady = False
466 self.__dataReady = False
462 self.byblock = byblock
467 self.byblock = byblock
463 self.stride = stride
468 self.stride = stride
464
469
465 if n == None and timeInterval == None:
470 if n == None and timeInterval == None:
466 raise ValueError("n or timeInterval should be specified ...")
471 raise ValueError("n or timeInterval should be specified ...")
467
472
468 if n != None:
473 if n != None:
469 self.n = n
474 self.n = n
470 self.__byTime = False
475 self.__byTime = False
471 else:
476 else:
472 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
477 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
473 self.n = 9999
478 self.n = 9999
474 self.__byTime = True
479 self.__byTime = True
475
480
476 if overlapping:
481 if overlapping:
477 self.__withOverlapping = True
482 self.__withOverlapping = True
478 self.__buffer = None
483 self.__buffer = None
479 else:
484 else:
480 self.__withOverlapping = False
485 self.__withOverlapping = False
481 self.__buffer = 0
486 self.__buffer = 0
482
487
483 self.__profIndex = 0
488 self.__profIndex = 0
484
489
485 def putData(self, data):
490 def putData(self, data):
486
491
487 """
492 """
488 Add a profile to the __buffer and increase in one the __profileIndex
493 Add a profile to the __buffer and increase in one the __profileIndex
489
494
490 """
495 """
491
496
492 if not self.__withOverlapping:
497 if not self.__withOverlapping:
493 self.__buffer += data.copy()
498 self.__buffer += data.copy()
494 self.__profIndex += 1
499 self.__profIndex += 1
495 return
500 return
496
501
497 #Overlapping data
502 #Overlapping data
498 nChannels, nHeis = data.shape
503 nChannels, nHeis = data.shape
499 data = numpy.reshape(data, (1, nChannels, nHeis))
504 data = numpy.reshape(data, (1, nChannels, nHeis))
500
505
501 #If the buffer is empty then it takes the data value
506 #If the buffer is empty then it takes the data value
502 if self.__buffer is None:
507 if self.__buffer is None:
503 self.__buffer = data
508 self.__buffer = data
504 self.__profIndex += 1
509 self.__profIndex += 1
505 return
510 return
506
511
507 #If the buffer length is lower than n then stakcing the data value
512 #If the buffer length is lower than n then stakcing the data value
508 if self.__profIndex < self.n:
513 if self.__profIndex < self.n:
509 self.__buffer = numpy.vstack((self.__buffer, data))
514 self.__buffer = numpy.vstack((self.__buffer, data))
510 self.__profIndex += 1
515 self.__profIndex += 1
511 return
516 return
512
517
513 #If the buffer length is equal to n then replacing the last buffer value with the data value
518 #If the buffer length is equal to n then replacing the last buffer value with the data value
514 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
519 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
515 self.__buffer[self.n-1] = data
520 self.__buffer[self.n-1] = data
516 self.__profIndex = self.n
521 self.__profIndex = self.n
517 return
522 return
518
523
519
524
520 def pushData(self):
525 def pushData(self):
521 """
526 """
522 Return the sum of the last profiles and the profiles used in the sum.
527 Return the sum of the last profiles and the profiles used in the sum.
523
528
524 Affected:
529 Affected:
525
530
526 self.__profileIndex
531 self.__profileIndex
527
532
528 """
533 """
529
534
530 if not self.__withOverlapping:
535 if not self.__withOverlapping:
531 data = self.__buffer
536 data = self.__buffer
532 n = self.__profIndex
537 n = self.__profIndex
533
538
534 self.__buffer = 0
539 self.__buffer = 0
535 self.__profIndex = 0
540 self.__profIndex = 0
536
541
537 return data, n
542 return data, n
538
543
539 #Integration with Overlapping
544 #Integration with Overlapping
540 data = numpy.sum(self.__buffer, axis=0)
545 data = numpy.sum(self.__buffer, axis=0)
541 # print data
546 # print data
542 # raise
547 # raise
543 n = self.__profIndex
548 n = self.__profIndex
544
549
545 return data, n
550 return data, n
546
551
547 def byProfiles(self, data):
552 def byProfiles(self, data):
548
553
549 self.__dataReady = False
554 self.__dataReady = False
550 avgdata = None
555 avgdata = None
551 # n = None
556 # n = None
552 # print data
557 # print data
553 # raise
558 # raise
554 self.putData(data)
559 self.putData(data)
555
560
556 if self.__profIndex == self.n:
561 if self.__profIndex == self.n:
557 avgdata, n = self.pushData()
562 avgdata, n = self.pushData()
558 self.__dataReady = True
563 self.__dataReady = True
559
564
560 return avgdata
565 return avgdata
561
566
562 def byTime(self, data, datatime):
567 def byTime(self, data, datatime):
563
568
564 self.__dataReady = False
569 self.__dataReady = False
565 avgdata = None
570 avgdata = None
566 n = None
571 n = None
567
572
568 self.putData(data)
573 self.putData(data)
569
574
570 if (datatime - self.__initime) >= self.__integrationtime:
575 if (datatime - self.__initime) >= self.__integrationtime:
571 avgdata, n = self.pushData()
576 avgdata, n = self.pushData()
572 self.n = n
577 self.n = n
573 self.__dataReady = True
578 self.__dataReady = True
574
579
575 return avgdata
580 return avgdata
576
581
577 def integrateByStride(self, data, datatime):
582 def integrateByStride(self, data, datatime):
578 # print data
583 # print data
579 if self.__profIndex == 0:
584 if self.__profIndex == 0:
580 self.__buffer = [[data.copy(), datatime]]
585 self.__buffer = [[data.copy(), datatime]]
581 else:
586 else:
582 self.__buffer.append([data.copy(),datatime])
587 self.__buffer.append([data.copy(),datatime])
583 self.__profIndex += 1
588 self.__profIndex += 1
584 self.__dataReady = False
589 self.__dataReady = False
585
590
586 if self.__profIndex == self.n * self.stride :
591 if self.__profIndex == self.n * self.stride :
587 self.__dataToPutStride = True
592 self.__dataToPutStride = True
588 self.__profIndexStride = 0
593 self.__profIndexStride = 0
589 self.__profIndex = 0
594 self.__profIndex = 0
590 self.__bufferStride = []
595 self.__bufferStride = []
591 for i in range(self.stride):
596 for i in range(self.stride):
592 current = self.__buffer[i::self.stride]
597 current = self.__buffer[i::self.stride]
593 data = numpy.sum([t[0] for t in current], axis=0)
598 data = numpy.sum([t[0] for t in current], axis=0)
594 avgdatatime = numpy.average([t[1] for t in current])
599 avgdatatime = numpy.average([t[1] for t in current])
595 # print data
600 # print data
596 self.__bufferStride.append((data, avgdatatime))
601 self.__bufferStride.append((data, avgdatatime))
597
602
598 if self.__dataToPutStride:
603 if self.__dataToPutStride:
599 self.__dataReady = True
604 self.__dataReady = True
600 self.__profIndexStride += 1
605 self.__profIndexStride += 1
601 if self.__profIndexStride == self.stride:
606 if self.__profIndexStride == self.stride:
602 self.__dataToPutStride = False
607 self.__dataToPutStride = False
603 # print self.__bufferStride[self.__profIndexStride - 1]
608 # print self.__bufferStride[self.__profIndexStride - 1]
604 # raise
609 # raise
605 return self.__bufferStride[self.__profIndexStride - 1]
610 return self.__bufferStride[self.__profIndexStride - 1]
606
611
607
612
608 return None, None
613 return None, None
609
614
610 def integrate(self, data, datatime=None):
615 def integrate(self, data, datatime=None):
611
616
612 if self.__initime == None:
617 if self.__initime == None:
613 self.__initime = datatime
618 self.__initime = datatime
614
619
615 if self.__byTime:
620 if self.__byTime:
616 avgdata = self.byTime(data, datatime)
621 avgdata = self.byTime(data, datatime)
617 else:
622 else:
618 avgdata = self.byProfiles(data)
623 avgdata = self.byProfiles(data)
619
624
620
625
621 self.__lastdatatime = datatime
626 self.__lastdatatime = datatime
622
627
623 if avgdata is None:
628 if avgdata is None:
624 return None, None
629 return None, None
625
630
626 avgdatatime = self.__initime
631 avgdatatime = self.__initime
627
632
628 deltatime = datatime - self.__lastdatatime
633 deltatime = datatime - self.__lastdatatime
629
634
630 if not self.__withOverlapping:
635 if not self.__withOverlapping:
631 self.__initime = datatime
636 self.__initime = datatime
632 else:
637 else:
633 self.__initime += deltatime
638 self.__initime += deltatime
634
639
635 return avgdata, avgdatatime
640 return avgdata, avgdatatime
636
641
637 def integrateByBlock(self, dataOut):
642 def integrateByBlock(self, dataOut):
638
643
639 times = int(dataOut.data.shape[1]/self.n)
644 times = int(dataOut.data.shape[1]/self.n)
640 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
645 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
641
646
642 id_min = 0
647 id_min = 0
643 id_max = self.n
648 id_max = self.n
644
649
645 for i in range(times):
650 for i in range(times):
646 junk = dataOut.data[:,id_min:id_max,:]
651 junk = dataOut.data[:,id_min:id_max,:]
647 avgdata[:,i,:] = junk.sum(axis=1)
652 avgdata[:,i,:] = junk.sum(axis=1)
648 id_min += self.n
653 id_min += self.n
649 id_max += self.n
654 id_max += self.n
650
655
651 timeInterval = dataOut.ippSeconds*self.n
656 timeInterval = dataOut.ippSeconds*self.n
652 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
657 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
653 self.__dataReady = True
658 self.__dataReady = True
654 return avgdata, avgdatatime
659 return avgdata, avgdatatime
655
660
656 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
661 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
657
662
658 if not self.isConfig:
663 if not self.isConfig:
659 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
664 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
660 self.isConfig = True
665 self.isConfig = True
661
666
662 if dataOut.flagDataAsBlock:
667 if dataOut.flagDataAsBlock:
663 """
668 """
664 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
669 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
665 """
670 """
666 avgdata, avgdatatime = self.integrateByBlock(dataOut)
671 avgdata, avgdatatime = self.integrateByBlock(dataOut)
667 dataOut.nProfiles /= self.n
672 dataOut.nProfiles /= self.n
668 else:
673 else:
669 if stride is None:
674 if stride is None:
670 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
675 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
671 else:
676 else:
672 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
677 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
673
678
674
679
675 # dataOut.timeInterval *= n
680 # dataOut.timeInterval *= n
676 dataOut.flagNoData = True
681 dataOut.flagNoData = True
677
682
678 if self.__dataReady:
683 if self.__dataReady:
679 dataOut.data = avgdata
684 dataOut.data = avgdata
680 if not dataOut.flagCohInt:
685 if not dataOut.flagCohInt:
681 dataOut.nCohInt *= self.n
686 dataOut.nCohInt *= self.n
682 dataOut.flagCohInt = True
687 dataOut.flagCohInt = True
683 dataOut.utctime = avgdatatime
688 dataOut.utctime = avgdatatime
684 # print avgdata, avgdatatime
689 # print avgdata, avgdatatime
685 # raise
690 # raise
686 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
691 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
687 dataOut.flagNoData = False
692 dataOut.flagNoData = False
688 return dataOut
693 return dataOut
689
694
690 class Decoder(Operation):
695 class Decoder(Operation):
691
696
692 isConfig = False
697 isConfig = False
693 __profIndex = 0
698 __profIndex = 0
694
699
695 code = None
700 code = None
696
701
697 nCode = None
702 nCode = None
698 nBaud = None
703 nBaud = None
699
704
700 def __init__(self, **kwargs):
705 def __init__(self, **kwargs):
701
706
702 Operation.__init__(self, **kwargs)
707 Operation.__init__(self, **kwargs)
703
708
704 self.times = None
709 self.times = None
705 self.osamp = None
710 self.osamp = None
706 # self.__setValues = False
711 # self.__setValues = False
707 self.isConfig = False
712 self.isConfig = False
708 self.setupReq = False
713 self.setupReq = False
709 def setup(self, code, osamp, dataOut):
714 def setup(self, code, osamp, dataOut):
710
715
711 self.__profIndex = 0
716 self.__profIndex = 0
712
717
713 self.code = code
718 self.code = code
714
719
715 self.nCode = len(code)
720 self.nCode = len(code)
716 self.nBaud = len(code[0])
721 self.nBaud = len(code[0])
717 if (osamp != None) and (osamp >1):
722 if (osamp != None) and (osamp >1):
718 self.osamp = osamp
723 self.osamp = osamp
719 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
724 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
720 self.nBaud = self.nBaud*self.osamp
725 self.nBaud = self.nBaud*self.osamp
721
726
722 self.__nChannels = dataOut.nChannels
727 self.__nChannels = dataOut.nChannels
723 self.__nProfiles = dataOut.nProfiles
728 self.__nProfiles = dataOut.nProfiles
724 self.__nHeis = dataOut.nHeights
729 self.__nHeis = dataOut.nHeights
725
730
726 if self.__nHeis < self.nBaud:
731 if self.__nHeis < self.nBaud:
727 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
732 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
728
733
729 #Frequency
734 #Frequency
730 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
735 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
731
736
732 __codeBuffer[:,0:self.nBaud] = self.code
737 __codeBuffer[:,0:self.nBaud] = self.code
733
738
734 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
739 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
735
740
736 if dataOut.flagDataAsBlock:
741 if dataOut.flagDataAsBlock:
737
742
738 self.ndatadec = self.__nHeis #- self.nBaud + 1
743 self.ndatadec = self.__nHeis #- self.nBaud + 1
739
744
740 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
745 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
741
746
742 else:
747 else:
743
748
744 #Time
749 #Time
745 self.ndatadec = self.__nHeis #- self.nBaud + 1
750 self.ndatadec = self.__nHeis #- self.nBaud + 1
746
751
747 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
752 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
748
753
749 def __convolutionInFreq(self, data):
754 def __convolutionInFreq(self, data):
750
755
751 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
756 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
752
757
753 fft_data = numpy.fft.fft(data, axis=1)
758 fft_data = numpy.fft.fft(data, axis=1)
754
759
755 conv = fft_data*fft_code
760 conv = fft_data*fft_code
756
761
757 data = numpy.fft.ifft(conv,axis=1)
762 data = numpy.fft.ifft(conv,axis=1)
758
763
759 return data
764 return data
760
765
761 def __convolutionInFreqOpt(self, data):
766 def __convolutionInFreqOpt(self, data):
762
767
763 raise NotImplementedError
768 raise NotImplementedError
764
769
765 def __convolutionInTime(self, data):
770 def __convolutionInTime(self, data):
766
771
767 code = self.code[self.__profIndex]
772 code = self.code[self.__profIndex]
768 for i in range(self.__nChannels):
773 for i in range(self.__nChannels):
769 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
774 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
770
775
771 return self.datadecTime
776 return self.datadecTime
772
777
773 def __convolutionByBlockInTime(self, data):
778 def __convolutionByBlockInTime(self, data):
774
779
775 repetitions = int(self.__nProfiles / self.nCode)
780 repetitions = int(self.__nProfiles / self.nCode)
776 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
781 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
777 junk = junk.flatten()
782 junk = junk.flatten()
778 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
783 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
779 profilesList = range(self.__nProfiles)
784 profilesList = range(self.__nProfiles)
780
785
781 for i in range(self.__nChannels):
786 for i in range(self.__nChannels):
782 for j in profilesList:
787 for j in profilesList:
783 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
788 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
784 return self.datadecTime
789 return self.datadecTime
785
790
786 def __convolutionByBlockInFreq(self, data):
791 def __convolutionByBlockInFreq(self, data):
787
792
788 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
793 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
789
794
790
795
791 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
796 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
792
797
793 fft_data = numpy.fft.fft(data, axis=2)
798 fft_data = numpy.fft.fft(data, axis=2)
794
799
795 conv = fft_data*fft_code
800 conv = fft_data*fft_code
796
801
797 data = numpy.fft.ifft(conv,axis=2)
802 data = numpy.fft.ifft(conv,axis=2)
798
803
799 return data
804 return data
800
805
801
806
802 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
807 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
803
808
804 if dataOut.flagDecodeData:
809 if dataOut.flagDecodeData:
805 print("This data is already decoded, recoding again ...")
810 print("This data is already decoded, recoding again ...")
806
811
807 if not self.isConfig:
812 if not self.isConfig:
808
813
809 if code is None:
814 if code is None:
810 if dataOut.code is None:
815 if dataOut.code is None:
811 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
816 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
812
817
813 code = dataOut.code
818 code = dataOut.code
814 else:
819 else:
815 code = numpy.array(code).reshape(nCode,nBaud)
820 code = numpy.array(code).reshape(nCode,nBaud)
816 self.setup(code, osamp, dataOut)
821 self.setup(code, osamp, dataOut)
817
822
818 self.isConfig = True
823 self.isConfig = True
819
824
820 if mode == 3:
825 if mode == 3:
821 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
826 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
822
827
823 if times != None:
828 if times != None:
824 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
829 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
825
830
826 if self.code is None:
831 if self.code is None:
827 print("Fail decoding: Code is not defined.")
832 print("Fail decoding: Code is not defined.")
828 return
833 return
829
834
830 self.__nProfiles = dataOut.nProfiles
835 self.__nProfiles = dataOut.nProfiles
831 datadec = None
836 datadec = None
832
837
833 if mode == 3:
838 if mode == 3:
834 mode = 0
839 mode = 0
835
840
836 if dataOut.flagDataAsBlock:
841 if dataOut.flagDataAsBlock:
837 """
842 """
838 Decoding when data have been read as block,
843 Decoding when data have been read as block,
839 """
844 """
840
845
841 if mode == 0:
846 if mode == 0:
842 datadec = self.__convolutionByBlockInTime(dataOut.data)
847 datadec = self.__convolutionByBlockInTime(dataOut.data)
843 if mode == 1:
848 if mode == 1:
844 datadec = self.__convolutionByBlockInFreq(dataOut.data)
849 datadec = self.__convolutionByBlockInFreq(dataOut.data)
845 else:
850 else:
846 """
851 """
847 Decoding when data have been read profile by profile
852 Decoding when data have been read profile by profile
848 """
853 """
849 if mode == 0:
854 if mode == 0:
850 datadec = self.__convolutionInTime(dataOut.data)
855 datadec = self.__convolutionInTime(dataOut.data)
851
856
852 if mode == 1:
857 if mode == 1:
853 datadec = self.__convolutionInFreq(dataOut.data)
858 datadec = self.__convolutionInFreq(dataOut.data)
854
859
855 if mode == 2:
860 if mode == 2:
856 datadec = self.__convolutionInFreqOpt(dataOut.data)
861 datadec = self.__convolutionInFreqOpt(dataOut.data)
857
862
858 if datadec is None:
863 if datadec is None:
859 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
864 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
860
865
861 dataOut.code = self.code
866 dataOut.code = self.code
862 dataOut.nCode = self.nCode
867 dataOut.nCode = self.nCode
863 dataOut.nBaud = self.nBaud
868 dataOut.nBaud = self.nBaud
864
869
865 dataOut.data = datadec
870 dataOut.data = datadec
866
871
867 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
872 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
868
873
869 dataOut.flagDecodeData = True #asumo q la data esta decodificada
874 dataOut.flagDecodeData = True #asumo q la data esta decodificada
870
875
871 if self.__profIndex == self.nCode-1:
876 if self.__profIndex == self.nCode-1:
872 self.__profIndex = 0
877 self.__profIndex = 0
873 return dataOut
878 return dataOut
874
879
875 self.__profIndex += 1
880 self.__profIndex += 1
876
881
877 return dataOut
882 return dataOut
878 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
883 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
879
884
880
885
881 class ProfileConcat(Operation):
886 class ProfileConcat(Operation):
882
887
883 isConfig = False
888 isConfig = False
884 buffer = None
889 buffer = None
885
890
886 def __init__(self, **kwargs):
891 def __init__(self, **kwargs):
887
892
888 Operation.__init__(self, **kwargs)
893 Operation.__init__(self, **kwargs)
889 self.profileIndex = 0
894 self.profileIndex = 0
890
895
891 def reset(self):
896 def reset(self):
892 self.buffer = numpy.zeros_like(self.buffer)
897 self.buffer = numpy.zeros_like(self.buffer)
893 self.start_index = 0
898 self.start_index = 0
894 self.times = 1
899 self.times = 1
895
900
896 def setup(self, data, m, n=1):
901 def setup(self, data, m, n=1):
897 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
902 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
898 self.nHeights = data.shape[1]#.nHeights
903 self.nHeights = data.shape[1]#.nHeights
899 self.start_index = 0
904 self.start_index = 0
900 self.times = 1
905 self.times = 1
901
906
902 def concat(self, data):
907 def concat(self, data):
903
908
904 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
909 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
905 self.start_index = self.start_index + self.nHeights
910 self.start_index = self.start_index + self.nHeights
906
911
907 def run(self, dataOut, m):
912 def run(self, dataOut, m):
908 dataOut.flagNoData = True
913 dataOut.flagNoData = True
909
914
910 if not self.isConfig:
915 if not self.isConfig:
911 self.setup(dataOut.data, m, 1)
916 self.setup(dataOut.data, m, 1)
912 self.isConfig = True
917 self.isConfig = True
913
918
914 if dataOut.flagDataAsBlock:
919 if dataOut.flagDataAsBlock:
915 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
920 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
916
921
917 else:
922 else:
918 self.concat(dataOut.data)
923 self.concat(dataOut.data)
919 self.times += 1
924 self.times += 1
920 if self.times > m:
925 if self.times > m:
921 dataOut.data = self.buffer
926 dataOut.data = self.buffer
922 self.reset()
927 self.reset()
923 dataOut.flagNoData = False
928 dataOut.flagNoData = False
924 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
929 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
925 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
930 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
926 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
931 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
927 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
932 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
928 dataOut.ippSeconds *= m
933 dataOut.ippSeconds *= m
929 return dataOut
934 return dataOut
930
935
931 class ProfileSelector(Operation):
936 class ProfileSelector(Operation):
932
937
933 profileIndex = None
938 profileIndex = None
934 # Tamanho total de los perfiles
939 # Tamanho total de los perfiles
935 nProfiles = None
940 nProfiles = None
936
941
937 def __init__(self, **kwargs):
942 def __init__(self, **kwargs):
938
943
939 Operation.__init__(self, **kwargs)
944 Operation.__init__(self, **kwargs)
940 self.profileIndex = 0
945 self.profileIndex = 0
941
946
942 def incProfileIndex(self):
947 def incProfileIndex(self):
943
948
944 self.profileIndex += 1
949 self.profileIndex += 1
945
950
946 if self.profileIndex >= self.nProfiles:
951 if self.profileIndex >= self.nProfiles:
947 self.profileIndex = 0
952 self.profileIndex = 0
948
953
949 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
954 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
950
955
951 if profileIndex < minIndex:
956 if profileIndex < minIndex:
952 return False
957 return False
953
958
954 if profileIndex > maxIndex:
959 if profileIndex > maxIndex:
955 return False
960 return False
956
961
957 return True
962 return True
958
963
959 def isThisProfileInList(self, profileIndex, profileList):
964 def isThisProfileInList(self, profileIndex, profileList):
960
965
961 if profileIndex not in profileList:
966 if profileIndex not in profileList:
962 return False
967 return False
963
968
964 return True
969 return True
965
970
966 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
971 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
967
972
968 """
973 """
969 ProfileSelector:
974 ProfileSelector:
970
975
971 Inputs:
976 Inputs:
972 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
977 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
973
978
974 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
979 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
975
980
976 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
981 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
977
982
978 """
983 """
979
984
980 if rangeList is not None:
985 if rangeList is not None:
981 if type(rangeList[0]) not in (tuple, list):
986 if type(rangeList[0]) not in (tuple, list):
982 rangeList = [rangeList]
987 rangeList = [rangeList]
983
988
984 dataOut.flagNoData = True
989 dataOut.flagNoData = True
985
990
986 if dataOut.flagDataAsBlock:
991 if dataOut.flagDataAsBlock:
987 """
992 """
988 data dimension = [nChannels, nProfiles, nHeis]
993 data dimension = [nChannels, nProfiles, nHeis]
989 """
994 """
990 if profileList != None:
995 if profileList != None:
991 dataOut.data = dataOut.data[:,profileList,:]
996 dataOut.data = dataOut.data[:,profileList,:]
992
997
993 if profileRangeList != None:
998 if profileRangeList != None:
994 minIndex = profileRangeList[0]
999 minIndex = profileRangeList[0]
995 maxIndex = profileRangeList[1]
1000 maxIndex = profileRangeList[1]
996 profileList = list(range(minIndex, maxIndex+1))
1001 profileList = list(range(minIndex, maxIndex+1))
997
1002
998 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
1003 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
999
1004
1000 if rangeList != None:
1005 if rangeList != None:
1001
1006
1002 profileList = []
1007 profileList = []
1003
1008
1004 for thisRange in rangeList:
1009 for thisRange in rangeList:
1005 minIndex = thisRange[0]
1010 minIndex = thisRange[0]
1006 maxIndex = thisRange[1]
1011 maxIndex = thisRange[1]
1007
1012
1008 profileList.extend(list(range(minIndex, maxIndex+1)))
1013 profileList.extend(list(range(minIndex, maxIndex+1)))
1009
1014
1010 dataOut.data = dataOut.data[:,profileList,:]
1015 dataOut.data = dataOut.data[:,profileList,:]
1011
1016
1012 dataOut.nProfiles = len(profileList)
1017 dataOut.nProfiles = len(profileList)
1013 dataOut.profileIndex = dataOut.nProfiles - 1
1018 dataOut.profileIndex = dataOut.nProfiles - 1
1014 dataOut.flagNoData = False
1019 dataOut.flagNoData = False
1015
1020
1016 return dataOut
1021 return dataOut
1017
1022
1018 """
1023 """
1019 data dimension = [nChannels, nHeis]
1024 data dimension = [nChannels, nHeis]
1020 """
1025 """
1021
1026
1022 if profileList != None:
1027 if profileList != None:
1023
1028
1024 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1029 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1025
1030
1026 self.nProfiles = len(profileList)
1031 self.nProfiles = len(profileList)
1027 dataOut.nProfiles = self.nProfiles
1032 dataOut.nProfiles = self.nProfiles
1028 dataOut.profileIndex = self.profileIndex
1033 dataOut.profileIndex = self.profileIndex
1029 dataOut.flagNoData = False
1034 dataOut.flagNoData = False
1030
1035
1031 self.incProfileIndex()
1036 self.incProfileIndex()
1032 return dataOut
1037 return dataOut
1033
1038
1034 if profileRangeList != None:
1039 if profileRangeList != None:
1035
1040
1036 minIndex = profileRangeList[0]
1041 minIndex = profileRangeList[0]
1037 maxIndex = profileRangeList[1]
1042 maxIndex = profileRangeList[1]
1038
1043
1039 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1044 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1040
1045
1041 self.nProfiles = maxIndex - minIndex + 1
1046 self.nProfiles = maxIndex - minIndex + 1
1042 dataOut.nProfiles = self.nProfiles
1047 dataOut.nProfiles = self.nProfiles
1043 dataOut.profileIndex = self.profileIndex
1048 dataOut.profileIndex = self.profileIndex
1044 dataOut.flagNoData = False
1049 dataOut.flagNoData = False
1045
1050
1046 self.incProfileIndex()
1051 self.incProfileIndex()
1047 return dataOut
1052 return dataOut
1048
1053
1049 if rangeList != None:
1054 if rangeList != None:
1050
1055
1051 nProfiles = 0
1056 nProfiles = 0
1052
1057
1053 for thisRange in rangeList:
1058 for thisRange in rangeList:
1054 minIndex = thisRange[0]
1059 minIndex = thisRange[0]
1055 maxIndex = thisRange[1]
1060 maxIndex = thisRange[1]
1056
1061
1057 nProfiles += maxIndex - minIndex + 1
1062 nProfiles += maxIndex - minIndex + 1
1058
1063
1059 for thisRange in rangeList:
1064 for thisRange in rangeList:
1060
1065
1061 minIndex = thisRange[0]
1066 minIndex = thisRange[0]
1062 maxIndex = thisRange[1]
1067 maxIndex = thisRange[1]
1063
1068
1064 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1069 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1065
1070
1066 self.nProfiles = nProfiles
1071 self.nProfiles = nProfiles
1067 dataOut.nProfiles = self.nProfiles
1072 dataOut.nProfiles = self.nProfiles
1068 dataOut.profileIndex = self.profileIndex
1073 dataOut.profileIndex = self.profileIndex
1069 dataOut.flagNoData = False
1074 dataOut.flagNoData = False
1070
1075
1071 self.incProfileIndex()
1076 self.incProfileIndex()
1072
1077
1073 break
1078 break
1074
1079
1075 return dataOut
1080 return dataOut
1076
1081
1077
1082
1078 if beam != None: #beam is only for AMISR data
1083 if beam != None: #beam is only for AMISR data
1079 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1084 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1080 dataOut.flagNoData = False
1085 dataOut.flagNoData = False
1081 dataOut.profileIndex = self.profileIndex
1086 dataOut.profileIndex = self.profileIndex
1082
1087
1083 self.incProfileIndex()
1088 self.incProfileIndex()
1084
1089
1085 return dataOut
1090 return dataOut
1086
1091
1087 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1092 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1088
1093
1089
1094
1090 class Reshaper(Operation):
1095 class Reshaper(Operation):
1091
1096
1092 def __init__(self, **kwargs):
1097 def __init__(self, **kwargs):
1093
1098
1094 Operation.__init__(self, **kwargs)
1099 Operation.__init__(self, **kwargs)
1095
1100
1096 self.__buffer = None
1101 self.__buffer = None
1097 self.__nitems = 0
1102 self.__nitems = 0
1098
1103
1099 def __appendProfile(self, dataOut, nTxs):
1104 def __appendProfile(self, dataOut, nTxs):
1100
1105
1101 if self.__buffer is None:
1106 if self.__buffer is None:
1102 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1107 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1103 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1108 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1104
1109
1105 ini = dataOut.nHeights * self.__nitems
1110 ini = dataOut.nHeights * self.__nitems
1106 end = ini + dataOut.nHeights
1111 end = ini + dataOut.nHeights
1107
1112
1108 self.__buffer[:, ini:end] = dataOut.data
1113 self.__buffer[:, ini:end] = dataOut.data
1109
1114
1110 self.__nitems += 1
1115 self.__nitems += 1
1111
1116
1112 return int(self.__nitems*nTxs)
1117 return int(self.__nitems*nTxs)
1113
1118
1114 def __getBuffer(self):
1119 def __getBuffer(self):
1115
1120
1116 if self.__nitems == int(1./self.__nTxs):
1121 if self.__nitems == int(1./self.__nTxs):
1117
1122
1118 self.__nitems = 0
1123 self.__nitems = 0
1119
1124
1120 return self.__buffer.copy()
1125 return self.__buffer.copy()
1121
1126
1122 return None
1127 return None
1123
1128
1124 def __checkInputs(self, dataOut, shape, nTxs):
1129 def __checkInputs(self, dataOut, shape, nTxs):
1125
1130
1126 if shape is None and nTxs is None:
1131 if shape is None and nTxs is None:
1127 raise ValueError("Reshaper: shape of factor should be defined")
1132 raise ValueError("Reshaper: shape of factor should be defined")
1128
1133
1129 if nTxs:
1134 if nTxs:
1130 if nTxs < 0:
1135 if nTxs < 0:
1131 raise ValueError("nTxs should be greater than 0")
1136 raise ValueError("nTxs should be greater than 0")
1132
1137
1133 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1138 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1134 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1139 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1135
1140
1136 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1141 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1137
1142
1138 return shape, nTxs
1143 return shape, nTxs
1139
1144
1140 if len(shape) != 2 and len(shape) != 3:
1145 if len(shape) != 2 and len(shape) != 3:
1141 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1146 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1142
1147
1143 if len(shape) == 2:
1148 if len(shape) == 2:
1144 shape_tuple = [dataOut.nChannels]
1149 shape_tuple = [dataOut.nChannels]
1145 shape_tuple.extend(shape)
1150 shape_tuple.extend(shape)
1146 else:
1151 else:
1147 shape_tuple = list(shape)
1152 shape_tuple = list(shape)
1148
1153
1149 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1154 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1150
1155
1151 return shape_tuple, nTxs
1156 return shape_tuple, nTxs
1152
1157
1153 def run(self, dataOut, shape=None, nTxs=None):
1158 def run(self, dataOut, shape=None, nTxs=None):
1154
1159
1155 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1160 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1156
1161
1157 dataOut.flagNoData = True
1162 dataOut.flagNoData = True
1158 profileIndex = None
1163 profileIndex = None
1159
1164
1160 if dataOut.flagDataAsBlock:
1165 if dataOut.flagDataAsBlock:
1161
1166
1162 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1167 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1163 dataOut.flagNoData = False
1168 dataOut.flagNoData = False
1164
1169
1165 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1170 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1166
1171
1167 else:
1172 else:
1168
1173
1169 if self.__nTxs < 1:
1174 if self.__nTxs < 1:
1170
1175
1171 self.__appendProfile(dataOut, self.__nTxs)
1176 self.__appendProfile(dataOut, self.__nTxs)
1172 new_data = self.__getBuffer()
1177 new_data = self.__getBuffer()
1173
1178
1174 if new_data is not None:
1179 if new_data is not None:
1175 dataOut.data = new_data
1180 dataOut.data = new_data
1176 dataOut.flagNoData = False
1181 dataOut.flagNoData = False
1177
1182
1178 profileIndex = dataOut.profileIndex*nTxs
1183 profileIndex = dataOut.profileIndex*nTxs
1179
1184
1180 else:
1185 else:
1181 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1186 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1182
1187
1183 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1188 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1184
1189
1185 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1190 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1186
1191
1187 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1192 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1188
1193
1189 dataOut.profileIndex = profileIndex
1194 dataOut.profileIndex = profileIndex
1190
1195
1191 dataOut.ippSeconds /= self.__nTxs
1196 dataOut.ippSeconds /= self.__nTxs
1192
1197
1193 return dataOut
1198 return dataOut
1194
1199
1195 class SplitProfiles(Operation):
1200 class SplitProfiles(Operation):
1196
1201
1197 def __init__(self, **kwargs):
1202 def __init__(self, **kwargs):
1198
1203
1199 Operation.__init__(self, **kwargs)
1204 Operation.__init__(self, **kwargs)
1200
1205
1201 def run(self, dataOut, n):
1206 def run(self, dataOut, n):
1202
1207
1203 dataOut.flagNoData = True
1208 dataOut.flagNoData = True
1204 profileIndex = None
1209 profileIndex = None
1205
1210
1206 if dataOut.flagDataAsBlock:
1211 if dataOut.flagDataAsBlock:
1207
1212
1208 #nchannels, nprofiles, nsamples
1213 #nchannels, nprofiles, nsamples
1209 shape = dataOut.data.shape
1214 shape = dataOut.data.shape
1210
1215
1211 if shape[2] % n != 0:
1216 if shape[2] % n != 0:
1212 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1217 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1213
1218
1214 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1219 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1215
1220
1216 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1221 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1217 dataOut.flagNoData = False
1222 dataOut.flagNoData = False
1218
1223
1219 profileIndex = int(dataOut.nProfiles/n) - 1
1224 profileIndex = int(dataOut.nProfiles/n) - 1
1220
1225
1221 else:
1226 else:
1222
1227
1223 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1228 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1224
1229
1225 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1230 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1226
1231
1227 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1232 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1228
1233
1229 dataOut.nProfiles = int(dataOut.nProfiles*n)
1234 dataOut.nProfiles = int(dataOut.nProfiles*n)
1230
1235
1231 dataOut.profileIndex = profileIndex
1236 dataOut.profileIndex = profileIndex
1232
1237
1233 dataOut.ippSeconds /= n
1238 dataOut.ippSeconds /= n
1234
1239
1235 return dataOut
1240 return dataOut
1236
1241
1237 class CombineProfiles(Operation):
1242 class CombineProfiles(Operation):
1238 def __init__(self, **kwargs):
1243 def __init__(self, **kwargs):
1239
1244
1240 Operation.__init__(self, **kwargs)
1245 Operation.__init__(self, **kwargs)
1241
1246
1242 self.__remData = None
1247 self.__remData = None
1243 self.__profileIndex = 0
1248 self.__profileIndex = 0
1244
1249
1245 def run(self, dataOut, n):
1250 def run(self, dataOut, n):
1246
1251
1247 dataOut.flagNoData = True
1252 dataOut.flagNoData = True
1248 profileIndex = None
1253 profileIndex = None
1249
1254
1250 if dataOut.flagDataAsBlock:
1255 if dataOut.flagDataAsBlock:
1251
1256
1252 #nchannels, nprofiles, nsamples
1257 #nchannels, nprofiles, nsamples
1253 shape = dataOut.data.shape
1258 shape = dataOut.data.shape
1254 new_shape = shape[0], shape[1]/n, shape[2]*n
1259 new_shape = shape[0], shape[1]/n, shape[2]*n
1255
1260
1256 if shape[1] % n != 0:
1261 if shape[1] % n != 0:
1257 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1262 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1258
1263
1259 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1264 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1260 dataOut.flagNoData = False
1265 dataOut.flagNoData = False
1261
1266
1262 profileIndex = int(dataOut.nProfiles*n) - 1
1267 profileIndex = int(dataOut.nProfiles*n) - 1
1263
1268
1264 else:
1269 else:
1265
1270
1266 #nchannels, nsamples
1271 #nchannels, nsamples
1267 if self.__remData is None:
1272 if self.__remData is None:
1268 newData = dataOut.data
1273 newData = dataOut.data
1269 else:
1274 else:
1270 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1275 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1271
1276
1272 self.__profileIndex += 1
1277 self.__profileIndex += 1
1273
1278
1274 if self.__profileIndex < n:
1279 if self.__profileIndex < n:
1275 self.__remData = newData
1280 self.__remData = newData
1276 #continue
1281 #continue
1277 return
1282 return
1278
1283
1279 self.__profileIndex = 0
1284 self.__profileIndex = 0
1280 self.__remData = None
1285 self.__remData = None
1281
1286
1282 dataOut.data = newData
1287 dataOut.data = newData
1283 dataOut.flagNoData = False
1288 dataOut.flagNoData = False
1284
1289
1285 profileIndex = dataOut.profileIndex/n
1290 profileIndex = dataOut.profileIndex/n
1286
1291
1287
1292
1288 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1293 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1289
1294
1290 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1295 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1291
1296
1292 dataOut.nProfiles = int(dataOut.nProfiles/n)
1297 dataOut.nProfiles = int(dataOut.nProfiles/n)
1293
1298
1294 dataOut.profileIndex = profileIndex
1299 dataOut.profileIndex = profileIndex
1295
1300
1296 dataOut.ippSeconds *= n
1301 dataOut.ippSeconds *= n
1297
1302
1298 return dataOut
1303 return dataOut
1299
1304
1300 class PulsePairVoltage(Operation):
1305 class PulsePairVoltage(Operation):
1301 '''
1306 '''
1302 Function PulsePair(Signal Power, Velocity)
1307 Function PulsePair(Signal Power, Velocity)
1303 The real component of Lag[0] provides Intensity Information
1308 The real component of Lag[0] provides Intensity Information
1304 The imag component of Lag[1] Phase provides Velocity Information
1309 The imag component of Lag[1] Phase provides Velocity Information
1305
1310
1306 Configuration Parameters:
1311 Configuration Parameters:
1307 nPRF = Number of Several PRF
1312 nPRF = Number of Several PRF
1308 theta = Degree Azimuth angel Boundaries
1313 theta = Degree Azimuth angel Boundaries
1309
1314
1310 Input:
1315 Input:
1311 self.dataOut
1316 self.dataOut
1312 lag[N]
1317 lag[N]
1313 Affected:
1318 Affected:
1314 self.dataOut.spc
1319 self.dataOut.spc
1315 '''
1320 '''
1316 isConfig = False
1321 isConfig = False
1317 __profIndex = 0
1322 __profIndex = 0
1318 __initime = None
1323 __initime = None
1319 __lastdatatime = None
1324 __lastdatatime = None
1320 __buffer = None
1325 __buffer = None
1321 noise = None
1326 noise = None
1322 __dataReady = False
1327 __dataReady = False
1323 n = None
1328 n = None
1324 __nch = 0
1329 __nch = 0
1325 __nHeis = 0
1330 __nHeis = 0
1326 removeDC = False
1331 removeDC = False
1327 ipp = None
1332 ipp = None
1328 lambda_ = 0
1333 lambda_ = 0
1329
1334
1330 def __init__(self,**kwargs):
1335 def __init__(self,**kwargs):
1331 Operation.__init__(self,**kwargs)
1336 Operation.__init__(self,**kwargs)
1332
1337
1333 def setup(self, dataOut, n = None, removeDC=False):
1338 def setup(self, dataOut, n = None, removeDC=False):
1334 '''
1339 '''
1335 n= Numero de PRF's de entrada
1340 n= Numero de PRF's de entrada
1336 '''
1341 '''
1337 self.__initime = None
1342 self.__initime = None
1338 self.__lastdatatime = 0
1343 self.__lastdatatime = 0
1339 self.__dataReady = False
1344 self.__dataReady = False
1340 self.__buffer = 0
1345 self.__buffer = 0
1341 self.__profIndex = 0
1346 self.__profIndex = 0
1342 self.noise = None
1347 self.noise = None
1343 self.__nch = dataOut.nChannels
1348 self.__nch = dataOut.nChannels
1344 self.__nHeis = dataOut.nHeights
1349 self.__nHeis = dataOut.nHeights
1345 self.removeDC = removeDC
1350 self.removeDC = removeDC
1346 self.lambda_ = 3.0e8/(9345.0e6)
1351 self.lambda_ = 3.0e8/(9345.0e6)
1347 self.ippSec = dataOut.ippSeconds
1352 self.ippSec = dataOut.ippSeconds
1348 self.nCohInt = dataOut.nCohInt
1353 self.nCohInt = dataOut.nCohInt
1349
1354
1350 if n == None:
1355 if n == None:
1351 raise ValueError("n should be specified.")
1356 raise ValueError("n should be specified.")
1352
1357
1353 if n != None:
1358 if n != None:
1354 if n<2:
1359 if n<2:
1355 raise ValueError("n should be greater than 2")
1360 raise ValueError("n should be greater than 2")
1356
1361
1357 self.n = n
1362 self.n = n
1358 self.__nProf = n
1363 self.__nProf = n
1359
1364
1360 self.__buffer = numpy.zeros((dataOut.nChannels,
1365 self.__buffer = numpy.zeros((dataOut.nChannels,
1361 n,
1366 n,
1362 dataOut.nHeights),
1367 dataOut.nHeights),
1363 dtype='complex')
1368 dtype='complex')
1364
1369
1365 def putData(self,data):
1370 def putData(self,data):
1366 '''
1371 '''
1367 Add a profile to he __buffer and increase in one the __profiel Index
1372 Add a profile to he __buffer and increase in one the __profiel Index
1368 '''
1373 '''
1369 self.__buffer[:,self.__profIndex,:]= data
1374 self.__buffer[:,self.__profIndex,:]= data
1370 self.__profIndex += 1
1375 self.__profIndex += 1
1371 return
1376 return
1372
1377
1373 def pushData(self,dataOut):
1378 def pushData(self,dataOut):
1374 '''
1379 '''
1375 Return the PULSEPAIR and the profiles used in the operation
1380 Return the PULSEPAIR and the profiles used in the operation
1376 Affected : self.__profileIndex
1381 Affected : self.__profileIndex
1377 '''
1382 '''
1378 #----------------- Remove DC-----------------------------------
1383 #----------------- Remove DC-----------------------------------
1379 if self.removeDC==True:
1384 if self.removeDC==True:
1380 mean = numpy.mean(self.__buffer,1)
1385 mean = numpy.mean(self.__buffer,1)
1381 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1386 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1382 dc= numpy.tile(tmp,[1,self.__nProf,1])
1387 dc= numpy.tile(tmp,[1,self.__nProf,1])
1383 self.__buffer = self.__buffer - dc
1388 self.__buffer = self.__buffer - dc
1384 #------------------Calculo de Potencia ------------------------
1389 #------------------Calculo de Potencia ------------------------
1385 pair0 = self.__buffer*numpy.conj(self.__buffer)
1390 pair0 = self.__buffer*numpy.conj(self.__buffer)
1386 pair0 = pair0.real
1391 pair0 = pair0.real
1387 lag_0 = numpy.sum(pair0,1)
1392 lag_0 = numpy.sum(pair0,1)
1388 #------------------Calculo de Ruido x canal--------------------
1393 #------------------Calculo de Ruido x canal--------------------
1389 self.noise = numpy.zeros(self.__nch)
1394 self.noise = numpy.zeros(self.__nch)
1390 for i in range(self.__nch):
1395 for i in range(self.__nch):
1391 daux = numpy.sort(pair0[i,:,:],axis= None)
1396 daux = numpy.sort(pair0[i,:,:],axis= None)
1392 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1397 self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
1393
1398
1394 self.noise = self.noise.reshape(self.__nch,1)
1399 self.noise = self.noise.reshape(self.__nch,1)
1395 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1400 self.noise = numpy.tile(self.noise,[1,self.__nHeis])
1396 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1401 noise_buffer = self.noise.reshape(self.__nch,1,self.__nHeis)
1397 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1402 noise_buffer = numpy.tile(noise_buffer,[1,self.__nProf,1])
1398 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1403 #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
1399 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1404 #------------------ P= S+N ,P=lag_0/N ---------------------------------
1400 #-------------------- Power --------------------------------------------------
1405 #-------------------- Power --------------------------------------------------
1401 data_power = lag_0/(self.n*self.nCohInt)
1406 data_power = lag_0/(self.n*self.nCohInt)
1402 #------------------ Senal ---------------------------------------------------
1407 #------------------ Senal ---------------------------------------------------
1403 data_intensity = pair0 - noise_buffer
1408 data_intensity = pair0 - noise_buffer
1404 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1409 data_intensity = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
1405 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1410 #data_intensity = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
1406 for i in range(self.__nch):
1411 for i in range(self.__nch):
1407 for j in range(self.__nHeis):
1412 for j in range(self.__nHeis):
1408 if data_intensity[i][j] < 0:
1413 if data_intensity[i][j] < 0:
1409 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1414 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
1410
1415
1411 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1416 #----------------- Calculo de Frecuencia y Velocidad doppler--------
1412 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1417 pair1 = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
1413 lag_1 = numpy.sum(pair1,1)
1418 lag_1 = numpy.sum(pair1,1)
1414 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1419 data_freq = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
1415 data_velocity = (self.lambda_/2.0)*data_freq
1420 data_velocity = (self.lambda_/2.0)*data_freq
1416
1421
1417 #---------------- Potencia promedio estimada de la Senal-----------
1422 #---------------- Potencia promedio estimada de la Senal-----------
1418 lag_0 = lag_0/self.n
1423 lag_0 = lag_0/self.n
1419 S = lag_0-self.noise
1424 S = lag_0-self.noise
1420
1425
1421 #---------------- Frecuencia Doppler promedio ---------------------
1426 #---------------- Frecuencia Doppler promedio ---------------------
1422 lag_1 = lag_1/(self.n-1)
1427 lag_1 = lag_1/(self.n-1)
1423 R1 = numpy.abs(lag_1)
1428 R1 = numpy.abs(lag_1)
1424
1429
1425 #---------------- Calculo del SNR----------------------------------
1430 #---------------- Calculo del SNR----------------------------------
1426 data_snrPP = S/self.noise
1431 data_snrPP = S/self.noise
1427 for i in range(self.__nch):
1432 for i in range(self.__nch):
1428 for j in range(self.__nHeis):
1433 for j in range(self.__nHeis):
1429 if data_snrPP[i][j] < 1.e-20:
1434 if data_snrPP[i][j] < 1.e-20:
1430 data_snrPP[i][j] = 1.e-20
1435 data_snrPP[i][j] = 1.e-20
1431
1436
1432 #----------------- Calculo del ancho espectral ----------------------
1437 #----------------- Calculo del ancho espectral ----------------------
1433 L = S/R1
1438 L = S/R1
1434 L = numpy.where(L<0,1,L)
1439 L = numpy.where(L<0,1,L)
1435 L = numpy.log(L)
1440 L = numpy.log(L)
1436 tmp = numpy.sqrt(numpy.absolute(L))
1441 tmp = numpy.sqrt(numpy.absolute(L))
1437 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1442 data_specwidth = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
1438 n = self.__profIndex
1443 n = self.__profIndex
1439
1444
1440 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1445 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1441 self.__profIndex = 0
1446 self.__profIndex = 0
1442 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n
1447 return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n
1443
1448
1444
1449
1445 def pulsePairbyProfiles(self,dataOut):
1450 def pulsePairbyProfiles(self,dataOut):
1446
1451
1447 self.__dataReady = False
1452 self.__dataReady = False
1448 data_power = None
1453 data_power = None
1449 data_intensity = None
1454 data_intensity = None
1450 data_velocity = None
1455 data_velocity = None
1451 data_specwidth = None
1456 data_specwidth = None
1452 data_snrPP = None
1457 data_snrPP = None
1453 self.putData(data=dataOut.data)
1458 self.putData(data=dataOut.data)
1454 if self.__profIndex == self.n:
1459 if self.__profIndex == self.n:
1455 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut)
1460 data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n = self.pushData(dataOut=dataOut)
1456 self.__dataReady = True
1461 self.__dataReady = True
1457
1462
1458 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth
1463 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth
1459
1464
1460
1465
1461 def pulsePairOp(self, dataOut, datatime= None):
1466 def pulsePairOp(self, dataOut, datatime= None):
1462
1467
1463 if self.__initime == None:
1468 if self.__initime == None:
1464 self.__initime = datatime
1469 self.__initime = datatime
1465 data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
1470 data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
1466 self.__lastdatatime = datatime
1471 self.__lastdatatime = datatime
1467
1472
1468 if data_power is None:
1473 if data_power is None:
1469 return None, None, None,None,None,None
1474 return None, None, None,None,None,None
1470
1475
1471 avgdatatime = self.__initime
1476 avgdatatime = self.__initime
1472 deltatime = datatime - self.__lastdatatime
1477 deltatime = datatime - self.__lastdatatime
1473 self.__initime = datatime
1478 self.__initime = datatime
1474
1479
1475 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime
1480 return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime
1476
1481
1477 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1482 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1478
1483
1479 if not self.isConfig:
1484 if not self.isConfig:
1480 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1485 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1481 self.isConfig = True
1486 self.isConfig = True
1482 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1487 data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
1483 dataOut.flagNoData = True
1488 dataOut.flagNoData = True
1484
1489
1485 if self.__dataReady:
1490 if self.__dataReady:
1486 dataOut.nCohInt *= self.n
1491 dataOut.nCohInt *= self.n
1487 dataOut.dataPP_POW = data_intensity # S
1492 dataOut.dataPP_POW = data_intensity # S
1488 dataOut.dataPP_POWER = data_power # P
1493 dataOut.dataPP_POWER = data_power # P
1489 dataOut.dataPP_DOP = data_velocity
1494 dataOut.dataPP_DOP = data_velocity
1490 dataOut.dataPP_SNR = data_snrPP
1495 dataOut.dataPP_SNR = data_snrPP
1491 dataOut.dataPP_WIDTH = data_specwidth
1496 dataOut.dataPP_WIDTH = data_specwidth
1492 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1497 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1493 dataOut.utctime = avgdatatime
1498 dataOut.utctime = avgdatatime
1494 dataOut.flagNoData = False
1499 dataOut.flagNoData = False
1495 return dataOut
1500 return dataOut
1496
1501
1497
1502
1498
1503
1499 # import collections
1504 # import collections
1500 # from scipy.stats import mode
1505 # from scipy.stats import mode
1501 #
1506 #
1502 # class Synchronize(Operation):
1507 # class Synchronize(Operation):
1503 #
1508 #
1504 # isConfig = False
1509 # isConfig = False
1505 # __profIndex = 0
1510 # __profIndex = 0
1506 #
1511 #
1507 # def __init__(self, **kwargs):
1512 # def __init__(self, **kwargs):
1508 #
1513 #
1509 # Operation.__init__(self, **kwargs)
1514 # Operation.__init__(self, **kwargs)
1510 # # self.isConfig = False
1515 # # self.isConfig = False
1511 # self.__powBuffer = None
1516 # self.__powBuffer = None
1512 # self.__startIndex = 0
1517 # self.__startIndex = 0
1513 # self.__pulseFound = False
1518 # self.__pulseFound = False
1514 #
1519 #
1515 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1520 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1516 #
1521 #
1517 # #Read data
1522 # #Read data
1518 #
1523 #
1519 # powerdB = dataOut.getPower(channel = channel)
1524 # powerdB = dataOut.getPower(channel = channel)
1520 # noisedB = dataOut.getNoise(channel = channel)[0]
1525 # noisedB = dataOut.getNoise(channel = channel)[0]
1521 #
1526 #
1522 # self.__powBuffer.extend(powerdB.flatten())
1527 # self.__powBuffer.extend(powerdB.flatten())
1523 #
1528 #
1524 # dataArray = numpy.array(self.__powBuffer)
1529 # dataArray = numpy.array(self.__powBuffer)
1525 #
1530 #
1526 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1531 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1527 #
1532 #
1528 # maxValue = numpy.nanmax(filteredPower)
1533 # maxValue = numpy.nanmax(filteredPower)
1529 #
1534 #
1530 # if maxValue < noisedB + 10:
1535 # if maxValue < noisedB + 10:
1531 # #No se encuentra ningun pulso de transmision
1536 # #No se encuentra ningun pulso de transmision
1532 # return None
1537 # return None
1533 #
1538 #
1534 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1539 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1535 #
1540 #
1536 # if len(maxValuesIndex) < 2:
1541 # if len(maxValuesIndex) < 2:
1537 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1542 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1538 # return None
1543 # return None
1539 #
1544 #
1540 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1545 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1541 #
1546 #
1542 # #Seleccionar solo valores con un espaciamiento de nSamples
1547 # #Seleccionar solo valores con un espaciamiento de nSamples
1543 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1548 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1544 #
1549 #
1545 # if len(pulseIndex) < 2:
1550 # if len(pulseIndex) < 2:
1546 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1551 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1547 # return None
1552 # return None
1548 #
1553 #
1549 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1554 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1550 #
1555 #
1551 # #remover senales que se distancien menos de 10 unidades o muestras
1556 # #remover senales que se distancien menos de 10 unidades o muestras
1552 # #(No deberian existir IPP menor a 10 unidades)
1557 # #(No deberian existir IPP menor a 10 unidades)
1553 #
1558 #
1554 # realIndex = numpy.where(spacing > 10 )[0]
1559 # realIndex = numpy.where(spacing > 10 )[0]
1555 #
1560 #
1556 # if len(realIndex) < 2:
1561 # if len(realIndex) < 2:
1557 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1562 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1558 # return None
1563 # return None
1559 #
1564 #
1560 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1565 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1561 # realPulseIndex = pulseIndex[realIndex]
1566 # realPulseIndex = pulseIndex[realIndex]
1562 #
1567 #
1563 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1568 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1564 #
1569 #
1565 # print "IPP = %d samples" %period
1570 # print "IPP = %d samples" %period
1566 #
1571 #
1567 # self.__newNSamples = dataOut.nHeights #int(period)
1572 # self.__newNSamples = dataOut.nHeights #int(period)
1568 # self.__startIndex = int(realPulseIndex[0])
1573 # self.__startIndex = int(realPulseIndex[0])
1569 #
1574 #
1570 # return 1
1575 # return 1
1571 #
1576 #
1572 #
1577 #
1573 # def setup(self, nSamples, nChannels, buffer_size = 4):
1578 # def setup(self, nSamples, nChannels, buffer_size = 4):
1574 #
1579 #
1575 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1580 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1576 # maxlen = buffer_size*nSamples)
1581 # maxlen = buffer_size*nSamples)
1577 #
1582 #
1578 # bufferList = []
1583 # bufferList = []
1579 #
1584 #
1580 # for i in range(nChannels):
1585 # for i in range(nChannels):
1581 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1586 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1582 # maxlen = buffer_size*nSamples)
1587 # maxlen = buffer_size*nSamples)
1583 #
1588 #
1584 # bufferList.append(bufferByChannel)
1589 # bufferList.append(bufferByChannel)
1585 #
1590 #
1586 # self.__nSamples = nSamples
1591 # self.__nSamples = nSamples
1587 # self.__nChannels = nChannels
1592 # self.__nChannels = nChannels
1588 # self.__bufferList = bufferList
1593 # self.__bufferList = bufferList
1589 #
1594 #
1590 # def run(self, dataOut, channel = 0):
1595 # def run(self, dataOut, channel = 0):
1591 #
1596 #
1592 # if not self.isConfig:
1597 # if not self.isConfig:
1593 # nSamples = dataOut.nHeights
1598 # nSamples = dataOut.nHeights
1594 # nChannels = dataOut.nChannels
1599 # nChannels = dataOut.nChannels
1595 # self.setup(nSamples, nChannels)
1600 # self.setup(nSamples, nChannels)
1596 # self.isConfig = True
1601 # self.isConfig = True
1597 #
1602 #
1598 # #Append new data to internal buffer
1603 # #Append new data to internal buffer
1599 # for thisChannel in range(self.__nChannels):
1604 # for thisChannel in range(self.__nChannels):
1600 # bufferByChannel = self.__bufferList[thisChannel]
1605 # bufferByChannel = self.__bufferList[thisChannel]
1601 # bufferByChannel.extend(dataOut.data[thisChannel])
1606 # bufferByChannel.extend(dataOut.data[thisChannel])
1602 #
1607 #
1603 # if self.__pulseFound:
1608 # if self.__pulseFound:
1604 # self.__startIndex -= self.__nSamples
1609 # self.__startIndex -= self.__nSamples
1605 #
1610 #
1606 # #Finding Tx Pulse
1611 # #Finding Tx Pulse
1607 # if not self.__pulseFound:
1612 # if not self.__pulseFound:
1608 # indexFound = self.__findTxPulse(dataOut, channel)
1613 # indexFound = self.__findTxPulse(dataOut, channel)
1609 #
1614 #
1610 # if indexFound == None:
1615 # if indexFound == None:
1611 # dataOut.flagNoData = True
1616 # dataOut.flagNoData = True
1612 # return
1617 # return
1613 #
1618 #
1614 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1619 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1615 # self.__pulseFound = True
1620 # self.__pulseFound = True
1616 # self.__startIndex = indexFound
1621 # self.__startIndex = indexFound
1617 #
1622 #
1618 # #If pulse was found ...
1623 # #If pulse was found ...
1619 # for thisChannel in range(self.__nChannels):
1624 # for thisChannel in range(self.__nChannels):
1620 # bufferByChannel = self.__bufferList[thisChannel]
1625 # bufferByChannel = self.__bufferList[thisChannel]
1621 # #print self.__startIndex
1626 # #print self.__startIndex
1622 # x = numpy.array(bufferByChannel)
1627 # x = numpy.array(bufferByChannel)
1623 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1628 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1624 #
1629 #
1625 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1630 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1626 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1631 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1627 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1632 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1628 #
1633 #
1629 # dataOut.data = self.__arrayBuffer
1634 # dataOut.data = self.__arrayBuffer
1630 #
1635 #
1631 # self.__startIndex += self.__newNSamples
1636 # self.__startIndex += self.__newNSamples
1632 #
1637 #
1633 # return
1638 # return
General Comments 0
You need to be logged in to leave comments. Login now