##// END OF EJS Templates
schain
RSS
LP Faraday update
rflores -
r1542:9f6529c5475c
parent child
Show More
Side by Side Unified
Context file:
r1542:9f6529c5475c Loading diff...:

The requested changes are too big and content was truncated. Show full diff

Show file before | Show file after | Hide comments
@@ -1,668 +1,669
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """API to create signal chain projects
6 6
7 7 The API is provide through class: Project
8 8 """
9 9
10 10 import re
11 11 import sys
12 12 import ast
13 13 import datetime
14 14 import traceback
15 15 import time
16 16 import multiprocessing
17 17 from multiprocessing import Process, Queue
18 18 from threading import Thread
19 19 from xml.etree.ElementTree import ElementTree, Element, SubElement
20 20
21 21 from schainpy.admin import Alarm, SchainWarning
22 22 from schainpy.model import *
23 23 from schainpy.utils import log
24 24
25 25 if 'darwin' in sys.platform and sys.version_info[0] == 3 and sys.version_info[1] > 7:
26 26 multiprocessing.set_start_method('fork')
27 27
28 28 class ConfBase():
29 29
30 30 def __init__(self):
31 31
32 32 self.id = '0'
33 33 self.name = None
34 34 self.priority = None
35 35 self.parameters = {}
36 36 self.object = None
37 37 self.operations = []
38 38
39 39 def getId(self):
40 40
41 41 return self.id
42 42
43 43 def getNewId(self):
44 44
45 45 return int(self.id) * 10 + len(self.operations) + 1
46 46
47 47 def updateId(self, new_id):
48 48
49 49 self.id = str(new_id)
50 50
51 51 n = 1
52 52 for conf in self.operations:
53 53 conf_id = str(int(new_id) * 10 + n)
54 54 conf.updateId(conf_id)
55 55 n += 1
56 56
57 57 def getKwargs(self):
58 58
59 59 params = {}
60 60
61 61 for key, value in self.parameters.items():
62 62 if value not in (None, '', ' '):
63 63 params[key] = value
64 64
65 65 return params
66 66
67 67 def update(self, **kwargs):
68 68
69 69 for key, value in kwargs.items():
70 70 self.addParameter(name=key, value=value)
71 71
72 72 def addParameter(self, name, value, format=None):
73 73 '''
74 74 '''
75 75
76 76 if isinstance(value, str) and re.search(r'(\d+/\d+/\d+)', value):
77 77 self.parameters[name] = datetime.date(*[int(x) for x in value.split('/')])
78 78 elif isinstance(value, str) and re.search(r'(\d+:\d+:\d+)', value):
79 79 self.parameters[name] = datetime.time(*[int(x) for x in value.split(':')])
80 80 else:
81 81 try:
82 82 self.parameters[name] = ast.literal_eval(value)
83 83 except:
84 84 if isinstance(value, str) and ',' in value:
85 85 self.parameters[name] = value.split(',')
86 86 else:
87 87 self.parameters[name] = value
88 88
89 89 def getParameters(self):
90 90
91 91 params = {}
92 92 for key, value in self.parameters.items():
93 93 s = type(value).__name__
94 94 if s == 'date':
95 95 params[key] = value.strftime('%Y/%m/%d')
96 96 elif s == 'time':
97 97 params[key] = value.strftime('%H:%M:%S')
98 98 else:
99 99 params[key] = str(value)
100 100
101 101 return params
102 102
103 103 def makeXml(self, element):
104 104
105 105 xml = SubElement(element, self.ELEMENTNAME)
106 106 for label in self.xml_labels:
107 107 xml.set(label, str(getattr(self, label)))
108 108
109 109 for key, value in self.getParameters().items():
110 110 xml_param = SubElement(xml, 'Parameter')
111 111 xml_param.set('name', key)
112 112 xml_param.set('value', value)
113 113
114 114 for conf in self.operations:
115 115 conf.makeXml(xml)
116 116
117 117 def __str__(self):
118 118
119 119 if self.ELEMENTNAME == 'Operation':
120 120 s = ' {}[id={}]\n'.format(self.name, self.id)
121 121 else:
122 122 s = '{}[id={}, inputId={}]\n'.format(self.name, self.id, self.inputId)
123 123
124 124 for key, value in self.parameters.items():
125 125 if self.ELEMENTNAME == 'Operation':
126 126 s += ' {}: {}\n'.format(key, value)
127 127 else:
128 128 s += ' {}: {}\n'.format(key, value)
129 129
130 130 for conf in self.operations:
131 131 s += str(conf)
132 132
133 133 return s
134 134
135 135 class OperationConf(ConfBase):
136 136
137 137 ELEMENTNAME = 'Operation'
138 138 xml_labels = ['id', 'name']
139 139
140 140 def setup(self, id, name, priority, project_id, err_queue):
141 141
142 142 self.id = str(id)
143 143 self.project_id = project_id
144 144 self.name = name
145 145 self.type = 'other'
146 146 self.err_queue = err_queue
147 147
148 148 def readXml(self, element, project_id, err_queue):
149 149
150 150 self.id = element.get('id')
151 151 self.name = element.get('name')
152 152 self.type = 'other'
153 153 self.project_id = str(project_id)
154 154 self.err_queue = err_queue
155 155
156 156 for elm in element.iter('Parameter'):
157 157 self.addParameter(elm.get('name'), elm.get('value'))
158 158
159 159 def createObject(self):
160 160
161 161 className = eval(self.name)
162 162
163 163 if 'Plot' in self.name or 'Writer' in self.name or 'Send' in self.name or 'print' in self.name:
164 164 kwargs = self.getKwargs()
165 165 opObj = className(self.id, self.id, self.project_id, self.err_queue, **kwargs)
166 166 opObj.start()
167 167 self.type = 'external'
168 168 else:
169 169 opObj = className()
170 170
171 171 self.object = opObj
172 172 return opObj
173 173
174 174 class ProcUnitConf(ConfBase):
175 175
176 176 ELEMENTNAME = 'ProcUnit'
177 177 xml_labels = ['id', 'inputId', 'name']
178 178
179 179 def setup(self, project_id, id, name, datatype, inputId, err_queue):
180 180 '''
181 181 '''
182 182
183 183 if datatype == None and name == None:
184 184 raise ValueError('datatype or name should be defined')
185 185
186 186 if name == None:
187 187 if 'Proc' in datatype:
188 188 name = datatype
189 189 else:
190 190 name = '%sProc' % (datatype)
191 191
192 192 if datatype == None:
193 193 datatype = name.replace('Proc', '')
194 194
195 195 self.id = str(id)
196 196 self.project_id = project_id
197 197 self.name = name
198 198 self.datatype = datatype
199 199 self.inputId = inputId
200 200 self.err_queue = err_queue
201 201 self.operations = []
202 202 self.parameters = {}
203 203
204 204 def removeOperation(self, id):
205 205
206 206 i = [1 if x.id == id else 0 for x in self.operations]
207 207 self.operations.pop(i.index(1))
208 208
209 209 def getOperation(self, id):
210 210
211 211 for conf in self.operations:
212 212 if conf.id == id:
213 213 return conf
214 214
215 215 def addOperation(self, name, optype='self'):
216 216 '''
217 217 '''
218 218
219 219 id = self.getNewId()
220 220 conf = OperationConf()
221 221 conf.setup(id, name=name, priority='0', project_id=self.project_id, err_queue=self.err_queue)
222 222 self.operations.append(conf)
223 223
224 224 return conf
225 225
226 226 def readXml(self, element, project_id, err_queue):
227 227
228 228 self.id = element.get('id')
229 229 self.name = element.get('name')
230 230 self.inputId = None if element.get('inputId') == 'None' else element.get('inputId')
231 231 self.datatype = element.get('datatype', self.name.replace(self.ELEMENTNAME.replace('Unit', ''), ''))
232 232 self.project_id = str(project_id)
233 233 self.err_queue = err_queue
234 234 self.operations = []
235 235 self.parameters = {}
236 236
237 237 for elm in element:
238 238 if elm.tag == 'Parameter':
239 239 self.addParameter(elm.get('name'), elm.get('value'))
240 240 elif elm.tag == 'Operation':
241 241 conf = OperationConf()
242 242 conf.readXml(elm, project_id, err_queue)
243 243 self.operations.append(conf)
244 244
245 245 def createObjects(self):
246 246 '''
247 247 Instancia de unidades de procesamiento.
248 248 '''
249 249
250 250 className = eval(self.name)
251 251 kwargs = self.getKwargs()
252 252 procUnitObj = className()
253 253 procUnitObj.name = self.name
254 254 log.success('creating process...', self.name)
255 255
256 256 for conf in self.operations:
257 257
258 258 opObj = conf.createObject()
259 259
260 260 log.success('adding operation: {}, type:{}'.format(
261 261 conf.name,
262 262 conf.type), self.name)
263 263
264 264 procUnitObj.addOperation(conf, opObj)
265 265
266 266 self.object = procUnitObj
267 267
268 268 def run(self):
269 269 '''
270 270 '''
271 271 #self.object.call(**self.getKwargs())
272 272
273 273 return self.object.call(**self.getKwargs())
274 274
275 275
276 276 class ReadUnitConf(ProcUnitConf):
277 277
278 278 ELEMENTNAME = 'ReadUnit'
279 279
280 280 def __init__(self):
281 281
282 282 self.id = None
283 283 self.datatype = None
284 284 self.name = None
285 285 self.inputId = None
286 286 self.operations = []
287 287 self.parameters = {}
288 288
289 289 def setup(self, project_id, id, name, datatype, err_queue, path='', startDate='', endDate='',
290 290 startTime='', endTime='', server=None, **kwargs):
291 291
292 292 if datatype == None and name == None:
293 293 raise ValueError('datatype or name should be defined')
294 294 if name == None:
295 295 if 'Reader' in datatype:
296 296 name = datatype
297 297 datatype = name.replace('Reader', '')
298 298 else:
299 299 name = '{}Reader'.format(datatype)
300 300 if datatype == None:
301 301 if 'Reader' in name:
302 302 datatype = name.replace('Reader', '')
303 303 else:
304 304 datatype = name
305 305 name = '{}Reader'.format(name)
306 306
307 307 self.id = id
308 308 self.project_id = project_id
309 309 self.name = name
310 310 self.datatype = datatype
311 311 self.err_queue = err_queue
312 312
313 313 self.addParameter(name='path', value=path)
314 314 self.addParameter(name='startDate', value=startDate)
315 315 self.addParameter(name='endDate', value=endDate)
316 316 self.addParameter(name='startTime', value=startTime)
317 317 self.addParameter(name='endTime', value=endTime)
318 318
319 319 for key, value in kwargs.items():
320 320 self.addParameter(name=key, value=value)
321 321
322 322
323 323 class Project(Process):
324 324 """API to create signal chain projects"""
325 325
326 326 ELEMENTNAME = 'Project'
327 327
328 328 def __init__(self, name=''):
329 329
330 330 Process.__init__(self)
331 331 self.id = '1'
332 332 if name:
333 333 self.name = '{} ({})'.format(Process.__name__, name)
334 334 self.filename = None
335 335 self.description = None
336 336 self.email = None
337 337 self.alarm = []
338 338 self.configurations = {}
339 339 # self.err_queue = Queue()
340 340 self.err_queue = None
341 341 self.started = False
342 342
343 343 def getNewId(self):
344 344
345 345 idList = list(self.configurations.keys())
346 346 id = int(self.id) * 10
347 347
348 348 while True:
349 349 id += 1
350 350
351 351 if str(id) in idList:
352 352 continue
353 353
354 354 break
355 355
356 356 return str(id)
357 357
358 358 def updateId(self, new_id):
359 359
360 360 self.id = str(new_id)
361 361
362 362 keyList = list(self.configurations.keys())
363 363 keyList.sort()
364 364
365 365 n = 1
366 366 new_confs = {}
367 367
368 368 for procKey in keyList:
369 369
370 370 conf = self.configurations[procKey]
371 371 idProcUnit = str(int(self.id) * 10 + n)
372 372 conf.updateId(idProcUnit)
373 373 new_confs[idProcUnit] = conf
374 374 n += 1
375 375
376 376 self.configurations = new_confs
377 377
378 378 def setup(self, id=1, name='', description='', email=None, alarm=[]):
379 379
380 380 self.id = str(id)
381 381 self.description = description
382 382 self.email = email
383 383 self.alarm = alarm
384 384 if name:
385 385 self.name = '{} ({})'.format(Process.__name__, name)
386 386
387 387 def update(self, **kwargs):
388 388
389 389 for key, value in kwargs.items():
390 390 setattr(self, key, value)
391 391
392 392 def clone(self):
393 393
394 394 p = Project()
395 395 p.id = self.id
396 396 p.name = self.name
397 397 p.description = self.description
398 398 p.configurations = self.configurations.copy()
399 399
400 400 return p
401 401
402 402 def addReadUnit(self, id=None, datatype=None, name=None, **kwargs):
403 403
404 404 '''
405 405 '''
406 406
407 407 if id is None:
408 408 idReadUnit = self.getNewId()
409 409 else:
410 410 idReadUnit = str(id)
411 411
412 412 conf = ReadUnitConf()
413 413 conf.setup(self.id, idReadUnit, name, datatype, self.err_queue, **kwargs)
414 414 self.configurations[conf.id] = conf
415 415
416 416 return conf
417 417
418 418 def addProcUnit(self, id=None, inputId='0', datatype=None, name=None):
419 419
420 420 '''
421 421 '''
422 422
423 423 if id is None:
424 424 idProcUnit = self.getNewId()
425 425 else:
426 426 idProcUnit = id
427 427
428 428 conf = ProcUnitConf()
429 429 conf.setup(self.id, idProcUnit, name, datatype, inputId, self.err_queue)
430 430 self.configurations[conf.id] = conf
431 431
432 432 return conf
433 433
434 434 def removeProcUnit(self, id):
435 435
436 436 if id in self.configurations:
437 437 self.configurations.pop(id)
438 438
439 439 def getReadUnit(self):
440 440
441 441 for obj in list(self.configurations.values()):
442 442 if obj.ELEMENTNAME == 'ReadUnit':
443 443 return obj
444 444
445 445 return None
446 446
447 447 def getProcUnit(self, id):
448 448
449 449 return self.configurations[id]
450 450
451 451 def getUnits(self):
452 452
453 453 keys = list(self.configurations)
454 454 keys.sort()
455 455
456 456 for key in keys:
457 457 yield self.configurations[key]
458 458
459 459 def updateUnit(self, id, **kwargs):
460 460
461 461 conf = self.configurations[id].update(**kwargs)
462 462
463 463 def makeXml(self):
464 464
465 465 xml = Element('Project')
466 466 xml.set('id', str(self.id))
467 467 xml.set('name', self.name)
468 468 xml.set('description', self.description)
469 469
470 470 for conf in self.configurations.values():
471 471 conf.makeXml(xml)
472 472
473 473 self.xml = xml
474 474
475 475 def writeXml(self, filename=None):
476 476
477 477 if filename == None:
478 478 if self.filename:
479 479 filename = self.filename
480 480 else:
481 481 filename = 'schain.xml'
482 482
483 483 if not filename:
484 484 print('filename has not been defined. Use setFilename(filename) for do it.')
485 485 return 0
486 486
487 487 abs_file = os.path.abspath(filename)
488 488
489 489 if not os.access(os.path.dirname(abs_file), os.W_OK):
490 490 print('No write permission on %s' % os.path.dirname(abs_file))
491 491 return 0
492 492
493 493 if os.path.isfile(abs_file) and not(os.access(abs_file, os.W_OK)):
494 494 print('File %s already exists and it could not be overwriten' % abs_file)
495 495 return 0
496 496
497 497 self.makeXml()
498 498
499 499 ElementTree(self.xml).write(abs_file, method='xml')
500 500
501 501 self.filename = abs_file
502 502
503 503 return 1
504 504
505 505 def readXml(self, filename):
506 506
507 507 abs_file = os.path.abspath(filename)
508 508
509 509 self.configurations = {}
510 510
511 511 try:
512 512 self.xml = ElementTree().parse(abs_file)
513 513 except:
514 514 log.error('Error reading %s, verify file format' % filename)
515 515 return 0
516 516
517 517 self.id = self.xml.get('id')
518 518 self.name = self.xml.get('name')
519 519 self.description = self.xml.get('description')
520 520
521 521 for element in self.xml:
522 522 if element.tag == 'ReadUnit':
523 523 conf = ReadUnitConf()
524 524 conf.readXml(element, self.id, self.err_queue)
525 525 self.configurations[conf.id] = conf
526 526 elif element.tag == 'ProcUnit':
527 527 conf = ProcUnitConf()
528 528 input_proc = self.configurations[element.get('inputId')]
529 529 conf.readXml(element, self.id, self.err_queue)
530 530 self.configurations[conf.id] = conf
531 531
532 532 self.filename = abs_file
533 533
534 534 return 1
535 535
536 536 def __str__(self):
537 537
538 538 text = '\nProject[id=%s, name=%s, description=%s]\n\n' % (
539 539 self.id,
540 540 self.name,
541 541 self.description,
542 542 )
543 543
544 544 for conf in self.configurations.values():
545 545 text += '{}'.format(conf)
546 546
547 547 return text
548 548
549 549 def createObjects(self):
550 550
551 551 keys = list(self.configurations.keys())
552 552 keys.sort()
553 553 for key in keys:
554 554 conf = self.configurations[key]
555 555 conf.createObjects()
556 556 if conf.inputId is not None:
557 557 if isinstance(conf.inputId, list):
558 558 conf.object.setInput([self.configurations[x].object for x in conf.inputId])
559 559 else:
560 560 conf.object.setInput([self.configurations[conf.inputId].object])
561 561
562 562 def monitor(self):
563 563
564 564 t = Thread(target=self._monitor, args=(self.err_queue, self.ctx))
565 565 t.start()
566 566
567 567 def _monitor(self, queue, ctx):
568 568
569 569 import socket
570 570
571 571 procs = 0
572 572 err_msg = ''
573 573
574 574 while True:
575 575 msg = queue.get()
576 576 if '#_start_#' in msg:
577 577 procs += 1
578 578 elif '#_end_#' in msg:
579 579 procs -= 1
580 580 else:
581 581 err_msg = msg
582 582
583 583 if procs == 0 or 'Traceback' in err_msg:
584 584 break
585 585 time.sleep(0.1)
586 586
587 587 if '|' in err_msg:
588 588 name, err = err_msg.split('|')
589 589 if 'SchainWarning' in err:
590 590 log.warning(err.split('SchainWarning:')[-1].split('\n')[0].strip(), name)
591 591 elif 'SchainError' in err:
592 592 log.error(err.split('SchainError:')[-1].split('\n')[0].strip(), name)
593 593 else:
594 594 log.error(err, name)
595 595 else:
596 596 name, err = self.name, err_msg
597 597
598 598 time.sleep(1)
599 599
600 600 ctx.term()
601 601
602 602 message = ''.join(err)
603 603
604 604 if err_msg:
605 605 subject = 'SChain v%s: Error running %s\n' % (
606 606 schainpy.__version__, self.name)
607 607
608 608 subtitle = 'Hostname: %s\n' % socket.gethostbyname(
609 609 socket.gethostname())
610 610 subtitle += 'Working directory: %s\n' % os.path.abspath('./')
611 611 subtitle += 'Configuration file: %s\n' % self.filename
612 612 subtitle += 'Time: %s\n' % str(datetime.datetime.now())
613 613
614 614 readUnitConfObj = self.getReadUnit()
615 615 if readUnitConfObj:
616 616 subtitle += '\nInput parameters:\n'
617 617 subtitle += '[Data path = %s]\n' % readUnitConfObj.parameters['path']
618 618 subtitle += '[Start date = %s]\n' % readUnitConfObj.parameters['startDate']
619 619 subtitle += '[End date = %s]\n' % readUnitConfObj.parameters['endDate']
620 620 subtitle += '[Start time = %s]\n' % readUnitConfObj.parameters['startTime']
621 621 subtitle += '[End time = %s]\n' % readUnitConfObj.parameters['endTime']
622 622
623 623 a = Alarm(
624 624 modes=self.alarm,
625 625 email=self.email,
626 626 message=message,
627 627 subject=subject,
628 628 subtitle=subtitle,
629 629 filename=self.filename
630 630 )
631 631
632 632 a.start()
633 633
634 634 def setFilename(self, filename):
635 635
636 636 self.filename = filename
637 637
638 638 def runProcs(self):
639 639
640 640 err = False
641 641 n = len(self.configurations)
642 642 #print(n)
643 643
644 644 while not err:
645 645 #print(self.getUnits())
646 646 for conf in self.getUnits():
647 647 #print(conf)
648 648 ok = conf.run()
649 649 #print("ok", ok)
650 650 if ok == 'Error':
651 651 n -= 1
652 652 continue
653 653 elif not ok:
654 654 break
655 655 #print("****************************************************end")
656 #exit(1)
656 657 if n == 0:
657 658 err = True
658 659
659 660 def run(self):
660 661
661 662 log.success('\nStarting Project {} [id={}]'.format(self.name, self.id), tag='')
662 663 self.started = True
663 664 self.start_time = time.time()
664 665 self.createObjects()
665 666 self.runProcs()
666 667 log.success('{} Done (Time: {:4.2f}s)'.format(
667 668 self.name,
668 669 time.time() - self.start_time), '')
Show file before | Show file after | Hide comments
@@ -1,1288 +1,1290
1 1 # Copyright (c) 2012-2021 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Classes to plot Spectra data
6 6
7 7 """
8 8
9 9 import os
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt, log
13 13
14 14 class SpectraPlot(Plot):
15 15 '''
16 16 Plot for Spectra data
17 17 '''
18 18
19 19 CODE = 'spc'
20 20 colormap = 'jet'
21 21 plot_type = 'pcolor'
22 22 buffering = False
23 23
24 24 def setup(self):
25 25
26 26 self.nplots = len(self.data.channels)
27 27 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
28 28 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
29 29 self.height = 2.6 * self.nrows
30 30 self.cb_label = 'dB'
31 31 if self.showprofile:
32 32 self.width = 4 * self.ncols
33 33 else:
34 34 self.width = 3.5 * self.ncols
35 35 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
36 36 self.ylabel = 'Range [km]'
37 37
38 38 def update(self, dataOut):
39 39
40 40 data = {}
41 41 meta = {}
42 42 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
43 43 data['spc'] = spc
44 44 data['rti'] = dataOut.getPower()
45 #print("NormFactor: ",dataOut.normFactor)
45 46 #data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
46 47 if hasattr(dataOut, 'LagPlot'): #Double Pulse
47 48 max_hei_id = dataOut.nHeights - 2*dataOut.LagPlot
48 49 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=46,ymax_index=max_hei_id)/dataOut.normFactor)
49 50 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=40,ymax_index=max_hei_id)/dataOut.normFactor)
50 51 data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=53,ymax_index=max_hei_id)/dataOut.normFactor)
51 52 data['noise'][0] = 10*numpy.log10(dataOut.getNoise(ymin_index=53)[0]/dataOut.normFactor)
52 53 #data['noise'][1] = 22.035507
53 54 else:
54 55 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
55 56 #data['noise'] = 10*numpy.log10(dataOut.getNoise(ymin_index=26,ymax_index=44)/dataOut.normFactor)
56 57 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
57 58
58 59 if self.CODE == 'spc_moments':
59 60 data['moments'] = dataOut.moments
60 61 if self.CODE == 'gaussian_fit':
61 62 data['gaussfit'] = dataOut.DGauFitParams
62 63
63 64 return data, meta
64 65
65 66 def plot(self):
66 67
67 68 if self.xaxis == "frequency":
68 69 x = self.data.xrange[0]
69 70 self.xlabel = "Frequency (kHz)"
70 71 elif self.xaxis == "time":
71 72 x = self.data.xrange[1]
72 73 self.xlabel = "Time (ms)"
73 74 else:
74 75 x = self.data.xrange[2]
75 76 self.xlabel = "Velocity (m/s)"
76 77
77 78 if (self.CODE == 'spc_moments') | (self.CODE == 'gaussian_fit'):
78 79 x = self.data.xrange[2]
79 80 self.xlabel = "Velocity (m/s)"
80 81
81 82 self.titles = []
82 83
83 84 y = self.data.yrange
84 85 self.y = y
85 86
86 87 data = self.data[-1]
87 88 z = data['spc']
88 89
89 90 self.CODE2 = 'spc_oblique'
90 91
91 92
92 93 for n, ax in enumerate(self.axes):
93 94 noise = data['noise'][n]
94 95 if self.CODE == 'spc_moments':
95 96 mean = data['moments'][n, 1]
96 97 if self.CODE == 'gaussian_fit':
97 98 gau0 = data['gaussfit'][n][2,:,0]
98 99 gau1 = data['gaussfit'][n][2,:,1]
99 100 if ax.firsttime:
100 101 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
101 102 self.xmin = self.xmin if self.xmin else numpy.nanmin(x)#-self.xmax
102 103 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
103 104 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
105
104 106 ax.plt = ax.pcolormesh(x, y, z[n].T,
105 107 vmin=self.zmin,
106 108 vmax=self.zmax,
107 cmap=plt.get_cmap(self.colormap)
109 cmap=plt.get_cmap(self.colormap),
108 110 )
109 111
110 112 if self.showprofile:
111 113 ax.plt_profile = self.pf_axes[n].plot(
112 114 data['rti'][n], y)[0]
113 115 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
114 116 color="k", linestyle="dashed", lw=1)[0]
115 117 if self.CODE == 'spc_moments':
116 118 ax.plt_mean = ax.plot(mean, y, color='k', lw=1)[0]
117 119 if self.CODE == 'gaussian_fit':
118 120 ax.plt_gau0 = ax.plot(gau0, y, color='r', lw=1)[0]
119 121 ax.plt_gau1 = ax.plot(gau1, y, color='y', lw=1)[0]
120 122 else:
121 123 ax.plt.set_array(z[n].T.ravel())
122 124 if self.showprofile:
123 125 ax.plt_profile.set_data(data['rti'][n], y)
124 126 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
125 127 if self.CODE == 'spc_moments':
126 128 ax.plt_mean.set_data(mean, y)
127 129 if self.CODE == 'gaussian_fit':
128 130 ax.plt_gau0.set_data(gau0, y)
129 131 ax.plt_gau1.set_data(gau1, y)
130 132 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
131 133
132 134 class SpectraObliquePlot(Plot):
133 135 '''
134 136 Plot for Spectra data
135 137 '''
136 138
137 139 CODE = 'spc_oblique'
138 140 colormap = 'jet'
139 141 plot_type = 'pcolor'
140 142
141 143 def setup(self):
142 144 self.xaxis = "oblique"
143 145 self.nplots = len(self.data.channels)
144 146 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
145 147 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
146 148 self.height = 2.6 * self.nrows
147 149 self.cb_label = 'dB'
148 150 if self.showprofile:
149 151 self.width = 4 * self.ncols
150 152 else:
151 153 self.width = 3.5 * self.ncols
152 154 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
153 155 self.ylabel = 'Range [km]'
154 156
155 157 def update(self, dataOut):
156 158
157 159 data = {}
158 160 meta = {}
159 161 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
160 162 data['spc'] = spc
161 163 data['rti'] = dataOut.getPower()
162 164 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
163 165 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
164 166
165 167 data['shift1'] = dataOut.Oblique_params[0][1]
166 168 data['shift2'] = dataOut.Oblique_params[0][4]
167 169 data['shift1_error'] = dataOut.Oblique_param_errors[0][1]
168 170 data['shift2_error'] = dataOut.Oblique_param_errors[0][4]
169 171
170 172 return data, meta
171 173
172 174 def plot(self):
173 175
174 176 if self.xaxis == "frequency":
175 177 x = self.data.xrange[0]
176 178 self.xlabel = "Frequency (kHz)"
177 179 elif self.xaxis == "time":
178 180 x = self.data.xrange[1]
179 181 self.xlabel = "Time (ms)"
180 182 else:
181 183 x = self.data.xrange[2]
182 184 self.xlabel = "Velocity (m/s)"
183 185
184 186 self.titles = []
185 187
186 188 y = self.data.yrange
187 189 self.y = y
188 190 z = self.data['spc']
189 191
190 192 for n, ax in enumerate(self.axes):
191 193 noise = self.data['noise'][n][-1]
192 194 shift1 = self.data['shift1']
193 195 shift2 = self.data['shift2']
194 196 err1 = self.data['shift1_error']
195 197 err2 = self.data['shift2_error']
196 198 if ax.firsttime:
197 199 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
198 200 self.xmin = self.xmin if self.xmin else -self.xmax
199 201 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
200 202 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
201 203 ax.plt = ax.pcolormesh(x, y, z[n].T,
202 204 vmin=self.zmin,
203 205 vmax=self.zmax,
204 206 cmap=plt.get_cmap(self.colormap)
205 207 )
206 208
207 209 if self.showprofile:
208 210 ax.plt_profile = self.pf_axes[n].plot(
209 211 self.data['rti'][n][-1], y)[0]
210 212 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
211 213 color="k", linestyle="dashed", lw=1)[0]
212 214
213 215 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^', elinewidth=0.2, marker='x', linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
214 216 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=0.2,marker='x',linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
215 217 else:
216 218 self.ploterr1.remove()
217 219 self.ploterr2.remove()
218 220 ax.plt.set_array(z[n].T.ravel())
219 221 if self.showprofile:
220 222 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
221 223 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
222 224 self.ploterr1 = ax.errorbar(shift1, y, xerr=err1, fmt='k^',elinewidth=0.2,marker='x',linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
223 225 self.ploterr2 = ax.errorbar(shift2, y, xerr=err2, fmt='m^',elinewidth=0.2,marker='x',linestyle='None',markersize=0.5,capsize=0.3,markeredgewidth=0.2)
224 226
225 227 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
226 228
227 229
228 230 class CrossSpectraPlot(Plot):
229 231
230 232 CODE = 'cspc'
231 233 colormap = 'jet'
232 234 plot_type = 'pcolor'
233 235 zmin_coh = None
234 236 zmax_coh = None
235 237 zmin_phase = None
236 238 zmax_phase = None
237 239
238 240 def setup(self):
239 241
240 242 self.ncols = 4
241 243 self.nplots = len(self.data.pairs) * 2
242 244 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
243 245 self.width = 3.1 * self.ncols
244 246 self.height = 5 * self.nrows
245 247 self.ylabel = 'Range [km]'
246 248 self.showprofile = False
247 249 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
248 250
249 251 def update(self, dataOut):
250 252
251 253 data = {}
252 254 meta = {}
253 255
254 256 spc = dataOut.data_spc
255 257 cspc = dataOut.data_cspc
256 258 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
257 259 meta['pairs'] = dataOut.pairsList
258 260
259 261 tmp = []
260 262
261 263 for n, pair in enumerate(meta['pairs']):
262 264 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
263 265 coh = numpy.abs(out)
264 266 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
265 267 tmp.append(coh)
266 268 tmp.append(phase)
267 269
268 270 data['cspc'] = numpy.array(tmp)
269 271
270 272 return data, meta
271 273
272 274 def plot(self):
273 275
274 276 if self.xaxis == "frequency":
275 277 x = self.data.xrange[0]
276 278 self.xlabel = "Frequency (kHz)"
277 279 elif self.xaxis == "time":
278 280 x = self.data.xrange[1]
279 281 self.xlabel = "Time (ms)"
280 282 else:
281 283 x = self.data.xrange[2]
282 284 self.xlabel = "Velocity (m/s)"
283 285
284 286 self.titles = []
285 287
286 288 y = self.data.yrange
287 289 self.y = y
288 290
289 291 data = self.data[-1]
290 292 cspc = data['cspc']
291 293
292 294 for n in range(len(self.data.pairs)):
293 295 pair = self.data.pairs[n]
294 296 coh = cspc[n*2]
295 297 phase = cspc[n*2+1]
296 298 ax = self.axes[2 * n]
297 299 if ax.firsttime:
298 300 ax.plt = ax.pcolormesh(x, y, coh.T,
299 301 vmin=0,
300 302 vmax=1,
301 303 cmap=plt.get_cmap(self.colormap_coh)
302 304 )
303 305 else:
304 306 ax.plt.set_array(coh.T.ravel())
305 307 self.titles.append(
306 308 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
307 309
308 310 ax = self.axes[2 * n + 1]
309 311 if ax.firsttime:
310 312 ax.plt = ax.pcolormesh(x, y, phase.T,
311 313 vmin=-180,
312 314 vmax=180,
313 315 cmap=plt.get_cmap(self.colormap_phase)
314 316 )
315 317 else:
316 318 ax.plt.set_array(phase.T.ravel())
317 319 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
318 320
319 321
320 322 class CrossSpectra4Plot(Plot):
321 323
322 324 CODE = 'cspc'
323 325 colormap = 'jet'
324 326 plot_type = 'pcolor'
325 327 zmin_coh = None
326 328 zmax_coh = None
327 329 zmin_phase = None
328 330 zmax_phase = None
329 331
330 332 def setup(self):
331 333
332 334 self.ncols = 4
333 335 self.nrows = len(self.data.pairs)
334 336 self.nplots = self.nrows * 4
335 337 self.width = 3.1 * self.ncols
336 338 self.height = 5 * self.nrows
337 339 self.ylabel = 'Range [km]'
338 340 self.showprofile = False
339 341 self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
340 342
341 343 def plot(self):
342 344
343 345 if self.xaxis == "frequency":
344 346 x = self.data.xrange[0]
345 347 self.xlabel = "Frequency (kHz)"
346 348 elif self.xaxis == "time":
347 349 x = self.data.xrange[1]
348 350 self.xlabel = "Time (ms)"
349 351 else:
350 352 x = self.data.xrange[2]
351 353 self.xlabel = "Velocity (m/s)"
352 354
353 355 self.titles = []
354 356
355 357
356 358 y = self.data.heights
357 359 self.y = y
358 360 nspc = self.data['spc']
359 361 #print(numpy.shape(self.data['spc']))
360 362 spc = self.data['cspc'][0]
361 363 #print(numpy.shape(nspc))
362 364 #exit()
363 365 #nspc[1,:,:] = numpy.flip(nspc[1,:,:],axis=0)
364 366 #print(numpy.shape(spc))
365 367 #exit()
366 368 cspc = self.data['cspc'][1]
367 369
368 370 #xflip=numpy.flip(x)
369 371 #print(numpy.shape(cspc))
370 372 #exit()
371 373
372 374 for n in range(self.nrows):
373 375 noise = self.data['noise'][:,-1]
374 376 pair = self.data.pairs[n]
375 377 #print(pair)
376 378 #exit()
377 379 ax = self.axes[4 * n]
378 380 if ax.firsttime:
379 381 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
380 382 self.xmin = self.xmin if self.xmin else -self.xmax
381 383 self.zmin = self.zmin if self.zmin else numpy.nanmin(nspc)
382 384 self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
383 385 ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
384 386 vmin=self.zmin,
385 387 vmax=self.zmax,
386 388 cmap=plt.get_cmap(self.colormap)
387 389 )
388 390 else:
389 391 #print(numpy.shape(nspc[pair[0]].T))
390 392 #exit()
391 393 ax.plt.set_array(nspc[pair[0]].T.ravel())
392 394 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
393 395
394 396 ax = self.axes[4 * n + 1]
395 397
396 398 if ax.firsttime:
397 399 ax.plt = ax.pcolormesh(x , y, numpy.flip(nspc[pair[1]],axis=0).T,
398 400 vmin=self.zmin,
399 401 vmax=self.zmax,
400 402 cmap=plt.get_cmap(self.colormap)
401 403 )
402 404 else:
403 405
404 406 ax.plt.set_array(numpy.flip(nspc[pair[1]],axis=0).T.ravel())
405 407 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
406 408
407 409 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
408 410 coh = numpy.abs(out)
409 411 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
410 412
411 413 ax = self.axes[4 * n + 2]
412 414 if ax.firsttime:
413 415 ax.plt = ax.pcolormesh(x, y, numpy.flip(coh,axis=0).T,
414 416 vmin=0,
415 417 vmax=1,
416 418 cmap=plt.get_cmap(self.colormap_coh)
417 419 )
418 420 else:
419 421 ax.plt.set_array(numpy.flip(coh,axis=0).T.ravel())
420 422 self.titles.append(
421 423 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
422 424
423 425 ax = self.axes[4 * n + 3]
424 426 if ax.firsttime:
425 427 ax.plt = ax.pcolormesh(x, y, numpy.flip(phase,axis=0).T,
426 428 vmin=-180,
427 429 vmax=180,
428 430 cmap=plt.get_cmap(self.colormap_phase)
429 431 )
430 432 else:
431 433 ax.plt.set_array(numpy.flip(phase,axis=0).T.ravel())
432 434 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
433 435
434 436
435 437 class CrossSpectra2Plot(Plot):
436 438
437 439 CODE = 'cspc'
438 440 colormap = 'jet'
439 441 plot_type = 'pcolor'
440 442 zmin_coh = None
441 443 zmax_coh = None
442 444 zmin_phase = None
443 445 zmax_phase = None
444 446
445 447 def setup(self):
446 448
447 449 self.ncols = 1
448 450 self.nrows = len(self.data.pairs)
449 451 self.nplots = self.nrows * 1
450 452 self.width = 3.1 * self.ncols
451 453 self.height = 5 * self.nrows
452 454 self.ylabel = 'Range [km]'
453 455 self.showprofile = False
454 456 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
455 457
456 458 def plot(self):
457 459
458 460 if self.xaxis == "frequency":
459 461 x = self.data.xrange[0]
460 462 self.xlabel = "Frequency (kHz)"
461 463 elif self.xaxis == "time":
462 464 x = self.data.xrange[1]
463 465 self.xlabel = "Time (ms)"
464 466 else:
465 467 x = self.data.xrange[2]
466 468 self.xlabel = "Velocity (m/s)"
467 469
468 470 self.titles = []
469 471
470 472
471 473 y = self.data.heights
472 474 self.y = y
473 475 #nspc = self.data['spc']
474 476 #print(numpy.shape(self.data['spc']))
475 477 #spc = self.data['cspc'][0]
476 478 #print(numpy.shape(spc))
477 479 #exit()
478 480 cspc = self.data['cspc'][1]
479 481 #print(numpy.shape(cspc))
480 482 #exit()
481 483
482 484 for n in range(self.nrows):
483 485 noise = self.data['noise'][:,-1]
484 486 pair = self.data.pairs[n]
485 487 #print(pair) #exit()
486 488
487 489
488 490
489 491 out = cspc[n]# / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
490 492
491 493 #print(out[:,53])
492 494 #exit()
493 495 cross = numpy.abs(out)
494 496 z = cross/self.data.nFactor
495 497 #print("here")
496 498 #print(dataOut.data_spc[0,0,0])
497 499 #exit()
498 500
499 501 cross = 10*numpy.log10(z)
500 502 #print(numpy.shape(cross))
501 503 #print(cross[0,:])
502 504 #print(self.data.nFactor)
503 505 #exit()
504 506 #phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
505 507
506 508 ax = self.axes[1 * n]
507 509 if ax.firsttime:
508 510 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
509 511 self.xmin = self.xmin if self.xmin else -self.xmax
510 512 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
511 513 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
512 514 ax.plt = ax.pcolormesh(x, y, cross.T,
513 515 vmin=self.zmin,
514 516 vmax=self.zmax,
515 517 cmap=plt.get_cmap(self.colormap)
516 518 )
517 519 else:
518 520 ax.plt.set_array(cross.T.ravel())
519 521 self.titles.append(
520 522 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
521 523
522 524
523 525 class CrossSpectra3Plot(Plot):
524 526
525 527 CODE = 'cspc'
526 528 colormap = 'jet'
527 529 plot_type = 'pcolor'
528 530 zmin_coh = None
529 531 zmax_coh = None
530 532 zmin_phase = None
531 533 zmax_phase = None
532 534
533 535 def setup(self):
534 536
535 537 self.ncols = 3
536 538 self.nrows = len(self.data.pairs)
537 539 self.nplots = self.nrows * 3
538 540 self.width = 3.1 * self.ncols
539 541 self.height = 5 * self.nrows
540 542 self.ylabel = 'Range [km]'
541 543 self.showprofile = False
542 544 self.plots_adjust.update({'left': 0.22, 'right': .90, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
543 545
544 546 def plot(self):
545 547
546 548 if self.xaxis == "frequency":
547 549 x = self.data.xrange[0]
548 550 self.xlabel = "Frequency (kHz)"
549 551 elif self.xaxis == "time":
550 552 x = self.data.xrange[1]
551 553 self.xlabel = "Time (ms)"
552 554 else:
553 555 x = self.data.xrange[2]
554 556 self.xlabel = "Velocity (m/s)"
555 557
556 558 self.titles = []
557 559
558 560
559 561 y = self.data.heights
560 562 self.y = y
561 563 #nspc = self.data['spc']
562 564 #print(numpy.shape(self.data['spc']))
563 565 #spc = self.data['cspc'][0]
564 566 #print(numpy.shape(spc))
565 567 #exit()
566 568 cspc = self.data['cspc'][1]
567 569 #print(numpy.shape(cspc))
568 570 #exit()
569 571
570 572 for n in range(self.nrows):
571 573 noise = self.data['noise'][:,-1]
572 574 pair = self.data.pairs[n]
573 575 #print(pair) #exit()
574 576
575 577
576 578
577 579 out = cspc[n]# / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
578 580
579 581 #print(out[:,53])
580 582 #exit()
581 583 cross = numpy.abs(out)
582 584 z = cross/self.data.nFactor
583 585 cross = 10*numpy.log10(z)
584 586
585 587 out_r= out.real/self.data.nFactor
586 588 #out_r = 10*numpy.log10(out_r)
587 589
588 590 out_i= out.imag/self.data.nFactor
589 591 #out_i = 10*numpy.log10(out_i)
590 592 #print(numpy.shape(cross))
591 593 #print(cross[0,:])
592 594 #print(self.data.nFactor)
593 595 #exit()
594 596 #phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
595 597
596 598 ax = self.axes[3 * n]
597 599 if ax.firsttime:
598 600 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
599 601 self.xmin = self.xmin if self.xmin else -self.xmax
600 602 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
601 603 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
602 604 ax.plt = ax.pcolormesh(x, y, cross.T,
603 605 vmin=self.zmin,
604 606 vmax=self.zmax,
605 607 cmap=plt.get_cmap(self.colormap)
606 608 )
607 609 else:
608 610 ax.plt.set_array(cross.T.ravel())
609 611 self.titles.append(
610 612 'Cross Spectra Power Ch{} * Ch{}'.format(pair[0], pair[1]))
611 613
612 614 ax = self.axes[3 * n + 1]
613 615 if ax.firsttime:
614 616 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
615 617 self.xmin = self.xmin if self.xmin else -self.xmax
616 618 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
617 619 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
618 620 ax.plt = ax.pcolormesh(x, y, out_r.T,
619 621 vmin=-1.e6,
620 622 vmax=0,
621 623 cmap=plt.get_cmap(self.colormap)
622 624 )
623 625 else:
624 626 ax.plt.set_array(out_r.T.ravel())
625 627 self.titles.append(
626 628 'Cross Spectra Real Ch{} * Ch{}'.format(pair[0], pair[1]))
627 629
628 630 ax = self.axes[3 * n + 2]
629 631
630 632
631 633 if ax.firsttime:
632 634 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
633 635 self.xmin = self.xmin if self.xmin else -self.xmax
634 636 self.zmin = self.zmin if self.zmin else numpy.nanmin(cross)
635 637 self.zmax = self.zmax if self.zmax else numpy.nanmax(cross)
636 638 ax.plt = ax.pcolormesh(x, y, out_i.T,
637 639 vmin=-1.e6,
638 640 vmax=1.e6,
639 641 cmap=plt.get_cmap(self.colormap)
640 642 )
641 643 else:
642 644 ax.plt.set_array(out_i.T.ravel())
643 645 self.titles.append(
644 646 'Cross Spectra Imag Ch{} * Ch{}'.format(pair[0], pair[1]))
645 647
646 648 class RTIPlot(Plot):
647 649 '''
648 650 Plot for RTI data
649 651 '''
650 652
651 653 CODE = 'rti'
652 654 colormap = 'jet'
653 655 plot_type = 'pcolorbuffer'
654 656
655 657 def setup(self):
656 658 self.xaxis = 'time'
657 659 self.ncols = 1
658 660 self.nrows = len(self.data.channels)
659 661 self.nplots = len(self.data.channels)
660 662 self.ylabel = 'Range [km]'
661 663 self.xlabel = 'Time'
662 664 self.cb_label = 'dB'
663 665 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
664 666 self.titles = ['{} Channel {}'.format(
665 667 self.CODE.upper(), x) for x in range(self.nrows)]
666 668
667 669 def update(self, dataOut):
668 670
669 671 data = {}
670 672 meta = {}
671 673 data['rti'] = dataOut.getPower()
672 674
673 675 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
674 676
675 677 return data, meta
676 678
677 679 def plot(self):
678 680
679 681 self.x = self.data.times
680 682 self.y = self.data.yrange
681 683 self.z = self.data[self.CODE]
682 684
683 685 self.z = numpy.ma.masked_invalid(self.z)
684 686
685 687 if self.decimation is None:
686 688 x, y, z = self.fill_gaps(self.x, self.y, self.z)
687 689 else:
688 690 x, y, z = self.fill_gaps(*self.decimate())
689 691
690 692 for n, ax in enumerate(self.axes):
691 693 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
692 694 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
693 695 if ax.firsttime:
694 696 ax.plt = ax.pcolormesh(x, y, z[n].T,
695 697 vmin=self.zmin,
696 698 vmax=self.zmax,
697 699 cmap=plt.get_cmap(self.colormap)
698 700 )
699 701 if self.showprofile:
700 702 ax.plot_profile = self.pf_axes[n].plot(
701 703 self.data['rti'][n][-1], self.y)[0]
702 704 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
703 705 color="k", linestyle="dashed", lw=1)[0]
704 706 else:
705 707 ax.collections.remove(ax.collections[0])
706 708 ax.plt = ax.pcolormesh(x, y, z[n].T,
707 709 vmin=self.zmin,
708 710 vmax=self.zmax,
709 711 cmap=plt.get_cmap(self.colormap)
710 712 )
711 713 if self.showprofile:
712 714 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
713 715 ax.plot_noise.set_data(numpy.repeat(
714 716 self.data['noise'][n][-1], len(self.y)), self.y)
715 717
716 718
717 719 class SpectrogramPlot(Plot):
718 720 '''
719 721 Plot for Spectrogram data
720 722 '''
721 723
722 724 CODE = 'Spectrogram_Profile'
723 725 colormap = 'binary'
724 726 plot_type = 'pcolorbuffer'
725 727
726 728 def setup(self):
727 729 self.xaxis = 'time'
728 730 self.ncols = 1
729 731 self.nrows = len(self.data.channels)
730 732 self.nplots = len(self.data.channels)
731 733 self.xlabel = 'Time'
732 734 #self.cb_label = 'dB'
733 735 self.plots_adjust.update({'hspace':1.2, 'left': 0.1, 'bottom': 0.12, 'right':0.95})
734 736 self.titles = []
735 737
736 738 #self.titles = ['{} Channel {} \n H = {} km ({} - {})'.format(
737 739 #self.CODE.upper(), x, self.data.heightList[self.data.hei], self.data.heightList[self.data.hei],self.data.heightList[self.data.hei]+(self.data.DH*self.data.nProfiles)) for x in range(self.nrows)]
738 740
739 741 self.titles = ['{} Channel {}'.format(
740 742 self.CODE.upper(), x) for x in range(self.nrows)]
741 743
742 744
743 745 def update(self, dataOut):
744 746 data = {}
745 747 meta = {}
746 748
747 749 maxHei = 1620#+12000
748 750 indb = numpy.where(dataOut.heightList <= maxHei)
749 751 hei = indb[0][-1]
750 752 #print(dataOut.heightList)
751 753
752 754 factor = dataOut.nIncohInt
753 755 z = dataOut.data_spc[:,:,hei] / factor
754 756 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
755 757 #buffer = 10 * numpy.log10(z)
756 758
757 759 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
758 760
759 761
760 762 #self.hei = hei
761 763 #self.heightList = dataOut.heightList
762 764 #self.DH = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
763 765 #self.nProfiles = dataOut.nProfiles
764 766
765 767 data['Spectrogram_Profile'] = 10 * numpy.log10(z)
766 768
767 769 data['hei'] = hei
768 770 data['DH'] = (dataOut.heightList[1] - dataOut.heightList[0])/dataOut.step
769 771 data['nProfiles'] = dataOut.nProfiles
770 772 #meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
771 773 '''
772 774 import matplotlib.pyplot as plt
773 775 plt.plot(10 * numpy.log10(z[0,:]))
774 776 plt.show()
775 777
776 778 from time import sleep
777 779 sleep(10)
778 780 '''
779 781 return data, meta
780 782
781 783 def plot(self):
782 784
783 785 self.x = self.data.times
784 786 self.z = self.data[self.CODE]
785 787 self.y = self.data.xrange[0]
786 788
787 789 hei = self.data['hei'][-1]
788 790 DH = self.data['DH'][-1]
789 791 nProfiles = self.data['nProfiles'][-1]
790 792
791 793 self.ylabel = "Frequency (kHz)"
792 794
793 795 self.z = numpy.ma.masked_invalid(self.z)
794 796
795 797 if self.decimation is None:
796 798 x, y, z = self.fill_gaps(self.x, self.y, self.z)
797 799 else:
798 800 x, y, z = self.fill_gaps(*self.decimate())
799 801
800 802 for n, ax in enumerate(self.axes):
801 803 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
802 804 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
803 805 data = self.data[-1]
804 806 if ax.firsttime:
805 807 ax.plt = ax.pcolormesh(x, y, z[n].T,
806 808 vmin=self.zmin,
807 809 vmax=self.zmax,
808 810 cmap=plt.get_cmap(self.colormap)
809 811 )
810 812 else:
811 813 ax.collections.remove(ax.collections[0])
812 814 ax.plt = ax.pcolormesh(x, y, z[n].T,
813 815 vmin=self.zmin,
814 816 vmax=self.zmax,
815 817 cmap=plt.get_cmap(self.colormap)
816 818 )
817 819
818 820 #self.titles.append('Spectrogram')
819 821
820 822 #self.titles.append('{} Channel {} \n H = {} km ({} - {})'.format(
821 823 #self.CODE.upper(), x, y[hei], y[hei],y[hei]+(DH*nProfiles)))
822 824
823 825
824 826
825 827
826 828 class CoherencePlot(RTIPlot):
827 829 '''
828 830 Plot for Coherence data
829 831 '''
830 832
831 833 CODE = 'coh'
832 834
833 835 def setup(self):
834 836 self.xaxis = 'time'
835 837 self.ncols = 1
836 838 self.nrows = len(self.data.pairs)
837 839 self.nplots = len(self.data.pairs)
838 840 self.ylabel = 'Range [km]'
839 841 self.xlabel = 'Time'
840 842 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
841 843 if self.CODE == 'coh':
842 844 self.cb_label = ''
843 845 self.titles = [
844 846 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
845 847 else:
846 848 self.cb_label = 'Degrees'
847 849 self.titles = [
848 850 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
849 851
850 852 def update(self, dataOut):
851 853
852 854 data = {}
853 855 meta = {}
854 856 data['coh'] = dataOut.getCoherence()
855 857 meta['pairs'] = dataOut.pairsList
856 858
857 859 return data, meta
858 860
859 861 class PhasePlot(CoherencePlot):
860 862 '''
861 863 Plot for Phase map data
862 864 '''
863 865
864 866 CODE = 'phase'
865 867 colormap = 'seismic'
866 868
867 869 def update(self, dataOut):
868 870
869 871 data = {}
870 872 meta = {}
871 873 data['phase'] = dataOut.getCoherence(phase=True)
872 874 meta['pairs'] = dataOut.pairsList
873 875
874 876 return data, meta
875 877
876 878 class NoisePlot(Plot):
877 879 '''
878 880 Plot for noise
879 881 '''
880 882
881 883 CODE = 'noise'
882 884 plot_type = 'scatterbuffer'
883 885
884 886 def setup(self):
885 887 self.xaxis = 'time'
886 888 self.ncols = 1
887 889 self.nrows = 1
888 890 self.nplots = 1
889 891 self.ylabel = 'Intensity [dB]'
890 892 self.xlabel = 'Time'
891 893 self.titles = ['Noise']
892 894 self.colorbar = False
893 895 self.plots_adjust.update({'right': 0.85 })
894 896
895 897 def update(self, dataOut):
896 898
897 899 data = {}
898 900 meta = {}
899 901 data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
900 902 meta['yrange'] = numpy.array([])
901 903
902 904 return data, meta
903 905
904 906 def plot(self):
905 907
906 908 x = self.data.times
907 909 xmin = self.data.min_time
908 910 xmax = xmin + self.xrange * 60 * 60
909 911 Y = self.data['noise']
910 912
911 913 if self.axes[0].firsttime:
912 914 self.ymin = numpy.nanmin(Y) - 5
913 915 self.ymax = numpy.nanmax(Y) + 5
914 916 for ch in self.data.channels:
915 917 y = Y[ch]
916 918 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
917 919 plt.legend(bbox_to_anchor=(1.18, 1.0))
918 920 else:
919 921 for ch in self.data.channels:
920 922 y = Y[ch]
921 923 self.axes[0].lines[ch].set_data(x, y)
922 924
923 925 self.ymin = numpy.nanmin(Y) - 5
924 926 self.ymax = numpy.nanmax(Y) + 10
925 927
926 928
927 929 class PowerProfilePlot(Plot):
928 930
929 931 CODE = 'pow_profile'
930 932 plot_type = 'scatter'
931 933
932 934 def setup(self):
933 935
934 936 self.ncols = 1
935 937 self.nrows = 1
936 938 self.nplots = 1
937 939 self.height = 4
938 940 self.width = 3
939 941 self.ylabel = 'Range [km]'
940 942 self.xlabel = 'Intensity [dB]'
941 943 self.titles = ['Power Profile']
942 944 self.colorbar = False
943 945
944 946 def update(self, dataOut):
945 947
946 948 data = {}
947 949 meta = {}
948 950 data[self.CODE] = dataOut.getPower()
949 951
950 952 return data, meta
951 953
952 954 def plot(self):
953 955
954 956 y = self.data.yrange
955 957 self.y = y
956 958
957 959 x = self.data[-1][self.CODE]
958 960
959 961 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
960 962 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
961 963
962 964 if self.axes[0].firsttime:
963 965 for ch in self.data.channels:
964 966 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
965 967 plt.legend()
966 968 else:
967 969 for ch in self.data.channels:
968 970 self.axes[0].lines[ch].set_data(x[ch], y)
969 971
970 972
971 973 class SpectraCutPlot(Plot):
972 974
973 975 CODE = 'spc_cut'
974 976 plot_type = 'scatter'
975 977 buffering = False
976 978
977 979 def setup(self):
978 980
979 981 self.nplots = len(self.data.channels)
980 982 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
981 983 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
982 984 self.width = 3.4 * self.ncols + 1.5
983 985 self.height = 3 * self.nrows
984 986 self.ylabel = 'Power [dB]'
985 987 self.colorbar = False
986 988 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
987 989
988 990 def update(self, dataOut):
989 991
990 992 data = {}
991 993 meta = {}
992 994 spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
993 995 data['spc'] = spc
994 996 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
995 997 if self.CODE == 'cut_gaussian_fit':
996 998 data['gauss_fit0'] = 10*numpy.log10(dataOut.GaussFit0/dataOut.normFactor)
997 999 data['gauss_fit1'] = 10*numpy.log10(dataOut.GaussFit1/dataOut.normFactor)
998 1000 return data, meta
999 1001
1000 1002 def plot(self):
1001 1003 if self.xaxis == "frequency":
1002 1004 x = self.data.xrange[0][1:]
1003 1005 self.xlabel = "Frequency (kHz)"
1004 1006 elif self.xaxis == "time":
1005 1007 x = self.data.xrange[1]
1006 1008 self.xlabel = "Time (ms)"
1007 1009 else:
1008 1010 x = self.data.xrange[2][:-1]
1009 1011 self.xlabel = "Velocity (m/s)"
1010 1012
1011 1013 if self.CODE == 'cut_gaussian_fit':
1012 1014 x = self.data.xrange[2][:-1]
1013 1015 self.xlabel = "Velocity (m/s)"
1014 1016
1015 1017 self.titles = []
1016 1018
1017 1019 y = self.data.yrange
1018 1020 data = self.data[-1]
1019 1021 z = data['spc']
1020 1022
1021 1023 if self.height_index:
1022 1024 index = numpy.array(self.height_index)
1023 1025 else:
1024 1026 index = numpy.arange(0, len(y), int((len(y))/9))
1025 1027
1026 1028 for n, ax in enumerate(self.axes):
1027 1029 if self.CODE == 'cut_gaussian_fit':
1028 1030 gau0 = data['gauss_fit0']
1029 1031 gau1 = data['gauss_fit1']
1030 1032 if ax.firsttime:
1031 1033 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
1032 1034 self.xmin = self.xmin if self.xmin else -self.xmax
1033 1035 self.ymin = self.ymin if self.ymin else numpy.nanmin(z[:,:,index])
1034 1036 self.ymax = self.ymax if self.ymax else numpy.nanmax(z[:,:,index])
1035 1037 #print(self.ymax)
1036 1038 #print(z[n, :, index])
1037 1039 ax.plt = ax.plot(x, z[n, :, index].T, lw=0.25)
1038 1040 if self.CODE == 'cut_gaussian_fit':
1039 1041 ax.plt_gau0 = ax.plot(x, gau0[n, :, index].T, lw=1, linestyle='-.')
1040 1042 for i, line in enumerate(ax.plt_gau0):
1041 1043 line.set_color(ax.plt[i].get_color())
1042 1044 ax.plt_gau1 = ax.plot(x, gau1[n, :, index].T, lw=1, linestyle='--')
1043 1045 for i, line in enumerate(ax.plt_gau1):
1044 1046 line.set_color(ax.plt[i].get_color())
1045 1047 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
1046 1048 self.figures[0].legend(ax.plt, labels, loc='center right')
1047 1049 else:
1048 1050 for i, line in enumerate(ax.plt):
1049 1051 line.set_data(x, z[n, :, index[i]].T)
1050 1052 for i, line in enumerate(ax.plt_gau0):
1051 1053 line.set_data(x, gau0[n, :, index[i]].T)
1052 1054 line.set_color(ax.plt[i].get_color())
1053 1055 for i, line in enumerate(ax.plt_gau1):
1054 1056 line.set_data(x, gau1[n, :, index[i]].T)
1055 1057 line.set_color(ax.plt[i].get_color())
1056 1058 self.titles.append('CH {}'.format(n))
1057 1059
1058 1060
1059 1061 class BeaconPhase(Plot):
1060 1062
1061 1063 __isConfig = None
1062 1064 __nsubplots = None
1063 1065
1064 1066 PREFIX = 'beacon_phase'
1065 1067
1066 1068 def __init__(self):
1067 1069 Plot.__init__(self)
1068 1070 self.timerange = 24*60*60
1069 1071 self.isConfig = False
1070 1072 self.__nsubplots = 1
1071 1073 self.counter_imagwr = 0
1072 1074 self.WIDTH = 800
1073 1075 self.HEIGHT = 400
1074 1076 self.WIDTHPROF = 120
1075 1077 self.HEIGHTPROF = 0
1076 1078 self.xdata = None
1077 1079 self.ydata = None
1078 1080
1079 1081 self.PLOT_CODE = BEACON_CODE
1080 1082
1081 1083 self.FTP_WEI = None
1082 1084 self.EXP_CODE = None
1083 1085 self.SUB_EXP_CODE = None
1084 1086 self.PLOT_POS = None
1085 1087
1086 1088 self.filename_phase = None
1087 1089
1088 1090 self.figfile = None
1089 1091
1090 1092 self.xmin = None
1091 1093 self.xmax = None
1092 1094
1093 1095 def getSubplots(self):
1094 1096
1095 1097 ncol = 1
1096 1098 nrow = 1
1097 1099
1098 1100 return nrow, ncol
1099 1101
1100 1102 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
1101 1103
1102 1104 self.__showprofile = showprofile
1103 1105 self.nplots = nplots
1104 1106
1105 1107 ncolspan = 7
1106 1108 colspan = 6
1107 1109 self.__nsubplots = 2
1108 1110
1109 1111 self.createFigure(id = id,
1110 1112 wintitle = wintitle,
1111 1113 widthplot = self.WIDTH+self.WIDTHPROF,
1112 1114 heightplot = self.HEIGHT+self.HEIGHTPROF,
1113 1115 show=show)
1114 1116
1115 1117 nrow, ncol = self.getSubplots()
1116 1118
1117 1119 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
1118 1120
1119 1121 def save_phase(self, filename_phase):
1120 1122 f = open(filename_phase,'w+')
1121 1123 f.write('\n\n')
1122 1124 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
1123 1125 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
1124 1126 f.close()
1125 1127
1126 1128 def save_data(self, filename_phase, data, data_datetime):
1127 1129 f=open(filename_phase,'a')
1128 1130 timetuple_data = data_datetime.timetuple()
1129 1131 day = str(timetuple_data.tm_mday)
1130 1132 month = str(timetuple_data.tm_mon)
1131 1133 year = str(timetuple_data.tm_year)
1132 1134 hour = str(timetuple_data.tm_hour)
1133 1135 minute = str(timetuple_data.tm_min)
1134 1136 second = str(timetuple_data.tm_sec)
1135 1137 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
1136 1138 f.close()
1137 1139
1138 1140 def plot(self):
1139 1141 log.warning('TODO: Not yet implemented...')
1140 1142
1141 1143 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
1142 1144 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
1143 1145 timerange=None,
1144 1146 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
1145 1147 server=None, folder=None, username=None, password=None,
1146 1148 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
1147 1149
1148 1150 if dataOut.flagNoData:
1149 1151 return dataOut
1150 1152
1151 1153 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
1152 1154 return
1153 1155
1154 1156 if pairsList == None:
1155 1157 pairsIndexList = dataOut.pairsIndexList[:10]
1156 1158 else:
1157 1159 pairsIndexList = []
1158 1160 for pair in pairsList:
1159 1161 if pair not in dataOut.pairsList:
1160 1162 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
1161 1163 pairsIndexList.append(dataOut.pairsList.index(pair))
1162 1164
1163 1165 if pairsIndexList == []:
1164 1166 return
1165 1167
1166 1168 # if len(pairsIndexList) > 4:
1167 1169 # pairsIndexList = pairsIndexList[0:4]
1168 1170
1169 1171 hmin_index = None
1170 1172 hmax_index = None
1171 1173
1172 1174 if hmin != None and hmax != None:
1173 1175 indexes = numpy.arange(dataOut.nHeights)
1174 1176 hmin_list = indexes[dataOut.heightList >= hmin]
1175 1177 hmax_list = indexes[dataOut.heightList <= hmax]
1176 1178
1177 1179 if hmin_list.any():
1178 1180 hmin_index = hmin_list[0]
1179 1181
1180 1182 if hmax_list.any():
1181 1183 hmax_index = hmax_list[-1]+1
1182 1184
1183 1185 x = dataOut.getTimeRange()
1184 1186
1185 1187 thisDatetime = dataOut.datatime
1186 1188
1187 1189 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
1188 1190 xlabel = "Local Time"
1189 1191 ylabel = "Phase (degrees)"
1190 1192
1191 1193 update_figfile = False
1192 1194
1193 1195 nplots = len(pairsIndexList)
1194 1196 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
1195 1197 phase_beacon = numpy.zeros(len(pairsIndexList))
1196 1198 for i in range(nplots):
1197 1199 pair = dataOut.pairsList[pairsIndexList[i]]
1198 1200 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
1199 1201 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
1200 1202 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
1201 1203 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
1202 1204 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
1203 1205
1204 1206 if dataOut.beacon_heiIndexList:
1205 1207 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
1206 1208 else:
1207 1209 phase_beacon[i] = numpy.average(phase)
1208 1210
1209 1211 if not self.isConfig:
1210 1212
1211 1213 nplots = len(pairsIndexList)
1212 1214
1213 1215 self.setup(id=id,
1214 1216 nplots=nplots,
1215 1217 wintitle=wintitle,
1216 1218 showprofile=showprofile,
1217 1219 show=show)
1218 1220
1219 1221 if timerange != None:
1220 1222 self.timerange = timerange
1221 1223
1222 1224 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
1223 1225
1224 1226 if ymin == None: ymin = 0
1225 1227 if ymax == None: ymax = 360
1226 1228
1227 1229 self.FTP_WEI = ftp_wei
1228 1230 self.EXP_CODE = exp_code
1229 1231 self.SUB_EXP_CODE = sub_exp_code
1230 1232 self.PLOT_POS = plot_pos
1231 1233
1232 1234 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
1233 1235 self.isConfig = True
1234 1236 self.figfile = figfile
1235 1237 self.xdata = numpy.array([])
1236 1238 self.ydata = numpy.array([])
1237 1239
1238 1240 update_figfile = True
1239 1241
1240 1242 #open file beacon phase
1241 1243 path = '%s%03d' %(self.PREFIX, self.id)
1242 1244 beacon_file = os.path.join(path,'%s.txt'%self.name)
1243 1245 self.filename_phase = os.path.join(figpath,beacon_file)
1244 1246 #self.save_phase(self.filename_phase)
1245 1247
1246 1248
1247 1249 #store data beacon phase
1248 1250 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
1249 1251
1250 1252 self.setWinTitle(title)
1251 1253
1252 1254
1253 1255 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
1254 1256
1255 1257 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
1256 1258
1257 1259 axes = self.axesList[0]
1258 1260
1259 1261 self.xdata = numpy.hstack((self.xdata, x[0:1]))
1260 1262
1261 1263 if len(self.ydata)==0:
1262 1264 self.ydata = phase_beacon.reshape(-1,1)
1263 1265 else:
1264 1266 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
1265 1267
1266 1268
1267 1269 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
1268 1270 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
1269 1271 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
1270 1272 XAxisAsTime=True, grid='both'
1271 1273 )
1272 1274
1273 1275 self.draw()
1274 1276
1275 1277 if dataOut.ltctime >= self.xmax:
1276 1278 self.counter_imagwr = wr_period
1277 1279 self.isConfig = False
1278 1280 update_figfile = True
1279 1281
1280 1282 self.save(figpath=figpath,
1281 1283 figfile=figfile,
1282 1284 save=save,
1283 1285 ftp=ftp,
1284 1286 wr_period=wr_period,
1285 1287 thisDatetime=thisDatetime,
1286 1288 update_figfile=update_figfile)
1287 1289
1288 1290 return dataOut
Show file before | Show file after | Hide comments
@@ -1,1287 +1,1285
1 1
2 2 import os
3 3 import time
4 4 import math
5 5 import datetime
6 6 import numpy
7 7 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator #YONG
8 8
9 9 from .jroplot_spectra import RTIPlot, NoisePlot
10 10
11 11 from schainpy.utils import log
12 12 from .plotting_codes import *
13 13
14 14 from schainpy.model.graphics.jroplot_base import Plot, plt
15 15
16 16 import matplotlib.pyplot as plt
17 17 import matplotlib.colors as colors
18 18 from matplotlib.ticker import MultipleLocator
19 19
20 20
21 21 class RTIDPPlot(RTIPlot):
22 22
23 23 '''Plot for RTI Double Pulse Experiment
24 24 '''
25 25
26 26 CODE = 'RTIDP'
27 27 colormap = 'jet'
28 28 plot_name = 'RTI'
29 29 plot_type = 'pcolorbuffer'
30 30
31 31 def setup(self):
32 32 self.xaxis = 'time'
33 33 self.ncols = 1
34 34 self.nrows = 3
35 35 self.nplots = self.nrows
36 36
37 37 self.ylabel = 'Range [km]'
38 38 self.xlabel = 'Time (LT)'
39 39
40 40 self.cb_label = 'Intensity (dB)'
41 41
42 42 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
43 43
44 44 self.titles = ['{} Channel {}'.format(
45 45 self.plot_name.upper(), '0x1'),'{} Channel {}'.format(
46 46 self.plot_name.upper(), '0'),'{} Channel {}'.format(
47 47 self.plot_name.upper(), '1')]
48 48
49 49 def update(self, dataOut):
50 50
51 51 data = {}
52 52 meta = {}
53 53 data['rti'] = dataOut.data_for_RTI_DP
54 54 data['NDP'] = dataOut.NDP
55 55
56 56 return data, meta
57 57
58 58 def plot(self):
59 59
60 60 NDP = self.data['NDP'][-1]
61 61 self.x = self.data.times
62 62 self.y = self.data.yrange[0:NDP]
63 63 self.z = self.data['rti']
64 64 self.z = numpy.ma.masked_invalid(self.z)
65 65
66 66 if self.decimation is None:
67 67 x, y, z = self.fill_gaps(self.x, self.y, self.z)
68 68 else:
69 69 x, y, z = self.fill_gaps(*self.decimate())
70 70
71 71 for n, ax in enumerate(self.axes):
72 72
73 73 self.zmax = self.zmax if self.zmax is not None else numpy.max(
74 74 self.z[1][0,12:40])
75 75 self.zmin = self.zmin if self.zmin is not None else numpy.min(
76 76 self.z[1][0,12:40])
77 77
78 78 if ax.firsttime:
79 79
80 80 if self.zlimits is not None:
81 81 self.zmin, self.zmax = self.zlimits[n]
82 82
83 83 ax.plt = ax.pcolormesh(x, y, z[n].T,
84 84 vmin=self.zmin,
85 85 vmax=self.zmax,
86 86 cmap=plt.get_cmap(self.colormap)
87 87 )
88 88 else:
89 89 #if self.zlimits is not None:
90 90 #self.zmin, self.zmax = self.zlimits[n]
91 91 ax.collections.remove(ax.collections[0])
92 92 ax.plt = ax.pcolormesh(x, y, z[n].T,
93 93 vmin=self.zmin,
94 94 vmax=self.zmax,
95 95 cmap=plt.get_cmap(self.colormap)
96 96 )
97 97
98 98
99 99 class RTILPPlot(RTIPlot):
100 100
101 101 '''
102 102 Plot for RTI Long Pulse
103 103 '''
104 104
105 105 CODE = 'RTILP'
106 106 colormap = 'jet'
107 107 plot_name = 'RTI LP'
108 108 plot_type = 'pcolorbuffer'
109 109
110 110 def setup(self):
111 111 self.xaxis = 'time'
112 112 self.ncols = 1
113 113 self.nrows = 4
114 114 self.nplots = self.nrows
115 115
116 116 self.ylabel = 'Range [km]'
117 117 self.xlabel = 'Time (LT)'
118 118
119 119 self.cb_label = 'Intensity (dB)'
120 120
121 121 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.1, 'right':0.95})
122 122
123 123
124 124 self.titles = ['{} Channel {}'.format(
125 125 self.plot_name.upper(), '0'),'{} Channel {}'.format(
126 126 self.plot_name.upper(), '1'),'{} Channel {}'.format(
127 127 self.plot_name.upper(), '2'),'{} Channel {}'.format(
128 128 self.plot_name.upper(), '3')]
129 129
130 130
131 131 def update(self, dataOut):
132 132
133 133 data = {}
134 134 meta = {}
135 135 data['rti'] = dataOut.data_for_RTI_LP
136 136 data['NRANGE'] = dataOut.NRANGE
137 137
138 138 return data, meta
139 139
140 140 def plot(self):
141 141
142 142 NRANGE = self.data['NRANGE'][-1]
143 143 self.x = self.data.times
144 144 self.y = self.data.yrange[0:NRANGE]
145 145
146 146 self.z = self.data['rti']
147 147
148 148 self.z = numpy.ma.masked_invalid(self.z)
149 149
150 150 if self.decimation is None:
151 151 x, y, z = self.fill_gaps(self.x, self.y, self.z)
152 152 else:
153 153 x, y, z = self.fill_gaps(*self.decimate())
154 154
155 155 for n, ax in enumerate(self.axes):
156 156
157 157 self.zmax = self.zmax if self.zmax is not None else numpy.max(
158 158 self.z[1][0,12:40])
159 159 self.zmin = self.zmin if self.zmin is not None else numpy.min(
160 160 self.z[1][0,12:40])
161 161
162 162 if ax.firsttime:
163 163
164 164 if self.zlimits is not None:
165 165 self.zmin, self.zmax = self.zlimits[n]
166 166
167 167
168 168 ax.plt = ax.pcolormesh(x, y, z[n].T,
169 169 vmin=self.zmin,
170 170 vmax=self.zmax,
171 171 cmap=plt.get_cmap(self.colormap)
172 172 )
173 173
174 174 else:
175 175 #if self.zlimits is not None:
176 176 #self.zmin, self.zmax = self.zlimits[n]
177 177 ax.collections.remove(ax.collections[0])
178 178 ax.plt = ax.pcolormesh(x, y, z[n].T,
179 179 vmin=self.zmin,
180 180 vmax=self.zmax,
181 181 cmap=plt.get_cmap(self.colormap)
182 182 )
183 183
184 184
185 185 class DenRTIPlot(RTIPlot):
186 186
187 187 '''
188 188 Plot for Den
189 189 '''
190 190
191 191 CODE = 'denrti'
192 192 colormap = 'jet'
193 193
194 194 def setup(self):
195 195 self.xaxis = 'time'
196 196 self.ncols = 1
197 197 self.nrows = self.data.shape(self.CODE)[0]
198 198 self.nplots = self.nrows
199 199
200 200 self.ylabel = 'Range [km]'
201 201 self.xlabel = 'Time (LT)'
202 202
203 203 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
204 204
205 205 if self.CODE == 'denrti':
206 206 self.cb_label = r'$\mathrm{N_e}$ Electron Density ($\mathrm{1/cm^3}$)'
207 207
208 208
209 209 self.titles = ['Electron Density RTI']
210 210
211 211 def update(self, dataOut):
212 212
213 213 data = {}
214 214 meta = {}
215 215
216 216 data['denrti'] = dataOut.DensityFinal
217 217
218 218 return data, meta
219 219
220 220 def plot(self):
221 221
222 222 self.x = self.data.times
223 223 self.y = self.data.yrange
224 224
225 225 self.z = self.data[self.CODE]
226 226
227 227 self.z = numpy.ma.masked_invalid(self.z)
228 228
229 229 if self.decimation is None:
230 230 x, y, z = self.fill_gaps(self.x, self.y, self.z)
231 231 else:
232 232 x, y, z = self.fill_gaps(*self.decimate())
233 233
234 234 for n, ax in enumerate(self.axes):
235 235
236 236 self.zmax = self.zmax if self.zmax is not None else numpy.max(
237 237 self.z[n])
238 238 self.zmin = self.zmin if self.zmin is not None else numpy.min(
239 239 self.z[n])
240 240
241 241 if ax.firsttime:
242 242
243 243 if self.zlimits is not None:
244 244 self.zmin, self.zmax = self.zlimits[n]
245 245 if numpy.log10(self.zmin)<0:
246 246 self.zmin=1
247 247 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
248 248 vmin=self.zmin,
249 249 vmax=self.zmax,
250 250 cmap=self.cmaps[n],
251 251 norm=colors.LogNorm()
252 252 )
253 253
254 254 else:
255 255 if self.zlimits is not None:
256 256 self.zmin, self.zmax = self.zlimits[n]
257 257 ax.collections.remove(ax.collections[0])
258 258 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
259 259 vmin=self.zmin,
260 260 vmax=self.zmax,
261 261 cmap=self.cmaps[n],
262 262 norm=colors.LogNorm()
263 263 )
264 264
265 265
266
267
268 266 class ETempRTIPlot(RTIPlot):
269 267
270 268 '''
271 269 Plot for Electron Temperature
272 270 '''
273 271
274 272 CODE = 'ETemp'
275 273 colormap = 'jet'
276 274
277 275 def setup(self):
278 276 self.xaxis = 'time'
279 277 self.ncols = 1
280 278 self.nrows = self.data.shape(self.CODE)[0]
281 279 self.nplots = self.nrows
282 280
283 281 self.ylabel = 'Range [km]'
284 282 self.xlabel = 'Time (LT)'
285 283 self.plots_adjust.update({'wspace': 0.8, 'hspace':0.2, 'left': 0.2, 'right': 0.9, 'bottom': 0.18})
286 284 if self.CODE == 'ETemp':
287 285 self.cb_label = 'Electron Temperature (K)'
288 286 self.titles = ['Electron Temperature RTI']
289 287 if self.CODE == 'ITemp':
290 288 self.cb_label = 'Ion Temperature (K)'
291 289 self.titles = ['Ion Temperature RTI']
292 290 if self.CODE == 'HeFracLP':
293 291 self.cb_label='He+ Fraction'
294 292 self.titles = ['He+ Fraction RTI']
295 293 self.zmax=0.16
296 294 if self.CODE== 'HFracLP':
297 295 self.cb_label='H+ Fraction'
298 296 self.titles = ['H+ Fraction RTI']
299 297
300 298 def update(self, dataOut):
301 299
302 300 data = {}
303 301 meta = {}
304 302
305 303 data['ETemp'] = dataOut.ElecTempFinal
306 304
307 305 return data, meta
308 306
309 307 def plot(self):
310 308
311 309 self.x = self.data.times
312 310 self.y = self.data.yrange
313 311
314 312
315 313 self.z = self.data[self.CODE]
316 314
317 315 self.z = numpy.ma.masked_invalid(self.z)
318 316
319 317 if self.decimation is None:
320 318 x, y, z = self.fill_gaps(self.x, self.y, self.z)
321 319 else:
322 320 x, y, z = self.fill_gaps(*self.decimate())
323 321
324 322 for n, ax in enumerate(self.axes):
325 323
326 324 self.zmax = self.zmax if self.zmax is not None else numpy.max(
327 325 self.z[n])
328 326 self.zmin = self.zmin if self.zmin is not None else numpy.min(
329 327 self.z[n])
330 328
331 329 if ax.firsttime:
332 330
333 331 if self.zlimits is not None:
334 332 self.zmin, self.zmax = self.zlimits[n]
335 333
336 334 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
337 335 vmin=self.zmin,
338 336 vmax=self.zmax,
339 337 cmap=self.cmaps[n]
340 338 )
341 339 #plt.tight_layout()
342 340
343 341 else:
344 342 if self.zlimits is not None:
345 343 self.zmin, self.zmax = self.zlimits[n]
346 344 ax.collections.remove(ax.collections[0])
347 345 ax.plt = ax.pcolormesh(x, y, z[n].T * self.factors[n],
348 346 vmin=self.zmin,
349 347 vmax=self.zmax,
350 348 cmap=self.cmaps[n]
351 349 )
352 350
353 351
354
355 352 class ITempRTIPlot(ETempRTIPlot):
356 353
357 354 '''
358 355 Plot for Ion Temperature
359 356 '''
360 357
361 358 CODE = 'ITemp'
362 359 colormap = 'jet'
363 360 plot_name = 'Ion Temperature'
364 361
365 362 def update(self, dataOut):
366 363
367 364 data = {}
368 365 meta = {}
369 366
370 367 data['ITemp'] = dataOut.IonTempFinal
371 368
372 369 return data, meta
373 370
374 371
375
376 372 class HFracRTIPlot(ETempRTIPlot):
377 373
378 374 '''
379 375 Plot for H+ LP
380 376 '''
381 377
382 378 CODE = 'HFracLP'
383 379 colormap = 'jet'
384 380 plot_name = 'H+ Frac'
385 381
386 382 def update(self, dataOut):
387 383
388 384 data = {}
389 385 meta = {}
390 386 data['HFracLP'] = dataOut.PhyFinal
391 387
392 388 return data, meta
393 389
394 390
395 391 class HeFracRTIPlot(ETempRTIPlot):
396 392
397 393 '''
398 394 Plot for He+ LP
399 395 '''
400 396
401 397 CODE = 'HeFracLP'
402 398 colormap = 'jet'
403 399 plot_name = 'He+ Frac'
404 400
405 401 def update(self, dataOut):
406 402
407 403 data = {}
408 404 meta = {}
409 405 data['HeFracLP'] = dataOut.PheFinal
410 406
411 407 return data, meta
412 408
413 409
414 410 class TempsDPPlot(Plot):
415 411 '''
416 412 Plot for Electron - Ion Temperatures
417 413 '''
418 414
419 415 CODE = 'tempsDP'
420 416 #plot_name = 'Temperatures'
421 417 plot_type = 'scatterbuffer'
422 418
423 419 def setup(self):
424 420
425 421 self.ncols = 1
426 422 self.nrows = 1
427 423 self.nplots = 1
428 424 self.ylabel = 'Range [km]'
429 425 self.xlabel = 'Temperature (K)'
430 426 self.titles = ['Electron/Ion Temperatures']
431 427 self.width = 3.5
432 428 self.height = 5.5
433 429 self.colorbar = False
434 430 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
435 431
436 432 def update(self, dataOut):
437 433 data = {}
438 434 meta = {}
439 435
440 436 data['Te'] = dataOut.te2
441 437 data['Ti'] = dataOut.ti2
442 438 data['Te_error'] = dataOut.ete2
443 439 data['Ti_error'] = dataOut.eti2
444 440
445 441 meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
446 442
447 443 return data, meta
448 444
449 445 def plot(self):
450 446
451 447 y = self.data.yrange
452 448
453 449 self.xmin = -100
454 450 self.xmax = 5000
455 451
456 452 ax = self.axes[0]
457 453
458 454 data = self.data[-1]
459 455
460 456 Te = data['Te']
461 457 Ti = data['Ti']
462 458 errTe = data['Te_error']
463 459 errTi = data['Ti_error']
464 460
465 461 if ax.firsttime:
466 462 ax.errorbar(Te, y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
467 463 ax.errorbar(Ti, y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
468 464 plt.legend(loc='lower right')
469 465 self.ystep_given = 50
470 466 ax.yaxis.set_minor_locator(MultipleLocator(15))
471 467 ax.grid(which='minor')
472 468
473 469 else:
474 470 self.clear_figures()
475 471 ax.errorbar(Te, y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
476 472 ax.errorbar(Ti, y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
477 473 plt.legend(loc='lower right')
478 474 ax.yaxis.set_minor_locator(MultipleLocator(15))
479 475
480 476
481 477 class TempsHPPlot(Plot):
482 478 '''
483 479 Plot for Temperatures Hybrid Experiment
484 480 '''
485 481
486 482 CODE = 'temps_LP'
487 483 #plot_name = 'Temperatures'
488 484 plot_type = 'scatterbuffer'
489 485
490 486
491 487 def setup(self):
492 488
493 489 self.ncols = 1
494 490 self.nrows = 1
495 491 self.nplots = 1
496 492 self.ylabel = 'Range [km]'
497 493 self.xlabel = 'Temperature (K)'
498 494 self.titles = ['Electron/Ion Temperatures']
499 495 self.width = 3.5
500 496 self.height = 6.5
501 497 self.colorbar = False
502 498 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
503 499
504 500 def update(self, dataOut):
505 501 data = {}
506 502 meta = {}
507 503
508 504
509 505 data['Te'] = numpy.concatenate((dataOut.te2[:dataOut.cut],dataOut.te[dataOut.cut:]))
510 506 data['Ti'] = numpy.concatenate((dataOut.ti2[:dataOut.cut],dataOut.ti[dataOut.cut:]))
511 507 data['Te_error'] = numpy.concatenate((dataOut.ete2[:dataOut.cut],dataOut.ete[dataOut.cut:]))
512 508 data['Ti_error'] = numpy.concatenate((dataOut.eti2[:dataOut.cut],dataOut.eti[dataOut.cut:]))
513 509
514 510 meta['yrange'] = dataOut.heightList[0:dataOut.NACF]
515 511
516 512 return data, meta
517 513
518 514 def plot(self):
519 515
520 516
521 517 self.y = self.data.yrange
522 518 self.xmin = -100
523 519 self.xmax = 4500
524 520 ax = self.axes[0]
525 521
526 522 data = self.data[-1]
527 523
528 524 Te = data['Te']
529 525 Ti = data['Ti']
530 526 errTe = data['Te_error']
531 527 errTi = data['Ti_error']
532 528
533 529 if ax.firsttime:
534 530
535 531 ax.errorbar(Te, self.y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
536 532 ax.errorbar(Ti, self.y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
537 533 plt.legend(loc='lower right')
538 534 self.ystep_given = 200
539 535 ax.yaxis.set_minor_locator(MultipleLocator(15))
540 536 ax.grid(which='minor')
541 537
542 538 else:
543 539 self.clear_figures()
544 540 ax.errorbar(Te, self.y, xerr=errTe, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='Te')
545 541 ax.errorbar(Ti, self.y, fmt='k^', xerr=errTi,elinewidth=1.0,color='b',linewidth=2.0, label='Ti')
546 542 plt.legend(loc='lower right')
547 543 ax.yaxis.set_minor_locator(MultipleLocator(15))
544 ax.grid(which='minor')
545
548 546
549 547 class FracsHPPlot(Plot):
550 548 '''
551 549 Plot for Composition LP
552 550 '''
553 551
554 552 CODE = 'fracs_LP'
555 553 plot_type = 'scatterbuffer'
556 554
557 555
558 556 def setup(self):
559 557
560 558 self.ncols = 1
561 559 self.nrows = 1
562 560 self.nplots = 1
563 561 self.ylabel = 'Range [km]'
564 562 self.xlabel = 'Frac'
565 563 self.titles = ['Composition']
566 564 self.width = 3.5
567 565 self.height = 6.5
568 566 self.colorbar = False
569 567 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
570 568
571 569 def update(self, dataOut):
572 570 data = {}
573 571 meta = {}
574 572
575 573 #aux_nan=numpy.zeros(dataOut.cut,'float32')
576 574 #aux_nan[:]=numpy.nan
577 575 #data['ph'] = numpy.concatenate((aux_nan,dataOut.ph[dataOut.cut:]))
578 576 #data['eph'] = numpy.concatenate((aux_nan,dataOut.eph[dataOut.cut:]))
579 577
580 578 data['ph'] = dataOut.ph[dataOut.cut:]
581 579 data['eph'] = dataOut.eph[dataOut.cut:]
582 580 data['phe'] = dataOut.phe[dataOut.cut:]
583 581 data['ephe'] = dataOut.ephe[dataOut.cut:]
584 582
585 583 data['cut'] = dataOut.cut
586 584
587 585 meta['yrange'] = dataOut.heightList[0:dataOut.NACF]
588 586
589 587
590 588 return data, meta
591 589
592 590 def plot(self):
593 591
594 592 data = self.data[-1]
595 593
596 594 ph = data['ph']
597 595 eph = data['eph']
598 596 phe = data['phe']
599 597 ephe = data['ephe']
600 598 cut = data['cut']
601 599 self.y = self.data.yrange
602 600
603 601 self.xmin = 0
604 602 self.xmax = 1
605 603 ax = self.axes[0]
606 604
607 605 if ax.firsttime:
608 606
609 607 ax.errorbar(ph, self.y[cut:], xerr=eph, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='H+')
610 608 ax.errorbar(phe, self.y[cut:], fmt='k^', xerr=ephe,elinewidth=1.0,color='b',linewidth=2.0, label='He+')
611 609 plt.legend(loc='lower right')
612 610 self.xstep_given = 0.2
613 611 self.ystep_given = 200
614 612 ax.yaxis.set_minor_locator(MultipleLocator(15))
615 613 ax.grid(which='minor')
616 614
617 615 else:
618 616 self.clear_figures()
619 617 ax.errorbar(ph, self.y[cut:], xerr=eph, fmt='r^',elinewidth=1.0,color='b',linewidth=2.0, label='H+')
620 618 ax.errorbar(phe, self.y[cut:], fmt='k^', xerr=ephe,elinewidth=1.0,color='b',linewidth=2.0, label='He+')
621 619 plt.legend(loc='lower right')
622 620 ax.yaxis.set_minor_locator(MultipleLocator(15))
623 621
624 622 class EDensityPlot(Plot):
625 623 '''
626 624 Plot for electron density
627 625 '''
628 626
629 627 CODE = 'den'
630 628 #plot_name = 'Electron Density'
631 629 plot_type = 'scatterbuffer'
632 630
633 631 def setup(self):
634 632
635 633 self.ncols = 1
636 634 self.nrows = 1
637 635 self.nplots = 1
638 636 self.ylabel = 'Range [km]'
639 637 self.xlabel = r'$\mathrm{N_e}$ Electron Density ($\mathrm{1/cm^3}$)'
640 638 self.titles = ['Electron Density']
641 639 self.width = 3.5
642 640 self.height = 5.5
643 641 self.colorbar = False
644 642 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
645 643
646 644 def update(self, dataOut):
647 645 data = {}
648 646 meta = {}
649 647
650 648 data['den_power'] = dataOut.ph2[:dataOut.NSHTS]
651 649 data['den_Faraday'] = dataOut.dphi[:dataOut.NSHTS]
652 650 data['den_error'] = dataOut.sdp2[:dataOut.NSHTS]
653 651 #data['err_Faraday'] = dataOut.sdn1[:dataOut.NSHTS]
654 652
655 653 data['NSHTS'] = dataOut.NSHTS
656 654
657 655 meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
658 656
659 657 return data, meta
660 658
661 659 def plot(self):
662 660
663 661 y = self.data.yrange
664 662
665 663 self.xmin = 1e3
666 664 self.xmax = 1e7
667 665
668 666 ax = self.axes[0]
669 667
670 668 data = self.data[-1]
671 669
672 670 DenPow = data['den_power']
673 671 DenFar = data['den_Faraday']
674 672 errDenPow = data['den_error']
675 673 #errFaraday = data['err_Faraday']
676 674
677 675 NSHTS = data['NSHTS']
678 676
679 677 if self.CODE == 'denLP':
680 678 DenPowLP = data['den_LP']
681 679 errDenPowLP = data['den_LP_error']
682 680 cut = data['cut']
683 681
684 682 if ax.firsttime:
685 683 self.autoxticks=False
686 684 #ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
687 685 ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday',markersize=2)
688 686 #ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
689 687 ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2)
690 688
691 689 if self.CODE=='denLP':
692 690 ax.errorbar(DenPowLP[cut:], y[cut:], xerr=errDenPowLP[cut:], fmt='r^-',elinewidth=1.0,color='r',linewidth=1.0, label='LP Profile',markersize=2)
693 691
694 692 plt.legend(loc='upper left',fontsize=8.5)
695 693 #plt.legend(loc='lower left',fontsize=8.5)
696 694 ax.set_xscale("log", nonposx='clip')
697 695 grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50)
698 696 self.ystep_given=100
699 697 if self.CODE=='denLP':
700 698 self.ystep_given=200
701 699 ax.set_yticks(grid_y_ticks,minor=True)
702 700 ax.grid(which='minor')
703 701
704 702 else:
705 703 dataBefore = self.data[-2]
706 704 DenPowBefore = dataBefore['den_power']
707 705 self.clear_figures()
708 706 #ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday Profile',markersize=2)
709 707 ax.errorbar(DenFar, y[:NSHTS], xerr=1, fmt='h-',elinewidth=1.0,color='g',linewidth=1.0, label='Faraday',markersize=2)
710 708 #ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power Profile',markersize=2)
711 709 ax.errorbar(DenPow, y[:NSHTS], fmt='k^-', xerr=errDenPow,elinewidth=1.0,color='b',linewidth=1.0, label='Power',markersize=2)
712 710 ax.errorbar(DenPowBefore, y[:NSHTS], elinewidth=1.0,color='r',linewidth=0.5,linestyle="dashed")
713 711
714 712 if self.CODE=='denLP':
715 713 ax.errorbar(DenPowLP[cut:], y[cut:], fmt='r^-', xerr=errDenPowLP[cut:],elinewidth=1.0,color='r',linewidth=1.0, label='LP Profile',markersize=2)
716 714
717 715 ax.set_xscale("log", nonposx='clip')
718 716 grid_y_ticks=numpy.arange(numpy.nanmin(y),numpy.nanmax(y),50)
719 717 ax.set_yticks(grid_y_ticks,minor=True)
720 718 ax.grid(which='minor')
721 719 plt.legend(loc='upper left',fontsize=8.5)
722 720 #plt.legend(loc='lower left',fontsize=8.5)
723 721
724 722 class FaradayAnglePlot(Plot):
725 723 '''
726 724 Plot for electron density
727 725 '''
728 726
729 727 CODE = 'angle'
730 728 plot_name = 'Faraday Angle'
731 729 plot_type = 'scatterbuffer'
732 730
733 731 def setup(self):
734 732
735 733 self.ncols = 1
736 734 self.nrows = 1
737 735 self.nplots = 1
738 736 self.ylabel = 'Range [km]'
739 737 self.xlabel = 'Faraday Angle (º)'
740 738 self.titles = ['Electron Density']
741 739 self.width = 3.5
742 740 self.height = 5.5
743 741 self.colorbar = False
744 742 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
745 743
746 744 def update(self, dataOut):
747 745 data = {}
748 746 meta = {}
749 747
750 748 data['angle'] = numpy.degrees(dataOut.phi)
751 749 #'''
752 750 print(dataOut.phi_uwrp)
753 751 print(data['angle'])
754 752 exit(1)
755 753 #'''
756 754 data['dphi'] = dataOut.dphi_uc*10
757 755 #print(dataOut.dphi)
758 756
759 757 #data['NSHTS'] = dataOut.NSHTS
760 758
761 759 #meta['yrange'] = dataOut.heightList[0:dataOut.NSHTS]
762 760
763 761 return data, meta
764 762
765 763 def plot(self):
766 764
767 765 data = self.data[-1]
768 766 self.x = data[self.CODE]
769 767 dphi = data['dphi']
770 768 self.y = self.data.yrange
771 769 self.xmin = -360#-180
772 770 self.xmax = 360#180
773 771 ax = self.axes[0]
774 772
775 773 if ax.firsttime:
776 774 self.autoxticks=False
777 775 #if self.CODE=='den':
778 776 ax.plot(self.x, self.y,marker='o',color='g',linewidth=1.0,markersize=2)
779 777 ax.plot(dphi, self.y,marker='o',color='blue',linewidth=1.0,markersize=2)
780 778
781 779 grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
782 780 self.ystep_given=100
783 781 if self.CODE=='denLP':
784 782 self.ystep_given=200
785 783 ax.set_yticks(grid_y_ticks,minor=True)
786 784 ax.grid(which='minor')
787 785 #plt.tight_layout()
788 786 else:
789 787
790 788 self.clear_figures()
791 789 #if self.CODE=='den':
792 790 #print(numpy.shape(self.x))
793 791 ax.plot(self.x, self.y, marker='o',color='g',linewidth=1.0, markersize=2)
794 792 ax.plot(dphi, self.y,marker='o',color='blue',linewidth=1.0,markersize=2)
795 793
796 794 grid_y_ticks=numpy.arange(numpy.nanmin(self.y),numpy.nanmax(self.y),50)
797 795 ax.set_yticks(grid_y_ticks,minor=True)
798 796 ax.grid(which='minor')
799 797
800 798 class EDensityHPPlot(EDensityPlot):
801 799
802 800 '''
803 801 Plot for Electron Density Hybrid Experiment
804 802 '''
805 803
806 804 CODE = 'denLP'
807 805 plot_name = 'Electron Density'
808 806 plot_type = 'scatterbuffer'
809 807
810 808 def update(self, dataOut):
811 809 data = {}
812 810 meta = {}
813 811
814 812 data['den_power'] = dataOut.ph2[:dataOut.NSHTS]
815 813 data['den_Faraday']=dataOut.dphi[:dataOut.NSHTS]
816 814 data['den_error']=dataOut.sdp2[:dataOut.NSHTS]
817 815 data['den_LP']=dataOut.ne[:dataOut.NACF]
818 816 data['den_LP_error']=dataOut.ene[:dataOut.NACF]*dataOut.ne[:dataOut.NACF]*0.434
819 817 #self.ene=10**dataOut.ene[:dataOut.NACF]
820 818 data['NSHTS']=dataOut.NSHTS
821 819 data['cut']=dataOut.cut
822 820
823 821 return data, meta
824 822
825 823
826 824 class ACFsPlot(Plot):
827 825 '''
828 826 Plot for ACFs Double Pulse Experiment
829 827 '''
830 828
831 829 CODE = 'acfs'
832 830 #plot_name = 'ACF'
833 831 plot_type = 'scatterbuffer'
834 832
835 833
836 834 def setup(self):
837 835 self.ncols = 1
838 836 self.nrows = 1
839 837 self.nplots = 1
840 838 self.ylabel = 'Range [km]'
841 839 self.xlabel = 'Lag (ms)'
842 840 self.titles = ['ACFs']
843 841 self.width = 3.5
844 842 self.height = 5.5
845 843 self.colorbar = False
846 844 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
847 845
848 846 def update(self, dataOut):
849 847 data = {}
850 848 meta = {}
851 849
852 850 data['ACFs'] = dataOut.acfs_to_plot
853 851 data['ACFs_error'] = dataOut.acfs_error_to_plot
854 852 data['lags'] = dataOut.lags_to_plot
855 853 data['Lag_contaminated_1'] = dataOut.x_igcej_to_plot
856 854 data['Lag_contaminated_2'] = dataOut.x_ibad_to_plot
857 855 data['Height_contaminated_1'] = dataOut.y_igcej_to_plot
858 856 data['Height_contaminated_2'] = dataOut.y_ibad_to_plot
859 857
860 858 meta['yrange'] = numpy.array([])
861 859 #meta['NSHTS'] = dataOut.NSHTS
862 860 #meta['DPL'] = dataOut.DPL
863 861 data['NSHTS'] = dataOut.NSHTS #This is metadata
864 862 data['DPL'] = dataOut.DPL #This is metadata
865 863
866 864 return data, meta
867 865
868 866 def plot(self):
869 867
870 868 data = self.data[-1]
871 869 #NSHTS = self.meta['NSHTS']
872 870 #DPL = self.meta['DPL']
873 871 NSHTS = data['NSHTS'] #This is metadata
874 872 DPL = data['DPL'] #This is metadata
875 873
876 874 lags = data['lags']
877 875 ACFs = data['ACFs']
878 876 errACFs = data['ACFs_error']
879 877 BadLag1 = data['Lag_contaminated_1']
880 878 BadLag2 = data['Lag_contaminated_2']
881 879 BadHei1 = data['Height_contaminated_1']
882 880 BadHei2 = data['Height_contaminated_2']
883 881
884 882 self.xmin = 0.0
885 883 self.xmax = 2.0
886 884 self.y = ACFs
887 885
888 886 ax = self.axes[0]
889 887
890 888 if ax.firsttime:
891 889
892 890 for i in range(NSHTS):
893 891 x_aux = numpy.isfinite(lags[i,:])
894 892 y_aux = numpy.isfinite(ACFs[i,:])
895 893 yerr_aux = numpy.isfinite(errACFs[i,:])
896 894 x_igcej_aux = numpy.isfinite(BadLag1[i,:])
897 895 y_igcej_aux = numpy.isfinite(BadHei1[i,:])
898 896 x_ibad_aux = numpy.isfinite(BadLag2[i,:])
899 897 y_ibad_aux = numpy.isfinite(BadHei2[i,:])
900 898 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
901 899 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',marker='o',linewidth=1.0,markersize=2)
902 900 ax.plot(BadLag1[i,x_igcej_aux],BadHei1[i,y_igcej_aux],'x',color='red',markersize=2)
903 901 ax.plot(BadLag2[i,x_ibad_aux],BadHei2[i,y_ibad_aux],'X',color='red',markersize=2)
904 902
905 903 self.xstep_given = (self.xmax-self.xmin)/(DPL-1)
906 904 self.ystep_given = 50
907 905 ax.yaxis.set_minor_locator(MultipleLocator(15))
908 906 ax.grid(which='minor')
909 907
910 908 else:
911 909 self.clear_figures()
912 910 for i in range(NSHTS):
913 911 x_aux = numpy.isfinite(lags[i,:])
914 912 y_aux = numpy.isfinite(ACFs[i,:])
915 913 yerr_aux = numpy.isfinite(errACFs[i,:])
916 914 x_igcej_aux = numpy.isfinite(BadLag1[i,:])
917 915 y_igcej_aux = numpy.isfinite(BadHei1[i,:])
918 916 x_ibad_aux = numpy.isfinite(BadLag2[i,:])
919 917 y_ibad_aux = numpy.isfinite(BadHei2[i,:])
920 918 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
921 919 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],linewidth=1.0,markersize=2,color='b',marker='o')
922 920 ax.plot(BadLag1[i,x_igcej_aux],BadHei1[i,y_igcej_aux],'x',color='red',markersize=2)
923 921 ax.plot(BadLag2[i,x_ibad_aux],BadHei2[i,y_ibad_aux],'X',color='red',markersize=2)
924 922 ax.yaxis.set_minor_locator(MultipleLocator(15))
925 923
926 924 class ACFsLPPlot(Plot):
927 925 '''
928 926 Plot for ACFs Double Pulse Experiment
929 927 '''
930 928
931 929 CODE = 'acfs_LP'
932 930 #plot_name = 'ACF'
933 931 plot_type = 'scatterbuffer'
934 932
935 933
936 934 def setup(self):
937 935 self.ncols = 1
938 936 self.nrows = 1
939 937 self.nplots = 1
940 938 self.ylabel = 'Range [km]'
941 939 self.xlabel = 'Lag (ms)'
942 940 self.titles = ['ACFs']
943 941 self.width = 3.5
944 942 self.height = 5.5
945 943 self.colorbar = False
946 944 self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
947 945
948 946 def update(self, dataOut):
949 947 data = {}
950 948 meta = {}
951 949
952 950 aux=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
953 951 errors=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
954 952 lags_LP_to_plot=numpy.zeros((dataOut.NACF,dataOut.IBITS),'float32')
955 953
956 954 for i in range(dataOut.NACF):
957 955 for j in range(dataOut.IBITS):
958 956 if numpy.abs(dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0])<1.0:
959 957 aux[i,j]=dataOut.output_LP_integrated.real[j,i,0]/dataOut.output_LP_integrated.real[0,i,0]
960 958 aux[i,j]=max(min(aux[i,j],1.0),-1.0)*dataOut.DH+dataOut.heightList[i]
961 959 lags_LP_to_plot[i,j]=dataOut.lags_LP[j]
962 960 errors[i,j]=dataOut.errors[j,i]/dataOut.output_LP_integrated.real[0,i,0]*dataOut.DH
963 961 else:
964 962 aux[i,j]=numpy.nan
965 963 lags_LP_to_plot[i,j]=numpy.nan
966 964 errors[i,j]=numpy.nan
967 965
968 966 data['ACFs'] = aux
969 967 data['ACFs_error'] = errors
970 968 data['lags'] = lags_LP_to_plot
971 969
972 970 meta['yrange'] = numpy.array([])
973 971 #meta['NACF'] = dataOut.NACF
974 972 #meta['NLAG'] = dataOut.NLAG
975 973 data['NACF'] = dataOut.NACF #This is metadata
976 974 data['NLAG'] = dataOut.NLAG #This is metadata
977 975
978 976 return data, meta
979 977
980 978 def plot(self):
981 979
982 980 data = self.data[-1]
983 981 #NACF = self.meta['NACF']
984 982 #NLAG = self.meta['NLAG']
985 983 NACF = data['NACF'] #This is metadata
986 984 NLAG = data['NLAG'] #This is metadata
987 985
988 986 lags = data['lags']
989 987 ACFs = data['ACFs']
990 988 errACFs = data['ACFs_error']
991 989
992 990 self.xmin = 0.0
993 991 self.xmax = 1.5
994 992
995 993 self.y = ACFs
996 994
997 995 ax = self.axes[0]
998 996
999 997 if ax.firsttime:
1000 998
1001 999 for i in range(NACF):
1002 1000 x_aux = numpy.isfinite(lags[i,:])
1003 1001 y_aux = numpy.isfinite(ACFs[i,:])
1004 1002 yerr_aux = numpy.isfinite(errACFs[i,:])
1005 1003
1006 1004 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
1007 1005 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',linewidth=1.0,markersize=2,ecolor='r')
1008 1006
1009 1007 #self.xstep_given = (self.xmax-self.xmin)/(self.data.NLAG-1)
1010 1008 self.xstep_given=0.3
1011 1009 self.ystep_given = 200
1012 1010 ax.yaxis.set_minor_locator(MultipleLocator(15))
1013 1011 ax.grid(which='minor')
1014 1012
1015 1013 else:
1016 1014 self.clear_figures()
1017 1015
1018 1016 for i in range(NACF):
1019 1017 x_aux = numpy.isfinite(lags[i,:])
1020 1018 y_aux = numpy.isfinite(ACFs[i,:])
1021 1019 yerr_aux = numpy.isfinite(errACFs[i,:])
1022 1020
1023 1021 if lags[i,:][~numpy.isnan(lags[i,:])].shape[0]>2:
1024 1022 ax.errorbar(lags[i,x_aux], ACFs[i,y_aux], yerr=errACFs[i,x_aux],color='b',linewidth=1.0,markersize=2,ecolor='r')
1025 1023
1026 1024 ax.yaxis.set_minor_locator(MultipleLocator(15))
1027 1025
1028 1026
1029 1027 class CrossProductsPlot(Plot):
1030 1028 '''
1031 1029 Plot for cross products
1032 1030 '''
1033 1031
1034 1032 CODE = 'crossprod'
1035 1033 plot_name = 'Cross Products'
1036 1034 plot_type = 'scatterbuffer'
1037 1035
1038 1036 def setup(self):
1039 1037
1040 1038 self.ncols = 3
1041 1039 self.nrows = 1
1042 1040 self.nplots = 3
1043 1041 self.ylabel = 'Range [km]'
1044 1042 self.titles = []
1045 1043 self.width = 3.5*self.nplots
1046 1044 self.height = 5.5
1047 1045 self.colorbar = False
1048 1046 self.plots_adjust.update({'wspace':.3, 'left': 0.12, 'right': 0.92, 'bottom': 0.1})
1049 1047
1050 1048
1051 1049 def update(self, dataOut):
1052 1050
1053 1051 data = {}
1054 1052 meta = {}
1055 1053
1056 1054 data['crossprod'] = dataOut.crossprods
1057 1055 data['NDP'] = dataOut.NDP
1058 1056
1059 1057 return data, meta
1060 1058
1061 1059 def plot(self):
1062 1060
1063 1061 NDP = self.data['NDP'][-1]
1064 1062 x = self.data['crossprod'][:,-1,:,:,:,:]
1065 1063 y = self.data.yrange[0:NDP]
1066 1064
1067 1065 for n, ax in enumerate(self.axes):
1068 1066
1069 1067 self.xmin=numpy.min(numpy.concatenate((x[n][0,20:30,0,0],x[n][1,20:30,0,0],x[n][2,20:30,0,0],x[n][3,20:30,0,0])))
1070 1068 self.xmax=numpy.max(numpy.concatenate((x[n][0,20:30,0,0],x[n][1,20:30,0,0],x[n][2,20:30,0,0],x[n][3,20:30,0,0])))
1071 1069
1072 1070 if ax.firsttime:
1073 1071
1074 1072 self.autoxticks=False
1075 1073 if n==0:
1076 1074 label1='kax'
1077 1075 label2='kay'
1078 1076 label3='kbx'
1079 1077 label4='kby'
1080 1078 self.xlimits=[(self.xmin,self.xmax)]
1081 1079 elif n==1:
1082 1080 label1='kax2'
1083 1081 label2='kay2'
1084 1082 label3='kbx2'
1085 1083 label4='kby2'
1086 1084 self.xlimits.append((self.xmin,self.xmax))
1087 1085 elif n==2:
1088 1086 label1='kaxay'
1089 1087 label2='kbxby'
1090 1088 label3='kaxbx'
1091 1089 label4='kaxby'
1092 1090 self.xlimits.append((self.xmin,self.xmax))
1093 1091
1094 1092 ax.plotline1 = ax.plot(x[n][0,:,0,0], y, color='r',linewidth=2.0, label=label1)
1095 1093 ax.plotline2 = ax.plot(x[n][1,:,0,0], y, color='k',linewidth=2.0, label=label2)
1096 1094 ax.plotline3 = ax.plot(x[n][2,:,0,0], y, color='b',linewidth=2.0, label=label3)
1097 1095 ax.plotline4 = ax.plot(x[n][3,:,0,0], y, color='m',linewidth=2.0, label=label4)
1098 1096 ax.legend(loc='upper right')
1099 1097 ax.set_xlim(self.xmin, self.xmax)
1100 1098 self.titles.append('{}'.format(self.plot_name.upper()))
1101 1099
1102 1100 else:
1103 1101
1104 1102 if n==0:
1105 1103 self.xlimits=[(self.xmin,self.xmax)]
1106 1104 else:
1107 1105 self.xlimits.append((self.xmin,self.xmax))
1108 1106
1109 1107 ax.set_xlim(self.xmin, self.xmax)
1110 1108
1111 1109 ax.plotline1[0].set_data(x[n][0,:,0,0],y)
1112 1110 ax.plotline2[0].set_data(x[n][1,:,0,0],y)
1113 1111 ax.plotline3[0].set_data(x[n][2,:,0,0],y)
1114 1112 ax.plotline4[0].set_data(x[n][3,:,0,0],y)
1115 1113 self.titles.append('{}'.format(self.plot_name.upper()))
1116 1114
1117 1115
1118 1116 class CrossProductsLPPlot(Plot):
1119 1117 '''
1120 1118 Plot for cross products LP
1121 1119 '''
1122 1120
1123 1121 CODE = 'crossprodslp'
1124 1122 plot_name = 'Cross Products LP'
1125 1123 plot_type = 'scatterbuffer'
1126 1124
1127 1125
1128 1126 def setup(self):
1129 1127
1130 1128 self.ncols = 2
1131 1129 self.nrows = 1
1132 1130 self.nplots = 2
1133 1131 self.ylabel = 'Range [km]'
1134 1132 self.xlabel = 'dB'
1135 1133 self.width = 3.5*self.nplots
1136 1134 self.height = 5.5
1137 1135 self.colorbar = False
1138 1136 self.titles = []
1139 1137 self.plots_adjust.update({'wspace': .8 ,'left': 0.17, 'right': 0.88, 'bottom': 0.1})
1140 1138
1141 1139 def update(self, dataOut):
1142 1140 data = {}
1143 1141 meta = {}
1144 1142
1145 1143 data['crossprodslp'] = 10*numpy.log10(numpy.abs(dataOut.output_LP))
1146 1144
1147 1145 data['NRANGE'] = dataOut.NRANGE #This is metadata
1148 1146 data['NLAG'] = dataOut.NLAG #This is metadata
1149 1147
1150 1148 return data, meta
1151 1149
1152 1150 def plot(self):
1153 1151
1154 1152 NRANGE = self.data['NRANGE'][-1]
1155 1153 NLAG = self.data['NLAG'][-1]
1156 1154
1157 1155 x = self.data[self.CODE][:,-1,:,:]
1158 1156 self.y = self.data.yrange[0:NRANGE]
1159 1157
1160 1158 label_array=numpy.array(['lag '+ str(x) for x in range(NLAG)])
1161 1159 color_array=['r','k','g','b','c','m','y','orange','steelblue','purple','peru','darksalmon','grey','limegreen','olive','midnightblue']
1162 1160
1163 1161
1164 1162 for n, ax in enumerate(self.axes):
1165 1163
1166 1164 self.xmin=28#30
1167 1165 self.xmax=70#70
1168 1166 #self.xmin=numpy.min(numpy.concatenate((self.x[0,:,n],self.x[1,:,n])))
1169 1167 #self.xmax=numpy.max(numpy.concatenate((self.x[0,:,n],self.x[1,:,n])))
1170 1168
1171 1169 if ax.firsttime:
1172 1170
1173 1171 self.autoxticks=False
1174 1172 if n == 0:
1175 1173 self.plotline_array=numpy.zeros((2,NLAG),dtype=object)
1176 1174
1177 1175 for i in range(NLAG):
1178 1176 self.plotline_array[n,i], = ax.plot(x[i,:,n], self.y, color=color_array[i],linewidth=1.0, label=label_array[i])
1179 1177
1180 1178 ax.legend(loc='upper right')
1181 1179 ax.set_xlim(self.xmin, self.xmax)
1182 1180 if n==0:
1183 1181 self.titles.append('{} CH0'.format(self.plot_name.upper()))
1184 1182 if n==1:
1185 1183 self.titles.append('{} CH1'.format(self.plot_name.upper()))
1186 1184 else:
1187 1185 for i in range(NLAG):
1188 1186 self.plotline_array[n,i].set_data(x[i,:,n],self.y)
1189 1187
1190 1188 if n==0:
1191 1189 self.titles.append('{} CH0'.format(self.plot_name.upper()))
1192 1190 if n==1:
1193 1191 self.titles.append('{} CH1'.format(self.plot_name.upper()))
1194 1192
1195 1193
1196 1194 class NoiseDPPlot(NoisePlot):
1197 1195 '''
1198 1196 Plot for noise Double Pulse
1199 1197 '''
1200 1198
1201 1199 CODE = 'noise'
1202 1200 #plot_name = 'Noise'
1203 1201 #plot_type = 'scatterbuffer'
1204 1202
1205 1203 def update(self, dataOut):
1206 1204
1207 1205 data = {}
1208 1206 meta = {}
1209 1207 data['noise'] = 10*numpy.log10(dataOut.noise_final)
1210 1208
1211 1209 return data, meta
1212 1210
1213 1211
1214 1212 class XmitWaveformPlot(Plot):
1215 1213 '''
1216 1214 Plot for xmit waveform
1217 1215 '''
1218 1216
1219 1217 CODE = 'xmit'
1220 1218 plot_name = 'Xmit Waveform'
1221 1219 plot_type = 'scatterbuffer'
1222 1220
1223 1221
1224 1222 def setup(self):
1225 1223
1226 1224 self.ncols = 1
1227 1225 self.nrows = 1
1228 1226 self.nplots = 1
1229 1227 self.ylabel = ''
1230 1228 self.xlabel = 'Number of Lag'
1231 1229 self.width = 5.5
1232 1230 self.height = 3.5
1233 1231 self.colorbar = False
1234 1232 self.plots_adjust.update({'right': 0.85 })
1235 1233 self.titles = [self.plot_name]
1236 1234 #self.plots_adjust.update({'left': 0.17, 'right': 0.88, 'bottom': 0.1})
1237 1235
1238 1236 #if not self.titles:
1239 1237 #self.titles = self.data.parameters \
1240 1238 #if self.data.parameters else ['{}'.format(self.plot_name.upper())]
1241 1239
1242 1240 def update(self, dataOut):
1243 1241
1244 1242 data = {}
1245 1243 meta = {}
1246 1244
1247 1245 y_1=numpy.arctan2(dataOut.output_LP[:,0,2].imag,dataOut.output_LP[:,0,2].real)* 180 / (numpy.pi*10)
1248 1246 y_2=numpy.abs(dataOut.output_LP[:,0,2])
1249 1247 norm=numpy.max(y_2)
1250 1248 norm=max(norm,0.1)
1251 1249 y_2=y_2/norm
1252 1250
1253 1251 meta['yrange'] = numpy.array([])
1254 1252
1255 1253 data['xmit'] = numpy.vstack((y_1,y_2))
1256 1254 data['NLAG'] = dataOut.NLAG
1257 1255
1258 1256 return data, meta
1259 1257
1260 1258 def plot(self):
1261 1259
1262 1260 data = self.data[-1]
1263 1261 NLAG = data['NLAG']
1264 1262 x = numpy.arange(0,NLAG,1,'float32')
1265 1263 y = data['xmit']
1266 1264
1267 1265 self.xmin = 0
1268 1266 self.xmax = NLAG-1
1269 1267 self.ymin = -1.0
1270 1268 self.ymax = 1.0
1271 1269 ax = self.axes[0]
1272 1270
1273 1271 if ax.firsttime:
1274 1272 ax.plotline0=ax.plot(x,y[0,:],color='blue')
1275 1273 ax.plotline1=ax.plot(x,y[1,:],color='red')
1276 1274 secax=ax.secondary_xaxis(location=0.5)
1277 1275 secax.xaxis.tick_bottom()
1278 1276 secax.tick_params( labelleft=False, labeltop=False,
1279 1277 labelright=False, labelbottom=False)
1280 1278
1281 1279 self.xstep_given = 3
1282 1280 self.ystep_given = .25
1283 1281 secax.set_xticks(numpy.linspace(self.xmin, self.xmax, 6)) #only works on matplotlib.version>3.2
1284 1282
1285 1283 else:
1286 1284 ax.plotline0[0].set_data(x,y[0,:])
1287 1285 ax.plotline1[0].set_data(x,y[1,:])
Show file before | Show file after | Hide comments
@@ -1,247 +1,251
1 1 '''
2 2 Base clases to create Processing units and operations, the MPDecorator
3 3 must be used in plotting and writing operations to allow to run as an
4 4 external process.
5 5 '''
6 6
7 7 import os
8 8 import inspect
9 9 import zmq
10 10 import time
11 11 import pickle
12 12 import traceback
13 13 from threading import Thread
14 14 from multiprocessing import Process, Queue
15 15 from schainpy.utils import log
16 16
17 17 import copy
18 18
19 19 QUEUE_SIZE = int(os.environ.get('QUEUE_MAX_SIZE', '10'))
20 20
21 21 class ProcessingUnit(object):
22 22 '''
23 23 Base class to create Signal Chain Units
24 24 '''
25 25
26 26 proc_type = 'processing'
27 27
28 28 def __init__(self):
29 29
30 30 self.dataIn = None
31 31 self.dataOut = None
32 32 self.isConfig = False
33 33 self.operations = []
34 34 self.name = 'Test'
35 35 self.inputs = []
36 36
37 37 def setInput(self, unit):
38 38
39 39 attr = 'dataIn'
40 40 for i, u in enumerate(unit):
41 41 if i==0:
42 42 #print(u.dataOut.flagNoData)
43 43 #exit(1)
44 44 self.dataIn = u.dataOut#.copy()
45 45 self.inputs.append('dataIn')
46 46 else:
47 47 setattr(self, 'dataIn{}'.format(i), u.dataOut)#.copy())
48 48 self.inputs.append('dataIn{}'.format(i))
49 49
50 50
51 51 def getAllowedArgs(self):
52 52 if hasattr(self, '__attrs__'):
53 53 return self.__attrs__
54 54 else:
55 55 return inspect.getargspec(self.run).args
56 56
57 57 def addOperation(self, conf, operation):
58 58 '''
59 59 '''
60 60
61 61 self.operations.append((operation, conf.type, conf.getKwargs()))
62 62
63 63 def getOperationObj(self, objId):
64 64
65 65 if objId not in list(self.operations.keys()):
66 66 return None
67 67
68 68 return self.operations[objId]
69 69
70 70 def call(self, **kwargs):
71 71 '''
72 72 '''
73
73 #print("call")
74 74 try:
75 #print("try")
76 75 if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error:
77 #print("Try")
78 #print(self.dataIn)
79 #print(self.dataIn.flagNoData)
76 #if self.dataIn is not None and self.dataIn.flagNoData and not self.dataIn.error and not self.dataIn.runNextUnit:
77 if self.dataIn.runNextUnit:
78 #print("SUCCESSSSSSS")
79 #exit(1)
80 return not self.dataIn.isReady()
81 else:
80 82 return self.dataIn.isReady()
81 83 elif self.dataIn is None or not self.dataIn.error:
82 84 #print([getattr(self, at) for at in self.inputs])
83 85 #print("Elif 1")
84 86 self.run(**kwargs)
85 87 elif self.dataIn.error:
86 88 #print("Elif 2")
87 89 self.dataOut.error = self.dataIn.error
88 90 self.dataOut.flagNoData = True
89 91 except:
90 92 #print("Except")
91 93 err = traceback.format_exc()
92 94 if 'SchainWarning' in err:
93 95 log.warning(err.split('SchainWarning:')[-1].split('\n')[0].strip(), self.name)
94 96 elif 'SchainError' in err:
95 97 log.error(err.split('SchainError:')[-1].split('\n')[0].strip(), self.name)
96 98 else:
97 99 log.error(err, self.name)
98 100 self.dataOut.error = True
99 101 #print("before op")
100 102 for op, optype, opkwargs in self.operations:
101 103 aux = self.dataOut.copy()
102 104 #aux = copy.deepcopy(self.dataOut)
103 105 #print("**********************Before",op)
104 106 if optype == 'other' and not self.dataOut.flagNoData:
105 107 #print("**********************Other",op)
106 108 #print(self.dataOut.flagNoData)
107 109 self.dataOut = op.run(self.dataOut, **opkwargs)
108 110 elif optype == 'external' and not self.dataOut.flagNoData:
109 111 op.queue.put(aux)
110 112 elif optype == 'external' and self.dataOut.error:
111 113 op.queue.put(aux)
112 114 #elif optype == 'external' and self.dataOut.isReady():
113 115 #op.queue.put(copy.deepcopy(self.dataOut))
114 116 #print(not self.dataOut.isReady())
117
115 118 try:
116 119 if self.dataOut.runNextUnit:
117 120 runNextUnit = self.dataOut.runNextUnit
118 121 #print(self.operations)
119 122 #print("Tru")
120 123
121 124 else:
122 125 runNextUnit = self.dataOut.isReady()
123 126 except:
124 127 runNextUnit = self.dataOut.isReady()
125 128 #if not self.dataOut.isReady():
126 129 #return 'Error' if self.dataOut.error else input()
127 130 #print("NexT",runNextUnit)
131 #print("error: ",self.dataOut.error)
128 132 return 'Error' if self.dataOut.error else runNextUnit# self.dataOut.isReady()
129 133
130 134 def setup(self):
131 135
132 136 raise NotImplementedError
133 137
134 138 def run(self):
135 139
136 140 raise NotImplementedError
137 141
138 142 def close(self):
139 143
140 144 return
141 145
142 146
143 147 class Operation(object):
144 148
145 149 '''
146 150 '''
147 151
148 152 proc_type = 'operation'
149 153
150 154 def __init__(self):
151 155
152 156 self.id = None
153 157 self.isConfig = False
154 158
155 159 if not hasattr(self, 'name'):
156 160 self.name = self.__class__.__name__
157 161
158 162 def getAllowedArgs(self):
159 163 if hasattr(self, '__attrs__'):
160 164 return self.__attrs__
161 165 else:
162 166 return inspect.getargspec(self.run).args
163 167
164 168 def setup(self):
165 169
166 170 self.isConfig = True
167 171
168 172 raise NotImplementedError
169 173
170 174 def run(self, dataIn, **kwargs):
171 175 """
172 176 Realiza las operaciones necesarias sobre la dataIn.data y actualiza los
173 177 atributos del objeto dataIn.
174 178
175 179 Input:
176 180
177 181 dataIn : objeto del tipo JROData
178 182
179 183 Return:
180 184
181 185 None
182 186
183 187 Affected:
184 188 __buffer : buffer de recepcion de datos.
185 189
186 190 """
187 191 if not self.isConfig:
188 192 self.setup(**kwargs)
189 193
190 194 raise NotImplementedError
191 195
192 196 def close(self):
193 197
194 198 return
195 199
196 200
197 201 def MPDecorator(BaseClass):
198 202 """
199 203 Multiprocessing class decorator
200 204
201 205 This function add multiprocessing features to a BaseClass.
202 206 """
203 207
204 208 class MPClass(BaseClass, Process):
205 209
206 210 def __init__(self, *args, **kwargs):
207 211 super(MPClass, self).__init__()
208 212 Process.__init__(self)
209 213
210 214 self.args = args
211 215 self.kwargs = kwargs
212 216 self.t = time.time()
213 217 self.op_type = 'external'
214 218 self.name = BaseClass.__name__
215 219 self.__doc__ = BaseClass.__doc__
216 220
217 221 if 'plot' in self.name.lower() and not self.name.endswith('_'):
218 222 self.name = '{}{}'.format(self.CODE.upper(), 'Plot')
219 223
220 224 self.start_time = time.time()
221 225 self.err_queue = args[3]
222 226 self.queue = Queue(maxsize=QUEUE_SIZE)
223 227 self.myrun = BaseClass.run
224 228
225 229 def run(self):
226 230
227 231 while True:
228 232
229 233 dataOut = self.queue.get()
230 234
231 235 if not dataOut.error:
232 236 try:
233 237 BaseClass.run(self, dataOut, **self.kwargs)
234 238 except:
235 239 err = traceback.format_exc()
236 240 log.error(err, self.name)
237 241 else:
238 242 break
239 243
240 244 self.close()
241 245
242 246 def close(self):
243 247
244 248 BaseClass.close(self)
245 249 log.success('Done...(Time:{:4.2f} secs)'.format(time.time() - self.start_time), self.name)
246 250
247 251 return MPClass
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
This diff has been collapsed as it changes many lines, (859 lines changed) Show them Hide them
@@ -1,796 +1,1645
1 1 import numpy
2 2
3 3 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
4 4 from schainpy.model.data.jrodata import Spectra
5 5 from schainpy.model.data.jrodata import hildebrand_sekhon
6 6
7 class SpectraAFCProc(ProcessingUnit):
7 class SpectraAFCProc_V0(ProcessingUnit):
8 8
9 9 def __init__(self, **kwargs):
10 10
11 11 ProcessingUnit.__init__(self, **kwargs)
12 12
13 13 self.buffer = None
14 14 self.firstdatatime = None
15 15 self.profIndex = 0
16 16 self.dataOut = Spectra()
17 17 self.id_min = None
18 18 self.id_max = None
19 19
20 20 def __updateSpecFromVoltage(self):
21 21
22 22 self.dataOut.plotting = "spectra_acf"
23 23
24 24 self.dataOut.timeZone = self.dataIn.timeZone
25 25 self.dataOut.dstFlag = self.dataIn.dstFlag
26 26 self.dataOut.errorCount = self.dataIn.errorCount
27 27 self.dataOut.useLocalTime = self.dataIn.useLocalTime
28 28
29 29 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
30 30 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
31 31 self.dataOut.ippSeconds = self.dataIn.getDeltaH()*(10**-6)/0.15
32 32
33 33 self.dataOut.channelList = self.dataIn.channelList
34 34 self.dataOut.heightList = self.dataIn.heightList
35 35 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
36 36
37 37 self.dataOut.nBaud = self.dataIn.nBaud
38 38 self.dataOut.nCode = self.dataIn.nCode
39 39 self.dataOut.code = self.dataIn.code
40 40 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
41 41
42 42 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
43 43 self.dataOut.utctime = self.firstdatatime
44 44 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
45 45 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
46 46 self.dataOut.flagShiftFFT = False
47 47
48 48 self.dataOut.nCohInt = self.dataIn.nCohInt
49 49 self.dataOut.nIncohInt = 1
50 50
51 51 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
52 52
53 53 self.dataOut.frequency = self.dataIn.frequency
54 54 self.dataOut.realtime = self.dataIn.realtime
55 55
56 56 self.dataOut.azimuth = self.dataIn.azimuth
57 57 self.dataOut.zenith = self.dataIn.zenith
58 58
59 59 self.dataOut.beam.codeList = self.dataIn.beam.codeList
60 60 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
61 61 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
62 62
63 63 def __decodeData(self, nProfiles, code):
64 64
65 65 if code is None:
66 66 return
67 67
68 68 for i in range(nProfiles):
69 69 self.buffer[:,i,:] = self.buffer[:,i,:]*code[0][i]
70 70
71 71 def __getFft(self):
72 72 """
73 73 Convierte valores de Voltaje a Spectra
74 74
75 75 Affected:
76 76 self.dataOut.data_spc
77 77 self.dataOut.data_cspc
78 78 self.dataOut.data_dc
79 79 self.dataOut.heightList
80 80 self.profIndex
81 81 self.buffer
82 82 self.dataOut.flagNoData
83 83 """
84 84 nsegments = self.dataOut.nHeights
85 85
86 86 _fft_buffer = numpy.zeros((self.dataOut.nChannels, self.dataOut.nProfiles, nsegments), dtype='complex')
87 87
88 88 for i in range(nsegments):
89 89 try:
90 90 _fft_buffer[:,:,i] = self.buffer[:,i:i+self.dataOut.nProfiles]
91 91
92 92 if self.code is not None:
93 93 _fft_buffer[:,:,i] = _fft_buffer[:,:,i]*self.code[0]
94 94 except:
95 95 pass
96 96
97 97 fft_volt = numpy.fft.fft(_fft_buffer, n=self.dataOut.nFFTPoints, axis=1)
98 98 fft_volt = fft_volt.astype(numpy.dtype('complex'))
99 99 dc = fft_volt[:,0,:]
100 100
101 101 #calculo de self-spectra
102 102 spc = fft_volt * numpy.conjugate(fft_volt)
103 103 data = numpy.fft.ifft(spc, axis=1)
104 104 data = numpy.fft.fftshift(data, axes=(1,))
105 105
106 106 spc = data
107 107
108 108 blocksize = 0
109 109 blocksize += dc.size
110 110 blocksize += spc.size
111 111
112 112 cspc = None
113 113 pairIndex = 0
114 114
115 115 if self.dataOut.pairsList != None:
116 116 #calculo de cross-spectra
117 117 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
118 118 for pair in self.dataOut.pairsList:
119 119 if pair[0] not in self.dataOut.channelList:
120 120 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList)))
121 121 if pair[1] not in self.dataOut.channelList:
122 122 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList)))
123 123
124 124 chan_index0 = self.dataOut.channelList.index(pair[0])
125 125 chan_index1 = self.dataOut.channelList.index(pair[1])
126 126
127 127 cspc[pairIndex,:,:] = fft_volt[chan_index0,:,:] * numpy.conjugate(fft_volt[chan_index1,:,:])
128 128 pairIndex += 1
129 129 blocksize += cspc.size
130 130
131 131 self.dataOut.data_spc = spc
132 132 self.dataOut.data_cspc = cspc
133 133 self.dataOut.data_dc = dc
134 134 self.dataOut.blockSize = blocksize
135 135 self.dataOut.flagShiftFFT = True
136 136
137 137 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], code=None, nCode=1, nBaud=1,real= None, imag=None):
138 138
139 139 self.dataOut.flagNoData = True
140 140
141 141 if self.dataIn.type == "Spectra":
142 142 self.dataOut.copy(self.dataIn)
143 #print(self.dataOut.data.shape)
144 #exit(1)
143 145 spc = self.dataOut.data_spc
144 146 data = numpy.fft.fftshift( spc, axes=(1,))
145 147 data = numpy.fft.ifft(data, axis=1)
146 #data = numpy.fft.fftshift( data, axes=(1,))
147 #acf = numpy.abs(data) # Autocorrelacion LLAMAR A ESTE VALOR ACF
148 acf = data
148 data = numpy.fft.fftshift( data, axes=(1,))
149 acf = numpy.abs(data) # Autocorrelacion LLAMAR A ESTE VALOR ACF
150 acf = data #Comentarlo?
151 print("acf",acf[0,:,150])
152 exit(1)
149 153 #'''
150 154 if real:
151 155 acf = data.real
152 156 if imag:
153 157 acf = data.imag
154 158 #'''
155 159 shape = acf.shape # nchannels, nprofiles, nsamples //nchannles, lags , alturas
156 160
157 161 '''
158 162 for j in range(shape[0]):
159 163 for i in range(shape[2]):
160 164 tmp = int(shape[1]/2)
161 165 #print(i,j,tmp)
162 166 value = (acf[j,:,i][tmp-1]+acf[j,:,i][tmp+1])/2.0
163 167 acf[j,:,i][tmp] = value
164 168 # Normalizando
165 169 for i in range(shape[0]):
166 170 for j in range(shape[2]):
167 171 acf[i,:,j]= acf[i,:,j] / numpy.max(numpy.abs(acf[i,:,j]))
168 172 '''
169 173 self.dataOut.data_acf = acf
170 self.dataOut.data_spc = acf
174 self.dataOut.data_spc = acf.real
171 175 #print(self.dataOut.data_acf[0,:,0])
172 176 #exit(1)
173 177 '''
174 178 shape = self.dataOut.data_acf.shape
175 179 resFactor = 5
176 180 z = self.dataOut.data_acf.copy()
177 181 min = numpy.min(z[0,:,0])
178 182 max =numpy.max(z[0,:,0])
179 183 deltaHeight = self.dataOut.heightList[1]-self.dataOut.heightList[0]
180 184 for i in range(shape[0]):
181 185 for j in range(shape[2]):
182 186 z[i,:,j]= (((z[i,:,j]-min)/(max-min))*deltaHeight*resFactor + j*deltaHeight)
183 187 #print(self.dataOut.data_spc.shape)
184 188 #print(self.dataOut.data_acf.shape)
185 189 '''
190 '''
186 191 import matplotlib.pyplot as plt
187 192 hei = 10
188 193 #print(self.dataOut.heightList)
189 194 print(self.dataOut.heightList[hei])
190 195 #plt.plot(z[0,0,:],self.dataOut.heightList)
191 196 aux = self.dataOut.data_acf[0,0,:]
192 197 power = aux*numpy.conjugate(aux)
193 198 print(power)
194 199 powerdb = numpy.log10(power)
195 200 plt.plot(powerdb,self.dataOut.heightList)
196 201 #plt.plot(self.dataOut.data_acf[0,:,1])
197 202 plt.ylim(0,1000)
198 203 plt.show()
199 204 exit(1)
205 '''
206 return True
207
208 if code is not None:
209 self.code = numpy.array(code).reshape(nCode,nBaud)
210 else:
211 self.code = None
212
213 if self.dataIn.type == "Voltage":
214
215 if nFFTPoints == None:
216 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
217
218 if nProfiles == None:
219 nProfiles = nFFTPoints
220
221 self.dataOut.ippFactor = 1
222
223 self.dataOut.nFFTPoints = nFFTPoints
224 self.dataOut.nProfiles = nProfiles
225 self.dataOut.pairsList = pairsList
226
227 # if self.buffer is None:
228 # self.buffer = numpy.zeros( (self.dataIn.nChannels, nProfiles, self.dataIn.nHeights),
229 # dtype='complex')
230
231 if not self.dataIn.flagDataAsBlock:
232 self.buffer = self.dataIn.data.copy()
233
234 # for i in range(self.dataIn.nHeights):
235 # self.buffer[:, self.profIndex, self.profIndex:] = voltage_data[:,:self.dataIn.nHeights - self.profIndex]
236 #
237 # self.profIndex += 1
238
239 else:
240 raise ValueError("")
241
242 self.firstdatatime = self.dataIn.utctime
243
244 self.profIndex == nProfiles
245
246 self.__updateSpecFromVoltage()
247
248 self.__getFft()
249
250 self.dataOut.flagNoData = False
251
252 return True
253
254 raise ValueError("The type of input object '%s' is not valid"%(self.dataIn.type))
255
256 def __selectPairs(self, pairsList):
257
258 if channelList == None:
259 return
260
261 pairsIndexListSelected = []
262
263 for thisPair in pairsList:
264
265 if thisPair not in self.dataOut.pairsList:
266 continue
267
268 pairIndex = self.dataOut.pairsList.index(thisPair)
269
270 pairsIndexListSelected.append(pairIndex)
271
272 if not pairsIndexListSelected:
273 self.dataOut.data_cspc = None
274 self.dataOut.pairsList = []
275 return
276
277 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
278 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
279
280 return
281
282 def __selectPairsByChannel(self, channelList=None):
283
284 if channelList == None:
285 return
286
287 pairsIndexListSelected = []
288 for pairIndex in self.dataOut.pairsIndexList:
289 #First pair
290 if self.dataOut.pairsList[pairIndex][0] not in channelList:
291 continue
292 #Second pair
293 if self.dataOut.pairsList[pairIndex][1] not in channelList:
294 continue
295
296 pairsIndexListSelected.append(pairIndex)
297
298 if not pairsIndexListSelected:
299 self.dataOut.data_cspc = None
300 self.dataOut.pairsList = []
301 return
302
303 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
304 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
305
306 return
307
308 def selectChannels(self, channelList):
309
310 channelIndexList = []
311
312 for channel in channelList:
313 if channel not in self.dataOut.channelList:
314 raise ValueError("Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList)))
315
316 index = self.dataOut.channelList.index(channel)
317 channelIndexList.append(index)
318
319 self.selectChannelsByIndex(channelIndexList)
320
321 def selectChannelsByIndex(self, channelIndexList):
322 """
323 Selecciona un bloque de datos en base a canales segun el channelIndexList
324
325 Input:
326 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
200 327
328 Affected:
329 self.dataOut.data_spc
330 self.dataOut.channelIndexList
331 self.dataOut.nChannels
332
333 Return:
334 None
335 """
336
337 for channelIndex in channelIndexList:
338 if channelIndex not in self.dataOut.channelIndexList:
339 raise ValueError("Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList))
340
341 # nChannels = len(channelIndexList)
342
343 data_spc = self.dataOut.data_spc[channelIndexList,:]
344 data_dc = self.dataOut.data_dc[channelIndexList,:]
345
346 self.dataOut.data_spc = data_spc
347 self.dataOut.data_dc = data_dc
348
349 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
350 # self.dataOut.nChannels = nChannels
351
352 self.__selectPairsByChannel(self.dataOut.channelList)
353
354 return 1
355
356 def selectHeights(self, minHei, maxHei):
357 """
358 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
359 minHei <= height <= maxHei
360
361 Input:
362 minHei : valor minimo de altura a considerar
363 maxHei : valor maximo de altura a considerar
364
365 Affected:
366 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
367
368 Return:
369 1 si el metodo se ejecuto con exito caso contrario devuelve 0
370 """
371
372 if (minHei > maxHei):
373 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei))
374
375 if (minHei < self.dataOut.heightList[0]):
376 minHei = self.dataOut.heightList[0]
377
378 if (maxHei > self.dataOut.heightList[-1]):
379 maxHei = self.dataOut.heightList[-1]
380
381 minIndex = 0
382 maxIndex = 0
383 heights = self.dataOut.heightList
384
385 inda = numpy.where(heights >= minHei)
386 indb = numpy.where(heights <= maxHei)
387
388 try:
389 minIndex = inda[0][0]
390 except:
391 minIndex = 0
392
393 try:
394 maxIndex = indb[0][-1]
395 except:
396 maxIndex = len(heights)
397
398 self.selectHeightsByIndex(minIndex, maxIndex)
399
400 return 1
401
402 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
403 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
404
405 if hei_ref != None:
406 newheis = numpy.where(self.dataOut.heightList>hei_ref)
407
408 minIndex = min(newheis[0])
409 maxIndex = max(newheis[0])
410 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
411 heightList = self.dataOut.heightList[minIndex:maxIndex+1]
412
413 # determina indices
414 nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
415 avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
416 beacon_dB = numpy.sort(avg_dB)[-nheis:]
417 beacon_heiIndexList = []
418 for val in avg_dB.tolist():
419 if val >= beacon_dB[0]:
420 beacon_heiIndexList.append(avg_dB.tolist().index(val))
421
422 #data_spc = data_spc[:,:,beacon_heiIndexList]
423 data_cspc = None
424 if self.dataOut.data_cspc is not None:
425 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
426 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
427
428 data_dc = None
429 if self.dataOut.data_dc is not None:
430 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
431 #data_dc = data_dc[:,beacon_heiIndexList]
432
433 self.dataOut.data_spc = data_spc
434 self.dataOut.data_cspc = data_cspc
435 self.dataOut.data_dc = data_dc
436 self.dataOut.heightList = heightList
437 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
438
439 return 1
440
441
442 def selectHeightsByIndex(self, minIndex, maxIndex):
443 """
444 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
445 minIndex <= index <= maxIndex
446
447 Input:
448 minIndex : valor de indice minimo de altura a considerar
449 maxIndex : valor de indice maximo de altura a considerar
450
451 Affected:
452 self.dataOut.data_spc
453 self.dataOut.data_cspc
454 self.dataOut.data_dc
455 self.dataOut.heightList
456
457 Return:
458 1 si el metodo se ejecuto con exito caso contrario devuelve 0
459 """
460
461 if (minIndex < 0) or (minIndex > maxIndex):
462 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
463
464 if (maxIndex >= self.dataOut.nHeights):
465 maxIndex = self.dataOut.nHeights-1
466
467 #Spectra
468 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
469
470 data_cspc = None
471 if self.dataOut.data_cspc is not None:
472 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
473
474 data_dc = None
475 if self.dataOut.data_dc is not None:
476 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
477
478 self.dataOut.data_spc = data_spc
479 self.dataOut.data_cspc = data_cspc
480 self.dataOut.data_dc = data_dc
481
482 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
483
484 return 1
485
486 def removeDC(self, mode = 2):
487 jspectra = self.dataOut.data_spc
488 jcspectra = self.dataOut.data_cspc
489
490
491 num_chan = jspectra.shape[0]
492 num_hei = jspectra.shape[2]
493
494 if jcspectra is not None:
495 jcspectraExist = True
496 num_pairs = jcspectra.shape[0]
497 else: jcspectraExist = False
498
499 freq_dc = jspectra.shape[1]/2
500 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
501
502 if ind_vel[0]<0:
503 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
504
505 if mode == 1:
506 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
507
508 if jcspectraExist:
509 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
510
511 if mode == 2:
512
513 vel = numpy.array([-2,-1,1,2])
514 xx = numpy.zeros([4,4])
515
516 for fil in range(4):
517 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
518
519 xx_inv = numpy.linalg.inv(xx)
520 xx_aux = xx_inv[0,:]
521
522 for ich in range(num_chan):
523 yy = jspectra[ich,ind_vel,:]
524 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
525
526 junkid = jspectra[ich,freq_dc,:]<=0
527 cjunkid = sum(junkid)
528
529 if cjunkid.any():
530 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
531
532 if jcspectraExist:
533 for ip in range(num_pairs):
534 yy = jcspectra[ip,ind_vel,:]
535 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
536
537
538 self.dataOut.data_spc = jspectra
539 self.dataOut.data_cspc = jcspectra
540
541 return 1
542
543 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
544
545 jspectra = self.dataOut.data_spc
546 jcspectra = self.dataOut.data_cspc
547 jnoise = self.dataOut.getNoise()
548 num_incoh = self.dataOut.nIncohInt
549
550 num_channel = jspectra.shape[0]
551 num_prof = jspectra.shape[1]
552 num_hei = jspectra.shape[2]
553
554 #hei_interf
555 if hei_interf is None:
556 count_hei = num_hei/2 #Como es entero no importa
557 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
558 hei_interf = numpy.asarray(hei_interf)[0]
559 #nhei_interf
560 if (nhei_interf == None):
561 nhei_interf = 5
562 if (nhei_interf < 1):
563 nhei_interf = 1
564 if (nhei_interf > count_hei):
565 nhei_interf = count_hei
566 if (offhei_interf == None):
567 offhei_interf = 0
568
569 ind_hei = range(num_hei)
570 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
571 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
572 mask_prof = numpy.asarray(range(num_prof))
573 num_mask_prof = mask_prof.size
574 comp_mask_prof = [0, num_prof/2]
575
576
577 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
578 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
579 jnoise = numpy.nan
580 noise_exist = jnoise[0] < numpy.Inf
581
582 #Subrutina de Remocion de la Interferencia
583 for ich in range(num_channel):
584 #Se ordena los espectros segun su potencia (menor a mayor)
585 power = jspectra[ich,mask_prof,:]
586 power = power[:,hei_interf]
587 power = power.sum(axis = 0)
588 psort = power.ravel().argsort()
589
590 #Se estima la interferencia promedio en los Espectros de Potencia empleando
591 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
592
593 if noise_exist:
594 # tmp_noise = jnoise[ich] / num_prof
595 tmp_noise = jnoise[ich]
596 junkspc_interf = junkspc_interf - tmp_noise
597 #junkspc_interf[:,comp_mask_prof] = 0
598
599 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
600 jspc_interf = jspc_interf.transpose()
601 #Calculando el espectro de interferencia promedio
602 noiseid = numpy.where(jspc_interf <= tmp_noise/ numpy.sqrt(num_incoh))
603 noiseid = noiseid[0]
604 cnoiseid = noiseid.size
605 interfid = numpy.where(jspc_interf > tmp_noise/ numpy.sqrt(num_incoh))
606 interfid = interfid[0]
607 cinterfid = interfid.size
608
609 if (cnoiseid > 0): jspc_interf[noiseid] = 0
610
611 #Expandiendo los perfiles a limpiar
612 if (cinterfid > 0):
613 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
614 new_interfid = numpy.asarray(new_interfid)
615 new_interfid = {x for x in new_interfid}
616 new_interfid = numpy.array(list(new_interfid))
617 new_cinterfid = new_interfid.size
618 else: new_cinterfid = 0
619
620 for ip in range(new_cinterfid):
621 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
622 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
623
624
625 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
626
627 #Removiendo la interferencia del punto de mayor interferencia
628 ListAux = jspc_interf[mask_prof].tolist()
629 maxid = ListAux.index(max(ListAux))
630
631
632 if cinterfid > 0:
633 for ip in range(cinterfid*(interf == 2) - 1):
634 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/numpy.sqrt(num_incoh))).nonzero()
635 cind = len(ind)
636
637 if (cind > 0):
638 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/numpy.sqrt(num_incoh))
639
640 ind = numpy.array([-2,-1,1,2])
641 xx = numpy.zeros([4,4])
642
643 for id1 in range(4):
644 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
645
646 xx_inv = numpy.linalg.inv(xx)
647 xx = xx_inv[:,0]
648 ind = (ind + maxid + num_mask_prof)%num_mask_prof
649 yy = jspectra[ich,mask_prof[ind],:]
650 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
651
652
653 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/numpy.sqrt(num_incoh))).nonzero()
654 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/numpy.sqrt(num_incoh))
655
656 #Remocion de Interferencia en el Cross Spectra
657 if jcspectra is None: return jspectra, jcspectra
658 num_pairs = jcspectra.size/(num_prof*num_hei)
659 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
660
661 for ip in range(num_pairs):
662
663 #-------------------------------------------
664
665 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
666 cspower = cspower[:,hei_interf]
667 cspower = cspower.sum(axis = 0)
668
669 cspsort = cspower.ravel().argsort()
670 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
671 junkcspc_interf = junkcspc_interf.transpose()
672 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
673
674 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
675
676 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
677 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
678 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
679
680 for iprof in range(num_prof):
681 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
682 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
683
684 #Removiendo la Interferencia
685 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
686
687 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
688 maxid = ListAux.index(max(ListAux))
689
690 ind = numpy.array([-2,-1,1,2])
691 xx = numpy.zeros([4,4])
692
693 for id1 in range(4):
694 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
695
696 xx_inv = numpy.linalg.inv(xx)
697 xx = xx_inv[:,0]
698
699 ind = (ind + maxid + num_mask_prof)%num_mask_prof
700 yy = jcspectra[ip,mask_prof[ind],:]
701 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
702
703 #Guardar Resultados
704 self.dataOut.data_spc = jspectra
705 self.dataOut.data_cspc = jcspectra
706
707 return 1
708
709 def setRadarFrequency(self, frequency=None):
710
711 if frequency != None:
712 self.dataOut.frequency = frequency
713
714 return 1
715
716 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
717 #validacion de rango
718 if minHei == None:
719 minHei = self.dataOut.heightList[0]
720
721 if maxHei == None:
722 maxHei = self.dataOut.heightList[-1]
723
724 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
725 print('minHei: %.2f is out of the heights range'%(minHei))
726 print('minHei is setting to %.2f'%(self.dataOut.heightList[0]))
727 minHei = self.dataOut.heightList[0]
728
729 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
730 print('maxHei: %.2f is out of the heights range'%(maxHei))
731 print('maxHei is setting to %.2f'%(self.dataOut.heightList[-1]))
732 maxHei = self.dataOut.heightList[-1]
733
734 # validacion de velocidades
735 velrange = self.dataOut.getVelRange(1)
736
737 if minVel == None:
738 minVel = velrange[0]
739
740 if maxVel == None:
741 maxVel = velrange[-1]
742
743 if (minVel < velrange[0]) or (minVel > maxVel):
744 print('minVel: %.2f is out of the velocity range'%(minVel))
745 print('minVel is setting to %.2f'%(velrange[0]))
746 minVel = velrange[0]
747
748 if (maxVel > velrange[-1]) or (maxVel < minVel):
749 print('maxVel: %.2f is out of the velocity range'%(maxVel))
750 print('maxVel is setting to %.2f'%(velrange[-1]))
751 maxVel = velrange[-1]
752
753 # seleccion de indices para rango
754 minIndex = 0
755 maxIndex = 0
756 heights = self.dataOut.heightList
757
758 inda = numpy.where(heights >= minHei)
759 indb = numpy.where(heights <= maxHei)
760
761 try:
762 minIndex = inda[0][0]
763 except:
764 minIndex = 0
765
766 try:
767 maxIndex = indb[0][-1]
768 except:
769 maxIndex = len(heights)
770
771 if (minIndex < 0) or (minIndex > maxIndex):
772 raise ValueError("some value in (%d,%d) is not valid" % (minIndex, maxIndex))
773
774 if (maxIndex >= self.dataOut.nHeights):
775 maxIndex = self.dataOut.nHeights-1
776
777 # seleccion de indices para velocidades
778 indminvel = numpy.where(velrange >= minVel)
779 indmaxvel = numpy.where(velrange <= maxVel)
780 try:
781 minIndexVel = indminvel[0][0]
782 except:
783 minIndexVel = 0
784
785 try:
786 maxIndexVel = indmaxvel[0][-1]
787 except:
788 maxIndexVel = len(velrange)
789
790 #seleccion del espectro
791 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
792 #estimacion de ruido
793 noise = numpy.zeros(self.dataOut.nChannels)
794
795 for channel in range(self.dataOut.nChannels):
796 daux = data_spc[channel,:,:]
797 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
798
799 self.dataOut.noise_estimation = noise.copy()
800
801 return 1
802
803 class SpectraAFCProc(ProcessingUnit):
804
805 def __init__(self, **kwargs):
806
807 ProcessingUnit.__init__(self, **kwargs)
808
809 self.buffer = None
810 self.firstdatatime = None
811 self.profIndex = 0
812 self.dataOut = Spectra()
813 self.id_min = None
814 self.id_max = None
815
816 def __updateSpecFromVoltage(self):
817
818 self.dataOut.plotting = "spectra_acf"
819
820 self.dataOut.timeZone = self.dataIn.timeZone
821 self.dataOut.dstFlag = self.dataIn.dstFlag
822 self.dataOut.errorCount = self.dataIn.errorCount
823 self.dataOut.useLocalTime = self.dataIn.useLocalTime
824
825 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
826 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
827 self.dataOut.ippSeconds = self.dataIn.getDeltaH()*(10**-6)/0.15
828
829 self.dataOut.channelList = self.dataIn.channelList
830 self.dataOut.heightList = self.dataIn.heightList
831 self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
832
833 self.dataOut.nBaud = self.dataIn.nBaud
834 self.dataOut.nCode = self.dataIn.nCode
835 self.dataOut.code = self.dataIn.code
836 # self.dataOut.nProfiles = self.dataOut.nFFTPoints
837
838 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
839 self.dataOut.utctime = self.firstdatatime
840 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
841 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
842 self.dataOut.flagShiftFFT = False
843
844 self.dataOut.nCohInt = self.dataIn.nCohInt
845 self.dataOut.nIncohInt = 1
846
847 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
848
849 self.dataOut.frequency = self.dataIn.frequency
850 self.dataOut.realtime = self.dataIn.realtime
851
852 self.dataOut.azimuth = self.dataIn.azimuth
853 self.dataOut.zenith = self.dataIn.zenith
854
855 self.dataOut.beam.codeList = self.dataIn.beam.codeList
856 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
857 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
858
859 def __decodeData(self, nProfiles, code):
860
861 if code is None:
862 return
863
864 for i in range(nProfiles):
865 self.buffer[:,i,:] = self.buffer[:,i,:]*code[0][i]
866
867 def __getFft(self):
868 """
869 Convierte valores de Voltaje a Spectra
870
871 Affected:
872 self.dataOut.data_spc
873 self.dataOut.data_cspc
874 self.dataOut.data_dc
875 self.dataOut.heightList
876 self.profIndex
877 self.buffer
878 self.dataOut.flagNoData
879 """
880 nsegments = self.dataOut.nHeights
881
882 _fft_buffer = numpy.zeros((self.dataOut.nChannels, self.dataOut.nProfiles, nsegments), dtype='complex')
883
884 for i in range(nsegments):
885 try:
886 _fft_buffer[:,:,i] = self.buffer[:,i:i+self.dataOut.nProfiles]
887
888 if self.code is not None:
889 _fft_buffer[:,:,i] = _fft_buffer[:,:,i]*self.code[0]
890 except:
891 pass
892
893 fft_volt = numpy.fft.fft(_fft_buffer, n=self.dataOut.nFFTPoints, axis=1)
894 fft_volt = fft_volt.astype(numpy.dtype('complex'))
895 dc = fft_volt[:,0,:]
896
897 #calculo de self-spectra
898 spc = fft_volt * numpy.conjugate(fft_volt)
899 data = numpy.fft.ifft(spc, axis=1)
900 data = numpy.fft.fftshift(data, axes=(1,))
901
902 spc = data
903
904 blocksize = 0
905 blocksize += dc.size
906 blocksize += spc.size
907
908 cspc = None
909 pairIndex = 0
910
911 if self.dataOut.pairsList != None:
912 #calculo de cross-spectra
913 cspc = numpy.zeros((self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
914 for pair in self.dataOut.pairsList:
915 if pair[0] not in self.dataOut.channelList:
916 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList)))
917 if pair[1] not in self.dataOut.channelList:
918 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" %(str(pair), str(self.dataOut.channelList)))
919
920 chan_index0 = self.dataOut.channelList.index(pair[0])
921 chan_index1 = self.dataOut.channelList.index(pair[1])
922
923 cspc[pairIndex,:,:] = fft_volt[chan_index0,:,:] * numpy.conjugate(fft_volt[chan_index1,:,:])
924 pairIndex += 1
925 blocksize += cspc.size
926
927 self.dataOut.data_spc = spc
928 self.dataOut.data_cspc = cspc
929 self.dataOut.data_dc = dc
930 self.dataOut.blockSize = blocksize
931 self.dataOut.flagShiftFFT = True
932
933 def run(self, nProfiles=None, nFFTPoints=None, pairsList=[], code=None, nCode=1, nBaud=1):
934
935 #self.dataOut.flagNoData = True
936
937 #self.dataIn.runNextUnit = runNextUnit
938 #print("here")
939 if self.dataIn.type == "Spectra":
940 self.dataOut.copy(self.dataIn)
941
942 spc = self.dataOut.data_spc
943 data = numpy.fft.fftshift( spc, axes=(1,))
944 data = numpy.fft.ifft(data, axis=1)
945 #data = numpy.fft.ifft(data, axis=1, n = 32)
946 #data = numpy.fft.fftshift( data, axes=(1,))
947 #acf = numpy.abs(data)
948 acf = data[:,:16,:]
949 #acf = data[:,16:,:]
950 #print("SUM: ",numpy.sum(acf))
951 #print(acf.shape)
952 '''
953 hei_id = 35
954
955 aux2 = numpy.fft.fft(self.dataOut.data[0,:,hei_id],n = 16)
956 aux2 = numpy.fft.fftshift(aux2)
957 aux2 = aux2*numpy.conjugate(aux2)
958 aux2 = aux2.real #Este valor es el que da SCh
959 print("spc 2: ",numpy.sum(aux2))
960 aux2 = numpy.fft.fftshift(aux2)
961 aux2 = numpy.fft.ifft(aux2)
962 print("Rate_Right?: ",aux2[0]/corr[0,0,hei_id])
963
964 print("AFC sum: ",numpy.sum(acf[0,:,hei_id]))
965 print("aux2 sum: ",numpy.sum(aux2))
966 print("Rate: ",acf[0,0,hei_id]/corr[0,0,hei_id])
967 print("Rate aux: ",acf[0,0,hei_id]/corr_aux[0,-1,hei_id])
968 '''
969 #print(acf[0,:,150])
970 #exit(1)
971 '''
972 if real:
973 acf = data.real
974 if imag:
975 acf = data.imag
976 '''
977 #shape = acf.shape # nchannels, nprofiles, nsamples //nchannles, lags , alturas
978 '''
979 for j in range(shape[0]):
980 for i in range(shape[2]):
981 tmp = int(shape[1]/2)
982 #print(i,j,tmp)
983 value = (acf[j,:,i][tmp-1]+acf[j,:,i][tmp+1])/2.0
984 acf[j,:,i][tmp] = value
985 '''
986 #acf = numpy.fft.fftshift( acf, axes=(1,))
987 '''
988 # Normalizando
989 for i in range(shape[0]):
990 for j in range(shape[2]):
991 acf[i,:,j]= acf[i,:,j] / numpy.max(numpy.abs(acf[i,:,j]))
992 '''
993
994 #self.dataOut.data_acf = acf[:,:16,:]#*2
995 #self.dataOut.data_spc = acf[:,:16,:].real#*2
996
997 self.dataOut.data_acf = acf
998 '''
999 self.dataOut.data_spc = data.imag
1000
1001 print("Real: ",data[0,:,26].real)
1002 print("Real dB:",10*numpy.log10(data[0,:,26].real))
1003 print("Imag: ",data[0,:,26].imag)
1004 print("Imag dB:",10*numpy.log10(data[0,:,26].imag))
1005 exit(1)
1006 '''
1007 #print("AFC",self.dataOut.flagNoData)
1008 #'''
1009 #print("acf 0: ", self.dataOut.data_acf[0,0,100])
1010 #print("spc: ",numpy.mean(self.dataOut.data_spc[0,:,100]))
1011 #print("spc 0: ",numpy.fft.fftshift(self.dataOut.data_spc[0,:,100])[0])
1012 #exit(1)
1013 #self.dataOut.data_spc = acf.real
1014 #print(self.dataOut.data_acf[0,:,0])
1015 #exit(1)
1016 #'''
1017 '''
1018 shape = self.dataOut.data_acf.shape
1019 resFactor = 5
1020 z = self.dataOut.data_acf.copy()
1021 min = numpy.min(z[0,:,0])
1022 max =numpy.max(z[0,:,0])
1023 deltaHeight = self.dataOut.heightList[1]-self.dataOut.heightList[0]
1024 for i in range(shape[0]):
1025 for j in range(shape[2]):
1026 z[i,:,j]= (((z[i,:,j]-min)/(max-min))*deltaHeight*resFactor + j*deltaHeight)
1027 #print(self.dataOut.data_spc.shape)
1028 #print(self.dataOut.data_acf.shape)
1029 '''
1030 '''
1031 import matplotlib.pyplot as plt
1032 #hei = 10
1033 #print(self.dataOut.heightList)
1034 #print(self.dataOut.heightList[hei])
1035 #plt.plot(z[0,0,:],self.dataOut.heightList)
1036 aux = self.dataOut.data_acf[0,:,100]
1037 power = aux*numpy.conjugate(aux)
1038 #print(power)
1039 powerdb = numpy.log10(power)
1040 #plt.plot(powerdb,self.dataOut.heightList)
1041 #plt.plot(self.dataOut.data_acf[0,:,33])
1042 #plt.plot(corr)
1043 #plt.imshow(corr.real[0])
1044 plt.imshow(self.dataOut.data_acf.real[0])
1045 #plt.ylim(0,1000)
1046 plt.grid()
1047 plt.show()
1048 exit(1)
1049 '''
201 1050 return True
202 1051
203 1052 if code is not None:
204 1053 self.code = numpy.array(code).reshape(nCode,nBaud)
205 1054 else:
206 1055 self.code = None
207 1056
208 1057 if self.dataIn.type == "Voltage":
209 1058
210 1059 if nFFTPoints == None:
211 1060 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
212 1061
213 1062 if nProfiles == None:
214 1063 nProfiles = nFFTPoints
215 1064
216 1065 self.dataOut.ippFactor = 1
217 1066
218 1067 self.dataOut.nFFTPoints = nFFTPoints
219 1068 self.dataOut.nProfiles = nProfiles
220 1069 self.dataOut.pairsList = pairsList
221 1070
222 1071 # if self.buffer is None:
223 1072 # self.buffer = numpy.zeros( (self.dataIn.nChannels, nProfiles, self.dataIn.nHeights),
224 1073 # dtype='complex')
225 1074
226 1075 if not self.dataIn.flagDataAsBlock:
227 1076 self.buffer = self.dataIn.data.copy()
228 1077
229 1078 # for i in range(self.dataIn.nHeights):
230 1079 # self.buffer[:, self.profIndex, self.profIndex:] = voltage_data[:,:self.dataIn.nHeights - self.profIndex]
231 1080 #
232 1081 # self.profIndex += 1
233 1082
234 1083 else:
235 1084 raise ValueError("")
236 1085
237 1086 self.firstdatatime = self.dataIn.utctime
238 1087
239 1088 self.profIndex == nProfiles
240 1089
241 1090 self.__updateSpecFromVoltage()
242 1091
243 1092 self.__getFft()
244 1093
245 1094 self.dataOut.flagNoData = False
246 1095
247 1096 return True
248 1097
249 1098 raise ValueError("The type of input object '%s' is not valid"%(self.dataIn.type))
250 1099
251 1100 def __selectPairs(self, pairsList):
252 1101
253 1102 if channelList == None:
254 1103 return
255 1104
256 1105 pairsIndexListSelected = []
257 1106
258 1107 for thisPair in pairsList:
259 1108
260 1109 if thisPair not in self.dataOut.pairsList:
261 1110 continue
262 1111
263 1112 pairIndex = self.dataOut.pairsList.index(thisPair)
264 1113
265 1114 pairsIndexListSelected.append(pairIndex)
266 1115
267 1116 if not pairsIndexListSelected:
268 1117 self.dataOut.data_cspc = None
269 1118 self.dataOut.pairsList = []
270 1119 return
271 1120
272 1121 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
273 1122 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
274 1123
275 1124 return
276 1125
277 1126 def __selectPairsByChannel(self, channelList=None):
278 1127
279 1128 if channelList == None:
280 1129 return
281 1130
282 1131 pairsIndexListSelected = []
283 1132 for pairIndex in self.dataOut.pairsIndexList:
284 1133 #First pair
285 1134 if self.dataOut.pairsList[pairIndex][0] not in channelList:
286 1135 continue
287 1136 #Second pair
288 1137 if self.dataOut.pairsList[pairIndex][1] not in channelList:
289 1138 continue
290 1139
291 1140 pairsIndexListSelected.append(pairIndex)
292 1141
293 1142 if not pairsIndexListSelected:
294 1143 self.dataOut.data_cspc = None
295 1144 self.dataOut.pairsList = []
296 1145 return
297 1146
298 1147 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
299 1148 self.dataOut.pairsList = [self.dataOut.pairsList[i] for i in pairsIndexListSelected]
300 1149
301 1150 return
302 1151
303 1152 def selectChannels(self, channelList):
304 1153
305 1154 channelIndexList = []
306 1155
307 1156 for channel in channelList:
308 1157 if channel not in self.dataOut.channelList:
309 1158 raise ValueError("Error selecting channels, Channel %d is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList)))
310 1159
311 1160 index = self.dataOut.channelList.index(channel)
312 1161 channelIndexList.append(index)
313 1162
314 1163 self.selectChannelsByIndex(channelIndexList)
315 1164
316 1165 def selectChannelsByIndex(self, channelIndexList):
317 1166 """
318 1167 Selecciona un bloque de datos en base a canales segun el channelIndexList
319 1168
320 1169 Input:
321 1170 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
322 1171
323 1172 Affected:
324 1173 self.dataOut.data_spc
325 1174 self.dataOut.channelIndexList
326 1175 self.dataOut.nChannels
327 1176
328 1177 Return:
329 1178 None
330 1179 """
331 1180
332 1181 for channelIndex in channelIndexList:
333 1182 if channelIndex not in self.dataOut.channelIndexList:
334 1183 raise ValueError("Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList))
335 1184
336 1185 # nChannels = len(channelIndexList)
337 1186
338 1187 data_spc = self.dataOut.data_spc[channelIndexList,:]
339 1188 data_dc = self.dataOut.data_dc[channelIndexList,:]
340 1189
341 1190 self.dataOut.data_spc = data_spc
342 1191 self.dataOut.data_dc = data_dc
343 1192
344 1193 self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
345 1194 # self.dataOut.nChannels = nChannels
346 1195
347 1196 self.__selectPairsByChannel(self.dataOut.channelList)
348 1197
349 1198 return 1
350 1199
351 1200 def selectHeights(self, minHei, maxHei):
352 1201 """
353 1202 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
354 1203 minHei <= height <= maxHei
355 1204
356 1205 Input:
357 1206 minHei : valor minimo de altura a considerar
358 1207 maxHei : valor maximo de altura a considerar
359 1208
360 1209 Affected:
361 1210 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
362 1211
363 1212 Return:
364 1213 1 si el metodo se ejecuto con exito caso contrario devuelve 0
365 1214 """
366 1215
367 1216 if (minHei > maxHei):
368 1217 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei))
369 1218
370 1219 if (minHei < self.dataOut.heightList[0]):
371 1220 minHei = self.dataOut.heightList[0]
372 1221
373 1222 if (maxHei > self.dataOut.heightList[-1]):
374 1223 maxHei = self.dataOut.heightList[-1]
375 1224
376 1225 minIndex = 0
377 1226 maxIndex = 0
378 1227 heights = self.dataOut.heightList
379 1228
380 1229 inda = numpy.where(heights >= minHei)
381 1230 indb = numpy.where(heights <= maxHei)
382 1231
383 1232 try:
384 1233 minIndex = inda[0][0]
385 1234 except:
386 1235 minIndex = 0
387 1236
388 1237 try:
389 1238 maxIndex = indb[0][-1]
390 1239 except:
391 1240 maxIndex = len(heights)
392 1241
393 1242 self.selectHeightsByIndex(minIndex, maxIndex)
394 1243
395 1244 return 1
396 1245
397 1246 def getBeaconSignal(self, tauindex = 0, channelindex = 0, hei_ref=None):
398 1247 newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
399 1248
400 1249 if hei_ref != None:
401 1250 newheis = numpy.where(self.dataOut.heightList>hei_ref)
402 1251
403 1252 minIndex = min(newheis[0])
404 1253 maxIndex = max(newheis[0])
405 1254 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
406 1255 heightList = self.dataOut.heightList[minIndex:maxIndex+1]
407 1256
408 1257 # determina indices
409 1258 nheis = int(self.dataOut.radarControllerHeaderObj.txB/(self.dataOut.heightList[1]-self.dataOut.heightList[0]))
410 1259 avg_dB = 10*numpy.log10(numpy.sum(data_spc[channelindex,:,:],axis=0))
411 1260 beacon_dB = numpy.sort(avg_dB)[-nheis:]
412 1261 beacon_heiIndexList = []
413 1262 for val in avg_dB.tolist():
414 1263 if val >= beacon_dB[0]:
415 1264 beacon_heiIndexList.append(avg_dB.tolist().index(val))
416 1265
417 1266 #data_spc = data_spc[:,:,beacon_heiIndexList]
418 1267 data_cspc = None
419 1268 if self.dataOut.data_cspc is not None:
420 1269 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
421 1270 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
422 1271
423 1272 data_dc = None
424 1273 if self.dataOut.data_dc is not None:
425 1274 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
426 1275 #data_dc = data_dc[:,beacon_heiIndexList]
427 1276
428 1277 self.dataOut.data_spc = data_spc
429 1278 self.dataOut.data_cspc = data_cspc
430 1279 self.dataOut.data_dc = data_dc
431 1280 self.dataOut.heightList = heightList
432 1281 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
433 1282
434 1283 return 1
435 1284
436 1285
437 1286 def selectHeightsByIndex(self, minIndex, maxIndex):
438 1287 """
439 1288 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
440 1289 minIndex <= index <= maxIndex
441 1290
442 1291 Input:
443 1292 minIndex : valor de indice minimo de altura a considerar
444 1293 maxIndex : valor de indice maximo de altura a considerar
445 1294
446 1295 Affected:
447 1296 self.dataOut.data_spc
448 1297 self.dataOut.data_cspc
449 1298 self.dataOut.data_dc
450 1299 self.dataOut.heightList
451 1300
452 1301 Return:
453 1302 1 si el metodo se ejecuto con exito caso contrario devuelve 0
454 1303 """
455 1304
456 1305 if (minIndex < 0) or (minIndex > maxIndex):
457 1306 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
458 1307
459 1308 if (maxIndex >= self.dataOut.nHeights):
460 1309 maxIndex = self.dataOut.nHeights-1
461 1310
462 1311 #Spectra
463 1312 data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
464 1313
465 1314 data_cspc = None
466 1315 if self.dataOut.data_cspc is not None:
467 1316 data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
468 1317
469 1318 data_dc = None
470 1319 if self.dataOut.data_dc is not None:
471 1320 data_dc = self.dataOut.data_dc[:,minIndex:maxIndex+1]
472 1321
473 1322 self.dataOut.data_spc = data_spc
474 1323 self.dataOut.data_cspc = data_cspc
475 1324 self.dataOut.data_dc = data_dc
476 1325
477 1326 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1]
478 1327
479 1328 return 1
480 1329
481 1330 def removeDC(self, mode = 2):
482 1331 jspectra = self.dataOut.data_spc
483 1332 jcspectra = self.dataOut.data_cspc
484 1333
485 1334
486 1335 num_chan = jspectra.shape[0]
487 1336 num_hei = jspectra.shape[2]
488 1337
489 1338 if jcspectra is not None:
490 1339 jcspectraExist = True
491 1340 num_pairs = jcspectra.shape[0]
492 1341 else: jcspectraExist = False
493 1342
494 1343 freq_dc = jspectra.shape[1]/2
495 1344 ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
496 1345
497 1346 if ind_vel[0]<0:
498 1347 ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
499 1348
500 1349 if mode == 1:
501 1350 jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
502 1351
503 1352 if jcspectraExist:
504 1353 jcspectra[:,freq_dc,:] = (jcspectra[:,ind_vel[1],:] + jcspectra[:,ind_vel[2],:])/2
505 1354
506 1355 if mode == 2:
507 1356
508 1357 vel = numpy.array([-2,-1,1,2])
509 1358 xx = numpy.zeros([4,4])
510 1359
511 1360 for fil in range(4):
512 1361 xx[fil,:] = vel[fil]**numpy.asarray(range(4))
513 1362
514 1363 xx_inv = numpy.linalg.inv(xx)
515 1364 xx_aux = xx_inv[0,:]
516 1365
517 1366 for ich in range(num_chan):
518 1367 yy = jspectra[ich,ind_vel,:]
519 1368 jspectra[ich,freq_dc,:] = numpy.dot(xx_aux,yy)
520 1369
521 1370 junkid = jspectra[ich,freq_dc,:]<=0
522 1371 cjunkid = sum(junkid)
523 1372
524 1373 if cjunkid.any():
525 1374 jspectra[ich,freq_dc,junkid.nonzero()] = (jspectra[ich,ind_vel[1],junkid] + jspectra[ich,ind_vel[2],junkid])/2
526 1375
527 1376 if jcspectraExist:
528 1377 for ip in range(num_pairs):
529 1378 yy = jcspectra[ip,ind_vel,:]
530 1379 jcspectra[ip,freq_dc,:] = numpy.dot(xx_aux,yy)
531 1380
532 1381
533 1382 self.dataOut.data_spc = jspectra
534 1383 self.dataOut.data_cspc = jcspectra
535 1384
536 1385 return 1
537 1386
538 1387 def removeInterference(self, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None):
539 1388
540 1389 jspectra = self.dataOut.data_spc
541 1390 jcspectra = self.dataOut.data_cspc
542 1391 jnoise = self.dataOut.getNoise()
543 1392 num_incoh = self.dataOut.nIncohInt
544 1393
545 1394 num_channel = jspectra.shape[0]
546 1395 num_prof = jspectra.shape[1]
547 1396 num_hei = jspectra.shape[2]
548 1397
549 1398 #hei_interf
550 1399 if hei_interf is None:
551 1400 count_hei = num_hei/2 #Como es entero no importa
552 1401 hei_interf = numpy.asmatrix(range(count_hei)) + num_hei - count_hei
553 1402 hei_interf = numpy.asarray(hei_interf)[0]
554 1403 #nhei_interf
555 1404 if (nhei_interf == None):
556 1405 nhei_interf = 5
557 1406 if (nhei_interf < 1):
558 1407 nhei_interf = 1
559 1408 if (nhei_interf > count_hei):
560 1409 nhei_interf = count_hei
561 1410 if (offhei_interf == None):
562 1411 offhei_interf = 0
563 1412
564 1413 ind_hei = range(num_hei)
565 1414 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
566 1415 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
567 1416 mask_prof = numpy.asarray(range(num_prof))
568 1417 num_mask_prof = mask_prof.size
569 1418 comp_mask_prof = [0, num_prof/2]
570 1419
571 1420
572 1421 #noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
573 1422 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
574 1423 jnoise = numpy.nan
575 1424 noise_exist = jnoise[0] < numpy.Inf
576 1425
577 1426 #Subrutina de Remocion de la Interferencia
578 1427 for ich in range(num_channel):
579 1428 #Se ordena los espectros segun su potencia (menor a mayor)
580 1429 power = jspectra[ich,mask_prof,:]
581 1430 power = power[:,hei_interf]
582 1431 power = power.sum(axis = 0)
583 1432 psort = power.ravel().argsort()
584 1433
585 1434 #Se estima la interferencia promedio en los Espectros de Potencia empleando
586 1435 junkspc_interf = jspectra[ich,:,hei_interf[psort[range(offhei_interf, nhei_interf + offhei_interf)]]]
587 1436
588 1437 if noise_exist:
589 1438 # tmp_noise = jnoise[ich] / num_prof
590 1439 tmp_noise = jnoise[ich]
591 1440 junkspc_interf = junkspc_interf - tmp_noise
592 1441 #junkspc_interf[:,comp_mask_prof] = 0
593 1442
594 1443 jspc_interf = junkspc_interf.sum(axis = 0) / nhei_interf
595 1444 jspc_interf = jspc_interf.transpose()
596 1445 #Calculando el espectro de interferencia promedio
597 1446 noiseid = numpy.where(jspc_interf <= tmp_noise/ numpy.sqrt(num_incoh))
598 1447 noiseid = noiseid[0]
599 1448 cnoiseid = noiseid.size
600 1449 interfid = numpy.where(jspc_interf > tmp_noise/ numpy.sqrt(num_incoh))
601 1450 interfid = interfid[0]
602 1451 cinterfid = interfid.size
603 1452
604 1453 if (cnoiseid > 0): jspc_interf[noiseid] = 0
605 1454
606 1455 #Expandiendo los perfiles a limpiar
607 1456 if (cinterfid > 0):
608 1457 new_interfid = (numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof)%num_prof
609 1458 new_interfid = numpy.asarray(new_interfid)
610 1459 new_interfid = {x for x in new_interfid}
611 1460 new_interfid = numpy.array(list(new_interfid))
612 1461 new_cinterfid = new_interfid.size
613 1462 else: new_cinterfid = 0
614 1463
615 1464 for ip in range(new_cinterfid):
616 1465 ind = junkspc_interf[:,new_interfid[ip]].ravel().argsort()
617 1466 jspc_interf[new_interfid[ip]] = junkspc_interf[ind[nhei_interf/2],new_interfid[ip]]
618 1467
619 1468
620 1469 jspectra[ich,:,ind_hei] = jspectra[ich,:,ind_hei] - jspc_interf #Corregir indices
621 1470
622 1471 #Removiendo la interferencia del punto de mayor interferencia
623 1472 ListAux = jspc_interf[mask_prof].tolist()
624 1473 maxid = ListAux.index(max(ListAux))
625 1474
626 1475
627 1476 if cinterfid > 0:
628 1477 for ip in range(cinterfid*(interf == 2) - 1):
629 1478 ind = (jspectra[ich,interfid[ip],:] < tmp_noise*(1 + 1/numpy.sqrt(num_incoh))).nonzero()
630 1479 cind = len(ind)
631 1480
632 1481 if (cind > 0):
633 1482 jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/numpy.sqrt(num_incoh))
634 1483
635 1484 ind = numpy.array([-2,-1,1,2])
636 1485 xx = numpy.zeros([4,4])
637 1486
638 1487 for id1 in range(4):
639 1488 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
640 1489
641 1490 xx_inv = numpy.linalg.inv(xx)
642 1491 xx = xx_inv[:,0]
643 1492 ind = (ind + maxid + num_mask_prof)%num_mask_prof
644 1493 yy = jspectra[ich,mask_prof[ind],:]
645 1494 jspectra[ich,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
646 1495
647 1496
648 1497 indAux = (jspectra[ich,:,:] < tmp_noise*(1-1/numpy.sqrt(num_incoh))).nonzero()
649 1498 jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/numpy.sqrt(num_incoh))
650 1499
651 1500 #Remocion de Interferencia en el Cross Spectra
652 1501 if jcspectra is None: return jspectra, jcspectra
653 1502 num_pairs = jcspectra.size/(num_prof*num_hei)
654 1503 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
655 1504
656 1505 for ip in range(num_pairs):
657 1506
658 1507 #-------------------------------------------
659 1508
660 1509 cspower = numpy.abs(jcspectra[ip,mask_prof,:])
661 1510 cspower = cspower[:,hei_interf]
662 1511 cspower = cspower.sum(axis = 0)
663 1512
664 1513 cspsort = cspower.ravel().argsort()
665 1514 junkcspc_interf = jcspectra[ip,:,hei_interf[cspsort[range(offhei_interf, nhei_interf + offhei_interf)]]]
666 1515 junkcspc_interf = junkcspc_interf.transpose()
667 1516 jcspc_interf = junkcspc_interf.sum(axis = 1)/nhei_interf
668 1517
669 1518 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
670 1519
671 1520 median_real = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
672 1521 median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
673 1522 junkcspc_interf[comp_mask_prof,:] = numpy.complex(median_real, median_imag)
674 1523
675 1524 for iprof in range(num_prof):
676 1525 ind = numpy.abs(junkcspc_interf[iprof,:]).ravel().argsort()
677 1526 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf/2]]
678 1527
679 1528 #Removiendo la Interferencia
680 1529 jcspectra[ip,:,ind_hei] = jcspectra[ip,:,ind_hei] - jcspc_interf
681 1530
682 1531 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
683 1532 maxid = ListAux.index(max(ListAux))
684 1533
685 1534 ind = numpy.array([-2,-1,1,2])
686 1535 xx = numpy.zeros([4,4])
687 1536
688 1537 for id1 in range(4):
689 1538 xx[:,id1] = ind[id1]**numpy.asarray(range(4))
690 1539
691 1540 xx_inv = numpy.linalg.inv(xx)
692 1541 xx = xx_inv[:,0]
693 1542
694 1543 ind = (ind + maxid + num_mask_prof)%num_mask_prof
695 1544 yy = jcspectra[ip,mask_prof[ind],:]
696 1545 jcspectra[ip,mask_prof[maxid],:] = numpy.dot(yy.transpose(),xx)
697 1546
698 1547 #Guardar Resultados
699 1548 self.dataOut.data_spc = jspectra
700 1549 self.dataOut.data_cspc = jcspectra
701 1550
702 1551 return 1
703 1552
704 1553 def setRadarFrequency(self, frequency=None):
705 1554
706 1555 if frequency != None:
707 1556 self.dataOut.frequency = frequency
708 1557
709 1558 return 1
710 1559
711 1560 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
712 1561 #validacion de rango
713 1562 if minHei == None:
714 1563 minHei = self.dataOut.heightList[0]
715 1564
716 1565 if maxHei == None:
717 1566 maxHei = self.dataOut.heightList[-1]
718 1567
719 1568 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
720 1569 print('minHei: %.2f is out of the heights range'%(minHei))
721 1570 print('minHei is setting to %.2f'%(self.dataOut.heightList[0]))
722 1571 minHei = self.dataOut.heightList[0]
723 1572
724 1573 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
725 1574 print('maxHei: %.2f is out of the heights range'%(maxHei))
726 1575 print('maxHei is setting to %.2f'%(self.dataOut.heightList[-1]))
727 1576 maxHei = self.dataOut.heightList[-1]
728 1577
729 1578 # validacion de velocidades
730 1579 velrange = self.dataOut.getVelRange(1)
731 1580
732 1581 if minVel == None:
733 1582 minVel = velrange[0]
734 1583
735 1584 if maxVel == None:
736 1585 maxVel = velrange[-1]
737 1586
738 1587 if (minVel < velrange[0]) or (minVel > maxVel):
739 1588 print('minVel: %.2f is out of the velocity range'%(minVel))
740 1589 print('minVel is setting to %.2f'%(velrange[0]))
741 1590 minVel = velrange[0]
742 1591
743 1592 if (maxVel > velrange[-1]) or (maxVel < minVel):
744 1593 print('maxVel: %.2f is out of the velocity range'%(maxVel))
745 1594 print('maxVel is setting to %.2f'%(velrange[-1]))
746 1595 maxVel = velrange[-1]
747 1596
748 1597 # seleccion de indices para rango
749 1598 minIndex = 0
750 1599 maxIndex = 0
751 1600 heights = self.dataOut.heightList
752 1601
753 1602 inda = numpy.where(heights >= minHei)
754 1603 indb = numpy.where(heights <= maxHei)
755 1604
756 1605 try:
757 1606 minIndex = inda[0][0]
758 1607 except:
759 1608 minIndex = 0
760 1609
761 1610 try:
762 1611 maxIndex = indb[0][-1]
763 1612 except:
764 1613 maxIndex = len(heights)
765 1614
766 1615 if (minIndex < 0) or (minIndex > maxIndex):
767 1616 raise ValueError("some value in (%d,%d) is not valid" % (minIndex, maxIndex))
768 1617
769 1618 if (maxIndex >= self.dataOut.nHeights):
770 1619 maxIndex = self.dataOut.nHeights-1
771 1620
772 1621 # seleccion de indices para velocidades
773 1622 indminvel = numpy.where(velrange >= minVel)
774 1623 indmaxvel = numpy.where(velrange <= maxVel)
775 1624 try:
776 1625 minIndexVel = indminvel[0][0]
777 1626 except:
778 1627 minIndexVel = 0
779 1628
780 1629 try:
781 1630 maxIndexVel = indmaxvel[0][-1]
782 1631 except:
783 1632 maxIndexVel = len(velrange)
784 1633
785 1634 #seleccion del espectro
786 1635 data_spc = self.dataOut.data_spc[:,minIndexVel:maxIndexVel+1,minIndex:maxIndex+1]
787 1636 #estimacion de ruido
788 1637 noise = numpy.zeros(self.dataOut.nChannels)
789 1638
790 1639 for channel in range(self.dataOut.nChannels):
791 1640 daux = data_spc[channel,:,:]
792 1641 noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
793 1642
794 1643 self.dataOut.noise_estimation = noise.copy()
795 1644
796 1645 return 1
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
Show file before | Show file after | Hide comments
1 NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff
General Comments 0
You need to be logged in to leave comments. Login now