From 83c0b1494611f9a4637d19e614e8c027fe9da7ae 2021-10-05 14:46:32 From: joabAM Date: 2021-10-05 14:46:32 Subject: [PATCH] Adicionada operación Clean Rayleigh a Spectra Proc --- diff --git a/schainpy/model/graphics/jroplot_heispectra.py b/schainpy/model/graphics/jroplot_heispectra.py index a98a5cd..da391ea 100644 --- a/schainpy/model/graphics/jroplot_heispectra.py +++ b/schainpy/model/graphics/jroplot_heispectra.py @@ -9,7 +9,7 @@ import numpy from schainpy.model.graphics.jroplot_base import Plot, plt - +import matplotlib.pyplot as plt class SpectraHeisPlot(Plot): @@ -33,17 +33,18 @@ class SpectraHeisPlot(Plot): meta = {} spc = 10*numpy.log10(dataOut.data_spc / dataOut.normFactor) data['spc_heis'] = spc - - return data, meta + + return data, meta def plot(self): c = 3E8 - deltaHeight = self.data.yrange[1] - self.data.yrange[0] + deltaHeight = self.data.yrange[1] - self.data.yrange[0] # yrange = heightList x = numpy.arange(-1*len(self.data.yrange)/2., len(self.data.yrange)/2.)*(c/(2*deltaHeight*len(self.data.yrange)*1000)) self.y = self.data[-1]['spc_heis'] self.titles = [] + Maintitle = "Range from %d km to %d km" %(int(self.data.yrange[0]),int(self.data.yrange[-1])) for n, ax in enumerate(self.axes): ychannel = self.y[n,:] if ax.firsttime: @@ -54,7 +55,7 @@ class SpectraHeisPlot(Plot): ax.plt.set_data(x, ychannel) self.titles.append("Channel {}: {:4.2f}dB".format(n, numpy.max(ychannel))) - + plt.suptitle(Maintitle) class RTIHeisPlot(Plot): @@ -80,8 +81,8 @@ class RTIHeisPlot(Plot): spc = dataOut.data_spc / dataOut.normFactor spc = 10*numpy.log10(numpy.average(spc, axis=1)) data['rti_heis'] = spc - - return data, meta + + return data, meta def plot(self): diff --git a/schainpy/model/graphics/jroplot_voltage.py b/schainpy/model/graphics/jroplot_voltage.py index 6faf42a..c0d2664 100644 --- a/schainpy/model/graphics/jroplot_voltage.py +++ b/schainpy/model/graphics/jroplot_voltage.py @@ -22,11 +22,11 @@ class ScopePlot(Plot): def setup(self): self.xaxis = 'Range (Km)' - self.ncols = 1 - self.nrows = 1 - self.nplots = 1 + self.nplots = len(self.data.channels) + self.nrows = int(numpy.ceil(self.nplots/2)) + self.ncols = int(numpy.ceil(self.nplots/self.nrows)) self.ylabel = 'Intensity [dB]' - self.titles = ['Scope'] + self.titles = ['Channel '+str(self.data.channels[i])+" " for i in self.data.channels] self.colorbar = False self.width = 6 self.height = 4 @@ -56,15 +56,13 @@ class ScopePlot(Plot): yreal = y[channelIndexList,:].real yimag = y[channelIndexList,:].imag - title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y")) + Maintitle = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y")) self.xlabel = "Range (Km)" self.ylabel = "Intensity - IQ" self.y = yreal self.x = x - self.titles[0] = title - for i,ax in enumerate(self.axes): title = "Channel %d" %(i) if ax.firsttime: @@ -75,6 +73,7 @@ class ScopePlot(Plot): else: ax.plt_r.set_data(x, yreal[i,:]) ax.plt_i.set_data(x, yimag[i,:]) + plt.suptitle(Maintitle) def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle): y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:]) @@ -99,6 +98,7 @@ class ScopePlot(Plot): else: ax.plt_r.set_data(x, ychannel) + def plot_weatherpower(self, x, y, channelIndexList, thisDatetime, wintitle): diff --git a/schainpy/model/io/jroIO_base.py b/schainpy/model/io/jroIO_base.py index 7d17366..f3fadec 100644 --- a/schainpy/model/io/jroIO_base.py +++ b/schainpy/model/io/jroIO_base.py @@ -384,7 +384,7 @@ def isRadarFolder(folder): def isRadarFile(file): - try: + try: year = int(file[1:5]) doy = int(file[5:8]) set = int(file[8:11]) @@ -395,10 +395,10 @@ def isRadarFile(file): def getDateFromRadarFile(file): - try: + try: year = int(file[1:5]) doy = int(file[5:8]) - set = int(file[8:11]) + set = int(file[8:11]) except: return None @@ -417,7 +417,7 @@ def getDateFromRadarFolder(folder): return thisDate def parse_format(s, fmt): - + for i in range(fmt.count('%')): x = fmt.index('%') d = DT_DIRECTIVES[fmt[x:x+2]] @@ -484,7 +484,7 @@ class Reader(object): def run(self): - raise NotImplementedError + raise NotImplementedError def getAllowedArgs(self): if hasattr(self, '__attrs__'): @@ -496,19 +496,19 @@ class Reader(object): for key, value in kwargs.items(): setattr(self, key, value) - + def find_folders(self, path, startDate, endDate, folderfmt, last=False): - folders = [x for f in path.split(',') + folders = [x for f in path.split(',') for x in os.listdir(f) if os.path.isdir(os.path.join(f, x))] folders.sort() if last: folders = [folders[-1]] - for folder in folders: - try: - dt = datetime.datetime.strptime(parse_format(folder, folderfmt), folderfmt).date() + for folder in folders: + try: + dt = datetime.datetime.strptime(parse_format(folder, folderfmt), folderfmt).date() if dt >= startDate and dt <= endDate: yield os.path.join(path, folder) else: @@ -517,38 +517,38 @@ class Reader(object): log.log('Skiping folder {}'.format(folder), self.name) continue return - - def find_files(self, folders, ext, filefmt, startDate=None, endDate=None, + + def find_files(self, folders, ext, filefmt, startDate=None, endDate=None, expLabel='', last=False): - - for path in folders: + + for path in folders: files = glob.glob1(path, '*{}'.format(ext)) files.sort() if last: - if files: + if files: fo = files[-1] - try: + try: dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date() - yield os.path.join(path, expLabel, fo) - except Exception as e: + yield os.path.join(path, expLabel, fo) + except Exception as e: pass return else: return for fo in files: - try: - dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date() + try: + dt = datetime.datetime.strptime(parse_format(fo, filefmt), filefmt).date() if dt >= startDate and dt <= endDate: yield os.path.join(path, expLabel, fo) else: log.log('Skiping file {}'.format(fo), self.name) except Exception as e: log.log('Skiping file {}'.format(fo), self.name) - continue + continue def searchFilesOffLine(self, path, startDate, endDate, - expLabel, ext, walk, + expLabel, ext, walk, filefmt, folderfmt): """Search files in offline mode for the given arguments @@ -561,12 +561,12 @@ class Reader(object): path, startDate, endDate, folderfmt) else: folders = path.split(',') - + return self.find_files( - folders, ext, filefmt, startDate, endDate, expLabel) + folders, ext, filefmt, startDate, endDate, expLabel) def searchFilesOnLine(self, path, startDate, endDate, - expLabel, ext, walk, + expLabel, ext, walk, filefmt, folderfmt): """Search for the last file of the last folder @@ -579,13 +579,13 @@ class Reader(object): Return: generator with the full path of last filename """ - + if walk: folders = self.find_folders( path, startDate, endDate, folderfmt, last=True) else: folders = path.split(',') - + return self.find_files( folders, ext, filefmt, startDate, endDate, expLabel, last=True) @@ -594,13 +594,13 @@ class Reader(object): while True: if self.fp != None: - self.fp.close() + self.fp.close() if self.online: newFile = self.setNextFileOnline() else: newFile = self.setNextFileOffline() - + if not(newFile): if self.online: raise schainpy.admin.SchainError('Time to wait for new files reach') @@ -609,10 +609,10 @@ class Reader(object): raise schainpy.admin.SchainWarning('No files found in the given path') else: raise schainpy.admin.SchainWarning('No more files to read') - + if self.verifyFile(self.filename): break - + log.log('Opening file: %s' % self.filename, self.name) self.readFirstHeader() @@ -625,7 +625,7 @@ class Reader(object): self.filename self.fp self.filesize - + Return: boolean @@ -633,7 +633,7 @@ class Reader(object): nextFile = True nextDay = False - for nFiles in range(self.nFiles+1): + for nFiles in range(self.nFiles+1): for nTries in range(self.nTries): fullfilename, filename = self.checkForRealPath(nextFile, nextDay) if fullfilename is not None: @@ -643,18 +643,18 @@ class Reader(object): self.name) time.sleep(self.delay) nextFile = False - continue - + continue + if fullfilename is not None: break - + self.nTries = 1 - nextFile = True + nextFile = True if nFiles == (self.nFiles - 1): log.log('Trying with next day...', self.name) nextDay = True - self.nTries = 3 + self.nTries = 3 if fullfilename: self.fileSize = os.path.getsize(fullfilename) @@ -666,18 +666,18 @@ class Reader(object): self.flagNoMoreFiles = 0 self.fileIndex += 1 return 1 - else: + else: return 0 - + def setNextFileOffline(self): """Open the next file to be readed in offline mode""" - + try: filename = next(self.filenameList) self.fileIndex +=1 except StopIteration: self.flagNoMoreFiles = 1 - return 0 + return 0 self.filename = filename self.fileSize = os.path.getsize(filename) @@ -685,22 +685,22 @@ class Reader(object): self.flagIsNewFile = 1 return 1 - + @staticmethod def isDateTimeInRange(dt, startDate, endDate, startTime, endTime): """Check if the given datetime is in range""" - - if startDate <= dt.date() <= endDate: - if startTime <= dt.time() <= endTime: - return True + startDateTime= datetime.datetime.combine(startDate,startTime) + endDateTime = datetime.datetime.combine(endDate,endTime) + if startDateTime <= dt <= endDateTime: + return True return False - + def verifyFile(self, filename): """Check for a valid file - + Arguments: filename -- full path filename - + Return: boolean """ @@ -711,7 +711,7 @@ class Reader(object): """Check if the next file to be readed exists""" raise NotImplementedError - + def readFirstHeader(self): """Parse the file header""" @@ -783,8 +783,8 @@ class JRODataReader(Reader): Return: str -- fullpath of the file """ - - + + if nextFile: self.set += 1 if nextDay: @@ -796,7 +796,7 @@ class JRODataReader(Reader): prefixFileList = ['d', 'D'] elif self.ext.lower() == ".pdata": # spectra prefixFileList = ['p', 'P'] - + # barrido por las combinaciones posibles for prefixDir in prefixDirList: thispath = self.path @@ -816,9 +816,9 @@ class JRODataReader(Reader): if os.path.exists(fullfilename): return fullfilename, filename - - return None, filename - + + return None, filename + def __waitNewBlock(self): """ Return 1 si se encontro un nuevo bloque de datos, 0 de otra forma. @@ -860,9 +860,9 @@ class JRODataReader(Reader): def __setNewBlock(self): if self.fp == None: - return 0 - - if self.flagIsNewFile: + return 0 + + if self.flagIsNewFile: self.lastUTTime = self.basicHeaderObj.utc return 1 @@ -875,12 +875,12 @@ class JRODataReader(Reader): currentSize = self.fileSize - self.fp.tell() neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize - + if (currentSize >= neededSize): self.basicHeaderObj.read(self.fp) self.lastUTTime = self.basicHeaderObj.utc return 1 - + if self.__waitNewBlock(): self.lastUTTime = self.basicHeaderObj.utc return 1 @@ -966,10 +966,10 @@ class JRODataReader(Reader): except IOError: log.error("File {} can't be opened".format(filename), self.name) return False - + if self.online and self.waitDataBlock(0): pass - + basicHeaderObj = BasicHeader(LOCALTIME) systemHeaderObj = SystemHeader() radarControllerHeaderObj = RadarControllerHeader() @@ -996,7 +996,7 @@ class JRODataReader(Reader): dt2 = basicHeaderObj.datatime if not self.isDateTimeInRange(dt1, self.startDate, self.endDate, self.startTime, self.endTime) and not \ self.isDateTimeInRange(dt2, self.startDate, self.endDate, self.startTime, self.endTime): - flag = False + flag = False fp.close() return flag @@ -1105,11 +1105,11 @@ class JRODataReader(Reader): return dateList def setup(self, **kwargs): - + self.set_kwargs(**kwargs) if not self.ext.startswith('.'): self.ext = '.{}'.format(self.ext) - + if self.server is not None: if 'tcp://' in self.server: address = server @@ -1131,36 +1131,36 @@ class JRODataReader(Reader): for nTries in range(self.nTries): fullpath = self.searchFilesOnLine(self.path, self.startDate, - self.endDate, self.expLabel, self.ext, self.walk, + self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt) try: fullpath = next(fullpath) except: fullpath = None - + if fullpath: break log.warning( 'Waiting {} sec for a valid file in {}: try {} ...'.format( - self.delay, self.path, nTries + 1), + self.delay, self.path, nTries + 1), self.name) time.sleep(self.delay) if not(fullpath): raise schainpy.admin.SchainError( - 'There isn\'t any valid file in {}'.format(self.path)) + 'There isn\'t any valid file in {}'.format(self.path)) pathname, filename = os.path.split(fullpath) self.year = int(filename[1:5]) self.doy = int(filename[5:8]) - self.set = int(filename[8:11]) - 1 + self.set = int(filename[8:11]) - 1 else: log.log("Searching files in {}".format(self.path), self.name) - self.filenameList = self.searchFilesOffLine(self.path, self.startDate, + self.filenameList = self.searchFilesOffLine(self.path, self.startDate, self.endDate, self.expLabel, self.ext, self.walk, self.filefmt, self.folderfmt) - + self.setNextFile() return @@ -1181,7 +1181,7 @@ class JRODataReader(Reader): self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs - + def getFirstHeader(self): raise NotImplementedError @@ -1214,8 +1214,8 @@ class JRODataReader(Reader): """ Arguments: - path : - startDate : + path : + startDate : endDate : startTime : endTime : @@ -1284,7 +1284,7 @@ class JRODataWriter(Reader): dtype_width = get_dtype_width(dtype_index) return dtype_width - + def getProcessFlags(self): processFlags = 0 @@ -1322,9 +1322,9 @@ class JRODataWriter(Reader): self.basicHeaderObj.size = self.basicHeaderSize # bytes self.basicHeaderObj.version = self.versionFile - self.basicHeaderObj.dataBlock = self.nTotalBlocks + self.basicHeaderObj.dataBlock = self.nTotalBlocks utc = numpy.floor(self.dataOut.utctime) - milisecond = (self.dataOut.utctime - utc) * 1000.0 + milisecond = (self.dataOut.utctime - utc) * 1000.0 self.basicHeaderObj.utc = utc self.basicHeaderObj.miliSecond = milisecond self.basicHeaderObj.timeZone = self.dataOut.timeZone @@ -1465,9 +1465,9 @@ class JRODataWriter(Reader): if self.dataOut.datatime.date() > self.fileDate: setFile = 0 self.nTotalBlocks = 0 - + filen = '{}{:04d}{:03d}{:03d}{}'.format( - self.optchar, timeTuple.tm_year, timeTuple.tm_yday, setFile, ext) + self.optchar, timeTuple.tm_year, timeTuple.tm_yday, setFile, ext) filename = os.path.join(path, subfolder, filen) @@ -1515,11 +1515,11 @@ class JRODataWriter(Reader): self.ext = ext.lower() self.path = path - + if set is None: self.setFile = -1 else: - self.setFile = set - 1 + self.setFile = set - 1 self.blocksPerFile = blocksPerFile self.profilesPerBlock = profilesPerBlock diff --git a/schainpy/model/io/jroIO_kamisr.py b/schainpy/model/io/jroIO_kamisr.py index e90a2b1..bb97e31 100644 --- a/schainpy/model/io/jroIO_kamisr.py +++ b/schainpy/model/io/jroIO_kamisr.py @@ -347,7 +347,7 @@ class AMISRReader(ProcessingUnit): else: #get the last file - 1 - self.filenameList = [self.filenameList[-1]] + self.filenameList = [self.filenameList[-2]] new_dirnameList = [] for dirname in self.dirnameList: junk = numpy.array([dirname in x for x in self.filenameList]) @@ -421,7 +421,7 @@ class AMISRReader(ProcessingUnit): self.__selectDataForTimes(online=True) filename = self.filenameList[0] wait = 0 - #self.__waitForNewFile=180 ## DEBUG: + self.__waitForNewFile=300 ## DEBUG: while self.__filename_online == filename: print('waiting %d seconds to get a new file...'%(self.__waitForNewFile)) if wait == 5: @@ -641,13 +641,14 @@ class AMISRReader(ProcessingUnit): #self.dataOut.utctime = self.timeset[self.profileIndex] + t_comp #print(self.dataOut.utctime) self.dataOut.profileIndex = self.profileIndex + #print("N profile:",self.profileIndex,self.newProfiles,self.nblocks,self.dataOut.utctime) self.dataOut.flagNoData = False # if indexprof == 0: # print self.dataOut.utctime self.profileIndex += 1 - return self.dataOut.data + #return self.dataOut.data def run(self, **kwargs): @@ -660,3 +661,4 @@ class AMISRReader(ProcessingUnit): self.isConfig = True self.getData() + return(self.dataOut.data) diff --git a/schainpy/model/proc/jroproc_heispectra.py b/schainpy/model/proc/jroproc_heispectra.py index be414c1..fc74dc2 100644 --- a/schainpy/model/proc/jroproc_heispectra.py +++ b/schainpy/model/proc/jroproc_heispectra.py @@ -5,7 +5,6 @@ from schainpy.model.data.jrodata import SpectraHeis from schainpy.utils import log - class SpectraHeisProc(ProcessingUnit): def __init__(self):#, **kwargs): @@ -89,7 +88,7 @@ class SpectraHeisProc(ProcessingUnit): if self.dataIn.type == "Fits": self.__updateObjFromFits() self.dataOut.flagNoData = False - return + return if self.dataIn.type == "SpectraHeis": self.dataOut.copy(self.dataIn) @@ -343,5 +342,5 @@ class IncohInt4SpectraHeis(Operation): # dataOut.timeInterval = dataOut.ippSeconds * dataOut.nIncohInt # dataOut.timeInterval = self.__timeInterval*self.n dataOut.flagNoData = False - - return dataOut \ No newline at end of file + + return dataOut diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index 192cafa..13b9546 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -12,12 +12,15 @@ import time import itertools import numpy +import math from schainpy.model.proc.jroproc_base import ProcessingUnit, MPDecorator, Operation from schainpy.model.data.jrodata import Spectra from schainpy.model.data.jrodata import hildebrand_sekhon from schainpy.utils import log +from scipy.optimize import curve_fit + class SpectraProc(ProcessingUnit): @@ -478,6 +481,562 @@ class removeDC(Operation): return self.dataOut +# import matplotlib.pyplot as plt + +def fit_func( x, a0, a1, a2): #, a3, a4, a5): + z = (x - a1) / a2 + y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2 + return y +class CleanRayleigh(Operation): + + def __init__(self): + + Operation.__init__(self) + self.i=0 + self.isConfig = False + self.__dataReady = False + self.__profIndex = 0 + self.byTime = False + self.byProfiles = False + + self.bloques = None + self.bloque0 = None + + self.index = 0 + + self.buffer = 0 + self.buffer2 = 0 + self.buffer3 = 0 + #self.min_hei = None + #self.max_hei = None + + def setup(self,dataOut,min_hei,max_hei,n, timeInterval,factor_stdv): + + self.nChannels = dataOut.nChannels + self.nProf = dataOut.nProfiles + self.nPairs = dataOut.data_cspc.shape[0] + self.pairsArray = numpy.array(dataOut.pairsList) + self.spectra = dataOut.data_spc + self.cspectra = dataOut.data_cspc + self.heights = dataOut.heightList #alturas totales + self.nHeights = len(self.heights) + self.min_hei = min_hei + self.max_hei = max_hei + if (self.min_hei == None): + self.min_hei = 0 + if (self.max_hei == None): + self.max_hei = dataOut.heightList[-1] + self.hval = ((self.max_hei>=self.heights) & (self.heights >= self.min_hei)).nonzero() + self.heightsClean = self.heights[self.hval] #alturas filtradas + self.hval = self.hval[0] # forma (N,), an solo N elementos -> Indices de alturas + self.nHeightsClean = len(self.heightsClean) + self.channels = dataOut.channelList + self.nChan = len(self.channels) + self.nIncohInt = dataOut.nIncohInt + self.__initime = dataOut.utctime + self.maxAltInd = self.hval[-1]+1 + self.minAltInd = self.hval[0] + + self.crosspairs = dataOut.pairsList + self.nPairs = len(self.crosspairs) + self.normFactor = dataOut.normFactor + self.nFFTPoints = dataOut.nFFTPoints + self.ippSeconds = dataOut.ippSeconds + self.currentTime = self.__initime + self.pairsArray = numpy.array(dataOut.pairsList) + self.factor_stdv = factor_stdv + + + + if n != None : + self.byProfiles = True + self.nIntProfiles = n + else: + self.__integrationtime = timeInterval + + self.__dataReady = False + self.isConfig = True + + + + def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5): + #print (dataOut.utctime) + if not self.isConfig : + #print("Setting config") + self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv) + #print("Config Done") + tini=dataOut.utctime + + if self.byProfiles: + if self.__profIndex == self.nIntProfiles: + self.__dataReady = True + else: + if (tini - self.__initime) >= self.__integrationtime: + #print(tini - self.__initime,self.__profIndex) + self.__dataReady = True + self.__initime = tini + + #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0): + + if self.__dataReady: + print("Data ready",self.__profIndex) + self.__profIndex = 0 + jspc = self.buffer + jcspc = self.buffer2 + #jnoise = self.buffer3 + self.buffer = dataOut.data_spc + self.buffer2 = dataOut.data_cspc + #self.buffer3 = dataOut.noise + self.currentTime = dataOut.utctime + if numpy.any(jspc) : + #print( jspc.shape, jcspc.shape) + jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights)) + jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights)) + self.__dataReady = False + #print( jspc.shape, jcspc.shape) + dataOut.flagNoData = False + else: + dataOut.flagNoData = True + self.__dataReady = False + return dataOut + else: + #print( len(self.buffer)) + if numpy.any(self.buffer): + self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0) + self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0) + self.buffer3 += dataOut.data_dc + else: + self.buffer = dataOut.data_spc + self.buffer2 = dataOut.data_cspc + self.buffer3 = dataOut.data_dc + #print self.index, self.fint + #print self.buffer2.shape + dataOut.flagNoData = True ## NOTE: ?? revisar LUEGO + self.__profIndex += 1 + return dataOut ## NOTE: REV + + # if self.index == 0 and self.fint == 1 : + # if jspc != None: + # print len(jspc), jspc.shape + # jspc= numpy.reshape(jspc,(4,128,63,len(jspc)/4)) + # print jspc.shape + # dataOut.flagNoData = True + # return dataOut + # if path != None: + # sys.path.append(path) + # self.library = importlib.import_module(file) + # + # To be inserted as a parameter + + #Set constants + #constants = self.library.setConstants(dataOut) + #dataOut.constants = constants + + #snrth= 20 + + #crosspairs = dataOut.groupList + #noise = dataOut.noise + #print( nProf,heights) + #print( jspc.shape, jspc.shape[0]) + #print noise + #print jnoise[len(jnoise)-1,:], numpy.nansum(jnoise,axis=0)/len(jnoise) + #jnoise = jnoise/N + #noise = numpy.nansum(jnoise,axis=0)#/len(jnoise) + #print( noise) + #power = numpy.sum(spectra, axis=1) + #print power[0,:] + #print("CROSSPAIRS",crosspairs) + #nPairs = len(crosspairs) + #print(numpy.shape(dataOut.data_spc)) + + #absc = dataOut.abscissaList[:-1] + + #print absc.shape + #nBlocks=149 + #print('spectra', spectra.shape) + #print('noise print', crosspairs) + #print('spectra', spectra.shape) + #print('cspectra', cspectra.shape) + #print numpy.array(dataOut.data_pre[1]).shape + #spec, cspec = self.__DiffCoherent(snrth, spectra, cspectra, nProf, heights,nChan, nHei, nPairs, channels, noise*nProf, crosspairs) + + + #index = tini.tm_hour*12+tini.tm_min/5 + # jspc = jspc/self.nFFTPoints/self.normFactor + # jcspc = jcspc/self.nFFTPoints/self.normFactor + + + + #dataOut.data_spc,dataOut.data_cspc = self.CleanRayleigh(dataOut,jspc,jcspc,crosspairs,heights,channels,nProf,nHei,nChan,nPairs,nIncohInt,nBlocks=nBlocks) + #tmp_spectra,tmp_cspectra,sat_spectra,sat_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.min_hei,self.max_hei) + tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv) + #jspectra = tmp_spectra*len(jspc[:,0,0,0]) + #jcspectra = tmp_cspectra*len(jspc[:,0,0,0]) + + dataOut.data_spc = tmp_spectra + dataOut.data_cspc = tmp_cspectra + dataOut.data_dc = self.buffer3 + dataOut.nIncohInt *= self.nIntProfiles + dataOut.utctime = self.currentTime #tiempo promediado + print("Time: ",time.localtime(dataOut.utctime)) + # dataOut.data_spc = sat_spectra + # dataOut.data_cspc = sat_cspectra + self.buffer = 0 + self.buffer2 = 0 + self.buffer3 = 0 + + return dataOut + + def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv): + print("OP cleanRayleigh") + #import matplotlib.pyplot as plt + #for k in range(149): + + rfunc = cspectra.copy() #self.bloques + val_spc = spectra*0.0 #self.bloque0*0.0 + val_cspc = cspectra*0.0 #self.bloques*0.0 + in_sat_spectra = spectra.copy() #self.bloque0 + in_sat_cspectra = cspectra.copy() #self.bloques + + raxs = math.ceil(math.sqrt(self.nPairs)) + caxs = math.ceil(self.nPairs/raxs) + + #print(self.hval) + #print numpy.absolute(rfunc[:,0,0,14]) + for ih in range(self.minAltInd,self.maxAltInd): + for ifreq in range(self.nFFTPoints): + # fig, axs = plt.subplots(raxs, caxs) + # fig2, axs2 = plt.subplots(raxs, caxs) + col_ax = 0 + row_ax = 0 + for ii in range(self.nPairs): #PARES DE CANALES SELF y CROSS + #print("ii: ",ii) + if (col_ax%caxs==0 and col_ax!=0): + col_ax = 0 + row_ax += 1 + func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia? + #print(func2clean.shape) + val = (numpy.isfinite(func2clean)==True).nonzero() + + if len(val)>0: #limitador + min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40) + if min_val <= -40 : + min_val = -40 + max_val = numpy.around(numpy.amax(func2clean)+2) #< 200 + if max_val >= 200 : + max_val = 200 + #print min_val, max_val + step = 1 + #print("Getting bins and the histogram") + x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step + y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step)) + #print(len(y_dist),len(binstep[:-1])) + #print(row_ax,col_ax, " ..") + #print(self.pairsArray[ii][0],self.pairsArray[ii][1]) + mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist) + sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist)) + parg = [numpy.amax(y_dist),mean,sigma] + gauss_fit, covariance = None, None + newY = None + try : + gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg) + mode = gauss_fit[1] + stdv = gauss_fit[2] + #print(" FIT OK",gauss_fit) + ''' + newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2]) + axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green') + axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red') + axs[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))''' + except: + mode = mean + stdv = sigma + #print("FIT FAIL") + + + #print(mode,stdv) + #Removing echoes greater than mode + 3*stdv + #factor_stdv = 2 + noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero() + #noval tiene los indices que se van a remover + #print("Pair ",ii," novals: ",len(noval[0])) + if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N) + novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero() + #print(novall) + #print(" ") + cross_pairs = self.pairsArray[ii] + #Getting coherent echoes which are removed. + if len(novall[0]) > 0: + + val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1 + val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1 + val_cspc[novall[0],ii,ifreq,ih] = 1 + #print("OUT NOVALL 1") + #Removing coherent from ISR data + + #print(spectra[:,ii,ifreq,ih]) + new_a = numpy.delete(cspectra[:,ii,ifreq,ih], noval[0]) + mean_cspc = numpy.mean(new_a) + new_b = numpy.delete(spectra[:,cross_pairs[0],ifreq,ih], noval[0]) + mean_spc0 = numpy.mean(new_b) + new_c = numpy.delete(spectra[:,cross_pairs[1],ifreq,ih], noval[0]) + mean_spc1 = numpy.mean(new_c) + spectra[noval,cross_pairs[0],ifreq,ih] = mean_spc0 + spectra[noval,cross_pairs[1],ifreq,ih] = mean_spc1 + cspectra[noval,ii,ifreq,ih] = mean_cspc + + ''' + func2clean = 10*numpy.log10(numpy.absolute(cspectra[:,ii,ifreq,ih])) + y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step)) + axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red') + axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green') + axs2[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii])) + ''' + + col_ax += 1 #contador de ploteo columnas + ##print(col_ax) + ''' + title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km" + title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED" + fig.suptitle(title) + fig2.suptitle(title2) + plt.show()''' + + ''' channels = channels + cross_pairs = cross_pairs + #print("OUT NOVALL 2") + + vcross0 = (cross_pairs[0] == channels[ii]).nonzero() + vcross1 = (cross_pairs[1] == channels[ii]).nonzero() + vcross = numpy.concatenate((vcross0,vcross1),axis=None) + #print('vcros =', vcross) + + #Getting coherent echoes which are removed. + if len(novall) > 0: + #val_spc[novall,ii,ifreq,ih] = 1 + val_spc[ii,ifreq,ih,novall] = 1 + if len(vcross) > 0: + val_cspc[vcross,ifreq,ih,novall] = 1 + + #Removing coherent from ISR data. + self.bloque0[ii,ifreq,ih,noval] = numpy.nan + if len(vcross) > 0: + self.bloques[vcross,ifreq,ih,noval] = numpy.nan + ''' + + print("Getting average of the spectra and cross-spectra from incoherent echoes.") + out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan + out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan + for ih in range(self.nHeights): + for ifreq in range(self.nFFTPoints): + for ich in range(self.nChan): + tmp = spectra[:,ich,ifreq,ih] + valid = (numpy.isfinite(tmp[:])==True).nonzero() +# if ich == 0 and ifreq == 0 and ih == 17 : +# print tmp +# print valid +# print len(valid[0]) + #print('TMP',tmp) + if len(valid[0]) >0 : + out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0]) + #for icr in range(nPairs): + for icr in range(self.nPairs): + tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih]) + valid = (numpy.isfinite(tmp)==True).nonzero() + if len(valid[0]) > 0: + out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0]) + ''' + # print('##########################################################') + print("Removing fake coherent echoes (at least 4 points around the point)") + + val_spectra = numpy.sum(val_spc,0) + val_cspectra = numpy.sum(val_cspc,0) + + val_spectra = self.REM_ISOLATED_POINTS(val_spectra,4) + val_cspectra = self.REM_ISOLATED_POINTS(val_cspectra,4) + + for i in range(nChan): + for j in range(nProf): + for k in range(nHeights): + if numpy.isfinite(val_spectra[i,j,k]) and val_spectra[i,j,k] < 1 : + val_spc[:,i,j,k] = 0.0 + for i in range(nPairs): + for j in range(nProf): + for k in range(nHeights): + if numpy.isfinite(val_cspectra[i,j,k]) and val_cspectra[i,j,k] < 1 : + val_cspc[:,i,j,k] = 0.0 + + # val_spc = numpy.reshape(val_spc, (len(spectra[:,0,0,0]),nProf*nHeights*nChan)) + # if numpy.isfinite(val_spectra)==str(True): + # noval = (val_spectra<1).nonzero() + # if len(noval) > 0: + # val_spc[:,noval] = 0.0 + # val_spc = numpy.reshape(val_spc, (149,nChan,nProf,nHeights)) + + #val_cspc = numpy.reshape(val_spc, (149,nChan*nHeights*nProf)) + #if numpy.isfinite(val_cspectra)==str(True): + # noval = (val_cspectra<1).nonzero() + # if len(noval) > 0: + # val_cspc[:,noval] = 0.0 + # val_cspc = numpy.reshape(val_cspc, (149,nChan,nProf,nHeights)) + tmp_sat_spectra = spectra.copy() + tmp_sat_spectra = tmp_sat_spectra*numpy.nan + tmp_sat_cspectra = cspectra.copy() + tmp_sat_cspectra = tmp_sat_cspectra*numpy.nan + ''' + # fig = plt.figure(figsize=(6,5)) + # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8 + # ax = fig.add_axes([left, bottom, width, height]) + # cp = ax.contour(10*numpy.log10(numpy.absolute(spectra[0,0,:,:]))) + # ax.clabel(cp, inline=True,fontsize=10) + # plt.show() + ''' + val = (val_spc > 0).nonzero() + if len(val[0]) > 0: + tmp_sat_spectra[val] = in_sat_spectra[val] + val = (val_cspc > 0).nonzero() + if len(val[0]) > 0: + tmp_sat_cspectra[val] = in_sat_cspectra[val] + + print("Getting average of the spectra and cross-spectra from incoherent echoes 2") + sat_spectra = numpy.zeros((nChan,nProf,nHeights), dtype=float) + sat_cspectra = numpy.zeros((nPairs,nProf,nHeights), dtype=complex) + for ih in range(nHeights): + for ifreq in range(nProf): + for ich in range(nChan): + tmp = numpy.squeeze(tmp_sat_spectra[:,ich,ifreq,ih]) + valid = (numpy.isfinite(tmp)).nonzero() + if len(valid[0]) > 0: + sat_spectra[ich,ifreq,ih] = numpy.nansum(tmp)/len(valid[0]) + + for icr in range(nPairs): + tmp = numpy.squeeze(tmp_sat_cspectra[:,icr,ifreq,ih]) + valid = (numpy.isfinite(tmp)).nonzero() + if len(valid[0]) > 0: + sat_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)/len(valid[0]) + ''' + #self.__dataReady= True + #sat_spectra, sat_cspectra= sat_spectra, sat_cspectra + #if not self.__dataReady: + #return None, None + #return out_spectra, out_cspectra ,sat_spectra,sat_cspectra + return out_spectra, out_cspectra + + def REM_ISOLATED_POINTS(self,array,rth): +# import matplotlib.pyplot as plt + if rth == None : + rth = 4 + print("REM ISO") + num_prof = len(array[0,:,0]) + num_hei = len(array[0,0,:]) + n2d = len(array[:,0,0]) + + for ii in range(n2d) : + #print ii,n2d + tmp = array[ii,:,:] + #print tmp.shape, array[ii,101,:],array[ii,102,:] + + # fig = plt.figure(figsize=(6,5)) + # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8 + # ax = fig.add_axes([left, bottom, width, height]) + # x = range(num_prof) + # y = range(num_hei) + # cp = ax.contour(y,x,tmp) + # ax.clabel(cp, inline=True,fontsize=10) + # plt.show() + + #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs) + tmp = numpy.reshape(tmp,num_prof*num_hei) + indxs1 = (numpy.isfinite(tmp)==True).nonzero() + indxs2 = (tmp > 0).nonzero() + + indxs1 = (indxs1[0]) + indxs2 = indxs2[0] + #indxs1 = numpy.array(indxs1[0]) + #indxs2 = numpy.array(indxs2[0]) + indxs = None + #print indxs1 , indxs2 + for iv in range(len(indxs2)): + indv = numpy.array((indxs1 == indxs2[iv]).nonzero()) + #print len(indxs2), indv + if len(indv[0]) > 0 : + indxs = numpy.concatenate((indxs,indxs2[iv]), axis=None) + # print indxs + indxs = indxs[1:] + #print(indxs, len(indxs)) + if len(indxs) < 4 : + array[ii,:,:] = 0. + return + + xpos = numpy.mod(indxs ,num_hei) + ypos = (indxs / num_hei) + sx = numpy.argsort(xpos) # Ordering respect to "x" (time) + #print sx + xpos = xpos[sx] + ypos = ypos[sx] + + # *********************************** Cleaning isolated points ********************************** + ic = 0 + while True : + r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2))) + #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh) + #plt.plot(r) + #plt.show() + no_coh1 = (numpy.isfinite(r)==True).nonzero() + no_coh2 = (r <= rth).nonzero() + #print r, no_coh1, no_coh2 + no_coh1 = numpy.array(no_coh1[0]) + no_coh2 = numpy.array(no_coh2[0]) + no_coh = None + #print valid1 , valid2 + for iv in range(len(no_coh2)): + indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero()) + if len(indv[0]) > 0 : + no_coh = numpy.concatenate((no_coh,no_coh2[iv]), axis=None) + no_coh = no_coh[1:] + #print len(no_coh), no_coh + if len(no_coh) < 4 : + #print xpos[ic], ypos[ic], ic + # plt.plot(r) + # plt.show() + xpos[ic] = numpy.nan + ypos[ic] = numpy.nan + + ic = ic + 1 + if (ic == len(indxs)) : + break + #print( xpos, ypos) + + indxs = (numpy.isfinite(list(xpos))==True).nonzero() + #print indxs[0] + if len(indxs[0]) < 4 : + array[ii,:,:] = 0. + return + + xpos = xpos[indxs[0]] + ypos = ypos[indxs[0]] + for i in range(0,len(ypos)): + ypos[i]=int(ypos[i]) + junk = tmp + tmp = junk*0.0 + + tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))] + array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei)) + + #print array.shape + #tmp = numpy.reshape(tmp,(num_prof,num_hei)) + #print tmp.shape + + # fig = plt.figure(figsize=(6,5)) + # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8 + # ax = fig.add_axes([left, bottom, width, height]) + # x = range(num_prof) + # y = range(num_hei) + # cp = ax.contour(y,x,array[ii,:,:]) + # ax.clabel(cp, inline=True,fontsize=10) + # plt.show() + return array + class removeInterference(Operation): def removeInterference2(self):