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