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