##// END OF EJS Templates
Adicionada operación Clean Rayleigh a Spectra Proc
joabAM -
r1385:83c0b1494611
parent child
Show More
@@ -9,7 +9,7
9 9 import numpy
10 10
11 11 from schainpy.model.graphics.jroplot_base import Plot, plt
12
12 import matplotlib.pyplot as plt
13 13
14 14 class SpectraHeisPlot(Plot):
15 15
@@ -33,17 +33,18 class SpectraHeisPlot(Plot):
33 33 meta = {}
34 34 spc = 10*numpy.log10(dataOut.data_spc / dataOut.normFactor)
35 35 data['spc_heis'] = spc
36
37 return data, meta
36
37 return data, meta
38 38
39 39 def plot(self):
40 40
41 41 c = 3E8
42 deltaHeight = self.data.yrange[1] - self.data.yrange[0]
42 deltaHeight = self.data.yrange[1] - self.data.yrange[0] # yrange = heightList
43 43 x = numpy.arange(-1*len(self.data.yrange)/2., len(self.data.yrange)/2.)*(c/(2*deltaHeight*len(self.data.yrange)*1000))
44 44 self.y = self.data[-1]['spc_heis']
45 45 self.titles = []
46 46
47 Maintitle = "Range from %d km to %d km" %(int(self.data.yrange[0]),int(self.data.yrange[-1]))
47 48 for n, ax in enumerate(self.axes):
48 49 ychannel = self.y[n,:]
49 50 if ax.firsttime:
@@ -54,7 +55,7 class SpectraHeisPlot(Plot):
54 55 ax.plt.set_data(x, ychannel)
55 56
56 57 self.titles.append("Channel {}: {:4.2f}dB".format(n, numpy.max(ychannel)))
57
58 plt.suptitle(Maintitle)
58 59
59 60 class RTIHeisPlot(Plot):
60 61
@@ -80,8 +81,8 class RTIHeisPlot(Plot):
80 81 spc = dataOut.data_spc / dataOut.normFactor
81 82 spc = 10*numpy.log10(numpy.average(spc, axis=1))
82 83 data['rti_heis'] = spc
83
84 return data, meta
84
85 return data, meta
85 86
86 87 def plot(self):
87 88
@@ -22,11 +22,11 class ScopePlot(Plot):
22 22 def setup(self):
23 23
24 24 self.xaxis = 'Range (Km)'
25 self.ncols = 1
26 self.nrows = 1
27 self.nplots = 1
25 self.nplots = len(self.data.channels)
26 self.nrows = int(numpy.ceil(self.nplots/2))
27 self.ncols = int(numpy.ceil(self.nplots/self.nrows))
28 28 self.ylabel = 'Intensity [dB]'
29 self.titles = ['Scope']
29 self.titles = ['Channel '+str(self.data.channels[i])+" " for i in self.data.channels]
30 30 self.colorbar = False
31 31 self.width = 6
32 32 self.height = 4
@@ -56,15 +56,13 class ScopePlot(Plot):
56 56
57 57 yreal = y[channelIndexList,:].real
58 58 yimag = y[channelIndexList,:].imag
59 title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
59 Maintitle = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
60 60 self.xlabel = "Range (Km)"
61 61 self.ylabel = "Intensity - IQ"
62 62
63 63 self.y = yreal
64 64 self.x = x
65 65
66 self.titles[0] = title
67
68 66 for i,ax in enumerate(self.axes):
69 67 title = "Channel %d" %(i)
70 68 if ax.firsttime:
@@ -75,6 +73,7 class ScopePlot(Plot):
75 73 else:
76 74 ax.plt_r.set_data(x, yreal[i,:])
77 75 ax.plt_i.set_data(x, yimag[i,:])
76 plt.suptitle(Maintitle)
78 77
79 78 def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle):
80 79 y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
@@ -99,6 +98,7 class ScopePlot(Plot):
99 98 else:
100 99 ax.plt_r.set_data(x, ychannel)
101 100
101
102 102 def plot_weatherpower(self, x, y, channelIndexList, thisDatetime, wintitle):
103 103
104 104
@@ -384,7 +384,7 def isRadarFolder(folder):
384 384
385 385
386 386 def isRadarFile(file):
387 try:
387 try:
388 388 year = int(file[1:5])
389 389 doy = int(file[5:8])
390 390 set = int(file[8:11])
@@ -395,10 +395,10 def isRadarFile(file):
395 395
396 396
397 397 def getDateFromRadarFile(file):
398 try:
398 try:
399 399 year = int(file[1:5])
400 400 doy = int(file[5:8])
401 set = int(file[8:11])
401 set = int(file[8:11])
402 402 except:
403 403 return None
404 404
@@ -417,7 +417,7 def getDateFromRadarFolder(folder):
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]]
@@ -484,7 +484,7 class Reader(object):
484 484
485 485 def run(self):
486 486
487 raise NotImplementedError
487 raise NotImplementedError
488 488
489 489 def getAllowedArgs(self):
490 490 if hasattr(self, '__attrs__'):
@@ -496,19 +496,19 class Reader(object):
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 folders = [x for f in path.split(',')
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 for folder in folders:
510 try:
511 dt = datetime.datetime.strptime(parse_format(folder, folderfmt), folderfmt).date()
509 for folder in folders:
510 try:
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:
@@ -517,38 +517,38 class Reader(object):
517 517 log.log('Skiping folder {}'.format(folder), self.name)
518 518 continue
519 519 return
520
521 def find_files(self, folders, ext, filefmt, startDate=None, endDate=None,
520
521 def find_files(self, folders, ext, filefmt, startDate=None, endDate=None,
522 522 expLabel='', last=False):
523
524 for path in folders:
523
524 for path in folders:
525 525 files = glob.glob1(path, '*{}'.format(ext))
526 526 files.sort()
527 527 if last:
528 if files:
528 if files:
529 529 fo = files[-1]
530 try:
530 try:
531 531 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
532 yield os.path.join(path, expLabel, fo)
533 except Exception as e:
532 yield os.path.join(path, expLabel, fo)
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 try:
541 dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date()
540 try:
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 continue
548 continue
549 549
550 550 def searchFilesOffLine(self, path, startDate, endDate,
551 expLabel, ext, walk,
551 expLabel, ext, walk,
552 552 filefmt, folderfmt):
553 553 """Search files in offline mode for the given arguments
554 554
@@ -561,12 +561,12 class Reader(object):
561 561 path, startDate, endDate, folderfmt)
562 562 else:
563 563 folders = path.split(',')
564
564
565 565 return self.find_files(
566 folders, ext, filefmt, startDate, endDate, expLabel)
566 folders, ext, filefmt, startDate, endDate, expLabel)
567 567
568 568 def searchFilesOnLine(self, path, startDate, endDate,
569 expLabel, ext, walk,
569 expLabel, ext, walk,
570 570 filefmt, folderfmt):
571 571 """Search for the last file of the last folder
572 572
@@ -579,13 +579,13 class Reader(object):
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
@@ -594,13 +594,13 class Reader(object):
594 594
595 595 while True:
596 596 if self.fp != None:
597 self.fp.close()
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')
@@ -609,10 +609,10 class Reader(object):
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()
@@ -625,7 +625,7 class Reader(object):
625 625 self.filename
626 626 self.fp
627 627 self.filesize
628
628
629 629 Return:
630 630 boolean
631 631
@@ -633,7 +633,7 class Reader(object):
633 633 nextFile = True
634 634 nextDay = False
635 635
636 for nFiles in range(self.nFiles+1):
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:
@@ -643,18 +643,18 class Reader(object):
643 643 self.name)
644 644 time.sleep(self.delay)
645 645 nextFile = False
646 continue
647
646 continue
647
648 648 if fullfilename is not None:
649 649 break
650
650
651 651 self.nTries = 1
652 nextFile = True
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 self.nTries = 3
657 self.nTries = 3
658 658
659 659 if fullfilename:
660 660 self.fileSize = os.path.getsize(fullfilename)
@@ -666,18 +666,18 class Reader(object):
666 666 self.flagNoMoreFiles = 0
667 667 self.fileIndex += 1
668 668 return 1
669 else:
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 return 0
680 return 0
681 681
682 682 self.filename = filename
683 683 self.fileSize = os.path.getsize(filename)
@@ -685,22 +685,22 class Reader(object):
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
693 if startDate <= dt.date() <= endDate:
694 if startTime <= dt.time() <= endTime:
695 return True
692 startDateTime= datetime.datetime.combine(startDate,startTime)
693 endDateTime = datetime.datetime.combine(endDate,endTime)
694 if startDateTime <= dt <= endDateTime:
695 return True
696 696 return False
697
697
698 698 def verifyFile(self, filename):
699 699 """Check for a valid file
700
700
701 701 Arguments:
702 702 filename -- full path filename
703
703
704 704 Return:
705 705 boolean
706 706 """
@@ -711,7 +711,7 class Reader(object):
711 711 """Check if the next file to be readed exists"""
712 712
713 713 raise NotImplementedError
714
714
715 715 def readFirstHeader(self):
716 716 """Parse the file header"""
717 717
@@ -783,8 +783,8 class JRODataReader(Reader):
783 783 Return:
784 784 str -- fullpath of the file
785 785 """
786
787
786
787
788 788 if nextFile:
789 789 self.set += 1
790 790 if nextDay:
@@ -796,7 +796,7 class JRODataReader(Reader):
796 796 prefixFileList = ['d', 'D']
797 797 elif self.ext.lower() == ".pdata": # spectra
798 798 prefixFileList = ['p', 'P']
799
799
800 800 # barrido por las combinaciones posibles
801 801 for prefixDir in prefixDirList:
802 802 thispath = self.path
@@ -816,9 +816,9 class JRODataReader(Reader):
816 816
817 817 if os.path.exists(fullfilename):
818 818 return fullfilename, filename
819
820 return None, filename
821
819
820 return None, filename
821
822 822 def __waitNewBlock(self):
823 823 """
824 824 Return 1 si se encontro un nuevo bloque de datos, 0 de otra forma.
@@ -860,9 +860,9 class JRODataReader(Reader):
860 860 def __setNewBlock(self):
861 861
862 862 if self.fp == None:
863 return 0
864
865 if self.flagIsNewFile:
863 return 0
864
865 if self.flagIsNewFile:
866 866 self.lastUTTime = self.basicHeaderObj.utc
867 867 return 1
868 868
@@ -875,12 +875,12 class JRODataReader(Reader):
875 875
876 876 currentSize = self.fileSize - self.fp.tell()
877 877 neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize
878
878
879 879 if (currentSize >= neededSize):
880 880 self.basicHeaderObj.read(self.fp)
881 881 self.lastUTTime = self.basicHeaderObj.utc
882 882 return 1
883
883
884 884 if self.__waitNewBlock():
885 885 self.lastUTTime = self.basicHeaderObj.utc
886 886 return 1
@@ -966,10 +966,10 class JRODataReader(Reader):
966 966 except IOError:
967 967 log.error("File {} can't be opened".format(filename), self.name)
968 968 return False
969
969
970 970 if self.online and self.waitDataBlock(0):
971 971 pass
972
972
973 973 basicHeaderObj = BasicHeader(LOCALTIME)
974 974 systemHeaderObj = SystemHeader()
975 975 radarControllerHeaderObj = RadarControllerHeader()
@@ -996,7 +996,7 class JRODataReader(Reader):
996 996 dt2 = basicHeaderObj.datatime
997 997 if not self.isDateTimeInRange(dt1, self.startDate, self.endDate, self.startTime, self.endTime) and not \
998 998 self.isDateTimeInRange(dt2, self.startDate, self.endDate, self.startTime, self.endTime):
999 flag = False
999 flag = False
1000 1000
1001 1001 fp.close()
1002 1002 return flag
@@ -1105,11 +1105,11 class JRODataReader(Reader):
1105 1105 return dateList
1106 1106
1107 1107 def setup(self, **kwargs):
1108
1108
1109 1109 self.set_kwargs(**kwargs)
1110 1110 if not self.ext.startswith('.'):
1111 1111 self.ext = '.{}'.format(self.ext)
1112
1112
1113 1113 if self.server is not None:
1114 1114 if 'tcp://' in self.server:
1115 1115 address = server
@@ -1131,36 +1131,36 class JRODataReader(Reader):
1131 1131
1132 1132 for nTries in range(self.nTries):
1133 1133 fullpath = self.searchFilesOnLine(self.path, self.startDate,
1134 self.endDate, self.expLabel, self.ext, self.walk,
1134 self.endDate, self.expLabel, self.ext, self.walk,
1135 1135 self.filefmt, self.folderfmt)
1136 1136
1137 1137 try:
1138 1138 fullpath = next(fullpath)
1139 1139 except:
1140 1140 fullpath = None
1141
1141
1142 1142 if fullpath:
1143 1143 break
1144 1144
1145 1145 log.warning(
1146 1146 'Waiting {} sec for a valid file in {}: try {} ...'.format(
1147 self.delay, self.path, nTries + 1),
1147 self.delay, self.path, nTries + 1),
1148 1148 self.name)
1149 1149 time.sleep(self.delay)
1150 1150
1151 1151 if not(fullpath):
1152 1152 raise schainpy.admin.SchainError(
1153 'There isn\'t any valid file in {}'.format(self.path))
1153 'There isn\'t any valid file in {}'.format(self.path))
1154 1154
1155 1155 pathname, filename = os.path.split(fullpath)
1156 1156 self.year = int(filename[1:5])
1157 1157 self.doy = int(filename[5:8])
1158 self.set = int(filename[8:11]) - 1
1158 self.set = int(filename[8:11]) - 1
1159 1159 else:
1160 1160 log.log("Searching files in {}".format(self.path), self.name)
1161 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
1161 self.filenameList = self.searchFilesOffLine(self.path, self.startDate,
1162 1162 self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt)
1163
1163
1164 1164 self.setNextFile()
1165 1165
1166 1166 return
@@ -1181,7 +1181,7 class JRODataReader(Reader):
1181 1181 self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
1182 1182
1183 1183 self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
1184
1184
1185 1185 def getFirstHeader(self):
1186 1186
1187 1187 raise NotImplementedError
@@ -1214,8 +1214,8 class JRODataReader(Reader):
1214 1214 """
1215 1215
1216 1216 Arguments:
1217 path :
1218 startDate :
1217 path :
1218 startDate :
1219 1219 endDate :
1220 1220 startTime :
1221 1221 endTime :
@@ -1284,7 +1284,7 class JRODataWriter(Reader):
1284 1284 dtype_width = get_dtype_width(dtype_index)
1285 1285
1286 1286 return dtype_width
1287
1287
1288 1288 def getProcessFlags(self):
1289 1289
1290 1290 processFlags = 0
@@ -1322,9 +1322,9 class JRODataWriter(Reader):
1322 1322
1323 1323 self.basicHeaderObj.size = self.basicHeaderSize # bytes
1324 1324 self.basicHeaderObj.version = self.versionFile
1325 self.basicHeaderObj.dataBlock = self.nTotalBlocks
1325 self.basicHeaderObj.dataBlock = self.nTotalBlocks
1326 1326 utc = numpy.floor(self.dataOut.utctime)
1327 milisecond = (self.dataOut.utctime - utc) * 1000.0
1327 milisecond = (self.dataOut.utctime - utc) * 1000.0
1328 1328 self.basicHeaderObj.utc = utc
1329 1329 self.basicHeaderObj.miliSecond = milisecond
1330 1330 self.basicHeaderObj.timeZone = self.dataOut.timeZone
@@ -1465,9 +1465,9 class JRODataWriter(Reader):
1465 1465 if self.dataOut.datatime.date() > self.fileDate:
1466 1466 setFile = 0
1467 1467 self.nTotalBlocks = 0
1468
1468
1469 1469 filen = '{}{:04d}{:03d}{:03d}{}'.format(
1470 self.optchar, timeTuple.tm_year, timeTuple.tm_yday, setFile, ext)
1470 self.optchar, timeTuple.tm_year, timeTuple.tm_yday, setFile, ext)
1471 1471
1472 1472 filename = os.path.join(path, subfolder, filen)
1473 1473
@@ -1515,11 +1515,11 class JRODataWriter(Reader):
1515 1515 self.ext = ext.lower()
1516 1516
1517 1517 self.path = path
1518
1518
1519 1519 if set is None:
1520 1520 self.setFile = -1
1521 1521 else:
1522 self.setFile = set - 1
1522 self.setFile = set - 1
1523 1523
1524 1524 self.blocksPerFile = blocksPerFile
1525 1525 self.profilesPerBlock = profilesPerBlock
@@ -347,7 +347,7 class AMISRReader(ProcessingUnit):
347 347
348 348 else:
349 349 #get the last file - 1
350 self.filenameList = [self.filenameList[-1]]
350 self.filenameList = [self.filenameList[-2]]
351 351 new_dirnameList = []
352 352 for dirname in self.dirnameList:
353 353 junk = numpy.array([dirname in x for x in self.filenameList])
@@ -421,7 +421,7 class AMISRReader(ProcessingUnit):
421 421 self.__selectDataForTimes(online=True)
422 422 filename = self.filenameList[0]
423 423 wait = 0
424 #self.__waitForNewFile=180 ## DEBUG:
424 self.__waitForNewFile=300 ## DEBUG:
425 425 while self.__filename_online == filename:
426 426 print('waiting %d seconds to get a new file...'%(self.__waitForNewFile))
427 427 if wait == 5:
@@ -641,13 +641,14 class AMISRReader(ProcessingUnit):
641 641 #self.dataOut.utctime = self.timeset[self.profileIndex] + t_comp
642 642 #print(self.dataOut.utctime)
643 643 self.dataOut.profileIndex = self.profileIndex
644 #print("N profile:",self.profileIndex,self.newProfiles,self.nblocks,self.dataOut.utctime)
644 645 self.dataOut.flagNoData = False
645 646 # if indexprof == 0:
646 647 # print self.dataOut.utctime
647 648
648 649 self.profileIndex += 1
649 650
650 return self.dataOut.data
651 #return self.dataOut.data
651 652
652 653
653 654 def run(self, **kwargs):
@@ -660,3 +661,4 class AMISRReader(ProcessingUnit):
660 661 self.isConfig = True
661 662
662 663 self.getData()
664 return(self.dataOut.data)
@@ -5,7 +5,6 from schainpy.model.data.jrodata import SpectraHeis
5 5 from schainpy.utils import log
6 6
7 7
8
9 8 class SpectraHeisProc(ProcessingUnit):
10 9
11 10 def __init__(self):#, **kwargs):
@@ -89,7 +88,7 class SpectraHeisProc(ProcessingUnit):
89 88 if self.dataIn.type == "Fits":
90 89 self.__updateObjFromFits()
91 90 self.dataOut.flagNoData = False
92 return
91 return
93 92
94 93 if self.dataIn.type == "SpectraHeis":
95 94 self.dataOut.copy(self.dataIn)
@@ -343,5 +342,5 class IncohInt4SpectraHeis(Operation):
343 342 # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nIncohInt
344 343 # dataOut.timeInterval = self.__timeInterval*self.n
345 344 dataOut.flagNoData = False
346
347 return dataOut No newline at end of file
345
346 return dataOut
This diff has been collapsed as it changes many lines, (559 lines changed) Show them Hide them
@@ -12,12 +12,15 import time
12 12 import itertools
13 13
14 14 import numpy
15 import math
15 16
16 17 from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation
17 18 from schainpy.model.data.jrodata import Spectra
18 19 from schainpy.model.data.jrodata import hildebrand_sekhon
19 20 from schainpy.utils import log
20 21
22 from scipy.optimize import curve_fit
23
21 24
22 25 class SpectraProc(ProcessingUnit):
23 26
@@ -478,6 +481,562 class removeDC(Operation):
478 481
479 482 return self.dataOut
480 483
484 # import matplotlib.pyplot as plt
485
486 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
487 z = (x - a1) / a2
488 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
489 return y
490 class CleanRayleigh(Operation):
491
492 def __init__(self):
493
494 Operation.__init__(self)
495 self.i=0
496 self.isConfig = False
497 self.__dataReady = False
498 self.__profIndex = 0
499 self.byTime = False
500 self.byProfiles = False
501
502 self.bloques = None
503 self.bloque0 = None
504
505 self.index = 0
506
507 self.buffer = 0
508 self.buffer2 = 0
509 self.buffer3 = 0
510 #self.min_hei = None
511 #self.max_hei = None
512
513 def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv):
514
515 self.nChannels = dataOut.nChannels
516 self.nProf = dataOut.nProfiles
517 self.nPairs = dataOut.data_cspc.shape[0]
518 self.pairsArray = numpy.array(dataOut.pairsList)
519 self.spectra = dataOut.data_spc
520 self.cspectra = dataOut.data_cspc
521 self.heights = dataOut.heightList #alturas totales
522 self.nHeights = len(self.heights)
523 self.min_hei = min_hei
524 self.max_hei = max_hei
525 if (self.min_hei == None):
526 self.min_hei = 0
527 if (self.max_hei == None):
528 self.max_hei = dataOut.heightList[-1]
529 self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero()
530 self.heightsClean = self.heights[self.hval] #alturas filtradas
531 self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas
532 self.nHeightsClean = len(self.heightsClean)
533 self.channels = dataOut.channelList
534 self.nChan = len(self.channels)
535 self.nIncohInt = dataOut.nIncohInt
536 self.__initime = dataOut.utctime
537 self.maxAltInd = self.hval[-1]+1
538 self.minAltInd = self.hval[0]
539
540 self.crosspairs = dataOut.pairsList
541 self.nPairs = len(self.crosspairs)
542 self.normFactor = dataOut.normFactor
543 self.nFFTPoints = dataOut.nFFTPoints
544 self.ippSeconds = dataOut.ippSeconds
545 self.currentTime = self.__initime
546 self.pairsArray = numpy.array(dataOut.pairsList)
547 self.factor_stdv = factor_stdv
548
549
550
551 if n != None :
552 self.byProfiles = True
553 self.nIntProfiles = n
554 else:
555 self.__integrationtime = timeInterval
556
557 self.__dataReady = False
558 self.isConfig = True
559
560
561
562 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
563 #print (dataOut.utctime)
564 if not self.isConfig :
565 #print("Setting config")
566 self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
567 #print("Config Done")
568 tini=dataOut.utctime
569
570 if self.byProfiles:
571 if self.__profIndex == self.nIntProfiles:
572 self.__dataReady = True
573 else:
574 if (tini - self.__initime) >= self.__integrationtime:
575 #print(tini - self.__initime,self.__profIndex)
576 self.__dataReady = True
577 self.__initime = tini
578
579 #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
580
581 if self.__dataReady:
582 print("Data ready",self.__profIndex)
583 self.__profIndex = 0
584 jspc = self.buffer
585 jcspc = self.buffer2
586 #jnoise = self.buffer3
587 self.buffer = dataOut.data_spc
588 self.buffer2 = dataOut.data_cspc
589 #self.buffer3 = dataOut.noise
590 self.currentTime = dataOut.utctime
591 if numpy.any(jspc) :
592 #print( jspc.shape, jcspc.shape)
593 jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
594 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
595 self.__dataReady = False
596 #print( jspc.shape, jcspc.shape)
597 dataOut.flagNoData = False
598 else:
599 dataOut.flagNoData = True
600 self.__dataReady = False
601 return dataOut
602 else:
603 #print( len(self.buffer))
604 if numpy.any(self.buffer):
605 self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
606 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
607 self.buffer3 += dataOut.data_dc
608 else:
609 self.buffer = dataOut.data_spc
610 self.buffer2 = dataOut.data_cspc
611 self.buffer3 = dataOut.data_dc
612 #print self.index, self.fint
613 #print self.buffer2.shape
614 dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO
615 self.__profIndex += 1
616 return dataOut ## NOTE: REV
617
618 # if self.index == 0 and self.fint == 1 :
619 # if jspc != None:
620 # print len(jspc), jspc.shape
621 # jspc= numpy.reshape(jspc,(4,128,63,len(jspc)/4))
622 # print jspc.shape
623 # dataOut.flagNoData = True
624 # return dataOut
625 # if path != None:
626 # sys.path.append(path)
627 # self.library = importlib.import_module(file)
628 #
629 # To be inserted as a parameter
630
631 #Set constants
632 #constants = self.library.setConstants(dataOut)
633 #dataOut.constants = constants
634
635 #snrth= 20
636
637 #crosspairs = dataOut.groupList
638 #noise = dataOut.noise
639 #print( nProf,heights)
640 #print( jspc.shape, jspc.shape[0])
641 #print noise
642 #print jnoise[len(jnoise)-1,:], numpy.nansum(jnoise,axis=0)/len(jnoise)
643 #jnoise = jnoise/N
644 #noise = numpy.nansum(jnoise,axis=0)#/len(jnoise)
645 #print( noise)
646 #power = numpy.sum(spectra, axis=1)
647 #print power[0,:]
648 #print("CROSSPAIRS",crosspairs)
649 #nPairs = len(crosspairs)
650 #print(numpy.shape(dataOut.data_spc))
651
652 #absc = dataOut.abscissaList[:-1]
653
654 #print absc.shape
655 #nBlocks=149
656 #print('spectra', spectra.shape)
657 #print('noise print', crosspairs)
658 #print('spectra', spectra.shape)
659 #print('cspectra', cspectra.shape)
660 #print numpy.array(dataOut.data_pre[1]).shape
661 #spec, cspec = self.__DiffCoherent(snrth, spectra, cspectra, nProf, heights,nChan, nHei, nPairs, channels, noise*nProf, crosspairs)
662
663
664 #index = tini.tm_hour*12+tini.tm_min/5
665 # jspc = jspc/self.nFFTPoints/self.normFactor
666 # jcspc = jcspc/self.nFFTPoints/self.normFactor
667
668
669
670 #dataOut.data_spc,dataOut.data_cspc = self.CleanRayleigh(dataOut,jspc,jcspc,crosspairs,heights,channels,nProf,nHei,nChan,nPairs,nIncohInt,nBlocks=nBlocks)
671 #tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.min_hei,self.max_hei)
672 tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
673 #jspectra = tmp_spectra*len(jspc[:,0,0,0])
674 #jcspectra = tmp_cspectra*len(jspc[:,0,0,0])
675
676 dataOut.data_spc = tmp_spectra
677 dataOut.data_cspc = tmp_cspectra
678 dataOut.data_dc = self.buffer3
679 dataOut.nIncohInt *= self.nIntProfiles
680 dataOut.utctime = self.currentTime #tiempo promediado
681 print("Time: ",time.localtime(dataOut.utctime))
682 # dataOut.data_spc = sat_spectra
683 # dataOut.data_cspc = sat_cspectra
684 self.buffer = 0
685 self.buffer2 = 0
686 self.buffer3 = 0
687
688 return dataOut
689
690 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
691 print("OP cleanRayleigh")
692 #import matplotlib.pyplot as plt
693 #for k in range(149):
694
695 rfunc = cspectra.copy() #self.bloques
696 val_spc = spectra*0.0 #self.bloque0*0.0
697 val_cspc = cspectra*0.0 #self.bloques*0.0
698 in_sat_spectra = spectra.copy() #self.bloque0
699 in_sat_cspectra = cspectra.copy() #self.bloques
700
701 raxs = math.ceil(math.sqrt(self.nPairs))
702 caxs = math.ceil(self.nPairs/raxs)
703
704 #print(self.hval)
705 #print numpy.absolute(rfunc[:,0,0,14])
706 for ih in range(self.minAltInd,self.maxAltInd):
707 for ifreq in range(self.nFFTPoints):
708 # fig, axs = plt.subplots(raxs, caxs)
709 # fig2, axs2 = plt.subplots(raxs, caxs)
710 col_ax = 0
711 row_ax = 0
712 for ii in range(self.nPairs): #PARES DE CANALES SELF y CROSS
713 #print("ii: ",ii)
714 if (col_ax%caxs==0 and col_ax!=0):
715 col_ax = 0
716 row_ax += 1
717 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
718 #print(func2clean.shape)
719 val = (numpy.isfinite(func2clean)==True).nonzero()
720
721 if len(val)>0: #limitador
722 min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
723 if min_val <= -40 :
724 min_val = -40
725 max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
726 if max_val >= 200 :
727 max_val = 200
728 #print min_val, max_val
729 step = 1
730 #print("Getting bins and the histogram")
731 x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
732 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
733 #print(len(y_dist),len(binstep[:-1]))
734 #print(row_ax,col_ax, " ..")
735 #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
736 mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
737 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
738 parg = [numpy.amax(y_dist),mean,sigma]
739 gauss_fit, covariance = None, None
740 newY = None
741 try :
742 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
743 mode = gauss_fit[1]
744 stdv = gauss_fit[2]
745 #print(" FIT OK",gauss_fit)
746 '''
747 newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
748 axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
749 axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
750 axs[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))'''
751 except:
752 mode = mean
753 stdv = sigma
754 #print("FIT FAIL")
755
756
757 #print(mode,stdv)
758 #Removing echoes greater than mode + 3*stdv
759 #factor_stdv = 2
760 noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
761 #noval tiene los indices que se van a remover
762 #print("Pair ",ii," novals: ",len(noval[0]))
763 if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
764 novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
765 #print(novall)
766 #print(" ")
767 cross_pairs = self.pairsArray[ii]
768 #Getting coherent echoes which are removed.
769 if len(novall[0]) > 0:
770
771 val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
772 val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
773 val_cspc[novall[0],ii,ifreq,ih] = 1
774 #print("OUT NOVALL 1")
775 #Removing coherent from ISR data
776
777 #print(spectra[:,ii,ifreq,ih])
778 new_a = numpy.delete(cspectra[:,ii,ifreq,ih], noval[0])
779 mean_cspc = numpy.mean(new_a)
780 new_b = numpy.delete(spectra[:,cross_pairs[0],ifreq,ih], noval[0])
781 mean_spc0 = numpy.mean(new_b)
782 new_c = numpy.delete(spectra[:,cross_pairs[1],ifreq,ih], noval[0])
783 mean_spc1 = numpy.mean(new_c)
784 spectra[noval,cross_pairs[0],ifreq,ih] = mean_spc0
785 spectra[noval,cross_pairs[1],ifreq,ih] = mean_spc1
786 cspectra[noval,ii,ifreq,ih] = mean_cspc
787
788 '''
789 func2clean = 10*numpy.log10(numpy.absolute(cspectra[:,ii,ifreq,ih]))
790 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
791 axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
792 axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
793 axs2[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))
794 '''
795
796 col_ax += 1 #contador de ploteo columnas
797 ##print(col_ax)
798 '''
799 title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
800 title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
801 fig.suptitle(title)
802 fig2.suptitle(title2)
803 plt.show()'''
804
805 ''' channels = channels
806 cross_pairs = cross_pairs
807 #print("OUT NOVALL 2")
808
809 vcross0 = (cross_pairs[0] == channels[ii]).nonzero()
810 vcross1 = (cross_pairs[1] == channels[ii]).nonzero()
811 vcross = numpy.concatenate((vcross0,vcross1),axis=None)
812 #print('vcros =', vcross)
813
814 #Getting coherent echoes which are removed.
815 if len(novall) > 0:
816 #val_spc[novall,ii,ifreq,ih] = 1
817 val_spc[ii,ifreq,ih,novall] = 1
818 if len(vcross) > 0:
819 val_cspc[vcross,ifreq,ih,novall] = 1
820
821 #Removing coherent from ISR data.
822 self.bloque0[ii,ifreq,ih,noval] = numpy.nan
823 if len(vcross) > 0:
824 self.bloques[vcross,ifreq,ih,noval] = numpy.nan
825 '''
826
827 print("Getting average of the spectra and cross-spectra from incoherent echoes.")
828 out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
829 out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
830 for ih in range(self.nHeights):
831 for ifreq in range(self.nFFTPoints):
832 for ich in range(self.nChan):
833 tmp = spectra[:,ich,ifreq,ih]
834 valid = (numpy.isfinite(tmp[:])==True).nonzero()
835 # if ich == 0 and ifreq == 0 and ih == 17 :
836 # print tmp
837 # print valid
838 # print len(valid[0])
839 #print('TMP',tmp)
840 if len(valid[0]) >0 :
841 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
842 #for icr in range(nPairs):
843 for icr in range(self.nPairs):
844 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
845 valid = (numpy.isfinite(tmp)==True).nonzero()
846 if len(valid[0]) > 0:
847 out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
848 '''
849 # print('##########################################################')
850 print("Removing fake coherent echoes (at least 4 points around the point)")
851
852 val_spectra = numpy.sum(val_spc,0)
853 val_cspectra = numpy.sum(val_cspc,0)
854
855 val_spectra = self.REM_ISOLATED_POINTS(val_spectra,4)
856 val_cspectra = self.REM_ISOLATED_POINTS(val_cspectra,4)
857
858 for i in range(nChan):
859 for j in range(nProf):
860 for k in range(nHeights):
861 if numpy.isfinite(val_spectra[i,j,k]) and val_spectra[i,j,k] < 1 :
862 val_spc[:,i,j,k] = 0.0
863 for i in range(nPairs):
864 for j in range(nProf):
865 for k in range(nHeights):
866 if numpy.isfinite(val_cspectra[i,j,k]) and val_cspectra[i,j,k] < 1 :
867 val_cspc[:,i,j,k] = 0.0
868
869 # val_spc = numpy.reshape(val_spc, (len(spectra[:,0,0,0]),nProf*nHeights*nChan))
870 # if numpy.isfinite(val_spectra)==str(True):
871 # noval = (val_spectra<1).nonzero()
872 # if len(noval) > 0:
873 # val_spc[:,noval] = 0.0
874 # val_spc = numpy.reshape(val_spc, (149,nChan,nProf,nHeights))
875
876 #val_cspc = numpy.reshape(val_spc, (149,nChan*nHeights*nProf))
877 #if numpy.isfinite(val_cspectra)==str(True):
878 # noval = (val_cspectra<1).nonzero()
879 # if len(noval) > 0:
880 # val_cspc[:,noval] = 0.0
881 # val_cspc = numpy.reshape(val_cspc, (149,nChan,nProf,nHeights))
882 tmp_sat_spectra = spectra.copy()
883 tmp_sat_spectra = tmp_sat_spectra*numpy.nan
884 tmp_sat_cspectra = cspectra.copy()
885 tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan
886 '''
887 # fig = plt.figure(figsize=(6,5))
888 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
889 # ax = fig.add_axes([left, bottom, width, height])
890 # cp = ax.contour(10*numpy.log10(numpy.absolute(spectra[0,0,:,:])))
891 # ax.clabel(cp, inline=True,fontsize=10)
892 # plt.show()
893 '''
894 val = (val_spc > 0).nonzero()
895 if len(val[0]) > 0:
896 tmp_sat_spectra[val] = in_sat_spectra[val]
897 val = (val_cspc > 0).nonzero()
898 if len(val[0]) > 0:
899 tmp_sat_cspectra[val] = in_sat_cspectra[val]
900
901 print("Getting average of the spectra and cross-spectra from incoherent echoes 2")
902 sat_spectra = numpy.zeros((nChan,nProf,nHeights), dtype=float)
903 sat_cspectra = numpy.zeros((nPairs,nProf,nHeights), dtype=complex)
904 for ih in range(nHeights):
905 for ifreq in range(nProf):
906 for ich in range(nChan):
907 tmp = numpy.squeeze(tmp_sat_spectra[:,ich,ifreq,ih])
908 valid = (numpy.isfinite(tmp)).nonzero()
909 if len(valid[0]) > 0:
910 sat_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
911
912 for icr in range(nPairs):
913 tmp = numpy.squeeze(tmp_sat_cspectra[:,icr,ifreq,ih])
914 valid = (numpy.isfinite(tmp)).nonzero()
915 if len(valid[0]) > 0:
916 sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0])
917 '''
918 #self.__dataReady= True
919 #sat_spectra, sat_cspectra= sat_spectra, sat_cspectra
920 #if not self.__dataReady:
921 #return None, None
922 #return out_spectra, out_cspectra ,sat_spectra,sat_cspectra
923 return out_spectra, out_cspectra
924
925 def REM_ISOLATED_POINTS(self,array,rth):
926 # import matplotlib.pyplot as plt
927 if rth == None :
928 rth = 4
929 print("REM ISO")
930 num_prof = len(array[0,:,0])
931 num_hei = len(array[0,0,:])
932 n2d = len(array[:,0,0])
933
934 for ii in range(n2d) :
935 #print ii,n2d
936 tmp = array[ii,:,:]
937 #print tmp.shape, array[ii,101,:],array[ii,102,:]
938
939 # fig = plt.figure(figsize=(6,5))
940 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
941 # ax = fig.add_axes([left, bottom, width, height])
942 # x = range(num_prof)
943 # y = range(num_hei)
944 # cp = ax.contour(y,x,tmp)
945 # ax.clabel(cp, inline=True,fontsize=10)
946 # plt.show()
947
948 #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
949 tmp = numpy.reshape(tmp,num_prof*num_hei)
950 indxs1 = (numpy.isfinite(tmp)==True).nonzero()
951 indxs2 = (tmp > 0).nonzero()
952
953 indxs1 = (indxs1[0])
954 indxs2 = indxs2[0]
955 #indxs1 = numpy.array(indxs1[0])
956 #indxs2 = numpy.array(indxs2[0])
957 indxs = None
958 #print indxs1 , indxs2
959 for iv in range(len(indxs2)):
960 indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
961 #print len(indxs2), indv
962 if len(indv[0]) > 0 :
963 indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None)
964 # print indxs
965 indxs = indxs[1:]
966 #print(indxs, len(indxs))
967 if len(indxs) < 4 :
968 array[ii,:,:] = 0.
969 return
970
971 xpos = numpy.mod(indxs ,num_hei)
972 ypos = (indxs / num_hei)
973 sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
974 #print sx
975 xpos = xpos[sx]
976 ypos = ypos[sx]
977
978 # *********************************** Cleaning isolated points **********************************
979 ic = 0
980 while True :
981 r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
982 #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
983 #plt.plot(r)
984 #plt.show()
985 no_coh1 = (numpy.isfinite(r)==True).nonzero()
986 no_coh2 = (r <= rth).nonzero()
987 #print r, no_coh1, no_coh2
988 no_coh1 = numpy.array(no_coh1[0])
989 no_coh2 = numpy.array(no_coh2[0])
990 no_coh = None
991 #print valid1 , valid2
992 for iv in range(len(no_coh2)):
993 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
994 if len(indv[0]) > 0 :
995 no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
996 no_coh = no_coh[1:]
997 #print len(no_coh), no_coh
998 if len(no_coh) < 4 :
999 #print xpos[ic], ypos[ic], ic
1000 # plt.plot(r)
1001 # plt.show()
1002 xpos[ic] = numpy.nan
1003 ypos[ic] = numpy.nan
1004
1005 ic = ic + 1
1006 if (ic == len(indxs)) :
1007 break
1008 #print( xpos, ypos)
1009
1010 indxs = (numpy.isfinite(list(xpos))==True).nonzero()
1011 #print indxs[0]
1012 if len(indxs[0]) < 4 :
1013 array[ii,:,:] = 0.
1014 return
1015
1016 xpos = xpos[indxs[0]]
1017 ypos = ypos[indxs[0]]
1018 for i in range(0,len(ypos)):
1019 ypos[i]=int(ypos[i])
1020 junk = tmp
1021 tmp = junk*0.0
1022
1023 tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
1024 array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
1025
1026 #print array.shape
1027 #tmp = numpy.reshape(tmp,(num_prof,num_hei))
1028 #print tmp.shape
1029
1030 # fig = plt.figure(figsize=(6,5))
1031 # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
1032 # ax = fig.add_axes([left, bottom, width, height])
1033 # x = range(num_prof)
1034 # y = range(num_hei)
1035 # cp = ax.contour(y,x,array[ii,:,:])
1036 # ax.clabel(cp, inline=True,fontsize=10)
1037 # plt.show()
1038 return array
1039
481 1040 class removeInterference(Operation):
482 1041
483 1042 def removeInterference2(self):
General Comments 0
You need to be logged in to leave comments. Login now