##// END OF EJS Templates
schain
RSS
Update version, fix kwargs for self operations (methods), Add SendToWeb...
Juan C. Espinoza -
r906:a52f011a763e
parent child
Show More
Show file before | Show file after | Hide comments
@@ -1,98 +1,102
1 1 # Byte-compiled / optimized / DLL files
2 2 __pycache__/
3 3 *.py[cod]
4 4 *$py.class
5 5
6 6 # C extensions
7 7 *.so
8 8
9 9 # Distribution / packaging
10 10 .Python
11 11 env/
12 12 build/
13 13 develop-eggs/
14 14 dist/
15 15 downloads/
16 16 eggs/
17 17 .eggs/
18 18 lib/
19 19 lib64/
20 20 parts/
21 21 sdist/
22 22 var/
23 23 wheels/
24 24 *.egg-info/
25 25 .installed.cfg
26 26 *.egg
27 27
28 28 # PyInstaller
29 29 # Usually these files are written by a python script from a template
30 30 # before PyInstaller builds the exe, so as to inject date/other infos into it.
31 31 *.manifest
32 32 *.spec
33 33
34 34 # Installer logs
35 35 pip-log.txt
36 36 pip-delete-this-directory.txt
37 37
38 38 # Unit test / coverage reports
39 39 htmlcov/
40 40 .tox/
41 41 .coverage
42 42 .coverage.*
43 43 .cache
44 44 nosetests.xml
45 45 coverage.xml
46 46 *,cover
47 47 .hypothesis/
48 48
49 49 # Translations
50 50 *.mo
51 51 *.pot
52 52
53 53 # Django stuff:
54 54 *.log
55 55 local_settings.py
56 56
57 57 # Flask stuff:
58 58 instance/
59 59 .webassets-cache
60 60
61 61 # Scrapy stuff:
62 62 .scrapy
63 63
64 64 # Sphinx documentation
65 65 docs/_build/
66 66
67 67 # PyBuilder
68 68 target/
69 69
70 70 # Jupyter Notebook
71 71 .ipynb_checkpoints
72 72
73 73 # pyenv
74 74 .python-version
75 75
76 76 # celery beat schedule file
77 77 celerybeat-schedule
78 78
79 79 # SageMath parsed files
80 80 *.sage.py
81 81
82 82 # dotenv
83 83 .env
84 84
85 85 # virtualenv
86 86 .venv
87 87 venv/
88 88 ENV/
89 89
90 90 # Spyder project settings
91 91 .spyderproject
92 92 .spyproject
93 93
94 94 # Rope project settings
95 95 .ropeproject
96 96
97 97 # mkdocs documentation
98 98 /site
99
100 # eclipse
101 .project
102 .pydevproject
Show file before | Show file after | Hide comments
@@ -1,7 +1,7
1 1 '''
2 2 Created on Feb 7, 2012
3 3
4 4 @author $Author$
5 5 @version $Id$
6 6 '''
7 __version__ = "2.2.5" No newline at end of file
7 __version__ = "2.3" No newline at end of file
Show file before | Show file after | Hide comments
@@ -1,1318 +1,1323
1 1 '''
2 2 Created on September , 2012
3 3 @author:
4 4 '''
5 5
6 6 import sys
7 7 import ast
8 8 import datetime
9 9 import traceback
10 10 import math
11 11 from multiprocessing import Process, Queue, cpu_count
12 12
13 13 import schainpy
14 14 import schainpy.admin
15 15
16 16 from xml.etree.ElementTree import ElementTree, Element, SubElement, tostring
17 17 from xml.dom import minidom
18 18
19 19 from schainpy.model import *
20 20 from time import sleep
21 21
22 22 def prettify(elem):
23 23 """Return a pretty-printed XML string for the Element.
24 24 """
25 25 rough_string = tostring(elem, 'utf-8')
26 26 reparsed = minidom.parseString(rough_string)
27 27 return reparsed.toprettyxml(indent=" ")
28 28
29 29 def multiSchain(child, nProcess=cpu_count(), startDate=None, endDate=None, receiver=None):
30 30 skip = 0
31 31 cursor = 0
32 32 nFiles = None
33 33 processes = []
34 34
35 35 dt1 = datetime.datetime.strptime(startDate, '%Y/%m/%d')
36 36 dt2 = datetime.datetime.strptime(endDate, '%Y/%m/%d')
37 37 days = (dt2 - dt1).days
38 38 print days
39 39 for day in range(days+1):
40 40 skip = 0
41 41 cursor = 0
42 42 q = Queue()
43 43 processes = []
44 44 dt = (dt1 + datetime.timedelta(day)).strftime('%Y/%m/%d')
45 45 firstProcess = Process(target=child, args=(cursor, skip, q, dt))
46 46 firstProcess.start()
47 47 nFiles = q.get()
48 48 firstProcess.terminate()
49 49 skip = int(math.ceil(nFiles/nProcess))
50 50 try:
51 51 while True:
52 52 processes.append(Process(target=child, args=(cursor, skip, q, dt)))
53 53 processes[cursor].start()
54 54 if nFiles < cursor*skip:
55 55 break
56 56 cursor += 1
57 57 except KeyboardInterrupt:
58 58 for process in processes:
59 59 process.terminate()
60 60 process.join()
61 61 for process in processes:
62 62 process.join()
63 63 # process.terminate()
64 64 sleep(3)
65 65
66 66 try:
67 67 while True:
68 68 pass
69 69 except KeyboardInterrupt:
70 70 for process in processes:
71 71 process.terminate()
72 72 process.join()
73 73
74 74 class ParameterConf():
75 75
76 76 id = None
77 77 name = None
78 78 value = None
79 79 format = None
80 80
81 81 __formated_value = None
82 82
83 83 ELEMENTNAME = 'Parameter'
84 84
85 85 def __init__(self):
86 86
87 87 self.format = 'str'
88 88
89 89 def getElementName(self):
90 90
91 91 return self.ELEMENTNAME
92 92
93 93 def getValue(self):
94 94
95 95 value = self.value
96 96 format = self.format
97 97
98 98 if self.__formated_value != None:
99 99
100 100 return self.__formated_value
101 101
102 102 if format == 'obj':
103 103 return value
104 104
105 105 if format == 'str':
106 106 self.__formated_value = str(value)
107 107 return self.__formated_value
108 108
109 109 if value == '':
110 110 raise ValueError, "%s: This parameter value is empty" %self.name
111 111
112 112 if format == 'list':
113 113 strList = value.split(',')
114 114
115 115 self.__formated_value = strList
116 116
117 117 return self.__formated_value
118 118
119 119 if format == 'intlist':
120 120 """
121 121 Example:
122 122 value = (0,1,2)
123 123 """
124 124
125 125 new_value = ast.literal_eval(value)
126 126
127 127 if type(new_value) not in (tuple, list):
128 128 new_value = [int(new_value)]
129 129
130 130 self.__formated_value = new_value
131 131
132 132 return self.__formated_value
133 133
134 134 if format == 'floatlist':
135 135 """
136 136 Example:
137 137 value = (0.5, 1.4, 2.7)
138 138 """
139 139
140 140 new_value = ast.literal_eval(value)
141 141
142 142 if type(new_value) not in (tuple, list):
143 143 new_value = [float(new_value)]
144 144
145 145 self.__formated_value = new_value
146 146
147 147 return self.__formated_value
148 148
149 149 if format == 'date':
150 150 strList = value.split('/')
151 151 intList = [int(x) for x in strList]
152 152 date = datetime.date(intList[0], intList[1], intList[2])
153 153
154 154 self.__formated_value = date
155 155
156 156 return self.__formated_value
157 157
158 158 if format == 'time':
159 159 strList = value.split(':')
160 160 intList = [int(x) for x in strList]
161 161 time = datetime.time(intList[0], intList[1], intList[2])
162 162
163 163 self.__formated_value = time
164 164
165 165 return self.__formated_value
166 166
167 167 if format == 'pairslist':
168 168 """
169 169 Example:
170 170 value = (0,1),(1,2)
171 171 """
172 172
173 173 new_value = ast.literal_eval(value)
174 174
175 175 if type(new_value) not in (tuple, list):
176 176 raise ValueError, "%s has to be a tuple or list of pairs" %value
177 177
178 178 if type(new_value[0]) not in (tuple, list):
179 179 if len(new_value) != 2:
180 180 raise ValueError, "%s has to be a tuple or list of pairs" %value
181 181 new_value = [new_value]
182 182
183 183 for thisPair in new_value:
184 184 if len(thisPair) != 2:
185 185 raise ValueError, "%s has to be a tuple or list of pairs" %value
186 186
187 187 self.__formated_value = new_value
188 188
189 189 return self.__formated_value
190 190
191 191 if format == 'multilist':
192 192 """
193 193 Example:
194 194 value = (0,1,2),(3,4,5)
195 195 """
196 196 multiList = ast.literal_eval(value)
197 197
198 198 if type(multiList[0]) == int:
199 199 multiList = ast.literal_eval("(" + value + ")")
200 200
201 201 self.__formated_value = multiList
202 202
203 203 return self.__formated_value
204 204
205 205 if format == 'bool':
206 206 value = int(value)
207 207
208 208 if format == 'int':
209 209 value = float(value)
210 210
211 211 format_func = eval(format)
212 212
213 213 self.__formated_value = format_func(value)
214 214
215 215 return self.__formated_value
216 216
217 217 def updateId(self, new_id):
218 218
219 219 self.id = str(new_id)
220 220
221 221 def setup(self, id, name, value, format='str'):
222 222
223 223 self.id = str(id)
224 224 self.name = name
225 225 if format == 'obj':
226 226 self.value = value
227 227 else:
228 228 self.value = str(value)
229 229 self.format = str.lower(format)
230 230
231 231 self.getValue()
232 232
233 233 return 1
234 234
235 235 def update(self, name, value, format='str'):
236 236
237 237 self.name = name
238 238 self.value = str(value)
239 239 self.format = format
240 240
241 241 def makeXml(self, opElement):
242 242 if self.name not in ('queue',):
243 243 parmElement = SubElement(opElement, self.ELEMENTNAME)
244 244 parmElement.set('id', str(self.id))
245 245 parmElement.set('name', self.name)
246 246 parmElement.set('value', self.value)
247 247 parmElement.set('format', self.format)
248 248
249 249 def readXml(self, parmElement):
250 250
251 251 self.id = parmElement.get('id')
252 252 self.name = parmElement.get('name')
253 253 self.value = parmElement.get('value')
254 254 self.format = str.lower(parmElement.get('format'))
255 255
256 256 #Compatible with old signal chain version
257 257 if self.format == 'int' and self.name == 'idfigure':
258 258 self.name = 'id'
259 259
260 260 def printattr(self):
261 261
262 262 print "Parameter[%s]: name = %s, value = %s, format = %s" %(self.id, self.name, self.value, self.format)
263 263
264 264 class OperationConf():
265 265
266 266 id = None
267 267 name = None
268 268 priority = None
269 269 type = None
270 270
271 271 parmConfObjList = []
272 272
273 273 ELEMENTNAME = 'Operation'
274 274
275 275 def __init__(self):
276 276
277 277 self.id = '0'
278 278 self.name = None
279 279 self.priority = None
280 280 self.type = 'self'
281 281
282 282
283 283 def __getNewId(self):
284 284
285 285 return int(self.id)*10 + len(self.parmConfObjList) + 1
286 286
287 287 def updateId(self, new_id):
288 288
289 289 self.id = str(new_id)
290 290
291 291 n = 1
292 292 for parmObj in self.parmConfObjList:
293 293
294 294 idParm = str(int(new_id)*10 + n)
295 295 parmObj.updateId(idParm)
296 296
297 297 n += 1
298 298
299 299 def getElementName(self):
300 300
301 301 return self.ELEMENTNAME
302 302
303 303 def getParameterObjList(self):
304 304
305 305 return self.parmConfObjList
306 306
307 307 def getParameterObj(self, parameterName):
308 308
309 309 for parmConfObj in self.parmConfObjList:
310 310
311 311 if parmConfObj.name != parameterName:
312 312 continue
313 313
314 314 return parmConfObj
315 315
316 316 return None
317 317
318 318 def getParameterObjfromValue(self, parameterValue):
319 319
320 320 for parmConfObj in self.parmConfObjList:
321 321
322 322 if parmConfObj.getValue() != parameterValue:
323 323 continue
324 324
325 325 return parmConfObj.getValue()
326 326
327 327 return None
328 328
329 329 def getParameterValue(self, parameterName):
330 330
331 331 parameterObj = self.getParameterObj(parameterName)
332 332
333 333 # if not parameterObj:
334 334 # return None
335 335
336 336 value = parameterObj.getValue()
337 337
338 338 return value
339 339
340 340
341 341 def getKwargs(self):
342 342
343 343 kwargs = {}
344 344
345 345 for parmConfObj in self.parmConfObjList:
346 346 if self.name == 'run' and parmConfObj.name == 'datatype':
347 347 continue
348 348
349 349 kwargs[parmConfObj.name] = parmConfObj.getValue()
350 350
351 351 return kwargs
352 352
353 353 def setup(self, id, name, priority, type):
354 354
355 355 self.id = str(id)
356 356 self.name = name
357 357 self.type = type
358 358 self.priority = priority
359 359
360 360 self.parmConfObjList = []
361 361
362 362 def removeParameters(self):
363 363
364 364 for obj in self.parmConfObjList:
365 365 del obj
366 366
367 367 self.parmConfObjList = []
368 368
369 369 def addParameter(self, name, value, format='str'):
370 370
371 371 id = self.__getNewId()
372 372
373 373 parmConfObj = ParameterConf()
374 374 if not parmConfObj.setup(id, name, value, format):
375 375 return None
376 376
377 377 self.parmConfObjList.append(parmConfObj)
378 378
379 379 return parmConfObj
380 380
381 381 def changeParameter(self, name, value, format='str'):
382 382
383 383 parmConfObj = self.getParameterObj(name)
384 384 parmConfObj.update(name, value, format)
385 385
386 386 return parmConfObj
387 387
388 388 def makeXml(self, procUnitElement):
389 389
390 390 opElement = SubElement(procUnitElement, self.ELEMENTNAME)
391 391 opElement.set('id', str(self.id))
392 392 opElement.set('name', self.name)
393 393 opElement.set('type', self.type)
394 394 opElement.set('priority', str(self.priority))
395 395
396 396 for parmConfObj in self.parmConfObjList:
397 397 parmConfObj.makeXml(opElement)
398 398
399 399 def readXml(self, opElement):
400 400
401 401 self.id = opElement.get('id')
402 402 self.name = opElement.get('name')
403 403 self.type = opElement.get('type')
404 404 self.priority = opElement.get('priority')
405 405
406 406 #Compatible with old signal chain version
407 407 #Use of 'run' method instead 'init'
408 408 if self.type == 'self' and self.name == 'init':
409 409 self.name = 'run'
410 410
411 411 self.parmConfObjList = []
412 412
413 413 parmElementList = opElement.iter(ParameterConf().getElementName())
414 414
415 415 for parmElement in parmElementList:
416 416 parmConfObj = ParameterConf()
417 417 parmConfObj.readXml(parmElement)
418 418
419 419 #Compatible with old signal chain version
420 420 #If an 'plot' OPERATION is found, changes name operation by the value of its type PARAMETER
421 421 if self.type != 'self' and self.name == 'Plot':
422 422 if parmConfObj.format == 'str' and parmConfObj.name == 'type':
423 423 self.name = parmConfObj.value
424 424 continue
425 425
426 426 self.parmConfObjList.append(parmConfObj)
427 427
428 428 def printattr(self):
429 429
430 430 print "%s[%s]: name = %s, type = %s, priority = %s" %(self.ELEMENTNAME,
431 431 self.id,
432 432 self.name,
433 433 self.type,
434 434 self.priority)
435 435
436 436 for parmConfObj in self.parmConfObjList:
437 437 parmConfObj.printattr()
438 438
439 439 def createObject(self, plotter_queue=None):
440 440
441 if self.type == 'self':
441
442 if self.type == 'self':
442 443 raise ValueError, "This operation type cannot be created"
443 444
444 445 if self.type == 'plotter':
445 446 #Plotter(plotter_name)
446 447 if not plotter_queue:
447 448 raise ValueError, "plotter_queue is not defined. Use:\nmyProject = Project()\nmyProject.setPlotterQueue(plotter_queue)"
448 449
449 450 opObj = Plotter(self.name, plotter_queue)
450 451
451 452 if self.type == 'external' or self.type == 'other':
452 print self.name
453
453 454 className = eval(self.name)
454 455 kwargs = self.getKwargs()
455 print kwargs
456
456 457 opObj = className(**kwargs)
457 458
458 459 return opObj
459 460
460 461
461 462 class ProcUnitConf():
462 463
463 464 id = None
464 465 name = None
465 466 datatype = None
466 467 inputId = None
467 468 parentId = None
468 469
469 470 opConfObjList = []
470 471
471 472 procUnitObj = None
472 473 opObjList = []
473 474
474 475 ELEMENTNAME = 'ProcUnit'
475 476
476 477 def __init__(self):
477 478
478 479 self.id = None
479 480 self.datatype = None
480 481 self.name = None
481 482 self.inputId = None
482 483
483 484 self.opConfObjList = []
484 485
485 486 self.procUnitObj = None
486 487 self.opObjDict = {}
487 488
488 489 def __getPriority(self):
489 490
490 491 return len(self.opConfObjList)+1
491 492
492 493 def __getNewId(self):
493 494
494 495 return int(self.id)*10 + len(self.opConfObjList) + 1
495 496
496 497 def getElementName(self):
497 498
498 499 return self.ELEMENTNAME
499 500
500 501 def getId(self):
501 502
502 503 return self.id
503 504
504 505 def updateId(self, new_id, parentId=parentId):
505 506
506 507
507 508 new_id = int(parentId)*10 + (int(self.id) % 10)
508 509 new_inputId = int(parentId)*10 + (int(self.inputId) % 10)
509 510
510 511 #If this proc unit has not inputs
511 512 if self.inputId == '0':
512 513 new_inputId = 0
513 514
514 515 n = 1
515 516 for opConfObj in self.opConfObjList:
516 517
517 518 idOp = str(int(new_id)*10 + n)
518 519 opConfObj.updateId(idOp)
519 520
520 521 n += 1
521 522
522 523 self.parentId = str(parentId)
523 524 self.id = str(new_id)
524 525 self.inputId = str(new_inputId)
525 526
526 527
527 528 def getInputId(self):
528 529
529 530 return self.inputId
530 531
531 532 def getOperationObjList(self):
532 533
533 534 return self.opConfObjList
534 535
535 536 def getOperationObj(self, name=None):
536 537
537 538 for opConfObj in self.opConfObjList:
538 539
539 540 if opConfObj.name != name:
540 541 continue
541 542
542 543 return opConfObj
543 544
544 545 return None
545 546
546 547 def getOpObjfromParamValue(self, value=None):
547 548
548 549 for opConfObj in self.opConfObjList:
549 550 if opConfObj.getParameterObjfromValue(parameterValue=value) != value:
550 551 continue
551 552 return opConfObj
552 553 return None
553 554
554 555 def getProcUnitObj(self):
555 556
556 557 return self.procUnitObj
557 558
558 559 def setup(self, id, name, datatype, inputId, parentId=None):
559 560
560 561 #Compatible with old signal chain version
561 562 if datatype==None and name==None:
562 563 raise ValueError, "datatype or name should be defined"
563 564
564 565 if name==None:
565 566 if 'Proc' in datatype:
566 567 name = datatype
567 568 else:
568 569 name = '%sProc' %(datatype)
569 570
570 571 if datatype==None:
571 572 datatype = name.replace('Proc','')
572 573
573 574 self.id = str(id)
574 575 self.name = name
575 576 self.datatype = datatype
576 577 self.inputId = inputId
577 578 self.parentId = parentId
578 579
579 580 self.opConfObjList = []
580 581
581 582 self.addOperation(name='run', optype='self')
582 583
583 584 def removeOperations(self):
584 585
585 586 for obj in self.opConfObjList:
586 587 del obj
587 588
588 589 self.opConfObjList = []
589 590 self.addOperation(name='run')
590 591
591 592 def addParameter(self, **kwargs):
592 593 '''
593 594 Add parameters to "run" operation
594 595 '''
595 596 opObj = self.opConfObjList[0]
596 597
597 598 opObj.addParameter(**kwargs)
598 599
599 600 return opObj
600 601
601 602 def addOperation(self, name, optype='self'):
602 603
603 604 id = self.__getNewId()
604 605 priority = self.__getPriority()
605 606
606 607 opConfObj = OperationConf()
607 608 opConfObj.setup(id, name=name, priority=priority, type=optype)
608 609
609 610 self.opConfObjList.append(opConfObj)
610 611
611 612 return opConfObj
612 613
613 614 def makeXml(self, projectElement):
614 615
615 616 procUnitElement = SubElement(projectElement, self.ELEMENTNAME)
616 617 procUnitElement.set('id', str(self.id))
617 618 procUnitElement.set('name', self.name)
618 619 procUnitElement.set('datatype', self.datatype)
619 620 procUnitElement.set('inputId', str(self.inputId))
620 621
621 622 for opConfObj in self.opConfObjList:
622 623 opConfObj.makeXml(procUnitElement)
623 624
624 625 def readXml(self, upElement):
625 626
626 627 self.id = upElement.get('id')
627 628 self.name = upElement.get('name')
628 629 self.datatype = upElement.get('datatype')
629 630 self.inputId = upElement.get('inputId')
630 631
631 632 if self.ELEMENTNAME == "ReadUnit":
632 633 self.datatype = self.datatype.replace("Reader", "")
633 634
634 635 if self.ELEMENTNAME == "ProcUnit":
635 636 self.datatype = self.datatype.replace("Proc", "")
636 637
637 638 if self.inputId == 'None':
638 639 self.inputId = '0'
639 640
640 641 self.opConfObjList = []
641 642
642 643 opElementList = upElement.iter(OperationConf().getElementName())
643 644
644 645 for opElement in opElementList:
645 646 opConfObj = OperationConf()
646 647 opConfObj.readXml(opElement)
647 648 self.opConfObjList.append(opConfObj)
648 649
649 650 def printattr(self):
650 651
651 652 print "%s[%s]: name = %s, datatype = %s, inputId = %s" %(self.ELEMENTNAME,
652 653 self.id,
653 654 self.name,
654 655 self.datatype,
655 656 self.inputId)
656 657
657 658 for opConfObj in self.opConfObjList:
658 659 opConfObj.printattr()
659 660
660 661
661 662 def getKwargs(self):
662 663
663 664 opObj = self.opConfObjList[0]
664 665 kwargs = opObj.getKwargs()
665 666
666 667 return kwargs
667 668
668 669 def createObjects(self, plotter_queue=None):
669 670
670 671 className = eval(self.name)
671 672 kwargs = self.getKwargs()
672 673 procUnitObj = className(**kwargs)
673 674
674 for opConfObj in self.opConfObjList:
675
676 if opConfObj.type == 'self':
675 for opConfObj in self.opConfObjList:
676
677 if opConfObj.type=='self' and self.name=='run':
678 continue
679 elif opConfObj.type=='self':
680 procUnitObj.addOperationKwargs(opConfObj.id, **opConfObj.getKwargs())
677 681 continue
678 682
679 683 opObj = opConfObj.createObject(plotter_queue)
680 684
681 685 self.opObjDict[opConfObj.id] = opObj
686
682 687 procUnitObj.addOperation(opObj, opConfObj.id)
683 688
684 689 self.procUnitObj = procUnitObj
685 690
686 691 return procUnitObj
687 692
688 693 def run(self):
689 694
690 695 is_ok = False
691 696
692 697 for opConfObj in self.opConfObjList:
693 698
694 699 kwargs = {}
695 700 for parmConfObj in opConfObj.getParameterObjList():
696 701 if opConfObj.name == 'run' and parmConfObj.name == 'datatype':
697 702 continue
698 703
699 704 kwargs[parmConfObj.name] = parmConfObj.getValue()
700 705
701 706 #ini = time.time()
702 707
703 708 #print "\tRunning the '%s' operation with %s" %(opConfObj.name, opConfObj.id)
704 709 sts = self.procUnitObj.call(opType = opConfObj.type,
705 710 opName = opConfObj.name,
706 711 opId = opConfObj.id,
707 712 )
708 713
709 714 # total_time = time.time() - ini
710 715 #
711 716 # if total_time > 0.002:
712 717 # print "%s::%s took %f seconds" %(self.name, opConfObj.name, total_time)
713 718
714 719 is_ok = is_ok or sts
715 720
716 721 return is_ok
717 722
718 723 def close(self):
719 724
720 725 for opConfObj in self.opConfObjList:
721 726 if opConfObj.type == 'self':
722 727 continue
723 728
724 729 opObj = self.procUnitObj.getOperationObj(opConfObj.id)
725 730 opObj.close()
726 731
727 732 self.procUnitObj.close()
728 733
729 734 return
730 735
731 736 class ReadUnitConf(ProcUnitConf):
732 737
733 738 path = None
734 739 startDate = None
735 740 endDate = None
736 741 startTime = None
737 742 endTime = None
738 743
739 744 ELEMENTNAME = 'ReadUnit'
740 745
741 746 def __init__(self):
742 747
743 748 self.id = None
744 749 self.datatype = None
745 750 self.name = None
746 751 self.inputId = None
747 752
748 753 self.parentId = None
749 754
750 755 self.opConfObjList = []
751 756 self.opObjList = []
752 757
753 758 def getElementName(self):
754 759
755 760 return self.ELEMENTNAME
756 761
757 762 def setup(self, id, name, datatype, path, startDate="", endDate="", startTime="", endTime="", parentId=None, queue=None, **kwargs):
758 763
759 764 #Compatible with old signal chain version
760 765 if datatype==None and name==None:
761 766 raise ValueError, "datatype or name should be defined"
762 767
763 768 if name==None:
764 769 if 'Reader' in datatype:
765 770 name = datatype
766 771 else:
767 772 name = '%sReader' %(datatype)
768 773
769 774 if datatype==None:
770 775 datatype = name.replace('Reader','')
771 776
772 777 self.id = id
773 778 self.name = name
774 779 self.datatype = datatype
775 780
776 781 self.path = os.path.abspath(path)
777 782 self.startDate = startDate
778 783 self.endDate = endDate
779 784 self.startTime = startTime
780 785 self.endTime = endTime
781 786
782 787 self.inputId = '0'
783 788 self.parentId = parentId
784 789 self.queue = queue
785 790 self.addRunOperation(**kwargs)
786 791
787 792 def update(self, datatype, path, startDate, endDate, startTime, endTime, parentId=None, name=None, **kwargs):
788 793
789 794 #Compatible with old signal chain version
790 795 if datatype==None and name==None:
791 796 raise ValueError, "datatype or name should be defined"
792 797
793 798 if name==None:
794 799 if 'Reader' in datatype:
795 800 name = datatype
796 801 else:
797 802 name = '%sReader' %(datatype)
798 803
799 804 if datatype==None:
800 805 datatype = name.replace('Reader','')
801 806
802 807 self.datatype = datatype
803 808 self.name = name
804 809 self.path = path
805 810 self.startDate = startDate
806 811 self.endDate = endDate
807 812 self.startTime = startTime
808 813 self.endTime = endTime
809 814
810 815 self.inputId = '0'
811 816 self.parentId = parentId
812 817
813 818 self.updateRunOperation(**kwargs)
814 819
815 820 def removeOperations(self):
816 821
817 822 for obj in self.opConfObjList:
818 823 del obj
819 824
820 825 self.opConfObjList = []
821 826
822 827 def addRunOperation(self, **kwargs):
823 828
824 829 opObj = self.addOperation(name = 'run', optype = 'self')
825 830
826 831 opObj.addParameter(name='datatype' , value=self.datatype, format='str')
827 832 opObj.addParameter(name='path' , value=self.path, format='str')
828 833 opObj.addParameter(name='startDate' , value=self.startDate, format='date')
829 834 opObj.addParameter(name='endDate' , value=self.endDate, format='date')
830 835 opObj.addParameter(name='startTime' , value=self.startTime, format='time')
831 836 opObj.addParameter(name='endTime' , value=self.endTime, format='time')
832 837 opObj.addParameter(name='queue' , value=self.queue, format='obj')
833 838
834 839 for key, value in kwargs.items():
835 840 opObj.addParameter(name=key, value=value, format=type(value).__name__)
836 841
837 842 return opObj
838 843
839 844 def updateRunOperation(self, **kwargs):
840 845
841 846 opObj = self.getOperationObj(name = 'run')
842 847 opObj.removeParameters()
843 848
844 849 opObj.addParameter(name='datatype' , value=self.datatype, format='str')
845 850 opObj.addParameter(name='path' , value=self.path, format='str')
846 851 opObj.addParameter(name='startDate' , value=self.startDate, format='date')
847 852 opObj.addParameter(name='endDate' , value=self.endDate, format='date')
848 853 opObj.addParameter(name='startTime' , value=self.startTime, format='time')
849 854 opObj.addParameter(name='endTime' , value=self.endTime, format='time')
850 855
851 856 for key, value in kwargs.items():
852 857 opObj.addParameter(name=key, value=value, format=type(value).__name__)
853 858
854 859 return opObj
855 860
856 861 # def makeXml(self, projectElement):
857 862 #
858 863 # procUnitElement = SubElement(projectElement, self.ELEMENTNAME)
859 864 # procUnitElement.set('id', str(self.id))
860 865 # procUnitElement.set('name', self.name)
861 866 # procUnitElement.set('datatype', self.datatype)
862 867 # procUnitElement.set('inputId', str(self.inputId))
863 868 #
864 869 # for opConfObj in self.opConfObjList:
865 870 # opConfObj.makeXml(procUnitElement)
866 871
867 872 def readXml(self, upElement):
868 873
869 874 self.id = upElement.get('id')
870 875 self.name = upElement.get('name')
871 876 self.datatype = upElement.get('datatype')
872 877 self.inputId = upElement.get('inputId')
873 878
874 879 if self.ELEMENTNAME == "ReadUnit":
875 880 self.datatype = self.datatype.replace("Reader", "")
876 881
877 882 if self.inputId == 'None':
878 883 self.inputId = '0'
879 884
880 885 self.opConfObjList = []
881 886
882 887 opElementList = upElement.iter(OperationConf().getElementName())
883 888
884 889 for opElement in opElementList:
885 890 opConfObj = OperationConf()
886 891 opConfObj.readXml(opElement)
887 892 self.opConfObjList.append(opConfObj)
888 893
889 894 if opConfObj.name == 'run':
890 895 self.path = opConfObj.getParameterValue('path')
891 896 self.startDate = opConfObj.getParameterValue('startDate')
892 897 self.endDate = opConfObj.getParameterValue('endDate')
893 898 self.startTime = opConfObj.getParameterValue('startTime')
894 899 self.endTime = opConfObj.getParameterValue('endTime')
895 900
896 901 class Project():
897 902
898 903 id = None
899 904 name = None
900 905 description = None
901 906 filename = None
902 907
903 908 procUnitConfObjDict = None
904 909
905 910 ELEMENTNAME = 'Project'
906 911
907 912 plotterQueue = None
908 913
909 914 def __init__(self, plotter_queue=None):
910 915
911 916 self.id = None
912 917 self.name = None
913 918 self.description = None
914 919
915 920 self.plotterQueue = plotter_queue
916 921
917 922 self.procUnitConfObjDict = {}
918 923
919 924 def __getNewId(self):
920 925
921 926 idList = self.procUnitConfObjDict.keys()
922 927
923 928 id = int(self.id)*10
924 929
925 930 while True:
926 931 id += 1
927 932
928 933 if str(id) in idList:
929 934 continue
930 935
931 936 break
932 937
933 938 return str(id)
934 939
935 940 def getElementName(self):
936 941
937 942 return self.ELEMENTNAME
938 943
939 944 def getId(self):
940 945
941 946 return self.id
942 947
943 948 def updateId(self, new_id):
944 949
945 950 self.id = str(new_id)
946 951
947 952 keyList = self.procUnitConfObjDict.keys()
948 953 keyList.sort()
949 954
950 955 n = 1
951 956 newProcUnitConfObjDict = {}
952 957
953 958 for procKey in keyList:
954 959
955 960 procUnitConfObj = self.procUnitConfObjDict[procKey]
956 961 idProcUnit = str(int(self.id)*10 + n)
957 962 procUnitConfObj.updateId(idProcUnit, parentId = self.id)
958 963
959 964 newProcUnitConfObjDict[idProcUnit] = procUnitConfObj
960 965 n += 1
961 966
962 967 self.procUnitConfObjDict = newProcUnitConfObjDict
963 968
964 969 def setup(self, id, name, description):
965 970
966 971 self.id = str(id)
967 972 self.name = name
968 973 self.description = description
969 974
970 975 def update(self, name, description):
971 976
972 977 self.name = name
973 978 self.description = description
974 979
975 980 def addReadUnit(self, id=None, datatype=None, name=None, **kwargs):
976 981
977 982 if id is None:
978 983 idReadUnit = self.__getNewId()
979 984 else:
980 985 idReadUnit = str(id)
981 986
982 987 readUnitConfObj = ReadUnitConf()
983 988 readUnitConfObj.setup(idReadUnit, name, datatype, parentId=self.id, **kwargs)
984 989
985 990 self.procUnitConfObjDict[readUnitConfObj.getId()] = readUnitConfObj
986 991
987 992 return readUnitConfObj
988 993
989 994 def addProcUnit(self, inputId='0', datatype=None, name=None):
990 995
991 996 idProcUnit = self.__getNewId()
992 997
993 998 procUnitConfObj = ProcUnitConf()
994 999 procUnitConfObj.setup(idProcUnit, name, datatype, inputId, parentId=self.id)
995 1000
996 1001 self.procUnitConfObjDict[procUnitConfObj.getId()] = procUnitConfObj
997 1002
998 1003 return procUnitConfObj
999 1004
1000 1005 def removeProcUnit(self, id):
1001 1006
1002 1007 if id in self.procUnitConfObjDict.keys():
1003 1008 self.procUnitConfObjDict.pop(id)
1004 1009
1005 1010 def getReadUnitId(self):
1006 1011
1007 1012 readUnitConfObj = self.getReadUnitObj()
1008 1013
1009 1014 return readUnitConfObj.id
1010 1015
1011 1016 def getReadUnitObj(self):
1012 1017
1013 1018 for obj in self.procUnitConfObjDict.values():
1014 1019 if obj.getElementName() == "ReadUnit":
1015 1020 return obj
1016 1021
1017 1022 return None
1018 1023
1019 1024 def getProcUnitObj(self, id=None, name=None):
1020 1025
1021 1026 if id != None:
1022 1027 return self.procUnitConfObjDict[id]
1023 1028
1024 1029 if name != None:
1025 1030 return self.getProcUnitObjByName(name)
1026 1031
1027 1032 return None
1028 1033
1029 1034 def getProcUnitObjByName(self, name):
1030 1035
1031 1036 for obj in self.procUnitConfObjDict.values():
1032 1037 if obj.name == name:
1033 1038 return obj
1034 1039
1035 1040 return None
1036 1041
1037 1042 def procUnitItems(self):
1038 1043
1039 1044 return self.procUnitConfObjDict.items()
1040 1045
1041 1046 def makeXml(self):
1042 1047
1043 1048 projectElement = Element('Project')
1044 1049 projectElement.set('id', str(self.id))
1045 1050 projectElement.set('name', self.name)
1046 1051 projectElement.set('description', self.description)
1047 1052
1048 1053 for procUnitConfObj in self.procUnitConfObjDict.values():
1049 1054 procUnitConfObj.makeXml(projectElement)
1050 1055
1051 1056 self.projectElement = projectElement
1052 1057
1053 1058 def writeXml(self, filename=None):
1054 1059
1055 1060 if filename == None:
1056 1061 if self.filename:
1057 1062 filename = self.filename
1058 1063 else:
1059 1064 filename = "schain.xml"
1060 1065
1061 1066 if not filename:
1062 1067 print "filename has not been defined. Use setFilename(filename) for do it."
1063 1068 return 0
1064 1069
1065 1070 abs_file = os.path.abspath(filename)
1066 1071
1067 1072 if not os.access(os.path.dirname(abs_file), os.W_OK):
1068 1073 print "No write permission on %s" %os.path.dirname(abs_file)
1069 1074 return 0
1070 1075
1071 1076 if os.path.isfile(abs_file) and not(os.access(abs_file, os.W_OK)):
1072 1077 print "File %s already exists and it could not be overwriten" %abs_file
1073 1078 return 0
1074 1079
1075 1080 self.makeXml()
1076 1081
1077 1082 ElementTree(self.projectElement).write(abs_file, method='xml')
1078 1083
1079 1084 self.filename = abs_file
1080 1085
1081 1086 return 1
1082 1087
1083 1088 def readXml(self, filename = None):
1084 1089
1085 1090 if not filename:
1086 1091 print "filename is not defined"
1087 1092 return 0
1088 1093
1089 1094 abs_file = os.path.abspath(filename)
1090 1095
1091 1096 if not os.path.isfile(abs_file):
1092 1097 print "%s file does not exist" %abs_file
1093 1098 return 0
1094 1099
1095 1100 self.projectElement = None
1096 1101 self.procUnitConfObjDict = {}
1097 1102
1098 1103 try:
1099 1104 self.projectElement = ElementTree().parse(abs_file)
1100 1105 except:
1101 1106 print "Error reading %s, verify file format" %filename
1102 1107 return 0
1103 1108
1104 1109 self.project = self.projectElement.tag
1105 1110
1106 1111 self.id = self.projectElement.get('id')
1107 1112 self.name = self.projectElement.get('name')
1108 1113 self.description = self.projectElement.get('description')
1109 1114
1110 1115 readUnitElementList = self.projectElement.iter(ReadUnitConf().getElementName())
1111 1116
1112 1117 for readUnitElement in readUnitElementList:
1113 1118 readUnitConfObj = ReadUnitConf()
1114 1119 readUnitConfObj.readXml(readUnitElement)
1115 1120
1116 1121 if readUnitConfObj.parentId == None:
1117 1122 readUnitConfObj.parentId = self.id
1118 1123
1119 1124 self.procUnitConfObjDict[readUnitConfObj.getId()] = readUnitConfObj
1120 1125
1121 1126 procUnitElementList = self.projectElement.iter(ProcUnitConf().getElementName())
1122 1127
1123 1128 for procUnitElement in procUnitElementList:
1124 1129 procUnitConfObj = ProcUnitConf()
1125 1130 procUnitConfObj.readXml(procUnitElement)
1126 1131
1127 1132 if procUnitConfObj.parentId == None:
1128 1133 procUnitConfObj.parentId = self.id
1129 1134
1130 1135 self.procUnitConfObjDict[procUnitConfObj.getId()] = procUnitConfObj
1131 1136
1132 1137 self.filename = abs_file
1133 1138
1134 1139 return 1
1135 1140
1136 1141 def printattr(self):
1137 1142
1138 1143 print "Project[%s]: name = %s, description = %s" %(self.id,
1139 1144 self.name,
1140 1145 self.description)
1141 1146
1142 1147 for procUnitConfObj in self.procUnitConfObjDict.values():
1143 1148 procUnitConfObj.printattr()
1144 1149
1145 1150 def createObjects(self):
1146 1151
1147 1152 for procUnitConfObj in self.procUnitConfObjDict.values():
1148 1153 procUnitConfObj.createObjects(self.plotterQueue)
1149 1154
1150 1155 def __connect(self, objIN, thisObj):
1151 1156
1152 1157 thisObj.setInput(objIN.getOutputObj())
1153 1158
1154 1159 def connectObjects(self):
1155 1160
1156 1161 for thisPUConfObj in self.procUnitConfObjDict.values():
1157 1162
1158 1163 inputId = thisPUConfObj.getInputId()
1159 1164
1160 1165 if int(inputId) == 0:
1161 1166 continue
1162 1167
1163 1168 #Get input object
1164 1169 puConfINObj = self.procUnitConfObjDict[inputId]
1165 1170 puObjIN = puConfINObj.getProcUnitObj()
1166 1171
1167 1172 #Get current object
1168 1173 thisPUObj = thisPUConfObj.getProcUnitObj()
1169 1174
1170 1175 self.__connect(puObjIN, thisPUObj)
1171 1176
1172 1177 def __handleError(self, procUnitConfObj, send_email=True):
1173 1178
1174 1179 import socket
1175 1180
1176 1181 err = traceback.format_exception(sys.exc_info()[0],
1177 1182 sys.exc_info()[1],
1178 1183 sys.exc_info()[2])
1179 1184
1180 1185 print "***** Error occurred in %s *****" %(procUnitConfObj.name)
1181 1186 print "***** %s" %err[-1]
1182 1187
1183 1188 message = "".join(err)
1184 1189
1185 1190 sys.stderr.write(message)
1186 1191
1187 1192 if not send_email:
1188 1193 return
1189 1194
1190 1195 subject = "SChain v%s: Error running %s\n" %(schainpy.__version__, procUnitConfObj.name)
1191 1196
1192 1197 subtitle = "%s: %s\n" %(procUnitConfObj.getElementName() ,procUnitConfObj.name)
1193 1198 subtitle += "Hostname: %s\n" %socket.gethostbyname(socket.gethostname())
1194 1199 subtitle += "Working directory: %s\n" %os.path.abspath("./")
1195 1200 subtitle += "Configuration file: %s\n" %self.filename
1196 1201 subtitle += "Time: %s\n" %str(datetime.datetime.now())
1197 1202
1198 1203 readUnitConfObj = self.getReadUnitObj()
1199 1204 if readUnitConfObj:
1200 1205 subtitle += "\nInput parameters:\n"
1201 1206 subtitle += "[Data path = %s]\n" %readUnitConfObj.path
1202 1207 subtitle += "[Data type = %s]\n" %readUnitConfObj.datatype
1203 1208 subtitle += "[Start date = %s]\n" %readUnitConfObj.startDate
1204 1209 subtitle += "[End date = %s]\n" %readUnitConfObj.endDate
1205 1210 subtitle += "[Start time = %s]\n" %readUnitConfObj.startTime
1206 1211 subtitle += "[End time = %s]\n" %readUnitConfObj.endTime
1207 1212
1208 1213 adminObj = schainpy.admin.SchainNotify()
1209 1214 adminObj.sendAlert(message=message,
1210 1215 subject=subject,
1211 1216 subtitle=subtitle,
1212 1217 filename=self.filename)
1213 1218
1214 1219 def isPaused(self):
1215 1220 return 0
1216 1221
1217 1222 def isStopped(self):
1218 1223 return 0
1219 1224
1220 1225 def runController(self):
1221 1226 """
1222 1227 returns 0 when this process has been stopped, 1 otherwise
1223 1228 """
1224 1229
1225 1230 if self.isPaused():
1226 1231 print "Process suspended"
1227 1232
1228 1233 while True:
1229 1234 sleep(0.1)
1230 1235
1231 1236 if not self.isPaused():
1232 1237 break
1233 1238
1234 1239 if self.isStopped():
1235 1240 break
1236 1241
1237 1242 print "Process reinitialized"
1238 1243
1239 1244 if self.isStopped():
1240 1245 print "Process stopped"
1241 1246 return 0
1242 1247
1243 1248 return 1
1244 1249
1245 1250 def setFilename(self, filename):
1246 1251
1247 1252 self.filename = filename
1248 1253
1249 1254 def setPlotterQueue(self, plotter_queue):
1250 1255
1251 1256 raise NotImplementedError, "Use schainpy.controller_api.ControllerThread instead Project class"
1252 1257
1253 1258 def getPlotterQueue(self):
1254 1259
1255 1260 raise NotImplementedError, "Use schainpy.controller_api.ControllerThread instead Project class"
1256 1261
1257 1262 def useExternalPlotter(self):
1258 1263
1259 1264 raise NotImplementedError, "Use schainpy.controller_api.ControllerThread instead Project class"
1260 1265
1261 1266 def run(self):
1262 1267
1263 1268 print
1264 1269 print "*"*60
1265 1270 print " Starting SIGNAL CHAIN PROCESSING v%s " %schainpy.__version__
1266 1271 print "*"*60
1267 1272 print
1268 1273
1269 1274 keyList = self.procUnitConfObjDict.keys()
1270 1275 keyList.sort()
1271 1276
1272 1277 while(True):
1273 1278
1274 1279 is_ok = False
1275 1280
1276 1281 for procKey in keyList:
1277 1282 # print "Running the '%s' process with %s" %(procUnitConfObj.name, procUnitConfObj.id)
1278 1283
1279 1284 procUnitConfObj = self.procUnitConfObjDict[procKey]
1280 1285
1281 1286 try:
1282 1287 sts = procUnitConfObj.run()
1283 1288 is_ok = is_ok or sts
1284 1289 except KeyboardInterrupt:
1285 1290 is_ok = False
1286 1291 break
1287 1292 except ValueError, e:
1288 1293 sleep(0.5)
1289 1294 self.__handleError(procUnitConfObj, send_email=True)
1290 1295 is_ok = False
1291 1296 break
1292 1297 except:
1293 1298 sleep(0.5)
1294 1299 self.__handleError(procUnitConfObj)
1295 1300 is_ok = False
1296 1301 break
1297 1302
1298 1303 #If every process unit finished so end process
1299 1304 if not(is_ok):
1300 1305 # print "Every process unit have finished"
1301 1306 break
1302 1307
1303 1308 if not self.runController():
1304 1309 break
1305 1310
1306 1311 #Closing every process
1307 1312 for procKey in keyList:
1308 1313 procUnitConfObj = self.procUnitConfObjDict[procKey]
1309 1314 procUnitConfObj.close()
1310 1315
1311 1316 print "Process finished"
1312 1317
1313 1318 def start(self):
1314 1319
1315 1320 self.writeXml()
1316 1321 self.createObjects()
1317 1322 self.connectObjects()
1318 1323 self.run()
Show file before | Show file after | Hide comments
@@ -1,377 +1,378
1 1
2 2 import os
3 3 import zmq
4 4 import time
5 5 import numpy
6 6 import datetime
7 7 import numpy as np
8 8 import matplotlib.pyplot as plt
9 9 from mpl_toolkits.axes_grid1 import make_axes_locatable
10 10 from matplotlib.ticker import FuncFormatter, LinearLocator
11 11 from multiprocessing import Process
12 12
13 13 from schainpy.model.proc.jroproc_base import Operation
14 14
15 15 #plt.ion()
16 16
17 17 func = lambda x, pos: ('%s') %(datetime.datetime.utcfromtimestamp(x).strftime('%H:%M'))
18 18
19 19 d1970 = datetime.datetime(1970,1,1)
20 20
21 21 class PlotData(Operation, Process):
22 22
23 23 CODE = 'Figure'
24 24 colormap = 'jet'
25 25 CONFLATE = True
26 26 __MAXNUMX = 80
27 27 __MAXNUMY = 80
28 28 __missing = 1E30
29 29
30 30 def __init__(self, **kwargs):
31 31
32 Operation.__init__(self, **kwargs)
32 Operation.__init__(self, plot=True, **kwargs)
33 33 Process.__init__(self)
34 self.kwargs['code'] = self.CODE
34 35 self.mp = False
35 36 self.dataOut = None
36 37 self.isConfig = False
37 38 self.figure = None
38 39 self.axes = []
39 40 self.localtime = kwargs.pop('localtime', True)
40 41 self.show = kwargs.get('show', True)
41 42 self.save = kwargs.get('save', False)
42 43 self.colormap = kwargs.get('colormap', self.colormap)
43 44 self.showprofile = kwargs.get('showprofile', False)
44 45 self.title = kwargs.get('wintitle', '')
45 46 self.xaxis = kwargs.get('xaxis', 'time')
46 47 self.zmin = kwargs.get('zmin', None)
47 48 self.zmax = kwargs.get('zmax', None)
48 49 self.xmin = kwargs.get('xmin', None)
49 50 self.xmax = kwargs.get('xmax', None)
50 51 self.xrange = kwargs.get('xrange', 24)
51 52 self.ymin = kwargs.get('ymin', None)
52 53 self.ymax = kwargs.get('ymax', None)
53 54 self.throttle_value = 1
54 55 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
55 56
56 57 if x_buffer.shape[0] < 2:
57 58 return x_buffer, y_buffer, z_buffer
58 59
59 60 deltas = x_buffer[1:] - x_buffer[0:-1]
60 61 x_median = np.median(deltas)
61 62
62 63 index = np.where(deltas > 5*x_median)
63 64
64 65 if len(index[0]) != 0:
65 66 z_buffer[::, index[0], ::] = self.__missing
66 67 z_buffer = np.ma.masked_inside(z_buffer,
67 68 0.99*self.__missing,
68 69 1.01*self.__missing)
69 70
70 71 return x_buffer, y_buffer, z_buffer
71 72
72 73 def decimate(self):
73 74
74 75 # dx = int(len(self.x)/self.__MAXNUMX) + 1
75 76 dy = int(len(self.y)/self.__MAXNUMY) + 1
76 77
77 78 # x = self.x[::dx]
78 79 x = self.x
79 80 y = self.y[::dy]
80 81 z = self.z[::, ::, ::dy]
81 82
82 83 return x, y, z
83 84
84 85 def __plot(self):
85 86
86 87 print 'plotting...{}'.format(self.CODE)
87 88
88 89 self.plot()
89 90 self.figure.suptitle('{} {} - Date:{}'.format(self.title, self.CODE.upper(),
90 91 datetime.datetime.utcfromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')))
91 92
92 93 if self.save:
93 94 figname = os.path.join(self.save, '{}_{}.png'.format(self.CODE,
94 95 datetime.datetime.utcfromtimestamp(self.times[0]).strftime('%y%m%d_%H%M%S')))
95 96 print 'Saving figure: {}'.format(figname)
96 97 self.figure.savefig(figname)
97 98
98 99 self.figure.canvas.draw()
99 100
100 101 def plot(self):
101 102
102 103 print 'plotting...{}'.format(self.CODE.upper())
103 104 return
104 105
105 106 def run(self):
106 107
107 108 print '[Starting] {}'.format(self.name)
108 109 context = zmq.Context()
109 110 receiver = context.socket(zmq.SUB)
110 111 receiver.setsockopt(zmq.SUBSCRIBE, '')
111 112 receiver.setsockopt(zmq.CONFLATE, self.CONFLATE)
112 113 receiver.connect("ipc:///tmp/zmq.plots")
113 114
114 115 while True:
115 116 try:
116 117 #if True:
117 118 self.data = receiver.recv_pyobj(flags=zmq.NOBLOCK)
118 119 self.dataOut = self.data['dataOut']
119 120 self.times = self.data['times']
120 121 self.times.sort()
121 122 self.throttle_value = self.data['throttle']
122 123 self.min_time = self.times[0]
123 124 self.max_time = self.times[-1]
124 125
125 126 if self.isConfig is False:
126 127 self.setup()
127 128 self.isConfig = True
128 129 self.__plot()
129 130
130 131 if self.data['ENDED'] is True:
131 132 # self.__plot()
132 133 self.isConfig = False
133 134
134 135 except zmq.Again as e:
135 136 print 'Waiting for data...'
136 137 plt.pause(self.throttle_value)
137 138 # time.sleep(3)
138 139
139 140 def close(self):
140 141 if self.dataOut:
141 142 self._plot()
142 143
143 144
144 145 class PlotSpectraData(PlotData):
145 146
146 147 CODE = 'spc'
147 148 colormap = 'jro'
148 149 CONFLATE = False
149 150 def setup(self):
150 151
151 152 ncolspan = 1
152 153 colspan = 1
153 154 self.ncols = int(numpy.sqrt(self.dataOut.nChannels)+0.9)
154 155 self.nrows = int(self.dataOut.nChannels*1./self.ncols + 0.9)
155 156 self.width = 3.6*self.ncols
156 157 self.height = 3.2*self.nrows
157 158 if self.showprofile:
158 159 ncolspan = 3
159 160 colspan = 2
160 161 self.width += 1.2*self.ncols
161 162
162 163 self.ylabel = 'Range [Km]'
163 164 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
164 165
165 166 if self.figure is None:
166 167 self.figure = plt.figure(figsize=(self.width, self.height),
167 168 edgecolor='k',
168 169 facecolor='w')
169 170 else:
170 171 self.figure.clf()
171 172
172 173 n = 0
173 174 for y in range(self.nrows):
174 175 for x in range(self.ncols):
175 176 if n >= self.dataOut.nChannels:
176 177 break
177 178 ax = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan), 1, colspan)
178 179 if self.showprofile:
179 180 ax.ax_profile = plt.subplot2grid((self.nrows, self.ncols*ncolspan), (y, x*ncolspan+colspan), 1, 1)
180 181
181 182 ax.firsttime = True
182 183 self.axes.append(ax)
183 184 n += 1
184 185
185 186 self.figure.subplots_adjust(wspace=0.9, hspace=0.5)
186 187 self.figure.show()
187 188
188 189 def plot(self):
189 190
190 191 if self.xaxis == "frequency":
191 192 x = self.dataOut.getFreqRange(1)/1000.
192 193 xlabel = "Frequency (kHz)"
193 194 elif self.xaxis == "time":
194 195 x = self.dataOut.getAcfRange(1)
195 196 xlabel = "Time (ms)"
196 197 else:
197 198 x = self.dataOut.getVelRange(1)
198 199 xlabel = "Velocity (m/s)"
199 200
200 201 y = self.dataOut.getHeiRange()
201 202 z = self.data[self.CODE]
202 203
203 204 for n, ax in enumerate(self.axes):
204 205
205 206 if ax.firsttime:
206 207 self.xmax = self.xmax if self.xmax else np.nanmax(x)
207 208 self.xmin = self.xmin if self.xmin else -self.xmax
208 209 self.ymin = self.ymin if self.ymin else np.nanmin(y)
209 210 self.ymax = self.ymax if self.ymax else np.nanmax(y)
210 211 self.zmin = self.zmin if self.zmin else np.nanmin(z)
211 212 self.zmax = self.zmax if self.zmax else np.nanmax(z)
212 213 ax.plot = ax.pcolormesh(x, y, z[n].T,
213 214 vmin=self.zmin,
214 215 vmax=self.zmax,
215 216 cmap=plt.get_cmap(self.colormap)
216 217 )
217 218 divider = make_axes_locatable(ax)
218 219 cax = divider.new_horizontal(size='3%', pad=0.05)
219 220 self.figure.add_axes(cax)
220 221 plt.colorbar(ax.plot, cax)
221 222
222 223 ax.set_xlim(self.xmin, self.xmax)
223 224 ax.set_ylim(self.ymin, self.ymax)
224 225
225 226 ax.xaxis.set_major_locator(LinearLocator(5))
226 227 #ax.yaxis.set_major_locator(LinearLocator(4))
227 228
228 229 ax.set_ylabel(self.ylabel)
229 230 ax.set_xlabel(xlabel)
230 231
231 232 ax.firsttime = False
232 233
233 234 if self.showprofile:
234 235 ax.plot_profile= ax.ax_profile.plot(self.data['rti'][self.max_time][n], y)[0]
235 236 ax.ax_profile.set_xlim(self.zmin, self.zmax)
236 237 ax.ax_profile.set_ylim(self.ymin, self.ymax)
237 238 ax.ax_profile.set_xlabel('dB')
238 239 ax.ax_profile.grid(b=True, axis='x')
239 240 ax.plot_noise = ax.ax_profile.plot(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y,
240 241 color="k", linestyle="dashed", lw=2)[0]
241 242 [tick.set_visible(False) for tick in ax.ax_profile.get_yticklabels()]
242 243 else:
243 244 ax.plot.set_array(z[n].T.ravel())
244 245 if self.showprofile:
245 246 ax.plot_profile.set_data(self.data['rti'][self.max_time][n], y)
246 247 ax.plot_noise.set_data(numpy.repeat(self.data['noise'][self.max_time][n], len(y)), y)
247 248
248 249 ax.set_title('{} - Noise: {:.2f} dB'.format(self.titles[n], self.data['noise'][self.max_time][n]),
249 250 size=8)
250 251
251 252 class PlotRTIData(PlotData):
252 253
253 254 CODE = 'rti'
254 255 colormap = 'jro'
255 256
256 257 def setup(self):
257 258 self.ncols = 1
258 259 self.nrows = self.dataOut.nChannels
259 260 self.width = 10
260 261 self.height = 2.2*self.nrows
261 262 self.ylabel = 'Range [Km]'
262 263 self.titles = ['Channel {}'.format(x) for x in self.dataOut.channelList]
263 264
264 265 if self.figure is None:
265 266 self.figure = plt.figure(figsize=(self.width, self.height),
266 267 edgecolor='k',
267 268 facecolor='w')
268 269 else:
269 270 self.figure.clf()
270 271 self.axes = []
271 272
272 273 for n in range(self.nrows):
273 274 ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
274 275 ax.firsttime = True
275 276 self.axes.append(ax)
276 277 self.figure.subplots_adjust(hspace=0.5)
277 278 self.figure.show()
278 279
279 280 def plot(self):
280 281
281 282 self.x = np.array(self.times)
282 283 self.y = self.dataOut.getHeiRange()
283 284 self.z = []
284 285
285 286 for ch in range(self.nrows):
286 287 self.z.append([self.data[self.CODE][t][ch] for t in self.times])
287 288
288 289 self.z = np.array(self.z)
289 290 for n, ax in enumerate(self.axes):
290 291
291 292 x, y, z = self.fill_gaps(*self.decimate())
292 293 xmin = self.min_time
293 294 xmax = xmin+self.xrange*60*60
294 295 if ax.firsttime:
295 296 self.ymin = self.ymin if self.ymin else np.nanmin(self.y)
296 297 self.ymax = self.ymax if self.ymax else np.nanmax(self.y)
297 298 self.zmin = self.zmin if self.zmin else np.nanmin(self.z)
298 299 self.zmax = self.zmax if self.zmax else np.nanmax(self.z)
299 300 plot = ax.pcolormesh(x, y, z[n].T,
300 301 vmin=self.zmin,
301 302 vmax=self.zmax,
302 303 cmap=plt.get_cmap(self.colormap)
303 304 )
304 305 divider = make_axes_locatable(ax)
305 306 cax = divider.new_horizontal(size='2%', pad=0.05)
306 307 self.figure.add_axes(cax)
307 308 plt.colorbar(plot, cax)
308 309 ax.set_ylim(self.ymin, self.ymax)
309 310 if self.xaxis == 'time':
310 311 ax.xaxis.set_major_formatter(FuncFormatter(func))
311 312 ax.xaxis.set_major_locator(LinearLocator(6))
312 313
313 314 # ax.yaxis.set_major_locator(LinearLocator(4))
314 315
315 316 ax.set_ylabel(self.ylabel)
316 317
317 318 # if self.xmin is None:
318 319 # xmin = self.min_time
319 320 # else:
320 321 # xmin = (datetime.datetime.combine(self.dataOut.datatime.date(),
321 322 # datetime.time(self.xmin, 0, 0))-d1970).total_seconds()
322 323
323 324 ax.set_xlim(xmin, xmax)
324 325 ax.firsttime = False
325 326 else:
326 327 ax.collections.remove(ax.collections[0])
327 328 ax.set_xlim(xmin, xmax)
328 329 plot = ax.pcolormesh(x, y, z[n].T,
329 330 vmin=self.zmin,
330 331 vmax=self.zmax,
331 332 cmap=plt.get_cmap(self.colormap)
332 333 )
333 334 ax.set_title('{} {}'.format(self.titles[n],
334 335 datetime.datetime.utcfromtimestamp(self.max_time).strftime('%y/%m/%d %H:%M:%S')),
335 336 size=8)
336 337
337 338
338 339 class PlotCOHData(PlotRTIData):
339 340
340 341 CODE = 'coh'
341 342
342 343 def setup(self):
343 344
344 345 self.ncols = 1
345 346 self.nrows = self.dataOut.nPairs
346 347 self.width = 10
347 348 self.height = 2.2*self.nrows
348 349 self.ylabel = 'Range [Km]'
349 350 self.titles = ['Channels {}'.format(x) for x in self.dataOut.pairsList]
350 351
351 352 if self.figure is None:
352 353 self.figure = plt.figure(figsize=(self.width, self.height),
353 354 edgecolor='k',
354 355 facecolor='w')
355 356 else:
356 357 self.figure.clf()
357 358
358 359 for n in range(self.nrows):
359 360 ax = self.figure.add_subplot(self.nrows, self.ncols, n+1)
360 361 ax.firsttime = True
361 362 self.axes.append(ax)
362 363
363 364 self.figure.subplots_adjust(hspace=0.5)
364 365 self.figure.show()
365 366
366 367 class PlotSNRData(PlotRTIData):
367 368
368 369 CODE = 'snr'
369 370
370 371 class PlotDOPData(PlotRTIData):
371 372 CODE = 'dop'
372 373 colormap = 'jet'
373 374
374 375 class PlotPHASEData(PlotCOHData):
375 376
376 377 CODE = 'phase'
377 378 colormap = 'seismic'
Show file before | Show file after | Hide comments
@@ -1,299 +1,315
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: jroproc_base.py 1 2012-11-12 18:56:07Z murco $
5 5 '''
6 6
7 7 class ProcessingUnit(object):
8 8
9 9 """
10 10 Esta es la clase base para el procesamiento de datos.
11 11
12 12 Contiene el metodo "call" para llamar operaciones. Las operaciones pueden ser:
13 13 - Metodos internos (callMethod)
14 14 - Objetos del tipo Operation (callObject). Antes de ser llamados, estos objetos
15 15 tienen que ser agreagados con el metodo "add".
16 16
17 17 """
18 18 # objeto de datos de entrada (Voltage, Spectra o Correlation)
19 19 dataIn = None
20 20 dataInList = []
21 21
22 22 # objeto de datos de entrada (Voltage, Spectra o Correlation)
23 23 dataOut = None
24 24
25 25 operations2RunDict = None
26 26
27 27 isConfig = False
28 28
29 29
30 30 def __init__(self, *args, **kwargs):
31 31
32 32 self.dataIn = None
33 33 self.dataInList = []
34 34
35 35 self.dataOut = None
36 36
37 37 self.operations2RunDict = {}
38 self.operationKwargs = {}
38 39
39 40 self.isConfig = False
40 41
41 42 self.args = args
42 43 self.kwargs = kwargs
43 44
45 def addOperationKwargs(self, objId, **kwargs):
46 '''
47 '''
48
49 self.operationKwargs[objId] = kwargs
50
51
44 52 def addOperation(self, opObj, objId):
45 53
46 54 """
47 55 Agrega un objeto del tipo "Operation" (opObj) a la lista de objetos "self.objectList" y retorna el
48 56 identificador asociado a este objeto.
49 57
50 58 Input:
51 59
52 60 object : objeto de la clase "Operation"
53 61
54 62 Return:
55 63
56 64 objId : identificador del objeto, necesario para ejecutar la operacion
57 65 """
58 66
59 67 self.operations2RunDict[objId] = opObj
60 68
61 69 return objId
62 70
63 71 def getOperationObj(self, objId):
64 72
65 73 if objId not in self.operations2RunDict.keys():
66 74 return None
67 75
68 76 return self.operations2RunDict[objId]
69 77
70 78 def operation(self, **kwargs):
71 79
72 80 """
73 81 Operacion directa sobre la data (dataOut.data). Es necesario actualizar los valores de los
74 82 atributos del objeto dataOut
75 83
76 84 Input:
77 85
78 86 **kwargs : Diccionario de argumentos de la funcion a ejecutar
79 87 """
80 88
81 89 raise NotImplementedError
82 90
83 def callMethod(self, name, **kwargs):
91 def callMethod(self, name, opId):
84 92
85 93 """
86 94 Ejecuta el metodo con el nombre "name" y con argumentos **kwargs de la propia clase.
87 95
88 96 Input:
89 97 name : nombre del metodo a ejecutar
90 98
91 99 **kwargs : diccionario con los nombres y valores de la funcion a ejecutar.
92 100
93 101 """
94 102
95 103 #Checking the inputs
96 104 if name == 'run':
97 105
98 106 if not self.checkInputs():
99 107 self.dataOut.flagNoData = True
100 108 return False
101 109 else:
102 110 #Si no es un metodo RUN la entrada es la misma dataOut (interna)
103 if self.dataOut.isEmpty():
111 if self.dataOut is not None and self.dataOut.isEmpty():
104 112 return False
105 113
106 114 #Getting the pointer to method
107 115 methodToCall = getattr(self, name)
108 116
109 117 #Executing the self method
110 118
111 119 if hasattr(self, 'mp'):
112 if self.mp is False:
113 self.mp = True
114 self.start()
120 if name=='run':
121 if self.mp is False:
122 self.mp = True
123 self.start()
124 else:
125 methodToCall(**self.operationKwargs[opId])
115 126 else:
116 methodToCall(**kwargs)
127 if name=='run':
128 methodToCall(**self.kwargs)
129 else:
130 methodToCall(**self.operationKwargs[opId])
117 131
118 132 if self.dataOut is None:
119 133 return False
120 134
121 135 if self.dataOut.isEmpty():
122 136 return False
123 137
124 138 return True
125 139
126 140 def callObject(self, objId):
127 141
128 142 """
129 143 Ejecuta la operacion asociada al identificador del objeto "objId"
130 144
131 145 Input:
132 146
133 147 objId : identificador del objeto a ejecutar
134 148
135 149 **kwargs : diccionario con los nombres y valores de la funcion a ejecutar.
136 150
137 151 Return:
138 152
139 153 None
140 154 """
141 155
142 156 if self.dataOut is not None and self.dataOut.isEmpty():
143 157 return False
144 158
145 159 externalProcObj = self.operations2RunDict[objId]
146 160
147 161 if hasattr(externalProcObj, 'mp'):
148 162 if externalProcObj.mp is False:
163 self.operationKwargs[objId] = externalProcObj.kwargs
149 164 externalProcObj.mp = True
150 165 externalProcObj.start()
151 166 else:
152 167 externalProcObj.run(self.dataOut, **externalProcObj.kwargs)
168 self.operationKwargs[objId] = externalProcObj.kwargs
153 169
154 170 return True
155 171
156 172 def call(self, opType, opName=None, opId=None):
157 173
158 174 """
159 175 Return True si ejecuta la operacion interna nombrada "opName" o la operacion externa
160 176 identificada con el id "opId"; con los argumentos "**kwargs".
161 177
162 178 False si la operacion no se ha ejecutado.
163 179
164 180 Input:
165 181
166 182 opType : Puede ser "self" o "external"
167 183
168 184 Depende del tipo de operacion para llamar a:callMethod or callObject:
169 185
170 186 1. If opType = "self": Llama a un metodo propio de esta clase:
171 187
172 188 name_method = getattr(self, name)
173 189 name_method(**kwargs)
174 190
175 191
176 192 2. If opType = "other" o"external": Llama al metodo "run()" de una instancia de la
177 193 clase "Operation" o de un derivado de ella:
178 194
179 195 instanceName = self.operationList[opId]
180 196 instanceName.run(**kwargs)
181 197
182 198 opName : Si la operacion es interna (opType = 'self'), entonces el "opName" sera
183 199 usada para llamar a un metodo interno de la clase Processing
184 200
185 201 opId : Si la operacion es externa (opType = 'other' o 'external), entonces el
186 202 "opId" sera usada para llamar al metodo "run" de la clase Operation
187 203 registrada anteriormente con ese Id
188 204
189 205 Exception:
190 206 Este objeto de tipo Operation debe de haber sido agregado antes con el metodo:
191 207 "addOperation" e identificado con el valor "opId" = el id de la operacion.
192 208 De lo contrario retornara un error del tipo ValueError
193 209
194 210 """
195 211
196 212 if opType == 'self':
197 213
198 214 if not opName:
199 215 raise ValueError, "opName parameter should be defined"
200 216
201 sts = self.callMethod(opName, **self.kwargs)
217 sts = self.callMethod(opName, opId)
202 218
203 219 elif opType == 'other' or opType == 'external' or opType == 'plotter':
204 220
205 221 if not opId:
206 222 raise ValueError, "opId parameter should be defined"
207 223
208 224 if opId not in self.operations2RunDict.keys():
209 225 raise ValueError, "Any operation with id=%s has been added" %str(opId)
210 226
211 227 sts = self.callObject(opId)
212 228
213 229 else:
214 230 raise ValueError, "opType should be 'self', 'external' or 'plotter'; and not '%s'" %opType
215 231
216 232 return sts
217 233
218 234 def setInput(self, dataIn):
219 235
220 236 self.dataIn = dataIn
221 237 self.dataInList.append(dataIn)
222 238
223 239 def getOutputObj(self):
224 240
225 241 return self.dataOut
226 242
227 243 def checkInputs(self):
228 244
229 245 for thisDataIn in self.dataInList:
230 246
231 247 if thisDataIn.isEmpty():
232 248 return False
233 249
234 250 return True
235 251
236 252 def setup(self):
237 253
238 254 raise NotImplementedError
239 255
240 256 def run(self):
241 257
242 258 raise NotImplementedError
243 259
244 260 def close(self):
245 261 #Close every thread, queue or any other object here is it is neccesary.
246 262 return
247 263
248 264 class Operation(object):
249 265
250 266 """
251 267 Clase base para definir las operaciones adicionales que se pueden agregar a la clase ProcessingUnit
252 268 y necesiten acumular informacion previa de los datos a procesar. De preferencia usar un buffer de
253 269 acumulacion dentro de esta clase
254 270
255 271 Ejemplo: Integraciones coherentes, necesita la informacion previa de los n perfiles anteriores (bufffer)
256 272
257 273 """
258 274
259 275 __buffer = None
260 276 isConfig = False
261 277
262 278 def __init__(self, **kwargs):
263 279
264 280 self.__buffer = None
265 281 self.isConfig = False
266 282 self.kwargs = kwargs
267 283
268 284 def setup(self):
269 285
270 286 self.isConfig = True
271 287
272 288 raise NotImplementedError
273 289
274 290 def run(self, dataIn, **kwargs):
275 291
276 292 """
277 293 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los
278 294 atributos del objeto dataIn.
279 295
280 296 Input:
281 297
282 298 dataIn : objeto del tipo JROData
283 299
284 300 Return:
285 301
286 302 None
287 303
288 304 Affected:
289 305 __buffer : buffer de recepcion de datos.
290 306
291 307 """
292 308 if not self.isConfig:
293 309 self.setup(**kwargs)
294 310
295 311 raise NotImplementedError
296 312
297 313 def close(self):
298 314
299 315 pass
Show file before | Show file after | Hide comments
@@ -1,2746 +1,2741
1 1 import numpy
2 2 import math
3 3 from scipy import optimize, interpolate, signal, stats, ndimage
4 4 import re
5 5 import datetime
6 6 import copy
7 7 import sys
8 8 import importlib
9 9 import itertools
10 10
11 11 from jroproc_base import ProcessingUnit, Operation
12 12 from schainpy.model.data.jrodata import Parameters, hildebrand_sekhon
13 13
14 14
15 15 class ParametersProc(ProcessingUnit):
16 16
17 17 nSeconds = None
18 18
19 19 def __init__(self):
20 20 ProcessingUnit.__init__(self)
21 21
22 22 # self.objectDict = {}
23 23 self.buffer = None
24 24 self.firstdatatime = None
25 25 self.profIndex = 0
26 26 self.dataOut = Parameters()
27 27
28 28 def __updateObjFromInput(self):
29 29
30 30 self.dataOut.inputUnit = self.dataIn.type
31 31
32 32 self.dataOut.timeZone = self.dataIn.timeZone
33 33 self.dataOut.dstFlag = self.dataIn.dstFlag
34 34 self.dataOut.errorCount = self.dataIn.errorCount
35 35 self.dataOut.useLocalTime = self.dataIn.useLocalTime
36 36
37 37 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
38 38 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 self.dataOut.heightList = self.dataIn.heightList
41 41 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
42 42 # self.dataOut.nHeights = self.dataIn.nHeights
43 43 # self.dataOut.nChannels = self.dataIn.nChannels
44 44 self.dataOut.nBaud = self.dataIn.nBaud
45 45 self.dataOut.nCode = self.dataIn.nCode
46 46 self.dataOut.code = self.dataIn.code
47 47 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
48 48 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
49 49 # self.dataOut.utctime = self.firstdatatime
50 50 self.dataOut.utctime = self.dataIn.utctime
51 51 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
52 52 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
53 53 self.dataOut.nCohInt = self.dataIn.nCohInt
54 54 # self.dataOut.nIncohInt = 1
55 55 self.dataOut.ippSeconds = self.dataIn.ippSeconds
56 56 # self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
57 57 self.dataOut.timeInterval = self.dataIn.timeInterval
58 58 self.dataOut.heightList = self.dataIn.getHeiRange()
59 59 self.dataOut.frequency = self.dataIn.frequency
60 60 self.dataOut.noise = self.dataIn.noise
61 61
62 62 def run(self):
63 63
64 64 #---------------------- Voltage Data ---------------------------
65 65
66 66 if self.dataIn.type == "Voltage":
67 67
68 68 self.__updateObjFromInput()
69 69 self.dataOut.data_pre = self.dataIn.data.copy()
70 70 self.dataOut.flagNoData = False
71 71 self.dataOut.utctimeInit = self.dataIn.utctime
72 72 self.dataOut.paramInterval = self.dataIn.nProfiles*self.dataIn.nCohInt*self.dataIn.ippSeconds
73 73 return
74 74
75 75 #---------------------- Spectra Data ---------------------------
76 76
77 77 if self.dataIn.type == "Spectra":
78 78
79 79 self.dataOut.data_pre = (self.dataIn.data_spc,self.dataIn.data_cspc)
80 80 self.dataOut.abscissaList = self.dataIn.getVelRange(1)
81 81 # self.dataOut.noise = self.dataIn.getNoise()
82 82 self.dataOut.normFactor = self.dataIn.normFactor
83 83 self.dataOut.groupList = self.dataIn.pairsList
84 84 self.dataOut.flagNoData = False
85 85
86 86 #---------------------- Correlation Data ---------------------------
87 87
88 88 if self.dataIn.type == "Correlation":
89 89 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.dataIn.splitFunctions()
90 90
91 91 self.dataOut.data_pre = (self.dataIn.data_cf[acf_ind,:], self.dataIn.data_cf[ccf_ind,:,:])
92 92 self.dataOut.normFactor = (self.dataIn.normFactor[acf_ind,:], self.dataIn.normFactor[ccf_ind,:])
93 93 self.dataOut.groupList = (acf_pairs, ccf_pairs)
94 94
95 95 self.dataOut.abscissaList = self.dataIn.lagRange
96 96 self.dataOut.noise = self.dataIn.noise
97 97 self.dataOut.data_SNR = self.dataIn.SNR
98 98 self.dataOut.flagNoData = False
99 99 self.dataOut.nAvg = self.dataIn.nAvg
100 100
101 101 #---------------------- Parameters Data ---------------------------
102 102
103 103 if self.dataIn.type == "Parameters":
104 104 self.dataOut.copy(self.dataIn)
105 105 self.dataOut.utctimeInit = self.dataIn.utctime
106 106 self.dataOut.flagNoData = False
107 107
108 108 return True
109 109
110 110 self.__updateObjFromInput()
111 111 self.dataOut.utctimeInit = self.dataIn.utctime
112 112 self.dataOut.paramInterval = self.dataIn.timeInterval
113 113
114 114 return
115 115
116 116 class SpectralMoments(Operation):
117 117
118 118 '''
119 119 Function SpectralMoments()
120 120
121 121 Calculates moments (power, mean, standard deviation) and SNR of the signal
122 122
123 123 Type of dataIn: Spectra
124 124
125 125 Configuration Parameters:
126 126
127 127 dirCosx : Cosine director in X axis
128 128 dirCosy : Cosine director in Y axis
129 129
130 130 elevation :
131 131 azimuth :
132 132
133 133 Input:
134 134 channelList : simple channel list to select e.g. [2,3,7]
135 135 self.dataOut.data_pre : Spectral data
136 136 self.dataOut.abscissaList : List of frequencies
137 137 self.dataOut.noise : Noise level per channel
138 138
139 139 Affected:
140 140 self.dataOut.data_param : Parameters per channel
141 141 self.dataOut.data_SNR : SNR per channel
142 142
143 143 '''
144 144
145 145 def run(self, dataOut):
146 146
147 147 dataOut.data_pre = dataOut.data_pre[0]
148 148 data = dataOut.data_pre
149 149 absc = dataOut.abscissaList[:-1]
150 150 noise = dataOut.noise
151 151 nChannel = data.shape[0]
152 152 data_param = numpy.zeros((nChannel, 4, data.shape[2]))
153 153
154 154 for ind in range(nChannel):
155 155 data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
156 156
157 157 dataOut.data_param = data_param[:,1:,:]
158 158 dataOut.data_SNR = data_param[:,0]
159 159 dataOut.data_DOP = data_param[:,1]
160 160 return
161 161
162 162 def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
163 163
164 164 if (nicoh is None): nicoh = 1
165 165 if (graph is None): graph = 0
166 166 if (smooth is None): smooth = 0
167 167 elif (self.smooth < 3): smooth = 0
168 168
169 169 if (type1 is None): type1 = 0
170 170 if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
171 171 if (snrth is None): snrth = -3
172 172 if (dc is None): dc = 0
173 173 if (aliasing is None): aliasing = 0
174 174 if (oldfd is None): oldfd = 0
175 175 if (wwauto is None): wwauto = 0
176 176
177 177 if (n0 < 1.e-20): n0 = 1.e-20
178 178
179 179 freq = oldfreq
180 180 vec_power = numpy.zeros(oldspec.shape[1])
181 181 vec_fd = numpy.zeros(oldspec.shape[1])
182 182 vec_w = numpy.zeros(oldspec.shape[1])
183 183 vec_snr = numpy.zeros(oldspec.shape[1])
184 184
185 185 for ind in range(oldspec.shape[1]):
186 186
187 187 spec = oldspec[:,ind]
188 188 aux = spec*fwindow
189 189 max_spec = aux.max()
190 190 m = list(aux).index(max_spec)
191 191
192 192 #Smooth
193 193 if (smooth == 0): spec2 = spec
194 194 else: spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
195 195
196 196 # Calculo de Momentos
197 197 bb = spec2[range(m,spec2.size)]
198 198 bb = (bb<n0).nonzero()
199 199 bb = bb[0]
200 200
201 201 ss = spec2[range(0,m + 1)]
202 202 ss = (ss<n0).nonzero()
203 203 ss = ss[0]
204 204
205 205 if (bb.size == 0):
206 206 bb0 = spec.size - 1 - m
207 207 else:
208 208 bb0 = bb[0] - 1
209 209 if (bb0 < 0):
210 210 bb0 = 0
211 211
212 212 if (ss.size == 0): ss1 = 1
213 213 else: ss1 = max(ss) + 1
214 214
215 215 if (ss1 > m): ss1 = m
216 216
217 217 valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
218 218 power = ((spec2[valid] - n0)*fwindow[valid]).sum()
219 219 fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
220 220 w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
221 221 snr = (spec2.mean()-n0)/n0
222 222
223 223 if (snr < 1.e-20) :
224 224 snr = 1.e-20
225 225
226 226 vec_power[ind] = power
227 227 vec_fd[ind] = fd
228 228 vec_w[ind] = w
229 229 vec_snr[ind] = snr
230 230
231 231 moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
232 232 return moments
233 233
234 234 #------------------ Get SA Parameters --------------------------
235 235
236 236 def GetSAParameters(self):
237 237 #SA en frecuencia
238 238 pairslist = self.dataOut.groupList
239 239 num_pairs = len(pairslist)
240 240
241 241 vel = self.dataOut.abscissaList
242 242 spectra = self.dataOut.data_pre
243 243 cspectra = self.dataIn.data_cspc
244 244 delta_v = vel[1] - vel[0]
245 245
246 246 #Calculating the power spectrum
247 247 spc_pow = numpy.sum(spectra, 3)*delta_v
248 248 #Normalizing Spectra
249 249 norm_spectra = spectra/spc_pow
250 250 #Calculating the norm_spectra at peak
251 251 max_spectra = numpy.max(norm_spectra, 3)
252 252
253 253 #Normalizing Cross Spectra
254 254 norm_cspectra = numpy.zeros(cspectra.shape)
255 255
256 256 for i in range(num_chan):
257 257 norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
258 258
259 259 max_cspectra = numpy.max(norm_cspectra,2)
260 260 max_cspectra_index = numpy.argmax(norm_cspectra, 2)
261 261
262 262 for i in range(num_pairs):
263 263 cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
264 264 #------------------- Get Lags ----------------------------------
265 265
266 266 class SALags(Operation):
267 267 '''
268 268 Function GetMoments()
269 269
270 270 Input:
271 271 self.dataOut.data_pre
272 272 self.dataOut.abscissaList
273 273 self.dataOut.noise
274 274 self.dataOut.normFactor
275 275 self.dataOut.data_SNR
276 276 self.dataOut.groupList
277 277 self.dataOut.nChannels
278 278
279 279 Affected:
280 280 self.dataOut.data_param
281 281
282 282 '''
283 283 def run(self, dataOut):
284 284 data_acf = dataOut.data_pre[0]
285 285 data_ccf = dataOut.data_pre[1]
286 286 normFactor_acf = dataOut.normFactor[0]
287 287 normFactor_ccf = dataOut.normFactor[1]
288 288 pairs_acf = dataOut.groupList[0]
289 289 pairs_ccf = dataOut.groupList[1]
290 290
291 291 nHeights = dataOut.nHeights
292 292 absc = dataOut.abscissaList
293 293 noise = dataOut.noise
294 294 SNR = dataOut.data_SNR
295 295 nChannels = dataOut.nChannels
296 296 # pairsList = dataOut.groupList
297 297 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
298 298
299 299 for l in range(len(pairs_acf)):
300 300 data_acf[l,:,:] = data_acf[l,:,:]/normFactor_acf[l,:]
301 301
302 302 for l in range(len(pairs_ccf)):
303 303 data_ccf[l,:,:] = data_ccf[l,:,:]/normFactor_ccf[l,:]
304 304
305 305 dataOut.data_param = numpy.zeros((len(pairs_ccf)*2 + 1, nHeights))
306 306 dataOut.data_param[:-1,:] = self.__calculateTaus(data_acf, data_ccf, absc)
307 307 dataOut.data_param[-1,:] = self.__calculateLag1Phase(data_acf, absc)
308 308 return
309 309
310 310 # def __getPairsAutoCorr(self, pairsList, nChannels):
311 311 #
312 312 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
313 313 #
314 314 # for l in range(len(pairsList)):
315 315 # firstChannel = pairsList[l][0]
316 316 # secondChannel = pairsList[l][1]
317 317 #
318 318 # #Obteniendo pares de Autocorrelacion
319 319 # if firstChannel == secondChannel:
320 320 # pairsAutoCorr[firstChannel] = int(l)
321 321 #
322 322 # pairsAutoCorr = pairsAutoCorr.astype(int)
323 323 #
324 324 # pairsCrossCorr = range(len(pairsList))
325 325 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
326 326 #
327 327 # return pairsAutoCorr, pairsCrossCorr
328 328
329 329 def __calculateTaus(self, data_acf, data_ccf, lagRange):
330 330
331 331 lag0 = data_acf.shape[1]/2
332 332 #Funcion de Autocorrelacion
333 333 mean_acf = stats.nanmean(data_acf, axis = 0)
334 334
335 335 #Obtencion Indice de TauCross
336 336 ind_ccf = data_ccf.argmax(axis = 1)
337 337 #Obtencion Indice de TauAuto
338 338 ind_acf = numpy.zeros(ind_ccf.shape,dtype = 'int')
339 339 ccf_lag0 = data_ccf[:,lag0,:]
340 340
341 341 for i in range(ccf_lag0.shape[0]):
342 342 ind_acf[i,:] = numpy.abs(mean_acf - ccf_lag0[i,:]).argmin(axis = 0)
343 343
344 344 #Obtencion de TauCross y TauAuto
345 345 tau_ccf = lagRange[ind_ccf]
346 346 tau_acf = lagRange[ind_acf]
347 347
348 348 Nan1, Nan2 = numpy.where(tau_ccf == lagRange[0])
349 349
350 350 tau_ccf[Nan1,Nan2] = numpy.nan
351 351 tau_acf[Nan1,Nan2] = numpy.nan
352 352 tau = numpy.vstack((tau_ccf,tau_acf))
353 353
354 354 return tau
355 355
356 356 def __calculateLag1Phase(self, data, lagTRange):
357 357 data1 = stats.nanmean(data, axis = 0)
358 358 lag1 = numpy.where(lagTRange == 0)[0][0] + 1
359 359
360 360 phase = numpy.angle(data1[lag1,:])
361 361
362 362 return phase
363 363
364 364 class SpectralFitting(Operation):
365 365 '''
366 366 Function GetMoments()
367 367
368 368 Input:
369 369 Output:
370 370 Variables modified:
371 371 '''
372 372
373 373 def run(self, dataOut, getSNR = True, path=None, file=None, groupList=None):
374 374
375 375
376 376 if path != None:
377 377 sys.path.append(path)
378 378 self.dataOut.library = importlib.import_module(file)
379 379
380 380 #To be inserted as a parameter
381 381 groupArray = numpy.array(groupList)
382 382 # groupArray = numpy.array([[0,1],[2,3]])
383 383 self.dataOut.groupList = groupArray
384 384
385 385 nGroups = groupArray.shape[0]
386 386 nChannels = self.dataIn.nChannels
387 387 nHeights=self.dataIn.heightList.size
388 388
389 389 #Parameters Array
390 390 self.dataOut.data_param = None
391 391
392 392 #Set constants
393 393 constants = self.dataOut.library.setConstants(self.dataIn)
394 394 self.dataOut.constants = constants
395 395 M = self.dataIn.normFactor
396 396 N = self.dataIn.nFFTPoints
397 397 ippSeconds = self.dataIn.ippSeconds
398 398 K = self.dataIn.nIncohInt
399 399 pairsArray = numpy.array(self.dataIn.pairsList)
400 400
401 401 #List of possible combinations
402 402 listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
403 403 indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
404 404
405 405 if getSNR:
406 406 listChannels = groupArray.reshape((groupArray.size))
407 407 listChannels.sort()
408 408 noise = self.dataIn.getNoise()
409 409 self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
410 410
411 411 for i in range(nGroups):
412 412 coord = groupArray[i,:]
413 413
414 414 #Input data array
415 415 data = self.dataIn.data_spc[coord,:,:]/(M*N)
416 416 data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
417 417
418 418 #Cross Spectra data array for Covariance Matrixes
419 419 ind = 0
420 420 for pairs in listComb:
421 421 pairsSel = numpy.array([coord[x],coord[y]])
422 422 indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
423 423 ind += 1
424 424 dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
425 425 dataCross = dataCross**2/K
426 426
427 427 for h in range(nHeights):
428 428 # print self.dataOut.heightList[h]
429 429
430 430 #Input
431 431 d = data[:,h]
432 432
433 433 #Covariance Matrix
434 434 D = numpy.diag(d**2/K)
435 435 ind = 0
436 436 for pairs in listComb:
437 437 #Coordinates in Covariance Matrix
438 438 x = pairs[0]
439 439 y = pairs[1]
440 440 #Channel Index
441 441 S12 = dataCross[ind,:,h]
442 442 D12 = numpy.diag(S12)
443 443 #Completing Covariance Matrix with Cross Spectras
444 444 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
445 445 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
446 446 ind += 1
447 447 Dinv=numpy.linalg.inv(D)
448 448 L=numpy.linalg.cholesky(Dinv)
449 449 LT=L.T
450 450
451 451 dp = numpy.dot(LT,d)
452 452
453 453 #Initial values
454 454 data_spc = self.dataIn.data_spc[coord,:,h]
455 455
456 456 if (h>0)and(error1[3]<5):
457 457 p0 = self.dataOut.data_param[i,:,h-1]
458 458 else:
459 459 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
460 460
461 461 try:
462 462 #Least Squares
463 463 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
464 464 # minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
465 465 #Chi square error
466 466 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
467 467 #Error with Jacobian
468 468 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
469 469 except:
470 470 minp = p0*numpy.nan
471 471 error0 = numpy.nan
472 472 error1 = p0*numpy.nan
473 473
474 474 #Save
475 475 if self.dataOut.data_param is None:
476 476 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
477 477 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
478 478
479 479 self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
480 480 self.dataOut.data_param[i,:,h] = minp
481 481 return
482 482
483 483 def __residFunction(self, p, dp, LT, constants):
484 484
485 485 fm = self.dataOut.library.modelFunction(p, constants)
486 486 fmp=numpy.dot(LT,fm)
487 487
488 488 return dp-fmp
489 489
490 490 def __getSNR(self, z, noise):
491 491
492 492 avg = numpy.average(z, axis=1)
493 493 SNR = (avg.T-noise)/noise
494 494 SNR = SNR.T
495 495 return SNR
496 496
497 497 def __chisq(p,chindex,hindex):
498 498 #similar to Resid but calculates CHI**2
499 499 [LT,d,fm]=setupLTdfm(p,chindex,hindex)
500 500 dp=numpy.dot(LT,d)
501 501 fmp=numpy.dot(LT,fm)
502 502 chisq=numpy.dot((dp-fmp).T,(dp-fmp))
503 503 return chisq
504 504
505 505 class WindProfiler(Operation):
506 506
507 507 __isConfig = False
508 508
509 509 __initime = None
510 510 __lastdatatime = None
511 511 __integrationtime = None
512 512
513 513 __buffer = None
514 514
515 515 __dataReady = False
516 516
517 517 __firstdata = None
518 518
519 519 n = None
520 520
521 def __init__(self):
522 Operation.__init__(self)
523
524 521 def __calculateCosDir(self, elev, azim):
525 522 zen = (90 - elev)*numpy.pi/180
526 523 azim = azim*numpy.pi/180
527 524 cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
528 525 cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
529 526
530 527 signX = numpy.sign(numpy.cos(azim))
531 528 signY = numpy.sign(numpy.sin(azim))
532 529
533 530 cosDirX = numpy.copysign(cosDirX, signX)
534 531 cosDirY = numpy.copysign(cosDirY, signY)
535 532 return cosDirX, cosDirY
536 533
537 534 def __calculateAngles(self, theta_x, theta_y, azimuth):
538 535
539 536 dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
540 537 zenith_arr = numpy.arccos(dir_cosw)
541 538 azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
542 539
543 540 dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
544 541 dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
545 542
546 543 return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
547 544
548 545 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
549 546
550 547 #
551 548 if horOnly:
552 549 A = numpy.c_[dir_cosu,dir_cosv]
553 550 else:
554 551 A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
555 552 A = numpy.asmatrix(A)
556 553 A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
557 554
558 555 return A1
559 556
560 557 def __correctValues(self, heiRang, phi, velRadial, SNR):
561 558 listPhi = phi.tolist()
562 559 maxid = listPhi.index(max(listPhi))
563 560 minid = listPhi.index(min(listPhi))
564 561
565 562 rango = range(len(phi))
566 563 # rango = numpy.delete(rango,maxid)
567 564
568 565 heiRang1 = heiRang*math.cos(phi[maxid])
569 566 heiRangAux = heiRang*math.cos(phi[minid])
570 567 indOut = (heiRang1 < heiRangAux[0]).nonzero()
571 568 heiRang1 = numpy.delete(heiRang1,indOut)
572 569
573 570 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
574 571 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
575 572
576 573 for i in rango:
577 574 x = heiRang*math.cos(phi[i])
578 575 y1 = velRadial[i,:]
579 576 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
580 577
581 578 x1 = heiRang1
582 579 y11 = f1(x1)
583 580
584 581 y2 = SNR[i,:]
585 582 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
586 583 y21 = f2(x1)
587 584
588 585 velRadial1[i,:] = y11
589 586 SNR1[i,:] = y21
590 587
591 588 return heiRang1, velRadial1, SNR1
592 589
593 590 def __calculateVelUVW(self, A, velRadial):
594 591
595 592 #Operacion Matricial
596 593 # velUVW = numpy.zeros((velRadial.shape[1],3))
597 594 # for ind in range(velRadial.shape[1]):
598 595 # velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
599 596 # velUVW = velUVW.transpose()
600 597 velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
601 598 velUVW[:,:] = numpy.dot(A,velRadial)
602 599
603 600
604 601 return velUVW
605 602
606 603 # def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
607 604
608 605 def techniqueDBS(self, kwargs):
609 606 """
610 607 Function that implements Doppler Beam Swinging (DBS) technique.
611 608
612 609 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
613 610 Direction correction (if necessary), Ranges and SNR
614 611
615 612 Output: Winds estimation (Zonal, Meridional and Vertical)
616 613
617 614 Parameters affected: Winds, height range, SNR
618 615 """
619 616 velRadial0 = kwargs['velRadial']
620 617 heiRang = kwargs['heightList']
621 618 SNR0 = kwargs['SNR']
622 619
623 620 if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
624 621 theta_x = numpy.array(kwargs['dirCosx'])
625 622 theta_y = numpy.array(kwargs['dirCosy'])
626 623 else:
627 624 elev = numpy.array(kwargs['elevation'])
628 625 azim = numpy.array(kwargs['azimuth'])
629 626 theta_x, theta_y = self.__calculateCosDir(elev, azim)
630 627 azimuth = kwargs['correctAzimuth']
631 628 if kwargs.has_key('horizontalOnly'):
632 629 horizontalOnly = kwargs['horizontalOnly']
633 630 else: horizontalOnly = False
634 631 if kwargs.has_key('correctFactor'):
635 632 correctFactor = kwargs['correctFactor']
636 633 else: correctFactor = 1
637 634 if kwargs.has_key('channelList'):
638 635 channelList = kwargs['channelList']
639 636 if len(channelList) == 2:
640 637 horizontalOnly = True
641 638 arrayChannel = numpy.array(channelList)
642 639 param = param[arrayChannel,:,:]
643 640 theta_x = theta_x[arrayChannel]
644 641 theta_y = theta_y[arrayChannel]
645 642
646 643 azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(theta_x, theta_y, azimuth)
647 644 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correctFactor*velRadial0, SNR0)
648 645 A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
649 646
650 647 #Calculo de Componentes de la velocidad con DBS
651 648 winds = self.__calculateVelUVW(A,velRadial1)
652 649
653 650 return winds, heiRang1, SNR1
654 651
655 652 def __calculateDistance(self, posx, posy, pairs_ccf, azimuth = None):
656 653
657 654 nPairs = len(pairs_ccf)
658 655 posx = numpy.asarray(posx)
659 656 posy = numpy.asarray(posy)
660 657
661 658 #Rotacion Inversa para alinear con el azimuth
662 659 if azimuth!= None:
663 660 azimuth = azimuth*math.pi/180
664 661 posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
665 662 posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
666 663 else:
667 664 posx1 = posx
668 665 posy1 = posy
669 666
670 667 #Calculo de Distancias
671 668 distx = numpy.zeros(nPairs)
672 669 disty = numpy.zeros(nPairs)
673 670 dist = numpy.zeros(nPairs)
674 671 ang = numpy.zeros(nPairs)
675 672
676 673 for i in range(nPairs):
677 674 distx[i] = posx1[pairs_ccf[i][1]] - posx1[pairs_ccf[i][0]]
678 675 disty[i] = posy1[pairs_ccf[i][1]] - posy1[pairs_ccf[i][0]]
679 676 dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
680 677 ang[i] = numpy.arctan2(disty[i],distx[i])
681 678
682 679 return distx, disty, dist, ang
683 680 #Calculo de Matrices
684 681 # nPairs = len(pairs)
685 682 # ang1 = numpy.zeros((nPairs, 2, 1))
686 683 # dist1 = numpy.zeros((nPairs, 2, 1))
687 684 #
688 685 # for j in range(nPairs):
689 686 # dist1[j,0,0] = dist[pairs[j][0]]
690 687 # dist1[j,1,0] = dist[pairs[j][1]]
691 688 # ang1[j,0,0] = ang[pairs[j][0]]
692 689 # ang1[j,1,0] = ang[pairs[j][1]]
693 690 #
694 691 # return distx,disty, dist1,ang1
695 692
696 693
697 694 def __calculateVelVer(self, phase, lagTRange, _lambda):
698 695
699 696 Ts = lagTRange[1] - lagTRange[0]
700 697 velW = -_lambda*phase/(4*math.pi*Ts)
701 698
702 699 return velW
703 700
704 701 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
705 702 nPairs = tau1.shape[0]
706 703 nHeights = tau1.shape[1]
707 704 vel = numpy.zeros((nPairs,3,nHeights))
708 705 dist1 = numpy.reshape(dist, (dist.size,1))
709 706
710 707 angCos = numpy.cos(ang)
711 708 angSin = numpy.sin(ang)
712 709
713 710 vel0 = dist1*tau1/(2*tau2**2)
714 711 vel[:,0,:] = (vel0*angCos).sum(axis = 1)
715 712 vel[:,1,:] = (vel0*angSin).sum(axis = 1)
716 713
717 714 ind = numpy.where(numpy.isinf(vel))
718 715 vel[ind] = numpy.nan
719 716
720 717 return vel
721 718
722 719 # def __getPairsAutoCorr(self, pairsList, nChannels):
723 720 #
724 721 # pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
725 722 #
726 723 # for l in range(len(pairsList)):
727 724 # firstChannel = pairsList[l][0]
728 725 # secondChannel = pairsList[l][1]
729 726 #
730 727 # #Obteniendo pares de Autocorrelacion
731 728 # if firstChannel == secondChannel:
732 729 # pairsAutoCorr[firstChannel] = int(l)
733 730 #
734 731 # pairsAutoCorr = pairsAutoCorr.astype(int)
735 732 #
736 733 # pairsCrossCorr = range(len(pairsList))
737 734 # pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
738 735 #
739 736 # return pairsAutoCorr, pairsCrossCorr
740 737
741 738 # def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
742 739 def techniqueSA(self, kwargs):
743 740
744 741 """
745 742 Function that implements Spaced Antenna (SA) technique.
746 743
747 744 Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
748 745 Direction correction (if necessary), Ranges and SNR
749 746
750 747 Output: Winds estimation (Zonal, Meridional and Vertical)
751 748
752 749 Parameters affected: Winds
753 750 """
754 751 position_x = kwargs['positionX']
755 752 position_y = kwargs['positionY']
756 753 azimuth = kwargs['azimuth']
757 754
758 755 if kwargs.has_key('correctFactor'):
759 756 correctFactor = kwargs['correctFactor']
760 757 else:
761 758 correctFactor = 1
762 759
763 760 groupList = kwargs['groupList']
764 761 pairs_ccf = groupList[1]
765 762 tau = kwargs['tau']
766 763 _lambda = kwargs['_lambda']
767 764
768 765 #Cross Correlation pairs obtained
769 766 # pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairssList, nChannels)
770 767 # pairsArray = numpy.array(pairsList)[pairsCrossCorr]
771 768 # pairsSelArray = numpy.array(pairsSelected)
772 769 # pairs = []
773 770 #
774 771 # #Wind estimation pairs obtained
775 772 # for i in range(pairsSelArray.shape[0]/2):
776 773 # ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
777 774 # ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
778 775 # pairs.append((ind1,ind2))
779 776
780 777 indtau = tau.shape[0]/2
781 778 tau1 = tau[:indtau,:]
782 779 tau2 = tau[indtau:-1,:]
783 780 # tau1 = tau1[pairs,:]
784 781 # tau2 = tau2[pairs,:]
785 782 phase1 = tau[-1,:]
786 783
787 784 #---------------------------------------------------------------------
788 785 #Metodo Directo
789 786 distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairs_ccf,azimuth)
790 787 winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
791 788 winds = stats.nanmean(winds, axis=0)
792 789 #---------------------------------------------------------------------
793 790 #Metodo General
794 791 # distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
795 792 # #Calculo Coeficientes de Funcion de Correlacion
796 793 # F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
797 794 # #Calculo de Velocidades
798 795 # winds = self.calculateVelUV(F,G,A,B,H)
799 796
800 797 #---------------------------------------------------------------------
801 798 winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
802 799 winds = correctFactor*winds
803 800 return winds
804 801
805 802 def __checkTime(self, currentTime, paramInterval, outputInterval):
806 803
807 804 dataTime = currentTime + paramInterval
808 805 deltaTime = dataTime - self.__initime
809 806
810 807 if deltaTime >= outputInterval or deltaTime < 0:
811 808 self.__dataReady = True
812 809 return
813 810
814 811 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax, binkm=2):
815 812 '''
816 813 Function that implements winds estimation technique with detected meteors.
817 814
818 815 Input: Detected meteors, Minimum meteor quantity to wind estimation
819 816
820 817 Output: Winds estimation (Zonal and Meridional)
821 818
822 819 Parameters affected: Winds
823 820 '''
824 821 # print arrayMeteor.shape
825 822 #Settings
826 823 nInt = (heightMax - heightMin)/binkm
827 824 # print nInt
828 825 nInt = int(nInt)
829 826 # print nInt
830 827 winds = numpy.zeros((2,nInt))*numpy.nan
831 828
832 829 #Filter errors
833 830 error = numpy.where(arrayMeteor[:,-1] == 0)[0]
834 831 finalMeteor = arrayMeteor[error,:]
835 832
836 833 #Meteor Histogram
837 834 finalHeights = finalMeteor[:,2]
838 835 hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
839 836 nMeteorsPerI = hist[0]
840 837 heightPerI = hist[1]
841 838
842 839 #Sort of meteors
843 840 indSort = finalHeights.argsort()
844 841 finalMeteor2 = finalMeteor[indSort,:]
845 842
846 843 # Calculating winds
847 844 ind1 = 0
848 845 ind2 = 0
849 846
850 847 for i in range(nInt):
851 848 nMet = nMeteorsPerI[i]
852 849 ind1 = ind2
853 850 ind2 = ind1 + nMet
854 851
855 852 meteorAux = finalMeteor2[ind1:ind2,:]
856 853
857 854 if meteorAux.shape[0] >= meteorThresh:
858 855 vel = meteorAux[:, 6]
859 856 zen = meteorAux[:, 4]*numpy.pi/180
860 857 azim = meteorAux[:, 3]*numpy.pi/180
861 858
862 859 n = numpy.cos(zen)
863 860 # m = (1 - n**2)/(1 - numpy.tan(azim)**2)
864 861 # l = m*numpy.tan(azim)
865 862 l = numpy.sin(zen)*numpy.sin(azim)
866 863 m = numpy.sin(zen)*numpy.cos(azim)
867 864
868 865 A = numpy.vstack((l, m)).transpose()
869 866 A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
870 867 windsAux = numpy.dot(A1, vel)
871 868
872 869 winds[0,i] = windsAux[0]
873 870 winds[1,i] = windsAux[1]
874 871
875 872 return winds, heightPerI[:-1]
876 873
877 874 def techniqueNSM_SA(self, **kwargs):
878 875 metArray = kwargs['metArray']
879 876 heightList = kwargs['heightList']
880 877 timeList = kwargs['timeList']
881 878
882 879 rx_location = kwargs['rx_location']
883 880 groupList = kwargs['groupList']
884 881 azimuth = kwargs['azimuth']
885 882 dfactor = kwargs['dfactor']
886 883 k = kwargs['k']
887 884
888 885 azimuth1, dist = self.__calculateAzimuth1(rx_location, groupList, azimuth)
889 886 d = dist*dfactor
890 887 #Phase calculation
891 888 metArray1 = self.__getPhaseSlope(metArray, heightList, timeList)
892 889
893 890 metArray1[:,-2] = metArray1[:,-2]*metArray1[:,2]*1000/(k*d[metArray1[:,1].astype(int)]) #angles into velocities
894 891
895 892 velEst = numpy.zeros((heightList.size,2))*numpy.nan
896 893 azimuth1 = azimuth1*numpy.pi/180
897 894
898 895 for i in range(heightList.size):
899 896 h = heightList[i]
900 897 indH = numpy.where((metArray1[:,2] == h)&(numpy.abs(metArray1[:,-2]) < 100))[0]
901 898 metHeight = metArray1[indH,:]
902 899 if metHeight.shape[0] >= 2:
903 900 velAux = numpy.asmatrix(metHeight[:,-2]).T #Radial Velocities
904 901 iazim = metHeight[:,1].astype(int)
905 902 azimAux = numpy.asmatrix(azimuth1[iazim]).T #Azimuths
906 903 A = numpy.hstack((numpy.cos(azimAux),numpy.sin(azimAux)))
907 904 A = numpy.asmatrix(A)
908 905 A1 = numpy.linalg.pinv(A.transpose()*A)*A.transpose()
909 906 velHor = numpy.dot(A1,velAux)
910 907
911 908 velEst[i,:] = numpy.squeeze(velHor)
912 909 return velEst
913 910
914 911 def __getPhaseSlope(self, metArray, heightList, timeList):
915 912 meteorList = []
916 913 #utctime sec1 height SNR velRad ph0 ph1 ph2 coh0 coh1 coh2
917 914 #Putting back together the meteor matrix
918 915 utctime = metArray[:,0]
919 916 uniqueTime = numpy.unique(utctime)
920 917
921 918 phaseDerThresh = 0.5
922 919 ippSeconds = timeList[1] - timeList[0]
923 920 sec = numpy.where(timeList>1)[0][0]
924 921 nPairs = metArray.shape[1] - 6
925 922 nHeights = len(heightList)
926 923
927 924 for t in uniqueTime:
928 925 metArray1 = metArray[utctime==t,:]
929 926 # phaseDerThresh = numpy.pi/4 #reducir Phase thresh
930 927 tmet = metArray1[:,1].astype(int)
931 928 hmet = metArray1[:,2].astype(int)
932 929
933 930 metPhase = numpy.zeros((nPairs, heightList.size, timeList.size - 1))
934 931 metPhase[:,:] = numpy.nan
935 932 metPhase[:,hmet,tmet] = metArray1[:,6:].T
936 933
937 934 #Delete short trails
938 935 metBool = ~numpy.isnan(metPhase[0,:,:])
939 936 heightVect = numpy.sum(metBool, axis = 1)
940 937 metBool[heightVect<sec,:] = False
941 938 metPhase[:,heightVect<sec,:] = numpy.nan
942 939
943 940 #Derivative
944 941 metDer = numpy.abs(metPhase[:,:,1:] - metPhase[:,:,:-1])
945 942 phDerAux = numpy.dstack((numpy.full((nPairs,nHeights,1), False, dtype=bool),metDer > phaseDerThresh))
946 943 metPhase[phDerAux] = numpy.nan
947 944
948 945 #--------------------------METEOR DETECTION -----------------------------------------
949 946 indMet = numpy.where(numpy.any(metBool,axis=1))[0]
950 947
951 948 for p in numpy.arange(nPairs):
952 949 phase = metPhase[p,:,:]
953 950 phDer = metDer[p,:,:]
954 951
955 952 for h in indMet:
956 953 height = heightList[h]
957 954 phase1 = phase[h,:] #82
958 955 phDer1 = phDer[h,:]
959 956
960 957 phase1[~numpy.isnan(phase1)] = numpy.unwrap(phase1[~numpy.isnan(phase1)]) #Unwrap
961 958
962 959 indValid = numpy.where(~numpy.isnan(phase1))[0]
963 960 initMet = indValid[0]
964 961 endMet = 0
965 962
966 963 for i in range(len(indValid)-1):
967 964
968 965 #Time difference
969 966 inow = indValid[i]
970 967 inext = indValid[i+1]
971 968 idiff = inext - inow
972 969 #Phase difference
973 970 phDiff = numpy.abs(phase1[inext] - phase1[inow])
974 971
975 972 if idiff>sec or phDiff>numpy.pi/4 or inext==indValid[-1]: #End of Meteor
976 973 sizeTrail = inow - initMet + 1
977 974 if sizeTrail>3*sec: #Too short meteors
978 975 x = numpy.arange(initMet,inow+1)*ippSeconds
979 976 y = phase1[initMet:inow+1]
980 977 ynnan = ~numpy.isnan(y)
981 978 x = x[ynnan]
982 979 y = y[ynnan]
983 980 slope, intercept, r_value, p_value, std_err = stats.linregress(x,y)
984 981 ylin = x*slope + intercept
985 982 rsq = r_value**2
986 983 if rsq > 0.5:
987 984 vel = slope#*height*1000/(k*d)
988 985 estAux = numpy.array([utctime,p,height, vel, rsq])
989 986 meteorList.append(estAux)
990 987 initMet = inext
991 988 metArray2 = numpy.array(meteorList)
992 989
993 990 return metArray2
994 991
995 992 def __calculateAzimuth1(self, rx_location, pairslist, azimuth0):
996 993
997 994 azimuth1 = numpy.zeros(len(pairslist))
998 995 dist = numpy.zeros(len(pairslist))
999 996
1000 997 for i in range(len(rx_location)):
1001 998 ch0 = pairslist[i][0]
1002 999 ch1 = pairslist[i][1]
1003 1000
1004 1001 diffX = rx_location[ch0][0] - rx_location[ch1][0]
1005 1002 diffY = rx_location[ch0][1] - rx_location[ch1][1]
1006 1003 azimuth1[i] = numpy.arctan2(diffY,diffX)*180/numpy.pi
1007 1004 dist[i] = numpy.sqrt(diffX**2 + diffY**2)
1008 1005
1009 1006 azimuth1 -= azimuth0
1010 1007 return azimuth1, dist
1011 1008
1012 1009 def techniqueNSM_DBS(self, **kwargs):
1013 1010 metArray = kwargs['metArray']
1014 1011 heightList = kwargs['heightList']
1015 1012 timeList = kwargs['timeList']
1016 1013 zenithList = kwargs['zenithList']
1017 1014 nChan = numpy.max(cmet) + 1
1018 1015 nHeights = len(heightList)
1019 1016
1020 1017 utctime = metArray[:,0]
1021 1018 cmet = metArray[:,1]
1022 1019 hmet = metArray1[:,3].astype(int)
1023 1020 h1met = heightList[hmet]*zenithList[cmet]
1024 1021 vmet = metArray1[:,5]
1025 1022
1026 1023 for i in range(nHeights - 1):
1027 1024 hmin = heightList[i]
1028 1025 hmax = heightList[i + 1]
1029 1026
1030 1027 vthisH = vmet[(h1met>=hmin) & (h1met<hmax)]
1031 1028
1032 1029
1033 1030
1034 1031 return data_output
1035 1032
1036 1033 def run(self, dataOut, technique, **kwargs):
1037 1034
1038 1035 param = dataOut.data_param
1039 1036 if dataOut.abscissaList != None:
1040 1037 absc = dataOut.abscissaList[:-1]
1041 1038 noise = dataOut.noise
1042 1039 heightList = dataOut.heightList
1043 1040 SNR = dataOut.data_SNR
1044 1041
1045 1042 if technique == 'DBS':
1046 1043
1047 1044 kwargs['velRadial'] = param[:,1,:] #Radial velocity
1048 1045 kwargs['heightList'] = heightList
1049 1046 kwargs['SNR'] = SNR
1050 1047
1051 1048 dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(kwargs) #DBS Function
1052 1049 dataOut.utctimeInit = dataOut.utctime
1053 1050 dataOut.outputInterval = dataOut.paramInterval
1054 1051
1055 1052 elif technique == 'SA':
1056 1053
1057 1054 #Parameters
1058 1055 # position_x = kwargs['positionX']
1059 1056 # position_y = kwargs['positionY']
1060 1057 # azimuth = kwargs['azimuth']
1061 1058 #
1062 1059 # if kwargs.has_key('crosspairsList'):
1063 1060 # pairs = kwargs['crosspairsList']
1064 1061 # else:
1065 1062 # pairs = None
1066 1063 #
1067 1064 # if kwargs.has_key('correctFactor'):
1068 1065 # correctFactor = kwargs['correctFactor']
1069 1066 # else:
1070 1067 # correctFactor = 1
1071 1068
1072 1069 # tau = dataOut.data_param
1073 1070 # _lambda = dataOut.C/dataOut.frequency
1074 1071 # pairsList = dataOut.groupList
1075 1072 # nChannels = dataOut.nChannels
1076 1073
1077 1074 kwargs['groupList'] = dataOut.groupList
1078 1075 kwargs['tau'] = dataOut.data_param
1079 1076 kwargs['_lambda'] = dataOut.C/dataOut.frequency
1080 1077 # dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
1081 1078 dataOut.data_output = self.techniqueSA(kwargs)
1082 1079 dataOut.utctimeInit = dataOut.utctime
1083 1080 dataOut.outputInterval = dataOut.timeInterval
1084 1081
1085 1082 elif technique == 'Meteors':
1086 1083 dataOut.flagNoData = True
1087 1084 self.__dataReady = False
1088 1085
1089 1086 if kwargs.has_key('nHours'):
1090 1087 nHours = kwargs['nHours']
1091 1088 else:
1092 1089 nHours = 1
1093 1090
1094 1091 if kwargs.has_key('meteorsPerBin'):
1095 1092 meteorThresh = kwargs['meteorsPerBin']
1096 1093 else:
1097 1094 meteorThresh = 6
1098 1095
1099 1096 if kwargs.has_key('hmin'):
1100 1097 hmin = kwargs['hmin']
1101 1098 else: hmin = 70
1102 1099 if kwargs.has_key('hmax'):
1103 1100 hmax = kwargs['hmax']
1104 1101 else: hmax = 110
1105 1102
1106 1103 if kwargs.has_key('BinKm'):
1107 1104 binkm = kwargs['BinKm']
1108 1105 else:
1109 1106 binkm = 2
1110 1107
1111 1108 dataOut.outputInterval = nHours*3600
1112 1109
1113 1110 if self.__isConfig == False:
1114 1111 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1115 1112 #Get Initial LTC time
1116 1113 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1117 1114 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1118 1115
1119 1116 self.__isConfig = True
1120 1117
1121 1118 if self.__buffer is None:
1122 1119 self.__buffer = dataOut.data_param
1123 1120 self.__firstdata = copy.copy(dataOut)
1124 1121
1125 1122 else:
1126 1123 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1127 1124
1128 1125 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1129 1126
1130 1127 if self.__dataReady:
1131 1128 dataOut.utctimeInit = self.__initime
1132 1129
1133 1130 self.__initime += dataOut.outputInterval #to erase time offset
1134 1131
1135 1132 dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax, binkm)
1136 1133 dataOut.flagNoData = False
1137 1134 self.__buffer = None
1138 1135
1139 1136 elif technique == 'Meteors1':
1140 1137 dataOut.flagNoData = True
1141 1138 self.__dataReady = False
1142 1139
1143 1140 if kwargs.has_key('nMins'):
1144 1141 nMins = kwargs['nMins']
1145 1142 else: nMins = 20
1146 1143 if kwargs.has_key('rx_location'):
1147 1144 rx_location = kwargs['rx_location']
1148 1145 else: rx_location = [(0,1),(1,1),(1,0)]
1149 1146 if kwargs.has_key('azimuth'):
1150 1147 azimuth = kwargs['azimuth']
1151 1148 else: azimuth = 51
1152 1149 if kwargs.has_key('dfactor'):
1153 1150 dfactor = kwargs['dfactor']
1154 1151 if kwargs.has_key('mode'):
1155 1152 mode = kwargs['mode']
1156 1153 else: mode = 'SA'
1157 1154
1158 1155 #Borrar luego esto
1159 1156 if dataOut.groupList is None:
1160 1157 dataOut.groupList = [(0,1),(0,2),(1,2)]
1161 1158 groupList = dataOut.groupList
1162 1159 C = 3e8
1163 1160 freq = 50e6
1164 1161 lamb = C/freq
1165 1162 k = 2*numpy.pi/lamb
1166 1163
1167 1164 timeList = dataOut.abscissaList
1168 1165 heightList = dataOut.heightList
1169 1166
1170 1167 if self.__isConfig == False:
1171 1168 dataOut.outputInterval = nMins*60
1172 1169 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
1173 1170 #Get Initial LTC time
1174 1171 initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
1175 1172 minuteAux = initime.minute
1176 1173 minuteNew = int(numpy.floor(minuteAux/nMins)*nMins)
1177 1174 self.__initime = (initime.replace(minute = minuteNew, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
1178 1175
1179 1176 self.__isConfig = True
1180 1177
1181 1178 if self.__buffer is None:
1182 1179 self.__buffer = dataOut.data_param
1183 1180 self.__firstdata = copy.copy(dataOut)
1184 1181
1185 1182 else:
1186 1183 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
1187 1184
1188 1185 self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
1189 1186
1190 1187 if self.__dataReady:
1191 1188 dataOut.utctimeInit = self.__initime
1192 1189 self.__initime += dataOut.outputInterval #to erase time offset
1193 1190
1194 1191 metArray = self.__buffer
1195 1192 if mode == 'SA':
1196 1193 dataOut.data_output = self.techniqueNSM_SA(rx_location=rx_location, groupList=groupList, azimuth=azimuth, dfactor=dfactor, k=k,metArray=metArray, heightList=heightList,timeList=timeList)
1197 1194 elif mode == 'DBS':
1198 1195 dataOut.data_output = self.techniqueNSM_DBS(metArray=metArray,heightList=heightList,timeList=timeList)
1199 1196 dataOut.data_output = dataOut.data_output.T
1200 1197 dataOut.flagNoData = False
1201 1198 self.__buffer = None
1202 1199
1203 1200 return
1204 1201
1205 1202 class EWDriftsEstimation(Operation):
1206 1203
1207 def __init__(self):
1208 Operation.__init__(self)
1209 1204
1210 1205 def __correctValues(self, heiRang, phi, velRadial, SNR):
1211 1206 listPhi = phi.tolist()
1212 1207 maxid = listPhi.index(max(listPhi))
1213 1208 minid = listPhi.index(min(listPhi))
1214 1209
1215 1210 rango = range(len(phi))
1216 1211 # rango = numpy.delete(rango,maxid)
1217 1212
1218 1213 heiRang1 = heiRang*math.cos(phi[maxid])
1219 1214 heiRangAux = heiRang*math.cos(phi[minid])
1220 1215 indOut = (heiRang1 < heiRangAux[0]).nonzero()
1221 1216 heiRang1 = numpy.delete(heiRang1,indOut)
1222 1217
1223 1218 velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
1224 1219 SNR1 = numpy.zeros([len(phi),len(heiRang1)])
1225 1220
1226 1221 for i in rango:
1227 1222 x = heiRang*math.cos(phi[i])
1228 1223 y1 = velRadial[i,:]
1229 1224 f1 = interpolate.interp1d(x,y1,kind = 'cubic')
1230 1225
1231 1226 x1 = heiRang1
1232 1227 y11 = f1(x1)
1233 1228
1234 1229 y2 = SNR[i,:]
1235 1230 f2 = interpolate.interp1d(x,y2,kind = 'cubic')
1236 1231 y21 = f2(x1)
1237 1232
1238 1233 velRadial1[i,:] = y11
1239 1234 SNR1[i,:] = y21
1240 1235
1241 1236 return heiRang1, velRadial1, SNR1
1242 1237
1243 1238 def run(self, dataOut, zenith, zenithCorrection):
1244 1239 heiRang = dataOut.heightList
1245 1240 velRadial = dataOut.data_param[:,3,:]
1246 1241 SNR = dataOut.data_SNR
1247 1242
1248 1243 zenith = numpy.array(zenith)
1249 1244 zenith -= zenithCorrection
1250 1245 zenith *= numpy.pi/180
1251 1246
1252 1247 heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
1253 1248
1254 1249 alp = zenith[0]
1255 1250 bet = zenith[1]
1256 1251
1257 1252 w_w = velRadial1[0,:]
1258 1253 w_e = velRadial1[1,:]
1259 1254
1260 1255 w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
1261 1256 u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
1262 1257
1263 1258 winds = numpy.vstack((u,w))
1264 1259
1265 1260 dataOut.heightList = heiRang1
1266 1261 dataOut.data_output = winds
1267 1262 dataOut.data_SNR = SNR1
1268 1263
1269 1264 dataOut.utctimeInit = dataOut.utctime
1270 1265 dataOut.outputInterval = dataOut.timeInterval
1271 1266 return
1272 1267
1273 1268 #--------------- Non Specular Meteor ----------------
1274 1269
1275 1270 class NonSpecularMeteorDetection(Operation):
1276 1271
1277 1272 def run(self, mode, SNRthresh=8, phaseDerThresh=0.5, cohThresh=0.8, allData = False):
1278 1273 data_acf = self.dataOut.data_pre[0]
1279 1274 data_ccf = self.dataOut.data_pre[1]
1280 1275
1281 1276 lamb = self.dataOut.C/self.dataOut.frequency
1282 1277 tSamp = self.dataOut.ippSeconds*self.dataOut.nCohInt
1283 1278 paramInterval = self.dataOut.paramInterval
1284 1279
1285 1280 nChannels = data_acf.shape[0]
1286 1281 nLags = data_acf.shape[1]
1287 1282 nProfiles = data_acf.shape[2]
1288 1283 nHeights = self.dataOut.nHeights
1289 1284 nCohInt = self.dataOut.nCohInt
1290 1285 sec = numpy.round(nProfiles/self.dataOut.paramInterval)
1291 1286 heightList = self.dataOut.heightList
1292 1287 ippSeconds = self.dataOut.ippSeconds*self.dataOut.nCohInt*self.dataOut.nAvg
1293 1288 utctime = self.dataOut.utctime
1294 1289
1295 1290 self.dataOut.abscissaList = numpy.arange(0,paramInterval+ippSeconds,ippSeconds)
1296 1291
1297 1292 #------------------------ SNR --------------------------------------
1298 1293 power = data_acf[:,0,:,:].real
1299 1294 noise = numpy.zeros(nChannels)
1300 1295 SNR = numpy.zeros(power.shape)
1301 1296 for i in range(nChannels):
1302 1297 noise[i] = hildebrand_sekhon(power[i,:], nCohInt)
1303 1298 SNR[i] = (power[i]-noise[i])/noise[i]
1304 1299 SNRm = numpy.nanmean(SNR, axis = 0)
1305 1300 SNRdB = 10*numpy.log10(SNR)
1306 1301
1307 1302 if mode == 'SA':
1308 1303 nPairs = data_ccf.shape[0]
1309 1304 #---------------------- Coherence and Phase --------------------------
1310 1305 phase = numpy.zeros(data_ccf[:,0,:,:].shape)
1311 1306 # phase1 = numpy.copy(phase)
1312 1307 coh1 = numpy.zeros(data_ccf[:,0,:,:].shape)
1313 1308
1314 1309 for p in range(nPairs):
1315 1310 ch0 = self.dataOut.groupList[p][0]
1316 1311 ch1 = self.dataOut.groupList[p][1]
1317 1312 ccf = data_ccf[p,0,:,:]/numpy.sqrt(data_acf[ch0,0,:,:]*data_acf[ch1,0,:,:])
1318 1313 phase[p,:,:] = ndimage.median_filter(numpy.angle(ccf), size = (5,1)) #median filter
1319 1314 # phase1[p,:,:] = numpy.angle(ccf) #median filter
1320 1315 coh1[p,:,:] = ndimage.median_filter(numpy.abs(ccf), 5) #median filter
1321 1316 # coh1[p,:,:] = numpy.abs(ccf) #median filter
1322 1317 coh = numpy.nanmax(coh1, axis = 0)
1323 1318 # struc = numpy.ones((5,1))
1324 1319 # coh = ndimage.morphology.grey_dilation(coh, size=(10,1))
1325 1320 #---------------------- Radial Velocity ----------------------------
1326 1321 phaseAux = numpy.mean(numpy.angle(data_acf[:,1,:,:]), axis = 0)
1327 1322 velRad = phaseAux*lamb/(4*numpy.pi*tSamp)
1328 1323
1329 1324 if allData:
1330 1325 boolMetFin = ~numpy.isnan(SNRm)
1331 1326 # coh[:-1,:] = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
1332 1327 else:
1333 1328 #------------------------ Meteor mask ---------------------------------
1334 1329 # #SNR mask
1335 1330 # boolMet = (SNRdB>SNRthresh)#|(~numpy.isnan(SNRdB))
1336 1331 #
1337 1332 # #Erase small objects
1338 1333 # boolMet1 = self.__erase_small(boolMet, 2*sec, 5)
1339 1334 #
1340 1335 # auxEEJ = numpy.sum(boolMet1,axis=0)
1341 1336 # indOver = auxEEJ>nProfiles*0.8 #Use this later
1342 1337 # indEEJ = numpy.where(indOver)[0]
1343 1338 # indNEEJ = numpy.where(~indOver)[0]
1344 1339 #
1345 1340 # boolMetFin = boolMet1
1346 1341 #
1347 1342 # if indEEJ.size > 0:
1348 1343 # boolMet1[:,indEEJ] = False #Erase heights with EEJ
1349 1344 #
1350 1345 # boolMet2 = coh > cohThresh
1351 1346 # boolMet2 = self.__erase_small(boolMet2, 2*sec,5)
1352 1347 #
1353 1348 # #Final Meteor mask
1354 1349 # boolMetFin = boolMet1|boolMet2
1355 1350
1356 1351 #Coherence mask
1357 1352 boolMet1 = coh > 0.75
1358 1353 struc = numpy.ones((30,1))
1359 1354 boolMet1 = ndimage.morphology.binary_dilation(boolMet1, structure=struc)
1360 1355
1361 1356 #Derivative mask
1362 1357 derPhase = numpy.nanmean(numpy.abs(phase[:,1:,:] - phase[:,:-1,:]),axis=0)
1363 1358 boolMet2 = derPhase < 0.2
1364 1359 # boolMet2 = ndimage.morphology.binary_opening(boolMet2)
1365 1360 # boolMet2 = ndimage.morphology.binary_closing(boolMet2, structure = numpy.ones((10,1)))
1366 1361 boolMet2 = ndimage.median_filter(boolMet2,size=5)
1367 1362 boolMet2 = numpy.vstack((boolMet2,numpy.full((1,nHeights), True, dtype=bool)))
1368 1363 # #Final mask
1369 1364 # boolMetFin = boolMet2
1370 1365 boolMetFin = boolMet1&boolMet2
1371 1366 # boolMetFin = ndimage.morphology.binary_dilation(boolMetFin)
1372 1367 #Creating data_param
1373 1368 coordMet = numpy.where(boolMetFin)
1374 1369
1375 1370 tmet = coordMet[0]
1376 1371 hmet = coordMet[1]
1377 1372
1378 1373 data_param = numpy.zeros((tmet.size, 6 + nPairs))
1379 1374 data_param[:,0] = utctime
1380 1375 data_param[:,1] = tmet
1381 1376 data_param[:,2] = hmet
1382 1377 data_param[:,3] = SNRm[tmet,hmet]
1383 1378 data_param[:,4] = velRad[tmet,hmet]
1384 1379 data_param[:,5] = coh[tmet,hmet]
1385 1380 data_param[:,6:] = phase[:,tmet,hmet].T
1386 1381
1387 1382 elif mode == 'DBS':
1388 1383 self.dataOut.groupList = numpy.arange(nChannels)
1389 1384
1390 1385 #Radial Velocities
1391 1386 # phase = numpy.angle(data_acf[:,1,:,:])
1392 1387 phase = ndimage.median_filter(numpy.angle(data_acf[:,1,:,:]), size = (1,5,1))
1393 1388 velRad = phase*lamb/(4*numpy.pi*tSamp)
1394 1389
1395 1390 #Spectral width
1396 1391 acf1 = ndimage.median_filter(numpy.abs(data_acf[:,1,:,:]), size = (1,5,1))
1397 1392 acf2 = ndimage.median_filter(numpy.abs(data_acf[:,2,:,:]), size = (1,5,1))
1398 1393
1399 1394 spcWidth = (lamb/(2*numpy.sqrt(6)*numpy.pi*tSamp))*numpy.sqrt(numpy.log(acf1/acf2))
1400 1395 # velRad = ndimage.median_filter(velRad, size = (1,5,1))
1401 1396 if allData:
1402 1397 boolMetFin = ~numpy.isnan(SNRdB)
1403 1398 else:
1404 1399 #SNR
1405 1400 boolMet1 = (SNRdB>SNRthresh) #SNR mask
1406 1401 boolMet1 = ndimage.median_filter(boolMet1, size=(1,5,5))
1407 1402
1408 1403 #Radial velocity
1409 1404 boolMet2 = numpy.abs(velRad) < 30
1410 1405 boolMet2 = ndimage.median_filter(boolMet2, (1,5,5))
1411 1406
1412 1407 #Spectral Width
1413 1408 boolMet3 = spcWidth < 30
1414 1409 boolMet3 = ndimage.median_filter(boolMet3, (1,5,5))
1415 1410 # boolMetFin = self.__erase_small(boolMet1, 10,5)
1416 1411 boolMetFin = boolMet1&boolMet2&boolMet3
1417 1412
1418 1413 #Creating data_param
1419 1414 coordMet = numpy.where(boolMetFin)
1420 1415
1421 1416 cmet = coordMet[0]
1422 1417 tmet = coordMet[1]
1423 1418 hmet = coordMet[2]
1424 1419
1425 1420 data_param = numpy.zeros((tmet.size, 7))
1426 1421 data_param[:,0] = utctime
1427 1422 data_param[:,1] = cmet
1428 1423 data_param[:,2] = tmet
1429 1424 data_param[:,3] = hmet
1430 1425 data_param[:,4] = SNR[cmet,tmet,hmet].T
1431 1426 data_param[:,5] = velRad[cmet,tmet,hmet].T
1432 1427 data_param[:,6] = spcWidth[cmet,tmet,hmet].T
1433 1428
1434 1429 # self.dataOut.data_param = data_int
1435 1430 if len(data_param) == 0:
1436 1431 self.dataOut.flagNoData = True
1437 1432 else:
1438 1433 self.dataOut.data_param = data_param
1439 1434
1440 1435 def __erase_small(self, binArray, threshX, threshY):
1441 1436 labarray, numfeat = ndimage.measurements.label(binArray)
1442 1437 binArray1 = numpy.copy(binArray)
1443 1438
1444 1439 for i in range(1,numfeat + 1):
1445 1440 auxBin = (labarray==i)
1446 1441 auxSize = auxBin.sum()
1447 1442
1448 1443 x,y = numpy.where(auxBin)
1449 1444 widthX = x.max() - x.min()
1450 1445 widthY = y.max() - y.min()
1451 1446
1452 1447 #width X: 3 seg -> 12.5*3
1453 1448 #width Y:
1454 1449
1455 1450 if (auxSize < 50) or (widthX < threshX) or (widthY < threshY):
1456 1451 binArray1[auxBin] = False
1457 1452
1458 1453 return binArray1
1459 1454
1460 1455 #--------------- Specular Meteor ----------------
1461 1456
1462 1457 class SMDetection(Operation):
1463 1458 '''
1464 1459 Function DetectMeteors()
1465 1460 Project developed with paper:
1466 1461 HOLDSWORTH ET AL. 2004
1467 1462
1468 1463 Input:
1469 1464 self.dataOut.data_pre
1470 1465
1471 1466 centerReceiverIndex: From the channels, which is the center receiver
1472 1467
1473 1468 hei_ref: Height reference for the Beacon signal extraction
1474 1469 tauindex:
1475 1470 predefinedPhaseShifts: Predefined phase offset for the voltge signals
1476 1471
1477 1472 cohDetection: Whether to user Coherent detection or not
1478 1473 cohDet_timeStep: Coherent Detection calculation time step
1479 1474 cohDet_thresh: Coherent Detection phase threshold to correct phases
1480 1475
1481 1476 noise_timeStep: Noise calculation time step
1482 1477 noise_multiple: Noise multiple to define signal threshold
1483 1478
1484 1479 multDet_timeLimit: Multiple Detection Removal time limit in seconds
1485 1480 multDet_rangeLimit: Multiple Detection Removal range limit in km
1486 1481
1487 1482 phaseThresh: Maximum phase difference between receiver to be consider a meteor
1488 1483 SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor
1489 1484
1490 1485 hmin: Minimum Height of the meteor to use it in the further wind estimations
1491 1486 hmax: Maximum Height of the meteor to use it in the further wind estimations
1492 1487 azimuth: Azimuth angle correction
1493 1488
1494 1489 Affected:
1495 1490 self.dataOut.data_param
1496 1491
1497 1492 Rejection Criteria (Errors):
1498 1493 0: No error; analysis OK
1499 1494 1: SNR < SNR threshold
1500 1495 2: angle of arrival (AOA) ambiguously determined
1501 1496 3: AOA estimate not feasible
1502 1497 4: Large difference in AOAs obtained from different antenna baselines
1503 1498 5: echo at start or end of time series
1504 1499 6: echo less than 5 examples long; too short for analysis
1505 1500 7: echo rise exceeds 0.3s
1506 1501 8: echo decay time less than twice rise time
1507 1502 9: large power level before echo
1508 1503 10: large power level after echo
1509 1504 11: poor fit to amplitude for estimation of decay time
1510 1505 12: poor fit to CCF phase variation for estimation of radial drift velocity
1511 1506 13: height unresolvable echo: not valid height within 70 to 110 km
1512 1507 14: height ambiguous echo: more then one possible height within 70 to 110 km
1513 1508 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s
1514 1509 16: oscilatory echo, indicating event most likely not an underdense echo
1515 1510
1516 1511 17: phase difference in meteor Reestimation
1517 1512
1518 1513 Data Storage:
1519 1514 Meteors for Wind Estimation (8):
1520 1515 Utc Time | Range Height
1521 1516 Azimuth Zenith errorCosDir
1522 1517 VelRad errorVelRad
1523 1518 Phase0 Phase1 Phase2 Phase3
1524 1519 TypeError
1525 1520
1526 1521 '''
1527 1522
1528 1523 def run(self, dataOut, hei_ref = None, tauindex = 0,
1529 1524 phaseOffsets = None,
1530 1525 cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
1531 1526 noise_timeStep = 4, noise_multiple = 4,
1532 1527 multDet_timeLimit = 1, multDet_rangeLimit = 3,
1533 1528 phaseThresh = 20, SNRThresh = 5,
1534 1529 hmin = 50, hmax=150, azimuth = 0,
1535 1530 channelPositions = None) :
1536 1531
1537 1532
1538 1533 #Getting Pairslist
1539 1534 if channelPositions is None:
1540 1535 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
1541 1536 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
1542 1537 meteorOps = SMOperations()
1543 1538 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
1544 1539 heiRang = dataOut.getHeiRange()
1545 1540 #Get Beacon signal - No Beacon signal anymore
1546 1541 # newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
1547 1542 #
1548 1543 # if hei_ref != None:
1549 1544 # newheis = numpy.where(self.dataOut.heightList>hei_ref)
1550 1545 #
1551 1546
1552 1547
1553 1548 #****************REMOVING HARDWARE PHASE DIFFERENCES***************
1554 1549 # see if the user put in pre defined phase shifts
1555 1550 voltsPShift = dataOut.data_pre.copy()
1556 1551
1557 1552 # if predefinedPhaseShifts != None:
1558 1553 # hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
1559 1554 #
1560 1555 # # elif beaconPhaseShifts:
1561 1556 # # #get hardware phase shifts using beacon signal
1562 1557 # # hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
1563 1558 # # hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
1564 1559 #
1565 1560 # else:
1566 1561 # hardwarePhaseShifts = numpy.zeros(5)
1567 1562 #
1568 1563 # voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
1569 1564 # for i in range(self.dataOut.data_pre.shape[0]):
1570 1565 # voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
1571 1566
1572 1567 #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
1573 1568
1574 1569 #Remove DC
1575 1570 voltsDC = numpy.mean(voltsPShift,1)
1576 1571 voltsDC = numpy.mean(voltsDC,1)
1577 1572 for i in range(voltsDC.shape[0]):
1578 1573 voltsPShift[i] = voltsPShift[i] - voltsDC[i]
1579 1574
1580 1575 #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
1581 1576 # voltsPShift = voltsPShift[:,:,:newheis[0][0]]
1582 1577
1583 1578 #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
1584 1579 #Coherent Detection
1585 1580 if cohDetection:
1586 1581 #use coherent detection to get the net power
1587 1582 cohDet_thresh = cohDet_thresh*numpy.pi/180
1588 1583 voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, dataOut.timeInterval, pairslist0, cohDet_thresh)
1589 1584
1590 1585 #Non-coherent detection!
1591 1586 powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
1592 1587 #********** END OF COH/NON-COH POWER CALCULATION**********************
1593 1588
1594 1589 #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
1595 1590 #Get noise
1596 1591 noise, noise1 = self.__getNoise(powerNet, noise_timeStep, dataOut.timeInterval)
1597 1592 # noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
1598 1593 #Get signal threshold
1599 1594 signalThresh = noise_multiple*noise
1600 1595 #Meteor echoes detection
1601 1596 listMeteors = self.__findMeteors(powerNet, signalThresh)
1602 1597 #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
1603 1598
1604 1599 #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
1605 1600 #Parameters
1606 1601 heiRange = dataOut.getHeiRange()
1607 1602 rangeInterval = heiRange[1] - heiRange[0]
1608 1603 rangeLimit = multDet_rangeLimit/rangeInterval
1609 1604 timeLimit = multDet_timeLimit/dataOut.timeInterval
1610 1605 #Multiple detection removals
1611 1606 listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
1612 1607 #************ END OF REMOVE MULTIPLE DETECTIONS **********************
1613 1608
1614 1609 #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ********************
1615 1610 #Parameters
1616 1611 phaseThresh = phaseThresh*numpy.pi/180
1617 1612 thresh = [phaseThresh, noise_multiple, SNRThresh]
1618 1613 #Meteor reestimation (Errors N 1, 6, 12, 17)
1619 1614 listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist0, thresh, noise, dataOut.timeInterval, dataOut.frequency)
1620 1615 # listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
1621 1616 #Estimation of decay times (Errors N 7, 8, 11)
1622 1617 listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, dataOut.timeInterval, dataOut.frequency)
1623 1618 #******************* END OF METEOR REESTIMATION *******************
1624 1619
1625 1620 #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) **************************
1626 1621 #Calculating Radial Velocity (Error N 15)
1627 1622 radialStdThresh = 10
1628 1623 listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist0, dataOut.timeInterval)
1629 1624
1630 1625 if len(listMeteors4) > 0:
1631 1626 #Setting New Array
1632 1627 date = dataOut.utctime
1633 1628 arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
1634 1629
1635 1630 #Correcting phase offset
1636 1631 if phaseOffsets != None:
1637 1632 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
1638 1633 arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
1639 1634
1640 1635 #Second Pairslist
1641 1636 pairsList = []
1642 1637 pairx = (0,1)
1643 1638 pairy = (2,3)
1644 1639 pairsList.append(pairx)
1645 1640 pairsList.append(pairy)
1646 1641
1647 1642 jph = numpy.array([0,0,0,0])
1648 1643 h = (hmin,hmax)
1649 1644 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
1650 1645
1651 1646 # #Calculate AOA (Error N 3, 4)
1652 1647 # #JONES ET AL. 1998
1653 1648 # error = arrayParameters[:,-1]
1654 1649 # AOAthresh = numpy.pi/8
1655 1650 # phases = -arrayParameters[:,9:13]
1656 1651 # arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
1657 1652 #
1658 1653 # #Calculate Heights (Error N 13 and 14)
1659 1654 # error = arrayParameters[:,-1]
1660 1655 # Ranges = arrayParameters[:,2]
1661 1656 # zenith = arrayParameters[:,5]
1662 1657 # arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
1663 1658 # error = arrayParameters[:,-1]
1664 1659 #********************* END OF PARAMETERS CALCULATION **************************
1665 1660
1666 1661 #***************************+ PASS DATA TO NEXT STEP **********************
1667 1662 # arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
1668 1663 dataOut.data_param = arrayParameters
1669 1664
1670 1665 if arrayParameters is None:
1671 1666 dataOut.flagNoData = True
1672 1667 else:
1673 1668 dataOut.flagNoData = True
1674 1669
1675 1670 return
1676 1671
1677 1672 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
1678 1673
1679 1674 minIndex = min(newheis[0])
1680 1675 maxIndex = max(newheis[0])
1681 1676
1682 1677 voltage = voltage0[:,:,minIndex:maxIndex+1]
1683 1678 nLength = voltage.shape[1]/n
1684 1679 nMin = 0
1685 1680 nMax = 0
1686 1681 phaseOffset = numpy.zeros((len(pairslist),n))
1687 1682
1688 1683 for i in range(n):
1689 1684 nMax += nLength
1690 1685 phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
1691 1686 phaseCCF = numpy.mean(phaseCCF, axis = 2)
1692 1687 phaseOffset[:,i] = phaseCCF.transpose()
1693 1688 nMin = nMax
1694 1689 # phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
1695 1690
1696 1691 #Remove Outliers
1697 1692 factor = 2
1698 1693 wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
1699 1694 dw = numpy.std(wt,axis = 1)
1700 1695 dw = dw.reshape((dw.size,1))
1701 1696 ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
1702 1697 phaseOffset[ind] = numpy.nan
1703 1698 phaseOffset = stats.nanmean(phaseOffset, axis=1)
1704 1699
1705 1700 return phaseOffset
1706 1701
1707 1702 def __shiftPhase(self, data, phaseShift):
1708 1703 #this will shift the phase of a complex number
1709 1704 dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
1710 1705 return dataShifted
1711 1706
1712 1707 def __estimatePhaseDifference(self, array, pairslist):
1713 1708 nChannel = array.shape[0]
1714 1709 nHeights = array.shape[2]
1715 1710 numPairs = len(pairslist)
1716 1711 # phaseCCF = numpy.zeros((nChannel, 5, nHeights))
1717 1712 phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
1718 1713
1719 1714 #Correct phases
1720 1715 derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
1721 1716 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
1722 1717
1723 1718 if indDer[0].shape[0] > 0:
1724 1719 for i in range(indDer[0].shape[0]):
1725 1720 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
1726 1721 phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
1727 1722
1728 1723 # for j in range(numSides):
1729 1724 # phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
1730 1725 # phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
1731 1726 #
1732 1727 #Linear
1733 1728 phaseInt = numpy.zeros((numPairs,1))
1734 1729 angAllCCF = phaseCCF[:,[0,1,3,4],0]
1735 1730 for j in range(numPairs):
1736 1731 fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
1737 1732 phaseInt[j] = fit[1]
1738 1733 #Phase Differences
1739 1734 phaseDiff = phaseInt - phaseCCF[:,2,:]
1740 1735 phaseArrival = phaseInt.reshape(phaseInt.size)
1741 1736
1742 1737 #Dealias
1743 1738 phaseArrival = numpy.angle(numpy.exp(1j*phaseArrival))
1744 1739 # indAlias = numpy.where(phaseArrival > numpy.pi)
1745 1740 # phaseArrival[indAlias] -= 2*numpy.pi
1746 1741 # indAlias = numpy.where(phaseArrival < -numpy.pi)
1747 1742 # phaseArrival[indAlias] += 2*numpy.pi
1748 1743
1749 1744 return phaseDiff, phaseArrival
1750 1745
1751 1746 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
1752 1747 #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
1753 1748 #find the phase shifts of each channel over 1 second intervals
1754 1749 #only look at ranges below the beacon signal
1755 1750 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
1756 1751 numBlocks = int(volts.shape[1]/numProfPerBlock)
1757 1752 numHeights = volts.shape[2]
1758 1753 nChannel = volts.shape[0]
1759 1754 voltsCohDet = volts.copy()
1760 1755
1761 1756 pairsarray = numpy.array(pairslist)
1762 1757 indSides = pairsarray[:,1]
1763 1758 # indSides = numpy.array(range(nChannel))
1764 1759 # indSides = numpy.delete(indSides, indCenter)
1765 1760 #
1766 1761 # listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
1767 1762 listBlocks = numpy.array_split(volts, numBlocks, 1)
1768 1763
1769 1764 startInd = 0
1770 1765 endInd = 0
1771 1766
1772 1767 for i in range(numBlocks):
1773 1768 startInd = endInd
1774 1769 endInd = endInd + listBlocks[i].shape[1]
1775 1770
1776 1771 arrayBlock = listBlocks[i]
1777 1772 # arrayBlockCenter = listCenter[i]
1778 1773
1779 1774 #Estimate the Phase Difference
1780 1775 phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
1781 1776 #Phase Difference RMS
1782 1777 arrayPhaseRMS = numpy.abs(phaseDiff)
1783 1778 phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
1784 1779 indPhase = numpy.where(phaseRMSaux==4)
1785 1780 #Shifting
1786 1781 if indPhase[0].shape[0] > 0:
1787 1782 for j in range(indSides.size):
1788 1783 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
1789 1784 voltsCohDet[:,startInd:endInd,:] = arrayBlock
1790 1785
1791 1786 return voltsCohDet
1792 1787
1793 1788 def __calculateCCF(self, volts, pairslist ,laglist):
1794 1789
1795 1790 nHeights = volts.shape[2]
1796 1791 nPoints = volts.shape[1]
1797 1792 voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
1798 1793
1799 1794 for i in range(len(pairslist)):
1800 1795 volts1 = volts[pairslist[i][0]]
1801 1796 volts2 = volts[pairslist[i][1]]
1802 1797
1803 1798 for t in range(len(laglist)):
1804 1799 idxT = laglist[t]
1805 1800 if idxT >= 0:
1806 1801 vStacked = numpy.vstack((volts2[idxT:,:],
1807 1802 numpy.zeros((idxT, nHeights),dtype='complex')))
1808 1803 else:
1809 1804 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
1810 1805 volts2[:(nPoints + idxT),:]))
1811 1806 voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
1812 1807
1813 1808 vStacked = None
1814 1809 return voltsCCF
1815 1810
1816 1811 def __getNoise(self, power, timeSegment, timeInterval):
1817 1812 numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
1818 1813 numBlocks = int(power.shape[0]/numProfPerBlock)
1819 1814 numHeights = power.shape[1]
1820 1815
1821 1816 listPower = numpy.array_split(power, numBlocks, 0)
1822 1817 noise = numpy.zeros((power.shape[0], power.shape[1]))
1823 1818 noise1 = numpy.zeros((power.shape[0], power.shape[1]))
1824 1819
1825 1820 startInd = 0
1826 1821 endInd = 0
1827 1822
1828 1823 for i in range(numBlocks): #split por canal
1829 1824 startInd = endInd
1830 1825 endInd = endInd + listPower[i].shape[0]
1831 1826
1832 1827 arrayBlock = listPower[i]
1833 1828 noiseAux = numpy.mean(arrayBlock, 0)
1834 1829 # noiseAux = numpy.median(noiseAux)
1835 1830 # noiseAux = numpy.mean(arrayBlock)
1836 1831 noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
1837 1832
1838 1833 noiseAux1 = numpy.mean(arrayBlock)
1839 1834 noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
1840 1835
1841 1836 return noise, noise1
1842 1837
1843 1838 def __findMeteors(self, power, thresh):
1844 1839 nProf = power.shape[0]
1845 1840 nHeights = power.shape[1]
1846 1841 listMeteors = []
1847 1842
1848 1843 for i in range(nHeights):
1849 1844 powerAux = power[:,i]
1850 1845 threshAux = thresh[:,i]
1851 1846
1852 1847 indUPthresh = numpy.where(powerAux > threshAux)[0]
1853 1848 indDNthresh = numpy.where(powerAux <= threshAux)[0]
1854 1849
1855 1850 j = 0
1856 1851
1857 1852 while (j < indUPthresh.size - 2):
1858 1853 if (indUPthresh[j + 2] == indUPthresh[j] + 2):
1859 1854 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
1860 1855 indDNthresh = indDNthresh[indDNAux]
1861 1856
1862 1857 if (indDNthresh.size > 0):
1863 1858 indEnd = indDNthresh[0] - 1
1864 1859 indInit = indUPthresh[j] if isinstance(indUPthresh[j], (int, float)) else indUPthresh[j][0] ##CHECK!!!!
1865 1860
1866 1861 meteor = powerAux[indInit:indEnd + 1]
1867 1862 indPeak = meteor.argmax() + indInit
1868 1863 FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
1869 1864
1870 1865 listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
1871 1866 j = numpy.where(indUPthresh == indEnd)[0] + 1
1872 1867 else: j+=1
1873 1868 else: j+=1
1874 1869
1875 1870 return listMeteors
1876 1871
1877 1872 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
1878 1873
1879 1874 arrayMeteors = numpy.asarray(listMeteors)
1880 1875 listMeteors1 = []
1881 1876
1882 1877 while arrayMeteors.shape[0] > 0:
1883 1878 FLAs = arrayMeteors[:,4]
1884 1879 maxFLA = FLAs.argmax()
1885 1880 listMeteors1.append(arrayMeteors[maxFLA,:])
1886 1881
1887 1882 MeteorInitTime = arrayMeteors[maxFLA,1]
1888 1883 MeteorEndTime = arrayMeteors[maxFLA,3]
1889 1884 MeteorHeight = arrayMeteors[maxFLA,0]
1890 1885
1891 1886 #Check neighborhood
1892 1887 maxHeightIndex = MeteorHeight + rangeLimit
1893 1888 minHeightIndex = MeteorHeight - rangeLimit
1894 1889 minTimeIndex = MeteorInitTime - timeLimit
1895 1890 maxTimeIndex = MeteorEndTime + timeLimit
1896 1891
1897 1892 #Check Heights
1898 1893 indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
1899 1894 indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
1900 1895 indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
1901 1896
1902 1897 arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
1903 1898
1904 1899 return listMeteors1
1905 1900
1906 1901 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
1907 1902 numHeights = volts.shape[2]
1908 1903 nChannel = volts.shape[0]
1909 1904
1910 1905 thresholdPhase = thresh[0]
1911 1906 thresholdNoise = thresh[1]
1912 1907 thresholdDB = float(thresh[2])
1913 1908
1914 1909 thresholdDB1 = 10**(thresholdDB/10)
1915 1910 pairsarray = numpy.array(pairslist)
1916 1911 indSides = pairsarray[:,1]
1917 1912
1918 1913 pairslist1 = list(pairslist)
1919 1914 pairslist1.append((0,4))
1920 1915 pairslist1.append((1,3))
1921 1916
1922 1917 listMeteors1 = []
1923 1918 listPowerSeries = []
1924 1919 listVoltageSeries = []
1925 1920 #volts has the war data
1926 1921
1927 1922 if frequency == 30.175e6:
1928 1923 timeLag = 45*10**-3
1929 1924 else:
1930 1925 timeLag = 15*10**-3
1931 1926 lag = int(numpy.ceil(timeLag/timeInterval))
1932 1927
1933 1928 for i in range(len(listMeteors)):
1934 1929
1935 1930 ###################### 3.6 - 3.7 PARAMETERS REESTIMATION #########################
1936 1931 meteorAux = numpy.zeros(16)
1937 1932
1938 1933 #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
1939 1934 mHeight = int(listMeteors[i][0])
1940 1935 mStart = int(listMeteors[i][1])
1941 1936 mPeak = int(listMeteors[i][2])
1942 1937 mEnd = int(listMeteors[i][3])
1943 1938
1944 1939 #get the volt data between the start and end times of the meteor
1945 1940 meteorVolts = volts[:,mStart:mEnd+1,mHeight]
1946 1941 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
1947 1942
1948 1943 #3.6. Phase Difference estimation
1949 1944 phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
1950 1945
1951 1946 #3.7. Phase difference removal & meteor start, peak and end times reestimated
1952 1947 #meteorVolts0.- all Channels, all Profiles
1953 1948 meteorVolts0 = volts[:,:,mHeight]
1954 1949 meteorThresh = noise[:,mHeight]*thresholdNoise
1955 1950 meteorNoise = noise[:,mHeight]
1956 1951 meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
1957 1952 powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power
1958 1953
1959 1954 #Times reestimation
1960 1955 mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
1961 1956 if mStart1.size > 0:
1962 1957 mStart1 = mStart1[-1] + 1
1963 1958
1964 1959 else:
1965 1960 mStart1 = mPeak
1966 1961
1967 1962 mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
1968 1963 mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
1969 1964 if mEndDecayTime1.size == 0:
1970 1965 mEndDecayTime1 = powerNet0.size
1971 1966 else:
1972 1967 mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
1973 1968 # mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
1974 1969
1975 1970 #meteorVolts1.- all Channels, from start to end
1976 1971 meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
1977 1972 meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
1978 1973 if meteorVolts2.shape[1] == 0:
1979 1974 meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
1980 1975 meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
1981 1976 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
1982 1977 ##################### END PARAMETERS REESTIMATION #########################
1983 1978
1984 1979 ##################### 3.8 PHASE DIFFERENCE REESTIMATION ########################
1985 1980 # if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis
1986 1981 if meteorVolts2.shape[1] > 0:
1987 1982 #Phase Difference re-estimation
1988 1983 phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation
1989 1984 # phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
1990 1985 meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
1991 1986 phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
1992 1987 meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting
1993 1988
1994 1989 #Phase Difference RMS
1995 1990 phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
1996 1991 powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
1997 1992 #Data from Meteor
1998 1993 mPeak1 = powerNet1.argmax() + mStart1
1999 1994 mPeakPower1 = powerNet1.max()
2000 1995 noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
2001 1996 mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
2002 1997 Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
2003 1998 Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
2004 1999 PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1]
2005 2000 #Vectorize
2006 2001 meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
2007 2002 meteorAux[7:11] = phaseDiffint[0:4]
2008 2003
2009 2004 #Rejection Criterions
2010 2005 if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation
2011 2006 meteorAux[-1] = 17
2012 2007 elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB
2013 2008 meteorAux[-1] = 1
2014 2009
2015 2010
2016 2011 else:
2017 2012 meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
2018 2013 meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
2019 2014 PowerSeries = 0
2020 2015
2021 2016 listMeteors1.append(meteorAux)
2022 2017 listPowerSeries.append(PowerSeries)
2023 2018 listVoltageSeries.append(meteorVolts1)
2024 2019
2025 2020 return listMeteors1, listPowerSeries, listVoltageSeries
2026 2021
2027 2022 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
2028 2023
2029 2024 threshError = 10
2030 2025 #Depending if it is 30 or 50 MHz
2031 2026 if frequency == 30.175e6:
2032 2027 timeLag = 45*10**-3
2033 2028 else:
2034 2029 timeLag = 15*10**-3
2035 2030 lag = int(numpy.ceil(timeLag/timeInterval))
2036 2031
2037 2032 listMeteors1 = []
2038 2033
2039 2034 for i in range(len(listMeteors)):
2040 2035 meteorPower = listPower[i]
2041 2036 meteorAux = listMeteors[i]
2042 2037
2043 2038 if meteorAux[-1] == 0:
2044 2039
2045 2040 try:
2046 2041 indmax = meteorPower.argmax()
2047 2042 indlag = indmax + lag
2048 2043
2049 2044 y = meteorPower[indlag:]
2050 2045 x = numpy.arange(0, y.size)*timeLag
2051 2046
2052 2047 #first guess
2053 2048 a = y[0]
2054 2049 tau = timeLag
2055 2050 #exponential fit
2056 2051 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
2057 2052 y1 = self.__exponential_function(x, *popt)
2058 2053 #error estimation
2059 2054 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
2060 2055
2061 2056 decayTime = popt[1]
2062 2057 riseTime = indmax*timeInterval
2063 2058 meteorAux[11:13] = [decayTime, error]
2064 2059
2065 2060 #Table items 7, 8 and 11
2066 2061 if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s
2067 2062 meteorAux[-1] = 7
2068 2063 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
2069 2064 meteorAux[-1] = 8
2070 2065 if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time
2071 2066 meteorAux[-1] = 11
2072 2067
2073 2068
2074 2069 except:
2075 2070 meteorAux[-1] = 11
2076 2071
2077 2072
2078 2073 listMeteors1.append(meteorAux)
2079 2074
2080 2075 return listMeteors1
2081 2076
2082 2077 #Exponential Function
2083 2078
2084 2079 def __exponential_function(self, x, a, tau):
2085 2080 y = a*numpy.exp(-x/tau)
2086 2081 return y
2087 2082
2088 2083 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval):
2089 2084
2090 2085 pairslist1 = list(pairslist)
2091 2086 pairslist1.append((0,4))
2092 2087 pairslist1.append((1,3))
2093 2088 numPairs = len(pairslist1)
2094 2089 #Time Lag
2095 2090 timeLag = 45*10**-3
2096 2091 c = 3e8
2097 2092 lag = numpy.ceil(timeLag/timeInterval)
2098 2093 freq = 30.175e6
2099 2094
2100 2095 listMeteors1 = []
2101 2096
2102 2097 for i in range(len(listMeteors)):
2103 2098 meteorAux = listMeteors[i]
2104 2099 if meteorAux[-1] == 0:
2105 2100 mStart = listMeteors[i][1]
2106 2101 mPeak = listMeteors[i][2]
2107 2102 mLag = mPeak - mStart + lag
2108 2103
2109 2104 #get the volt data between the start and end times of the meteor
2110 2105 meteorVolts = listVolts[i]
2111 2106 meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
2112 2107
2113 2108 #Get CCF
2114 2109 allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
2115 2110
2116 2111 #Method 2
2117 2112 slopes = numpy.zeros(numPairs)
2118 2113 time = numpy.array([-2,-1,1,2])*timeInterval
2119 2114 angAllCCF = numpy.angle(allCCFs[:,[0,4,2,3],0])
2120 2115
2121 2116 #Correct phases
2122 2117 derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
2123 2118 indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
2124 2119
2125 2120 if indDer[0].shape[0] > 0:
2126 2121 for i in range(indDer[0].shape[0]):
2127 2122 signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
2128 2123 angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
2129 2124
2130 2125 # fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
2131 2126 for j in range(numPairs):
2132 2127 fit = stats.linregress(time, angAllCCF[j,:])
2133 2128 slopes[j] = fit[0]
2134 2129
2135 2130 #Remove Outlier
2136 2131 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
2137 2132 # slopes = numpy.delete(slopes,indOut)
2138 2133 # indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
2139 2134 # slopes = numpy.delete(slopes,indOut)
2140 2135
2141 2136 radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
2142 2137 radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
2143 2138 meteorAux[-2] = radialError
2144 2139 meteorAux[-3] = radialVelocity
2145 2140
2146 2141 #Setting Error
2147 2142 #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
2148 2143 if numpy.abs(radialVelocity) > 200:
2149 2144 meteorAux[-1] = 15
2150 2145 #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
2151 2146 elif radialError > radialStdThresh:
2152 2147 meteorAux[-1] = 12
2153 2148
2154 2149 listMeteors1.append(meteorAux)
2155 2150 return listMeteors1
2156 2151
2157 2152 def __setNewArrays(self, listMeteors, date, heiRang):
2158 2153
2159 2154 #New arrays
2160 2155 arrayMeteors = numpy.array(listMeteors)
2161 2156 arrayParameters = numpy.zeros((len(listMeteors), 13))
2162 2157
2163 2158 #Date inclusion
2164 2159 # date = re.findall(r'\((.*?)\)', date)
2165 2160 # date = date[0].split(',')
2166 2161 # date = map(int, date)
2167 2162 #
2168 2163 # if len(date)<6:
2169 2164 # date.append(0)
2170 2165 #
2171 2166 # date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
2172 2167 # arrayDate = numpy.tile(date, (len(listMeteors), 1))
2173 2168 arrayDate = numpy.tile(date, (len(listMeteors)))
2174 2169
2175 2170 #Meteor array
2176 2171 # arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
2177 2172 # arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
2178 2173
2179 2174 #Parameters Array
2180 2175 arrayParameters[:,0] = arrayDate #Date
2181 2176 arrayParameters[:,1] = heiRang[arrayMeteors[:,0].astype(int)] #Range
2182 2177 arrayParameters[:,6:8] = arrayMeteors[:,-3:-1] #Radial velocity and its error
2183 2178 arrayParameters[:,8:12] = arrayMeteors[:,7:11] #Phases
2184 2179 arrayParameters[:,-1] = arrayMeteors[:,-1] #Error
2185 2180
2186 2181
2187 2182 return arrayParameters
2188 2183
2189 2184 class CorrectSMPhases(Operation):
2190 2185
2191 2186 def run(self, dataOut, phaseOffsets, hmin = 50, hmax = 150, azimuth = 45, channelPositions = None):
2192 2187
2193 2188 arrayParameters = dataOut.data_param
2194 2189 pairsList = []
2195 2190 pairx = (0,1)
2196 2191 pairy = (2,3)
2197 2192 pairsList.append(pairx)
2198 2193 pairsList.append(pairy)
2199 2194 jph = numpy.zeros(4)
2200 2195
2201 2196 phaseOffsets = numpy.array(phaseOffsets)*numpy.pi/180
2202 2197 # arrayParameters[:,8:12] = numpy.unwrap(arrayParameters[:,8:12] + phaseOffsets)
2203 2198 arrayParameters[:,8:12] = numpy.angle(numpy.exp(1j*(arrayParameters[:,8:12] + phaseOffsets)))
2204 2199
2205 2200 meteorOps = SMOperations()
2206 2201 if channelPositions is None:
2207 2202 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2208 2203 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2209 2204
2210 2205 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2211 2206 h = (hmin,hmax)
2212 2207
2213 2208 arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, distances, jph)
2214 2209
2215 2210 dataOut.data_param = arrayParameters
2216 2211 return
2217 2212
2218 2213 class SMPhaseCalibration(Operation):
2219 2214
2220 2215 __buffer = None
2221 2216
2222 2217 __initime = None
2223 2218
2224 2219 __dataReady = False
2225 2220
2226 2221 __isConfig = False
2227 2222
2228 2223 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
2229 2224
2230 2225 dataTime = currentTime + paramInterval
2231 2226 deltaTime = dataTime - initTime
2232 2227
2233 2228 if deltaTime >= outputInterval or deltaTime < 0:
2234 2229 return True
2235 2230
2236 2231 return False
2237 2232
2238 2233 def __getGammas(self, pairs, d, phases):
2239 2234 gammas = numpy.zeros(2)
2240 2235
2241 2236 for i in range(len(pairs)):
2242 2237
2243 2238 pairi = pairs[i]
2244 2239
2245 2240 phip3 = phases[:,pairi[1]]
2246 2241 d3 = d[pairi[1]]
2247 2242 phip2 = phases[:,pairi[0]]
2248 2243 d2 = d[pairi[0]]
2249 2244 #Calculating gamma
2250 2245 # jdcos = alp1/(k*d1)
2251 2246 # jgamma = numpy.angle(numpy.exp(1j*(d0*alp1/d1 - alp0)))
2252 2247 jgamma = -phip2*d3/d2 - phip3
2253 2248 jgamma = numpy.angle(numpy.exp(1j*jgamma))
2254 2249 # jgamma[jgamma>numpy.pi] -= 2*numpy.pi
2255 2250 # jgamma[jgamma<-numpy.pi] += 2*numpy.pi
2256 2251
2257 2252 #Revised distribution
2258 2253 jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
2259 2254
2260 2255 #Histogram
2261 2256 nBins = 64.0
2262 2257 rmin = -0.5*numpy.pi
2263 2258 rmax = 0.5*numpy.pi
2264 2259 phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
2265 2260
2266 2261 meteorsY = phaseHisto[0]
2267 2262 phasesX = phaseHisto[1][:-1]
2268 2263 width = phasesX[1] - phasesX[0]
2269 2264 phasesX += width/2
2270 2265
2271 2266 #Gaussian aproximation
2272 2267 bpeak = meteorsY.argmax()
2273 2268 peak = meteorsY.max()
2274 2269 jmin = bpeak - 5
2275 2270 jmax = bpeak + 5 + 1
2276 2271
2277 2272 if jmin<0:
2278 2273 jmin = 0
2279 2274 jmax = 6
2280 2275 elif jmax > meteorsY.size:
2281 2276 jmin = meteorsY.size - 6
2282 2277 jmax = meteorsY.size
2283 2278
2284 2279 x0 = numpy.array([peak,bpeak,50])
2285 2280 coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
2286 2281
2287 2282 #Gammas
2288 2283 gammas[i] = coeff[0][1]
2289 2284
2290 2285 return gammas
2291 2286
2292 2287 def __residualFunction(self, coeffs, y, t):
2293 2288
2294 2289 return y - self.__gauss_function(t, coeffs)
2295 2290
2296 2291 def __gauss_function(self, t, coeffs):
2297 2292
2298 2293 return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
2299 2294
2300 2295 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
2301 2296 meteorOps = SMOperations()
2302 2297 nchan = 4
2303 2298 pairx = pairsList[0]
2304 2299 pairy = pairsList[1]
2305 2300 center_xangle = 0
2306 2301 center_yangle = 0
2307 2302 range_angle = numpy.array([10*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
2308 2303 ntimes = len(range_angle)
2309 2304
2310 2305 nstepsx = 20.0
2311 2306 nstepsy = 20.0
2312 2307
2313 2308 for iz in range(ntimes):
2314 2309 min_xangle = -range_angle[iz]/2 + center_xangle
2315 2310 max_xangle = range_angle[iz]/2 + center_xangle
2316 2311 min_yangle = -range_angle[iz]/2 + center_yangle
2317 2312 max_yangle = range_angle[iz]/2 + center_yangle
2318 2313
2319 2314 inc_x = (max_xangle-min_xangle)/nstepsx
2320 2315 inc_y = (max_yangle-min_yangle)/nstepsy
2321 2316
2322 2317 alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
2323 2318 alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
2324 2319 penalty = numpy.zeros((nstepsx,nstepsy))
2325 2320 jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
2326 2321 jph = numpy.zeros(nchan)
2327 2322
2328 2323 # Iterations looking for the offset
2329 2324 for iy in range(int(nstepsy)):
2330 2325 for ix in range(int(nstepsx)):
2331 2326 jph[pairy[1]] = alpha_y[iy]
2332 2327 jph[pairy[0]] = -gammas[1] - alpha_y[iy]*d[pairy[1]]/d[pairy[0]]
2333 2328
2334 2329 jph[pairx[1]] = alpha_x[ix]
2335 2330 jph[pairx[0]] = -gammas[0] - alpha_x[ix]*d[pairx[1]]/d[pairx[0]]
2336 2331
2337 2332 jph_array[:,ix,iy] = jph
2338 2333
2339 2334 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, d, jph)
2340 2335 error = meteorsArray1[:,-1]
2341 2336 ind1 = numpy.where(error==0)[0]
2342 2337 penalty[ix,iy] = ind1.size
2343 2338
2344 2339 i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
2345 2340 phOffset = jph_array[:,i,j]
2346 2341
2347 2342 center_xangle = phOffset[pairx[1]]
2348 2343 center_yangle = phOffset[pairy[1]]
2349 2344
2350 2345 phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
2351 2346 phOffset = phOffset*180/numpy.pi
2352 2347 return phOffset
2353 2348
2354 2349
2355 2350 def run(self, dataOut, hmin, hmax, channelPositions=None, nHours = 1):
2356 2351
2357 2352 dataOut.flagNoData = True
2358 2353 self.__dataReady = False
2359 2354 dataOut.outputInterval = nHours*3600
2360 2355
2361 2356 if self.__isConfig == False:
2362 2357 # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
2363 2358 #Get Initial LTC time
2364 2359 self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
2365 2360 self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
2366 2361
2367 2362 self.__isConfig = True
2368 2363
2369 2364 if self.__buffer is None:
2370 2365 self.__buffer = dataOut.data_param.copy()
2371 2366
2372 2367 else:
2373 2368 self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param))
2374 2369
2375 2370 self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
2376 2371
2377 2372 if self.__dataReady:
2378 2373 dataOut.utctimeInit = self.__initime
2379 2374 self.__initime += dataOut.outputInterval #to erase time offset
2380 2375
2381 2376 freq = dataOut.frequency
2382 2377 c = dataOut.C #m/s
2383 2378 lamb = c/freq
2384 2379 k = 2*numpy.pi/lamb
2385 2380 azimuth = 0
2386 2381 h = (hmin, hmax)
2387 2382 pairs = ((0,1),(2,3))
2388 2383
2389 2384 if channelPositions is None:
2390 2385 # channelPositions = [(2.5,0), (0,2.5), (0,0), (0,4.5), (-2,0)] #T
2391 2386 channelPositions = [(4.5,2), (2,4.5), (2,2), (2,0), (0,2)] #Estrella
2392 2387 meteorOps = SMOperations()
2393 2388 pairslist0, distances = meteorOps.getPhasePairs(channelPositions)
2394 2389
2395 2390 # distances1 = [-distances[0]*lamb, distances[1]*lamb, -distances[2]*lamb, distances[3]*lamb]
2396 2391
2397 2392 meteorsArray = self.__buffer
2398 2393 error = meteorsArray[:,-1]
2399 2394 boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
2400 2395 ind1 = numpy.where(boolError)[0]
2401 2396 meteorsArray = meteorsArray[ind1,:]
2402 2397 meteorsArray[:,-1] = 0
2403 2398 phases = meteorsArray[:,8:12]
2404 2399
2405 2400 #Calculate Gammas
2406 2401 gammas = self.__getGammas(pairs, distances, phases)
2407 2402 # gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
2408 2403 #Calculate Phases
2409 2404 phasesOff = self.__getPhases(azimuth, h, pairs, distances, gammas, meteorsArray)
2410 2405 phasesOff = phasesOff.reshape((1,phasesOff.size))
2411 2406 dataOut.data_output = -phasesOff
2412 2407 dataOut.flagNoData = False
2413 2408 dataOut.channelList = pairslist0
2414 2409 self.__buffer = None
2415 2410
2416 2411
2417 2412 return
2418 2413
2419 2414 class SMOperations():
2420 2415
2421 2416 def __init__(self):
2422 2417
2423 2418 return
2424 2419
2425 2420 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, distances, jph):
2426 2421
2427 2422 arrayParameters = arrayParameters0.copy()
2428 2423 hmin = h[0]
2429 2424 hmax = h[1]
2430 2425
2431 2426 #Calculate AOA (Error N 3, 4)
2432 2427 #JONES ET AL. 1998
2433 2428 AOAthresh = numpy.pi/8
2434 2429 error = arrayParameters[:,-1]
2435 2430 phases = -arrayParameters[:,8:12] + jph
2436 2431 # phases = numpy.unwrap(phases)
2437 2432 arrayParameters[:,3:6], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, distances, error, AOAthresh, azimuth)
2438 2433
2439 2434 #Calculate Heights (Error N 13 and 14)
2440 2435 error = arrayParameters[:,-1]
2441 2436 Ranges = arrayParameters[:,1]
2442 2437 zenith = arrayParameters[:,4]
2443 2438 arrayParameters[:,2], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
2444 2439
2445 2440 #----------------------- Get Final data ------------------------------------
2446 2441 # error = arrayParameters[:,-1]
2447 2442 # ind1 = numpy.where(error==0)[0]
2448 2443 # arrayParameters = arrayParameters[ind1,:]
2449 2444
2450 2445 return arrayParameters
2451 2446
2452 2447 def __getAOA(self, phases, pairsList, directions, error, AOAthresh, azimuth):
2453 2448
2454 2449 arrayAOA = numpy.zeros((phases.shape[0],3))
2455 2450 cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList,directions)
2456 2451
2457 2452 arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2458 2453 cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2459 2454 arrayAOA[:,2] = cosDirError
2460 2455
2461 2456 azimuthAngle = arrayAOA[:,0]
2462 2457 zenithAngle = arrayAOA[:,1]
2463 2458
2464 2459 #Setting Error
2465 2460 indError = numpy.where(numpy.logical_or(error == 3, error == 4))[0]
2466 2461 error[indError] = 0
2467 2462 #Number 3: AOA not fesible
2468 2463 indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
2469 2464 error[indInvalid] = 3
2470 2465 #Number 4: Large difference in AOAs obtained from different antenna baselines
2471 2466 indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
2472 2467 error[indInvalid] = 4
2473 2468 return arrayAOA, error
2474 2469
2475 2470 def __getDirectionCosines(self, arrayPhase, pairsList, distances):
2476 2471
2477 2472 #Initializing some variables
2478 2473 ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
2479 2474 ang_aux = ang_aux.reshape(1,ang_aux.size)
2480 2475
2481 2476 cosdir = numpy.zeros((arrayPhase.shape[0],2))
2482 2477 cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
2483 2478
2484 2479
2485 2480 for i in range(2):
2486 2481 ph0 = arrayPhase[:,pairsList[i][0]]
2487 2482 ph1 = arrayPhase[:,pairsList[i][1]]
2488 2483 d0 = distances[pairsList[i][0]]
2489 2484 d1 = distances[pairsList[i][1]]
2490 2485
2491 2486 ph0_aux = ph0 + ph1
2492 2487 ph0_aux = numpy.angle(numpy.exp(1j*ph0_aux))
2493 2488 # ph0_aux[ph0_aux > numpy.pi] -= 2*numpy.pi
2494 2489 # ph0_aux[ph0_aux < -numpy.pi] += 2*numpy.pi
2495 2490 #First Estimation
2496 2491 cosdir0[:,i] = (ph0_aux)/(2*numpy.pi*(d0 - d1))
2497 2492
2498 2493 #Most-Accurate Second Estimation
2499 2494 phi1_aux = ph0 - ph1
2500 2495 phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2501 2496 #Direction Cosine 1
2502 2497 cosdir1 = (phi1_aux + ang_aux)/(2*numpy.pi*(d0 + d1))
2503 2498
2504 2499 #Searching the correct Direction Cosine
2505 2500 cosdir0_aux = cosdir0[:,i]
2506 2501 cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
2507 2502 #Minimum Distance
2508 2503 cosDiff = (cosdir1 - cosdir0_aux)**2
2509 2504 indcos = cosDiff.argmin(axis = 1)
2510 2505 #Saving Value obtained
2511 2506 cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
2512 2507
2513 2508 return cosdir0, cosdir
2514 2509
2515 2510 def __calculateAOA(self, cosdir, azimuth):
2516 2511 cosdirX = cosdir[:,0]
2517 2512 cosdirY = cosdir[:,1]
2518 2513
2519 2514 zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2520 2515 azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth#0 deg north, 90 deg east
2521 2516 angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2522 2517
2523 2518 return angles
2524 2519
2525 2520 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
2526 2521
2527 2522 Ramb = 375 #Ramb = c/(2*PRF)
2528 2523 Re = 6371 #Earth Radius
2529 2524 heights = numpy.zeros(Ranges.shape)
2530 2525
2531 2526 R_aux = numpy.array([0,1,2])*Ramb
2532 2527 R_aux = R_aux.reshape(1,R_aux.size)
2533 2528
2534 2529 Ranges = Ranges.reshape(Ranges.size,1)
2535 2530
2536 2531 Ri = Ranges + R_aux
2537 2532 hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
2538 2533
2539 2534 #Check if there is a height between 70 and 110 km
2540 2535 h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
2541 2536 ind_h = numpy.where(h_bool == 1)[0]
2542 2537
2543 2538 hCorr = hi[ind_h, :]
2544 2539 ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
2545 2540
2546 2541 hCorr = hi[ind_hCorr]
2547 2542 heights[ind_h] = hCorr
2548 2543
2549 2544 #Setting Error
2550 2545 #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
2551 2546 #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
2552 2547 indError = numpy.where(numpy.logical_or(error == 13, error == 14))[0]
2553 2548 error[indError] = 0
2554 2549 indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
2555 2550 error[indInvalid2] = 14
2556 2551 indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2557 2552 error[indInvalid1] = 13
2558 2553
2559 2554 return heights, error
2560 2555
2561 2556 def getPhasePairs(self, channelPositions):
2562 2557 chanPos = numpy.array(channelPositions)
2563 2558 listOper = list(itertools.combinations(range(5),2))
2564 2559
2565 2560 distances = numpy.zeros(4)
2566 2561 axisX = []
2567 2562 axisY = []
2568 2563 distX = numpy.zeros(3)
2569 2564 distY = numpy.zeros(3)
2570 2565 ix = 0
2571 2566 iy = 0
2572 2567
2573 2568 pairX = numpy.zeros((2,2))
2574 2569 pairY = numpy.zeros((2,2))
2575 2570
2576 2571 for i in range(len(listOper)):
2577 2572 pairi = listOper[i]
2578 2573
2579 2574 posDif = numpy.abs(chanPos[pairi[0],:] - chanPos[pairi[1],:])
2580 2575
2581 2576 if posDif[0] == 0:
2582 2577 axisY.append(pairi)
2583 2578 distY[iy] = posDif[1]
2584 2579 iy += 1
2585 2580 elif posDif[1] == 0:
2586 2581 axisX.append(pairi)
2587 2582 distX[ix] = posDif[0]
2588 2583 ix += 1
2589 2584
2590 2585 for i in range(2):
2591 2586 if i==0:
2592 2587 dist0 = distX
2593 2588 axis0 = axisX
2594 2589 else:
2595 2590 dist0 = distY
2596 2591 axis0 = axisY
2597 2592
2598 2593 side = numpy.argsort(dist0)[:-1]
2599 2594 axis0 = numpy.array(axis0)[side,:]
2600 2595 chanC = int(numpy.intersect1d(axis0[0,:], axis0[1,:])[0])
2601 2596 axis1 = numpy.unique(numpy.reshape(axis0,4))
2602 2597 side = axis1[axis1 != chanC]
2603 2598 diff1 = chanPos[chanC,i] - chanPos[side[0],i]
2604 2599 diff2 = chanPos[chanC,i] - chanPos[side[1],i]
2605 2600 if diff1<0:
2606 2601 chan2 = side[0]
2607 2602 d2 = numpy.abs(diff1)
2608 2603 chan1 = side[1]
2609 2604 d1 = numpy.abs(diff2)
2610 2605 else:
2611 2606 chan2 = side[1]
2612 2607 d2 = numpy.abs(diff2)
2613 2608 chan1 = side[0]
2614 2609 d1 = numpy.abs(diff1)
2615 2610
2616 2611 if i==0:
2617 2612 chanCX = chanC
2618 2613 chan1X = chan1
2619 2614 chan2X = chan2
2620 2615 distances[0:2] = numpy.array([d1,d2])
2621 2616 else:
2622 2617 chanCY = chanC
2623 2618 chan1Y = chan1
2624 2619 chan2Y = chan2
2625 2620 distances[2:4] = numpy.array([d1,d2])
2626 2621 # axisXsides = numpy.reshape(axisX[ix,:],4)
2627 2622 #
2628 2623 # channelCentX = int(numpy.intersect1d(pairX[0,:], pairX[1,:])[0])
2629 2624 # channelCentY = int(numpy.intersect1d(pairY[0,:], pairY[1,:])[0])
2630 2625 #
2631 2626 # ind25X = numpy.where(pairX[0,:] != channelCentX)[0][0]
2632 2627 # ind20X = numpy.where(pairX[1,:] != channelCentX)[0][0]
2633 2628 # channel25X = int(pairX[0,ind25X])
2634 2629 # channel20X = int(pairX[1,ind20X])
2635 2630 # ind25Y = numpy.where(pairY[0,:] != channelCentY)[0][0]
2636 2631 # ind20Y = numpy.where(pairY[1,:] != channelCentY)[0][0]
2637 2632 # channel25Y = int(pairY[0,ind25Y])
2638 2633 # channel20Y = int(pairY[1,ind20Y])
2639 2634
2640 2635 # pairslist = [(channelCentX, channel25X),(channelCentX, channel20X),(channelCentY,channel25Y),(channelCentY, channel20Y)]
2641 2636 pairslist = [(chanCX, chan1X),(chanCX, chan2X),(chanCY,chan1Y),(chanCY, chan2Y)]
2642 2637
2643 2638 return pairslist, distances
2644 2639 # def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
2645 2640 #
2646 2641 # arrayAOA = numpy.zeros((phases.shape[0],3))
2647 2642 # cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
2648 2643 #
2649 2644 # arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
2650 2645 # cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
2651 2646 # arrayAOA[:,2] = cosDirError
2652 2647 #
2653 2648 # azimuthAngle = arrayAOA[:,0]
2654 2649 # zenithAngle = arrayAOA[:,1]
2655 2650 #
2656 2651 # #Setting Error
2657 2652 # #Number 3: AOA not fesible
2658 2653 # indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
2659 2654 # error[indInvalid] = 3
2660 2655 # #Number 4: Large difference in AOAs obtained from different antenna baselines
2661 2656 # indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
2662 2657 # error[indInvalid] = 4
2663 2658 # return arrayAOA, error
2664 2659 #
2665 2660 # def __getDirectionCosines(self, arrayPhase, pairsList):
2666 2661 #
2667 2662 # #Initializing some variables
2668 2663 # ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
2669 2664 # ang_aux = ang_aux.reshape(1,ang_aux.size)
2670 2665 #
2671 2666 # cosdir = numpy.zeros((arrayPhase.shape[0],2))
2672 2667 # cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
2673 2668 #
2674 2669 #
2675 2670 # for i in range(2):
2676 2671 # #First Estimation
2677 2672 # phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
2678 2673 # #Dealias
2679 2674 # indcsi = numpy.where(phi0_aux > numpy.pi)
2680 2675 # phi0_aux[indcsi] -= 2*numpy.pi
2681 2676 # indcsi = numpy.where(phi0_aux < -numpy.pi)
2682 2677 # phi0_aux[indcsi] += 2*numpy.pi
2683 2678 # #Direction Cosine 0
2684 2679 # cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
2685 2680 #
2686 2681 # #Most-Accurate Second Estimation
2687 2682 # phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
2688 2683 # phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
2689 2684 # #Direction Cosine 1
2690 2685 # cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
2691 2686 #
2692 2687 # #Searching the correct Direction Cosine
2693 2688 # cosdir0_aux = cosdir0[:,i]
2694 2689 # cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
2695 2690 # #Minimum Distance
2696 2691 # cosDiff = (cosdir1 - cosdir0_aux)**2
2697 2692 # indcos = cosDiff.argmin(axis = 1)
2698 2693 # #Saving Value obtained
2699 2694 # cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
2700 2695 #
2701 2696 # return cosdir0, cosdir
2702 2697 #
2703 2698 # def __calculateAOA(self, cosdir, azimuth):
2704 2699 # cosdirX = cosdir[:,0]
2705 2700 # cosdirY = cosdir[:,1]
2706 2701 #
2707 2702 # zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
2708 2703 # azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
2709 2704 # angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
2710 2705 #
2711 2706 # return angles
2712 2707 #
2713 2708 # def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
2714 2709 #
2715 2710 # Ramb = 375 #Ramb = c/(2*PRF)
2716 2711 # Re = 6371 #Earth Radius
2717 2712 # heights = numpy.zeros(Ranges.shape)
2718 2713 #
2719 2714 # R_aux = numpy.array([0,1,2])*Ramb
2720 2715 # R_aux = R_aux.reshape(1,R_aux.size)
2721 2716 #
2722 2717 # Ranges = Ranges.reshape(Ranges.size,1)
2723 2718 #
2724 2719 # Ri = Ranges + R_aux
2725 2720 # hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
2726 2721 #
2727 2722 # #Check if there is a height between 70 and 110 km
2728 2723 # h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
2729 2724 # ind_h = numpy.where(h_bool == 1)[0]
2730 2725 #
2731 2726 # hCorr = hi[ind_h, :]
2732 2727 # ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
2733 2728 #
2734 2729 # hCorr = hi[ind_hCorr]
2735 2730 # heights[ind_h] = hCorr
2736 2731 #
2737 2732 # #Setting Error
2738 2733 # #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
2739 2734 # #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
2740 2735 #
2741 2736 # indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
2742 2737 # error[indInvalid2] = 14
2743 2738 # indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
2744 2739 # error[indInvalid1] = 13
2745 2740 #
2746 2741 # return heights, error
Show file before | Show file after | Hide comments
@@ -1,393 +1,420
1 1 '''
2 2 @author: Juan C. Espinoza
3 3 '''
4 4
5 5 import time
6 6 import json
7 7 import numpy
8 8 import paho.mqtt.client as mqtt
9 9 import zmq
10 10 import cPickle as pickle
11 11 import datetime
12 12 from zmq.utils.monitor import recv_monitor_message
13 13 from functools import wraps
14 14 from threading import Thread
15 15 from multiprocessing import Process
16 16
17 17 from schainpy.model.proc.jroproc_base import Operation, ProcessingUnit
18 18
19 19 MAXNUMX = 100
20 20 MAXNUMY = 100
21 21
22 22 class PrettyFloat(float):
23 23 def __repr__(self):
24 24 return '%.2f' % self
25 25
26 26 def roundFloats(obj):
27 27 if isinstance(obj, list):
28 28 return map(roundFloats, obj)
29 29 elif isinstance(obj, float):
30 30 return round(obj, 2)
31 31
32 32 def decimate(z):
33 33 # dx = int(len(self.x)/self.__MAXNUMX) + 1
34 34 dy = int(len(z[0])/MAXNUMY) + 1
35 35 return z[::, ::dy]
36 36
37 37 class throttle(object):
38 38 """Decorator that prevents a function from being called more than once every
39 39 time period.
40 40 To create a function that cannot be called more than once a minute, but
41 41 will sleep until it can be called:
42 42 @throttle(minutes=1)
43 43 def foo():
44 44 pass
45 45
46 46 for i in range(10):
47 47 foo()
48 48 print "This function has run %s times." % i
49 49 """
50 50
51 51 def __init__(self, seconds=0, minutes=0, hours=0):
52 52 self.throttle_period = datetime.timedelta(
53 53 seconds=seconds, minutes=minutes, hours=hours
54 54 )
55 55
56 56 self.time_of_last_call = datetime.datetime.min
57 57
58 58 def __call__(self, fn):
59 59 @wraps(fn)
60 60 def wrapper(*args, **kwargs):
61 61 now = datetime.datetime.now()
62 62 time_since_last_call = now - self.time_of_last_call
63 63 time_left = self.throttle_period - time_since_last_call
64 64
65 65 if time_left > datetime.timedelta(seconds=0):
66 66 return
67 67
68 68 self.time_of_last_call = datetime.datetime.now()
69 69 return fn(*args, **kwargs)
70 70
71 71 return wrapper
72 72
73 73
74 74 class PublishData(Operation):
75 75 """Clase publish."""
76 76
77 77 def __init__(self, **kwargs):
78 78 """Inicio."""
79 79 Operation.__init__(self, **kwargs)
80 80 self.isConfig = False
81 81 self.client = None
82 82 self.zeromq = None
83 83 self.mqtt = None
84 84
85 85 def on_disconnect(self, client, userdata, rc):
86 86 if rc != 0:
87 87 print("Unexpected disconnection.")
88 88 self.connect()
89 89
90 90 def connect(self):
91 91 print 'trying to connect'
92 92 try:
93 93 self.client.connect(
94 94 host=self.host,
95 95 port=self.port,
96 96 keepalive=60*10,
97 97 bind_address='')
98 98 self.client.loop_start()
99 99 # self.client.publish(
100 100 # self.topic + 'SETUP',
101 101 # json.dumps(setup),
102 102 # retain=True
103 103 # )
104 104 except:
105 105 print "MQTT Conection error."
106 106 self.client = False
107 107
108 108 def setup(self, port=1883, username=None, password=None, clientId="user", zeromq=1, **kwargs):
109 109 self.counter = 0
110 110 self.topic = kwargs.get('topic', 'schain')
111 111 self.delay = kwargs.get('delay', 0)
112 112 self.plottype = kwargs.get('plottype', 'spectra')
113 113 self.host = kwargs.get('host', "10.10.10.82")
114 114 self.port = kwargs.get('port', 3000)
115 115 self.clientId = clientId
116 116 self.cnt = 0
117 117 self.zeromq = zeromq
118 118 self.mqtt = kwargs.get('plottype', 0)
119 119 self.client = None
120 120 setup = []
121 121 if mqtt is 1:
122 122 self.client = mqtt.Client(
123 123 client_id=self.clientId + self.topic + 'SCHAIN',
124 124 clean_session=True)
125 125 self.client.on_disconnect = self.on_disconnect
126 126 self.connect()
127 127 for plot in self.plottype:
128 128 setup.append({
129 129 'plot': plot,
130 130 'topic': self.topic + plot,
131 131 'title': getattr(self, plot + '_' + 'title', False),
132 132 'xlabel': getattr(self, plot + '_' + 'xlabel', False),
133 133 'ylabel': getattr(self, plot + '_' + 'ylabel', False),
134 134 'xrange': getattr(self, plot + '_' + 'xrange', False),
135 135 'yrange': getattr(self, plot + '_' + 'yrange', False),
136 136 'zrange': getattr(self, plot + '_' + 'zrange', False),
137 137 })
138 138 if zeromq is 1:
139 139 context = zmq.Context()
140 140 self.zmq_socket = context.socket(zmq.PUSH)
141 141 server = kwargs.get('server', 'zmq.pipe')
142 142
143 143 if 'tcp://' in server:
144 144 address = server
145 145 else:
146 146 address = 'ipc:///tmp/%s' % server
147 147
148 148 self.zmq_socket.connect(address)
149 149 time.sleep(1)
150 150
151 151
152 152
153 153 def publish_data(self):
154 154 self.dataOut.finished = False
155 155 if self.mqtt is 1:
156 156 yData = self.dataOut.heightList[:2].tolist()
157 157 if self.plottype == 'spectra':
158 158 data = getattr(self.dataOut, 'data_spc')
159 159 z = data/self.dataOut.normFactor
160 160 zdB = 10*numpy.log10(z)
161 161 xlen, ylen = zdB[0].shape
162 162 dx = int(xlen/MAXNUMX) + 1
163 163 dy = int(ylen/MAXNUMY) + 1
164 164 Z = [0 for i in self.dataOut.channelList]
165 165 for i in self.dataOut.channelList:
166 166 Z[i] = zdB[i][::dx, ::dy].tolist()
167 167 payload = {
168 168 'timestamp': self.dataOut.utctime,
169 169 'data': roundFloats(Z),
170 170 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
171 171 'interval': self.dataOut.getTimeInterval(),
172 172 'type': self.plottype,
173 173 'yData': yData
174 174 }
175 175 # print payload
176 176
177 177 elif self.plottype in ('rti', 'power'):
178 178 data = getattr(self.dataOut, 'data_spc')
179 179 z = data/self.dataOut.normFactor
180 180 avg = numpy.average(z, axis=1)
181 181 avgdB = 10*numpy.log10(avg)
182 182 xlen, ylen = z[0].shape
183 183 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
184 184 AVG = [0 for i in self.dataOut.channelList]
185 185 for i in self.dataOut.channelList:
186 186 AVG[i] = avgdB[i][::dy].tolist()
187 187 payload = {
188 188 'timestamp': self.dataOut.utctime,
189 189 'data': roundFloats(AVG),
190 190 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
191 191 'interval': self.dataOut.getTimeInterval(),
192 192 'type': self.plottype,
193 193 'yData': yData
194 194 }
195 195 elif self.plottype == 'noise':
196 196 noise = self.dataOut.getNoise()/self.dataOut.normFactor
197 197 noisedB = 10*numpy.log10(noise)
198 198 payload = {
199 199 'timestamp': self.dataOut.utctime,
200 200 'data': roundFloats(noisedB.reshape(-1, 1).tolist()),
201 201 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
202 202 'interval': self.dataOut.getTimeInterval(),
203 203 'type': self.plottype,
204 204 'yData': yData
205 205 }
206 206 elif self.plottype == 'snr':
207 207 data = getattr(self.dataOut, 'data_SNR')
208 208 avgdB = 10*numpy.log10(data)
209 209
210 210 ylen = data[0].size
211 211 dy = numpy.floor(ylen/self.__MAXNUMY) + 1
212 212 AVG = [0 for i in self.dataOut.channelList]
213 213 for i in self.dataOut.channelList:
214 214 AVG[i] = avgdB[i][::dy].tolist()
215 215 payload = {
216 216 'timestamp': self.dataOut.utctime,
217 217 'data': roundFloats(AVG),
218 218 'channels': ['Ch %s' % ch for ch in self.dataOut.channelList],
219 219 'type': self.plottype,
220 220 'yData': yData
221 221 }
222 222 else:
223 223 print "Tipo de grafico invalido"
224 224 payload = {
225 225 'data': 'None',
226 226 'timestamp': 'None',
227 227 'type': None
228 228 }
229 229 # print 'Publishing data to {}'.format(self.host)
230 230 self.client.publish(self.topic + self.plottype, json.dumps(payload), qos=0)
231 231
232 232 if self.zeromq is 1:
233 233 print '[Sending] {} - {}'.format(self.dataOut.type, self.dataOut.datatime)
234 234 self.zmq_socket.send_pyobj(self.dataOut)
235 235
236 236 def run(self, dataOut, **kwargs):
237 237 self.dataOut = dataOut
238 238 if not self.isConfig:
239 239 self.setup(**kwargs)
240 240 self.isConfig = True
241 241
242 242 self.publish_data()
243 243 time.sleep(self.delay)
244 244
245 245 def close(self):
246 246 if self.zeromq is 1:
247 247 self.dataOut.finished = True
248 248 self.zmq_socket.send_pyobj(self.dataOut)
249 249
250 250 if self.client:
251 251 self.client.loop_stop()
252 252 self.client.disconnect()
253 253
254 254
255 255 class ReceiverData(ProcessingUnit, Process):
256 256
257 257 throttle_value = 5
258 258
259 259 def __init__(self, **kwargs):
260 260
261 261 ProcessingUnit.__init__(self, **kwargs)
262 262 Process.__init__(self)
263 263 self.mp = False
264 264 self.isConfig = False
265 self.isWebConfig = False
265 266 self.plottypes =[]
266 267 self.connections = 0
267 268 server = kwargs.get('server', 'zmq.pipe')
269 plot_server = kwargs.get('plot_server', 'zmq.web')
268 270 if 'tcp://' in server:
269 271 address = server
270 272 else:
271 273 address = 'ipc:///tmp/%s' % server
272 274
275 if 'tcp://' in plot_server:
276 plot_address = plot_server
277 else:
278 plot_address = 'ipc:///tmp/%s' % plot_server
279
273 280 self.address = address
281 self.plot_address = plot_address
274 282 self.plottypes = [s.strip() for s in kwargs.get('plottypes', 'rti').split(',')]
275 283 self.realtime = kwargs.get('realtime', False)
276 284 self.throttle_value = kwargs.get('throttle', 10)
277 285 self.sendData = self.initThrottle(self.throttle_value)
278 286 self.setup()
279 287
280 288 def setup(self):
281 289
282 290 self.data = {}
283 291 self.data['times'] = []
284 292 for plottype in self.plottypes:
285 293 self.data[plottype] = {}
286 294 self.data['noise'] = {}
287 295 self.data['throttle'] = self.throttle_value
288 296 self.data['ENDED'] = False
289 297 self.isConfig = True
290 298 self.data_web = {}
291 299
292 300 def event_monitor(self, monitor):
293 301
294 302 events = {}
295 303
296 304 for name in dir(zmq):
297 305 if name.startswith('EVENT_'):
298 306 value = getattr(zmq, name)
299 307 events[value] = name
300 308
301 309 while monitor.poll():
302 310 evt = recv_monitor_message(monitor)
303 311 if evt['event'] == 32:
304 312 self.connections += 1
305 313 if evt['event'] == 512:
306 314 pass
307 315 if self.connections == 0 and self.started is True:
308 316 self.ended = True
309 317 # send('ENDED')
310 318 evt.update({'description': events[evt['event']]})
311 319
312 320 if evt['event'] == zmq.EVENT_MONITOR_STOPPED:
313 321 break
314 322 monitor.close()
315 323 print("event monitor thread done!")
316 324
317 325 def initThrottle(self, throttle_value):
318 326
319 327 @throttle(seconds=throttle_value)
320 328 def sendDataThrottled(fn_sender, data):
321 329 fn_sender(data)
322 330
323 331 return sendDataThrottled
324 332
325 333 def send(self, data):
326 334 # print '[sending] data=%s size=%s' % (data.keys(), len(data['times']))
327 335 self.sender.send_pyobj(data)
328 336
329 337 def update(self):
330 338
331 339 t = self.dataOut.ltctime
332 340 self.data['times'].append(t)
333 341 self.data['dataOut'] = self.dataOut
334 342 for plottype in self.plottypes:
335 343 if plottype == 'spc':
336 344 z = self.dataOut.data_spc/self.dataOut.normFactor
337 345 self.data[plottype] = 10*numpy.log10(z)
338 346 self.data['noise'][t] = 10*numpy.log10(self.dataOut.getNoise()/self.dataOut.normFactor)
339 347 if plottype == 'rti':
340 348 self.data[plottype][t] = self.dataOut.getPower()
341 349 if plottype == 'snr':
342 350 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_SNR)
343 351 if plottype == 'dop':
344 352 self.data[plottype][t] = 10*numpy.log10(self.dataOut.data_DOP)
345 353 if plottype == 'coh':
346 354 self.data[plottype][t] = self.dataOut.getCoherence()
347 355 if plottype == 'phase':
348 356 self.data[plottype][t] = self.dataOut.getCoherence(phase=True)
349 357 if self.realtime:
350 358 self.data_web[plottype] = roundFloats(decimate(self.data[plottype][t]).tolist())
351 359 self.data_web['time'] = t
352 360 def run(self):
353 361
354 362 print '[Starting] {} from {}'.format(self.name, self.address)
355 363
356 364 self.context = zmq.Context()
357 365 self.receiver = self.context.socket(zmq.PULL)
358 366 self.receiver.bind(self.address)
359 367 monitor = self.receiver.get_monitor_socket()
360 368 self.sender = self.context.socket(zmq.PUB)
361 if self.realtime:
362 self.sender_web = self.context.socket(zmq.PUB)
363 # self.sender_web.setsockopt(zmq.PUBLISH, 'realtime')
364 self.sender_web.bind("ipc:///tmp/zmq.web")
369 if self.realtime:
370 self.sender_web = self.context.socket(zmq.PUB)
371 self.sender_web.bind(self.plot_address)
365 372 self.sender.bind("ipc:///tmp/zmq.plots")
366 373
367 374 t = Thread(target=self.event_monitor, args=(monitor,))
368 375 t.start()
369 376
370 377 while True:
371 378 self.dataOut = self.receiver.recv_pyobj()
372 379 # print '[Receiving] {} - {}'.format(self.dataOut.type,
373 380 # self.dataOut.datatime.ctime())
374 381
375 382 self.update()
376 383
377 384 if self.dataOut.finished is True:
378 385 self.send(self.data)
379 386 self.connections -= 1
380 387 if self.connections == 0 and self.started:
381 388 self.ended = True
382 389 self.data['ENDED'] = True
383 390 self.send(self.data)
384 391 self.setup()
385 392 else:
386 393 if self.realtime:
387 394 self.send(self.data)
388 395 self.sender_web.send_string(json.dumps(self.data_web))
389 396 else:
390 397 self.sendData(self.send, self.data)
391 398 self.started = True
392 399
393 400 return
401
402 def sendToWeb(self):
403
404 if not self.isWebConfig:
405 context = zmq.Context()
406 sender_web_config = context.socket(zmq.PUB)
407 if 'tcp://' in self.plot_address:
408 print self.plot_address
409 dum, address, port = self.plot_address.split(':')
410 conf_address = '{}:{}:{}'.format(dum, address, int(port)+1)
411 else:
412 conf_address = self.plot_address + '.config'
413 sender_web_config.bind(conf_address)
414
415 for kwargs in self.operationKwargs.values():
416 if 'plot' in kwargs:
417 sender_web_config.send_string(json.dumps(kwargs))
418 print kwargs
419 self.isWebConfig = True
420
General Comments 0
You need to be logged in to leave comments. Login now