##// END OF EJS Templates
schain
RSS
añadido FaradayIntegration para limpieza de datos, restringida la operación al funcionamiento con la pdata regular
joabAM -
r1399:101813202a42
parent child
Show More
Side by Side Unified
Context file:
r1399:101813202a42 Loading diff...:
Show file before | Show file after | Hide comments
@@ -1,1576 +1,1577
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 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 startDateTime= datetime.datetime.combine(startDate,startTime)
693 693 endDateTime = datetime.datetime.combine(endDate,endTime)
694 #print("dt eval: ", dt, startDateTime,endDateTime)
694 695 if startDateTime <= dt <= endDateTime:
695 696 return True
696 697 return False
697 698
698 699 def verifyFile(self, filename):
699 700 """Check for a valid file
700 701
701 702 Arguments:
702 703 filename -- full path filename
703 704
704 705 Return:
705 706 boolean
706 707 """
707 708
708 709 return True
709 710
710 711 def checkForRealPath(self, nextFile, nextDay):
711 712 """Check if the next file to be readed exists"""
712 713
713 714 raise NotImplementedError
714 715
715 716 def readFirstHeader(self):
716 717 """Parse the file header"""
717 718
718 719 pass
719 720
720 721 def waitDataBlock(self, pointer_location, blocksize=None):
721 722 """
722 723 """
723 724
724 725 currentPointer = pointer_location
725 726 if blocksize is None:
726 727 neededSize = self.processingHeaderObj.blockSize # + self.basicHeaderSize
727 728 else:
728 729 neededSize = blocksize
729 730
730 731 for nTries in range(self.nTries):
731 732 self.fp.close()
732 733 self.fp = open(self.filename, 'rb')
733 734 self.fp.seek(currentPointer)
734 735
735 736 self.fileSize = os.path.getsize(self.filename)
736 737 currentSize = self.fileSize - currentPointer
737 738
738 739 if (currentSize >= neededSize):
739 740 return 1
740 741
741 742 log.warning(
742 743 "Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1),
743 744 self.name
744 745 )
745 746 time.sleep(self.delay)
746 747
747 748 return 0
748 749
749 750 class JRODataReader(Reader):
750 751
751 752 utc = 0
752 753 nReadBlocks = 0
753 754 foldercounter = 0
754 755 firstHeaderSize = 0
755 756 basicHeaderSize = 24
756 757 __isFirstTimeOnline = 1
757 758 filefmt = "*%Y%j***"
758 759 folderfmt = "*%Y%j"
759 760 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'online', 'delay', 'walk']
760 761
761 762 def getDtypeWidth(self):
762 763
763 764 dtype_index = get_dtype_index(self.dtype)
764 765 dtype_width = get_dtype_width(dtype_index)
765 766
766 767 return dtype_width
767 768
768 769 def checkForRealPath(self, nextFile, nextDay):
769 770 """Check if the next file to be readed exists.
770 771
771 772 Example :
772 773 nombre correcto del file es .../.../D2009307/P2009307367.ext
773 774
774 775 Entonces la funcion prueba con las siguientes combinaciones
775 776 .../.../y2009307367.ext
776 777 .../.../Y2009307367.ext
777 778 .../.../x2009307/y2009307367.ext
778 779 .../.../x2009307/Y2009307367.ext
779 780 .../.../X2009307/y2009307367.ext
780 781 .../.../X2009307/Y2009307367.ext
781 782 siendo para este caso, la ultima combinacion de letras, identica al file buscado
782 783
783 784 Return:
784 785 str -- fullpath of the file
785 786 """
786 787
787 788
788 789 if nextFile:
789 790 self.set += 1
790 791 if nextDay:
791 792 self.set = 0
792 793 self.doy += 1
793 794 foldercounter = 0
794 795 prefixDirList = [None, 'd', 'D']
795 796 if self.ext.lower() == ".r": # voltage
796 797 prefixFileList = ['d', 'D']
797 798 elif self.ext.lower() == ".pdata": # spectra
798 799 prefixFileList = ['p', 'P']
799 800
800 801 # barrido por las combinaciones posibles
801 802 for prefixDir in prefixDirList:
802 803 thispath = self.path
803 804 if prefixDir != None:
804 805 # formo el nombre del directorio xYYYYDDD (x=d o x=D)
805 806 if foldercounter == 0:
806 807 thispath = os.path.join(self.path, "%s%04d%03d" %
807 808 (prefixDir, self.year, self.doy))
808 809 else:
809 810 thispath = os.path.join(self.path, "%s%04d%03d_%02d" % (
810 811 prefixDir, self.year, self.doy, foldercounter))
811 812 for prefixFile in prefixFileList: # barrido por las dos combinaciones posibles de "D"
812 813 # formo el nombre del file xYYYYDDDSSS.ext
813 814 filename = "%s%04d%03d%03d%s" % (prefixFile, self.year, self.doy, self.set, self.ext)
814 815 fullfilename = os.path.join(
815 816 thispath, filename)
816 817
817 818 if os.path.exists(fullfilename):
818 819 return fullfilename, filename
819 820
820 821 return None, filename
821 822
822 823 def __waitNewBlock(self):
823 824 """
824 825 Return 1 si se encontro un nuevo bloque de datos, 0 de otra forma.
825 826
826 827 Si el modo de lectura es OffLine siempre retorn 0
827 828 """
828 829 if not self.online:
829 830 return 0
830 831
831 832 if (self.nReadBlocks >= self.processingHeaderObj.dataBlocksPerFile):
832 833 return 0
833 834
834 835 currentPointer = self.fp.tell()
835 836
836 837 neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
837 838
838 839 for nTries in range(self.nTries):
839 840
840 841 self.fp.close()
841 842 self.fp = open(self.filename, 'rb')
842 843 self.fp.seek(currentPointer)
843 844
844 845 self.fileSize = os.path.getsize(self.filename)
845 846 currentSize = self.fileSize - currentPointer
846 847
847 848 if (currentSize >= neededSize):
848 849 self.basicHeaderObj.read(self.fp)
849 850 return 1
850 851
851 852 if self.fileSize == self.fileSizeByHeader:
852 853 # self.flagEoF = True
853 854 return 0
854 855
855 856 print("[Reading] Waiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries + 1))
856 857 time.sleep(self.delay)
857 858
858 859 return 0
859 860
860 861 def __setNewBlock(self):
861 862
862 863 if self.fp == None:
863 864 return 0
864 865
865 866 if self.flagIsNewFile:
866 867 self.lastUTTime = self.basicHeaderObj.utc
867 868 return 1
868 869
869 870 if self.realtime:
870 871 self.flagDiscontinuousBlock = 1
871 872 if not(self.setNextFile()):
872 873 return 0
873 874 else:
874 875 return 1
875 876
876 877 currentSize = self.fileSize - self.fp.tell()
877 878 neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
878 879
879 880 if (currentSize >= neededSize):
880 881 self.basicHeaderObj.read(self.fp)
881 882 self.lastUTTime = self.basicHeaderObj.utc
882 883 return 1
883 884
884 885 if self.__waitNewBlock():
885 886 self.lastUTTime = self.basicHeaderObj.utc
886 887 return 1
887 888
888 889 if not(self.setNextFile()):
889 890 return 0
890 891
891 892 deltaTime = self.basicHeaderObj.utc - self.lastUTTime
892 893 self.lastUTTime = self.basicHeaderObj.utc
893 894
894 895 self.flagDiscontinuousBlock = 0
895 896
896 897 if deltaTime > self.maxTimeStep:
897 898 self.flagDiscontinuousBlock = 1
898 899
899 900 return 1
900 901
901 902 def readNextBlock(self):
902 903
903 904 while True:
904 905 if not(self.__setNewBlock()):
905 906 continue
906 907
907 908 if not(self.readBlock()):
908 909 return 0
909 910
910 911 self.getBasicHeader()
911 912
912 913 if not self.isDateTimeInRange(self.dataOut.datatime, self.startDate, self.endDate, self.startTime, self.endTime):
913 914 print("[Reading] Block No. %d/%d -> %s [Skipping]" % (self.nReadBlocks,
914 915 self.processingHeaderObj.dataBlocksPerFile,
915 916 self.dataOut.datatime.ctime()))
916 917 continue
917 918
918 919 break
919 920
920 921 if self.verbose:
921 922 print("[Reading] Block No. %d/%d -> %s" % (self.nReadBlocks,
922 923 self.processingHeaderObj.dataBlocksPerFile,
923 924 self.dataOut.datatime.ctime()))
924 925 return 1
925 926
926 927 def readFirstHeader(self):
927 928
928 929 self.basicHeaderObj.read(self.fp)
929 930 self.systemHeaderObj.read(self.fp)
930 931 self.radarControllerHeaderObj.read(self.fp)
931 932 self.processingHeaderObj.read(self.fp)
932 933 self.firstHeaderSize = self.basicHeaderObj.size
933 934
934 935 datatype = int(numpy.log2((self.processingHeaderObj.processFlags &
935 936 PROCFLAG.DATATYPE_MASK)) - numpy.log2(PROCFLAG.DATATYPE_CHAR))
936 937 if datatype == 0:
937 938 datatype_str = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
938 939 elif datatype == 1:
939 940 datatype_str = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
940 941 elif datatype == 2:
941 942 datatype_str = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
942 943 elif datatype == 3:
943 944 datatype_str = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
944 945 elif datatype == 4:
945 946 datatype_str = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
946 947 elif datatype == 5:
947 948 datatype_str = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
948 949 else:
949 950 raise ValueError('Data type was not defined')
950 951
951 952 self.dtype = datatype_str
952 953 #self.ippSeconds = 2 * 1000 * self.radarControllerHeaderObj.ipp / self.c
953 954 self.fileSizeByHeader = self.processingHeaderObj.dataBlocksPerFile * self.processingHeaderObj.blockSize + \
954 955 self.firstHeaderSize + self.basicHeaderSize * \
955 956 (self.processingHeaderObj.dataBlocksPerFile - 1)
956 957 # self.dataOut.channelList = numpy.arange(self.systemHeaderObj.numChannels)
957 958 # self.dataOut.channelIndexList = numpy.arange(self.systemHeaderObj.numChannels)
958 959 self.getBlockDimension()
959 960
960 961 def verifyFile(self, filename):
961 962
962 963 flag = True
963 964
964 965 try:
965 966 fp = open(filename, 'rb')
966 967 except IOError:
967 968 log.error("File {} can't be opened".format(filename), self.name)
968 969 return False
969 970
970 971 if self.online and self.waitDataBlock(0):
971 972 pass
972 973
973 974 basicHeaderObj = BasicHeader(LOCALTIME)
974 975 systemHeaderObj = SystemHeader()
975 976 radarControllerHeaderObj = RadarControllerHeader()
976 977 processingHeaderObj = ProcessingHeader()
977 978
978 979 if not(basicHeaderObj.read(fp)):
979 980 flag = False
980 981 if not(systemHeaderObj.read(fp)):
981 982 flag = False
982 983 if not(radarControllerHeaderObj.read(fp)):
983 984 flag = False
984 985 if not(processingHeaderObj.read(fp)):
985 986 flag = False
986 987 if not self.online:
987 988 dt1 = basicHeaderObj.datatime
988 989 pos = self.fileSize-processingHeaderObj.blockSize-24
989 990 if pos<0:
990 991 flag = False
991 992 log.error('Invalid size for file: {}'.format(self.filename), self.name)
992 993 else:
993 994 fp.seek(pos)
994 995 if not(basicHeaderObj.read(fp)):
995 996 flag = False
996 997 dt2 = basicHeaderObj.datatime
997 998 if not self.isDateTimeInRange(dt1, self.startDate, self.endDate, self.startTime, self.endTime) and not \
998 999 self.isDateTimeInRange(dt2, self.startDate, self.endDate, self.startTime, self.endTime):
999 1000 flag = False
1000 1001
1001 1002 fp.close()
1002 1003 return flag
1003 1004
1004 1005 def findDatafiles(self, path, startDate=None, endDate=None, expLabel='', ext='.r', walk=True, include_path=False):
1005 1006
1006 1007 path_empty = True
1007 1008
1008 1009 dateList = []
1009 1010 pathList = []
1010 1011
1011 1012 multi_path = path.split(',')
1012 1013
1013 1014 if not walk:
1014 1015
1015 1016 for single_path in multi_path:
1016 1017
1017 1018 if not os.path.isdir(single_path):
1018 1019 continue
1019 1020
1020 1021 fileList = glob.glob1(single_path, "*" + ext)
1021 1022
1022 1023 if not fileList:
1023 1024 continue
1024 1025
1025 1026 path_empty = False
1026 1027
1027 1028 fileList.sort()
1028 1029
1029 1030 for thisFile in fileList:
1030 1031
1031 1032 if not os.path.isfile(os.path.join(single_path, thisFile)):
1032 1033 continue
1033 1034
1034 1035 if not isRadarFile(thisFile):
1035 1036 continue
1036 1037
1037 1038 if not isFileInDateRange(thisFile, startDate, endDate):
1038 1039 continue
1039 1040
1040 1041 thisDate = getDateFromRadarFile(thisFile)
1041 1042
1042 1043 if thisDate in dateList or single_path in pathList:
1043 1044 continue
1044 1045
1045 1046 dateList.append(thisDate)
1046 1047 pathList.append(single_path)
1047 1048
1048 1049 else:
1049 1050 for single_path in multi_path:
1050 1051
1051 1052 if not os.path.isdir(single_path):
1052 1053 continue
1053 1054
1054 1055 dirList = []
1055 1056
1056 1057 for thisPath in os.listdir(single_path):
1057 1058
1058 1059 if not os.path.isdir(os.path.join(single_path, thisPath)):
1059 1060 continue
1060 1061
1061 1062 if not isRadarFolder(thisPath):
1062 1063 continue
1063 1064
1064 1065 if not isFolderInDateRange(thisPath, startDate, endDate):
1065 1066 continue
1066 1067
1067 1068 dirList.append(thisPath)
1068 1069
1069 1070 if not dirList:
1070 1071 continue
1071 1072
1072 1073 dirList.sort()
1073 1074
1074 1075 for thisDir in dirList:
1075 1076
1076 1077 datapath = os.path.join(single_path, thisDir, expLabel)
1077 1078 fileList = glob.glob1(datapath, "*" + ext)
1078 1079
1079 1080 if not fileList:
1080 1081 continue
1081 1082
1082 1083 path_empty = False
1083 1084
1084 1085 thisDate = getDateFromRadarFolder(thisDir)
1085 1086
1086 1087 pathList.append(datapath)
1087 1088 dateList.append(thisDate)
1088 1089
1089 1090 dateList.sort()
1090 1091
1091 1092 if walk:
1092 1093 pattern_path = os.path.join(multi_path[0], "[dYYYYDDD]", expLabel)
1093 1094 else:
1094 1095 pattern_path = multi_path[0]
1095 1096
1096 1097 if path_empty:
1097 1098 raise schainpy.admin.SchainError("[Reading] No *%s files in %s for %s to %s" % (ext, pattern_path, startDate, endDate))
1098 1099 else:
1099 1100 if not dateList:
1100 1101 raise schainpy.admin.SchainError("[Reading] Date range selected invalid [%s - %s]: No *%s files in %s)" % (startDate, endDate, ext, path))
1101 1102
1102 1103 if include_path:
1103 1104 return dateList, pathList
1104 1105
1105 1106 return dateList
1106 1107
1107 1108 def setup(self, **kwargs):
1108 1109
1109 1110 self.set_kwargs(**kwargs)
1110 1111 if not self.ext.startswith('.'):
1111 1112 self.ext = '.{}'.format(self.ext)
1112 1113
1113 1114 if self.server is not None:
1114 1115 if 'tcp://' in self.server:
1115 1116 address = server
1116 1117 else:
1117 1118 address = 'ipc:///tmp/%s' % self.server
1118 1119 self.server = address
1119 1120 self.context = zmq.Context()
1120 1121 self.receiver = self.context.socket(zmq.PULL)
1121 1122 self.receiver.connect(self.server)
1122 1123 time.sleep(0.5)
1123 1124 print('[Starting] ReceiverData from {}'.format(self.server))
1124 1125 else:
1125 1126 self.server = None
1126 1127 if self.path == None:
1127 1128 raise ValueError("[Reading] The path is not valid")
1128 1129
1129 1130 if self.online:
1130 1131 log.log("[Reading] Searching files in online mode...", self.name)
1131 1132
1132 1133 for nTries in range(self.nTries):
1133 1134 fullpath = self.searchFilesOnLine(self.path, self.startDate,
1134 1135 self.endDate, self.expLabel, self.ext, self.walk,
1135 1136 self.filefmt, self.folderfmt)
1136 1137
1137 1138 try:
1138 1139 fullpath = next(fullpath)
1139 1140 except:
1140 1141 fullpath = None
1141 1142
1142 1143 if fullpath:
1143 1144 break
1144 1145
1145 1146 log.warning(
1146 1147 'Waiting {} sec for a valid file in {}: try {} ...'.format(
1147 1148 self.delay, self.path, nTries + 1),
1148 1149 self.name)
1149 1150 time.sleep(self.delay)
1150 1151
1151 1152 if not(fullpath):
1152 1153 raise schainpy.admin.SchainError(
1153 1154 'There isn\'t any valid file in {}'.format(self.path))
1154 1155
1155 1156 pathname, filename = os.path.split(fullpath)
1156 1157 self.year = int(filename[1:5])
1157 1158 self.doy = int(filename[5:8])
1158 1159 self.set = int(filename[8:11]) - 1
1159 1160 else:
1160 1161 log.log("Searching files in {}".format(self.path), self.name)
1161 1162 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
1162 1163 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
1163 1164
1164 1165 self.setNextFile()
1165 1166
1166 1167 return
1167 1168
1168 1169 def getBasicHeader(self):
1169 1170
1170 1171 self.dataOut.utctime = self.basicHeaderObj.utc + self.basicHeaderObj.miliSecond / \
1171 1172 1000. + self.profileIndex * self.radarControllerHeaderObj.ippSeconds
1172 1173
1173 1174 self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
1174 1175
1175 1176 self.dataOut.timeZone = self.basicHeaderObj.timeZone
1176 1177
1177 1178 self.dataOut.dstFlag = self.basicHeaderObj.dstFlag
1178 1179
1179 1180 self.dataOut.errorCount = self.basicHeaderObj.errorCount
1180 1181
1181 1182 self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
1182 1183
1183 1184 self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
1184 1185
1185 1186 def getFirstHeader(self):
1186 1187
1187 1188 raise NotImplementedError
1188 1189
1189 1190 def getData(self):
1190 1191
1191 1192 raise NotImplementedError
1192 1193
1193 1194 def hasNotDataInBuffer(self):
1194 1195
1195 1196 raise NotImplementedError
1196 1197
1197 1198 def readBlock(self):
1198 1199
1199 1200 raise NotImplementedError
1200 1201
1201 1202 def isEndProcess(self):
1202 1203
1203 1204 return self.flagNoMoreFiles
1204 1205
1205 1206 def printReadBlocks(self):
1206 1207
1207 1208 print("[Reading] Number of read blocks per file %04d" % self.nReadBlocks)
1208 1209
1209 1210 def printTotalBlocks(self):
1210 1211
1211 1212 print("[Reading] Number of read blocks %04d" % self.nTotalBlocks)
1212 1213
1213 1214 def run(self, **kwargs):
1214 1215 """
1215 1216
1216 1217 Arguments:
1217 1218 path :
1218 1219 startDate :
1219 1220 endDate :
1220 1221 startTime :
1221 1222 endTime :
1222 1223 set :
1223 1224 expLabel :
1224 1225 ext :
1225 1226 online :
1226 1227 delay :
1227 1228 walk :
1228 1229 getblock :
1229 1230 nTxs :
1230 1231 realtime :
1231 1232 blocksize :
1232 1233 blocktime :
1233 1234 skip :
1234 1235 cursor :
1235 1236 warnings :
1236 1237 server :
1237 1238 verbose :
1238 1239 format :
1239 1240 oneDDict :
1240 1241 twoDDict :
1241 1242 independentParam :
1242 1243 """
1243 1244
1244 1245 if not(self.isConfig):
1245 1246 self.setup(**kwargs)
1246 1247 self.isConfig = True
1247 1248 if self.server is None:
1248 1249 self.getData()
1249 1250 else:
1250 1251 self.getFromServer()
1251 1252
1252 1253
1253 1254 class JRODataWriter(Reader):
1254 1255
1255 1256 """
1256 1257 Esta clase permite escribir datos a archivos procesados (.r o ,pdata). La escritura
1257 1258 de los datos siempre se realiza por bloques.
1258 1259 """
1259 1260
1260 1261 setFile = None
1261 1262 profilesPerBlock = None
1262 1263 blocksPerFile = None
1263 1264 nWriteBlocks = 0
1264 1265 fileDate = None
1265 1266
1266 1267 def __init__(self, dataOut=None):
1267 1268 raise NotImplementedError
1268 1269
1269 1270 def hasAllDataInBuffer(self):
1270 1271 raise NotImplementedError
1271 1272
1272 1273 def setBlockDimension(self):
1273 1274 raise NotImplementedError
1274 1275
1275 1276 def writeBlock(self):
1276 1277 raise NotImplementedError
1277 1278
1278 1279 def putData(self):
1279 1280 raise NotImplementedError
1280 1281
1281 1282 def getDtypeWidth(self):
1282 1283
1283 1284 dtype_index = get_dtype_index(self.dtype)
1284 1285 dtype_width = get_dtype_width(dtype_index)
1285 1286
1286 1287 return dtype_width
1287 1288
1288 1289 def getProcessFlags(self):
1289 1290
1290 1291 processFlags = 0
1291 1292
1292 1293 dtype_index = get_dtype_index(self.dtype)
1293 1294 procflag_dtype = get_procflag_dtype(dtype_index)
1294 1295
1295 1296 processFlags += procflag_dtype
1296 1297
1297 1298 if self.dataOut.flagDecodeData:
1298 1299 processFlags += PROCFLAG.DECODE_DATA
1299 1300
1300 1301 if self.dataOut.flagDeflipData:
1301 1302 processFlags += PROCFLAG.DEFLIP_DATA
1302 1303
1303 1304 if self.dataOut.code is not None:
1304 1305 processFlags += PROCFLAG.DEFINE_PROCESS_CODE
1305 1306
1306 1307 if self.dataOut.nCohInt > 1:
1307 1308 processFlags += PROCFLAG.COHERENT_INTEGRATION
1308 1309
1309 1310 if self.dataOut.type == "Spectra":
1310 1311 if self.dataOut.nIncohInt > 1:
1311 1312 processFlags += PROCFLAG.INCOHERENT_INTEGRATION
1312 1313
1313 1314 if self.dataOut.data_dc is not None:
1314 1315 processFlags += PROCFLAG.SAVE_CHANNELS_DC
1315 1316
1316 1317 if self.dataOut.flagShiftFFT:
1317 1318 processFlags += PROCFLAG.SHIFT_FFT_DATA
1318 1319
1319 1320 return processFlags
1320 1321
1321 1322 def setBasicHeader(self):
1322 1323
1323 1324 self.basicHeaderObj.size = self.basicHeaderSize # bytes
1324 1325 self.basicHeaderObj.version = self.versionFile
1325 1326 self.basicHeaderObj.dataBlock = self.nTotalBlocks
1326 1327 utc = numpy.floor(self.dataOut.utctime)
1327 1328 milisecond = (self.dataOut.utctime - utc) * 1000.0
1328 1329 self.basicHeaderObj.utc = utc
1329 1330 self.basicHeaderObj.miliSecond = milisecond
1330 1331 self.basicHeaderObj.timeZone = self.dataOut.timeZone
1331 1332 self.basicHeaderObj.dstFlag = self.dataOut.dstFlag
1332 1333 self.basicHeaderObj.errorCount = self.dataOut.errorCount
1333 1334
1334 1335 def setFirstHeader(self):
1335 1336 """
1336 1337 Obtiene una copia del First Header
1337 1338
1338 1339 Affected:
1339 1340
1340 1341 self.basicHeaderObj
1341 1342 self.systemHeaderObj
1342 1343 self.radarControllerHeaderObj
1343 1344 self.processingHeaderObj self.
1344 1345
1345 1346 Return:
1346 1347 None
1347 1348 """
1348 1349
1349 1350 raise NotImplementedError
1350 1351
1351 1352 def __writeFirstHeader(self):
1352 1353 """
1353 1354 Escribe el primer header del file es decir el Basic header y el Long header (SystemHeader, RadarControllerHeader, ProcessingHeader)
1354 1355
1355 1356 Affected:
1356 1357 __dataType
1357 1358
1358 1359 Return:
1359 1360 None
1360 1361 """
1361 1362
1362 1363 # CALCULAR PARAMETROS
1363 1364
1364 1365 sizeLongHeader = self.systemHeaderObj.size + \
1365 1366 self.radarControllerHeaderObj.size + self.processingHeaderObj.size
1366 1367 self.basicHeaderObj.size = self.basicHeaderSize + sizeLongHeader
1367 1368
1368 1369 self.basicHeaderObj.write(self.fp)
1369 1370 self.systemHeaderObj.write(self.fp)
1370 1371 self.radarControllerHeaderObj.write(self.fp)
1371 1372 self.processingHeaderObj.write(self.fp)
1372 1373
1373 1374 def __setNewBlock(self):
1374 1375 """
1375 1376 Si es un nuevo file escribe el First Header caso contrario escribe solo el Basic Header
1376 1377
1377 1378 Return:
1378 1379 0 : si no pudo escribir nada
1379 1380 1 : Si escribio el Basic el First Header
1380 1381 """
1381 1382 if self.fp == None:
1382 1383 self.setNextFile()
1383 1384
1384 1385 if self.flagIsNewFile:
1385 1386 return 1
1386 1387
1387 1388 if self.blockIndex < self.processingHeaderObj.dataBlocksPerFile:
1388 1389 self.basicHeaderObj.write(self.fp)
1389 1390 return 1
1390 1391
1391 1392 if not(self.setNextFile()):
1392 1393 return 0
1393 1394
1394 1395 return 1
1395 1396
1396 1397 def writeNextBlock(self):
1397 1398 """
1398 1399 Selecciona el bloque siguiente de datos y los escribe en un file
1399 1400
1400 1401 Return:
1401 1402 0 : Si no hizo pudo escribir el bloque de datos
1402 1403 1 : Si no pudo escribir el bloque de datos
1403 1404 """
1404 1405 if not(self.__setNewBlock()):
1405 1406 return 0
1406 1407
1407 1408 self.writeBlock()
1408 1409
1409 1410 print("[Writing] Block No. %d/%d" % (self.blockIndex,
1410 1411 self.processingHeaderObj.dataBlocksPerFile))
1411 1412
1412 1413 return 1
1413 1414
1414 1415 def setNextFile(self):
1415 1416 """Determina el siguiente file que sera escrito
1416 1417
1417 1418 Affected:
1418 1419 self.filename
1419 1420 self.subfolder
1420 1421 self.fp
1421 1422 self.setFile
1422 1423 self.flagIsNewFile
1423 1424
1424 1425 Return:
1425 1426 0 : Si el archivo no puede ser escrito
1426 1427 1 : Si el archivo esta listo para ser escrito
1427 1428 """
1428 1429 ext = self.ext
1429 1430 path = self.path
1430 1431
1431 1432 if self.fp != None:
1432 1433 self.fp.close()
1433 1434
1434 1435
1435 1436 if not os.path.exists(path):
1436 1437 os.mkdir(path)
1437 1438
1438 1439 timeTuple = time.localtime(self.dataOut.utctime)
1439 1440 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year, timeTuple.tm_yday)
1440 1441
1441 1442 fullpath = os.path.join(path, subfolder)
1442 1443 setFile = self.setFile
1443 1444
1444 1445 if not(os.path.exists(fullpath)):
1445 1446 os.mkdir(fullpath)
1446 1447 setFile = -1 # inicializo mi contador de seteo
1447 1448 else:
1448 1449 filesList = os.listdir(fullpath)
1449 1450 if len(filesList) > 0:
1450 1451 filesList = sorted(filesList, key=str.lower)
1451 1452 filen = filesList[-1]
1452 1453 # el filename debera tener el siguiente formato
1453 1454 # 0 1234 567 89A BCDE (hex)
1454 1455 # x YYYY DDD SSS .ext
1455 1456 if isNumber(filen[8:11]):
1456 1457 # inicializo mi contador de seteo al seteo del ultimo file
1457 1458 setFile = int(filen[8:11])
1458 1459 else:
1459 1460 setFile = -1
1460 1461 else:
1461 1462 setFile = -1 # inicializo mi contador de seteo
1462 1463
1463 1464 setFile += 1
1464 1465
1465 1466 # If this is a new day it resets some values
1466 1467 if self.dataOut.datatime.date() > self.fileDate:
1467 1468 setFile = 0
1468 1469 self.nTotalBlocks = 0
1469 1470
1470 1471 filen = '{}{:04d}{:03d}{:03d}{}'.format(
1471 1472 self.optchar, timeTuple.tm_year, timeTuple.tm_yday, setFile, ext)
1472 1473
1473 1474 filename = os.path.join(path, subfolder, filen)
1474 1475
1475 1476 fp = open(filename, 'wb')
1476 1477
1477 1478 self.blockIndex = 0
1478 1479 self.filename = filename
1479 1480 self.subfolder = subfolder
1480 1481 self.fp = fp
1481 1482 self.setFile = setFile
1482 1483 self.flagIsNewFile = 1
1483 1484 self.fileDate = self.dataOut.datatime.date()
1484 1485 self.setFirstHeader()
1485 1486
1486 1487 print('[Writing] Opening file: %s' % self.filename)
1487 1488
1488 1489 self.__writeFirstHeader()
1489 1490
1490 1491 return 1
1491 1492
1492 1493 def setup(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=None, ext=None, datatype=4):
1493 1494 """
1494 1495 Setea el tipo de formato en la cual sera guardada la data y escribe el First Header
1495 1496
1496 1497 Inputs:
1497 1498 path : directory where data will be saved
1498 1499 profilesPerBlock : number of profiles per block
1499 1500 set : initial file set
1500 1501 datatype : An integer number that defines data type:
1501 1502 0 : int8 (1 byte)
1502 1503 1 : int16 (2 bytes)
1503 1504 2 : int32 (4 bytes)
1504 1505 3 : int64 (8 bytes)
1505 1506 4 : float32 (4 bytes)
1506 1507 5 : double64 (8 bytes)
1507 1508
1508 1509 Return:
1509 1510 0 : Si no realizo un buen seteo
1510 1511 1 : Si realizo un buen seteo
1511 1512 """
1512 1513
1513 1514 if ext == None:
1514 1515 ext = self.ext
1515 1516
1516 1517 self.ext = ext.lower()
1517 1518
1518 1519 self.path = path
1519 1520
1520 1521 if set is None:
1521 1522 self.setFile = -1
1522 1523 else:
1523 1524 self.setFile = set - 1
1524 1525
1525 1526 self.blocksPerFile = blocksPerFile
1526 1527 self.profilesPerBlock = profilesPerBlock
1527 1528 self.dataOut = dataOut
1528 1529 self.fileDate = self.dataOut.datatime.date()
1529 1530 self.dtype = self.dataOut.dtype
1530 1531
1531 1532 if datatype is not None:
1532 1533 self.dtype = get_numpy_dtype(datatype)
1533 1534
1534 1535 if not(self.setNextFile()):
1535 1536 print("[Writing] There isn't a next file")
1536 1537 return 0
1537 1538
1538 1539 self.setBlockDimension()
1539 1540
1540 1541 return 1
1541 1542
1542 1543 def run(self, dataOut, path, blocksPerFile=100, profilesPerBlock=64, set=None, ext=None, datatype=4, **kwargs):
1543 1544
1544 1545 if not(self.isConfig):
1545 1546
1546 1547 self.setup(dataOut, path, blocksPerFile, profilesPerBlock=profilesPerBlock,
1547 1548 set=set, ext=ext, datatype=datatype, **kwargs)
1548 1549 self.isConfig = True
1549 1550
1550 1551 self.dataOut = dataOut
1551 1552 self.putData()
1552 1553 return self.dataOut
1553 1554
1554 1555 @MPDecorator
1555 1556 class printInfo(Operation):
1556 1557
1557 1558 def __init__(self):
1558 1559
1559 1560 Operation.__init__(self)
1560 1561 self.__printInfo = True
1561 1562
1562 1563 def run(self, dataOut, headers = ['systemHeaderObj', 'radarControllerHeaderObj', 'processingHeaderObj']):
1563 1564 if self.__printInfo == False:
1564 1565 return
1565 1566
1566 1567 for header in headers:
1567 1568 if hasattr(dataOut, header):
1568 1569 obj = getattr(dataOut, header)
1569 1570 if hasattr(obj, 'printInfo'):
1570 1571 obj.printInfo()
1571 1572 else:
1572 1573 print(obj)
1573 1574 else:
1574 1575 log.warning('Header {} Not found in object'.format(header))
1575 1576
1576 1577 self.__printInfo = False
Show file before | Show file after | Hide comments
@@ -1,651 +1,652
1 1 import os
2 2 import time
3 3 import datetime
4 4
5 5 import numpy
6 6 import h5py
7 7
8 8 import schainpy.admin
9 9 from schainpy.model.data.jrodata import *
10 10 from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
11 11 from schainpy.model.io.jroIO_base import *
12 12 from schainpy.utils import log
13 13
14 14
15 15 class HDFReader(Reader, ProcessingUnit):
16 16 """Processing unit to read HDF5 format files
17 17
18 18 This unit reads HDF5 files created with `HDFWriter` operation contains
19 19 by default two groups Data and Metadata all variables would be saved as `dataOut`
20 20 attributes.
21 21 It is possible to read any HDF5 file by given the structure in the `description`
22 22 parameter, also you can add extra values to metadata with the parameter `extras`.
23 23
24 24 Parameters:
25 25 -----------
26 26 path : str
27 27 Path where files are located.
28 28 startDate : date
29 29 Start date of the files
30 30 endDate : list
31 31 End date of the files
32 32 startTime : time
33 33 Start time of the files
34 34 endTime : time
35 35 End time of the files
36 36 description : dict, optional
37 37 Dictionary with the description of the HDF5 file
38 38 extras : dict, optional
39 39 Dictionary with extra metadata to be be added to `dataOut`
40 40
41 41 Examples
42 42 --------
43 43
44 44 desc = {
45 45 'Data': {
46 46 'data_output': ['u', 'v', 'w'],
47 47 'utctime': 'timestamps',
48 48 } ,
49 49 'Metadata': {
50 50 'heightList': 'heights'
51 51 }
52 52 }
53 53
54 54 desc = {
55 55 'Data': {
56 56 'data_output': 'winds',
57 57 'utctime': 'timestamps'
58 58 },
59 59 'Metadata': {
60 60 'heightList': 'heights'
61 61 }
62 62 }
63 63
64 64 extras = {
65 65 'timeZone': 300
66 66 }
67 67
68 68 reader = project.addReadUnit(
69 69 name='HDFReader',
70 70 path='/path/to/files',
71 71 startDate='2019/01/01',
72 72 endDate='2019/01/31',
73 73 startTime='00:00:00',
74 74 endTime='23:59:59',
75 75 # description=json.dumps(desc),
76 76 # extras=json.dumps(extras),
77 77 )
78 78
79 79 """
80 80
81 81 __attrs__ = ['path', 'startDate', 'endDate', 'startTime', 'endTime', 'description', 'extras']
82 82
83 83 def __init__(self):
84 84 ProcessingUnit.__init__(self)
85 85 self.dataOut = Parameters()
86 86 self.ext = ".hdf5"
87 87 self.optchar = "D"
88 88 self.meta = {}
89 89 self.data = {}
90 90 self.open_file = h5py.File
91 91 self.open_mode = 'r'
92 92 self.description = {}
93 93 self.extras = {}
94 94 self.filefmt = "*%Y%j***"
95 95 self.folderfmt = "*%Y%j"
96 96 self.utcoffset = 0
97 97
98 98 def setup(self, **kwargs):
99 99
100 100 self.set_kwargs(**kwargs)
101 101 if not self.ext.startswith('.'):
102 102 self.ext = '.{}'.format(self.ext)
103 103
104 104 if self.online:
105 105 log.log("Searching files in online mode...", self.name)
106 106
107 107 for nTries in range(self.nTries):
108 108 fullpath = self.searchFilesOnLine(self.path, self.startDate,
109 109 self.endDate, self.expLabel, self.ext, self.walk,
110 110 self.filefmt, self.folderfmt)
111 111 pathname, filename = os.path.split(fullpath)
112 112 #print(pathname,filename)
113 113 try:
114 114 fullpath = next(fullpath)
115 115
116 116 except:
117 117 fullpath = None
118 118
119 119 if fullpath:
120 120 break
121 121
122 122 log.warning(
123 123 'Waiting {} sec for a valid file in {}: try {} ...'.format(
124 124 self.delay, self.path, nTries + 1),
125 125 self.name)
126 126 time.sleep(self.delay)
127 127
128 128 if not(fullpath):
129 129 raise schainpy.admin.SchainError(
130 130 'There isn\'t any valid file in {}'.format(self.path))
131 131
132 132 pathname, filename = os.path.split(fullpath)
133 133 self.year = int(filename[1:5])
134 134 self.doy = int(filename[5:8])
135 135 self.set = int(filename[8:11]) - 1
136 136 else:
137 137 log.log("Searching files in {}".format(self.path), self.name)
138 138 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
139 139 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
140 140
141 141 self.setNextFile()
142 142
143 143 return
144 144
145 145
146 146 def readFirstHeader(self):
147 147 '''Read metadata and data'''
148 148
149 149 self.__readMetadata()
150 150 self.__readData()
151 151 self.__setBlockList()
152 152
153 153 if 'type' in self.meta:
154 154 self.dataOut = eval(self.meta['type'])()
155 155
156 156 for attr in self.meta:
157 print("attr: ", attr)
157 #print("attr: ", attr)
158 158 setattr(self.dataOut, attr, self.meta[attr])
159 159
160 160
161 161 self.blockIndex = 0
162 162
163 163 return
164 164
165 165 def __setBlockList(self):
166 166 '''
167 167 Selects the data within the times defined
168 168
169 169 self.fp
170 170 self.startTime
171 171 self.endTime
172 172 self.blockList
173 173 self.blocksPerFile
174 174
175 175 '''
176 176
177 177 startTime = self.startTime
178 178 endTime = self.endTime
179 179 thisUtcTime = self.data['utctime'] + self.utcoffset
180 180 self.interval = numpy.min(thisUtcTime[1:] - thisUtcTime[:-1])
181 181 thisDatetime = datetime.datetime.utcfromtimestamp(thisUtcTime[0])
182 182 self.startFileDatetime = thisDatetime
183 183 thisDate = thisDatetime.date()
184 184 thisTime = thisDatetime.time()
185 185
186 186 startUtcTime = (datetime.datetime.combine(thisDate, startTime) - datetime.datetime(1970, 1, 1)).total_seconds()
187 187 endUtcTime = (datetime.datetime.combine(thisDate, endTime) - datetime.datetime(1970, 1, 1)).total_seconds()
188 188
189 189 ind = numpy.where(numpy.logical_and(thisUtcTime >= startUtcTime, thisUtcTime < endUtcTime))[0]
190 190
191 191 self.blockList = ind
192 192 self.blocksPerFile = len(ind)
193 193 self.blocksPerFile = len(thisUtcTime)
194 194 return
195 195
196 196 def __readMetadata(self):
197 197 '''
198 198 Reads Metadata
199 199 '''
200 200
201 201 meta = {}
202 202
203 203 if self.description:
204 204 for key, value in self.description['Metadata'].items():
205 205 meta[key] = self.fp[value][()]
206 206 else:
207 207 grp = self.fp['Metadata']
208 208 for name in grp:
209 209 meta[name] = grp[name][()]
210 210
211 211 if self.extras:
212 212 for key, value in self.extras.items():
213 213 meta[key] = value
214 214 self.meta = meta
215 215
216 216 return
217 217
218 218
219 219
220 220 def checkForRealPath(self, nextFile, nextDay):
221 221
222 222 # print("check FRP")
223 223 # dt = self.startFileDatetime + datetime.timedelta(1)
224 224 # filename = '{}.{}{}'.format(self.path, dt.strftime('%Y%m%d'), self.ext)
225 225 # fullfilename = os.path.join(self.path, filename)
226 226 # print("check Path ",fullfilename,filename)
227 227 # if os.path.exists(fullfilename):
228 228 # return fullfilename, filename
229 229 # return None, filename
230 230 return None,None
231 231
232 232 def __readData(self):
233 233
234 234 data = {}
235 235
236 236 if self.description:
237 237 for key, value in self.description['Data'].items():
238 238 if isinstance(value, str):
239 239 if isinstance(self.fp[value], h5py.Dataset):
240 240 data[key] = self.fp[value][()]
241 241 elif isinstance(self.fp[value], h5py.Group):
242 242 array = []
243 243 for ch in self.fp[value]:
244 244 array.append(self.fp[value][ch][()])
245 245 data[key] = numpy.array(array)
246 246 elif isinstance(value, list):
247 247 array = []
248 248 for ch in value:
249 249 array.append(self.fp[ch][()])
250 250 data[key] = numpy.array(array)
251 251 else:
252 252 grp = self.fp['Data']
253 253 for name in grp:
254 254 if isinstance(grp[name], h5py.Dataset):
255 255 array = grp[name][()]
256 256 elif isinstance(grp[name], h5py.Group):
257 257 array = []
258 258 for ch in grp[name]:
259 259 array.append(grp[name][ch][()])
260 260 array = numpy.array(array)
261 261 else:
262 262 log.warning('Unknown type: {}'.format(name))
263 263
264 264 if name in self.description:
265 265 key = self.description[name]
266 266 else:
267 267 key = name
268 268 data[key] = array
269 269
270 270 self.data = data
271 271 return
272 272
273 273 def getData(self):
274 274 if not self.isDateTimeInRange(self.startFileDatetime, self.startDate, self.endDate, self.startTime, self.endTime):
275 275 self.dataOut.flagNoData = True
276 self.dataOut.error = True
276 self.blockIndex = self.blocksPerFile
277 #self.dataOut.error = True TERMINA EL PROGRAMA, removido
277 278 return
278 279 for attr in self.data:
279 280 if self.data[attr].ndim == 1:
280 281 setattr(self.dataOut, attr, self.data[attr][self.blockIndex])
281 282 else:
282 283 setattr(self.dataOut, attr, self.data[attr][:, self.blockIndex])
283 284
284 285 self.dataOut.flagNoData = False
285 286 self.blockIndex += 1
286 287
287 288 log.log("Block No. {}/{} -> {}".format(
288 289 self.blockIndex,
289 290 self.blocksPerFile,
290 291 self.dataOut.datatime.ctime()), self.name)
291 292
292 293 return
293 294
294 295 def run(self, **kwargs):
295 296
296 297 if not(self.isConfig):
297 298 self.setup(**kwargs)
298 299 self.isConfig = True
299 300
300 301 if self.blockIndex == self.blocksPerFile:
301 302 self.setNextFile()
302 303
303 304 self.getData()
304 305
305 306 return
306 307
307 308 @MPDecorator
308 309 class HDFWriter(Operation):
309 310 """Operation to write HDF5 files.
310 311
311 312 The HDF5 file contains by default two groups Data and Metadata where
312 313 you can save any `dataOut` attribute specified by `dataList` and `metadataList`
313 314 parameters, data attributes are normaly time dependent where the metadata
314 315 are not.
315 316 It is possible to customize the structure of the HDF5 file with the
316 317 optional description parameter see the examples.
317 318
318 319 Parameters:
319 320 -----------
320 321 path : str
321 322 Path where files will be saved.
322 323 blocksPerFile : int
323 324 Number of blocks per file
324 325 metadataList : list
325 326 List of the dataOut attributes that will be saved as metadata
326 327 dataList : int
327 328 List of the dataOut attributes that will be saved as data
328 329 setType : bool
329 330 If True the name of the files corresponds to the timestamp of the data
330 331 description : dict, optional
331 332 Dictionary with the desired description of the HDF5 file
332 333
333 334 Examples
334 335 --------
335 336
336 337 desc = {
337 338 'data_output': {'winds': ['z', 'w', 'v']},
338 339 'utctime': 'timestamps',
339 340 'heightList': 'heights'
340 341 }
341 342 desc = {
342 343 'data_output': ['z', 'w', 'v'],
343 344 'utctime': 'timestamps',
344 345 'heightList': 'heights'
345 346 }
346 347 desc = {
347 348 'Data': {
348 349 'data_output': 'winds',
349 350 'utctime': 'timestamps'
350 351 },
351 352 'Metadata': {
352 353 'heightList': 'heights'
353 354 }
354 355 }
355 356
356 357 writer = proc_unit.addOperation(name='HDFWriter')
357 358 writer.addParameter(name='path', value='/path/to/file')
358 359 writer.addParameter(name='blocksPerFile', value='32')
359 360 writer.addParameter(name='metadataList', value='heightList,timeZone')
360 361 writer.addParameter(name='dataList',value='data_output,utctime')
361 362 # writer.addParameter(name='description',value=json.dumps(desc))
362 363
363 364 """
364 365
365 366 ext = ".hdf5"
366 367 optchar = "D"
367 368 filename = None
368 369 path = None
369 370 setFile = None
370 371 fp = None
371 372 firsttime = True
372 373 #Configurations
373 374 blocksPerFile = None
374 375 blockIndex = None
375 376 dataOut = None
376 377 #Data Arrays
377 378 dataList = None
378 379 metadataList = None
379 380 currentDay = None
380 381 lastTime = None
381 382
382 383 def __init__(self):
383 384
384 385 Operation.__init__(self)
385 386 return
386 387
387 388 def setup(self, path=None, blocksPerFile=10, metadataList=None, dataList=None, setType=None, description=None):
388 389 self.path = path
389 390 self.blocksPerFile = blocksPerFile
390 391 self.metadataList = metadataList
391 392 self.dataList = [s.strip() for s in dataList]
392 393 self.setType = setType
393 394 self.description = description
394 395
395 396 if self.metadataList is None:
396 397 self.metadataList = self.dataOut.metadata_list
397 398
398 399 tableList = []
399 400 dsList = []
400 401
401 402 for i in range(len(self.dataList)):
402 403 dsDict = {}
403 404 if hasattr(self.dataOut, self.dataList[i]):
404 405 dataAux = getattr(self.dataOut, self.dataList[i])
405 406 dsDict['variable'] = self.dataList[i]
406 407 else:
407 408 log.warning('Attribute {} not found in dataOut', self.name)
408 409 continue
409 410
410 411 if dataAux is None:
411 412 continue
412 413 elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)):
413 414 dsDict['nDim'] = 0
414 415 else:
415 416 dsDict['nDim'] = len(dataAux.shape)
416 417 dsDict['shape'] = dataAux.shape
417 418 dsDict['dsNumber'] = dataAux.shape[0]
418 419 dsDict['dtype'] = dataAux.dtype
419 420
420 421 dsList.append(dsDict)
421 422
422 423 self.dsList = dsList
423 424 self.currentDay = self.dataOut.datatime.date()
424 425
425 426 def timeFlag(self):
426 427 currentTime = self.dataOut.utctime
427 428 timeTuple = time.localtime(currentTime)
428 429 dataDay = timeTuple.tm_yday
429 430 #print("time UTC: ",currentTime, self.dataOut.datatime)
430 431 if self.lastTime is None:
431 432 self.lastTime = currentTime
432 433 self.currentDay = dataDay
433 434 return False
434 435
435 436 timeDiff = currentTime - self.lastTime
436 437
437 438 #Si el dia es diferente o si la diferencia entre un dato y otro supera la hora
438 439 if dataDay != self.currentDay:
439 440 self.currentDay = dataDay
440 441 return True
441 442 elif timeDiff > 3*60*60:
442 443 self.lastTime = currentTime
443 444 return True
444 445 else:
445 446 self.lastTime = currentTime
446 447 return False
447 448
448 449 def run(self, dataOut, path, blocksPerFile=10, metadataList=None,
449 450 dataList=[], setType=None, description={}):
450 451
451 452 self.dataOut = dataOut
452 453 if not(self.isConfig):
453 454 self.setup(path=path, blocksPerFile=blocksPerFile,
454 455 metadataList=metadataList, dataList=dataList,
455 456 setType=setType, description=description)
456 457
457 458 self.isConfig = True
458 459 self.setNextFile()
459 460
460 461 self.putData()
461 462 return
462 463
463 464 def setNextFile(self):
464 465
465 466 ext = self.ext
466 467 path = self.path
467 468 setFile = self.setFile
468 469
469 470 timeTuple = time.gmtime(self.dataOut.utctime)
470 471 #print("path: ",timeTuple)
471 472 subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday)
472 473 fullpath = os.path.join(path, subfolder)
473 474
474 475 if os.path.exists(fullpath):
475 476 filesList = os.listdir(fullpath)
476 477 filesList = [k for k in filesList if k.startswith(self.optchar)]
477 478 if len( filesList ) > 0:
478 479 filesList = sorted(filesList, key=str.lower)
479 480 filen = filesList[-1]
480 481 # el filename debera tener el siguiente formato
481 482 # 0 1234 567 89A BCDE (hex)
482 483 # x YYYY DDD SSS .ext
483 484 if isNumber(filen[8:11]):
484 485 setFile = int(filen[8:11]) #inicializo mi contador de seteo al seteo del ultimo file
485 486 else:
486 487 setFile = -1
487 488 else:
488 489 setFile = -1 #inicializo mi contador de seteo
489 490 else:
490 491 os.makedirs(fullpath)
491 492 setFile = -1 #inicializo mi contador de seteo
492 493
493 494 if self.setType is None:
494 495 setFile += 1
495 496 file = '%s%4.4d%3.3d%03d%s' % (self.optchar,
496 497 timeTuple.tm_year,
497 498 timeTuple.tm_yday,
498 499 setFile,
499 500 ext )
500 501 else:
501 502 setFile = timeTuple.tm_hour*60+timeTuple.tm_min
502 503 file = '%s%4.4d%3.3d%04d%s' % (self.optchar,
503 504 timeTuple.tm_year,
504 505 timeTuple.tm_yday,
505 506 setFile,
506 507 ext )
507 508
508 509 self.filename = os.path.join( path, subfolder, file )
509 510
510 511 #Setting HDF5 File
511 512 self.fp = h5py.File(self.filename, 'w')
512 513 #write metadata
513 514 self.writeMetadata(self.fp)
514 515 #Write data
515 516 self.writeData(self.fp)
516 517
517 518 def getLabel(self, name, x=None):
518 519
519 520 if x is None:
520 521 if 'Data' in self.description:
521 522 data = self.description['Data']
522 523 if 'Metadata' in self.description:
523 524 data.update(self.description['Metadata'])
524 525 else:
525 526 data = self.description
526 527 if name in data:
527 528 if isinstance(data[name], str):
528 529 return data[name]
529 530 elif isinstance(data[name], list):
530 531 return None
531 532 elif isinstance(data[name], dict):
532 533 for key, value in data[name].items():
533 534 return key
534 535 return name
535 536 else:
536 537 if 'Metadata' in self.description:
537 538 meta = self.description['Metadata']
538 539 else:
539 540 meta = self.description
540 541 if name in meta:
541 542 if isinstance(meta[name], list):
542 543 return meta[name][x]
543 544 elif isinstance(meta[name], dict):
544 545 for key, value in meta[name].items():
545 546 return value[x]
546 547 if 'cspc' in name:
547 548 return 'pair{:02d}'.format(x)
548 549 else:
549 550 return 'channel{:02d}'.format(x)
550 551
551 552 def writeMetadata(self, fp):
552 553
553 554 if self.description:
554 555 if 'Metadata' in self.description:
555 556 grp = fp.create_group('Metadata')
556 557 else:
557 558 grp = fp
558 559 else:
559 560 grp = fp.create_group('Metadata')
560 561
561 562 for i in range(len(self.metadataList)):
562 563 if not hasattr(self.dataOut, self.metadataList[i]):
563 564 log.warning('Metadata: `{}` not found'.format(self.metadataList[i]), self.name)
564 565 continue
565 566 value = getattr(self.dataOut, self.metadataList[i])
566 567 if isinstance(value, bool):
567 568 if value is True:
568 569 value = 1
569 570 else:
570 571 value = 0
571 572 grp.create_dataset(self.getLabel(self.metadataList[i]), data=value)
572 573 return
573 574
574 575 def writeData(self, fp):
575 576
576 577 if self.description:
577 578 if 'Data' in self.description:
578 579 grp = fp.create_group('Data')
579 580 else:
580 581 grp = fp
581 582 else:
582 583 grp = fp.create_group('Data')
583 584
584 585 dtsets = []
585 586 data = []
586 587
587 588 for dsInfo in self.dsList:
588 589 if dsInfo['nDim'] == 0:
589 590 ds = grp.create_dataset(
590 591 self.getLabel(dsInfo['variable']),
591 592 (self.blocksPerFile, ),
592 593 chunks=True,
593 594 dtype=numpy.float64)
594 595 dtsets.append(ds)
595 596 data.append((dsInfo['variable'], -1))
596 597 else:
597 598 label = self.getLabel(dsInfo['variable'])
598 599 if label is not None:
599 600 sgrp = grp.create_group(label)
600 601 else:
601 602 sgrp = grp
602 603 for i in range(dsInfo['dsNumber']):
603 604 ds = sgrp.create_dataset(
604 605 self.getLabel(dsInfo['variable'], i),
605 606 (self.blocksPerFile, ) + dsInfo['shape'][1:],
606 607 chunks=True,
607 608 dtype=dsInfo['dtype'])
608 609 dtsets.append(ds)
609 610 data.append((dsInfo['variable'], i))
610 611 fp.flush()
611 612
612 613 log.log('Creating file: {}'.format(fp.filename), self.name)
613 614
614 615 self.ds = dtsets
615 616 self.data = data
616 617 self.firsttime = True
617 618 self.blockIndex = 0
618 619 return
619 620
620 621 def putData(self):
621 622
622 623 if (self.blockIndex == self.blocksPerFile) or self.timeFlag():
623 624 self.closeFile()
624 625 self.setNextFile()
625 626
626 627 for i, ds in enumerate(self.ds):
627 628 attr, ch = self.data[i]
628 629 if ch == -1:
629 630 ds[self.blockIndex] = getattr(self.dataOut, attr)
630 631 else:
631 632 ds[self.blockIndex] = getattr(self.dataOut, attr)[ch]
632 633
633 634 self.fp.flush()
634 635 self.blockIndex += 1
635 636 log.log('Block No. {}/{}'.format(self.blockIndex, self.blocksPerFile), self.name)
636 637
637 638 return
638 639
639 640 def closeFile(self):
640 641
641 642 if self.blockIndex != self.blocksPerFile:
642 643 for ds in self.ds:
643 644 ds.resize(self.blockIndex, axis=0)
644 645
645 646 if self.fp:
646 647 self.fp.flush()
647 648 self.fp.close()
648 649
649 650 def close(self):
650 651
651 652 self.closeFile()
Show file before | Show file after | Hide comments
@@ -1,1357 +1,1683
1 1 # Copyright (c) 2012-2020 Jicamarca Radio Observatory
2 2 # All rights reserved.
3 3 #
4 4 # Distributed under the terms of the BSD 3-clause license.
5 5 """Spectra processing Unit and operations
6 6
7 7 Here you will find the processing unit `SpectraProc` and several operations
8 8 to work with Spectra data type
9 9 """
10 10
11 11 import time
12 12 import itertools
13 13
14 14 import numpy
15 15 import math
16 16
17 17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
18 18 from schainpy.model.data.jrodata import Spectra
19 19 from schainpy.model.data.jrodata import hildebrand_sekhon
20 20 from schainpy.utils import log
21 21
22 22 from scipy.optimize import curve_fit
23 23
24 24
25 25 class SpectraProc(ProcessingUnit):
26 26
27 27 def __init__(self):
28 28
29 29 ProcessingUnit.__init__(self)
30 30
31 31 self.buffer = None
32 32 self.firstdatatime = None
33 33 self.profIndex = 0
34 34 self.dataOut = Spectra()
35 35 self.id_min = None
36 36 self.id_max = None
37 37 self.setupReq = False #Agregar a todas las unidades de proc
38 38
39 39 def __updateSpecFromVoltage(self):
40 40
41 41 self.dataOut.timeZone = self.dataIn.timeZone
42 42 self.dataOut.dstFlag = self.dataIn.dstFlag
43 43 self.dataOut.errorCount = self.dataIn.errorCount
44 44 self.dataOut.useLocalTime = self.dataIn.useLocalTime
45 45 try:
46 46 self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
47 47 except:
48 48 pass
49 49 self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
50 50 self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
51 51 self.dataOut.channelList = self.dataIn.channelList
52 52 self.dataOut.heightList = self.dataIn.heightList
53 53 self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
54 54 self.dataOut.nProfiles = self.dataOut.nFFTPoints
55 55 self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
56 56 self.dataOut.utctime = self.firstdatatime
57 57 self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
58 58 self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
59 59 self.dataOut.flagShiftFFT = False
60 60 self.dataOut.nCohInt = self.dataIn.nCohInt
61 61 self.dataOut.nIncohInt = 1
62 62 self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
63 63 self.dataOut.frequency = self.dataIn.frequency
64 64 self.dataOut.realtime = self.dataIn.realtime
65 65 self.dataOut.azimuth = self.dataIn.azimuth
66 66 self.dataOut.zenith = self.dataIn.zenith
67 67 self.dataOut.codeList = self.dataIn.codeList
68 68 self.dataOut.azimuthList = self.dataIn.azimuthList
69 69 self.dataOut.elevationList = self.dataIn.elevationList
70 70
71 71 def __getFft(self):
72 72 """
73 73 Convierte valores de Voltaje a Spectra
74 74
75 75 Affected:
76 76 self.dataOut.data_spc
77 77 self.dataOut.data_cspc
78 78 self.dataOut.data_dc
79 79 self.dataOut.heightList
80 80 self.profIndex
81 81 self.buffer
82 82 self.dataOut.flagNoData
83 83 """
84 84 fft_volt = numpy.fft.fft(
85 85 self.buffer, n=self.dataOut.nFFTPoints, axis=1)
86 86 fft_volt = fft_volt.astype(numpy.dtype('complex'))
87 87 dc = fft_volt[:, 0, :]
88 88
89 89 # calculo de self-spectra
90 90 fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
91 91 spc = fft_volt * numpy.conjugate(fft_volt)
92 92 spc = spc.real
93 93
94 94 blocksize = 0
95 95 blocksize += dc.size
96 96 blocksize += spc.size
97 97
98 98 cspc = None
99 99 pairIndex = 0
100 100 if self.dataOut.pairsList != None:
101 101 # calculo de cross-spectra
102 102 cspc = numpy.zeros(
103 103 (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
104 104 for pair in self.dataOut.pairsList:
105 105 if pair[0] not in self.dataOut.channelList:
106 106 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
107 107 str(pair), str(self.dataOut.channelList)))
108 108 if pair[1] not in self.dataOut.channelList:
109 109 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
110 110 str(pair), str(self.dataOut.channelList)))
111 111
112 112 cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
113 113 numpy.conjugate(fft_volt[pair[1], :, :])
114 114 pairIndex += 1
115 115 blocksize += cspc.size
116 116
117 117 self.dataOut.data_spc = spc
118 118 self.dataOut.data_cspc = cspc
119 119 self.dataOut.data_dc = dc
120 120 self.dataOut.blockSize = blocksize
121 121 self.dataOut.flagShiftFFT = False
122 122
123 123 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
124 124
125 125 if self.dataIn.type == "Spectra":
126 126 self.dataOut.copy(self.dataIn)
127 127 if shift_fft:
128 128 #desplaza a la derecha en el eje 2 determinadas posiciones
129 129 shift = int(self.dataOut.nFFTPoints/2)
130 130 self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
131 131
132 132 if self.dataOut.data_cspc is not None:
133 133 #desplaza a la derecha en el eje 2 determinadas posiciones
134 134 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
135 135 if pairsList:
136 136 self.__selectPairs(pairsList)
137 137
138 138 elif self.dataIn.type == "Voltage":
139 139
140 140 self.dataOut.flagNoData = True
141 141
142 142 if nFFTPoints == None:
143 143 raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
144 144
145 145 if nProfiles == None:
146 146 nProfiles = nFFTPoints
147 147
148 148 if ippFactor == None:
149 149 self.dataOut.ippFactor = 1
150 150
151 151 self.dataOut.nFFTPoints = nFFTPoints
152 152
153 153 if self.buffer is None:
154 154 self.buffer = numpy.zeros((self.dataIn.nChannels,
155 155 nProfiles,
156 156 self.dataIn.nHeights),
157 157 dtype='complex')
158 158
159 159 if self.dataIn.flagDataAsBlock:
160 160 nVoltProfiles = self.dataIn.data.shape[1]
161 161
162 162 if nVoltProfiles == nProfiles:
163 163 self.buffer = self.dataIn.data.copy()
164 164 self.profIndex = nVoltProfiles
165 165
166 166 elif nVoltProfiles < nProfiles:
167 167
168 168 if self.profIndex == 0:
169 169 self.id_min = 0
170 170 self.id_max = nVoltProfiles
171 171
172 172 self.buffer[:, self.id_min:self.id_max,
173 173 :] = self.dataIn.data
174 174 self.profIndex += nVoltProfiles
175 175 self.id_min += nVoltProfiles
176 176 self.id_max += nVoltProfiles
177 177 else:
178 178 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
179 179 self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
180 180 self.dataOut.flagNoData = True
181 181 else:
182 182 self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
183 183 self.profIndex += 1
184 184
185 185 if self.firstdatatime == None:
186 186 self.firstdatatime = self.dataIn.utctime
187 187
188 188 if self.profIndex == nProfiles:
189 189 self.__updateSpecFromVoltage()
190 190 if pairsList == None:
191 191 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
192 192 else:
193 193 self.dataOut.pairsList = pairsList
194 194 self.__getFft()
195 195 self.dataOut.flagNoData = False
196 196 self.firstdatatime = None
197 197 self.profIndex = 0
198 198 else:
199 199 raise ValueError("The type of input object '%s' is not valid".format(
200 200 self.dataIn.type))
201 201
202 202 def __selectPairs(self, pairsList):
203 203
204 204 if not pairsList:
205 205 return
206 206
207 207 pairs = []
208 208 pairsIndex = []
209 209
210 210 for pair in pairsList:
211 211 if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
212 212 continue
213 213 pairs.append(pair)
214 214 pairsIndex.append(pairs.index(pair))
215 215
216 216 self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
217 217 self.dataOut.pairsList = pairs
218 218
219 219 return
220 220
221 221 def selectFFTs(self, minFFT, maxFFT ):
222 222 """
223 223 Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
224 224 minFFT<= FFT <= maxFFT
225 225 """
226 226
227 227 if (minFFT > maxFFT):
228 228 raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
229 229
230 230 if (minFFT < self.dataOut.getFreqRange()[0]):
231 231 minFFT = self.dataOut.getFreqRange()[0]
232 232
233 233 if (maxFFT > self.dataOut.getFreqRange()[-1]):
234 234 maxFFT = self.dataOut.getFreqRange()[-1]
235 235
236 236 minIndex = 0
237 237 maxIndex = 0
238 238 FFTs = self.dataOut.getFreqRange()
239 239
240 240 inda = numpy.where(FFTs >= minFFT)
241 241 indb = numpy.where(FFTs <= maxFFT)
242 242
243 243 try:
244 244 minIndex = inda[0][0]
245 245 except:
246 246 minIndex = 0
247 247
248 248 try:
249 249 maxIndex = indb[0][-1]
250 250 except:
251 251 maxIndex = len(FFTs)
252 252
253 253 self.selectFFTsByIndex(minIndex, maxIndex)
254 254
255 255 return 1
256 256
257 257 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
258 258 newheis = numpy.where(
259 259 self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
260 260
261 261 if hei_ref != None:
262 262 newheis = numpy.where(self.dataOut.heightList > hei_ref)
263 263
264 264 minIndex = min(newheis[0])
265 265 maxIndex = max(newheis[0])
266 266 data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
267 267 heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
268 268
269 269 # determina indices
270 270 nheis = int(self.dataOut.radarControllerHeaderObj.txB /
271 271 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
272 272 avg_dB = 10 * \
273 273 numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
274 274 beacon_dB = numpy.sort(avg_dB)[-nheis:]
275 275 beacon_heiIndexList = []
276 276 for val in avg_dB.tolist():
277 277 if val >= beacon_dB[0]:
278 278 beacon_heiIndexList.append(avg_dB.tolist().index(val))
279 279
280 280 #data_spc = data_spc[:,:,beacon_heiIndexList]
281 281 data_cspc = None
282 282 if self.dataOut.data_cspc is not None:
283 283 data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
284 284 #data_cspc = data_cspc[:,:,beacon_heiIndexList]
285 285
286 286 data_dc = None
287 287 if self.dataOut.data_dc is not None:
288 288 data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
289 289 #data_dc = data_dc[:,beacon_heiIndexList]
290 290
291 291 self.dataOut.data_spc = data_spc
292 292 self.dataOut.data_cspc = data_cspc
293 293 self.dataOut.data_dc = data_dc
294 294 self.dataOut.heightList = heightList
295 295 self.dataOut.beacon_heiIndexList = beacon_heiIndexList
296 296
297 297 return 1
298 298
299 299 def selectFFTsByIndex(self, minIndex, maxIndex):
300 300 """
301 301
302 302 """
303 303
304 304 if (minIndex < 0) or (minIndex > maxIndex):
305 305 raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
306 306
307 307 if (maxIndex >= self.dataOut.nProfiles):
308 308 maxIndex = self.dataOut.nProfiles-1
309 309
310 310 #Spectra
311 311 data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
312 312
313 313 data_cspc = None
314 314 if self.dataOut.data_cspc is not None:
315 315 data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
316 316
317 317 data_dc = None
318 318 if self.dataOut.data_dc is not None:
319 319 data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
320 320
321 321 self.dataOut.data_spc = data_spc
322 322 self.dataOut.data_cspc = data_cspc
323 323 self.dataOut.data_dc = data_dc
324 324
325 325 self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
326 326 self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
327 327 self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
328 328
329 329 return 1
330 330
331 331 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
332 332 # validacion de rango
333 333 if minHei == None:
334 334 minHei = self.dataOut.heightList[0]
335 335
336 336 if maxHei == None:
337 337 maxHei = self.dataOut.heightList[-1]
338 338
339 339 if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
340 340 print('minHei: %.2f is out of the heights range' % (minHei))
341 341 print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
342 342 minHei = self.dataOut.heightList[0]
343 343
344 344 if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
345 345 print('maxHei: %.2f is out of the heights range' % (maxHei))
346 346 print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
347 347 maxHei = self.dataOut.heightList[-1]
348 348
349 349 # validacion de velocidades
350 350 velrange = self.dataOut.getVelRange(1)
351 351
352 352 if minVel == None:
353 353 minVel = velrange[0]
354 354
355 355 if maxVel == None:
356 356 maxVel = velrange[-1]
357 357
358 358 if (minVel < velrange[0]) or (minVel > maxVel):
359 359 print('minVel: %.2f is out of the velocity range' % (minVel))
360 360 print('minVel is setting to %.2f' % (velrange[0]))
361 361 minVel = velrange[0]
362 362
363 363 if (maxVel > velrange[-1]) or (maxVel < minVel):
364 364 print('maxVel: %.2f is out of the velocity range' % (maxVel))
365 365 print('maxVel is setting to %.2f' % (velrange[-1]))
366 366 maxVel = velrange[-1]
367 367
368 368 # seleccion de indices para rango
369 369 minIndex = 0
370 370 maxIndex = 0
371 371 heights = self.dataOut.heightList
372 372
373 373 inda = numpy.where(heights >= minHei)
374 374 indb = numpy.where(heights <= maxHei)
375 375
376 376 try:
377 377 minIndex = inda[0][0]
378 378 except:
379 379 minIndex = 0
380 380
381 381 try:
382 382 maxIndex = indb[0][-1]
383 383 except:
384 384 maxIndex = len(heights)
385 385
386 386 if (minIndex < 0) or (minIndex > maxIndex):
387 387 raise ValueError("some value in (%d,%d) is not valid" % (
388 388 minIndex, maxIndex))
389 389
390 390 if (maxIndex >= self.dataOut.nHeights):
391 391 maxIndex = self.dataOut.nHeights - 1
392 392
393 393 # seleccion de indices para velocidades
394 394 indminvel = numpy.where(velrange >= minVel)
395 395 indmaxvel = numpy.where(velrange <= maxVel)
396 396 try:
397 397 minIndexVel = indminvel[0][0]
398 398 except:
399 399 minIndexVel = 0
400 400
401 401 try:
402 402 maxIndexVel = indmaxvel[0][-1]
403 403 except:
404 404 maxIndexVel = len(velrange)
405 405
406 406 # seleccion del espectro
407 407 data_spc = self.dataOut.data_spc[:,
408 408 minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
409 409 # estimacion de ruido
410 410 noise = numpy.zeros(self.dataOut.nChannels)
411 411
412 412 for channel in range(self.dataOut.nChannels):
413 413 daux = data_spc[channel, :, :]
414 414 sortdata = numpy.sort(daux, axis=None)
415 415 noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
416 416
417 417 self.dataOut.noise_estimation = noise.copy()
418 418
419 419 return 1
420 420
421 421 class removeDC(Operation):
422 422
423 423 def run(self, dataOut, mode=2):
424 424 self.dataOut = dataOut
425 425 jspectra = self.dataOut.data_spc
426 426 jcspectra = self.dataOut.data_cspc
427 427
428 428 num_chan = jspectra.shape[0]
429 429 num_hei = jspectra.shape[2]
430 430
431 431 if jcspectra is not None:
432 432 jcspectraExist = True
433 433 num_pairs = jcspectra.shape[0]
434 434 else:
435 435 jcspectraExist = False
436 436
437 437 freq_dc = int(jspectra.shape[1] / 2)
438 438 ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
439 439 ind_vel = ind_vel.astype(int)
440 440
441 441 if ind_vel[0] < 0:
442 442 ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
443 443
444 444 if mode == 1:
445 445 jspectra[:, freq_dc, :] = (
446 446 jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2 # CORRECCION
447 447
448 448 if jcspectraExist:
449 449 jcspectra[:, freq_dc, :] = (
450 450 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
451 451
452 452 if mode == 2:
453 453
454 454 vel = numpy.array([-2, -1, 1, 2])
455 455 xx = numpy.zeros([4, 4])
456 456
457 457 for fil in range(4):
458 458 xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
459 459
460 460 xx_inv = numpy.linalg.inv(xx)
461 461 xx_aux = xx_inv[0, :]
462 462
463 463 for ich in range(num_chan):
464 464 yy = jspectra[ich, ind_vel, :]
465 465 jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
466 466
467 467 junkid = jspectra[ich, freq_dc, :] <= 0
468 468 cjunkid = sum(junkid)
469 469
470 470 if cjunkid.any():
471 471 jspectra[ich, freq_dc, junkid.nonzero()] = (
472 472 jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
473 473
474 474 if jcspectraExist:
475 475 for ip in range(num_pairs):
476 476 yy = jcspectra[ip, ind_vel, :]
477 477 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
478 478
479 479 self.dataOut.data_spc = jspectra
480 480 self.dataOut.data_cspc = jcspectra
481 481
482 482 return self.dataOut
483 483
484 484 # import matplotlib.pyplot as plt
485 485
486 486 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
487 487 z = (x - a1) / a2
488 488 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
489 489 return y
490 490
491 491
492 492 class CleanRayleigh(Operation):
493 493
494 494 def __init__(self):
495 495
496 496 Operation.__init__(self)
497 497 self.i=0
498 498 self.isConfig = False
499 499 self.__dataReady = False
500 500 self.__profIndex = 0
501 501 self.byTime = False
502 502 self.byProfiles = False
503 503
504 504 self.bloques = None
505 505 self.bloque0 = None
506 506
507 507 self.index = 0
508 508
509 509 self.buffer = 0
510 510 self.buffer2 = 0
511 511 self.buffer3 = 0
512 512
513 513
514 514 def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
515 515
516 516 self.nChannels = dataOut.nChannels
517 517 self.nProf = dataOut.nProfiles
518 518 self.nPairs = dataOut.data_cspc.shape[0]
519 519 self.pairsArray = numpy.array(dataOut.pairsList)
520 520 self.spectra = dataOut.data_spc
521 521 self.cspectra = dataOut.data_cspc
522 522 self.heights = dataOut.heightList #alturas totales
523 523 self.nHeights = len(self.heights)
524 524 self.min_hei = min_hei
525 525 self.max_hei = max_hei
526 526 if (self.min_hei == None):
527 527 self.min_hei = 0
528 528 if (self.max_hei == None):
529 529 self.max_hei = dataOut.heightList[-1]
530 530 self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
531 531 self.heightsClean = self.heights[self.hval] #alturas filtradas
532 532 self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
533 533 self.nHeightsClean = len(self.heightsClean)
534 534 self.channels = dataOut.channelList
535 535 self.nChan = len(self.channels)
536 536 self.nIncohInt = dataOut.nIncohInt
537 537 self.__initime = dataOut.utctime
538 538 self.maxAltInd = self.hval[-1]+1
539 539 self.minAltInd = self.hval[0]
540 540
541 541 self.crosspairs = dataOut.pairsList
542 542 self.nPairs = len(self.crosspairs)
543 543 self.normFactor = dataOut.normFactor
544 544 self.nFFTPoints = dataOut.nFFTPoints
545 545 self.ippSeconds = dataOut.ippSeconds
546 546 self.currentTime = self.__initime
547 547 self.pairsArray = numpy.array(dataOut.pairsList)
548 548 self.factor_stdv = factor_stdv
549 549 #print("CHANNELS: ",[x for x in self.channels])
550 550
551 551 if n != None :
552 552 self.byProfiles = True
553 553 self.nIntProfiles = n
554 554 else:
555 555 self.__integrationtime = timeInterval
556 556
557 557 self.__dataReady = False
558 558 self.isConfig = True
559 559
560 560
561 561
562 562 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
563 563 #print (dataOut.utctime)
564 564 if not self.isConfig :
565 565 #print("Setting config")
566 566 self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
567 567 #print("Config Done")
568 568 tini=dataOut.utctime
569 569
570 570 if self.byProfiles:
571 571 if self.__profIndex == self.nIntProfiles:
572 572 self.__dataReady = True
573 573 else:
574 574 if (tini - self.__initime) >= self.__integrationtime:
575 575 #print(tini - self.__initime,self.__profIndex)
576 576 self.__dataReady = True
577 577 self.__initime = tini
578 578
579 579 #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
580 580
581 581 if self.__dataReady:
582 582 #print("Data ready",self.__profIndex)
583 583 self.__profIndex = 0
584 584 jspc = self.buffer
585 585 jcspc = self.buffer2
586 586 #jnoise = self.buffer3
587 587 self.buffer = dataOut.data_spc
588 588 self.buffer2 = dataOut.data_cspc
589 589 #self.buffer3 = dataOut.noise
590 590 self.currentTime = dataOut.utctime
591 591 if numpy.any(jspc) :
592 592 #print( jspc.shape, jcspc.shape)
593 593 jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
594 594 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
595 595 self.__dataReady = False
596 596 #print( jspc.shape, jcspc.shape)
597 597 dataOut.flagNoData = False
598 598 else:
599 599 dataOut.flagNoData = True
600 600 self.__dataReady = False
601 601 return dataOut
602 602 else:
603 603 #print( len(self.buffer))
604 604 if numpy.any(self.buffer):
605 605 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
606 606 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
607 607 self.buffer3 += dataOut.data_dc
608 608 else:
609 609 self.buffer = dataOut.data_spc
610 610 self.buffer2 = dataOut.data_cspc
611 611 self.buffer3 = dataOut.data_dc
612 612 #print self.index, self.fint
613 613 #print self.buffer2.shape
614 614 dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
615 615 self.__profIndex += 1
616 616 return dataOut ## NOTE: REV
617 617
618 618
619 619 #index = tini.tm_hour*12+tini.tm_min/5
620 620 '''REVISAR'''
621 621 # jspc = jspc/self.nFFTPoints/self.normFactor
622 622 # jcspc = jcspc/self.nFFTPoints/self.normFactor
623 623
624 624
625 625
626 626 tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
627 627 dataOut.data_spc = tmp_spectra
628 628 dataOut.data_cspc = tmp_cspectra
629 629
630 630 #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
631 631
632 632 dataOut.data_dc = self.buffer3
633 633 dataOut.nIncohInt *= self.nIntProfiles
634 634 dataOut.utctime = self.currentTime #tiempo promediado
635 635 #print("Time: ",time.localtime(dataOut.utctime))
636 636 # dataOut.data_spc = sat_spectra
637 637 # dataOut.data_cspc = sat_cspectra
638 638 self.buffer = 0
639 639 self.buffer2 = 0
640 640 self.buffer3 = 0
641 641
642 642 return dataOut
643 643
644 644 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
645 645 #print("OP cleanRayleigh")
646 646 #import matplotlib.pyplot as plt
647 647 #for k in range(149):
648 648 #channelsProcssd = []
649 649 #channelA_ok = False
650 650 #rfunc = cspectra.copy() #self.bloques
651 651 rfunc = spectra.copy()
652 652 #rfunc = cspectra
653 653 #val_spc = spectra*0.0 #self.bloque0*0.0
654 654 #val_cspc = cspectra*0.0 #self.bloques*0.0
655 655 #in_sat_spectra = spectra.copy() #self.bloque0
656 656 #in_sat_cspectra = cspectra.copy() #self.bloques
657 657
658 658
659 659 ###ONLY FOR TEST:
660 660 raxs = math.ceil(math.sqrt(self.nPairs))
661 661 caxs = math.ceil(self.nPairs/raxs)
662 662 if self.nPairs <4:
663 663 raxs = 2
664 664 caxs = 2
665 665 #print(raxs, caxs)
666 666 fft_rev = 14 #nFFT to plot
667 667 hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot
668 668 hei_rev = hei_rev[0]
669 669 #print(hei_rev)
670 670
671 671 #print numpy.absolute(rfunc[:,0,0,14])
672 672
673 673 gauss_fit, covariance = None, None
674 674 for ih in range(self.minAltInd,self.maxAltInd):
675 675 for ifreq in range(self.nFFTPoints):
676 676 '''
677 677 ###ONLY FOR TEST:
678 678 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
679 679 fig, axs = plt.subplots(raxs, caxs)
680 680 fig2, axs2 = plt.subplots(raxs, caxs)
681 681 col_ax = 0
682 682 row_ax = 0
683 683 '''
684 684 #print(self.nPairs)
685 685 for ii in range(self.nChan): #PARES DE CANALES SELF y CROSS
686 686 # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
687 687 # continue
688 688 # if not self.crosspairs[ii][0] in channelsProcssd:
689 689 # channelA_ok = True
690 690 #print("pair: ",self.crosspairs[ii])
691 691 '''
692 692 ###ONLY FOR TEST:
693 693 if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1):
694 694 col_ax = 0
695 695 row_ax += 1
696 696 '''
697 697 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
698 698 #print(func2clean.shape)
699 699 val = (numpy.isfinite(func2clean)==True).nonzero()
700 700
701 701 if len(val)>0: #limitador
702 702 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
703 703 if min_val <= -40 :
704 704 min_val = -40
705 705 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
706 706 if max_val >= 200 :
707 707 max_val = 200
708 708 #print min_val, max_val
709 709 step = 1
710 710 #print("Getting bins and the histogram")
711 711 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
712 712 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
713 713 #print(len(y_dist),len(binstep[:-1]))
714 714 #print(row_ax,col_ax, " ..")
715 715 #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
716 716 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
717 717 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
718 718 parg = [numpy.amax(y_dist),mean,sigma]
719 719
720 720 newY = None
721 721
722 722 try :
723 723 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
724 724 mode = gauss_fit[1]
725 725 stdv = gauss_fit[2]
726 726 #print(" FIT OK",gauss_fit)
727 727 '''
728 728 ###ONLY FOR TEST:
729 729 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
730 730 newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
731 731 axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
732 732 axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
733 733 axs[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
734 734 '''
735 735 except:
736 736 mode = mean
737 737 stdv = sigma
738 738 #print("FIT FAIL")
739 739 #continue
740 740
741 741
742 742 #print(mode,stdv)
743 743 #Removing echoes greater than mode + std_factor*stdv
744 744 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
745 745 #noval tiene los indices que se van a remover
746 746 #print("Chan ",ii," novals: ",len(noval[0]))
747 747 if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
748 748 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
749 749 #print(novall)
750 750 #print(" ",self.pairsArray[ii])
751 751 #cross_pairs = self.pairsArray[ii]
752 752 #Getting coherent echoes which are removed.
753 753 # if len(novall[0]) > 0:
754 754 #
755 755 # val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
756 756 # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
757 757 # val_cspc[novall[0],ii,ifreq,ih] = 1
758 758 #print("OUT NOVALL 1")
759 759 try:
760 760 pair = (self.channels[ii],self.channels[ii + 1])
761 761 except:
762 762 pair = (99,99)
763 763 #print("par ", pair)
764 764 if ( pair in self.crosspairs):
765 765 q = self.crosspairs.index(pair)
766 766 #print("está aqui: ", q, (ii,ii + 1))
767 767 new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0])
768 768 cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
769 769
770 770 #if channelA_ok:
771 771 #chA = self.channels.index(cross_pairs[0])
772 772 new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0])
773 773 spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
774 774 #channelA_ok = False
775 775
776 776 # chB = self.channels.index(cross_pairs[1])
777 777 # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
778 778 # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
779 779 #
780 780 # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A
781 781 # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B
782 782 '''
783 783 ###ONLY FOR TEST:
784 784 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
785 785 func2clean = 10*numpy.log10(numpy.absolute(spectra[:,ii,ifreq,ih]))
786 786 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
787 787 axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
788 788 axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
789 789 axs2[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
790 790 '''
791 791 '''
792 792 ###ONLY FOR TEST:
793 793 col_ax += 1 #contador de ploteo columnas
794 794 ##print(col_ax)
795 795 ###ONLY FOR TEST:
796 796 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
797 797 title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
798 798 title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
799 799 fig.suptitle(title)
800 800 fig2.suptitle(title2)
801 801 plt.show()
802 802 '''
803 803 ##################################################################################################
804 804
805 805 #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
806 806 out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
807 807 out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
808 808 for ih in range(self.nHeights):
809 809 for ifreq in range(self.nFFTPoints):
810 810 for ich in range(self.nChan):
811 811 tmp = spectra[:,ich,ifreq,ih]
812 812 valid = (numpy.isfinite(tmp[:])==True).nonzero()
813 813
814 814 if len(valid[0]) >0 :
815 815 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
816 816
817 817 for icr in range(self.nPairs):
818 818 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
819 819 valid = (numpy.isfinite(tmp)==True).nonzero()
820 820 if len(valid[0]) > 0:
821 821 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
822 822
823 823 return out_spectra, out_cspectra
824 824
825 825 def REM_ISOLATED_POINTS(self,array,rth):
826 826 # import matplotlib.pyplot as plt
827 827 if rth == None :
828 828 rth = 4
829 829 #print("REM ISO")
830 830 num_prof = len(array[0,:,0])
831 831 num_hei = len(array[0,0,:])
832 832 n2d = len(array[:,0,0])
833 833
834 834 for ii in range(n2d) :
835 835 #print ii,n2d
836 836 tmp = array[ii,:,:]
837 837 #print tmp.shape, array[ii,101,:],array[ii,102,:]
838 838
839 839 # fig = plt.figure(figsize=(6,5))
840 840 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
841 841 # ax = fig.add_axes([left, bottom, width, height])
842 842 # x = range(num_prof)
843 843 # y = range(num_hei)
844 844 # cp = ax.contour(y,x,tmp)
845 845 # ax.clabel(cp, inline=True,fontsize=10)
846 846 # plt.show()
847 847
848 848 #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
849 849 tmp = numpy.reshape(tmp,num_prof*num_hei)
850 850 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
851 851 indxs2 = (tmp > 0).nonzero()
852 852
853 853 indxs1 = (indxs1[0])
854 854 indxs2 = indxs2[0]
855 855 #indxs1 = numpy.array(indxs1[0])
856 856 #indxs2 = numpy.array(indxs2[0])
857 857 indxs = None
858 858 #print indxs1 , indxs2
859 859 for iv in range(len(indxs2)):
860 860 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
861 861 #print len(indxs2), indv
862 862 if len(indv[0]) > 0 :
863 863 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
864 864 # print indxs
865 865 indxs = indxs[1:]
866 866 #print(indxs, len(indxs))
867 867 if len(indxs) < 4 :
868 868 array[ii,:,:] = 0.
869 869 return
870 870
871 871 xpos = numpy.mod(indxs ,num_hei)
872 872 ypos = (indxs / num_hei)
873 873 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
874 874 #print sx
875 875 xpos = xpos[sx]
876 876 ypos = ypos[sx]
877 877
878 878 # *********************************** Cleaning isolated points **********************************
879 879 ic = 0
880 880 while True :
881 881 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
882 882 #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
883 883 #plt.plot(r)
884 884 #plt.show()
885 885 no_coh1 = (numpy.isfinite(r)==True).nonzero()
886 886 no_coh2 = (r <= rth).nonzero()
887 887 #print r, no_coh1, no_coh2
888 888 no_coh1 = numpy.array(no_coh1[0])
889 889 no_coh2 = numpy.array(no_coh2[0])
890 890 no_coh = None
891 891 #print valid1 , valid2
892 892 for iv in range(len(no_coh2)):
893 893 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
894 894 if len(indv[0]) > 0 :
895 895 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
896 896 no_coh = no_coh[1:]
897 897 #print len(no_coh), no_coh
898 898 if len(no_coh) < 4 :
899 899 #print xpos[ic], ypos[ic], ic
900 900 # plt.plot(r)
901 901 # plt.show()
902 902 xpos[ic] = numpy.nan
903 903 ypos[ic] = numpy.nan
904 904
905 905 ic = ic + 1
906 906 if (ic == len(indxs)) :
907 907 break
908 908 #print( xpos, ypos)
909 909
910 910 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
911 911 #print indxs[0]
912 912 if len(indxs[0]) < 4 :
913 913 array[ii,:,:] = 0.
914 914 return
915 915
916 916 xpos = xpos[indxs[0]]
917 917 ypos = ypos[indxs[0]]
918 918 for i in range(0,len(ypos)):
919 919 ypos[i]=int(ypos[i])
920 920 junk = tmp
921 921 tmp = junk*0.0
922 922
923 923 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
924 924 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
925 925
926 926 #print array.shape
927 927 #tmp = numpy.reshape(tmp,(num_prof,num_hei))
928 928 #print tmp.shape
929 929
930 930 # fig = plt.figure(figsize=(6,5))
931 931 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
932 932 # ax = fig.add_axes([left, bottom, width, height])
933 933 # x = range(num_prof)
934 934 # y = range(num_hei)
935 935 # cp = ax.contour(y,x,array[ii,:,:])
936 936 # ax.clabel(cp, inline=True,fontsize=10)
937 937 # plt.show()
938 938 return array
939 939
940
941 class IntegrationFaradaySpectra(Operation):
942
943 __profIndex = 0
944 __withOverapping = False
945
946 __byTime = False
947 __initime = None
948 __lastdatatime = None
949 __integrationtime = None
950
951 __buffer_spc = None
952 __buffer_cspc = None
953 __buffer_dc = None
954
955 __dataReady = False
956
957 __timeInterval = None
958
959 n = None
960
961 def __init__(self):
962
963 Operation.__init__(self)
964
965 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None):
966 """
967 Set the parameters of the integration class.
968
969 Inputs:
970
971 n : Number of coherent integrations
972 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
973 overlapping :
974
975 """
976
977 self.__initime = None
978 self.__lastdatatime = 0
979
980 self.__buffer_spc = []
981 self.__buffer_cspc = []
982 self.__buffer_dc = 0
983
984 self.__profIndex = 0
985 self.__dataReady = False
986 self.__byTime = False
987
988 #self.ByLags = dataOut.ByLags ###REDEFINIR
989 self.ByLags = False
990
991 if DPL != None:
992 self.DPL=DPL
993 else:
994 #self.DPL=dataOut.DPL ###REDEFINIR
995 self.DPL=0
996
997 if n is None and timeInterval is None:
998 raise ValueError("n or timeInterval should be specified ...")
999
1000 if n is not None:
1001 self.n = int(n)
1002 else:
1003
1004 self.__integrationtime = int(timeInterval)
1005 self.n = None
1006 self.__byTime = True
1007
1008 def putData(self, data_spc, data_cspc, data_dc):
1009 """
1010 Add a profile to the __buffer_spc and increase in one the __profileIndex
1011
1012 """
1013
1014 self.__buffer_spc.append(data_spc)
1015
1016 if data_cspc is None:
1017 self.__buffer_cspc = None
1018 else:
1019 self.__buffer_cspc.append(data_cspc)
1020
1021 if data_dc is None:
1022 self.__buffer_dc = None
1023 else:
1024 self.__buffer_dc += data_dc
1025
1026 self.__profIndex += 1
1027
1028 return
1029
1030 def hildebrand_sekhon_Integration(self,data,navg):
1031
1032 sortdata = numpy.sort(data, axis=None)
1033 sortID=data.argsort()
1034 lenOfData = len(sortdata)
1035 nums_min = lenOfData*0.75
1036 if nums_min <= 5:
1037 nums_min = 5
1038 sump = 0.
1039 sumq = 0.
1040 j = 0
1041 cont = 1
1042 while((cont == 1)and(j < lenOfData)):
1043 sump += sortdata[j]
1044 sumq += sortdata[j]**2
1045 if j > nums_min:
1046 rtest = float(j)/(j-1) + 1.0/navg
1047 if ((sumq*j) > (rtest*sump**2)):
1048 j = j - 1
1049 sump = sump - sortdata[j]
1050 sumq = sumq - sortdata[j]**2
1051 cont = 0
1052 j += 1
1053 #lnoise = sump / j
1054
1055 return j,sortID
1056
1057 def pushData(self):
1058 """
1059 Return the sum of the last profiles and the profiles used in the sum.
1060
1061 Affected:
1062
1063 self.__profileIndex
1064
1065 """
1066 bufferH=None
1067 buffer=None
1068 buffer1=None
1069 buffer_cspc=None
1070 self.__buffer_spc=numpy.array(self.__buffer_spc)
1071 self.__buffer_cspc=numpy.array(self.__buffer_cspc)
1072 freq_dc = int(self.__buffer_spc.shape[2] / 2)
1073 #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
1074 for k in range(7,self.nHeights):
1075 buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
1076 outliers_IDs_cspc=[]
1077 cspc_outliers_exist=False
1078 #print("AQUIII")
1079 for i in range(self.nChannels):#dataOut.nChannels):
1080
1081 buffer1=numpy.copy(self.__buffer_spc[:,i,:,k])
1082 indexes=[]
1083 #sortIDs=[]
1084 outliers_IDs=[]
1085
1086 for j in range(self.nProfiles):
1087 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
1088 # continue
1089 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
1090 # continue
1091 buffer=buffer1[:,j]
1092 index,sortID=self.hildebrand_sekhon_Integration(buffer,1)
1093
1094 indexes.append(index)
1095 #sortIDs.append(sortID)
1096 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
1097
1098 outliers_IDs=numpy.array(outliers_IDs)
1099 outliers_IDs=outliers_IDs.ravel()
1100 outliers_IDs=numpy.unique(outliers_IDs)
1101 outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
1102 indexes=numpy.array(indexes)
1103 indexmin=numpy.min(indexes)
1104
1105 if indexmin != buffer1.shape[0]:
1106 cspc_outliers_exist=True
1107 ###sortdata=numpy.sort(buffer1,axis=0)
1108 ###avg2=numpy.mean(sortdata[:indexmin,:],axis=0)
1109 lt=outliers_IDs
1110 avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
1111
1112 for p in list(outliers_IDs):
1113 buffer1[p,:]=avg
1114
1115 self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
1116 ###cspc IDs
1117 #indexmin_cspc+=indexmin_cspc
1118 outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
1119
1120 #if not breakFlag:
1121 outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
1122 if cspc_outliers_exist:
1123 #sortdata=numpy.sort(buffer_cspc,axis=0)
1124 #avg=numpy.mean(sortdata[:indexmin_cpsc,:],axis=0)
1125 lt=outliers_IDs_cspc
1126
1127 avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
1128 for p in list(outliers_IDs_cspc):
1129 buffer_cspc[p,:]=avg
1130
1131 self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
1132 #else:
1133 #break
1134
1135
1136
1137
1138 buffer=None
1139 bufferH=None
1140 buffer1=None
1141 buffer_cspc=None
1142
1143 #print("cpsc",self.__buffer_cspc[:,0,0,0,0])
1144 #print(self.__profIndex)
1145 #exit()
1146
1147 buffer=None
1148 #print(self.__buffer_spc[:,1,3,20,0])
1149 #print(self.__buffer_spc[:,1,5,37,0])
1150 data_spc = numpy.sum(self.__buffer_spc,axis=0)
1151 data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
1152
1153 #print(numpy.shape(data_spc))
1154 #data_spc[1,4,20,0]=numpy.nan
1155
1156 #data_cspc = self.__buffer_cspc
1157 data_dc = self.__buffer_dc
1158 n = self.__profIndex
1159
1160 self.__buffer_spc = []
1161 self.__buffer_cspc = []
1162 self.__buffer_dc = 0
1163 self.__profIndex = 0
1164
1165 return data_spc, data_cspc, data_dc, n
1166
1167 def byProfiles(self, *args):
1168
1169 self.__dataReady = False
1170 avgdata_spc = None
1171 avgdata_cspc = None
1172 avgdata_dc = None
1173
1174 self.putData(*args)
1175
1176 if self.__profIndex == self.n:
1177
1178 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1179 self.n = n
1180 self.__dataReady = True
1181
1182 return avgdata_spc, avgdata_cspc, avgdata_dc
1183
1184 def byTime(self, datatime, *args):
1185
1186 self.__dataReady = False
1187 avgdata_spc = None
1188 avgdata_cspc = None
1189 avgdata_dc = None
1190
1191 self.putData(*args)
1192
1193 if (datatime - self.__initime) >= self.__integrationtime:
1194 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1195 self.n = n
1196 self.__dataReady = True
1197
1198 return avgdata_spc, avgdata_cspc, avgdata_dc
1199
1200 def integrate(self, datatime, *args):
1201
1202 if self.__profIndex == 0:
1203 self.__initime = datatime
1204
1205 if self.__byTime:
1206 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1207 datatime, *args)
1208 else:
1209 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1210
1211 if not self.__dataReady:
1212 return None, None, None, None
1213
1214 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1215
1216 def run(self, dataOut, n=None, DPL = None,timeInterval=None, overlapping=False):
1217 if n == 1:
1218 return dataOut
1219
1220 dataOut.flagNoData = True
1221
1222 if not self.isConfig:
1223 self.setup(dataOut, n, timeInterval, overlapping,DPL )
1224 self.isConfig = True
1225
1226 if not self.ByLags:
1227 self.nProfiles=dataOut.nProfiles
1228 self.nChannels=dataOut.nChannels
1229 self.nHeights=dataOut.nHeights
1230 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1231 dataOut.data_spc,
1232 dataOut.data_cspc,
1233 dataOut.data_dc)
1234 else:
1235 self.nProfiles=dataOut.nProfiles
1236 self.nChannels=dataOut.nChannels
1237 self.nHeights=dataOut.nHeights
1238 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1239 dataOut.dataLag_spc,
1240 dataOut.dataLag_cspc,
1241 dataOut.dataLag_dc)
1242
1243 if self.__dataReady:
1244
1245 if not self.ByLags:
1246
1247 dataOut.data_spc = numpy.squeeze(avgdata_spc)
1248 dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
1249 dataOut.data_dc = avgdata_dc
1250 else:
1251 dataOut.dataLag_spc = avgdata_spc
1252 dataOut.dataLag_cspc = avgdata_cspc
1253 dataOut.dataLag_dc = avgdata_dc
1254
1255 dataOut.data_spc=dataOut.dataLag_spc[:,:,:,dataOut.LagPlot]
1256 dataOut.data_cspc=dataOut.dataLag_cspc[:,:,:,dataOut.LagPlot]
1257 dataOut.data_dc=dataOut.dataLag_dc[:,:,dataOut.LagPlot]
1258
1259
1260 dataOut.nIncohInt *= self.n
1261 dataOut.utctime = avgdatatime
1262 dataOut.flagNoData = False
1263
1264 return dataOut
1265
940 1266 class removeInterference(Operation):
941 1267
942 1268 def removeInterference2(self):
943 1269
944 1270 cspc = self.dataOut.data_cspc
945 1271 spc = self.dataOut.data_spc
946 1272 Heights = numpy.arange(cspc.shape[2])
947 1273 realCspc = numpy.abs(cspc)
948 1274
949 1275 for i in range(cspc.shape[0]):
950 1276 LinePower= numpy.sum(realCspc[i], axis=0)
951 1277 Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
952 1278 SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
953 1279 InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
954 1280 InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
955 1281 InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
956 1282
957 1283
958 1284 InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
959 1285 #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
960 1286 if len(InterferenceRange)<int(cspc.shape[1]*0.3):
961 1287 cspc[i,InterferenceRange,:] = numpy.NaN
962 1288
963 1289 self.dataOut.data_cspc = cspc
964 1290
965 1291 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
966 1292
967 1293 jspectra = self.dataOut.data_spc
968 1294 jcspectra = self.dataOut.data_cspc
969 1295 jnoise = self.dataOut.getNoise()
970 1296 num_incoh = self.dataOut.nIncohInt
971 1297
972 1298 num_channel = jspectra.shape[0]
973 1299 num_prof = jspectra.shape[1]
974 1300 num_hei = jspectra.shape[2]
975 1301
976 1302 # hei_interf
977 1303 if hei_interf is None:
978 1304 count_hei = int(num_hei / 2)
979 1305 hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
980 1306 hei_interf = numpy.asarray(hei_interf)[0]
981 1307 # nhei_interf
982 1308 if (nhei_interf == None):
983 1309 nhei_interf = 5
984 1310 if (nhei_interf < 1):
985 1311 nhei_interf = 1
986 1312 if (nhei_interf > count_hei):
987 1313 nhei_interf = count_hei
988 1314 if (offhei_interf == None):
989 1315 offhei_interf = 0
990 1316
991 1317 ind_hei = list(range(num_hei))
992 1318 # mask_prof = numpy.asarray(range(num_prof - 2)) + 1
993 1319 # mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
994 1320 mask_prof = numpy.asarray(list(range(num_prof)))
995 1321 num_mask_prof = mask_prof.size
996 1322 comp_mask_prof = [0, num_prof / 2]
997 1323
998 1324 # noise_exist: Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
999 1325 if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
1000 1326 jnoise = numpy.nan
1001 1327 noise_exist = jnoise[0] < numpy.Inf
1002 1328
1003 1329 # Subrutina de Remocion de la Interferencia
1004 1330 for ich in range(num_channel):
1005 1331 # Se ordena los espectros segun su potencia (menor a mayor)
1006 1332 power = jspectra[ich, mask_prof, :]
1007 1333 power = power[:, hei_interf]
1008 1334 power = power.sum(axis=0)
1009 1335 psort = power.ravel().argsort()
1010 1336
1011 1337 # Se estima la interferencia promedio en los Espectros de Potencia empleando
1012 1338 junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
1013 1339 offhei_interf, nhei_interf + offhei_interf))]]]
1014 1340
1015 1341 if noise_exist:
1016 1342 # tmp_noise = jnoise[ich] / num_prof
1017 1343 tmp_noise = jnoise[ich]
1018 1344 junkspc_interf = junkspc_interf - tmp_noise
1019 1345 #junkspc_interf[:,comp_mask_prof] = 0
1020 1346
1021 1347 jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
1022 1348 jspc_interf = jspc_interf.transpose()
1023 1349 # Calculando el espectro de interferencia promedio
1024 1350 noiseid = numpy.where(
1025 1351 jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
1026 1352 noiseid = noiseid[0]
1027 1353 cnoiseid = noiseid.size
1028 1354 interfid = numpy.where(
1029 1355 jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
1030 1356 interfid = interfid[0]
1031 1357 cinterfid = interfid.size
1032 1358
1033 1359 if (cnoiseid > 0):
1034 1360 jspc_interf[noiseid] = 0
1035 1361
1036 1362 # Expandiendo los perfiles a limpiar
1037 1363 if (cinterfid > 0):
1038 1364 new_interfid = (
1039 1365 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
1040 1366 new_interfid = numpy.asarray(new_interfid)
1041 1367 new_interfid = {x for x in new_interfid}
1042 1368 new_interfid = numpy.array(list(new_interfid))
1043 1369 new_cinterfid = new_interfid.size
1044 1370 else:
1045 1371 new_cinterfid = 0
1046 1372
1047 1373 for ip in range(new_cinterfid):
1048 1374 ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
1049 1375 jspc_interf[new_interfid[ip]
1050 1376 ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
1051 1377
1052 1378 jspectra[ich, :, ind_hei] = jspectra[ich, :,
1053 1379 ind_hei] - jspc_interf # Corregir indices
1054 1380
1055 1381 # Removiendo la interferencia del punto de mayor interferencia
1056 1382 ListAux = jspc_interf[mask_prof].tolist()
1057 1383 maxid = ListAux.index(max(ListAux))
1058 1384
1059 1385 if cinterfid > 0:
1060 1386 for ip in range(cinterfid * (interf == 2) - 1):
1061 1387 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
1062 1388 (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
1063 1389 cind = len(ind)
1064 1390
1065 1391 if (cind > 0):
1066 1392 jspectra[ich, interfid[ip], ind] = tmp_noise * \
1067 1393 (1 + (numpy.random.uniform(cind) - 0.5) /
1068 1394 numpy.sqrt(num_incoh))
1069 1395
1070 1396 ind = numpy.array([-2, -1, 1, 2])
1071 1397 xx = numpy.zeros([4, 4])
1072 1398
1073 1399 for id1 in range(4):
1074 1400 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1075 1401
1076 1402 xx_inv = numpy.linalg.inv(xx)
1077 1403 xx = xx_inv[:, 0]
1078 1404 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1079 1405 yy = jspectra[ich, mask_prof[ind], :]
1080 1406 jspectra[ich, mask_prof[maxid], :] = numpy.dot(
1081 1407 yy.transpose(), xx)
1082 1408
1083 1409 indAux = (jspectra[ich, :, :] < tmp_noise *
1084 1410 (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
1085 1411 jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
1086 1412 (1 - 1 / numpy.sqrt(num_incoh))
1087 1413
1088 1414 # Remocion de Interferencia en el Cross Spectra
1089 1415 if jcspectra is None:
1090 1416 return jspectra, jcspectra
1091 1417 num_pairs = int(jcspectra.size / (num_prof * num_hei))
1092 1418 jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
1093 1419
1094 1420 for ip in range(num_pairs):
1095 1421
1096 1422 #-------------------------------------------
1097 1423
1098 1424 cspower = numpy.abs(jcspectra[ip, mask_prof, :])
1099 1425 cspower = cspower[:, hei_interf]
1100 1426 cspower = cspower.sum(axis=0)
1101 1427
1102 1428 cspsort = cspower.ravel().argsort()
1103 1429 junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
1104 1430 offhei_interf, nhei_interf + offhei_interf))]]]
1105 1431 junkcspc_interf = junkcspc_interf.transpose()
1106 1432 jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
1107 1433
1108 1434 ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
1109 1435
1110 1436 median_real = int(numpy.median(numpy.real(
1111 1437 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1112 1438 median_imag = int(numpy.median(numpy.imag(
1113 1439 junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
1114 1440 comp_mask_prof = [int(e) for e in comp_mask_prof]
1115 1441 junkcspc_interf[comp_mask_prof, :] = numpy.complex(
1116 1442 median_real, median_imag)
1117 1443
1118 1444 for iprof in range(num_prof):
1119 1445 ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
1120 1446 jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
1121 1447
1122 1448 # Removiendo la Interferencia
1123 1449 jcspectra[ip, :, ind_hei] = jcspectra[ip,
1124 1450 :, ind_hei] - jcspc_interf
1125 1451
1126 1452 ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
1127 1453 maxid = ListAux.index(max(ListAux))
1128 1454
1129 1455 ind = numpy.array([-2, -1, 1, 2])
1130 1456 xx = numpy.zeros([4, 4])
1131 1457
1132 1458 for id1 in range(4):
1133 1459 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
1134 1460
1135 1461 xx_inv = numpy.linalg.inv(xx)
1136 1462 xx = xx_inv[:, 0]
1137 1463
1138 1464 ind = (ind + maxid + num_mask_prof) % num_mask_prof
1139 1465 yy = jcspectra[ip, mask_prof[ind], :]
1140 1466 jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
1141 1467
1142 1468 # Guardar Resultados
1143 1469 self.dataOut.data_spc = jspectra
1144 1470 self.dataOut.data_cspc = jcspectra
1145 1471
1146 1472 return 1
1147 1473
1148 1474 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
1149 1475
1150 1476 self.dataOut = dataOut
1151 1477
1152 1478 if mode == 1:
1153 1479 self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
1154 1480 elif mode == 2:
1155 1481 self.removeInterference2()
1156 1482
1157 1483 return self.dataOut
1158 1484
1159 1485
1160 1486 class IncohInt(Operation):
1161 1487
1162 1488 __profIndex = 0
1163 1489 __withOverapping = False
1164 1490
1165 1491 __byTime = False
1166 1492 __initime = None
1167 1493 __lastdatatime = None
1168 1494 __integrationtime = None
1169 1495
1170 1496 __buffer_spc = None
1171 1497 __buffer_cspc = None
1172 1498 __buffer_dc = None
1173 1499
1174 1500 __dataReady = False
1175 1501
1176 1502 __timeInterval = None
1177 1503
1178 1504 n = None
1179 1505
1180 1506 def __init__(self):
1181 1507
1182 1508 Operation.__init__(self)
1183 1509
1184 1510 def setup(self, n=None, timeInterval=None, overlapping=False):
1185 1511 """
1186 1512 Set the parameters of the integration class.
1187 1513
1188 1514 Inputs:
1189 1515
1190 1516 n : Number of coherent integrations
1191 1517 timeInterval : Time of integration. If the parameter "n" is selected this one does not work
1192 1518 overlapping :
1193 1519
1194 1520 """
1195 1521
1196 1522 self.__initime = None
1197 1523 self.__lastdatatime = 0
1198 1524
1199 1525 self.__buffer_spc = 0
1200 1526 self.__buffer_cspc = 0
1201 1527 self.__buffer_dc = 0
1202 1528
1203 1529 self.__profIndex = 0
1204 1530 self.__dataReady = False
1205 1531 self.__byTime = False
1206 1532
1207 1533 if n is None and timeInterval is None:
1208 1534 raise ValueError("n or timeInterval should be specified ...")
1209 1535
1210 1536 if n is not None:
1211 1537 self.n = int(n)
1212 1538 else:
1213 1539
1214 1540 self.__integrationtime = int(timeInterval)
1215 1541 self.n = None
1216 1542 self.__byTime = True
1217 1543
1218 1544 def putData(self, data_spc, data_cspc, data_dc):
1219 1545 """
1220 1546 Add a profile to the __buffer_spc and increase in one the __profileIndex
1221 1547
1222 1548 """
1223 1549
1224 1550 self.__buffer_spc += data_spc
1225 1551
1226 1552 if data_cspc is None:
1227 1553 self.__buffer_cspc = None
1228 1554 else:
1229 1555 self.__buffer_cspc += data_cspc
1230 1556
1231 1557 if data_dc is None:
1232 1558 self.__buffer_dc = None
1233 1559 else:
1234 1560 self.__buffer_dc += data_dc
1235 1561
1236 1562 self.__profIndex += 1
1237 1563
1238 1564 return
1239 1565
1240 1566 def pushData(self):
1241 1567 """
1242 1568 Return the sum of the last profiles and the profiles used in the sum.
1243 1569
1244 1570 Affected:
1245 1571
1246 1572 self.__profileIndex
1247 1573
1248 1574 """
1249 1575
1250 1576 data_spc = self.__buffer_spc
1251 1577 data_cspc = self.__buffer_cspc
1252 1578 data_dc = self.__buffer_dc
1253 1579 n = self.__profIndex
1254 1580
1255 1581 self.__buffer_spc = 0
1256 1582 self.__buffer_cspc = 0
1257 1583 self.__buffer_dc = 0
1258 1584 self.__profIndex = 0
1259 1585
1260 1586 return data_spc, data_cspc, data_dc, n
1261 1587
1262 1588 def byProfiles(self, *args):
1263 1589
1264 1590 self.__dataReady = False
1265 1591 avgdata_spc = None
1266 1592 avgdata_cspc = None
1267 1593 avgdata_dc = None
1268 1594
1269 1595 self.putData(*args)
1270 1596
1271 1597 if self.__profIndex == self.n:
1272 1598
1273 1599 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1274 1600 self.n = n
1275 1601 self.__dataReady = True
1276 1602
1277 1603 return avgdata_spc, avgdata_cspc, avgdata_dc
1278 1604
1279 1605 def byTime(self, datatime, *args):
1280 1606
1281 1607 self.__dataReady = False
1282 1608 avgdata_spc = None
1283 1609 avgdata_cspc = None
1284 1610 avgdata_dc = None
1285 1611
1286 1612 self.putData(*args)
1287 1613
1288 1614 if (datatime - self.__initime) >= self.__integrationtime:
1289 1615 avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
1290 1616 self.n = n
1291 1617 self.__dataReady = True
1292 1618
1293 1619 return avgdata_spc, avgdata_cspc, avgdata_dc
1294 1620
1295 1621 def integrate(self, datatime, *args):
1296 1622
1297 1623 if self.__profIndex == 0:
1298 1624 self.__initime = datatime
1299 1625
1300 1626 if self.__byTime:
1301 1627 avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
1302 1628 datatime, *args)
1303 1629 else:
1304 1630 avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
1305 1631
1306 1632 if not self.__dataReady:
1307 1633 return None, None, None, None
1308 1634
1309 1635 return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
1310 1636
1311 1637 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
1312 1638 if n == 1:
1313 1639 return dataOut
1314 1640
1315 1641 dataOut.flagNoData = True
1316 1642
1317 1643 if not self.isConfig:
1318 1644 self.setup(n, timeInterval, overlapping)
1319 1645 self.isConfig = True
1320 1646
1321 1647 avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
1322 1648 dataOut.data_spc,
1323 1649 dataOut.data_cspc,
1324 1650 dataOut.data_dc)
1325 1651
1326 1652 if self.__dataReady:
1327 1653
1328 1654 dataOut.data_spc = avgdata_spc
1329 1655 dataOut.data_cspc = avgdata_cspc
1330 1656 dataOut.data_dc = avgdata_dc
1331 1657 dataOut.nIncohInt *= self.n
1332 1658 dataOut.utctime = avgdatatime
1333 1659 dataOut.flagNoData = False
1334 1660
1335 1661 return dataOut
1336 1662
1337 1663 class dopplerFlip(Operation):
1338 1664
1339 1665 def run(self, dataOut):
1340 1666 # arreglo 1: (num_chan, num_profiles, num_heights)
1341 1667 self.dataOut = dataOut
1342 1668 # JULIA-oblicua, indice 2
1343 1669 # arreglo 2: (num_profiles, num_heights)
1344 1670 jspectra = self.dataOut.data_spc[2]
1345 1671 jspectra_tmp = numpy.zeros(jspectra.shape)
1346 1672 num_profiles = jspectra.shape[0]
1347 1673 freq_dc = int(num_profiles / 2)
1348 1674 # Flip con for
1349 1675 for j in range(num_profiles):
1350 1676 jspectra_tmp[num_profiles-j-1]= jspectra[j]
1351 1677 # Intercambio perfil de DC con perfil inmediato anterior
1352 1678 jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
1353 1679 jspectra_tmp[freq_dc]= jspectra[freq_dc]
1354 1680 # canal modificado es re-escrito en el arreglo de canales
1355 1681 self.dataOut.data_spc[2] = jspectra_tmp
1356 1682
1357 1683 return self.dataOut
General Comments 0
You need to be logged in to leave comments. Login now