##// END OF EJS Templates
Update plots modules
Juan C. Espinoza -
r1308:fe4d0955eff6
parent child
Show More
@@ -1,648 +1,648
1 1 '''
2 2 Main routines to create a Signal Chain project
3 3 '''
4 4
5 5 import re
6 6 import sys
7 7 import ast
8 8 import datetime
9 9 import traceback
10 10 import time
11 11 from multiprocessing import Process, Queue
12 12 from threading import Thread
13 13 from xml.etree.ElementTree import ElementTree, Element, SubElement
14 14
15 15 from schainpy.admin import Alarm, SchainWarning
16 16 from schainpy.model import *
17 17 from schainpy.utils import log
18 18
19 19
20 20 class ConfBase():
21 21
22 22 def __init__(self):
23 23
24 24 self.id = '0'
25 25 self.name = None
26 26 self.priority = None
27 27 self.parameters = {}
28 28 self.object = None
29 29 self.operations = []
30 30
31 31 def getId(self):
32 32
33 33 return self.id
34 34
35 35 def getNewId(self):
36 36
37 37 return int(self.id) * 10 + len(self.operations) + 1
38 38
39 39 def updateId(self, new_id):
40 40
41 41 self.id = str(new_id)
42 42
43 43 n = 1
44 44 for conf in self.operations:
45 45 conf_id = str(int(new_id) * 10 + n)
46 46 conf.updateId(conf_id)
47 47 n += 1
48 48
49 49 def getKwargs(self):
50 50
51 51 params = {}
52 52
53 53 for key, value in self.parameters.items():
54 54 if value not in (None, '', ' '):
55 55 params[key] = value
56 56
57 57 return params
58 58
59 59 def update(self, **kwargs):
60 60
61 61 for key, value in kwargs.items():
62 62 self.addParameter(name=key, value=value)
63 63
64 64 def addParameter(self, name, value, format=None):
65 65 '''
66 66 '''
67 67
68 68 if isinstance(value, str) and re.search(r'(\d+/\d+/\d+)', value):
69 69 self.parameters[name] = datetime.date(*[int(x) for x in value.split('/')])
70 70 elif isinstance(value, str) and re.search(r'(\d+:\d+:\d+)', value):
71 71 self.parameters[name] = datetime.time(*[int(x) for x in value.split(':')])
72 72 else:
73 73 try:
74 74 self.parameters[name] = ast.literal_eval(value)
75 75 except:
76 76 if isinstance(value, str) and ',' in value:
77 77 self.parameters[name] = value.split(',')
78 78 else:
79 79 self.parameters[name] = value
80 80
81 81 def getParameters(self):
82 82
83 83 params = {}
84 84 for key, value in self.parameters.items():
85 85 s = type(value).__name__
86 86 if s == 'date':
87 87 params[key] = value.strftime('%Y/%m/%d')
88 88 elif s == 'time':
89 89 params[key] = value.strftime('%H:%M:%S')
90 90 else:
91 91 params[key] = str(value)
92 92
93 93 return params
94 94
95 95 def makeXml(self, element):
96 96
97 97 xml = SubElement(element, self.ELEMENTNAME)
98 98 for label in self.xml_labels:
99 99 xml.set(label, str(getattr(self, label)))
100 100
101 101 for key, value in self.getParameters().items():
102 102 xml_param = SubElement(xml, 'Parameter')
103 103 xml_param.set('name', key)
104 104 xml_param.set('value', value)
105 105
106 106 for conf in self.operations:
107 107 conf.makeXml(xml)
108 108
109 109 def __str__(self):
110 110
111 111 if self.ELEMENTNAME == 'Operation':
112 112 s = ' {}[id={}]\n'.format(self.name, self.id)
113 113 else:
114 114 s = '{}[id={}, inputId={}]\n'.format(self.name, self.id, self.inputId)
115 115
116 116 for key, value in self.parameters.items():
117 117 if self.ELEMENTNAME == 'Operation':
118 118 s += ' {}: {}\n'.format(key, value)
119 119 else:
120 120 s += ' {}: {}\n'.format(key, value)
121 121
122 122 for conf in self.operations:
123 123 s += str(conf)
124 124
125 125 return s
126 126
127 127 class OperationConf(ConfBase):
128 128
129 129 ELEMENTNAME = 'Operation'
130 130 xml_labels = ['id', 'name']
131 131
132 132 def setup(self, id, name, priority, project_id, err_queue):
133 133
134 134 self.id = str(id)
135 135 self.project_id = project_id
136 136 self.name = name
137 137 self.type = 'other'
138 138 self.err_queue = err_queue
139 139
140 140 def readXml(self, element, project_id, err_queue):
141 141
142 142 self.id = element.get('id')
143 143 self.name = element.get('name')
144 144 self.type = 'other'
145 145 self.project_id = str(project_id)
146 146 self.err_queue = err_queue
147 147
148 148 for elm in element.iter('Parameter'):
149 149 self.addParameter(elm.get('name'), elm.get('value'))
150 150
151 151 def createObject(self):
152 152
153 153 className = eval(self.name)
154 154
155 if 'Plot' in self.name or 'Writer' in self.name or 'Send' in self.name:
155 if 'Plot' in self.name or 'Writer' in self.name or 'Send' in self.name or 'print' in self.name:
156 156 kwargs = self.getKwargs()
157 157 opObj = className(self.id, self.id, self.project_id, self.err_queue, **kwargs)
158 158 opObj.start()
159 159 self.type = 'external'
160 160 else:
161 161 opObj = className()
162 162
163 163 self.object = opObj
164 164 return opObj
165 165
166 166 class ProcUnitConf(ConfBase):
167 167
168 168 ELEMENTNAME = 'ProcUnit'
169 169 xml_labels = ['id', 'inputId', 'name']
170 170
171 171 def setup(self, project_id, id, name, datatype, inputId, err_queue):
172 172 '''
173 173 '''
174 174
175 175 if datatype == None and name == None:
176 176 raise ValueError('datatype or name should be defined')
177 177
178 178 if name == None:
179 179 if 'Proc' in datatype:
180 180 name = datatype
181 181 else:
182 182 name = '%sProc' % (datatype)
183 183
184 184 if datatype == None:
185 185 datatype = name.replace('Proc', '')
186 186
187 187 self.id = str(id)
188 188 self.project_id = project_id
189 189 self.name = name
190 190 self.datatype = datatype
191 191 self.inputId = inputId
192 192 self.err_queue = err_queue
193 193 self.operations = []
194 194 self.parameters = {}
195 195
196 196 def removeOperation(self, id):
197 197
198 198 i = [1 if x.id==id else 0 for x in self.operations]
199 199 self.operations.pop(i.index(1))
200 200
201 201 def getOperation(self, id):
202 202
203 203 for conf in self.operations:
204 204 if conf.id == id:
205 205 return conf
206 206
207 207 def addOperation(self, name, optype='self'):
208 208 '''
209 209 '''
210 210
211 211 id = self.getNewId()
212 212 conf = OperationConf()
213 213 conf.setup(id, name=name, priority='0', project_id=self.project_id, err_queue=self.err_queue)
214 214 self.operations.append(conf)
215 215
216 216 return conf
217 217
218 218 def readXml(self, element, project_id, err_queue):
219 219
220 220 self.id = element.get('id')
221 221 self.name = element.get('name')
222 222 self.inputId = None if element.get('inputId') == 'None' else element.get('inputId')
223 223 self.datatype = element.get('datatype', self.name.replace(self.ELEMENTNAME.replace('Unit', ''), ''))
224 224 self.project_id = str(project_id)
225 225 self.err_queue = err_queue
226 226 self.operations = []
227 227 self.parameters = {}
228 228
229 229 for elm in element:
230 230 if elm.tag == 'Parameter':
231 231 self.addParameter(elm.get('name'), elm.get('value'))
232 232 elif elm.tag == 'Operation':
233 233 conf = OperationConf()
234 234 conf.readXml(elm, project_id, err_queue)
235 235 self.operations.append(conf)
236 236
237 237 def createObjects(self):
238 238 '''
239 239 Instancia de unidades de procesamiento.
240 240 '''
241 241
242 242 className = eval(self.name)
243 243 kwargs = self.getKwargs()
244 244 procUnitObj = className()
245 245 procUnitObj.name = self.name
246 246 log.success('creating process...', self.name)
247 247
248 248 for conf in self.operations:
249 249
250 250 opObj = conf.createObject()
251 251
252 252 log.success('adding operation: {}, type:{}'.format(
253 253 conf.name,
254 254 conf.type), self.name)
255 255
256 256 procUnitObj.addOperation(conf, opObj)
257 257
258 258 self.object = procUnitObj
259 259
260 260 def run(self):
261 261 '''
262 262 '''
263 263
264 264 return self.object.call(**self.getKwargs())
265 265
266 266
267 267 class ReadUnitConf(ProcUnitConf):
268 268
269 269 ELEMENTNAME = 'ReadUnit'
270 270
271 271 def __init__(self):
272 272
273 273 self.id = None
274 274 self.datatype = None
275 275 self.name = None
276 276 self.inputId = None
277 277 self.operations = []
278 278 self.parameters = {}
279 279
280 280 def setup(self, project_id, id, name, datatype, err_queue, path='', startDate='', endDate='',
281 281 startTime='', endTime='', server=None, **kwargs):
282 282
283 283 if datatype == None and name == None:
284 284 raise ValueError('datatype or name should be defined')
285 285 if name == None:
286 286 if 'Reader' in datatype:
287 287 name = datatype
288 288 datatype = name.replace('Reader','')
289 289 else:
290 290 name = '{}Reader'.format(datatype)
291 291 if datatype == None:
292 292 if 'Reader' in name:
293 293 datatype = name.replace('Reader','')
294 294 else:
295 295 datatype = name
296 296 name = '{}Reader'.format(name)
297 297
298 298 self.id = id
299 299 self.project_id = project_id
300 300 self.name = name
301 301 self.datatype = datatype
302 302 self.err_queue = err_queue
303 303
304 304 self.addParameter(name='path', value=path)
305 305 self.addParameter(name='startDate', value=startDate)
306 306 self.addParameter(name='endDate', value=endDate)
307 307 self.addParameter(name='startTime', value=startTime)
308 308 self.addParameter(name='endTime', value=endTime)
309 309
310 310 for key, value in kwargs.items():
311 311 self.addParameter(name=key, value=value)
312 312
313 313
314 314 class Project(Process):
315 315
316 316 ELEMENTNAME = 'Project'
317 317
318 318 def __init__(self):
319 319
320 320 Process.__init__(self)
321 321 self.id = None
322 322 self.filename = None
323 323 self.description = None
324 324 self.email = None
325 325 self.alarm = []
326 326 self.configurations = {}
327 327 # self.err_queue = Queue()
328 328 self.err_queue = None
329 329 self.started = False
330 330
331 331 def getNewId(self):
332 332
333 333 idList = list(self.configurations.keys())
334 334 id = int(self.id) * 10
335 335
336 336 while True:
337 337 id += 1
338 338
339 339 if str(id) in idList:
340 340 continue
341 341
342 342 break
343 343
344 344 return str(id)
345 345
346 346 def updateId(self, new_id):
347 347
348 348 self.id = str(new_id)
349 349
350 350 keyList = list(self.configurations.keys())
351 351 keyList.sort()
352 352
353 353 n = 1
354 354 new_confs = {}
355 355
356 356 for procKey in keyList:
357 357
358 358 conf = self.configurations[procKey]
359 359 idProcUnit = str(int(self.id) * 10 + n)
360 360 conf.updateId(idProcUnit)
361 361 new_confs[idProcUnit] = conf
362 362 n += 1
363 363
364 364 self.configurations = new_confs
365 365
366 366 def setup(self, id=1, name='', description='', email=None, alarm=[]):
367 367
368 368 self.id = str(id)
369 369 self.description = description
370 370 self.email = email
371 371 self.alarm = alarm
372 372 if name:
373 373 self.name = '{} ({})'.format(Process.__name__, name)
374 374
375 375 def update(self, **kwargs):
376 376
377 377 for key, value in kwargs.items():
378 378 setattr(self, key, value)
379 379
380 380 def clone(self):
381 381
382 382 p = Project()
383 383 p.id = self.id
384 384 p.name = self.name
385 385 p.description = self.description
386 386 p.configurations = self.configurations.copy()
387 387
388 388 return p
389 389
390 390 def addReadUnit(self, id=None, datatype=None, name=None, **kwargs):
391 391
392 392 '''
393 393 '''
394 394
395 395 if id is None:
396 396 idReadUnit = self.getNewId()
397 397 else:
398 398 idReadUnit = str(id)
399 399
400 400 conf = ReadUnitConf()
401 401 conf.setup(self.id, idReadUnit, name, datatype, self.err_queue, **kwargs)
402 402 self.configurations[conf.id] = conf
403 403
404 404 return conf
405 405
406 406 def addProcUnit(self, id=None, inputId='0', datatype=None, name=None):
407 407
408 408 '''
409 409 '''
410 410
411 411 if id is None:
412 412 idProcUnit = self.getNewId()
413 413 else:
414 414 idProcUnit = id
415 415
416 416 conf = ProcUnitConf()
417 417 conf.setup(self.id, idProcUnit, name, datatype, inputId, self.err_queue)
418 418 self.configurations[conf.id] = conf
419 419
420 420 return conf
421 421
422 422 def removeProcUnit(self, id):
423 423
424 424 if id in self.configurations:
425 425 self.configurations.pop(id)
426 426
427 427 def getReadUnit(self):
428 428
429 429 for obj in list(self.configurations.values()):
430 430 if obj.ELEMENTNAME == 'ReadUnit':
431 431 return obj
432 432
433 433 return None
434 434
435 435 def getProcUnit(self, id):
436 436
437 437 return self.configurations[id]
438 438
439 439 def getUnits(self):
440 440
441 441 keys = list(self.configurations)
442 442 keys.sort()
443 443
444 444 for key in keys:
445 445 yield self.configurations[key]
446 446
447 447 def updateUnit(self, id, **kwargs):
448 448
449 449 conf = self.configurations[id].update(**kwargs)
450 450
451 451 def makeXml(self):
452 452
453 453 xml = Element('Project')
454 454 xml.set('id', str(self.id))
455 455 xml.set('name', self.name)
456 456 xml.set('description', self.description)
457 457
458 458 for conf in self.configurations.values():
459 459 conf.makeXml(xml)
460 460
461 461 self.xml = xml
462 462
463 463 def writeXml(self, filename=None):
464 464
465 465 if filename == None:
466 466 if self.filename:
467 467 filename = self.filename
468 468 else:
469 469 filename = 'schain.xml'
470 470
471 471 if not filename:
472 472 print('filename has not been defined. Use setFilename(filename) for do it.')
473 473 return 0
474 474
475 475 abs_file = os.path.abspath(filename)
476 476
477 477 if not os.access(os.path.dirname(abs_file), os.W_OK):
478 478 print('No write permission on %s' % os.path.dirname(abs_file))
479 479 return 0
480 480
481 481 if os.path.isfile(abs_file) and not(os.access(abs_file, os.W_OK)):
482 482 print('File %s already exists and it could not be overwriten' % abs_file)
483 483 return 0
484 484
485 485 self.makeXml()
486 486
487 487 ElementTree(self.xml).write(abs_file, method='xml')
488 488
489 489 self.filename = abs_file
490 490
491 491 return 1
492 492
493 493 def readXml(self, filename):
494 494
495 495 abs_file = os.path.abspath(filename)
496 496
497 497 self.configurations = {}
498 498
499 499 try:
500 500 self.xml = ElementTree().parse(abs_file)
501 501 except:
502 502 log.error('Error reading %s, verify file format' % filename)
503 503 return 0
504 504
505 505 self.id = self.xml.get('id')
506 506 self.name = self.xml.get('name')
507 507 self.description = self.xml.get('description')
508 508
509 509 for element in self.xml:
510 510 if element.tag == 'ReadUnit':
511 511 conf = ReadUnitConf()
512 512 conf.readXml(element, self.id, self.err_queue)
513 513 self.configurations[conf.id] = conf
514 514 elif element.tag == 'ProcUnit':
515 515 conf = ProcUnitConf()
516 516 input_proc = self.configurations[element.get('inputId')]
517 517 conf.readXml(element, self.id, self.err_queue)
518 518 self.configurations[conf.id] = conf
519 519
520 520 self.filename = abs_file
521 521
522 522 return 1
523 523
524 524 def __str__(self):
525 525
526 526 text = '\nProject[id=%s, name=%s, description=%s]\n\n' % (
527 527 self.id,
528 528 self.name,
529 529 self.description,
530 530 )
531 531
532 532 for conf in self.configurations.values():
533 533 text += '{}'.format(conf)
534 534
535 535 return text
536 536
537 537 def createObjects(self):
538 538
539 539 keys = list(self.configurations.keys())
540 540 keys.sort()
541 541 for key in keys:
542 542 conf = self.configurations[key]
543 543 conf.createObjects()
544 544 if conf.inputId is not None:
545 545 conf.object.setInput(self.configurations[conf.inputId].object)
546 546
547 547 def monitor(self):
548 548
549 549 t = Thread(target=self._monitor, args=(self.err_queue, self.ctx))
550 550 t.start()
551 551
552 552 def _monitor(self, queue, ctx):
553 553
554 554 import socket
555 555
556 556 procs = 0
557 557 err_msg = ''
558 558
559 559 while True:
560 560 msg = queue.get()
561 561 if '#_start_#' in msg:
562 562 procs += 1
563 563 elif '#_end_#' in msg:
564 564 procs -=1
565 565 else:
566 566 err_msg = msg
567 567
568 568 if procs == 0 or 'Traceback' in err_msg:
569 569 break
570 570 time.sleep(0.1)
571 571
572 572 if '|' in err_msg:
573 573 name, err = err_msg.split('|')
574 574 if 'SchainWarning' in err:
575 575 log.warning(err.split('SchainWarning:')[-1].split('\n')[0].strip(), name)
576 576 elif 'SchainError' in err:
577 577 log.error(err.split('SchainError:')[-1].split('\n')[0].strip(), name)
578 578 else:
579 579 log.error(err, name)
580 580 else:
581 581 name, err = self.name, err_msg
582 582
583 583 time.sleep(1)
584 584
585 585 ctx.term()
586 586
587 587 message = ''.join(err)
588 588
589 589 if err_msg:
590 590 subject = 'SChain v%s: Error running %s\n' % (
591 591 schainpy.__version__, self.name)
592 592
593 593 subtitle = 'Hostname: %s\n' % socket.gethostbyname(
594 594 socket.gethostname())
595 595 subtitle += 'Working directory: %s\n' % os.path.abspath('./')
596 596 subtitle += 'Configuration file: %s\n' % self.filename
597 597 subtitle += 'Time: %s\n' % str(datetime.datetime.now())
598 598
599 599 readUnitConfObj = self.getReadUnit()
600 600 if readUnitConfObj:
601 601 subtitle += '\nInput parameters:\n'
602 602 subtitle += '[Data path = %s]\n' % readUnitConfObj.parameters['path']
603 603 subtitle += '[Start date = %s]\n' % readUnitConfObj.parameters['startDate']
604 604 subtitle += '[End date = %s]\n' % readUnitConfObj.parameters['endDate']
605 605 subtitle += '[Start time = %s]\n' % readUnitConfObj.parameters['startTime']
606 606 subtitle += '[End time = %s]\n' % readUnitConfObj.parameters['endTime']
607 607
608 608 a = Alarm(
609 609 modes=self.alarm,
610 610 email=self.email,
611 611 message=message,
612 612 subject=subject,
613 613 subtitle=subtitle,
614 614 filename=self.filename
615 615 )
616 616
617 617 a.start()
618 618
619 619 def setFilename(self, filename):
620 620
621 621 self.filename = filename
622 622
623 623 def runProcs(self):
624 624
625 625 err = False
626 626 n = len(self.configurations)
627 627
628 628 while not err:
629 629 for conf in self.getUnits():
630 630 ok = conf.run()
631 631 if ok is 'Error':
632 632 n -= 1
633 633 continue
634 634 elif not ok:
635 635 break
636 636 if n == 0:
637 637 err = True
638 638
639 639 def run(self):
640 640
641 641 log.success('\nStarting Project {} [id={}]'.format(self.name, self.id), tag='')
642 642 self.started = True
643 643 self.start_time = time.time()
644 644 self.createObjects()
645 645 self.runProcs()
646 646 log.success('{} Done (Time: {:4.2f}s)'.format(
647 647 self.name,
648 648 time.time()-self.start_time), '')
@@ -1,1394 +1,1391
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
5 5 '''
6 6
7 7 import copy
8 8 import numpy
9 9 import datetime
10 10 import json
11 11
12 12 import schainpy.admin
13 13 from schainpy.utils import log
14 14 from .jroheaderIO import SystemHeader, RadarControllerHeader
15 15 from schainpy.model.data import _noise
16 16
17 17
18 18 def getNumpyDtype(dataTypeCode):
19 19
20 20 if dataTypeCode == 0:
21 21 numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
22 22 elif dataTypeCode == 1:
23 23 numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
24 24 elif dataTypeCode == 2:
25 25 numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
26 26 elif dataTypeCode == 3:
27 27 numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
28 28 elif dataTypeCode == 4:
29 29 numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
30 30 elif dataTypeCode == 5:
31 31 numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
32 32 else:
33 33 raise ValueError('dataTypeCode was not defined')
34 34
35 35 return numpyDtype
36 36
37 37
38 38 def getDataTypeCode(numpyDtype):
39 39
40 40 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
41 41 datatype = 0
42 42 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
43 43 datatype = 1
44 44 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
45 45 datatype = 2
46 46 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
47 47 datatype = 3
48 48 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
49 49 datatype = 4
50 50 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
51 51 datatype = 5
52 52 else:
53 53 datatype = None
54 54
55 55 return datatype
56 56
57 57
58 58 def hildebrand_sekhon(data, navg):
59 59 """
60 60 This method is for the objective determination of the noise level in Doppler spectra. This
61 61 implementation technique is based on the fact that the standard deviation of the spectral
62 62 densities is equal to the mean spectral density for white Gaussian noise
63 63
64 64 Inputs:
65 65 Data : heights
66 66 navg : numbers of averages
67 67
68 68 Return:
69 69 mean : noise's level
70 70 """
71 71
72 72 sortdata = numpy.sort(data, axis=None)
73 73 '''
74 74 lenOfData = len(sortdata)
75 75 nums_min = lenOfData*0.2
76 76
77 77 if nums_min <= 5:
78 78
79 79 nums_min = 5
80 80
81 81 sump = 0.
82 82 sumq = 0.
83 83
84 84 j = 0
85 85 cont = 1
86 86
87 87 while((cont == 1)and(j < lenOfData)):
88 88
89 89 sump += sortdata[j]
90 90 sumq += sortdata[j]**2
91 91
92 92 if j > nums_min:
93 93 rtest = float(j)/(j-1) + 1.0/navg
94 94 if ((sumq*j) > (rtest*sump**2)):
95 95 j = j - 1
96 96 sump = sump - sortdata[j]
97 97 sumq = sumq - sortdata[j]**2
98 98 cont = 0
99 99
100 100 j += 1
101 101
102 102 lnoise = sump / j
103 103 '''
104 104 return _noise.hildebrand_sekhon(sortdata, navg)
105 105
106 106
107 107 class Beam:
108 108
109 109 def __init__(self):
110 110 self.codeList = []
111 111 self.azimuthList = []
112 112 self.zenithList = []
113 113
114 114
115 115 class GenericData(object):
116 116
117 117 flagNoData = True
118 118
119 119 def copy(self, inputObj=None):
120 120
121 121 if inputObj == None:
122 122 return copy.deepcopy(self)
123 123
124 124 for key in list(inputObj.__dict__.keys()):
125 125
126 126 attribute = inputObj.__dict__[key]
127 127
128 128 # If this attribute is a tuple or list
129 129 if type(inputObj.__dict__[key]) in (tuple, list):
130 130 self.__dict__[key] = attribute[:]
131 131 continue
132 132
133 133 # If this attribute is another object or instance
134 134 if hasattr(attribute, '__dict__'):
135 135 self.__dict__[key] = attribute.copy()
136 136 continue
137 137
138 138 self.__dict__[key] = inputObj.__dict__[key]
139 139
140 140 def deepcopy(self):
141 141
142 142 return copy.deepcopy(self)
143 143
144 144 def isEmpty(self):
145 145
146 146 return self.flagNoData
147 147
148 148 def isReady(self):
149 149
150 150 return not self.flagNoData
151 151
152 152
153 153 class JROData(GenericData):
154 154
155 155 # m_BasicHeader = BasicHeader()
156 156 # m_ProcessingHeader = ProcessingHeader()
157 157
158 158 systemHeaderObj = SystemHeader()
159 159 radarControllerHeaderObj = RadarControllerHeader()
160 160 # data = None
161 161 type = None
162 162 datatype = None # dtype but in string
163 163 # dtype = None
164 164 # nChannels = None
165 165 # nHeights = None
166 166 nProfiles = None
167 167 heightList = None
168 168 channelList = None
169 169 flagDiscontinuousBlock = False
170 170 useLocalTime = False
171 171 utctime = None
172 172 timeZone = None
173 173 dstFlag = None
174 174 errorCount = None
175 175 blocksize = None
176 176 # nCode = None
177 177 # nBaud = None
178 178 # code = None
179 179 flagDecodeData = False # asumo q la data no esta decodificada
180 180 flagDeflipData = False # asumo q la data no esta sin flip
181 181 flagShiftFFT = False
182 182 # ippSeconds = None
183 183 # timeInterval = None
184 184 nCohInt = None
185 185 # noise = None
186 186 windowOfFilter = 1
187 187 # Speed of ligth
188 188 C = 3e8
189 189 frequency = 49.92e6
190 190 realtime = False
191 191 beacon_heiIndexList = None
192 192 last_block = None
193 193 blocknow = None
194 194 azimuth = None
195 195 zenith = None
196 196 beam = Beam()
197 197 profileIndex = None
198 198 error = None
199 199 data = None
200 200 nmodes = None
201 201
202 202 def __str__(self):
203 203
204 204 return '{} - {}'.format(self.type, self.getDatatime())
205 205
206 206 def getNoise(self):
207 207
208 208 raise NotImplementedError
209 209
210 210 def getNChannels(self):
211 211
212 212 return len(self.channelList)
213 213
214 214 def getChannelIndexList(self):
215 215
216 216 return list(range(self.nChannels))
217 217
218 218 def getNHeights(self):
219 219
220 220 return len(self.heightList)
221 221
222 222 def getHeiRange(self, extrapoints=0):
223 223
224 224 heis = self.heightList
225 225 # deltah = self.heightList[1] - self.heightList[0]
226 226 #
227 227 # heis.append(self.heightList[-1])
228 228
229 229 return heis
230 230
231 231 def getDeltaH(self):
232 232
233 233 delta = self.heightList[1] - self.heightList[0]
234 234
235 235 return delta
236 236
237 237 def getltctime(self):
238 238
239 239 if self.useLocalTime:
240 240 return self.utctime - self.timeZone * 60
241 241
242 242 return self.utctime
243 243
244 244 def getDatatime(self):
245 245
246 246 datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
247 247 return datatimeValue
248 248
249 249 def getTimeRange(self):
250 250
251 251 datatime = []
252 252
253 253 datatime.append(self.ltctime)
254 254 datatime.append(self.ltctime + self.timeInterval + 1)
255 255
256 256 datatime = numpy.array(datatime)
257 257
258 258 return datatime
259 259
260 260 def getFmaxTimeResponse(self):
261 261
262 262 period = (10**-6) * self.getDeltaH() / (0.15)
263 263
264 264 PRF = 1. / (period * self.nCohInt)
265 265
266 266 fmax = PRF
267 267
268 268 return fmax
269 269
270 270 def getFmax(self):
271 271 PRF = 1. / (self.ippSeconds * self.nCohInt)
272 272
273 273 fmax = PRF
274 274 return fmax
275 275
276 276 def getVmax(self):
277 277
278 278 _lambda = self.C / self.frequency
279 279
280 280 vmax = self.getFmax() * _lambda / 2
281 281
282 282 return vmax
283 283
284 284 def get_ippSeconds(self):
285 285 '''
286 286 '''
287 287 return self.radarControllerHeaderObj.ippSeconds
288 288
289 289 def set_ippSeconds(self, ippSeconds):
290 290 '''
291 291 '''
292 292
293 293 self.radarControllerHeaderObj.ippSeconds = ippSeconds
294 294
295 295 return
296 296
297 297 def get_dtype(self):
298 298 '''
299 299 '''
300 300 return getNumpyDtype(self.datatype)
301 301
302 302 def set_dtype(self, numpyDtype):
303 303 '''
304 304 '''
305 305
306 306 self.datatype = getDataTypeCode(numpyDtype)
307 307
308 308 def get_code(self):
309 309 '''
310 310 '''
311 311 return self.radarControllerHeaderObj.code
312 312
313 313 def set_code(self, code):
314 314 '''
315 315 '''
316 316 self.radarControllerHeaderObj.code = code
317 317
318 318 return
319 319
320 320 def get_ncode(self):
321 321 '''
322 322 '''
323 323 return self.radarControllerHeaderObj.nCode
324 324
325 325 def set_ncode(self, nCode):
326 326 '''
327 327 '''
328 328 self.radarControllerHeaderObj.nCode = nCode
329 329
330 330 return
331 331
332 332 def get_nbaud(self):
333 333 '''
334 334 '''
335 335 return self.radarControllerHeaderObj.nBaud
336 336
337 337 def set_nbaud(self, nBaud):
338 338 '''
339 339 '''
340 340 self.radarControllerHeaderObj.nBaud = nBaud
341 341
342 342 return
343 343
344 344 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
345 345 channelIndexList = property(
346 346 getChannelIndexList, "I'm the 'channelIndexList' property.")
347 347 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
348 348 #noise = property(getNoise, "I'm the 'nHeights' property.")
349 349 datatime = property(getDatatime, "I'm the 'datatime' property")
350 350 ltctime = property(getltctime, "I'm the 'ltctime' property")
351 351 ippSeconds = property(get_ippSeconds, set_ippSeconds)
352 352 dtype = property(get_dtype, set_dtype)
353 353 # timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
354 354 code = property(get_code, set_code)
355 355 nCode = property(get_ncode, set_ncode)
356 356 nBaud = property(get_nbaud, set_nbaud)
357 357
358 358
359 359 class Voltage(JROData):
360 360
361 361 # data es un numpy array de 2 dmensiones (canales, alturas)
362 362 data = None
363 363 data_intensity = None
364 364 data_velocity = None
365 365 def __init__(self):
366 366 '''
367 367 Constructor
368 368 '''
369 369
370 370 self.useLocalTime = True
371 371 self.radarControllerHeaderObj = RadarControllerHeader()
372 372 self.systemHeaderObj = SystemHeader()
373 373 self.type = "Voltage"
374 374 self.data = None
375 375 # self.dtype = None
376 376 # self.nChannels = 0
377 377 # self.nHeights = 0
378 378 self.nProfiles = None
379 379 self.heightList = None
380 380 self.channelList = None
381 381 # self.channelIndexList = None
382 382 self.flagNoData = True
383 383 self.flagDiscontinuousBlock = False
384 384 self.utctime = None
385 385 self.timeZone = None
386 386 self.dstFlag = None
387 387 self.errorCount = None
388 388 self.nCohInt = None
389 389 self.blocksize = None
390 390 self.flagDecodeData = False # asumo q la data no esta decodificada
391 391 self.flagDeflipData = False # asumo q la data no esta sin flip
392 392 self.flagShiftFFT = False
393 393 self.flagDataAsBlock = False # Asumo que la data es leida perfil a perfil
394 394 self.profileIndex = 0
395 395
396 396 def getNoisebyHildebrand(self, channel=None):
397 397 """
398 398 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
399 399
400 400 Return:
401 401 noiselevel
402 402 """
403 403
404 404 if channel != None:
405 405 data = self.data[channel]
406 406 nChannels = 1
407 407 else:
408 408 data = self.data
409 409 nChannels = self.nChannels
410 410
411 411 noise = numpy.zeros(nChannels)
412 412 power = data * numpy.conjugate(data)
413 413
414 414 for thisChannel in range(nChannels):
415 415 if nChannels == 1:
416 416 daux = power[:].real
417 417 else:
418 418 daux = power[thisChannel, :].real
419 419 noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
420 420
421 421 return noise
422 422
423 423 def getNoise(self, type=1, channel=None):
424 424
425 425 if type == 1:
426 426 noise = self.getNoisebyHildebrand(channel)
427 427
428 428 return noise
429 429
430 430 def getPower(self, channel=None):
431 431
432 432 if channel != None:
433 433 data = self.data[channel]
434 434 else:
435 435 data = self.data
436 436
437 437 power = data * numpy.conjugate(data)
438 438 powerdB = 10 * numpy.log10(power.real)
439 439 powerdB = numpy.squeeze(powerdB)
440 440
441 441 return powerdB
442 442
443 443 def getTimeInterval(self):
444 444
445 445 timeInterval = self.ippSeconds * self.nCohInt
446 446
447 447 return timeInterval
448 448
449 449 noise = property(getNoise, "I'm the 'nHeights' property.")
450 450 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
451 451
452 452
453 453 class Spectra(JROData):
454 454
455 455 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
456 456 data_spc = None
457 457 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
458 458 data_cspc = None
459 459 # data dc es un numpy array de 2 dmensiones (canales, alturas)
460 460 data_dc = None
461 461 # data power
462 462 data_pwr = None
463 463 nFFTPoints = None
464 464 # nPairs = None
465 465 pairsList = None
466 466 nIncohInt = None
467 467 wavelength = None # Necesario para cacular el rango de velocidad desde la frecuencia
468 468 nCohInt = None # se requiere para determinar el valor de timeInterval
469 469 ippFactor = None
470 470 profileIndex = 0
471 471 plotting = "spectra"
472 472
473 473 def __init__(self):
474 474 '''
475 475 Constructor
476 476 '''
477 477
478 478 self.useLocalTime = True
479 479 self.radarControllerHeaderObj = RadarControllerHeader()
480 480 self.systemHeaderObj = SystemHeader()
481 481 self.type = "Spectra"
482 482 # self.data = None
483 483 # self.dtype = None
484 484 # self.nChannels = 0
485 485 # self.nHeights = 0
486 486 self.nProfiles = None
487 487 self.heightList = None
488 488 self.channelList = None
489 489 # self.channelIndexList = None
490 490 self.pairsList = None
491 491 self.flagNoData = True
492 492 self.flagDiscontinuousBlock = False
493 493 self.utctime = None
494 494 self.nCohInt = None
495 495 self.nIncohInt = None
496 496 self.blocksize = None
497 497 self.nFFTPoints = None
498 498 self.wavelength = None
499 499 self.flagDecodeData = False # asumo q la data no esta decodificada
500 500 self.flagDeflipData = False # asumo q la data no esta sin flip
501 501 self.flagShiftFFT = False
502 502 self.ippFactor = 1
503 503 #self.noise = None
504 504 self.beacon_heiIndexList = []
505 505 self.noise_estimation = None
506 506
507 507 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
508 508 """
509 509 Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
510 510
511 511 Return:
512 512 noiselevel
513 513 """
514 514
515 515 noise = numpy.zeros(self.nChannels)
516 516
517 517 for channel in range(self.nChannels):
518 518 daux = self.data_spc[channel,
519 519 xmin_index:xmax_index, ymin_index:ymax_index]
520 520 noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
521 521
522 522 return noise
523 523
524 524 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
525 525
526 526 if self.noise_estimation is not None:
527 527 # this was estimated by getNoise Operation defined in jroproc_spectra.py
528 528 return self.noise_estimation
529 529 else:
530 530 noise = self.getNoisebyHildebrand(
531 531 xmin_index, xmax_index, ymin_index, ymax_index)
532 532 return noise
533 533
534 534 def getFreqRangeTimeResponse(self, extrapoints=0):
535 535
536 536 deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
537 537 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
538 538
539 539 return freqrange
540 540
541 541 def getAcfRange(self, extrapoints=0):
542 542
543 543 deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
544 544 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
545 545
546 546 return freqrange
547 547
548 548 def getFreqRange(self, extrapoints=0):
549 549
550 550 deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
551 551 freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
552 552
553 553 return freqrange
554 554
555 555 def getVelRange(self, extrapoints=0):
556 556
557 557 deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
558 558 velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
559 559
560 560 if self.nmodes:
561 561 return velrange/self.nmodes
562 562 else:
563 563 return velrange
564 564
565 565 def getNPairs(self):
566 566
567 567 return len(self.pairsList)
568 568
569 569 def getPairsIndexList(self):
570 570
571 571 return list(range(self.nPairs))
572 572
573 573 def getNormFactor(self):
574 574
575 575 pwcode = 1
576 576
577 577 if self.flagDecodeData:
578 578 pwcode = numpy.sum(self.code[0]**2)
579 579 #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
580 580 normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
581 581
582 582 return normFactor
583 583
584 584 def getFlagCspc(self):
585 585
586 586 if self.data_cspc is None:
587 587 return True
588 588
589 589 return False
590 590
591 591 def getFlagDc(self):
592 592
593 593 if self.data_dc is None:
594 594 return True
595 595
596 596 return False
597 597
598 598 def getTimeInterval(self):
599 599
600 600 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
601 601 if self.nmodes:
602 602 return self.nmodes*timeInterval
603 603 else:
604 604 return timeInterval
605 605
606 606 def getPower(self):
607 607
608 608 factor = self.normFactor
609 609 z = self.data_spc / factor
610 610 z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
611 611 avg = numpy.average(z, axis=1)
612 612
613 613 return 10 * numpy.log10(avg)
614 614
615 615 def getCoherence(self, pairsList=None, phase=False):
616 616
617 617 z = []
618 618 if pairsList is None:
619 619 pairsIndexList = self.pairsIndexList
620 620 else:
621 621 pairsIndexList = []
622 622 for pair in pairsList:
623 623 if pair not in self.pairsList:
624 624 raise ValueError("Pair %s is not in dataOut.pairsList" % (
625 625 pair))
626 626 pairsIndexList.append(self.pairsList.index(pair))
627 627 for i in range(len(pairsIndexList)):
628 628 pair = self.pairsList[pairsIndexList[i]]
629 629 ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
630 630 powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
631 631 powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
632 632 avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
633 633 if phase:
634 634 data = numpy.arctan2(avgcoherenceComplex.imag,
635 635 avgcoherenceComplex.real) * 180 / numpy.pi
636 636 else:
637 637 data = numpy.abs(avgcoherenceComplex)
638 638
639 639 z.append(data)
640 640
641 641 return numpy.array(z)
642 642
643 643 def setValue(self, value):
644 644
645 645 print("This property should not be initialized")
646 646
647 647 return
648 648
649 649 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
650 650 pairsIndexList = property(
651 651 getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
652 652 normFactor = property(getNormFactor, setValue,
653 653 "I'm the 'getNormFactor' property.")
654 654 flag_cspc = property(getFlagCspc, setValue)
655 655 flag_dc = property(getFlagDc, setValue)
656 656 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
657 657 timeInterval = property(getTimeInterval, setValue,
658 658 "I'm the 'timeInterval' property")
659 659
660 660
661 661 class SpectraHeis(Spectra):
662 662
663 663 data_spc = None
664 664 data_cspc = None
665 665 data_dc = None
666 666 nFFTPoints = None
667 667 # nPairs = None
668 668 pairsList = None
669 669 nCohInt = None
670 670 nIncohInt = None
671 671
672 672 def __init__(self):
673 673
674 674 self.radarControllerHeaderObj = RadarControllerHeader()
675 675
676 676 self.systemHeaderObj = SystemHeader()
677 677
678 678 self.type = "SpectraHeis"
679 679
680 680 # self.dtype = None
681 681
682 682 # self.nChannels = 0
683 683
684 684 # self.nHeights = 0
685 685
686 686 self.nProfiles = None
687 687
688 688 self.heightList = None
689 689
690 690 self.channelList = None
691 691
692 692 # self.channelIndexList = None
693 693
694 694 self.flagNoData = True
695 695
696 696 self.flagDiscontinuousBlock = False
697 697
698 698 # self.nPairs = 0
699 699
700 700 self.utctime = None
701 701
702 702 self.blocksize = None
703 703
704 704 self.profileIndex = 0
705 705
706 706 self.nCohInt = 1
707 707
708 708 self.nIncohInt = 1
709 709
710 710 def getNormFactor(self):
711 711 pwcode = 1
712 712 if self.flagDecodeData:
713 713 pwcode = numpy.sum(self.code[0]**2)
714 714
715 715 normFactor = self.nIncohInt * self.nCohInt * pwcode
716 716
717 717 return normFactor
718 718
719 719 def getTimeInterval(self):
720 720
721 721 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
722 722
723 723 return timeInterval
724 724
725 725 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
726 726 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
727 727
728 728
729 729 class Fits(JROData):
730 730
731 731 heightList = None
732 732 channelList = None
733 733 flagNoData = True
734 734 flagDiscontinuousBlock = False
735 735 useLocalTime = False
736 736 utctime = None
737 737 timeZone = None
738 738 # ippSeconds = None
739 739 # timeInterval = None
740 740 nCohInt = None
741 741 nIncohInt = None
742 742 noise = None
743 743 windowOfFilter = 1
744 744 # Speed of ligth
745 745 C = 3e8
746 746 frequency = 49.92e6
747 747 realtime = False
748 748
749 749 def __init__(self):
750 750
751 751 self.type = "Fits"
752 752
753 753 self.nProfiles = None
754 754
755 755 self.heightList = None
756 756
757 757 self.channelList = None
758 758
759 759 # self.channelIndexList = None
760 760
761 761 self.flagNoData = True
762 762
763 763 self.utctime = None
764 764
765 765 self.nCohInt = 1
766 766
767 767 self.nIncohInt = 1
768 768
769 769 self.useLocalTime = True
770 770
771 771 self.profileIndex = 0
772 772
773 773 # self.utctime = None
774 774 # self.timeZone = None
775 775 # self.ltctime = None
776 776 # self.timeInterval = None
777 777 # self.header = None
778 778 # self.data_header = None
779 779 # self.data = None
780 780 # self.datatime = None
781 781 # self.flagNoData = False
782 782 # self.expName = ''
783 783 # self.nChannels = None
784 784 # self.nSamples = None
785 785 # self.dataBlocksPerFile = None
786 786 # self.comments = ''
787 787 #
788 788
789 789 def getltctime(self):
790 790
791 791 if self.useLocalTime:
792 792 return self.utctime - self.timeZone * 60
793 793
794 794 return self.utctime
795 795
796 796 def getDatatime(self):
797 797
798 798 datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
799 799 return datatime
800 800
801 801 def getTimeRange(self):
802 802
803 803 datatime = []
804 804
805 805 datatime.append(self.ltctime)
806 806 datatime.append(self.ltctime + self.timeInterval)
807 807
808 808 datatime = numpy.array(datatime)
809 809
810 810 return datatime
811 811
812 812 def getHeiRange(self):
813 813
814 814 heis = self.heightList
815 815
816 816 return heis
817 817
818 818 def getNHeights(self):
819 819
820 820 return len(self.heightList)
821 821
822 822 def getNChannels(self):
823 823
824 824 return len(self.channelList)
825 825
826 826 def getChannelIndexList(self):
827 827
828 828 return list(range(self.nChannels))
829 829
830 830 def getNoise(self, type=1):
831 831
832 832 #noise = numpy.zeros(self.nChannels)
833 833
834 834 if type == 1:
835 835 noise = self.getNoisebyHildebrand()
836 836
837 837 if type == 2:
838 838 noise = self.getNoisebySort()
839 839
840 840 if type == 3:
841 841 noise = self.getNoisebyWindow()
842 842
843 843 return noise
844 844
845 845 def getTimeInterval(self):
846 846
847 847 timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
848 848
849 849 return timeInterval
850 850
851 851 def get_ippSeconds(self):
852 852 '''
853 853 '''
854 854 return self.ipp_sec
855 855
856 856
857 857 datatime = property(getDatatime, "I'm the 'datatime' property")
858 858 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
859 859 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
860 860 channelIndexList = property(
861 861 getChannelIndexList, "I'm the 'channelIndexList' property.")
862 862 noise = property(getNoise, "I'm the 'nHeights' property.")
863 863
864 864 ltctime = property(getltctime, "I'm the 'ltctime' property")
865 865 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
866 866 ippSeconds = property(get_ippSeconds, '')
867 867
868 868 class Correlation(JROData):
869 869
870 870 noise = None
871 871 SNR = None
872 872 #--------------------------------------------------
873 873 mode = None
874 874 split = False
875 875 data_cf = None
876 876 lags = None
877 877 lagRange = None
878 878 pairsList = None
879 879 normFactor = None
880 880 #--------------------------------------------------
881 881 # calculateVelocity = None
882 882 nLags = None
883 883 nPairs = None
884 884 nAvg = None
885 885
886 886 def __init__(self):
887 887 '''
888 888 Constructor
889 889 '''
890 890 self.radarControllerHeaderObj = RadarControllerHeader()
891 891
892 892 self.systemHeaderObj = SystemHeader()
893 893
894 894 self.type = "Correlation"
895 895
896 896 self.data = None
897 897
898 898 self.dtype = None
899 899
900 900 self.nProfiles = None
901 901
902 902 self.heightList = None
903 903
904 904 self.channelList = None
905 905
906 906 self.flagNoData = True
907 907
908 908 self.flagDiscontinuousBlock = False
909 909
910 910 self.utctime = None
911 911
912 912 self.timeZone = None
913 913
914 914 self.dstFlag = None
915 915
916 916 self.errorCount = None
917 917
918 918 self.blocksize = None
919 919
920 920 self.flagDecodeData = False # asumo q la data no esta decodificada
921 921
922 922 self.flagDeflipData = False # asumo q la data no esta sin flip
923 923
924 924 self.pairsList = None
925 925
926 926 self.nPoints = None
927 927
928 928 def getPairsList(self):
929 929
930 930 return self.pairsList
931 931
932 932 def getNoise(self, mode=2):
933 933
934 934 indR = numpy.where(self.lagR == 0)[0][0]
935 935 indT = numpy.where(self.lagT == 0)[0][0]
936 936
937 937 jspectra0 = self.data_corr[:, :, indR, :]
938 938 jspectra = copy.copy(jspectra0)
939 939
940 940 num_chan = jspectra.shape[0]
941 941 num_hei = jspectra.shape[2]
942 942
943 943 freq_dc = jspectra.shape[1] / 2
944 944 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
945 945
946 946 if ind_vel[0] < 0:
947 947 ind_vel[list(range(0, 1))] = ind_vel[list(
948 948 range(0, 1))] + self.num_prof
949 949
950 950 if mode == 1:
951 951 jspectra[:, freq_dc, :] = (
952 952 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
953 953
954 954 if mode == 2:
955 955
956 956 vel = numpy.array([-2, -1, 1, 2])
957 957 xx = numpy.zeros([4, 4])
958 958
959 959 for fil in range(4):
960 960 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
961 961
962 962 xx_inv = numpy.linalg.inv(xx)
963 963 xx_aux = xx_inv[0, :]
964 964
965 965 for ich in range(num_chan):
966 966 yy = jspectra[ich, ind_vel, :]
967 967 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
968 968
969 969 junkid = jspectra[ich, freq_dc, :] <= 0
970 970 cjunkid = sum(junkid)
971 971
972 972 if cjunkid.any():
973 973 jspectra[ich, freq_dc, junkid.nonzero()] = (
974 974 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
975 975
976 976 noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
977 977
978 978 return noise
979 979
980 980 def getTimeInterval(self):
981 981
982 982 timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
983 983
984 984 return timeInterval
985 985
986 986 def splitFunctions(self):
987 987
988 988 pairsList = self.pairsList
989 989 ccf_pairs = []
990 990 acf_pairs = []
991 991 ccf_ind = []
992 992 acf_ind = []
993 993 for l in range(len(pairsList)):
994 994 chan0 = pairsList[l][0]
995 995 chan1 = pairsList[l][1]
996 996
997 997 # Obteniendo pares de Autocorrelacion
998 998 if chan0 == chan1:
999 999 acf_pairs.append(chan0)
1000 1000 acf_ind.append(l)
1001 1001 else:
1002 1002 ccf_pairs.append(pairsList[l])
1003 1003 ccf_ind.append(l)
1004 1004
1005 1005 data_acf = self.data_cf[acf_ind]
1006 1006 data_ccf = self.data_cf[ccf_ind]
1007 1007
1008 1008 return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
1009 1009
1010 1010 def getNormFactor(self):
1011 1011 acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
1012 1012 acf_pairs = numpy.array(acf_pairs)
1013 1013 normFactor = numpy.zeros((self.nPairs, self.nHeights))
1014 1014
1015 1015 for p in range(self.nPairs):
1016 1016 pair = self.pairsList[p]
1017 1017
1018 1018 ch0 = pair[0]
1019 1019 ch1 = pair[1]
1020 1020
1021 1021 ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
1022 1022 ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
1023 1023 normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
1024 1024
1025 1025 return normFactor
1026 1026
1027 1027 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
1028 1028 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
1029 1029
1030 1030
1031 1031 class Parameters(Spectra):
1032 1032
1033 1033 experimentInfo = None # Information about the experiment
1034 1034 # Information from previous data
1035 1035 inputUnit = None # Type of data to be processed
1036 1036 operation = None # Type of operation to parametrize
1037 1037 # normFactor = None #Normalization Factor
1038 1038 groupList = None # List of Pairs, Groups, etc
1039 1039 # Parameters
1040 1040 data_param = None # Parameters obtained
1041 1041 data_pre = None # Data Pre Parametrization
1042 1042 data_SNR = None # Signal to Noise Ratio
1043 1043 # heightRange = None #Heights
1044 1044 abscissaList = None # Abscissa, can be velocities, lags or time
1045 1045 # noise = None #Noise Potency
1046 1046 utctimeInit = None # Initial UTC time
1047 1047 paramInterval = None # Time interval to calculate Parameters in seconds
1048 1048 useLocalTime = True
1049 1049 # Fitting
1050 1050 data_error = None # Error of the estimation
1051 1051 constants = None
1052 1052 library = None
1053 1053 # Output signal
1054 1054 outputInterval = None # Time interval to calculate output signal in seconds
1055 1055 data_output = None # Out signal
1056 1056 nAvg = None
1057 1057 noise_estimation = None
1058 1058 GauSPC = None # Fit gaussian SPC
1059 1059
1060 1060 def __init__(self):
1061 1061 '''
1062 1062 Constructor
1063 1063 '''
1064 1064 self.radarControllerHeaderObj = RadarControllerHeader()
1065 1065
1066 1066 self.systemHeaderObj = SystemHeader()
1067 1067
1068 1068 self.type = "Parameters"
1069 1069
1070 1070 def getTimeRange1(self, interval):
1071 1071
1072 1072 datatime = []
1073 1073
1074 1074 if self.useLocalTime:
1075 1075 time1 = self.utctimeInit - self.timeZone * 60
1076 1076 else:
1077 1077 time1 = self.utctimeInit
1078 1078
1079 1079 datatime.append(time1)
1080 1080 datatime.append(time1 + interval)
1081 1081 datatime = numpy.array(datatime)
1082 1082
1083 1083 return datatime
1084 1084
1085 1085 def getTimeInterval(self):
1086 1086
1087 1087 if hasattr(self, 'timeInterval1'):
1088 1088 return self.timeInterval1
1089 1089 else:
1090 1090 return self.paramInterval
1091 1091
1092 1092 def setValue(self, value):
1093 1093
1094 1094 print("This property should not be initialized")
1095 1095
1096 1096 return
1097 1097
1098 1098 def getNoise(self):
1099 1099
1100 1100 return self.spc_noise
1101 1101
1102 1102 timeInterval = property(getTimeInterval)
1103 1103 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
1104 1104
1105 1105
1106 1106 class PlotterData(object):
1107 1107 '''
1108 1108 Object to hold data to be plotted
1109 1109 '''
1110 1110
1111 1111 MAXNUMX = 200
1112 1112 MAXNUMY = 200
1113 1113
1114 1114 def __init__(self, code, throttle_value, exp_code, buffering=True, snr=False):
1115 1115
1116 1116 self.key = code
1117 1117 self.throttle = throttle_value
1118 1118 self.exp_code = exp_code
1119 1119 self.buffering = buffering
1120 1120 self.ready = False
1121 1121 self.flagNoData = False
1122 1122 self.localtime = False
1123 1123 self.data = {}
1124 1124 self.meta = {}
1125 self.__times = []
1126 1125 self.__heights = []
1127 1126
1128 1127 if 'snr' in code:
1129 1128 self.plottypes = ['snr']
1130 1129 elif code == 'spc':
1131 1130 self.plottypes = ['spc', 'noise', 'rti']
1131 elif code == 'cspc':
1132 self.plottypes = ['cspc', 'spc', 'noise', 'rti']
1132 1133 elif code == 'rti':
1133 1134 self.plottypes = ['noise', 'rti']
1134 1135 else:
1135 1136 self.plottypes = [code]
1136 1137
1137 1138 if 'snr' not in self.plottypes and snr:
1138 1139 self.plottypes.append('snr')
1139 1140
1140 1141 for plot in self.plottypes:
1141 1142 self.data[plot] = {}
1142 1143
1143 1144
1144 1145 def __str__(self):
1145 1146 dum = ['{}{}'.format(key, self.shape(key)) for key in self.data]
1146 return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times))
1147 return 'Data[{}][{}]'.format(';'.join(dum), len(self.times))
1147 1148
1148 1149 def __len__(self):
1149 return len(self.__times)
1150 return len(self.data[self.key])
1150 1151
1151 1152 def __getitem__(self, key):
1152 1153
1153 1154 if key not in self.data:
1154 1155 raise KeyError(log.error('Missing key: {}'.format(key)))
1155 1156 if 'spc' in key or not self.buffering:
1156 ret = self.data[key]
1157 ret = self.data[key][self.tm]
1157 1158 elif 'scope' in key:
1158 1159 ret = numpy.array(self.data[key][float(self.tm)])
1159 1160 else:
1160 1161 ret = numpy.array([self.data[key][x] for x in self.times])
1161 1162 if ret.ndim > 1:
1162 1163 ret = numpy.swapaxes(ret, 0, 1)
1163 1164 return ret
1164 1165
1165 1166 def __contains__(self, key):
1166 1167 return key in self.data
1167 1168
1168 1169 def setup(self):
1169 1170 '''
1170 1171 Configure object
1171 1172 '''
1172 1173 self.type = ''
1173 1174 self.ready = False
1175 del self.data
1174 1176 self.data = {}
1175 self.__times = []
1176 1177 self.__heights = []
1177 1178 self.__all_heights = set()
1178 1179 for plot in self.plottypes:
1179 1180 if 'snr' in plot:
1180 1181 plot = 'snr'
1181 1182 elif 'spc_moments' == plot:
1182 1183 plot = 'moments'
1183 1184 self.data[plot] = {}
1184 1185
1185 1186 if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data:
1186 1187 self.data['noise'] = {}
1187 1188 self.data['rti'] = {}
1188 1189 if 'noise' not in self.plottypes:
1189 1190 self.plottypes.append('noise')
1190 1191 if 'rti' not in self.plottypes:
1191 1192 self.plottypes.append('rti')
1192 1193
1193 1194 def shape(self, key):
1194 1195 '''
1195 1196 Get the shape of the one-element data for the given key
1196 1197 '''
1197 1198
1198 1199 if len(self.data[key]):
1199 1200 if 'spc' in key or not self.buffering:
1200 1201 return self.data[key].shape
1201 return self.data[key][self.__times[0]].shape
1202 return self.data[key][self.times[0]].shape
1202 1203 return (0,)
1203 1204
1204 1205 def update(self, dataOut, tm):
1205 1206 '''
1206 1207 Update data object with new dataOut
1207 1208 '''
1208 if tm in self.__times:
1209 return
1209
1210 1210 self.profileIndex = dataOut.profileIndex
1211 1211 self.tm = tm
1212 1212 self.type = dataOut.type
1213 1213 self.parameters = getattr(dataOut, 'parameters', [])
1214 1214
1215 1215 if hasattr(dataOut, 'meta'):
1216 1216 self.meta.update(dataOut.meta)
1217 1217
1218 1218 if hasattr(dataOut, 'pairsList'):
1219 1219 self.pairs = dataOut.pairsList
1220 1220
1221 1221 self.interval = dataOut.getTimeInterval()
1222 1222 self.localtime = dataOut.useLocalTime
1223 1223 if True in ['spc' in ptype for ptype in self.plottypes]:
1224 1224 self.xrange = (dataOut.getFreqRange(1)/1000.,
1225 1225 dataOut.getAcfRange(1), dataOut.getVelRange(1))
1226 self.factor = dataOut.normFactor
1227 1226 self.__heights.append(dataOut.heightList)
1228 1227 self.__all_heights.update(dataOut.heightList)
1229 self.__times.append(tm)
1228
1230 1229 for plot in self.plottypes:
1231 1230 if plot in ('spc', 'spc_moments', 'spc_cut'):
1232 1231 z = dataOut.data_spc/dataOut.normFactor
1233 1232 buffer = 10*numpy.log10(z)
1234 1233 if plot == 'cspc':
1235 z = dataOut.data_spc/dataOut.normFactor
1236 1234 buffer = (dataOut.data_spc, dataOut.data_cspc)
1237 1235 if plot == 'noise':
1238 1236 buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
1239 1237 if plot in ('rti', 'spcprofile'):
1240 1238 buffer = dataOut.getPower()
1241 1239 if plot == 'snr_db':
1242 1240 buffer = dataOut.data_SNR
1243 1241 if plot == 'snr':
1244 1242 buffer = 10*numpy.log10(dataOut.data_SNR)
1245 1243 if plot == 'dop':
1246 1244 buffer = dataOut.data_DOP
1247 1245 if plot == 'pow':
1248 1246 buffer = 10*numpy.log10(dataOut.data_POW)
1249 1247 if plot == 'width':
1250 1248 buffer = dataOut.data_WIDTH
1251 1249 if plot == 'coh':
1252 1250 buffer = dataOut.getCoherence()
1253 1251 if plot == 'phase':
1254 1252 buffer = dataOut.getCoherence(phase=True)
1255 1253 if plot == 'output':
1256 1254 buffer = dataOut.data_output
1257 1255 if plot == 'param':
1258 1256 buffer = dataOut.data_param
1259 1257 if plot == 'scope':
1260 1258 buffer = dataOut.data
1261 1259 self.flagDataAsBlock = dataOut.flagDataAsBlock
1262 1260 self.nProfiles = dataOut.nProfiles
1263 1261 if plot == 'pp_power':
1264 1262 buffer = dataOut.data_intensity
1265 1263 self.flagDataAsBlock = dataOut.flagDataAsBlock
1266 1264 self.nProfiles = dataOut.nProfiles
1267 1265 if plot == 'pp_velocity':
1268 1266 buffer = dataOut.data_velocity
1269 1267 self.flagDataAsBlock = dataOut.flagDataAsBlock
1270 1268 self.nProfiles = dataOut.nProfiles
1271 1269
1272 1270 if plot == 'spc':
1273 self.data['spc'] = buffer
1271 self.data['spc'][tm] = buffer
1274 1272 elif plot == 'cspc':
1275 self.data['spc'] = buffer[0]
1276 self.data['cspc'] = buffer[1]
1273 self.data['cspc'][tm] = buffer
1277 1274 elif plot == 'spc_moments':
1278 self.data['spc'] = buffer
1275 self.data['spc'][tm] = buffer
1279 1276 self.data['moments'][tm] = dataOut.moments
1280 1277 else:
1281 1278 if self.buffering:
1282 1279 self.data[plot][tm] = buffer
1283 1280 else:
1284 self.data[plot] = buffer
1281 self.data[plot][tm] = buffer
1285 1282
1286 1283 if dataOut.channelList is None:
1287 1284 self.channels = range(buffer.shape[0])
1288 1285 else:
1289 1286 self.channels = dataOut.channelList
1290 1287
1291 1288 if buffer is None:
1292 1289 self.flagNoData = True
1293 1290 raise schainpy.admin.SchainWarning('Attribute data_{} is empty'.format(self.key))
1294 1291
1295 1292 def normalize_heights(self):
1296 1293 '''
1297 1294 Ensure same-dimension of the data for different heighList
1298 1295 '''
1299 1296
1300 1297 H = numpy.array(list(self.__all_heights))
1301 1298 H.sort()
1302 1299 for key in self.data:
1303 1300 shape = self.shape(key)[:-1] + H.shape
1304 1301 for tm, obj in list(self.data[key].items()):
1305 h = self.__heights[self.__times.index(tm)]
1302 h = self.__heights[self.times.index(tm)]
1306 1303 if H.size == h.size:
1307 1304 continue
1308 1305 index = numpy.where(numpy.in1d(H, h))[0]
1309 1306 dummy = numpy.zeros(shape) + numpy.nan
1310 1307 if len(shape) == 2:
1311 1308 dummy[:, index] = obj
1312 1309 else:
1313 1310 dummy[index] = obj
1314 1311 self.data[key][tm] = dummy
1315 1312
1316 self.__heights = [H for tm in self.__times]
1313 self.__heights = [H for tm in self.times]
1317 1314
1318 def jsonify(self, plot_name, plot_type, decimate=False):
1315 def jsonify(self, tm, plot_name, plot_type, decimate=False):
1319 1316 '''
1320 1317 Convert data to json
1321 1318 '''
1322 1319
1323 tm = self.times[-1]
1324 1320 dy = int(self.heights.size/self.MAXNUMY) + 1
1325 if self.key in ('spc', 'cspc') or not self.buffering:
1326 dx = int(self.data[self.key].shape[1]/self.MAXNUMX) + 1
1321 if self.key in ('spc', 'cspc'):
1322 dx = int(self.data[self.key][tm].shape[1]/self.MAXNUMX) + 1
1327 1323 data = self.roundFloats(
1328 self.data[self.key][::, ::dx, ::dy].tolist())
1324 self.data[self.key][tm][::, ::dx, ::dy].tolist())
1329 1325 else:
1330 1326 if self.key is 'noise':
1331 1327 data = [[x] for x in self.roundFloats(self.data[self.key][tm].tolist())]
1332 1328 else:
1333 1329 data = self.roundFloats(self.data[self.key][tm][::, ::dy].tolist())
1334 1330
1335 1331 meta = {}
1336 1332 ret = {
1337 1333 'plot': plot_name,
1338 1334 'code': self.exp_code,
1339 1335 'time': float(tm),
1340 1336 'data': data,
1341 1337 }
1342 1338 meta['type'] = plot_type
1343 1339 meta['interval'] = float(self.interval)
1344 1340 meta['localtime'] = self.localtime
1345 1341 meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
1346 1342 if 'spc' in self.data or 'cspc' in self.data:
1347 1343 meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
1348 1344 else:
1349 1345 meta['xrange'] = []
1350 1346
1351 1347 meta.update(self.meta)
1352 1348 ret['metadata'] = meta
1353 1349 return json.dumps(ret)
1354 1350
1355 1351 @property
1356 1352 def times(self):
1357 1353 '''
1358 1354 Return the list of times of the current data
1359 1355 '''
1360 1356
1361 ret = numpy.array(self.__times)
1362 ret.sort()
1357 ret = numpy.array([*self.data[self.key]])
1358 if self:
1359 ret.sort()
1363 1360 return ret
1364 1361
1365 1362 @property
1366 1363 def min_time(self):
1367 1364 '''
1368 1365 Return the minimun time value
1369 1366 '''
1370 1367
1371 1368 return self.times[0]
1372 1369
1373 1370 @property
1374 1371 def max_time(self):
1375 1372 '''
1376 1373 Return the maximun time value
1377 1374 '''
1378 1375
1379 1376 return self.times[-1]
1380 1377
1381 1378 @property
1382 1379 def heights(self):
1383 1380 '''
1384 1381 Return the list of heights of the current data
1385 1382 '''
1386 1383
1387 1384 return numpy.array(self.__heights[-1])
1388 1385
1389 1386 @staticmethod
1390 1387 def roundFloats(obj):
1391 1388 if isinstance(obj, list):
1392 1389 return list(map(PlotterData.roundFloats, obj))
1393 1390 elif isinstance(obj, float):
1394 1391 return round(obj, 2)
@@ -1,705 +1,717
1 1
2 2 import os
3 3 import sys
4 4 import zmq
5 5 import time
6 6 import numpy
7 7 import datetime
8 8 from queue import Queue
9 9 from functools import wraps
10 10 from threading import Thread
11 11 import matplotlib
12 12
13 13 if 'BACKEND' in os.environ:
14 14 matplotlib.use(os.environ['BACKEND'])
15 15 elif 'linux' in sys.platform:
16 16 matplotlib.use("TkAgg")
17 17 elif 'darwin' in sys.platform:
18 18 matplotlib.use('WxAgg')
19 19 else:
20 20 from schainpy.utils import log
21 21 log.warning('Using default Backend="Agg"', 'INFO')
22 22 matplotlib.use('Agg')
23 23
24 24 import matplotlib.pyplot as plt
25 25 from matplotlib.patches import Polygon
26 26 from mpl_toolkits.axes_grid1 import make_axes_locatable
27 27 from matplotlib.ticker import FuncFormatter, LinearLocator, MultipleLocator
28 28
29 29 from schainpy.model.data.jrodata import PlotterData
30 30 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
31 31 from schainpy.utils import log
32 32
33 33 jet_values = matplotlib.pyplot.get_cmap('jet', 100)(numpy.arange(100))[10:90]
34 34 blu_values = matplotlib.pyplot.get_cmap(
35 35 'seismic_r', 20)(numpy.arange(20))[10:15]
36 36 ncmap = matplotlib.colors.LinearSegmentedColormap.from_list(
37 37 'jro', numpy.vstack((blu_values, jet_values)))
38 38 matplotlib.pyplot.register_cmap(cmap=ncmap)
39 39
40 40 CMAPS = [plt.get_cmap(s) for s in ('jro', 'jet', 'viridis',
41 41 'plasma', 'inferno', 'Greys', 'seismic', 'bwr', 'coolwarm')]
42 42
43 43 EARTH_RADIUS = 6.3710e3
44 44
45 45 def ll2xy(lat1, lon1, lat2, lon2):
46 46
47 47 p = 0.017453292519943295
48 48 a = 0.5 - numpy.cos((lat2 - lat1) * p)/2 + numpy.cos(lat1 * p) * \
49 49 numpy.cos(lat2 * p) * (1 - numpy.cos((lon2 - lon1) * p)) / 2
50 50 r = 12742 * numpy.arcsin(numpy.sqrt(a))
51 51 theta = numpy.arctan2(numpy.sin((lon2-lon1)*p)*numpy.cos(lat2*p), numpy.cos(lat1*p)
52 52 * numpy.sin(lat2*p)-numpy.sin(lat1*p)*numpy.cos(lat2*p)*numpy.cos((lon2-lon1)*p))
53 53 theta = -theta + numpy.pi/2
54 54 return r*numpy.cos(theta), r*numpy.sin(theta)
55 55
56 56
57 57 def km2deg(km):
58 58 '''
59 59 Convert distance in km to degrees
60 60 '''
61 61
62 62 return numpy.rad2deg(km/EARTH_RADIUS)
63 63
64 64
65 65 def figpause(interval):
66 66 backend = plt.rcParams['backend']
67 67 if backend in matplotlib.rcsetup.interactive_bk:
68 68 figManager = matplotlib._pylab_helpers.Gcf.get_active()
69 69 if figManager is not None:
70 70 canvas = figManager.canvas
71 71 if canvas.figure.stale:
72 72 canvas.draw()
73 73 try:
74 74 canvas.start_event_loop(interval)
75 75 except:
76 76 pass
77 77 return
78 78
79 79
80 80 def popup(message):
81 81 '''
82 82 '''
83 83
84 84 fig = plt.figure(figsize=(12, 8), facecolor='r')
85 85 text = '\n'.join([s.strip() for s in message.split(':')])
86 86 fig.text(0.01, 0.5, text, ha='left', va='center',
87 87 size='20', weight='heavy', color='w')
88 88 fig.show()
89 89 figpause(1000)
90 90
91 91
92 92 class Throttle(object):
93 93 '''
94 94 Decorator that prevents a function from being called more than once every
95 95 time period.
96 96 To create a function that cannot be called more than once a minute, but
97 97 will sleep until it can be called:
98 98 @Throttle(minutes=1)
99 99 def foo():
100 100 pass
101 101
102 102 for i in range(10):
103 103 foo()
104 104 print "This function has run %s times." % i
105 105 '''
106 106
107 107 def __init__(self, seconds=0, minutes=0, hours=0):
108 108 self.throttle_period = datetime.timedelta(
109 109 seconds=seconds, minutes=minutes, hours=hours
110 110 )
111 111
112 112 self.time_of_last_call = datetime.datetime.min
113 113
114 114 def __call__(self, fn):
115 115 @wraps(fn)
116 116 def wrapper(*args, **kwargs):
117 117 coerce = kwargs.pop('coerce', None)
118 118 if coerce:
119 119 self.time_of_last_call = datetime.datetime.now()
120 120 return fn(*args, **kwargs)
121 121 else:
122 122 now = datetime.datetime.now()
123 123 time_since_last_call = now - self.time_of_last_call
124 124 time_left = self.throttle_period - time_since_last_call
125 125
126 126 if time_left > datetime.timedelta(seconds=0):
127 127 return
128 128
129 129 self.time_of_last_call = datetime.datetime.now()
130 130 return fn(*args, **kwargs)
131 131
132 132 return wrapper
133 133
134 134 def apply_throttle(value):
135 135
136 136 @Throttle(seconds=value)
137 137 def fnThrottled(fn):
138 138 fn()
139 139
140 140 return fnThrottled
141 141
142 142
143 143 @MPDecorator
144 144 class Plot(Operation):
145 145 '''
146 146 Base class for Schain plotting operations
147 147 '''
148 148
149 149 CODE = 'Figure'
150 150 colormap = 'jet'
151 151 bgcolor = 'white'
152 152 buffering = True
153 153 __missing = 1E30
154 154
155 155 __attrs__ = ['show', 'save', 'ymin', 'ymax', 'zmin', 'zmax', 'title',
156 156 'showprofile']
157 157
158 158 def __init__(self):
159 159
160 160 Operation.__init__(self)
161 161 self.isConfig = False
162 162 self.isPlotConfig = False
163 163 self.save_counter = 1
164 164 self.sender_time = 0
165 165 self.data = None
166 166 self.firsttime = True
167 167 self.sender_queue = Queue(maxsize=10)
168 168 self.plots_adjust = {'left': 0.125, 'right': 0.9, 'bottom': 0.15, 'top': 0.9, 'wspace': 0.2, 'hspace': 0.2}
169 169
170 170 def __fmtTime(self, x, pos):
171 171 '''
172 172 '''
173 173
174 174 return '{}'.format(self.getDateTime(x).strftime('%H:%M'))
175 175
176 176 def __setup(self, **kwargs):
177 177 '''
178 178 Initialize variables
179 179 '''
180 180
181 181 self.figures = []
182 182 self.axes = []
183 183 self.cb_axes = []
184 184 self.localtime = kwargs.pop('localtime', True)
185 185 self.show = kwargs.get('show', True)
186 186 self.save = kwargs.get('save', False)
187 187 self.save_period = kwargs.get('save_period', 1)
188 188 self.ftp = kwargs.get('ftp', False)
189 189 self.colormap = kwargs.get('colormap', self.colormap)
190 190 self.colormap_coh = kwargs.get('colormap_coh', 'jet')
191 191 self.colormap_phase = kwargs.get('colormap_phase', 'RdBu_r')
192 192 self.colormaps = kwargs.get('colormaps', None)
193 193 self.bgcolor = kwargs.get('bgcolor', self.bgcolor)
194 194 self.showprofile = kwargs.get('showprofile', False)
195 195 self.title = kwargs.get('wintitle', self.CODE.upper())
196 196 self.cb_label = kwargs.get('cb_label', None)
197 197 self.cb_labels = kwargs.get('cb_labels', None)
198 198 self.labels = kwargs.get('labels', None)
199 199 self.xaxis = kwargs.get('xaxis', 'frequency')
200 200 self.zmin = kwargs.get('zmin', None)
201 201 self.zmax = kwargs.get('zmax', None)
202 202 self.zlimits = kwargs.get('zlimits', None)
203 203 self.xmin = kwargs.get('xmin', None)
204 204 self.xmax = kwargs.get('xmax', None)
205 self.xrange = kwargs.get('xrange', 24)
205 self.xrange = kwargs.get('xrange', 12)
206 206 self.xscale = kwargs.get('xscale', None)
207 207 self.ymin = kwargs.get('ymin', None)
208 208 self.ymax = kwargs.get('ymax', None)
209 209 self.yscale = kwargs.get('yscale', None)
210 210 self.xlabel = kwargs.get('xlabel', None)
211 211 self.attr_time = kwargs.get('attr_time', 'utctime')
212 212 self.decimation = kwargs.get('decimation', None)
213 213 self.showSNR = kwargs.get('showSNR', False)
214 214 self.oneFigure = kwargs.get('oneFigure', True)
215 215 self.width = kwargs.get('width', None)
216 216 self.height = kwargs.get('height', None)
217 217 self.colorbar = kwargs.get('colorbar', True)
218 218 self.factors = kwargs.get('factors', [1, 1, 1, 1, 1, 1, 1, 1])
219 219 self.channels = kwargs.get('channels', None)
220 220 self.titles = kwargs.get('titles', [])
221 221 self.polar = False
222 222 self.type = kwargs.get('type', 'iq')
223 223 self.grid = kwargs.get('grid', False)
224 224 self.pause = kwargs.get('pause', False)
225 225 self.save_code = kwargs.get('save_code', None)
226 226 self.realtime = kwargs.get('realtime', True)
227 227 self.throttle = kwargs.get('throttle', 0)
228 228 self.exp_code = kwargs.get('exp_code', None)
229 229 self.plot_server = kwargs.get('plot_server', False)
230 230 self.sender_period = kwargs.get('sender_period', 60)
231 self.tag = kwargs.get('tag', '')
231 232 self.height_index = kwargs.get('height_index', None)
232 233 self.__throttle_plot = apply_throttle(self.throttle)
233 234 self.data = PlotterData(
234 235 self.CODE, self.throttle, self.exp_code, self.buffering, snr=self.showSNR)
235 236
236 237 if self.plot_server:
237 238 if not self.plot_server.startswith('tcp://'):
238 239 self.plot_server = 'tcp://{}'.format(self.plot_server)
239 240 log.success(
240 241 'Sending to server: {}'.format(self.plot_server),
241 242 self.name
242 243 )
243 244 if 'plot_name' in kwargs:
244 245 self.plot_name = kwargs['plot_name']
245 246
246 247 def __setup_plot(self):
247 248 '''
248 249 Common setup for all figures, here figures and axes are created
249 250 '''
250 251
251 252 self.setup()
252 253
253 254 self.time_label = 'LT' if self.localtime else 'UTC'
254 255
255 256 if self.width is None:
256 257 self.width = 8
257 258
258 259 self.figures = []
259 260 self.axes = []
260 261 self.cb_axes = []
261 262 self.pf_axes = []
262 263 self.cmaps = []
263 264
264 265 size = '15%' if self.ncols == 1 else '30%'
265 266 pad = '4%' if self.ncols == 1 else '8%'
266 267
267 268 if self.oneFigure:
268 269 if self.height is None:
269 270 self.height = 1.4 * self.nrows + 1
270 271 fig = plt.figure(figsize=(self.width, self.height),
271 272 edgecolor='k',
272 273 facecolor='w')
273 274 self.figures.append(fig)
274 275 for n in range(self.nplots):
275 276 ax = fig.add_subplot(self.nrows, self.ncols,
276 277 n + 1, polar=self.polar)
277 278 ax.tick_params(labelsize=8)
278 279 ax.firsttime = True
279 280 ax.index = 0
280 281 ax.press = None
281 282 self.axes.append(ax)
282 283 if self.showprofile:
283 284 cax = self.__add_axes(ax, size=size, pad=pad)
284 285 cax.tick_params(labelsize=8)
285 286 self.pf_axes.append(cax)
286 287 else:
287 288 if self.height is None:
288 289 self.height = 3
289 290 for n in range(self.nplots):
290 291 fig = plt.figure(figsize=(self.width, self.height),
291 292 edgecolor='k',
292 293 facecolor='w')
293 294 ax = fig.add_subplot(1, 1, 1, polar=self.polar)
294 295 ax.tick_params(labelsize=8)
295 296 ax.firsttime = True
296 297 ax.index = 0
297 298 ax.press = None
298 299 self.figures.append(fig)
299 300 self.axes.append(ax)
300 301 if self.showprofile:
301 302 cax = self.__add_axes(ax, size=size, pad=pad)
302 303 cax.tick_params(labelsize=8)
303 304 self.pf_axes.append(cax)
304 305
305 306 for n in range(self.nrows):
306 307 if self.colormaps is not None:
307 308 cmap = plt.get_cmap(self.colormaps[n])
308 309 else:
309 310 cmap = plt.get_cmap(self.colormap)
310 311 cmap.set_bad(self.bgcolor, 1.)
311 312 self.cmaps.append(cmap)
312 313
313 314 def __add_axes(self, ax, size='30%', pad='8%'):
314 315 '''
315 316 Add new axes to the given figure
316 317 '''
317 318 divider = make_axes_locatable(ax)
318 319 nax = divider.new_horizontal(size=size, pad=pad)
319 320 ax.figure.add_axes(nax)
320 321 return nax
321 322
322 323 def fill_gaps(self, x_buffer, y_buffer, z_buffer):
323 324 '''
324 325 Create a masked array for missing data
325 326 '''
326 327 if x_buffer.shape[0] < 2:
327 328 return x_buffer, y_buffer, z_buffer
328 329
329 330 deltas = x_buffer[1:] - x_buffer[0:-1]
330 331 x_median = numpy.median(deltas)
331 332
332 333 index = numpy.where(deltas > 5 * x_median)
333 334
334 335 if len(index[0]) != 0:
335 336 z_buffer[::, index[0], ::] = self.__missing
336 337 z_buffer = numpy.ma.masked_inside(z_buffer,
337 338 0.99 * self.__missing,
338 339 1.01 * self.__missing)
339 340
340 341 return x_buffer, y_buffer, z_buffer
341 342
342 343 def decimate(self):
343 344
344 345 # dx = int(len(self.x)/self.__MAXNUMX) + 1
345 346 dy = int(len(self.y) / self.decimation) + 1
346 347
347 348 # x = self.x[::dx]
348 349 x = self.x
349 350 y = self.y[::dy]
350 351 z = self.z[::, ::, ::dy]
351 352
352 353 return x, y, z
353 354
354 355 def format(self):
355 356 '''
356 357 Set min and max values, labels, ticks and titles
357 358 '''
358 359
359 360 if self.xmin is None:
360 361 xmin = self.data.min_time
361 362 else:
362 363 if self.xaxis is 'time':
363 364 dt = self.getDateTime(self.data.min_time)
364 365 xmin = (dt.replace(hour=int(self.xmin), minute=0, second=0) -
365 366 datetime.datetime(1970, 1, 1)).total_seconds()
366 367 if self.data.localtime:
367 368 xmin += time.timezone
368 369 else:
369 370 xmin = self.xmin
370 371
371 372 if self.xmax is None:
372 373 xmax = xmin + self.xrange * 60 * 60
373 374 else:
374 375 if self.xaxis is 'time':
375 376 dt = self.getDateTime(self.data.max_time)
376 377 xmax = self.xmax - 1
377 378 xmax = (dt.replace(hour=int(xmax), minute=59, second=59) -
378 379 datetime.datetime(1970, 1, 1) + datetime.timedelta(seconds=1)).total_seconds()
379 380 if self.data.localtime:
380 381 xmax += time.timezone
381 382 else:
382 383 xmax = self.xmax
383 384
384 385 ymin = self.ymin if self.ymin else numpy.nanmin(self.y)
385 386 ymax = self.ymax if self.ymax else numpy.nanmax(self.y)
386 387
387 388 for n, ax in enumerate(self.axes):
388 389 if ax.firsttime:
389 390
390 391 dig = int(numpy.log10(ymax))
391 392 if dig == 0:
392 393 digD = len(str(ymax)) - 2
393 394 ydec = ymax*(10**digD)
394 395
395 396 dig = int(numpy.log10(ydec))
396 397 ystep = ((ydec + (10**(dig)))//10**(dig))*(10**(dig))
397 398 ystep = ystep/5
398 399 ystep = ystep/(10**digD)
399 400
400 401 else:
401 402 ystep = ((ymax + (10**(dig)))//10**(dig))*(10**(dig))
402 403 ystep = ystep/5
403 404
404 405 if self.xaxis is not 'time':
405 406
406 407 dig = int(numpy.log10(xmax))
407 408
408 409 if dig <= 0:
409 410 digD = len(str(xmax)) - 2
410 411 xdec = xmax*(10**digD)
411 412
412 413 dig = int(numpy.log10(xdec))
413 414 xstep = ((xdec + (10**(dig)))//10**(dig))*(10**(dig))
414 415 xstep = xstep*0.5
415 416 xstep = xstep/(10**digD)
416 417
417 418 else:
418 419 xstep = ((xmax + (10**(dig)))//10**(dig))*(10**(dig))
419 420 xstep = xstep/5
420 421
421 422 ax.set_facecolor(self.bgcolor)
422 423 ax.yaxis.set_major_locator(MultipleLocator(ystep))
423 424 if self.xscale:
424 425 ax.xaxis.set_major_formatter(FuncFormatter(
425 426 lambda x, pos: '{0:g}'.format(x*self.xscale)))
426 427 if self.xscale:
427 428 ax.yaxis.set_major_formatter(FuncFormatter(
428 429 lambda x, pos: '{0:g}'.format(x*self.yscale)))
429 430 if self.xaxis is 'time':
430 431 ax.xaxis.set_major_formatter(FuncFormatter(self.__fmtTime))
431 432 ax.xaxis.set_major_locator(LinearLocator(9))
432 433 else:
433 434 ax.xaxis.set_major_locator(MultipleLocator(xstep))
434 435 if self.xlabel is not None:
435 436 ax.set_xlabel(self.xlabel)
436 437 ax.set_ylabel(self.ylabel)
437 438 ax.firsttime = False
438 439 if self.showprofile:
439 440 self.pf_axes[n].set_ylim(ymin, ymax)
440 441 self.pf_axes[n].set_xlim(self.zmin, self.zmax)
441 442 self.pf_axes[n].set_xlabel('dB')
442 443 self.pf_axes[n].grid(b=True, axis='x')
443 444 [tick.set_visible(False)
444 445 for tick in self.pf_axes[n].get_yticklabels()]
445 446 if self.colorbar:
446 447 ax.cbar = plt.colorbar(
447 448 ax.plt, ax=ax, fraction=0.05, pad=0.02, aspect=10)
448 449 ax.cbar.ax.tick_params(labelsize=8)
449 450 ax.cbar.ax.press = None
450 451 if self.cb_label:
451 452 ax.cbar.set_label(self.cb_label, size=8)
452 453 elif self.cb_labels:
453 454 ax.cbar.set_label(self.cb_labels[n], size=8)
454 455 else:
455 456 ax.cbar = None
456 457 if self.grid:
457 458 ax.grid(True)
458 459
459 460 if not self.polar:
460 461 ax.set_xlim(xmin, xmax)
461 462 ax.set_ylim(ymin, ymax)
462 463 ax.set_title('{} {} {}'.format(
463 464 self.titles[n],
464 465 self.getDateTime(self.data.max_time).strftime(
465 466 '%Y-%m-%d %H:%M:%S'),
466 467 self.time_label),
467 468 size=8)
468 469 else:
469 470 ax.set_title('{}'.format(self.titles[n]), size=8)
470 471 ax.set_ylim(0, 90)
471 472 ax.set_yticks(numpy.arange(0, 90, 20))
472 473 ax.yaxis.labelpad = 40
473 474
474 475 if self.firsttime:
475 476 for n, fig in enumerate(self.figures):
476 477 fig.subplots_adjust(**self.plots_adjust)
477 478 self.firsttime = False
478 479
479 480 def clear_figures(self):
480 481 '''
481 482 Reset axes for redraw plots
482 483 '''
483 484
484 485 for ax in self.axes+self.pf_axes+self.cb_axes:
485 486 ax.clear()
486 487 ax.firsttime = True
487 488 if ax.cbar:
488 489 ax.cbar.remove()
489 490
490 491 def __plot(self):
491 492 '''
492 493 Main function to plot, format and save figures
493 494 '''
494 495
495 496 self.plot()
496 497 self.format()
497 498
498 499 for n, fig in enumerate(self.figures):
499 500 if self.nrows == 0 or self.nplots == 0:
500 501 log.warning('No data', self.name)
501 502 fig.text(0.5, 0.5, 'No Data', fontsize='large', ha='center')
502 503 fig.canvas.manager.set_window_title(self.CODE)
503 504 continue
504 505
505 506 fig.canvas.manager.set_window_title('{} - {}'.format(self.title,
506 507 self.getDateTime(self.data.max_time).strftime('%Y/%m/%d')))
507 508 fig.canvas.draw()
508 509 if self.show:
509 510 fig.show()
510 511 figpause(0.01)
511 512
512 513 if self.save:
513 514 self.save_figure(n)
514 515
515 516 if self.plot_server:
516 517 self.send_to_server()
517 518
518 519 def save_figure(self, n):
519 520 '''
520 521 '''
521 522
522 523 if self.save_counter < self.save_period:
523 524 self.save_counter += 1
524 525 return
525 526
526 527 self.save_counter = 1
527 528
528 529 fig = self.figures[n]
529 530
530 531 if self.save_code:
531 532 if isinstance(self.save_code, str):
532 533 labels = [self.save_code for x in self.figures]
533 534 else:
534 535 labels = self.save_code
535 536 else:
536 537 labels = [self.CODE for x in self.figures]
537 538
538 539 figname = os.path.join(
539 540 self.save,
540 541 labels[n],
541 542 '{}_{}.png'.format(
542 543 labels[n],
543 544 self.getDateTime(self.data.max_time).strftime(
544 545 '%Y%m%d_%H%M%S'
545 546 ),
546 547 )
547 548 )
548 549 log.log('Saving figure: {}'.format(figname), self.name)
549 550 if not os.path.isdir(os.path.dirname(figname)):
550 551 os.makedirs(os.path.dirname(figname))
551 552 fig.savefig(figname)
552 553
553 554 if self.throttle == 0:
554 555 figname = os.path.join(
555 556 self.save,
556 557 '{}_{}.png'.format(
557 558 labels[n],
558 559 self.getDateTime(self.data.min_time).strftime(
559 560 '%Y%m%d'
560 561 ),
561 562 )
562 563 )
563 564 fig.savefig(figname)
564 565
565 566 def send_to_server(self):
566 567 '''
567 568 '''
568 569
569 570 interval = self.data.tm - self.sender_time
570 571 if interval < self.sender_period:
571 572 return
572 573
573 574 self.sender_time = self.data.tm
574 575
575 attrs = ['titles', 'zmin', 'zmax', 'colormap']
576 attrs = ['titles', 'zmin', 'zmax', 'tag', 'ymin', 'ymax']
576 577 for attr in attrs:
577 578 value = getattr(self, attr)
578 if value is not None:
579 self.data.meta[attr] = getattr(self, attr)
579 if value:
580 if isinstance(value, (numpy.float32, numpy.float64)):
581 value = round(float(value), 2)
582 self.data.meta[attr] = value
583 if self.colormap == 'jet':
584 self.data.meta['colormap'] = 'Jet'
585 elif 'RdBu' in self.colormap:
586 self.data.meta['colormap'] = 'RdBu'
587 else:
588 self.data.meta['colormap'] = 'Viridis'
580 589 self.data.meta['interval'] = int(interval)
581 msg = self.data.jsonify(self.plot_name, self.plot_type)
582 self.sender_queue.put(msg)
590 # msg = self.data.jsonify(self.data.tm, self.plot_name, self.plot_type)
591 self.sender_queue.put(self.data.tm)
583 592
584 593 while True:
585 594 if self.sender_queue.empty():
586 595 break
587 self.socket.send_string(self.sender_queue.get())
596 tm = self.sender_queue.get()
597 msg = self.data.jsonify(tm, self.plot_name, self.plot_type)
598 self.socket.send_string(msg)
588 599 socks = dict(self.poll.poll(5000))
589 600 if socks.get(self.socket) == zmq.POLLIN:
590 601 reply = self.socket.recv_string()
591 602 if reply == 'ok':
592 603 log.log("Response from server ok", self.name)
593 604 time.sleep(0.1)
594 605 continue
595 606 else:
596 607 log.warning(
597 608 "Malformed reply from server: {}".format(reply), self.name)
598 609 else:
599 610 log.warning(
600 611 "No response from server, retrying...", self.name)
601 self.sender_queue.put(msg)
612 self.sender_queue.put(self.data.tm)
602 613 self.socket.setsockopt(zmq.LINGER, 0)
603 614 self.socket.close()
604 615 self.poll.unregister(self.socket)
605 616 time.sleep(0.1)
606 617 self.socket = self.context.socket(zmq.REQ)
607 618 self.socket.connect(self.plot_server)
608 619 self.poll.register(self.socket, zmq.POLLIN)
609 620 break
610 621
611 622 def setup(self):
612 623 '''
613 624 This method should be implemented in the child class, the following
614 625 attributes should be set:
615 626
616 627 self.nrows: number of rows
617 628 self.ncols: number of cols
618 629 self.nplots: number of plots (channels or pairs)
619 630 self.ylabel: label for Y axes
620 631 self.titles: list of axes title
621 632
622 633 '''
623 634 raise NotImplementedError
624 635
625 636 def plot(self):
626 637 '''
627 638 Must be defined in the child class
628 639 '''
629 640 raise NotImplementedError
630 641
631 642 def run(self, dataOut, **kwargs):
632 643 '''
633 644 Main plotting routine
634 645 '''
635 646
636 647 if self.isConfig is False:
637 648 self.__setup(**kwargs)
638 649
639 650 t = getattr(dataOut, self.attr_time)
640 651
641 652 if dataOut.useLocalTime:
642 653 self.getDateTime = datetime.datetime.fromtimestamp
643 654 if not self.localtime:
644 655 t += time.timezone
645 656 else:
646 657 self.getDateTime = datetime.datetime.utcfromtimestamp
647 658 if self.localtime:
648 659 t -= time.timezone
649 660
650 if 'buffer' in self.plot_type:
651 if self.xmin is None:
652 self.tmin = t
661 if self.xmin is None:
662 self.tmin = t
663 if 'buffer' in self.plot_type:
653 664 self.xmin = self.getDateTime(t).hour
654 else:
655 self.tmin = (
656 self.getDateTime(t).replace(
657 hour=int(self.xmin),
658 minute=0,
659 second=0) - self.getDateTime(0)).total_seconds()
665 else:
666 self.tmin = (
667 self.getDateTime(t).replace(
668 hour=int(self.xmin),
669 minute=0,
670 second=0) - self.getDateTime(0)).total_seconds()
660 671
661 672 self.data.setup()
662 673 self.isConfig = True
663 674 if self.plot_server:
664 675 self.context = zmq.Context()
665 676 self.socket = self.context.socket(zmq.REQ)
666 677 self.socket.connect(self.plot_server)
667 678 self.poll = zmq.Poller()
668 679 self.poll.register(self.socket, zmq.POLLIN)
669 680
670 681 tm = getattr(dataOut, self.attr_time)
671 682
672 683 if not dataOut.useLocalTime and self.localtime:
673 684 tm -= time.timezone
674 685 if dataOut.useLocalTime and not self.localtime:
675 686 tm += time.timezone
676 687
677 if self.xaxis is 'time' and self.data and (tm - self.tmin) >= self.xrange*60*60:
688 if self.data and (tm - self.tmin) >= self.xrange*60*60:
678 689 self.save_counter = self.save_period
679 690 self.__plot()
680 self.xmin += self.xrange
681 if self.xmin >= 24:
682 self.xmin -= 24
691 if 'time' in self.xaxis:
692 self.xmin += self.xrange
693 if self.xmin >= 24:
694 self.xmin -= 24
683 695 self.tmin += self.xrange*60*60
684 696 self.data.setup()
685 697 self.clear_figures()
686 698
687 699 self.data.update(dataOut, tm)
688 700
689 701 if self.isPlotConfig is False:
690 702 self.__setup_plot()
691 703 self.isPlotConfig = True
692 704
693 705 if self.throttle == 0:
694 706 self.__plot()
695 707 else:
696 708 self.__throttle_plot(self.__plot)#, coerce=coerce)
697 709
698 710 def close(self):
699 711
700 712 if self.data and not self.data.flagNoData:
701 713 self.save_counter = self.save_period
702 714 self.__plot()
703 715 if self.data and not self.data.flagNoData and self.pause:
704 716 figpause(10)
705 717
@@ -1,651 +1,651
1 1 '''
2 2 Created on Jul 9, 2014
3 3 Modified on May 10, 2020
4 4
5 5 @author: Juan C. Espinoza
6 6 '''
7 7
8 8 import os
9 9 import datetime
10 10 import numpy
11 11
12 12 from schainpy.model.graphics.jroplot_base import Plot, plt
13 13
14 14
15 15 class SpectraPlot(Plot):
16 16 '''
17 17 Plot for Spectra data
18 18 '''
19 19
20 20 CODE = 'spc'
21 21 colormap = 'jet'
22 22 plot_name = 'Spectra'
23 23 plot_type = 'pcolor'
24 24
25 25 def setup(self):
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 = 3 * 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.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
36 36 self.ylabel = 'Range [km]'
37 37
38 38 def plot(self):
39 39 if self.xaxis == "frequency":
40 40 x = self.data.xrange[0]
41 41 self.xlabel = "Frequency (kHz)"
42 42 elif self.xaxis == "time":
43 43 x = self.data.xrange[1]
44 44 self.xlabel = "Time (ms)"
45 45 else:
46 46 x = self.data.xrange[2]
47 47 self.xlabel = "Velocity (m/s)"
48 48
49 49 if self.CODE == 'spc_moments':
50 50 x = self.data.xrange[2]
51 51 self.xlabel = "Velocity (m/s)"
52 52
53 53 self.titles = []
54 54
55 55 y = self.data.heights
56 56 self.y = y
57 57 z = self.data['spc']
58 58
59 59 for n, ax in enumerate(self.axes):
60 60 noise = self.data['noise'][n][-1]
61 61 if self.CODE == 'spc_moments':
62 62 mean = self.data['moments'][n, :, 1, :][-1]
63 63 if ax.firsttime:
64 64 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
65 65 self.xmin = self.xmin if self.xmin else -self.xmax
66 66 self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
67 67 self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
68 68 ax.plt = ax.pcolormesh(x, y, z[n].T,
69 69 vmin=self.zmin,
70 70 vmax=self.zmax,
71 71 cmap=plt.get_cmap(self.colormap)
72 72 )
73 73
74 74 if self.showprofile:
75 75 ax.plt_profile = self.pf_axes[n].plot(
76 76 self.data['rti'][n][-1], y)[0]
77 77 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
78 78 color="k", linestyle="dashed", lw=1)[0]
79 79 if self.CODE == 'spc_moments':
80 80 ax.plt_mean = ax.plot(mean, y, color='k')[0]
81 81 else:
82 82 ax.plt.set_array(z[n].T.ravel())
83 83 if self.showprofile:
84 84 ax.plt_profile.set_data(self.data['rti'][n][-1], y)
85 85 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
86 86 if self.CODE == 'spc_moments':
87 87 ax.plt_mean.set_data(mean, y)
88 88 self.titles.append('CH {}: {:3.2f}dB'.format(n, noise))
89 89
90 90
91 91 class CrossSpectraPlot(Plot):
92 92
93 93 CODE = 'cspc'
94 94 colormap = 'jet'
95 95 plot_name = 'CrossSpectra'
96 96 plot_type = 'pcolor'
97 97 zmin_coh = None
98 98 zmax_coh = None
99 99 zmin_phase = None
100 100 zmax_phase = None
101 101
102 102 def setup(self):
103 103
104 104 self.ncols = 4
105 105 self.nrows = len(self.data.pairs)
106 106 self.nplots = self.nrows * 4
107 107 self.width = 3.4 * self.ncols
108 108 self.height = 3 * self.nrows
109 109 self.ylabel = 'Range [km]'
110 110 self.showprofile = False
111 111 self.plots_adjust.update({'bottom': 0.08})
112 112
113 113 def plot(self):
114 114
115 115 if self.xaxis == "frequency":
116 116 x = self.data.xrange[0]
117 117 self.xlabel = "Frequency (kHz)"
118 118 elif self.xaxis == "time":
119 119 x = self.data.xrange[1]
120 120 self.xlabel = "Time (ms)"
121 121 else:
122 122 x = self.data.xrange[2]
123 123 self.xlabel = "Velocity (m/s)"
124 124
125 125 self.titles = []
126 126
127 127 y = self.data.heights
128 128 self.y = y
129 spc = self.data['spc']
130 cspc = self.data['cspc']
129 nspc = self.data['spc']
130 spc = self.data['cspc'][0]
131 cspc = self.data['cspc'][1]
131 132
132 133 for n in range(self.nrows):
133 noise = self.data['noise'][n][-1]
134 noise = self.data['noise'][:,-1]
134 135 pair = self.data.pairs[n]
135 136 ax = self.axes[4 * n]
136 spc0 = 10.*numpy.log10(spc[pair[0]]/self.data.factor)
137 137 if ax.firsttime:
138 138 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
139 139 self.xmin = self.xmin if self.xmin else -self.xmax
140 self.zmin = self.zmin if self.zmin else numpy.nanmin(spc)
141 self.zmax = self.zmax if self.zmax else numpy.nanmax(spc)
142 ax.plt = ax.pcolormesh(x , y , spc0.T,
140 self.zmin = self.zmin if self.zmin else numpy.nanmin(nspc)
141 self.zmax = self.zmax if self.zmax else numpy.nanmax(nspc)
142 ax.plt = ax.pcolormesh(x , y , nspc[pair[0]].T,
143 143 vmin=self.zmin,
144 144 vmax=self.zmax,
145 145 cmap=plt.get_cmap(self.colormap)
146 146 )
147 147 else:
148 ax.plt.set_array(spc0.T.ravel())
149 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise))
148 ax.plt.set_array(nspc[pair[0]].T.ravel())
149 self.titles.append('CH {}: {:3.2f}dB'.format(pair[0], noise[pair[0]]))
150 150
151 151 ax = self.axes[4 * n + 1]
152 spc1 = 10.*numpy.log10(spc[pair[1]]/self.data.factor)
153 152 if ax.firsttime:
154 ax.plt = ax.pcolormesh(x , y, spc1.T,
153 ax.plt = ax.pcolormesh(x , y, nspc[pair[1]].T,
155 154 vmin=self.zmin,
156 155 vmax=self.zmax,
157 156 cmap=plt.get_cmap(self.colormap)
158 157 )
159 158 else:
160 ax.plt.set_array(spc1.T.ravel())
161 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise))
159 ax.plt.set_array(nspc[pair[1]].T.ravel())
160 self.titles.append('CH {}: {:3.2f}dB'.format(pair[1], noise[pair[1]]))
162 161
163 162 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
164 163 coh = numpy.abs(out)
165 164 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
166 165
167 166 ax = self.axes[4 * n + 2]
168 167 if ax.firsttime:
169 168 ax.plt = ax.pcolormesh(x, y, coh.T,
170 169 vmin=0,
171 170 vmax=1,
172 171 cmap=plt.get_cmap(self.colormap_coh)
173 172 )
174 173 else:
175 174 ax.plt.set_array(coh.T.ravel())
176 175 self.titles.append(
177 176 'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
178 177
179 178 ax = self.axes[4 * n + 3]
180 179 if ax.firsttime:
181 180 ax.plt = ax.pcolormesh(x, y, phase.T,
182 181 vmin=-180,
183 182 vmax=180,
184 183 cmap=plt.get_cmap(self.colormap_phase)
185 184 )
186 185 else:
187 186 ax.plt.set_array(phase.T.ravel())
188 187 self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
189 188
190 189
191 190 class RTIPlot(Plot):
192 191 '''
193 192 Plot for RTI data
194 193 '''
195 194
196 195 CODE = 'rti'
197 196 colormap = 'jet'
198 197 plot_name = 'RTI'
199 198 plot_type = 'pcolorbuffer'
200 199
201 200 def setup(self):
202 201 self.xaxis = 'time'
203 202 self.ncols = 1
204 203 self.nrows = len(self.data.channels)
205 204 self.nplots = len(self.data.channels)
206 205 self.ylabel = 'Range [km]'
207 206 self.xlabel = 'Time'
208 207 self.cb_label = 'dB'
209 208 self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
210 209 self.titles = ['{} Channel {}'.format(
211 210 self.CODE.upper(), x) for x in range(self.nrows)]
212 211
213 212 def plot(self):
214 213 self.x = self.data.times
215 214 self.y = self.data.heights
216 215 self.z = self.data[self.CODE]
217 216 self.z = numpy.ma.masked_invalid(self.z)
218 217
219 218 if self.decimation is None:
220 219 x, y, z = self.fill_gaps(self.x, self.y, self.z)
221 220 else:
222 221 x, y, z = self.fill_gaps(*self.decimate())
223 222
224 223 for n, ax in enumerate(self.axes):
225 224 self.zmin = self.zmin if self.zmin else numpy.min(self.z)
226 225 self.zmax = self.zmax if self.zmax else numpy.max(self.z)
227 226 if ax.firsttime:
228 227 ax.plt = ax.pcolormesh(x, y, z[n].T,
229 228 vmin=self.zmin,
230 229 vmax=self.zmax,
231 230 cmap=plt.get_cmap(self.colormap)
232 231 )
233 232 if self.showprofile:
234 233 ax.plot_profile = self.pf_axes[n].plot(
235 234 self.data['rti'][n][-1], self.y)[0]
236 235 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(self.data['noise'][n][-1], len(self.y)), self.y,
237 236 color="k", linestyle="dashed", lw=1)[0]
238 237 else:
239 238 ax.collections.remove(ax.collections[0])
240 239 ax.plt = ax.pcolormesh(x, y, z[n].T,
241 240 vmin=self.zmin,
242 241 vmax=self.zmax,
243 242 cmap=plt.get_cmap(self.colormap)
244 243 )
245 244 if self.showprofile:
246 245 ax.plot_profile.set_data(self.data['rti'][n][-1], self.y)
247 246 ax.plot_noise.set_data(numpy.repeat(
248 247 self.data['noise'][n][-1], len(self.y)), self.y)
249 248
250 249
251 250 class CoherencePlot(RTIPlot):
252 251 '''
253 252 Plot for Coherence data
254 253 '''
255 254
256 255 CODE = 'coh'
257 256 plot_name = 'Coherence'
258 257
259 258 def setup(self):
260 259 self.xaxis = 'time'
261 260 self.ncols = 1
262 261 self.nrows = len(self.data.pairs)
263 262 self.nplots = len(self.data.pairs)
264 263 self.ylabel = 'Range [km]'
265 264 self.xlabel = 'Time'
266 265 self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
267 266 if self.CODE == 'coh':
268 267 self.cb_label = ''
269 268 self.titles = [
270 269 'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
271 270 else:
272 271 self.cb_label = 'Degrees'
273 272 self.titles = [
274 273 'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
275 274
276 275
277 276 class PhasePlot(CoherencePlot):
278 277 '''
279 278 Plot for Phase map data
280 279 '''
281 280
282 281 CODE = 'phase'
283 282 colormap = 'seismic'
284 283 plot_name = 'Phase'
285 284
286 285
287 286 class NoisePlot(Plot):
288 287 '''
289 288 Plot for noise
290 289 '''
291 290
292 291 CODE = 'noise'
293 292 plot_name = 'Noise'
294 293 plot_type = 'scatterbuffer'
295 294
296 295
297 296 def setup(self):
298 297 self.xaxis = 'time'
299 298 self.ncols = 1
300 299 self.nrows = 1
301 300 self.nplots = 1
302 301 self.ylabel = 'Intensity [dB]'
302 self.xlabel = 'Time'
303 303 self.titles = ['Noise']
304 304 self.colorbar = False
305 305
306 306 def plot(self):
307 307
308 308 x = self.data.times
309 309 xmin = self.data.min_time
310 310 xmax = xmin + self.xrange * 60 * 60
311 311 Y = self.data[self.CODE]
312 312
313 313 if self.axes[0].firsttime:
314 314 for ch in self.data.channels:
315 315 y = Y[ch]
316 316 self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
317 317 plt.legend()
318 318 else:
319 319 for ch in self.data.channels:
320 320 y = Y[ch]
321 321 self.axes[0].lines[ch].set_data(x, y)
322 322
323 323 self.ymin = numpy.nanmin(Y) - 5
324 324 self.ymax = numpy.nanmax(Y) + 5
325 325
326 326
327 327 class PowerProfilePlot(Plot):
328 328
329 329 CODE = 'spcprofile'
330 330 plot_name = 'Power Profile'
331 331 plot_type = 'scatter'
332 332 buffering = False
333 333
334 334 def setup(self):
335 335
336 336 self.ncols = 1
337 337 self.nrows = 1
338 338 self.nplots = 1
339 339 self.height = 4
340 340 self.width = 3
341 341 self.ylabel = 'Range [km]'
342 342 self.xlabel = 'Intensity [dB]'
343 343 self.titles = ['Power Profile']
344 344 self.colorbar = False
345 345
346 346 def plot(self):
347 347
348 348 y = self.data.heights
349 349 self.y = y
350 350
351 351 x = self.data['spcprofile']
352 352
353 353 if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
354 354 if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
355 355
356 356 if self.axes[0].firsttime:
357 357 for ch in self.data.channels:
358 358 self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
359 359 plt.legend()
360 360 else:
361 361 for ch in self.data.channels:
362 362 self.axes[0].lines[ch].set_data(x[ch], y)
363 363
364 364
365 365 class SpectraCutPlot(Plot):
366 366
367 367 CODE = 'spc_cut'
368 368 plot_name = 'Spectra Cut'
369 369 plot_type = 'scatter'
370 370 buffering = False
371 371
372 372 def setup(self):
373 373
374 374 self.nplots = len(self.data.channels)
375 375 self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
376 376 self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
377 377 self.width = 3.4 * self.ncols + 1.5
378 378 self.height = 3 * self.nrows
379 379 self.ylabel = 'Power [dB]'
380 380 self.colorbar = False
381 381 self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
382 382
383 383 def plot(self):
384 384 if self.xaxis == "frequency":
385 385 x = self.data.xrange[0][1:]
386 386 self.xlabel = "Frequency (kHz)"
387 387 elif self.xaxis == "time":
388 388 x = self.data.xrange[1]
389 389 self.xlabel = "Time (ms)"
390 390 else:
391 391 x = self.data.xrange[2]
392 392 self.xlabel = "Velocity (m/s)"
393 393
394 394 self.titles = []
395 395
396 396 y = self.data.heights
397 397 #self.y = y
398 398 z = self.data['spc_cut']
399 399
400 400 if self.height_index:
401 401 index = numpy.array(self.height_index)
402 402 else:
403 403 index = numpy.arange(0, len(y), int((len(y))/9))
404 404
405 405 for n, ax in enumerate(self.axes):
406 406 if ax.firsttime:
407 407 self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
408 408 self.xmin = self.xmin if self.xmin else -self.xmax
409 409 self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
410 410 self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
411 411 ax.plt = ax.plot(x, z[n, :, index].T)
412 412 labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
413 413 self.figures[0].legend(ax.plt, labels, loc='center right')
414 414 else:
415 415 for i, line in enumerate(ax.plt):
416 416 line.set_data(x, z[n, :, i])
417 417 self.titles.append('CH {}'.format(n))
418 418
419 419
420 420 class BeaconPhase(Plot):
421 421
422 422 __isConfig = None
423 423 __nsubplots = None
424 424
425 425 PREFIX = 'beacon_phase'
426 426
427 427 def __init__(self):
428 428 Plot.__init__(self)
429 429 self.timerange = 24*60*60
430 430 self.isConfig = False
431 431 self.__nsubplots = 1
432 432 self.counter_imagwr = 0
433 433 self.WIDTH = 800
434 434 self.HEIGHT = 400
435 435 self.WIDTHPROF = 120
436 436 self.HEIGHTPROF = 0
437 437 self.xdata = None
438 438 self.ydata = None
439 439
440 440 self.PLOT_CODE = BEACON_CODE
441 441
442 442 self.FTP_WEI = None
443 443 self.EXP_CODE = None
444 444 self.SUB_EXP_CODE = None
445 445 self.PLOT_POS = None
446 446
447 447 self.filename_phase = None
448 448
449 449 self.figfile = None
450 450
451 451 self.xmin = None
452 452 self.xmax = None
453 453
454 454 def getSubplots(self):
455 455
456 456 ncol = 1
457 457 nrow = 1
458 458
459 459 return nrow, ncol
460 460
461 461 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
462 462
463 463 self.__showprofile = showprofile
464 464 self.nplots = nplots
465 465
466 466 ncolspan = 7
467 467 colspan = 6
468 468 self.__nsubplots = 2
469 469
470 470 self.createFigure(id = id,
471 471 wintitle = wintitle,
472 472 widthplot = self.WIDTH+self.WIDTHPROF,
473 473 heightplot = self.HEIGHT+self.HEIGHTPROF,
474 474 show=show)
475 475
476 476 nrow, ncol = self.getSubplots()
477 477
478 478 self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
479 479
480 480 def save_phase(self, filename_phase):
481 481 f = open(filename_phase,'w+')
482 482 f.write('\n\n')
483 483 f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
484 484 f.write('DD MM YYYY HH MM SS pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
485 485 f.close()
486 486
487 487 def save_data(self, filename_phase, data, data_datetime):
488 488 f=open(filename_phase,'a')
489 489 timetuple_data = data_datetime.timetuple()
490 490 day = str(timetuple_data.tm_mday)
491 491 month = str(timetuple_data.tm_mon)
492 492 year = str(timetuple_data.tm_year)
493 493 hour = str(timetuple_data.tm_hour)
494 494 minute = str(timetuple_data.tm_min)
495 495 second = str(timetuple_data.tm_sec)
496 496 f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n')
497 497 f.close()
498 498
499 499 def plot(self):
500 500 log.warning('TODO: Not yet implemented...')
501 501
502 502 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
503 503 xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
504 504 timerange=None,
505 505 save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
506 506 server=None, folder=None, username=None, password=None,
507 507 ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
508 508
509 509 if dataOut.flagNoData:
510 510 return dataOut
511 511
512 512 if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
513 513 return
514 514
515 515 if pairsList == None:
516 516 pairsIndexList = dataOut.pairsIndexList[:10]
517 517 else:
518 518 pairsIndexList = []
519 519 for pair in pairsList:
520 520 if pair not in dataOut.pairsList:
521 521 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
522 522 pairsIndexList.append(dataOut.pairsList.index(pair))
523 523
524 524 if pairsIndexList == []:
525 525 return
526 526
527 527 # if len(pairsIndexList) > 4:
528 528 # pairsIndexList = pairsIndexList[0:4]
529 529
530 530 hmin_index = None
531 531 hmax_index = None
532 532
533 533 if hmin != None and hmax != None:
534 534 indexes = numpy.arange(dataOut.nHeights)
535 535 hmin_list = indexes[dataOut.heightList >= hmin]
536 536 hmax_list = indexes[dataOut.heightList <= hmax]
537 537
538 538 if hmin_list.any():
539 539 hmin_index = hmin_list[0]
540 540
541 541 if hmax_list.any():
542 542 hmax_index = hmax_list[-1]+1
543 543
544 544 x = dataOut.getTimeRange()
545 545 #y = dataOut.getHeiRange()
546 546
547 547
548 548 thisDatetime = dataOut.datatime
549 549
550 550 title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
551 551 xlabel = "Local Time"
552 552 ylabel = "Phase (degrees)"
553 553
554 554 update_figfile = False
555 555
556 556 nplots = len(pairsIndexList)
557 557 #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
558 558 phase_beacon = numpy.zeros(len(pairsIndexList))
559 559 for i in range(nplots):
560 560 pair = dataOut.pairsList[pairsIndexList[i]]
561 561 ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
562 562 powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
563 563 powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
564 564 avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
565 565 phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
566 566
567 567 if dataOut.beacon_heiIndexList:
568 568 phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
569 569 else:
570 570 phase_beacon[i] = numpy.average(phase)
571 571
572 572 if not self.isConfig:
573 573
574 574 nplots = len(pairsIndexList)
575 575
576 576 self.setup(id=id,
577 577 nplots=nplots,
578 578 wintitle=wintitle,
579 579 showprofile=showprofile,
580 580 show=show)
581 581
582 582 if timerange != None:
583 583 self.timerange = timerange
584 584
585 585 self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
586 586
587 587 if ymin == None: ymin = 0
588 588 if ymax == None: ymax = 360
589 589
590 590 self.FTP_WEI = ftp_wei
591 591 self.EXP_CODE = exp_code
592 592 self.SUB_EXP_CODE = sub_exp_code
593 593 self.PLOT_POS = plot_pos
594 594
595 595 self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
596 596 self.isConfig = True
597 597 self.figfile = figfile
598 598 self.xdata = numpy.array([])
599 599 self.ydata = numpy.array([])
600 600
601 601 update_figfile = True
602 602
603 603 #open file beacon phase
604 604 path = '%s%03d' %(self.PREFIX, self.id)
605 605 beacon_file = os.path.join(path,'%s.txt'%self.name)
606 606 self.filename_phase = os.path.join(figpath,beacon_file)
607 607 #self.save_phase(self.filename_phase)
608 608
609 609
610 610 #store data beacon phase
611 611 #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
612 612
613 613 self.setWinTitle(title)
614 614
615 615
616 616 title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
617 617
618 618 legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
619 619
620 620 axes = self.axesList[0]
621 621
622 622 self.xdata = numpy.hstack((self.xdata, x[0:1]))
623 623
624 624 if len(self.ydata)==0:
625 625 self.ydata = phase_beacon.reshape(-1,1)
626 626 else:
627 627 self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
628 628
629 629
630 630 axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
631 631 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
632 632 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
633 633 XAxisAsTime=True, grid='both'
634 634 )
635 635
636 636 self.draw()
637 637
638 638 if dataOut.ltctime >= self.xmax:
639 639 self.counter_imagwr = wr_period
640 640 self.isConfig = False
641 641 update_figfile = True
642 642
643 643 self.save(figpath=figpath,
644 644 figfile=figfile,
645 645 save=save,
646 646 ftp=ftp,
647 647 wr_period=wr_period,
648 648 thisDatetime=thisDatetime,
649 649 update_figfile=update_figfile)
650 650
651 651 return dataOut No newline at end of file
@@ -1,23 +1,24
1 1 '''
2 2
3 3 $Author: murco $
4 4 $Id: JRODataIO.py 169 2012-11-19 21:57:03Z murco $
5 5 '''
6 6
7 7 from .jroIO_voltage import *
8 8 from .jroIO_spectra import *
9 9 from .jroIO_heispectra import *
10 10 from .jroIO_usrp import *
11 11 from .jroIO_digitalRF import *
12 12 from .jroIO_kamisr import *
13 13 from .jroIO_param import *
14 14 from .jroIO_hf import *
15 15
16 16 from .jroIO_madrigal import *
17 17
18 18 from .bltrIO_param import *
19 19 from .bltrIO_spectra import *
20 20 from .jroIO_mira35c import *
21 21 from .julIO_param import *
22 22
23 from .pxIO_param import * No newline at end of file
23 from .pxIO_param import *
24 from .jroIO_simulator import * No newline at end of file
@@ -1,1580 +1,1580
1 1 """
2 2 Created on Jul 2, 2014
3 3
4 4 @author: roj-idl71
5 5 """
6 6 import os
7 7 import sys
8 8 import glob
9 9 import time
10 10 import numpy
11 11 import fnmatch
12 12 import inspect
13 13 import time
14 14 import datetime
15 15 import zmq
16 16
17 from schainpy.model.proc.jroproc_base import Operation
17 from schainpy.model.proc.jroproc_base import Operation, MPDecorator
18 18 from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader, SystemHeader, RadarControllerHeader, ProcessingHeader
19 19 from schainpy.model.data.jroheaderIO import get_dtype_index, get_numpy_dtype, get_procflag_dtype, get_dtype_width
20 20 from schainpy.utils import log
21 21 import schainpy.admin
22 22
23 23 LOCALTIME = True
24 24 DT_DIRECTIVES = {
25 25 '%Y': 4,
26 26 '%y': 2,
27 27 '%m': 2,
28 28 '%d': 2,
29 29 '%j': 3,
30 30 '%H': 2,
31 31 '%M': 2,
32 32 '%S': 2,
33 33 '%f': 6
34 34 }
35 35
36 36
37 37 def isNumber(cad):
38 38 """
39 39 Chequea si el conjunto de caracteres que componen un string puede ser convertidos a un numero.
40 40
41 41 Excepciones:
42 42 Si un determinado string no puede ser convertido a numero
43 43 Input:
44 44 str, string al cual se le analiza para determinar si convertible a un numero o no
45 45
46 46 Return:
47 47 True : si el string es uno numerico
48 48 False : no es un string numerico
49 49 """
50 50 try:
51 51 float(cad)
52 52 return True
53 53 except:
54 54 return False
55 55
56 56
57 57 def isFileInEpoch(filename, startUTSeconds, endUTSeconds):
58 58 """
59 59 Esta funcion determina si un archivo de datos se encuentra o no dentro del rango de fecha especificado.
60 60
61 61 Inputs:
62 62 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
63 63
64 64 startUTSeconds : fecha inicial del rango seleccionado. La fecha esta dada en
65 65 segundos contados desde 01/01/1970.
66 66 endUTSeconds : fecha final del rango seleccionado. La fecha esta dada en
67 67 segundos contados desde 01/01/1970.
68 68
69 69 Return:
70 70 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
71 71 fecha especificado, de lo contrario retorna False.
72 72
73 73 Excepciones:
74 74 Si el archivo no existe o no puede ser abierto
75 75 Si la cabecera no puede ser leida.
76 76
77 77 """
78 78 basicHeaderObj = BasicHeader(LOCALTIME)
79 79
80 80 try:
81 81 fp = open(filename, 'rb')
82 82 except IOError:
83 83 print("The file %s can't be opened" % (filename))
84 84 return 0
85 85
86 86 sts = basicHeaderObj.read(fp)
87 87 fp.close()
88 88
89 89 if not(sts):
90 90 print("Skipping the file %s because it has not a valid header" % (filename))
91 91 return 0
92 92
93 93 if not ((startUTSeconds <= basicHeaderObj.utc) and (endUTSeconds > basicHeaderObj.utc)):
94 94 return 0
95 95
96 96 return 1
97 97
98 98
99 99 def isTimeInRange(thisTime, startTime, endTime):
100 100 if endTime >= startTime:
101 101 if (thisTime < startTime) or (thisTime > endTime):
102 102 return 0
103 103 return 1
104 104 else:
105 105 if (thisTime < startTime) and (thisTime > endTime):
106 106 return 0
107 107 return 1
108 108
109 109
110 110 def isFileInTimeRange(filename, startDate, endDate, startTime, endTime):
111 111 """
112 112 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
113 113
114 114 Inputs:
115 115 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
116 116
117 117 startDate : fecha inicial del rango seleccionado en formato datetime.date
118 118
119 119 endDate : fecha final del rango seleccionado en formato datetime.date
120 120
121 121 startTime : tiempo inicial del rango seleccionado en formato datetime.time
122 122
123 123 endTime : tiempo final del rango seleccionado en formato datetime.time
124 124
125 125 Return:
126 126 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
127 127 fecha especificado, de lo contrario retorna False.
128 128
129 129 Excepciones:
130 130 Si el archivo no existe o no puede ser abierto
131 131 Si la cabecera no puede ser leida.
132 132
133 133 """
134 134
135 135 try:
136 136 fp = open(filename, 'rb')
137 137 except IOError:
138 138 print("The file %s can't be opened" % (filename))
139 139 return None
140 140
141 141 firstBasicHeaderObj = BasicHeader(LOCALTIME)
142 142 systemHeaderObj = SystemHeader()
143 143 radarControllerHeaderObj = RadarControllerHeader()
144 144 processingHeaderObj = ProcessingHeader()
145 145
146 146 lastBasicHeaderObj = BasicHeader(LOCALTIME)
147 147
148 148 sts = firstBasicHeaderObj.read(fp)
149 149
150 150 if not(sts):
151 151 print("[Reading] Skipping the file %s because it has not a valid header" % (filename))
152 152 return None
153 153
154 154 if not systemHeaderObj.read(fp):
155 155 return None
156 156
157 157 if not radarControllerHeaderObj.read(fp):
158 158 return None
159 159
160 160 if not processingHeaderObj.read(fp):
161 161 return None
162 162
163 163 filesize = os.path.getsize(filename)
164 164
165 165 offset = processingHeaderObj.blockSize + 24 # header size
166 166
167 167 if filesize <= offset:
168 168 print("[Reading] %s: This file has not enough data" % filename)
169 169 return None
170 170
171 171 fp.seek(-offset, 2)
172 172
173 173 sts = lastBasicHeaderObj.read(fp)
174 174
175 175 fp.close()
176 176
177 177 thisDatetime = lastBasicHeaderObj.datatime
178 178 thisTime_last_block = thisDatetime.time()
179 179
180 180 thisDatetime = firstBasicHeaderObj.datatime
181 181 thisDate = thisDatetime.date()
182 182 thisTime_first_block = thisDatetime.time()
183 183
184 184 # General case
185 185 # o>>>>>>>>>>>>>><<<<<<<<<<<<<<o
186 186 #-----------o----------------------------o-----------
187 187 # startTime endTime
188 188
189 189 if endTime >= startTime:
190 190 if (thisTime_last_block < startTime) or (thisTime_first_block > endTime):
191 191 return None
192 192
193 193 return thisDatetime
194 194
195 195 # If endTime < startTime then endTime belongs to the next day
196 196
197 197 #<<<<<<<<<<<o o>>>>>>>>>>>
198 198 #-----------o----------------------------o-----------
199 199 # endTime startTime
200 200
201 201 if (thisDate == startDate) and (thisTime_last_block < startTime):
202 202 return None
203 203
204 204 if (thisDate == endDate) and (thisTime_first_block > endTime):
205 205 return None
206 206
207 207 if (thisTime_last_block < startTime) and (thisTime_first_block > endTime):
208 208 return None
209 209
210 210 return thisDatetime
211 211
212 212
213 213 def isFolderInDateRange(folder, startDate=None, endDate=None):
214 214 """
215 215 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
216 216
217 217 Inputs:
218 218 folder : nombre completo del directorio.
219 219 Su formato deberia ser "/path_root/?YYYYDDD"
220 220
221 221 siendo:
222 222 YYYY : Anio (ejemplo 2015)
223 223 DDD : Dia del anio (ejemplo 305)
224 224
225 225 startDate : fecha inicial del rango seleccionado en formato datetime.date
226 226
227 227 endDate : fecha final del rango seleccionado en formato datetime.date
228 228
229 229 Return:
230 230 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
231 231 fecha especificado, de lo contrario retorna False.
232 232 Excepciones:
233 233 Si el directorio no tiene el formato adecuado
234 234 """
235 235
236 236 basename = os.path.basename(folder)
237 237
238 238 if not isRadarFolder(basename):
239 239 print("The folder %s has not the rigth format" % folder)
240 240 return 0
241 241
242 242 if startDate and endDate:
243 243 thisDate = getDateFromRadarFolder(basename)
244 244
245 245 if thisDate < startDate:
246 246 return 0
247 247
248 248 if thisDate > endDate:
249 249 return 0
250 250
251 251 return 1
252 252
253 253
254 254 def isFileInDateRange(filename, startDate=None, endDate=None):
255 255 """
256 256 Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado.
257 257
258 258 Inputs:
259 259 filename : nombre completo del archivo de datos en formato Jicamarca (.r)
260 260
261 261 Su formato deberia ser "?YYYYDDDsss"
262 262
263 263 siendo:
264 264 YYYY : Anio (ejemplo 2015)
265 265 DDD : Dia del anio (ejemplo 305)
266 266 sss : set
267 267
268 268 startDate : fecha inicial del rango seleccionado en formato datetime.date
269 269
270 270 endDate : fecha final del rango seleccionado en formato datetime.date
271 271
272 272 Return:
273 273 Boolean : Retorna True si el archivo de datos contiene datos en el rango de
274 274 fecha especificado, de lo contrario retorna False.
275 275 Excepciones:
276 276 Si el archivo no tiene el formato adecuado
277 277 """
278 278
279 279 basename = os.path.basename(filename)
280 280
281 281 if not isRadarFile(basename):
282 282 print("The filename %s has not the rigth format" % filename)
283 283 return 0
284 284
285 285 if startDate and endDate:
286 286 thisDate = getDateFromRadarFile(basename)
287 287
288 288 if thisDate < startDate:
289 289 return 0
290 290
291 291 if thisDate > endDate:
292 292 return 0
293 293
294 294 return 1
295 295
296 296
297 297 def getFileFromSet(path, ext, set):
298 298 validFilelist = []
299 299 fileList = os.listdir(path)
300 300
301 301 # 0 1234 567 89A BCDE
302 302 # H YYYY DDD SSS .ext
303 303
304 304 for thisFile in fileList:
305 305 try:
306 306 year = int(thisFile[1:5])
307 307 doy = int(thisFile[5:8])
308 308 except:
309 309 continue
310 310
311 311 if (os.path.splitext(thisFile)[-1].lower() != ext.lower()):
312 312 continue
313 313
314 314 validFilelist.append(thisFile)
315 315
316 316 myfile = fnmatch.filter(
317 317 validFilelist, '*%4.4d%3.3d%3.3d*' % (year, doy, set))
318 318
319 319 if len(myfile) != 0:
320 320 return myfile[0]
321 321 else:
322 322 filename = '*%4.4d%3.3d%3.3d%s' % (year, doy, set, ext.lower())
323 323 print('the filename %s does not exist' % filename)
324 324 print('...going to the last file: ')
325 325
326 326 if validFilelist:
327 327 validFilelist = sorted(validFilelist, key=str.lower)
328 328 return validFilelist[-1]
329 329
330 330 return None
331 331
332 332
333 333 def getlastFileFromPath(path, ext):
334 334 """
335 335 Depura el fileList dejando solo los que cumplan el formato de "PYYYYDDDSSS.ext"
336 336 al final de la depuracion devuelve el ultimo file de la lista que quedo.
337 337
338 338 Input:
339 339 fileList : lista conteniendo todos los files (sin path) que componen una determinada carpeta
340 340 ext : extension de los files contenidos en una carpeta
341 341
342 342 Return:
343 343 El ultimo file de una determinada carpeta, no se considera el path.
344 344 """
345 345 validFilelist = []
346 346 fileList = os.listdir(path)
347 347
348 348 # 0 1234 567 89A BCDE
349 349 # H YYYY DDD SSS .ext
350 350
351 351 for thisFile in fileList:
352 352
353 353 year = thisFile[1:5]
354 354 if not isNumber(year):
355 355 continue
356 356
357 357 doy = thisFile[5:8]
358 358 if not isNumber(doy):
359 359 continue
360 360
361 361 year = int(year)
362 362 doy = int(doy)
363 363
364 364 if (os.path.splitext(thisFile)[-1].lower() != ext.lower()):
365 365 continue
366 366
367 367 validFilelist.append(thisFile)
368 368
369 369 if validFilelist:
370 370 validFilelist = sorted(validFilelist, key=str.lower)
371 371 return validFilelist[-1]
372 372
373 373 return None
374 374
375 375
376 376 def isRadarFolder(folder):
377 377 try:
378 378 year = int(folder[1:5])
379 379 doy = int(folder[5:8])
380 380 except:
381 381 return 0
382 382
383 383 return 1
384 384
385 385
386 386 def isRadarFile(file):
387 387 try:
388 388 year = int(file[1:5])
389 389 doy = int(file[5:8])
390 390 set = int(file[8:11])
391 391 except:
392 392 return 0
393 393
394 394 return 1
395 395
396 396
397 397 def getDateFromRadarFile(file):
398 398 try:
399 399 year = int(file[1:5])
400 400 doy = int(file[5:8])
401 401 set = int(file[8:11])
402 402 except:
403 403 return None
404 404
405 405 thisDate = datetime.date(year, 1, 1) + datetime.timedelta(doy - 1)
406 406 return thisDate
407 407
408 408
409 409 def getDateFromRadarFolder(folder):
410 410 try:
411 411 year = int(folder[1:5])
412 412 doy = int(folder[5:8])
413 413 except:
414 414 return None
415 415
416 416 thisDate = datetime.date(year, 1, 1) + datetime.timedelta(doy - 1)
417 417 return thisDate
418 418
419 419 def parse_format(s, fmt):
420 420
421 421 for i in range(fmt.count('%')):
422 422 x = fmt.index('%')
423 423 d = DT_DIRECTIVES[fmt[x:x+2]]
424 424 fmt = fmt.replace(fmt[x:x+2], s[x:x+d])
425 425 return fmt
426 426
427 427 class Reader(object):
428 428
429 429 c = 3E8
430 430 isConfig = False
431 431 dtype = None
432 432 pathList = []
433 433 filenameList = []
434 434 datetimeList = []
435 435 filename = None
436 436 ext = None
437 437 flagIsNewFile = 1
438 438 flagDiscontinuousBlock = 0
439 439 flagIsNewBlock = 0
440 440 flagNoMoreFiles = 0
441 441 fp = None
442 442 firstHeaderSize = 0
443 443 basicHeaderSize = 24
444 444 versionFile = 1103
445 445 fileSize = None
446 446 fileSizeByHeader = None
447 447 fileIndex = -1
448 448 profileIndex = None
449 449 blockIndex = 0
450 450 nTotalBlocks = 0
451 451 maxTimeStep = 30
452 452 lastUTTime = None
453 453 datablock = None
454 454 dataOut = None
455 455 getByBlock = False
456 456 path = None
457 457 startDate = None
458 458 endDate = None
459 459 startTime = datetime.time(0, 0, 0)
460 460 endTime = datetime.time(23, 59, 59)
461 461 set = None
462 462 expLabel = ""
463 463 online = False
464 464 delay = 60
465 465 nTries = 3 # quantity tries
466 466 nFiles = 3 # number of files for searching
467 467 walk = True
468 468 getblock = False
469 469 nTxs = 1
470 470 realtime = False
471 471 blocksize = 0
472 472 blocktime = None
473 473 warnings = True
474 474 verbose = True
475 475 server = None
476 476 format = None
477 477 oneDDict = None
478 478 twoDDict = None
479 479 independentParam = None
480 480 filefmt = None
481 481 folderfmt = None
482 482 open_file = open
483 483 open_mode = 'rb'
484 484
485 485 def run(self):
486 486
487 487 raise NotImplementedError
488 488
489 489 def getAllowedArgs(self):
490 490 if hasattr(self, '__attrs__'):
491 491 return self.__attrs__
492 492 else:
493 493 return inspect.getargspec(self.run).args
494 494
495 495 def set_kwargs(self, **kwargs):
496 496
497 497 for key, value in kwargs.items():
498 498 setattr(self, key, value)
499 499
500 500 def find_folders(self, path, startDate, endDate, folderfmt, last=False):
501 501
502 502 folders = [x for f in path.split(',')
503 503 for x in os.listdir(f) if os.path.isdir(os.path.join(f, x))]
504 504 folders.sort()
505 505
506 506 if last:
507 507 folders = [folders[-1]]
508 508
509 509 for folder in folders:
510 510 try:
511 511 dt = datetime.datetime.strptime(parse_format(folder, folderfmt), folderfmt).date()
512 512 if dt >= startDate and dt <= endDate:
513 513 yield os.path.join(path, folder)
514 514 else:
515 515 log.log('Skiping folder {}'.format(folder), self.name)
516 516 except Exception as e:
517 517 log.log('Skiping folder {}'.format(folder), self.name)
518 518 continue
519 519 return
520 520
521 521 def find_files(self, folders, ext, filefmt, startDate=None, endDate=None,
522 522 expLabel='', last=False):
523 523
524 524 for path in folders:
525 525 files = glob.glob1(path, '*{}'.format(ext))
526 526 files.sort()
527 527 if last:
528 528 if files:
529 529 fo = files[-1]
530 530 try:
531 531 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
532 532 yield os.path.join(path, expLabel, fo)
533 533 except Exception as e:
534 534 pass
535 535 return
536 536 else:
537 537 return
538 538
539 539 for fo in files:
540 540 try:
541 541 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
542 542 if dt >= startDate and dt <= endDate:
543 543 yield os.path.join(path, expLabel, fo)
544 544 else:
545 545 log.log('Skiping file {}'.format(fo), self.name)
546 546 except Exception as e:
547 547 log.log('Skiping file {}'.format(fo), self.name)
548 548 continue
549 549
550 550 def searchFilesOffLine(self, path, startDate, endDate,
551 551 expLabel, ext, walk,
552 552 filefmt, folderfmt):
553 553 """Search files in offline mode for the given arguments
554 554
555 555 Return:
556 556 Generator of files
557 557 """
558 558
559 559 if walk:
560 560 folders = self.find_folders(
561 561 path, startDate, endDate, folderfmt)
562 562 else:
563 563 folders = path.split(',')
564 564
565 565 return self.find_files(
566 566 folders, ext, filefmt, startDate, endDate, expLabel)
567 567
568 568 def searchFilesOnLine(self, path, startDate, endDate,
569 569 expLabel, ext, walk,
570 570 filefmt, folderfmt):
571 571 """Search for the last file of the last folder
572 572
573 573 Arguments:
574 574 path : carpeta donde estan contenidos los files que contiene data
575 575 expLabel : Nombre del subexperimento (subfolder)
576 576 ext : extension de los files
577 577 walk : Si es habilitado no realiza busquedas dentro de los ubdirectorios (doypath)
578 578
579 579 Return:
580 580 generator with the full path of last filename
581 581 """
582 582
583 583 if walk:
584 584 folders = self.find_folders(
585 585 path, startDate, endDate, folderfmt, last=True)
586 586 else:
587 587 folders = path.split(',')
588 588
589 589 return self.find_files(
590 590 folders, ext, filefmt, startDate, endDate, expLabel, last=True)
591 591
592 592 def setNextFile(self):
593 593 """Set the next file to be readed open it and parse de file header"""
594 594
595 595 while True:
596 596 if self.fp != None:
597 597 self.fp.close()
598 598
599 599 if self.online:
600 600 newFile = self.setNextFileOnline()
601 601 else:
602 602 newFile = self.setNextFileOffline()
603 603
604 604 if not(newFile):
605 605 if self.online:
606 606 raise schainpy.admin.SchainError('Time to wait for new files reach')
607 607 else:
608 608 if self.fileIndex == -1:
609 609 raise schainpy.admin.SchainWarning('No files found in the given path')
610 610 else:
611 611 raise schainpy.admin.SchainWarning('No more files to read')
612 612
613 613 if self.verifyFile(self.filename):
614 614 break
615 615
616 616 log.log('Opening file: %s' % self.filename, self.name)
617 617
618 618 self.readFirstHeader()
619 619 self.nReadBlocks = 0
620 620
621 621 def setNextFileOnline(self):
622 622 """Check for the next file to be readed in online mode.
623 623
624 624 Set:
625 625 self.filename
626 626 self.fp
627 627 self.filesize
628 628
629 629 Return:
630 630 boolean
631 631
632 632 """
633 633 nextFile = True
634 634 nextDay = False
635 635
636 636 for nFiles in range(self.nFiles+1):
637 637 for nTries in range(self.nTries):
638 638 fullfilename, filename = self.checkForRealPath(nextFile, nextDay)
639 639 if fullfilename is not None:
640 640 break
641 641 log.warning(
642 642 "Waiting %0.2f sec for the next file: \"%s\" , try %02d ..." % (self.delay, filename, nTries + 1),
643 643 self.name)
644 644 time.sleep(self.delay)
645 645 nextFile = False
646 646 continue
647 647
648 648 if fullfilename is not None:
649 649 break
650 650
651 651 self.nTries = 1
652 652 nextFile = True
653 653
654 654 if nFiles == (self.nFiles - 1):
655 655 log.log('Trying with next day...', self.name)
656 656 nextDay = True
657 657 self.nTries = 3
658 658
659 659 if fullfilename:
660 660 self.fileSize = os.path.getsize(fullfilename)
661 661 self.filename = fullfilename
662 662 self.flagIsNewFile = 1
663 663 if self.fp != None:
664 664 self.fp.close()
665 665 self.fp = self.open_file(fullfilename, self.open_mode)
666 666 self.flagNoMoreFiles = 0
667 667 self.fileIndex += 1
668 668 return 1
669 669 else:
670 670 return 0
671 671
672 672 def setNextFileOffline(self):
673 673 """Open the next file to be readed in offline mode"""
674 674
675 675 try:
676 676 filename = next(self.filenameList)
677 677 self.fileIndex +=1
678 678 except StopIteration:
679 679 self.flagNoMoreFiles = 1
680 680 return 0
681 681
682 682 self.filename = filename
683 683 self.fileSize = os.path.getsize(filename)
684 684 self.fp = self.open_file(filename, self.open_mode)
685 685 self.flagIsNewFile = 1
686 686
687 687 return 1
688 688
689 689 @staticmethod
690 690 def isDateTimeInRange(dt, startDate, endDate, startTime, endTime):
691 691 """Check if the given datetime is in range"""
692 692
693 693 if startDate <= dt.date() <= endDate:
694 694 if startTime <= dt.time() <= endTime:
695 695 return True
696 696 return False
697 697
698 698 def verifyFile(self, filename):
699 699 """Check for a valid file
700 700
701 701 Arguments:
702 702 filename -- full path filename
703 703
704 704 Return:
705 705 boolean
706 706 """
707 707
708 708 return True
709 709
710 710 def checkForRealPath(self, nextFile, nextDay):
711 711 """Check if the next file to be readed exists"""
712 712
713 713 raise NotImplementedError
714 714
715 715 def readFirstHeader(self):
716 716 """Parse the file header"""
717 717
718 718 pass
719 719
720 720 class JRODataReader(Reader):
721 721
722 722 utc = 0
723 723 nReadBlocks = 0
724 724 foldercounter = 0
725 725 firstHeaderSize = 0
726 726 basicHeaderSize = 24
727 727 __isFirstTimeOnline = 1
728 728 filefmt = "*%Y%j***"
729 729 folderfmt = "*%Y%j"
730 730 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'online', 'delay', 'walk']
731 731
732 732 def getDtypeWidth(self):
733 733
734 734 dtype_index = get_dtype_index(self.dtype)
735 735 dtype_width = get_dtype_width(dtype_index)
736 736
737 737 return dtype_width
738 738
739 739 def checkForRealPath(self, nextFile, nextDay):
740 740 """Check if the next file to be readed exists.
741 741
742 742 Example :
743 743 nombre correcto del file es .../.../D2009307/P2009307367.ext
744 744
745 745 Entonces la funcion prueba con las siguientes combinaciones
746 746 .../.../y2009307367.ext
747 747 .../.../Y2009307367.ext
748 748 .../.../x2009307/y2009307367.ext
749 749 .../.../x2009307/Y2009307367.ext
750 750 .../.../X2009307/y2009307367.ext
751 751 .../.../X2009307/Y2009307367.ext
752 752 siendo para este caso, la ultima combinacion de letras, identica al file buscado
753 753
754 754 Return:
755 755 str -- fullpath of the file
756 756 """
757 757
758 758
759 759 if nextFile:
760 760 self.set += 1
761 761 if nextDay:
762 762 self.set = 0
763 763 self.doy += 1
764 764 foldercounter = 0
765 765 prefixDirList = [None, 'd', 'D']
766 766 if self.ext.lower() == ".r": # voltage
767 767 prefixFileList = ['d', 'D']
768 768 elif self.ext.lower() == ".pdata": # spectra
769 769 prefixFileList = ['p', 'P']
770 770
771 771 # barrido por las combinaciones posibles
772 772 for prefixDir in prefixDirList:
773 773 thispath = self.path
774 774 if prefixDir != None:
775 775 # formo el nombre del directorio xYYYYDDD (x=d o x=D)
776 776 if foldercounter == 0:
777 777 thispath = os.path.join(self.path, "%s%04d%03d" %
778 778 (prefixDir, self.year, self.doy))
779 779 else:
780 780 thispath = os.path.join(self.path, "%s%04d%03d_%02d" % (
781 781 prefixDir, self.year, self.doy, foldercounter))
782 782 for prefixFile in prefixFileList: # barrido por las dos combinaciones posibles de "D"
783 783 # formo el nombre del file xYYYYDDDSSS.ext
784 784 filename = "%s%04d%03d%03d%s" % (prefixFile, self.year, self.doy, self.set, self.ext)
785 785 fullfilename = os.path.join(
786 786 thispath, filename)
787 787
788 788 if os.path.exists(fullfilename):
789 789 return fullfilename, filename
790 790
791 791 return None, filename
792 792
793 793 def __waitNewBlock(self):
794 794 """
795 795 Return 1 si se encontro un nuevo bloque de datos, 0 de otra forma.
796 796
797 797 Si el modo de lectura es OffLine siempre retorn 0
798 798 """
799 799 if not self.online:
800 800 return 0
801 801
802 802 if (self.nReadBlocks >= self.processingHeaderObj.dataBlocksPerFile):
803 803 return 0
804 804
805 805 currentPointer = self.fp.tell()
806 806
807 807 neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
808 808
809 809 for nTries in range(self.nTries):
810 810
811 811 self.fp.close()
812 812 self.fp = open(self.filename, 'rb')
813 813 self.fp.seek(currentPointer)
814 814
815 815 self.fileSize = os.path.getsize(self.filename)
816 816 currentSize = self.fileSize - currentPointer
817 817
818 818 if (currentSize >= neededSize):
819 819 self.basicHeaderObj.read(self.fp)
820 820 return 1
821 821
822 822 if self.fileSize == self.fileSizeByHeader:
823 823 # self.flagEoF = True
824 824 return 0
825 825
826 826 print("[Reading] Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1))
827 827 time.sleep(self.delay)
828 828
829 829 return 0
830 830
831 831 def waitDataBlock(self, pointer_location, blocksize=None):
832 832
833 833 currentPointer = pointer_location
834 834 if blocksize is None:
835 835 neededSize = self.processingHeaderObj.blockSize # + self.basicHeaderSize
836 836 else:
837 837 neededSize = blocksize
838 838
839 839 for nTries in range(self.nTries):
840 840 self.fp.close()
841 841 self.fp = open(self.filename, 'rb')
842 842 self.fp.seek(currentPointer)
843 843
844 844 self.fileSize = os.path.getsize(self.filename)
845 845 currentSize = self.fileSize - currentPointer
846 846
847 847 if (currentSize >= neededSize):
848 848 return 1
849 849
850 850 log.warning(
851 851 "Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1),
852 852 self.name
853 853 )
854 854 time.sleep(self.delay)
855 855
856 856 return 0
857 857
858 858 def __setNewBlock(self):
859 859
860 860 if self.fp == None:
861 861 return 0
862 862
863 863 if self.flagIsNewFile:
864 864 self.lastUTTime = self.basicHeaderObj.utc
865 865 return 1
866 866
867 867 if self.realtime:
868 868 self.flagDiscontinuousBlock = 1
869 869 if not(self.setNextFile()):
870 870 return 0
871 871 else:
872 872 return 1
873 873
874 874 currentSize = self.fileSize - self.fp.tell()
875 875 neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
876 876
877 877 if (currentSize >= neededSize):
878 878 self.basicHeaderObj.read(self.fp)
879 879 self.lastUTTime = self.basicHeaderObj.utc
880 880 return 1
881 881
882 882 if self.__waitNewBlock():
883 883 self.lastUTTime = self.basicHeaderObj.utc
884 884 return 1
885 885
886 886 if not(self.setNextFile()):
887 887 return 0
888 888
889 889 deltaTime = self.basicHeaderObj.utc - self.lastUTTime
890 890 self.lastUTTime = self.basicHeaderObj.utc
891 891
892 892 self.flagDiscontinuousBlock = 0
893 893
894 894 if deltaTime > self.maxTimeStep:
895 895 self.flagDiscontinuousBlock = 1
896 896
897 897 return 1
898 898
899 899 def readNextBlock(self):
900 900
901 901 while True:
902 902 self.__setNewBlock()
903 903
904 904 if not(self.readBlock()):
905 905 return 0
906 906
907 907 self.getBasicHeader()
908 908
909 909 if not self.isDateTimeInRange(self.dataOut.datatime, self.startDate, self.endDate, self.startTime, self.endTime):
910 910 print("[Reading] Block No. %d/%d -> %s [Skipping]" % (self.nReadBlocks,
911 911 self.processingHeaderObj.dataBlocksPerFile,
912 912 self.dataOut.datatime.ctime()))
913 913 continue
914 914
915 915 break
916 916
917 917 if self.verbose:
918 918 print("[Reading] Block No. %d/%d -> %s" % (self.nReadBlocks,
919 919 self.processingHeaderObj.dataBlocksPerFile,
920 920 self.dataOut.datatime.ctime()))
921 921 return 1
922 922
923 923 def readFirstHeader(self):
924 924
925 925 self.basicHeaderObj.read(self.fp)
926 926 self.systemHeaderObj.read(self.fp)
927 927 self.radarControllerHeaderObj.read(self.fp)
928 928 self.processingHeaderObj.read(self.fp)
929 929 self.firstHeaderSize = self.basicHeaderObj.size
930 930
931 931 datatype = int(numpy.log2((self.processingHeaderObj.processFlags &
932 932 PROCFLAG.DATATYPE_MASK)) - numpy.log2(PROCFLAG.DATATYPE_CHAR))
933 933 if datatype == 0:
934 934 datatype_str = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
935 935 elif datatype == 1:
936 936 datatype_str = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
937 937 elif datatype == 2:
938 938 datatype_str = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
939 939 elif datatype == 3:
940 940 datatype_str = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
941 941 elif datatype == 4:
942 942 datatype_str = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
943 943 elif datatype == 5:
944 944 datatype_str = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
945 945 else:
946 946 raise ValueError('Data type was not defined')
947 947
948 948 self.dtype = datatype_str
949 949 #self.ippSeconds = 2 * 1000 * self.radarControllerHeaderObj.ipp / self.c
950 950 self.fileSizeByHeader = self.processingHeaderObj.dataBlocksPerFile * self.processingHeaderObj.blockSize + \
951 951 self.firstHeaderSize + self.basicHeaderSize * \
952 952 (self.processingHeaderObj.dataBlocksPerFile - 1)
953 953 # self.dataOut.channelList = numpy.arange(self.systemHeaderObj.numChannels)
954 954 # self.dataOut.channelIndexList = numpy.arange(self.systemHeaderObj.numChannels)
955 955 self.getBlockDimension()
956 956
957 957 def verifyFile(self, filename, msgFlag=True):
958 958
959 959 msg = None
960 960
961 961 try:
962 962 fp = open(filename, 'rb')
963 963 except IOError:
964 964
965 965 if msgFlag:
966 966 print("[Reading] File %s can't be opened" % (filename))
967 967
968 968 return False
969 969
970 970 if self.waitDataBlock(0):
971 971 basicHeaderObj = BasicHeader(LOCALTIME)
972 972 systemHeaderObj = SystemHeader()
973 973 radarControllerHeaderObj = RadarControllerHeader()
974 974 processingHeaderObj = ProcessingHeader()
975 975
976 976 if not(basicHeaderObj.read(fp)):
977 977 fp.close()
978 978 return False
979 979
980 980 if not(systemHeaderObj.read(fp)):
981 981 fp.close()
982 982 return False
983 983
984 984 if not(radarControllerHeaderObj.read(fp)):
985 985 fp.close()
986 986 return False
987 987
988 988 if not(processingHeaderObj.read(fp)):
989 989 fp.close()
990 990 return False
991 991
992 992 if not self.online:
993 993 dt1 = basicHeaderObj.datatime
994 994 fp.seek(self.fileSize-processingHeaderObj.blockSize-24)
995 995 if not(basicHeaderObj.read(fp)):
996 996 fp.close()
997 997 return False
998 998 dt2 = basicHeaderObj.datatime
999 999 if not self.isDateTimeInRange(dt1, self.startDate, self.endDate, self.startTime, self.endTime) and not \
1000 1000 self.isDateTimeInRange(dt2, self.startDate, self.endDate, self.startTime, self.endTime):
1001 1001 return False
1002 1002
1003 1003 fp.close()
1004 1004
1005 1005 return True
1006 1006
1007 1007 def findDatafiles(self, path, startDate=None, endDate=None, expLabel='', ext='.r', walk=True, include_path=False):
1008 1008
1009 1009 path_empty = True
1010 1010
1011 1011 dateList = []
1012 1012 pathList = []
1013 1013
1014 1014 multi_path = path.split(',')
1015 1015
1016 1016 if not walk:
1017 1017
1018 1018 for single_path in multi_path:
1019 1019
1020 1020 if not os.path.isdir(single_path):
1021 1021 continue
1022 1022
1023 1023 fileList = glob.glob1(single_path, "*" + ext)
1024 1024
1025 1025 if not fileList:
1026 1026 continue
1027 1027
1028 1028 path_empty = False
1029 1029
1030 1030 fileList.sort()
1031 1031
1032 1032 for thisFile in fileList:
1033 1033
1034 1034 if not os.path.isfile(os.path.join(single_path, thisFile)):
1035 1035 continue
1036 1036
1037 1037 if not isRadarFile(thisFile):
1038 1038 continue
1039 1039
1040 1040 if not isFileInDateRange(thisFile, startDate, endDate):
1041 1041 continue
1042 1042
1043 1043 thisDate = getDateFromRadarFile(thisFile)
1044 1044
1045 1045 if thisDate in dateList or single_path in pathList:
1046 1046 continue
1047 1047
1048 1048 dateList.append(thisDate)
1049 1049 pathList.append(single_path)
1050 1050
1051 1051 else:
1052 1052 for single_path in multi_path:
1053 1053
1054 1054 if not os.path.isdir(single_path):
1055 1055 continue
1056 1056
1057 1057 dirList = []
1058 1058
1059 1059 for thisPath in os.listdir(single_path):
1060 1060
1061 1061 if not os.path.isdir(os.path.join(single_path, thisPath)):
1062 1062 continue
1063 1063
1064 1064 if not isRadarFolder(thisPath):
1065 1065 continue
1066 1066
1067 1067 if not isFolderInDateRange(thisPath, startDate, endDate):
1068 1068 continue
1069 1069
1070 1070 dirList.append(thisPath)
1071 1071
1072 1072 if not dirList:
1073 1073 continue
1074 1074
1075 1075 dirList.sort()
1076 1076
1077 1077 for thisDir in dirList:
1078 1078
1079 1079 datapath = os.path.join(single_path, thisDir, expLabel)
1080 1080 fileList = glob.glob1(datapath, "*" + ext)
1081 1081
1082 1082 if not fileList:
1083 1083 continue
1084 1084
1085 1085 path_empty = False
1086 1086
1087 1087 thisDate = getDateFromRadarFolder(thisDir)
1088 1088
1089 1089 pathList.append(datapath)
1090 1090 dateList.append(thisDate)
1091 1091
1092 1092 dateList.sort()
1093 1093
1094 1094 if walk:
1095 1095 pattern_path = os.path.join(multi_path[0], "[dYYYYDDD]", expLabel)
1096 1096 else:
1097 1097 pattern_path = multi_path[0]
1098 1098
1099 1099 if path_empty:
1100 1100 raise schainpy.admin.SchainError("[Reading] No *%s files in %s for %s to %s" % (ext, pattern_path, startDate, endDate))
1101 1101 else:
1102 1102 if not dateList:
1103 1103 raise schainpy.admin.SchainError("[Reading] Date range selected invalid [%s - %s]: No *%s files in %s)" % (startDate, endDate, ext, path))
1104 1104
1105 1105 if include_path:
1106 1106 return dateList, pathList
1107 1107
1108 1108 return dateList
1109 1109
1110 1110 def setup(self, **kwargs):
1111 1111
1112 1112 self.set_kwargs(**kwargs)
1113 1113 if not self.ext.startswith('.'):
1114 1114 self.ext = '.{}'.format(self.ext)
1115 1115
1116 1116 if self.server is not None:
1117 1117 if 'tcp://' in self.server:
1118 1118 address = server
1119 1119 else:
1120 1120 address = 'ipc:///tmp/%s' % self.server
1121 1121 self.server = address
1122 1122 self.context = zmq.Context()
1123 1123 self.receiver = self.context.socket(zmq.PULL)
1124 1124 self.receiver.connect(self.server)
1125 1125 time.sleep(0.5)
1126 1126 print('[Starting] ReceiverData from {}'.format(self.server))
1127 1127 else:
1128 1128 self.server = None
1129 1129 if self.path == None:
1130 1130 raise ValueError("[Reading] The path is not valid")
1131 1131
1132 1132 if self.online:
1133 1133 log.log("[Reading] Searching files in online mode...", self.name)
1134 1134
1135 1135 for nTries in range(self.nTries):
1136 1136 fullpath = self.searchFilesOnLine(self.path, self.startDate,
1137 1137 self.endDate, self.expLabel, self.ext, self.walk,
1138 1138 self.filefmt, self.folderfmt)
1139 1139
1140 1140 try:
1141 1141 fullpath = next(fullpath)
1142 1142 except:
1143 1143 fullpath = None
1144 1144
1145 1145 if fullpath:
1146 1146 break
1147 1147
1148 1148 log.warning(
1149 1149 'Waiting {} sec for a valid file in {}: try {} ...'.format(
1150 1150 self.delay, self.path, nTries + 1),
1151 1151 self.name)
1152 1152 time.sleep(self.delay)
1153 1153
1154 1154 if not(fullpath):
1155 1155 raise schainpy.admin.SchainError(
1156 1156 'There isn\'t any valid file in {}'.format(self.path))
1157 1157
1158 1158 pathname, filename = os.path.split(fullpath)
1159 1159 self.year = int(filename[1:5])
1160 1160 self.doy = int(filename[5:8])
1161 1161 self.set = int(filename[8:11]) - 1
1162 1162 else:
1163 1163 log.log("Searching files in {}".format(self.path), self.name)
1164 1164 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
1165 1165 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
1166 1166
1167 1167 self.setNextFile()
1168 1168
1169 1169 return
1170 1170
1171 1171 def getBasicHeader(self):
1172 1172
1173 1173 self.dataOut.utctime = self.basicHeaderObj.utc + self.basicHeaderObj.miliSecond / \
1174 1174 1000. + self.profileIndex * self.radarControllerHeaderObj.ippSeconds
1175 1175
1176 1176 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
1177 1177
1178 1178 self.dataOut.timeZone = self.basicHeaderObj.timeZone
1179 1179
1180 1180 self.dataOut.dstFlag = self.basicHeaderObj.dstFlag
1181 1181
1182 1182 self.dataOut.errorCount = self.basicHeaderObj.errorCount
1183 1183
1184 1184 self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
1185 1185
1186 1186 self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
1187 1187
1188 1188 # self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock*self.nTxs
1189 1189
1190 1190 def getFirstHeader(self):
1191 1191
1192 1192 raise NotImplementedError
1193 1193
1194 1194 def getData(self):
1195 1195
1196 1196 raise NotImplementedError
1197 1197
1198 1198 def hasNotDataInBuffer(self):
1199 1199
1200 1200 raise NotImplementedError
1201 1201
1202 1202 def readBlock(self):
1203 1203
1204 1204 raise NotImplementedError
1205 1205
1206 1206 def isEndProcess(self):
1207 1207
1208 1208 return self.flagNoMoreFiles
1209 1209
1210 1210 def printReadBlocks(self):
1211 1211
1212 1212 print("[Reading] Number of read blocks per file %04d" % self.nReadBlocks)
1213 1213
1214 1214 def printTotalBlocks(self):
1215 1215
1216 1216 print("[Reading] Number of read blocks %04d" % self.nTotalBlocks)
1217 1217
1218 1218 def run(self, **kwargs):
1219 1219 """
1220 1220
1221 1221 Arguments:
1222 1222 path :
1223 1223 startDate :
1224 1224 endDate :
1225 1225 startTime :
1226 1226 endTime :
1227 1227 set :
1228 1228 expLabel :
1229 1229 ext :
1230 1230 online :
1231 1231 delay :
1232 1232 walk :
1233 1233 getblock :
1234 1234 nTxs :
1235 1235 realtime :
1236 1236 blocksize :
1237 1237 blocktime :
1238 1238 skip :
1239 1239 cursor :
1240 1240 warnings :
1241 1241 server :
1242 1242 verbose :
1243 1243 format :
1244 1244 oneDDict :
1245 1245 twoDDict :
1246 1246 independentParam :
1247 1247 """
1248 1248
1249 1249 if not(self.isConfig):
1250 1250 self.setup(**kwargs)
1251 1251 self.isConfig = True
1252 1252 if self.server is None:
1253 1253 self.getData()
1254 1254 else:
1255 1255 self.getFromServer()
1256 1256
1257 1257
1258 1258 class JRODataWriter(Reader):
1259 1259
1260 1260 """
1261 1261 Esta clase permite escribir datos a archivos procesados (.r o ,pdata). La escritura
1262 1262 de los datos siempre se realiza por bloques.
1263 1263 """
1264 1264
1265 1265 setFile = None
1266 1266 profilesPerBlock = None
1267 1267 blocksPerFile = None
1268 1268 nWriteBlocks = 0
1269 1269 fileDate = None
1270 1270
1271 1271 def __init__(self, dataOut=None):
1272 1272 raise NotImplementedError
1273 1273
1274 1274 def hasAllDataInBuffer(self):
1275 1275 raise NotImplementedError
1276 1276
1277 1277 def setBlockDimension(self):
1278 1278 raise NotImplementedError
1279 1279
1280 1280 def writeBlock(self):
1281 1281 raise NotImplementedError
1282 1282
1283 1283 def putData(self):
1284 1284 raise NotImplementedError
1285 1285
1286 1286 def getDtypeWidth(self):
1287 1287
1288 1288 dtype_index = get_dtype_index(self.dtype)
1289 1289 dtype_width = get_dtype_width(dtype_index)
1290 1290
1291 1291 return dtype_width
1292 1292
1293 1293 def getProcessFlags(self):
1294 1294
1295 1295 processFlags = 0
1296 1296
1297 1297 dtype_index = get_dtype_index(self.dtype)
1298 1298 procflag_dtype = get_procflag_dtype(dtype_index)
1299 1299
1300 1300 processFlags += procflag_dtype
1301 1301
1302 1302 if self.dataOut.flagDecodeData:
1303 1303 processFlags += PROCFLAG.DECODE_DATA
1304 1304
1305 1305 if self.dataOut.flagDeflipData:
1306 1306 processFlags += PROCFLAG.DEFLIP_DATA
1307 1307
1308 1308 if self.dataOut.code is not None:
1309 1309 processFlags += PROCFLAG.DEFINE_PROCESS_CODE
1310 1310
1311 1311 if self.dataOut.nCohInt > 1:
1312 1312 processFlags += PROCFLAG.COHERENT_INTEGRATION
1313 1313
1314 1314 if self.dataOut.type == "Spectra":
1315 1315 if self.dataOut.nIncohInt > 1:
1316 1316 processFlags += PROCFLAG.INCOHERENT_INTEGRATION
1317 1317
1318 1318 if self.dataOut.data_dc is not None:
1319 1319 processFlags += PROCFLAG.SAVE_CHANNELS_DC
1320 1320
1321 1321 if self.dataOut.flagShiftFFT:
1322 1322 processFlags += PROCFLAG.SHIFT_FFT_DATA
1323 1323
1324 1324 return processFlags
1325 1325
1326 1326 def setBasicHeader(self):
1327 1327
1328 1328 self.basicHeaderObj.size = self.basicHeaderSize # bytes
1329 1329 self.basicHeaderObj.version = self.versionFile
1330 1330 self.basicHeaderObj.dataBlock = self.nTotalBlocks
1331 1331 utc = numpy.floor(self.dataOut.utctime)
1332 1332 milisecond = (self.dataOut.utctime - utc) * 1000.0
1333 1333 self.basicHeaderObj.utc = utc
1334 1334 self.basicHeaderObj.miliSecond = milisecond
1335 1335 self.basicHeaderObj.timeZone = self.dataOut.timeZone
1336 1336 self.basicHeaderObj.dstFlag = self.dataOut.dstFlag
1337 1337 self.basicHeaderObj.errorCount = self.dataOut.errorCount
1338 1338
1339 1339 def setFirstHeader(self):
1340 1340 """
1341 1341 Obtiene una copia del First Header
1342 1342
1343 1343 Affected:
1344 1344
1345 1345 self.basicHeaderObj
1346 1346 self.systemHeaderObj
1347 1347 self.radarControllerHeaderObj
1348 1348 self.processingHeaderObj self.
1349 1349
1350 1350 Return:
1351 1351 None
1352 1352 """
1353 1353
1354 1354 raise NotImplementedError
1355 1355
1356 1356 def __writeFirstHeader(self):
1357 1357 """
1358 1358 Escribe el primer header del file es decir el Basic header y el Long header (SystemHeader, RadarControllerHeader, ProcessingHeader)
1359 1359
1360 1360 Affected:
1361 1361 __dataType
1362 1362
1363 1363 Return:
1364 1364 None
1365 1365 """
1366 1366
1367 1367 # CALCULAR PARAMETROS
1368 1368
1369 1369 sizeLongHeader = self.systemHeaderObj.size + \
1370 1370 self.radarControllerHeaderObj.size + self.processingHeaderObj.size
1371 1371 self.basicHeaderObj.size = self.basicHeaderSize + sizeLongHeader
1372 1372
1373 1373 self.basicHeaderObj.write(self.fp)
1374 1374 self.systemHeaderObj.write(self.fp)
1375 1375 self.radarControllerHeaderObj.write(self.fp)
1376 1376 self.processingHeaderObj.write(self.fp)
1377 1377
1378 1378 def __setNewBlock(self):
1379 1379 """
1380 1380 Si es un nuevo file escribe el First Header caso contrario escribe solo el Basic Header
1381 1381
1382 1382 Return:
1383 1383 0 : si no pudo escribir nada
1384 1384 1 : Si escribio el Basic el First Header
1385 1385 """
1386 1386 if self.fp == None:
1387 1387 self.setNextFile()
1388 1388
1389 1389 if self.flagIsNewFile:
1390 1390 return 1
1391 1391
1392 1392 if self.blockIndex < self.processingHeaderObj.dataBlocksPerFile:
1393 1393 self.basicHeaderObj.write(self.fp)
1394 1394 return 1
1395 1395
1396 1396 if not(self.setNextFile()):
1397 1397 return 0
1398 1398
1399 1399 return 1
1400 1400
1401 1401 def writeNextBlock(self):
1402 1402 """
1403 1403 Selecciona el bloque siguiente de datos y los escribe en un file
1404 1404
1405 1405 Return:
1406 1406 0 : Si no hizo pudo escribir el bloque de datos
1407 1407 1 : Si no pudo escribir el bloque de datos
1408 1408 """
1409 1409 if not(self.__setNewBlock()):
1410 1410 return 0
1411 1411
1412 1412 self.writeBlock()
1413 1413
1414 1414 print("[Writing] Block No. %d/%d" % (self.blockIndex,
1415 1415 self.processingHeaderObj.dataBlocksPerFile))
1416 1416
1417 1417 return 1
1418 1418
1419 1419 def setNextFile(self):
1420 1420 """Determina el siguiente file que sera escrito
1421 1421
1422 1422 Affected:
1423 1423 self.filename
1424 1424 self.subfolder
1425 1425 self.fp
1426 1426 self.setFile
1427 1427 self.flagIsNewFile
1428 1428
1429 1429 Return:
1430 1430 0 : Si el archivo no puede ser escrito
1431 1431 1 : Si el archivo esta listo para ser escrito
1432 1432 """
1433 1433 ext = self.ext
1434 1434 path = self.path
1435 1435
1436 1436 if self.fp != None:
1437 1437 self.fp.close()
1438 1438
1439 1439 if not os.path.exists(path):
1440 1440 os.mkdir(path)
1441 1441
1442 1442 timeTuple = time.localtime(self.dataOut.utctime)
1443 1443 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year, timeTuple.tm_yday)
1444 1444
1445 1445 fullpath = os.path.join(path, subfolder)
1446 1446 setFile = self.setFile
1447 1447
1448 1448 if not(os.path.exists(fullpath)):
1449 1449 os.mkdir(fullpath)
1450 1450 setFile = -1 # inicializo mi contador de seteo
1451 1451 else:
1452 1452 filesList = os.listdir(fullpath)
1453 1453 if len(filesList) > 0:
1454 1454 filesList = sorted(filesList, key=str.lower)
1455 1455 filen = filesList[-1]
1456 1456 # el filename debera tener el siguiente formato
1457 1457 # 0 1234 567 89A BCDE (hex)
1458 1458 # x YYYY DDD SSS .ext
1459 1459 if isNumber(filen[8:11]):
1460 1460 # inicializo mi contador de seteo al seteo del ultimo file
1461 1461 setFile = int(filen[8:11])
1462 1462 else:
1463 1463 setFile = -1
1464 1464 else:
1465 1465 setFile = -1 # inicializo mi contador de seteo
1466 1466
1467 1467 setFile += 1
1468 1468
1469 1469 # If this is a new day it resets some values
1470 1470 if self.dataOut.datatime.date() > self.fileDate:
1471 1471 setFile = 0
1472 1472 self.nTotalBlocks = 0
1473 1473
1474 1474 filen = '{}{:04d}{:03d}{:03d}{}'.format(
1475 1475 self.optchar, timeTuple.tm_year, timeTuple.tm_yday, setFile, ext)
1476 1476
1477 1477 filename = os.path.join(path, subfolder, filen)
1478 1478
1479 1479 fp = open(filename, 'wb')
1480 1480
1481 1481 self.blockIndex = 0
1482 1482 self.filename = filename
1483 1483 self.subfolder = subfolder
1484 1484 self.fp = fp
1485 1485 self.setFile = setFile
1486 1486 self.flagIsNewFile = 1
1487 1487 self.fileDate = self.dataOut.datatime.date()
1488 1488 self.setFirstHeader()
1489 1489
1490 1490 print('[Writing] Opening file: %s' % self.filename)
1491 1491
1492 1492 self.__writeFirstHeader()
1493 1493
1494 1494 return 1
1495 1495
1496 1496 def setup(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=None, ext=None, datatype=4):
1497 1497 """
1498 1498 Setea el tipo de formato en la cual sera guardada la data y escribe el First Header
1499 1499
1500 1500 Inputs:
1501 1501 path : directory where data will be saved
1502 1502 profilesPerBlock : number of profiles per block
1503 1503 set : initial file set
1504 1504 datatype : An integer number that defines data type:
1505 1505 0 : int8 (1 byte)
1506 1506 1 : int16 (2 bytes)
1507 1507 2 : int32 (4 bytes)
1508 1508 3 : int64 (8 bytes)
1509 1509 4 : float32 (4 bytes)
1510 1510 5 : double64 (8 bytes)
1511 1511
1512 1512 Return:
1513 1513 0 : Si no realizo un buen seteo
1514 1514 1 : Si realizo un buen seteo
1515 1515 """
1516 1516
1517 1517 if ext == None:
1518 1518 ext = self.ext
1519 1519
1520 1520 self.ext = ext.lower()
1521 1521
1522 1522 self.path = path
1523 1523
1524 1524 if set is None:
1525 1525 self.setFile = -1
1526 1526 else:
1527 1527 self.setFile = set - 1
1528 1528
1529 1529 self.blocksPerFile = blocksPerFile
1530 1530 self.profilesPerBlock = profilesPerBlock
1531 1531 self.dataOut = dataOut
1532 1532 self.fileDate = self.dataOut.datatime.date()
1533 1533 self.dtype = self.dataOut.dtype
1534 1534
1535 1535 if datatype is not None:
1536 1536 self.dtype = get_numpy_dtype(datatype)
1537 1537
1538 1538 if not(self.setNextFile()):
1539 1539 print("[Writing] There isn't a next file")
1540 1540 return 0
1541 1541
1542 1542 self.setBlockDimension()
1543 1543
1544 1544 return 1
1545 1545
1546 1546 def run(self, dataOut, path, blocksPerFile=100, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs):
1547 1547
1548 1548 if not(self.isConfig):
1549 1549
1550 1550 self.setup(dataOut, path, blocksPerFile, profilesPerBlock=profilesPerBlock,
1551 1551 set=set, ext=ext, datatype=datatype, **kwargs)
1552 1552 self.isConfig = True
1553 1553
1554 1554 self.dataOut = dataOut
1555 1555 self.putData()
1556 1556 return self.dataOut
1557 1557
1558 @MPDecorator
1558 1559 class printInfo(Operation):
1559 1560
1560 1561 def __init__(self):
1561 1562
1562 1563 Operation.__init__(self)
1563 1564 self.__printInfo = True
1564 1565
1565 1566 def run(self, dataOut, headers = ['systemHeaderObj', 'radarControllerHeaderObj', 'processingHeaderObj']):
1566 1567 if self.__printInfo == False:
1567 return dataOut
1568 return
1568 1569
1569 1570 for header in headers:
1570 1571 if hasattr(dataOut, header):
1571 1572 obj = getattr(dataOut, header)
1572 1573 if hasattr(obj, 'printInfo'):
1573 1574 obj.printInfo()
1574 1575 else:
1575 1576 print(obj)
1576 1577 else:
1577 1578 log.warning('Header {} Not found in object'.format(header))
1578 1579
1579 1580 self.__printInfo = False
1580 return dataOut No newline at end of file
@@ -1,1557 +1,1574
1 1 import sys
2 2 import numpy,math
3 3 from scipy import interpolate
4 4 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
5 5 from schainpy.model.data.jrodata import Voltage
6 6 from schainpy.utils import log
7 7 from time import time
8 8
9 9
10 10
11 11 class VoltageProc(ProcessingUnit):
12 12
13 13 def __init__(self):
14 14
15 15 ProcessingUnit.__init__(self)
16 16
17 17 self.dataOut = Voltage()
18 18 self.flip = 1
19 19 self.setupReq = False
20 20
21 21 def run(self):
22 22
23 23 if self.dataIn.type == 'AMISR':
24 24 self.__updateObjFromAmisrInput()
25 25
26 26 if self.dataIn.type == 'Voltage':
27 27 self.dataOut.copy(self.dataIn)
28 28
29 29 def __updateObjFromAmisrInput(self):
30 30
31 31 self.dataOut.timeZone = self.dataIn.timeZone
32 32 self.dataOut.dstFlag = self.dataIn.dstFlag
33 33 self.dataOut.errorCount = self.dataIn.errorCount
34 34 self.dataOut.useLocalTime = self.dataIn.useLocalTime
35 35
36 36 self.dataOut.flagNoData = self.dataIn.flagNoData
37 37 self.dataOut.data = self.dataIn.data
38 38 self.dataOut.utctime = self.dataIn.utctime
39 39 self.dataOut.channelList = self.dataIn.channelList
40 40 #self.dataOut.timeInterval = self.dataIn.timeInterval
41 41 self.dataOut.heightList = self.dataIn.heightList
42 42 self.dataOut.nProfiles = self.dataIn.nProfiles
43 43
44 44 self.dataOut.nCohInt = self.dataIn.nCohInt
45 45 self.dataOut.ippSeconds = self.dataIn.ippSeconds
46 46 self.dataOut.frequency = self.dataIn.frequency
47 47
48 48 self.dataOut.azimuth = self.dataIn.azimuth
49 49 self.dataOut.zenith = self.dataIn.zenith
50 50
51 51 self.dataOut.beam.codeList = self.dataIn.beam.codeList
52 52 self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
53 53 self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
54 54
55 55
56 56 class selectChannels(Operation):
57 57
58 58 def run(self, dataOut, channelList):
59 59
60 60 channelIndexList = []
61 61 self.dataOut = dataOut
62 62 for channel in channelList:
63 63 if channel not in self.dataOut.channelList:
64 64 raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList)))
65 65
66 66 index = self.dataOut.channelList.index(channel)
67 67 channelIndexList.append(index)
68 68 self.selectChannelsByIndex(channelIndexList)
69 69 return self.dataOut
70 70
71 71 def selectChannelsByIndex(self, channelIndexList):
72 72 """
73 73 Selecciona un bloque de datos en base a canales segun el channelIndexList
74 74
75 75 Input:
76 76 channelIndexList : lista sencilla de canales a seleccionar por ej. [2,3,7]
77 77
78 78 Affected:
79 79 self.dataOut.data
80 80 self.dataOut.channelIndexList
81 81 self.dataOut.nChannels
82 82 self.dataOut.m_ProcessingHeader.totalSpectra
83 83 self.dataOut.systemHeaderObj.numChannels
84 84 self.dataOut.m_ProcessingHeader.blockSize
85 85
86 86 Return:
87 87 None
88 88 """
89 89
90 90 for channelIndex in channelIndexList:
91 91 if channelIndex not in self.dataOut.channelIndexList:
92 92 raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
93 93
94 94 if self.dataOut.type == 'Voltage':
95 95 if self.dataOut.flagDataAsBlock:
96 96 """
97 97 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
98 98 """
99 99 data = self.dataOut.data[channelIndexList,:,:]
100 100 else:
101 101 data = self.dataOut.data[channelIndexList,:]
102 102
103 103 self.dataOut.data = data
104 104 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
105 105 self.dataOut.channelList = range(len(channelIndexList))
106 106
107 107 elif self.dataOut.type == 'Spectra':
108 108 data_spc = self.dataOut.data_spc[channelIndexList, :]
109 109 data_dc = self.dataOut.data_dc[channelIndexList, :]
110 110
111 111 self.dataOut.data_spc = data_spc
112 112 self.dataOut.data_dc = data_dc
113 113
114 114 # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
115 115 self.dataOut.channelList = range(len(channelIndexList))
116 116 self.__selectPairsByChannel(channelIndexList)
117 117
118 118 return 1
119 119
120 120 def __selectPairsByChannel(self, channelList=None):
121 121
122 122 if channelList == None:
123 123 return
124 124
125 125 pairsIndexListSelected = []
126 126 for pairIndex in self.dataOut.pairsIndexList:
127 127 # First pair
128 128 if self.dataOut.pairsList[pairIndex][0] not in channelList:
129 129 continue
130 130 # Second pair
131 131 if self.dataOut.pairsList[pairIndex][1] not in channelList:
132 132 continue
133 133
134 134 pairsIndexListSelected.append(pairIndex)
135 135
136 136 if not pairsIndexListSelected:
137 137 self.dataOut.data_cspc = None
138 138 self.dataOut.pairsList = []
139 139 return
140 140
141 141 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
142 142 self.dataOut.pairsList = [self.dataOut.pairsList[i]
143 143 for i in pairsIndexListSelected]
144 144
145 145 return
146 146
147 147 class selectHeights(Operation):
148 148
149 149 def run(self, dataOut, minHei=None, maxHei=None):
150 150 """
151 151 Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
152 152 minHei <= height <= maxHei
153 153
154 154 Input:
155 155 minHei : valor minimo de altura a considerar
156 156 maxHei : valor maximo de altura a considerar
157 157
158 158 Affected:
159 159 Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
160 160
161 161 Return:
162 162 1 si el metodo se ejecuto con exito caso contrario devuelve 0
163 163 """
164 164
165 165 self.dataOut = dataOut
166 166
167 167 if minHei == None:
168 168 minHei = self.dataOut.heightList[0]
169 169
170 170 if maxHei == None:
171 171 maxHei = self.dataOut.heightList[-1]
172 172
173 173 if (minHei < self.dataOut.heightList[0]):
174 174 minHei = self.dataOut.heightList[0]
175 175
176 176 if (maxHei > self.dataOut.heightList[-1]):
177 177 maxHei = self.dataOut.heightList[-1]
178 178
179 179 minIndex = 0
180 180 maxIndex = 0
181 181 heights = self.dataOut.heightList
182 182
183 183 inda = numpy.where(heights >= minHei)
184 184 indb = numpy.where(heights <= maxHei)
185 185
186 186 try:
187 187 minIndex = inda[0][0]
188 188 except:
189 189 minIndex = 0
190 190
191 191 try:
192 192 maxIndex = indb[0][-1]
193 193 except:
194 194 maxIndex = len(heights)
195 195
196 196 self.selectHeightsByIndex(minIndex, maxIndex)
197 197
198 198 return self.dataOut
199 199
200 200 def selectHeightsByIndex(self, minIndex, maxIndex):
201 201 """
202 202 Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
203 203 minIndex <= index <= maxIndex
204 204
205 205 Input:
206 206 minIndex : valor de indice minimo de altura a considerar
207 207 maxIndex : valor de indice maximo de altura a considerar
208 208
209 209 Affected:
210 210 self.dataOut.data
211 211 self.dataOut.heightList
212 212
213 213 Return:
214 214 1 si el metodo se ejecuto con exito caso contrario devuelve 0
215 215 """
216 216
217 217 if self.dataOut.type == 'Voltage':
218 218 if (minIndex < 0) or (minIndex > maxIndex):
219 219 raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
220 220
221 221 if (maxIndex >= self.dataOut.nHeights):
222 222 maxIndex = self.dataOut.nHeights
223 223
224 224 #voltage
225 225 if self.dataOut.flagDataAsBlock:
226 226 """
227 227 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
228 228 """
229 229 data = self.dataOut.data[:,:, minIndex:maxIndex]
230 230 else:
231 231 data = self.dataOut.data[:, minIndex:maxIndex]
232 232
233 233 # firstHeight = self.dataOut.heightList[minIndex]
234 234
235 235 self.dataOut.data = data
236 236 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
237 237
238 238 if self.dataOut.nHeights <= 1:
239 239 raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
240 240 elif self.dataOut.type == 'Spectra':
241 241 if (minIndex < 0) or (minIndex > maxIndex):
242 242 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
243 243 minIndex, maxIndex))
244 244
245 245 if (maxIndex >= self.dataOut.nHeights):
246 246 maxIndex = self.dataOut.nHeights - 1
247 247
248 248 # Spectra
249 249 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
250 250
251 251 data_cspc = None
252 252 if self.dataOut.data_cspc is not None:
253 253 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
254 254
255 255 data_dc = None
256 256 if self.dataOut.data_dc is not None:
257 257 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
258 258
259 259 self.dataOut.data_spc = data_spc
260 260 self.dataOut.data_cspc = data_cspc
261 261 self.dataOut.data_dc = data_dc
262 262
263 263 self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
264 264
265 265 return 1
266 266
267 267
268 268 class filterByHeights(Operation):
269 269
270 270 def run(self, dataOut, window):
271 271
272 272 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
273 273
274 274 if window == None:
275 275 window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
276 276
277 277 newdelta = deltaHeight * window
278 278 r = dataOut.nHeights % window
279 279 newheights = (dataOut.nHeights-r)/window
280 280
281 281 if newheights <= 1:
282 282 raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
283 283
284 284 if dataOut.flagDataAsBlock:
285 285 """
286 286 Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
287 287 """
288 288 buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
289 289 buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
290 290 buffer = numpy.sum(buffer,3)
291 291
292 292 else:
293 293 buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
294 294 buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
295 295 buffer = numpy.sum(buffer,2)
296 296
297 297 dataOut.data = buffer
298 298 dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
299 299 dataOut.windowOfFilter = window
300 300
301 301 return dataOut
302 302
303 303
304 304 class setH0(Operation):
305 305
306 306 def run(self, dataOut, h0, deltaHeight = None):
307 307
308 308 if not deltaHeight:
309 309 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
310 310
311 311 nHeights = dataOut.nHeights
312 312
313 313 newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
314 314
315 315 dataOut.heightList = newHeiRange
316 316
317 317 return dataOut
318 318
319 319
320 320 class deFlip(Operation):
321 321
322 322 def run(self, dataOut, channelList = []):
323 323
324 324 data = dataOut.data.copy()
325 325
326 326 if dataOut.flagDataAsBlock:
327 327 flip = self.flip
328 328 profileList = list(range(dataOut.nProfiles))
329 329
330 330 if not channelList:
331 331 for thisProfile in profileList:
332 332 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
333 333 flip *= -1.0
334 334 else:
335 335 for thisChannel in channelList:
336 336 if thisChannel not in dataOut.channelList:
337 337 continue
338 338
339 339 for thisProfile in profileList:
340 340 data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
341 341 flip *= -1.0
342 342
343 343 self.flip = flip
344 344
345 345 else:
346 346 if not channelList:
347 347 data[:,:] = data[:,:]*self.flip
348 348 else:
349 349 for thisChannel in channelList:
350 350 if thisChannel not in dataOut.channelList:
351 351 continue
352 352
353 353 data[thisChannel,:] = data[thisChannel,:]*self.flip
354 354
355 355 self.flip *= -1.
356 356
357 357 dataOut.data = data
358 358
359 359 return dataOut
360 360
361 361
362 362 class setAttribute(Operation):
363 363 '''
364 364 Set an arbitrary attribute(s) to dataOut
365 365 '''
366 366
367 367 def __init__(self):
368 368
369 369 Operation.__init__(self)
370 370 self._ready = False
371 371
372 372 def run(self, dataOut, **kwargs):
373 373
374 374 for key, value in kwargs.items():
375 375 setattr(dataOut, key, value)
376 376
377 377 return dataOut
378 378
379 379
380 @MPDecorator
381 class printAttribute(Operation):
382 '''
383 Print an arbitrary attribute of dataOut
384 '''
385
386 def __init__(self):
387
388 Operation.__init__(self)
389
390 def run(self, dataOut, attributes):
391
392 for attr in attributes:
393 if hasattr(dataOut, attr):
394 log.log(getattr(dataOut, attr), attr)
395
396
380 397 class interpolateHeights(Operation):
381 398
382 399 def run(self, dataOut, topLim, botLim):
383 400 #69 al 72 para julia
384 401 #82-84 para meteoros
385 402 if len(numpy.shape(dataOut.data))==2:
386 403 sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
387 404 sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
388 405 #dataOut.data[:,botLim:limSup+1] = sampInterp
389 406 dataOut.data[:,botLim:topLim+1] = sampInterp
390 407 else:
391 408 nHeights = dataOut.data.shape[2]
392 409 x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
393 410 y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
394 411 f = interpolate.interp1d(x, y, axis = 2)
395 412 xnew = numpy.arange(botLim,topLim+1)
396 413 ynew = f(xnew)
397 414 dataOut.data[:,:,botLim:topLim+1] = ynew
398 415
399 416 return dataOut
400 417
401 418
402 419 class CohInt(Operation):
403 420
404 421 isConfig = False
405 422 __profIndex = 0
406 423 __byTime = False
407 424 __initime = None
408 425 __lastdatatime = None
409 426 __integrationtime = None
410 427 __buffer = None
411 428 __bufferStride = []
412 429 __dataReady = False
413 430 __profIndexStride = 0
414 431 __dataToPutStride = False
415 432 n = None
416 433
417 434 def __init__(self, **kwargs):
418 435
419 436 Operation.__init__(self, **kwargs)
420 437
421 438 # self.isConfig = False
422 439
423 440 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
424 441 """
425 442 Set the parameters of the integration class.
426 443
427 444 Inputs:
428 445
429 446 n : Number of coherent integrations
430 447 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
431 448 overlapping :
432 449 """
433 450
434 451 self.__initime = None
435 452 self.__lastdatatime = 0
436 453 self.__buffer = None
437 454 self.__dataReady = False
438 455 self.byblock = byblock
439 456 self.stride = stride
440 457
441 458 if n == None and timeInterval == None:
442 459 raise ValueError("n or timeInterval should be specified ...")
443 460
444 461 if n != None:
445 462 self.n = n
446 463 self.__byTime = False
447 464 else:
448 465 self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
449 466 self.n = 9999
450 467 self.__byTime = True
451 468
452 469 if overlapping:
453 470 self.__withOverlapping = True
454 471 self.__buffer = None
455 472 else:
456 473 self.__withOverlapping = False
457 474 self.__buffer = 0
458 475
459 476 self.__profIndex = 0
460 477
461 478 def putData(self, data):
462 479
463 480 """
464 481 Add a profile to the __buffer and increase in one the __profileIndex
465 482
466 483 """
467 484
468 485 if not self.__withOverlapping:
469 486 self.__buffer += data.copy()
470 487 self.__profIndex += 1
471 488 return
472 489
473 490 #Overlapping data
474 491 nChannels, nHeis = data.shape
475 492 data = numpy.reshape(data, (1, nChannels, nHeis))
476 493
477 494 #If the buffer is empty then it takes the data value
478 495 if self.__buffer is None:
479 496 self.__buffer = data
480 497 self.__profIndex += 1
481 498 return
482 499
483 500 #If the buffer length is lower than n then stakcing the data value
484 501 if self.__profIndex < self.n:
485 502 self.__buffer = numpy.vstack((self.__buffer, data))
486 503 self.__profIndex += 1
487 504 return
488 505
489 506 #If the buffer length is equal to n then replacing the last buffer value with the data value
490 507 self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
491 508 self.__buffer[self.n-1] = data
492 509 self.__profIndex = self.n
493 510 return
494 511
495 512
496 513 def pushData(self):
497 514 """
498 515 Return the sum of the last profiles and the profiles used in the sum.
499 516
500 517 Affected:
501 518
502 519 self.__profileIndex
503 520
504 521 """
505 522
506 523 if not self.__withOverlapping:
507 524 data = self.__buffer
508 525 n = self.__profIndex
509 526
510 527 self.__buffer = 0
511 528 self.__profIndex = 0
512 529
513 530 return data, n
514 531
515 532 #Integration with Overlapping
516 533 data = numpy.sum(self.__buffer, axis=0)
517 534 # print data
518 535 # raise
519 536 n = self.__profIndex
520 537
521 538 return data, n
522 539
523 540 def byProfiles(self, data):
524 541
525 542 self.__dataReady = False
526 543 avgdata = None
527 544 # n = None
528 545 # print data
529 546 # raise
530 547 self.putData(data)
531 548
532 549 if self.__profIndex == self.n:
533 550 avgdata, n = self.pushData()
534 551 self.__dataReady = True
535 552
536 553 return avgdata
537 554
538 555 def byTime(self, data, datatime):
539 556
540 557 self.__dataReady = False
541 558 avgdata = None
542 559 n = None
543 560
544 561 self.putData(data)
545 562
546 563 if (datatime - self.__initime) >= self.__integrationtime:
547 564 avgdata, n = self.pushData()
548 565 self.n = n
549 566 self.__dataReady = True
550 567
551 568 return avgdata
552 569
553 570 def integrateByStride(self, data, datatime):
554 571 # print data
555 572 if self.__profIndex == 0:
556 573 self.__buffer = [[data.copy(), datatime]]
557 574 else:
558 575 self.__buffer.append([data.copy(),datatime])
559 576 self.__profIndex += 1
560 577 self.__dataReady = False
561 578
562 579 if self.__profIndex == self.n * self.stride :
563 580 self.__dataToPutStride = True
564 581 self.__profIndexStride = 0
565 582 self.__profIndex = 0
566 583 self.__bufferStride = []
567 584 for i in range(self.stride):
568 585 current = self.__buffer[i::self.stride]
569 586 data = numpy.sum([t[0] for t in current], axis=0)
570 587 avgdatatime = numpy.average([t[1] for t in current])
571 588 # print data
572 589 self.__bufferStride.append((data, avgdatatime))
573 590
574 591 if self.__dataToPutStride:
575 592 self.__dataReady = True
576 593 self.__profIndexStride += 1
577 594 if self.__profIndexStride == self.stride:
578 595 self.__dataToPutStride = False
579 596 # print self.__bufferStride[self.__profIndexStride - 1]
580 597 # raise
581 598 return self.__bufferStride[self.__profIndexStride - 1]
582 599
583 600
584 601 return None, None
585 602
586 603 def integrate(self, data, datatime=None):
587 604
588 605 if self.__initime == None:
589 606 self.__initime = datatime
590 607
591 608 if self.__byTime:
592 609 avgdata = self.byTime(data, datatime)
593 610 else:
594 611 avgdata = self.byProfiles(data)
595 612
596 613
597 614 self.__lastdatatime = datatime
598 615
599 616 if avgdata is None:
600 617 return None, None
601 618
602 619 avgdatatime = self.__initime
603 620
604 621 deltatime = datatime - self.__lastdatatime
605 622
606 623 if not self.__withOverlapping:
607 624 self.__initime = datatime
608 625 else:
609 626 self.__initime += deltatime
610 627
611 628 return avgdata, avgdatatime
612 629
613 630 def integrateByBlock(self, dataOut):
614 631
615 632 times = int(dataOut.data.shape[1]/self.n)
616 633 avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
617 634
618 635 id_min = 0
619 636 id_max = self.n
620 637
621 638 for i in range(times):
622 639 junk = dataOut.data[:,id_min:id_max,:]
623 640 avgdata[:,i,:] = junk.sum(axis=1)
624 641 id_min += self.n
625 642 id_max += self.n
626 643
627 644 timeInterval = dataOut.ippSeconds*self.n
628 645 avgdatatime = (times - 1) * timeInterval + dataOut.utctime
629 646 self.__dataReady = True
630 647 return avgdata, avgdatatime
631 648
632 649 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
633 650
634 651 if not self.isConfig:
635 652 self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
636 653 self.isConfig = True
637 654
638 655 if dataOut.flagDataAsBlock:
639 656 """
640 657 Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
641 658 """
642 659 avgdata, avgdatatime = self.integrateByBlock(dataOut)
643 660 dataOut.nProfiles /= self.n
644 661 else:
645 662 if stride is None:
646 663 avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
647 664 else:
648 665 avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
649 666
650 667
651 668 # dataOut.timeInterval *= n
652 669 dataOut.flagNoData = True
653 670
654 671 if self.__dataReady:
655 672 dataOut.data = avgdata
656 673 dataOut.nCohInt *= self.n
657 674 dataOut.utctime = avgdatatime
658 675 # print avgdata, avgdatatime
659 676 # raise
660 677 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
661 678 dataOut.flagNoData = False
662 679 return dataOut
663 680
664 681 class Decoder(Operation):
665 682
666 683 isConfig = False
667 684 __profIndex = 0
668 685
669 686 code = None
670 687
671 688 nCode = None
672 689 nBaud = None
673 690
674 691 def __init__(self, **kwargs):
675 692
676 693 Operation.__init__(self, **kwargs)
677 694
678 695 self.times = None
679 696 self.osamp = None
680 697 # self.__setValues = False
681 698 self.isConfig = False
682 699 self.setupReq = False
683 700 def setup(self, code, osamp, dataOut):
684 701
685 702 self.__profIndex = 0
686 703
687 704 self.code = code
688 705
689 706 self.nCode = len(code)
690 707 self.nBaud = len(code[0])
691 708
692 709 if (osamp != None) and (osamp >1):
693 710 self.osamp = osamp
694 711 self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
695 712 self.nBaud = self.nBaud*self.osamp
696 713
697 714 self.__nChannels = dataOut.nChannels
698 715 self.__nProfiles = dataOut.nProfiles
699 716 self.__nHeis = dataOut.nHeights
700 717
701 718 if self.__nHeis < self.nBaud:
702 719 raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
703 720
704 721 #Frequency
705 722 __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
706 723
707 724 __codeBuffer[:,0:self.nBaud] = self.code
708 725
709 726 self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
710 727
711 728 if dataOut.flagDataAsBlock:
712 729
713 730 self.ndatadec = self.__nHeis #- self.nBaud + 1
714 731
715 732 self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
716 733
717 734 else:
718 735
719 736 #Time
720 737 self.ndatadec = self.__nHeis #- self.nBaud + 1
721 738
722 739 self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
723 740
724 741 def __convolutionInFreq(self, data):
725 742
726 743 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
727 744
728 745 fft_data = numpy.fft.fft(data, axis=1)
729 746
730 747 conv = fft_data*fft_code
731 748
732 749 data = numpy.fft.ifft(conv,axis=1)
733 750
734 751 return data
735 752
736 753 def __convolutionInFreqOpt(self, data):
737 754
738 755 raise NotImplementedError
739 756
740 757 def __convolutionInTime(self, data):
741 758
742 759 code = self.code[self.__profIndex]
743 760 for i in range(self.__nChannels):
744 761 self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
745 762
746 763 return self.datadecTime
747 764
748 765 def __convolutionByBlockInTime(self, data):
749 766
750 767 repetitions = int(self.__nProfiles / self.nCode)
751 768 junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
752 769 junk = junk.flatten()
753 770 code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
754 771 profilesList = range(self.__nProfiles)
755 772
756 773 for i in range(self.__nChannels):
757 774 for j in profilesList:
758 775 self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
759 776 return self.datadecTime
760 777
761 778 def __convolutionByBlockInFreq(self, data):
762 779
763 780 raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
764 781
765 782
766 783 fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
767 784
768 785 fft_data = numpy.fft.fft(data, axis=2)
769 786
770 787 conv = fft_data*fft_code
771 788
772 789 data = numpy.fft.ifft(conv,axis=2)
773 790
774 791 return data
775 792
776 793
777 794 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
778 795
779 796 if dataOut.flagDecodeData:
780 797 print("This data is already decoded, recoding again ...")
781 798
782 799 if not self.isConfig:
783 800
784 801 if code is None:
785 802 if dataOut.code is None:
786 803 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
787 804
788 805 code = dataOut.code
789 806 else:
790 807 code = numpy.array(code).reshape(nCode,nBaud)
791 808 self.setup(code, osamp, dataOut)
792 809
793 810 self.isConfig = True
794 811
795 812 if mode == 3:
796 813 sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
797 814
798 815 if times != None:
799 816 sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
800 817
801 818 if self.code is None:
802 819 print("Fail decoding: Code is not defined.")
803 820 return
804 821
805 822 self.__nProfiles = dataOut.nProfiles
806 823 datadec = None
807 824
808 825 if mode == 3:
809 826 mode = 0
810 827
811 828 if dataOut.flagDataAsBlock:
812 829 """
813 830 Decoding when data have been read as block,
814 831 """
815 832
816 833 if mode == 0:
817 834 datadec = self.__convolutionByBlockInTime(dataOut.data)
818 835 if mode == 1:
819 836 datadec = self.__convolutionByBlockInFreq(dataOut.data)
820 837 else:
821 838 """
822 839 Decoding when data have been read profile by profile
823 840 """
824 841 if mode == 0:
825 842 datadec = self.__convolutionInTime(dataOut.data)
826 843
827 844 if mode == 1:
828 845 datadec = self.__convolutionInFreq(dataOut.data)
829 846
830 847 if mode == 2:
831 848 datadec = self.__convolutionInFreqOpt(dataOut.data)
832 849
833 850 if datadec is None:
834 851 raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
835 852
836 853 dataOut.code = self.code
837 854 dataOut.nCode = self.nCode
838 855 dataOut.nBaud = self.nBaud
839 856
840 857 dataOut.data = datadec
841 858
842 859 dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
843 860
844 861 dataOut.flagDecodeData = True #asumo q la data esta decodificada
845 862
846 863 if self.__profIndex == self.nCode-1:
847 864 self.__profIndex = 0
848 865 return dataOut
849 866
850 867 self.__profIndex += 1
851 868
852 869 return dataOut
853 870 # dataOut.flagDeflipData = True #asumo q la data no esta sin flip
854 871
855 872
856 873 class ProfileConcat(Operation):
857 874
858 875 isConfig = False
859 876 buffer = None
860 877
861 878 def __init__(self, **kwargs):
862 879
863 880 Operation.__init__(self, **kwargs)
864 881 self.profileIndex = 0
865 882
866 883 def reset(self):
867 884 self.buffer = numpy.zeros_like(self.buffer)
868 885 self.start_index = 0
869 886 self.times = 1
870 887
871 888 def setup(self, data, m, n=1):
872 889 self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
873 890 self.nHeights = data.shape[1]#.nHeights
874 891 self.start_index = 0
875 892 self.times = 1
876 893
877 894 def concat(self, data):
878 895
879 896 self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
880 897 self.start_index = self.start_index + self.nHeights
881 898
882 899 def run(self, dataOut, m):
883 900 dataOut.flagNoData = True
884 901
885 902 if not self.isConfig:
886 903 self.setup(dataOut.data, m, 1)
887 904 self.isConfig = True
888 905
889 906 if dataOut.flagDataAsBlock:
890 907 raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
891 908
892 909 else:
893 910 self.concat(dataOut.data)
894 911 self.times += 1
895 912 if self.times > m:
896 913 dataOut.data = self.buffer
897 914 self.reset()
898 915 dataOut.flagNoData = False
899 916 # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
900 917 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
901 918 xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
902 919 dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
903 920 dataOut.ippSeconds *= m
904 921 return dataOut
905 922
906 923 class ProfileSelector(Operation):
907 924
908 925 profileIndex = None
909 926 # Tamanho total de los perfiles
910 927 nProfiles = None
911 928
912 929 def __init__(self, **kwargs):
913 930
914 931 Operation.__init__(self, **kwargs)
915 932 self.profileIndex = 0
916 933
917 934 def incProfileIndex(self):
918 935
919 936 self.profileIndex += 1
920 937
921 938 if self.profileIndex >= self.nProfiles:
922 939 self.profileIndex = 0
923 940
924 941 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
925 942
926 943 if profileIndex < minIndex:
927 944 return False
928 945
929 946 if profileIndex > maxIndex:
930 947 return False
931 948
932 949 return True
933 950
934 951 def isThisProfileInList(self, profileIndex, profileList):
935 952
936 953 if profileIndex not in profileList:
937 954 return False
938 955
939 956 return True
940 957
941 958 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
942 959
943 960 """
944 961 ProfileSelector:
945 962
946 963 Inputs:
947 964 profileList : Index of profiles selected. Example: profileList = (0,1,2,7,8)
948 965
949 966 profileRangeList : Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
950 967
951 968 rangeList : List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
952 969
953 970 """
954 971
955 972 if rangeList is not None:
956 973 if type(rangeList[0]) not in (tuple, list):
957 974 rangeList = [rangeList]
958 975
959 976 dataOut.flagNoData = True
960 977
961 978 if dataOut.flagDataAsBlock:
962 979 """
963 980 data dimension = [nChannels, nProfiles, nHeis]
964 981 """
965 982 if profileList != None:
966 983 dataOut.data = dataOut.data[:,profileList,:]
967 984
968 985 if profileRangeList != None:
969 986 minIndex = profileRangeList[0]
970 987 maxIndex = profileRangeList[1]
971 988 profileList = list(range(minIndex, maxIndex+1))
972 989
973 990 dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
974 991
975 992 if rangeList != None:
976 993
977 994 profileList = []
978 995
979 996 for thisRange in rangeList:
980 997 minIndex = thisRange[0]
981 998 maxIndex = thisRange[1]
982 999
983 1000 profileList.extend(list(range(minIndex, maxIndex+1)))
984 1001
985 1002 dataOut.data = dataOut.data[:,profileList,:]
986 1003
987 1004 dataOut.nProfiles = len(profileList)
988 1005 dataOut.profileIndex = dataOut.nProfiles - 1
989 1006 dataOut.flagNoData = False
990 1007
991 1008 return dataOut
992 1009
993 1010 """
994 1011 data dimension = [nChannels, nHeis]
995 1012 """
996 1013
997 1014 if profileList != None:
998 1015
999 1016 if self.isThisProfileInList(dataOut.profileIndex, profileList):
1000 1017
1001 1018 self.nProfiles = len(profileList)
1002 1019 dataOut.nProfiles = self.nProfiles
1003 1020 dataOut.profileIndex = self.profileIndex
1004 1021 dataOut.flagNoData = False
1005 1022
1006 1023 self.incProfileIndex()
1007 1024 return dataOut
1008 1025
1009 1026 if profileRangeList != None:
1010 1027
1011 1028 minIndex = profileRangeList[0]
1012 1029 maxIndex = profileRangeList[1]
1013 1030
1014 1031 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1015 1032
1016 1033 self.nProfiles = maxIndex - minIndex + 1
1017 1034 dataOut.nProfiles = self.nProfiles
1018 1035 dataOut.profileIndex = self.profileIndex
1019 1036 dataOut.flagNoData = False
1020 1037
1021 1038 self.incProfileIndex()
1022 1039 return dataOut
1023 1040
1024 1041 if rangeList != None:
1025 1042
1026 1043 nProfiles = 0
1027 1044
1028 1045 for thisRange in rangeList:
1029 1046 minIndex = thisRange[0]
1030 1047 maxIndex = thisRange[1]
1031 1048
1032 1049 nProfiles += maxIndex - minIndex + 1
1033 1050
1034 1051 for thisRange in rangeList:
1035 1052
1036 1053 minIndex = thisRange[0]
1037 1054 maxIndex = thisRange[1]
1038 1055
1039 1056 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
1040 1057
1041 1058 self.nProfiles = nProfiles
1042 1059 dataOut.nProfiles = self.nProfiles
1043 1060 dataOut.profileIndex = self.profileIndex
1044 1061 dataOut.flagNoData = False
1045 1062
1046 1063 self.incProfileIndex()
1047 1064
1048 1065 break
1049 1066
1050 1067 return dataOut
1051 1068
1052 1069
1053 1070 if beam != None: #beam is only for AMISR data
1054 1071 if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
1055 1072 dataOut.flagNoData = False
1056 1073 dataOut.profileIndex = self.profileIndex
1057 1074
1058 1075 self.incProfileIndex()
1059 1076
1060 1077 return dataOut
1061 1078
1062 1079 raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
1063 1080
1064 1081
1065 1082 class Reshaper(Operation):
1066 1083
1067 1084 def __init__(self, **kwargs):
1068 1085
1069 1086 Operation.__init__(self, **kwargs)
1070 1087
1071 1088 self.__buffer = None
1072 1089 self.__nitems = 0
1073 1090
1074 1091 def __appendProfile(self, dataOut, nTxs):
1075 1092
1076 1093 if self.__buffer is None:
1077 1094 shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
1078 1095 self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
1079 1096
1080 1097 ini = dataOut.nHeights * self.__nitems
1081 1098 end = ini + dataOut.nHeights
1082 1099
1083 1100 self.__buffer[:, ini:end] = dataOut.data
1084 1101
1085 1102 self.__nitems += 1
1086 1103
1087 1104 return int(self.__nitems*nTxs)
1088 1105
1089 1106 def __getBuffer(self):
1090 1107
1091 1108 if self.__nitems == int(1./self.__nTxs):
1092 1109
1093 1110 self.__nitems = 0
1094 1111
1095 1112 return self.__buffer.copy()
1096 1113
1097 1114 return None
1098 1115
1099 1116 def __checkInputs(self, dataOut, shape, nTxs):
1100 1117
1101 1118 if shape is None and nTxs is None:
1102 1119 raise ValueError("Reshaper: shape of factor should be defined")
1103 1120
1104 1121 if nTxs:
1105 1122 if nTxs < 0:
1106 1123 raise ValueError("nTxs should be greater than 0")
1107 1124
1108 1125 if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
1109 1126 raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
1110 1127
1111 1128 shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
1112 1129
1113 1130 return shape, nTxs
1114 1131
1115 1132 if len(shape) != 2 and len(shape) != 3:
1116 1133 raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
1117 1134
1118 1135 if len(shape) == 2:
1119 1136 shape_tuple = [dataOut.nChannels]
1120 1137 shape_tuple.extend(shape)
1121 1138 else:
1122 1139 shape_tuple = list(shape)
1123 1140
1124 1141 nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
1125 1142
1126 1143 return shape_tuple, nTxs
1127 1144
1128 1145 def run(self, dataOut, shape=None, nTxs=None):
1129 1146
1130 1147 shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
1131 1148
1132 1149 dataOut.flagNoData = True
1133 1150 profileIndex = None
1134 1151
1135 1152 if dataOut.flagDataAsBlock:
1136 1153
1137 1154 dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
1138 1155 dataOut.flagNoData = False
1139 1156
1140 1157 profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
1141 1158
1142 1159 else:
1143 1160
1144 1161 if self.__nTxs < 1:
1145 1162
1146 1163 self.__appendProfile(dataOut, self.__nTxs)
1147 1164 new_data = self.__getBuffer()
1148 1165
1149 1166 if new_data is not None:
1150 1167 dataOut.data = new_data
1151 1168 dataOut.flagNoData = False
1152 1169
1153 1170 profileIndex = dataOut.profileIndex*nTxs
1154 1171
1155 1172 else:
1156 1173 raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
1157 1174
1158 1175 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1159 1176
1160 1177 dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
1161 1178
1162 1179 dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
1163 1180
1164 1181 dataOut.profileIndex = profileIndex
1165 1182
1166 1183 dataOut.ippSeconds /= self.__nTxs
1167 1184
1168 1185 return dataOut
1169 1186
1170 1187 class SplitProfiles(Operation):
1171 1188
1172 1189 def __init__(self, **kwargs):
1173 1190
1174 1191 Operation.__init__(self, **kwargs)
1175 1192
1176 1193 def run(self, dataOut, n):
1177 1194
1178 1195 dataOut.flagNoData = True
1179 1196 profileIndex = None
1180 1197
1181 1198 if dataOut.flagDataAsBlock:
1182 1199
1183 1200 #nchannels, nprofiles, nsamples
1184 1201 shape = dataOut.data.shape
1185 1202
1186 1203 if shape[2] % n != 0:
1187 1204 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
1188 1205
1189 1206 new_shape = shape[0], shape[1]*n, int(shape[2]/n)
1190 1207
1191 1208 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1192 1209 dataOut.flagNoData = False
1193 1210
1194 1211 profileIndex = int(dataOut.nProfiles/n) - 1
1195 1212
1196 1213 else:
1197 1214
1198 1215 raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
1199 1216
1200 1217 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1201 1218
1202 1219 dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
1203 1220
1204 1221 dataOut.nProfiles = int(dataOut.nProfiles*n)
1205 1222
1206 1223 dataOut.profileIndex = profileIndex
1207 1224
1208 1225 dataOut.ippSeconds /= n
1209 1226
1210 1227 return dataOut
1211 1228
1212 1229 class CombineProfiles(Operation):
1213 1230 def __init__(self, **kwargs):
1214 1231
1215 1232 Operation.__init__(self, **kwargs)
1216 1233
1217 1234 self.__remData = None
1218 1235 self.__profileIndex = 0
1219 1236
1220 1237 def run(self, dataOut, n):
1221 1238
1222 1239 dataOut.flagNoData = True
1223 1240 profileIndex = None
1224 1241
1225 1242 if dataOut.flagDataAsBlock:
1226 1243
1227 1244 #nchannels, nprofiles, nsamples
1228 1245 shape = dataOut.data.shape
1229 1246 new_shape = shape[0], shape[1]/n, shape[2]*n
1230 1247
1231 1248 if shape[1] % n != 0:
1232 1249 raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
1233 1250
1234 1251 dataOut.data = numpy.reshape(dataOut.data, new_shape)
1235 1252 dataOut.flagNoData = False
1236 1253
1237 1254 profileIndex = int(dataOut.nProfiles*n) - 1
1238 1255
1239 1256 else:
1240 1257
1241 1258 #nchannels, nsamples
1242 1259 if self.__remData is None:
1243 1260 newData = dataOut.data
1244 1261 else:
1245 1262 newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
1246 1263
1247 1264 self.__profileIndex += 1
1248 1265
1249 1266 if self.__profileIndex < n:
1250 1267 self.__remData = newData
1251 1268 #continue
1252 1269 return
1253 1270
1254 1271 self.__profileIndex = 0
1255 1272 self.__remData = None
1256 1273
1257 1274 dataOut.data = newData
1258 1275 dataOut.flagNoData = False
1259 1276
1260 1277 profileIndex = dataOut.profileIndex/n
1261 1278
1262 1279
1263 1280 deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1264 1281
1265 1282 dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
1266 1283
1267 1284 dataOut.nProfiles = int(dataOut.nProfiles/n)
1268 1285
1269 1286 dataOut.profileIndex = profileIndex
1270 1287
1271 1288 dataOut.ippSeconds *= n
1272 1289
1273 1290 return dataOut
1274 1291
1275 1292 class PulsePairVoltage(Operation):
1276 1293 '''
1277 1294 Function PulsePair(Signal Power, Velocity)
1278 1295 The real component of Lag[0] provides Intensity Information
1279 1296 The imag component of Lag[1] Phase provides Velocity Information
1280 1297
1281 1298 Configuration Parameters:
1282 1299 nPRF = Number of Several PRF
1283 1300 theta = Degree Azimuth angel Boundaries
1284 1301
1285 1302 Input:
1286 1303 self.dataOut
1287 1304 lag[N]
1288 1305 Affected:
1289 1306 self.dataOut.spc
1290 1307 '''
1291 1308 isConfig = False
1292 1309 __profIndex = 0
1293 1310 __initime = None
1294 1311 __lastdatatime = None
1295 1312 __buffer = None
1296 1313 __buffer2 = []
1297 1314 __buffer3 = None
1298 1315 __dataReady = False
1299 1316 n = None
1300 1317 __nch = 0
1301 1318 __nHeis = 0
1302 1319 removeDC = False
1303 1320 ipp = None
1304 1321 lambda_ = 0
1305 1322
1306 1323 def __init__(self,**kwargs):
1307 1324 Operation.__init__(self,**kwargs)
1308 1325
1309 1326 def setup(self, dataOut, n = None, removeDC=False):
1310 1327 '''
1311 1328 n= Numero de PRF's de entrada
1312 1329 '''
1313 1330 self.__initime = None
1314 1331 self.__lastdatatime = 0
1315 1332 self.__dataReady = False
1316 1333 self.__buffer = 0
1317 1334 self.__buffer2 = []
1318 1335 self.__buffer3 = 0
1319 1336 self.__profIndex = 0
1320 1337
1321 1338 self.__nch = dataOut.nChannels
1322 1339 self.__nHeis = dataOut.nHeights
1323 1340 self.removeDC = removeDC
1324 1341 self.lambda_ = 3.0e8/(9345.0e6)
1325 1342 self.ippSec = dataOut.ippSeconds
1326 1343 self.nCohInt = dataOut.nCohInt
1327 1344 print("IPPseconds",dataOut.ippSeconds)
1328 1345
1329 1346 print("ELVALOR DE n es:", n)
1330 1347 if n == None:
1331 1348 raise ValueError("n should be specified.")
1332 1349
1333 1350 if n != None:
1334 1351 if n<2:
1335 1352 raise ValueError("n should be greater than 2")
1336 1353
1337 1354 self.n = n
1338 1355 self.__nProf = n
1339 1356
1340 1357 self.__buffer = numpy.zeros((dataOut.nChannels,
1341 1358 n,
1342 1359 dataOut.nHeights),
1343 1360 dtype='complex')
1344 1361
1345 1362 def putData(self,data):
1346 1363 '''
1347 1364 Add a profile to he __buffer and increase in one the __profiel Index
1348 1365 '''
1349 1366 self.__buffer[:,self.__profIndex,:]= data
1350 1367 self.__profIndex += 1
1351 1368 return
1352 1369
1353 1370 def pushData(self):
1354 1371 '''
1355 1372 Return the PULSEPAIR and the profiles used in the operation
1356 1373 Affected : self.__profileIndex
1357 1374 '''
1358 1375
1359 1376 if self.removeDC==True:
1360 1377 mean = numpy.mean(self.__buffer,1)
1361 1378 tmp = mean.reshape(self.__nch,1,self.__nHeis)
1362 1379 dc= numpy.tile(tmp,[1,self.__nProf,1])
1363 1380 self.__buffer = self.__buffer - dc
1364 1381
1365 1382 data_intensity = numpy.sum(self.__buffer*numpy.conj(self.__buffer),1)/(self.n*self.nCohInt)#*self.nCohInt)
1366 1383 pair1 = self.__buffer[:,1:,:]*numpy.conjugate(self.__buffer[:,:-1,:])
1367 1384 angle = numpy.angle(numpy.sum(pair1,1))*180/(math.pi)
1368 1385 #print(angle.shape)#print("__ANGLE__") #print("angle",angle[:,:10])
1369 1386 data_velocity = (self.lambda_/(4*math.pi*self.ippSec))*numpy.angle(numpy.sum(pair1,1))#self.ippSec*self.nCohInt
1370 1387 n = self.__profIndex
1371 1388
1372 1389 self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex')
1373 1390 self.__profIndex = 0
1374 1391 return data_intensity,data_velocity,n
1375 1392
1376 1393 def pulsePairbyProfiles(self,data):
1377 1394
1378 1395 self.__dataReady = False
1379 1396 data_intensity = None
1380 1397 data_velocity = None
1381 1398 self.putData(data)
1382 1399 if self.__profIndex == self.n:
1383 1400 data_intensity, data_velocity, n = self.pushData()
1384 1401 self.__dataReady = True
1385 1402
1386 1403 return data_intensity, data_velocity
1387 1404
1388 1405 def pulsePairOp(self, data, datatime= None):
1389 1406
1390 1407 if self.__initime == None:
1391 1408 self.__initime = datatime
1392 1409
1393 1410 data_intensity, data_velocity = self.pulsePairbyProfiles(data)
1394 1411 self.__lastdatatime = datatime
1395 1412
1396 1413 if data_intensity is None:
1397 1414 return None, None, None
1398 1415
1399 1416 avgdatatime = self.__initime
1400 1417 deltatime = datatime - self.__lastdatatime
1401 1418 self.__initime = datatime
1402 1419
1403 1420 return data_intensity, data_velocity, avgdatatime
1404 1421
1405 1422 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
1406 1423
1407 1424 if not self.isConfig:
1408 1425 self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs)
1409 1426 self.isConfig = True
1410 1427 data_intensity, data_velocity, avgdatatime = self.pulsePairOp(dataOut.data, dataOut.utctime)
1411 1428 dataOut.flagNoData = True
1412 1429
1413 1430 if self.__dataReady:
1414 1431 dataOut.nCohInt *= self.n
1415 1432 dataOut.data_intensity = data_intensity #valor para intensidad
1416 1433 dataOut.data_velocity = data_velocity #valor para velocidad
1417 1434 dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo.
1418 1435 dataOut.utctime = avgdatatime
1419 1436 dataOut.flagNoData = False
1420 1437 return dataOut
1421 1438
1422 1439
1423 1440 # import collections
1424 1441 # from scipy.stats import mode
1425 1442 #
1426 1443 # class Synchronize(Operation):
1427 1444 #
1428 1445 # isConfig = False
1429 1446 # __profIndex = 0
1430 1447 #
1431 1448 # def __init__(self, **kwargs):
1432 1449 #
1433 1450 # Operation.__init__(self, **kwargs)
1434 1451 # # self.isConfig = False
1435 1452 # self.__powBuffer = None
1436 1453 # self.__startIndex = 0
1437 1454 # self.__pulseFound = False
1438 1455 #
1439 1456 # def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
1440 1457 #
1441 1458 # #Read data
1442 1459 #
1443 1460 # powerdB = dataOut.getPower(channel = channel)
1444 1461 # noisedB = dataOut.getNoise(channel = channel)[0]
1445 1462 #
1446 1463 # self.__powBuffer.extend(powerdB.flatten())
1447 1464 #
1448 1465 # dataArray = numpy.array(self.__powBuffer)
1449 1466 #
1450 1467 # filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
1451 1468 #
1452 1469 # maxValue = numpy.nanmax(filteredPower)
1453 1470 #
1454 1471 # if maxValue < noisedB + 10:
1455 1472 # #No se encuentra ningun pulso de transmision
1456 1473 # return None
1457 1474 #
1458 1475 # maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
1459 1476 #
1460 1477 # if len(maxValuesIndex) < 2:
1461 1478 # #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
1462 1479 # return None
1463 1480 #
1464 1481 # phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
1465 1482 #
1466 1483 # #Seleccionar solo valores con un espaciamiento de nSamples
1467 1484 # pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
1468 1485 #
1469 1486 # if len(pulseIndex) < 2:
1470 1487 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1471 1488 # return None
1472 1489 #
1473 1490 # spacing = pulseIndex[1:] - pulseIndex[:-1]
1474 1491 #
1475 1492 # #remover senales que se distancien menos de 10 unidades o muestras
1476 1493 # #(No deberian existir IPP menor a 10 unidades)
1477 1494 #
1478 1495 # realIndex = numpy.where(spacing > 10 )[0]
1479 1496 #
1480 1497 # if len(realIndex) < 2:
1481 1498 # #Solo se encontro un pulso de transmision con ancho mayor a 1
1482 1499 # return None
1483 1500 #
1484 1501 # #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
1485 1502 # realPulseIndex = pulseIndex[realIndex]
1486 1503 #
1487 1504 # period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
1488 1505 #
1489 1506 # print "IPP = %d samples" %period
1490 1507 #
1491 1508 # self.__newNSamples = dataOut.nHeights #int(period)
1492 1509 # self.__startIndex = int(realPulseIndex[0])
1493 1510 #
1494 1511 # return 1
1495 1512 #
1496 1513 #
1497 1514 # def setup(self, nSamples, nChannels, buffer_size = 4):
1498 1515 #
1499 1516 # self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
1500 1517 # maxlen = buffer_size*nSamples)
1501 1518 #
1502 1519 # bufferList = []
1503 1520 #
1504 1521 # for i in range(nChannels):
1505 1522 # bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) + numpy.NAN,
1506 1523 # maxlen = buffer_size*nSamples)
1507 1524 #
1508 1525 # bufferList.append(bufferByChannel)
1509 1526 #
1510 1527 # self.__nSamples = nSamples
1511 1528 # self.__nChannels = nChannels
1512 1529 # self.__bufferList = bufferList
1513 1530 #
1514 1531 # def run(self, dataOut, channel = 0):
1515 1532 #
1516 1533 # if not self.isConfig:
1517 1534 # nSamples = dataOut.nHeights
1518 1535 # nChannels = dataOut.nChannels
1519 1536 # self.setup(nSamples, nChannels)
1520 1537 # self.isConfig = True
1521 1538 #
1522 1539 # #Append new data to internal buffer
1523 1540 # for thisChannel in range(self.__nChannels):
1524 1541 # bufferByChannel = self.__bufferList[thisChannel]
1525 1542 # bufferByChannel.extend(dataOut.data[thisChannel])
1526 1543 #
1527 1544 # if self.__pulseFound:
1528 1545 # self.__startIndex -= self.__nSamples
1529 1546 #
1530 1547 # #Finding Tx Pulse
1531 1548 # if not self.__pulseFound:
1532 1549 # indexFound = self.__findTxPulse(dataOut, channel)
1533 1550 #
1534 1551 # if indexFound == None:
1535 1552 # dataOut.flagNoData = True
1536 1553 # return
1537 1554 #
1538 1555 # self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
1539 1556 # self.__pulseFound = True
1540 1557 # self.__startIndex = indexFound
1541 1558 #
1542 1559 # #If pulse was found ...
1543 1560 # for thisChannel in range(self.__nChannels):
1544 1561 # bufferByChannel = self.__bufferList[thisChannel]
1545 1562 # #print self.__startIndex
1546 1563 # x = numpy.array(bufferByChannel)
1547 1564 # self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
1548 1565 #
1549 1566 # deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
1550 1567 # dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
1551 1568 # # dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
1552 1569 #
1553 1570 # dataOut.data = self.__arrayBuffer
1554 1571 #
1555 1572 # self.__startIndex += self.__newNSamples
1556 1573 #
1557 1574 # return
General Comments 0
You need to be logged in to leave comments. Login now