diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index 8ffee71..4221383 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -204,7 +204,7 @@ class JROData(GenericData): sampled_heightsFFT = None pulseLength_TxA = None deltaHeight = None - + def __str__(self): return '{} - {}'.format(self.type, self.datatime()) @@ -443,7 +443,7 @@ class Voltage(JROData): class Spectra(JROData): - data_outlier = 0 + data_outlier = None def __init__(self): ''' @@ -483,7 +483,6 @@ class Spectra(JROData): 'code', 'nCode', 'nBaud', 'ippSeconds', 'ipp','nIncohInt', 'nFFTPoints', 'nProfiles'] - self.max_nIncohInt = 1 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None): @@ -511,6 +510,7 @@ class Spectra(JROData): # # noise[channel] = hildebrand_sekhon(daux, self.nIncohInt) noise = numpy.zeros(self.nChannels) + for channel in range(self.nChannels): daux = self.data_spc[channel,xmin_index:xmax_index, ymin_index:ymax_index] @@ -576,6 +576,7 @@ class Spectra(JROData): if self.flagDecodeData: pwcode = numpy.sum(self.code[0]**2) + #print(self.flagDecodeData, pwcode) #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter @@ -612,15 +613,36 @@ class Spectra(JROData): factor = self.normFactor power = numpy.zeros( (self.nChannels,self.nHeights) ) for ch in range(self.nChannels): + z = None if hasattr(factor,'shape'): - z = numpy.divide(self.data_spc[ch],factor[ch]) + if factor.ndim > 1: + z = self.data_spc[ch]/factor[ch] + else: + z = self.data_spc[ch]/factor else: - z = numpy.divide(self.data_spc[ch],factor) + z = self.data_spc[ch]/factor z = numpy.where(numpy.isfinite(z), z, numpy.NAN) avg = numpy.average(z, axis=0) - power[ch,:] = 10 * numpy.log10(avg) + power[ch] = 10 * numpy.log10(avg) return power + @property + def max_nIncohInt(self): + + ints = numpy.zeros(self.nChannels) + for ch in range(self.nChannels): + if hasattr(self.nIncohInt,'shape'): + if self.nIncohInt.ndim > 1: + ints[ch,] = self.nIncohInt[ch].max() + else: + ints[ch,] = self.nIncohInt + self.nIncohInt = int(self.nIncohInt) + else: + ints[ch,] = self.nIncohInt + + return ints + + def getCoherence(self, pairsList=None, phase=False): z = [] @@ -901,7 +923,7 @@ class Parameters(Spectra): nAvg = None noise_estimation = None GauSPC = None # Fit gaussian SPC - max_nIncohInt = 1 + def __init__(self): ''' Constructor diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index 8a2e8dc..e31b912 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -56,12 +56,25 @@ class SpectraPlot(Plot): data = {} meta = {} - data['rti'] = dataOut.getPower() + #data['rti'] = dataOut.getPower() norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter - noise = 10*numpy.log10(dataOut.getNoise()/float(norm)) - data['noise'] = noise - spc = 10*numpy.log10(dataOut.data_spc/norm) + noise = 10*numpy.log10(dataOut.getNoise()/norm) + + + z = [] + for ch in range(dataOut.nChannels): + if hasattr(dataOut.normFactor,'shape'): + z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch])) + else: + z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor)) + + z = numpy.asarray(z) + z = numpy.where(numpy.isfinite(z), z, numpy.NAN) + spc = 10*numpy.log10(z) + data['spc'] = spc + data['rti'] = spc.mean(axis=1) + data['noise'] = noise meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) if self.CODE == 'spc_moments': data['moments'] = dataOut.moments @@ -270,10 +283,11 @@ class RTIPlot(Plot): self.update_list(dataOut) data = {} meta = {} + data['rti'] = dataOut.getPower() norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter - noise = 10*numpy.log10(dataOut.getNoise()/float(norm)) + noise = 10*numpy.log10(dataOut.getNoise()/norm) data['noise'] = noise return data, meta @@ -529,12 +543,23 @@ class SpectraCutPlot(Plot): data = {} meta = {} - norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter - n0 = 10*numpy.log10(dataOut.getNoise()/float(norm)) - - spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) + norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter#*dataOut.nFFTPoints + n0 = 10*numpy.log10(dataOut.getNoise()/norm) noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights) + + z = [] + for ch in range(dataOut.nChannels): + if hasattr(dataOut.normFactor,'shape'): + z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch])) + else: + z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor)) + + z = numpy.asarray(z) + z = numpy.where(numpy.isfinite(z), z, numpy.NAN) + spc = 10*numpy.log10(z) + + data['spc'] = spc - noise meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) @@ -860,18 +885,27 @@ class NoiselessSpectraPlot(Plot): data = {} meta = {} - norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter - n0 = 10*numpy.log10(dataOut.getNoise()/float(norm)) + norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter#*dataOut.nFFTPoints + n0 = 10*numpy.log10(dataOut.getNoise()/norm) + noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights) + z = [] + for ch in range(dataOut.nChannels): + if hasattr(dataOut.normFactor,'shape'): + z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor[ch])) + else: + z.append(numpy.divide(dataOut.data_spc[ch],dataOut.normFactor)) - #spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor) - spc = 10*numpy.log10(dataOut.data_spc/norm) + z = numpy.asarray(z) + z = numpy.where(numpy.isfinite(z), z, numpy.NAN) + spc = 10*numpy.log10(z) - noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights) - data['rti'] = dataOut.getPower() - noise - noise = numpy.repeat(n0,(dataOut.nFFTPoints*dataOut.nHeights)).reshape(dataOut.nChannels,dataOut.nFFTPoints,dataOut.nHeights) data['spc'] = spc - noise + #print(spc.shape) + data['rti'] = spc.mean(axis=1) + data['noise'] = noise + # data['noise'] = noise @@ -965,12 +999,12 @@ class NoiselessRTIPlot(Plot): self.update_list(dataOut) data = {} meta = {} + #print(dataOut.max_nIncohInt, dataOut.nIncohInt) + #print(dataOut.windowOfFilter,dataOut.nCohInt,dataOut.nProfiles,dataOut.max_nIncohInt,dataOut.nIncohInt) + norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter - norm = dataOut.nProfiles * dataOut.max_nIncohInt * dataOut.nCohInt * dataOut.windowOfFilter - #print("Norm: ", norm) - #print(dataOut.nProfiles , dataOut.max_nIncohInt ,dataOut.nCohInt, dataOut.windowOfFilter) - n0 = 10*numpy.log10(dataOut.getNoise()/float(norm)) + n0 = 10*numpy.log10(dataOut.getNoise()/norm) data['noise'] = n0 noise = numpy.repeat(n0,dataOut.nHeights).reshape(dataOut.nChannels,dataOut.nHeights) diff --git a/schainpy/model/io/utilsIO.py b/schainpy/model/io/utilsIO.py index 3813cae..2c806d9 100644 --- a/schainpy/model/io/utilsIO.py +++ b/schainpy/model/io/utilsIO.py @@ -274,6 +274,7 @@ class MergeH5(object): elif isinstance(meta[name], dict): for key, value in meta[name].items(): return value[x] + if 'cspc' in name: return 'pair{:02d}'.format(x) else: @@ -315,12 +316,25 @@ class MergeH5(object): self.dataOut.utctime = time ints = [data.nIncohInt for data in self.ch_dataIn] self.dataOut.nIncohInt = numpy.stack(ints, axis=1) + + print("nIncohInt 1: ",self.dataOut.nIncohInt.shape) + + if self.dataOut.nIncohInt.ndim > 3: + aux = self.dataOut.nIncohInt + self.dataOut.nIncohInt = None + self.dataOut.nIncohInt = aux[0] + + if self.dataOut.nIncohInt.ndim < 3: nIncohInt = numpy.repeat(self.dataOut.nIncohInt, self.dataOut.nHeights).reshape(self.blocksPerFile,self.nChannels, self.dataOut.nHeights) #nIncohInt = numpy.reshape(nIncohInt, (self.blocksPerFile,self.nChannels, self.dataOut.nHeights)) self.dataOut.nIncohInt = None self.dataOut.nIncohInt = nIncohInt - #print("nIncohInt: ", self.dataOut.nIncohInt.shape) + + if (self.dataOut.nIncohInt.shape)[0]==self.nChannels: ## ch,blocks, hei + self.dataOut.nIncohInt = numpy.swapaxes(self.dataOut.nIncohInt, 0, 1) ## blocks,ch, hei + + print("nIncohInt 2: ", self.dataOut.nIncohInt.shape) #print("utcTime: ", time.shape) #print("data_spc ",self.dataOut.data_spc.shape) pairsList = [pair for pair in itertools.combinations(self.channelList, 2)] @@ -371,8 +385,10 @@ class MergeH5(object): elif isinstance(dataAux, (int, float, numpy.integer, numpy.float)): dsDict['nDim'] = 0 else: + dsDict['nDim'] = len(dataAux.shape) -1 dsDict['shape'] = dataAux.shape + if len(dsDict['shape'])>=2: dsDict['dsNumber'] = dataAux.shape[1] else: @@ -399,6 +415,7 @@ class MergeH5(object): self.ch_dataIn[ch].utctime = None self.ch_dataIn[ch].nIncohInt = None self.meta ={} + self.blocksPerFile = None def writeData(self, outFilename): @@ -422,15 +439,17 @@ class MergeH5(object): sgrp = grp.create_group(label) else: sgrp = grp + k = -1*(dsInfo['nDim'] - 1) + #print(k, dsInfo['shape'], dsInfo['shape'][k:]) for i in range(dsInfo['dsNumber']): ds = sgrp.create_dataset( - self.getLabel(dsInfo['variable'], i), - (self.blocksPerFile, ) + dsInfo['shape'][2:], + self.getLabel(dsInfo['variable'], i),(self.blocksPerFile, ) + dsInfo['shape'][k:], chunks=True, dtype=dsInfo['dtype']) dtsets.append(ds) data.append((dsInfo['variable'], i)) + #print("\n",dtsets) print('Creating merged file: {}'.format(fp.filename)) @@ -442,6 +461,7 @@ class MergeH5(object): #print(ds, getattr(self.dataOut, attr)[ch].shape) aux = getattr(self.dataOut, attr)# block, ch, ... aux = numpy.swapaxes(aux,0,1) # ch, blocks, ... + #print(ds.shape, aux.shape) #ds[:] = getattr(self.dataOut, attr)[ch] ds[:] = aux[ch] @@ -468,7 +488,7 @@ class MergeH5(object): self.readFile(fp,ch) fp.close() self.getDataOut() - name = name[-16:-1] + name = name[-16:] #print("Final name out: ", name) outFile = os.path.join(self.pathOut, name) #print("Outfile: ", outFile) diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index cb71436..1b2268f 100755 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -141,7 +141,6 @@ class ParametersProc(ProcessingUnit): self.dataOut.nIncohInt = self.dataIn.nIncohInt self.dataOut.nFFTPoints = self.dataIn.nFFTPoints self.dataOut.ippFactor = self.dataIn.ippFactor - self.dataOut.max_nIncohInt = self.dataIn.max_nIncohInt self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy() self.dataOut.radarControllerHeaderTxt = self.dataOut.radarControllerHeaderObj.toString() self.dataOut.ipp = self.dataIn.ipp diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index 7847d9c..74f1aa6 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -63,7 +63,6 @@ class SpectraProc(ProcessingUnit): self.dataOut.flagShiftFFT = False self.dataOut.nCohInt = self.dataIn.nCohInt self.dataOut.nIncohInt = 1 - self.dataOut.max_nIncohInt = 1 self.dataOut.radar_ipp = self.dataIn.radar_ipp self.dataOut.sampled_heightsFFT = self.dataIn.sampled_heightsFFT self.dataOut.pulseLength_TxA = self.dataIn.pulseLength_TxA @@ -690,11 +689,12 @@ class getNoiseB(Operation): self.dataOut.noise_estimation = None noise = None + #print("data type: ",self.dataOut.type, self.dataOut.nIncohInt, self.dataOut.max_nIncohInt) if self.dataOut.type == 'Voltage': noise = self.dataOut.getNoise(ymin_index=self.minIndex, ymax_index=self.maxIndex) #print(minIndex, maxIndex,minIndexVel, maxIndexVel) elif self.dataOut.type == 'Spectra': - #print(self.dataOut.nChannels, self.minIndex, self.maxIndex,self.minIndexFFT, self.maxIndexFFT, self.dataOut.max_nIncohInt, self.dataOut.nIncohInt.shape) + #print(self.dataOut.nChannels, self.minIndex, self.maxIndex,self.minIndexFFT, self.maxIndexFFT, self.dataOut.max_nIncohInt, self.dataOut.nIncohInt) noise = numpy.zeros( self.dataOut.nChannels) norm = 1 @@ -702,8 +702,8 @@ class getNoiseB(Operation): if not hasattr(self.dataOut.nIncohInt,'__len__'): norm = 1 else: - norm = self.dataOut.max_nIncohInt/self.dataOut.nIncohInt[channel, self.minIndex:self.maxIndex] - #print("norm nIncoh: ", norm ,self.dataOut.data_spc.shape) + norm = self.dataOut.max_nIncohInt[channel]/self.dataOut.nIncohInt[channel, self.minIndex:self.maxIndex] + #print("norm nIncoh: ", norm ,self.dataOut.data_spc.shape, self.dataOut.max_nIncohInt) daux = self.dataOut.data_spc[channel,self.minIndexFFT:self.maxIndexFFT, self.minIndex:self.maxIndex] daux = numpy.multiply(daux, norm) #print("offset: ", self.offset, 10*numpy.log10(self.offset)) @@ -711,7 +711,8 @@ class getNoiseB(Operation): #print(daux.shape, daux) #noise[channel] = self.getNoiseByHS(daux, self.dataOut.max_nIncohInt)/self.offset sortdata = numpy.sort(daux, axis=None) - noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt)/self.offset + + noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt[channel])/self.offset #noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex) @@ -721,7 +722,7 @@ class getNoiseB(Operation): #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64)) #print("2: ",self.dataOut.noise_estimation) #print(self.dataOut.flagNoData) - # print("getNoise Done", noise, self.dataOut.nProfiles ,self.dataOut.ippFactor) + #print("getNoise Done", noise, self.dataOut.nProfiles ,self.dataOut.ippFactor) return self.dataOut def getNoiseByMean(self,data): @@ -798,464 +799,464 @@ def fit_func( x, a0, a1, a2): #, a3, a4, a5): 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 - - - 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("runing cleanRayleigh") - if not self.isConfig : - - self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv) - - tini=dataOut.utctime - - if self.byProfiles: - if self.__profIndex == self.nIntProfiles: - self.__dataReady = True - else: - if (tini - self.__initime) >= self.__integrationtime: - - self.__dataReady = True - self.__initime = tini - - #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0): - - if self.__dataReady: - - 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)) - try: - jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights)) - except: - print("no cspc") - 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) - try: - self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0) - self.buffer3 += dataOut.data_dc - except: - pass - 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 - - - #index = tini.tm_hour*12+tini.tm_min/5 - ''' - #REVISAR - ''' - # jspc = jspc/self.nFFTPoints/self.normFactor - # jcspc = jcspc/self.nFFTPoints/self.normFactor - - - - tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv) - dataOut.data_spc = tmp_spectra - dataOut.data_cspc = tmp_cspectra - - #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv) - - dataOut.data_dc = self.buffer3 - dataOut.nIncohInt *= self.nIntProfiles - dataOut.max_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): - #channelsProcssd = [] - #channelA_ok = False - #rfunc = cspectra.copy() #self.bloques - rfunc = spectra.copy() - #rfunc = cspectra - #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 - - - ###ONLY FOR TEST: - raxs = math.ceil(math.sqrt(self.nPairs)) - if raxs == 0: - raxs = 1 - caxs = math.ceil(self.nPairs/raxs) - if self.nPairs <4: - raxs = 2 - caxs = 2 - #print(raxs, caxs) - fft_rev = 14 #nFFT to plot - hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot - hei_rev = hei_rev[0] - #print(hei_rev) - - #print numpy.absolute(rfunc[:,0,0,14]) - - gauss_fit, covariance = None, None - for ih in range(self.minAltInd,self.maxAltInd): - for ifreq in range(self.nFFTPoints): - ''' - ###ONLY FOR TEST: - if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY - fig, axs = plt.subplots(raxs, caxs) - fig2, axs2 = plt.subplots(raxs, caxs) - col_ax = 0 - row_ax = 0 - ''' - #print(self.nPairs) - for ii in range(self.nChan): #PARES DE CANALES SELF y CROSS - # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS - # continue - # if not self.crosspairs[ii][0] in channelsProcssd: - # channelA_ok = True - #print("pair: ",self.crosspairs[ii]) - ''' - ###ONLY FOR TEST: - if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1): - 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] - - 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) - ''' - ###ONLY FOR TEST: - if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY - 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("CH "+str(self.channels[ii])) - ''' - except: - mode = mean - stdv = sigma - #print("FIT FAIL") - #continue - - - #print(mode,stdv) - #Removing echoes greater than mode + std_factor*stdv - noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero() - #noval tiene los indices que se van a remover - #print("Chan ",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(" ",self.pairsArray[ii]) - #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") - try: - pair = (self.channels[ii],self.channels[ii + 1]) - except: - pair = (99,99) - #print("par ", pair) - if ( pair in self.crosspairs): - q = self.crosspairs.index(pair) - #print("está aqui: ", q, (ii,ii + 1)) - new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0]) - cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra - - #if channelA_ok: - #chA = self.channels.index(cross_pairs[0]) - new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0]) - spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A - #channelA_ok = False - - # chB = self.channels.index(cross_pairs[1]) - # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0]) - # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B - # - # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A - # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B - ''' - ###ONLY FOR TEST: - if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY - func2clean = 10*numpy.log10(numpy.absolute(spectra[:,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("CH "+str(self.channels[ii])) - ''' - ''' - ###ONLY FOR TEST: - col_ax += 1 #contador de ploteo columnas - ##print(col_ax) - ###ONLY FOR TEST: - if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY - 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() - ''' - ################################################################################################## - - #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 len(valid[0]) >0 : - out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0]) - - 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]) - - 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 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 +# +# +# 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("runing cleanRayleigh") +# if not self.isConfig : +# +# self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv) +# +# tini=dataOut.utctime +# +# if self.byProfiles: +# if self.__profIndex == self.nIntProfiles: +# self.__dataReady = True +# else: +# if (tini - self.__initime) >= self.__integrationtime: +# +# self.__dataReady = True +# self.__initime = tini +# +# #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0): +# +# if self.__dataReady: +# +# 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)) +# try: +# jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights)) +# except: +# print("no cspc") +# 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) +# try: +# self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0) +# self.buffer3 += dataOut.data_dc +# except: +# pass +# 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 +# +# +# #index = tini.tm_hour*12+tini.tm_min/5 +# ''' +# #REVISAR +# ''' +# # jspc = jspc/self.nFFTPoints/self.normFactor +# # jcspc = jcspc/self.nFFTPoints/self.normFactor +# +# +# +# tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv) +# dataOut.data_spc = tmp_spectra +# dataOut.data_cspc = tmp_cspectra +# +# #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv) +# +# dataOut.data_dc = self.buffer3 +# dataOut.nIncohInt *= self.nIntProfiles +# dataOut.max_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): +# #channelsProcssd = [] +# #channelA_ok = False +# #rfunc = cspectra.copy() #self.bloques +# rfunc = spectra.copy() +# #rfunc = cspectra +# #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 +# +# +# ###ONLY FOR TEST: +# raxs = math.ceil(math.sqrt(self.nPairs)) +# if raxs == 0: +# raxs = 1 +# caxs = math.ceil(self.nPairs/raxs) +# if self.nPairs <4: +# raxs = 2 +# caxs = 2 +# #print(raxs, caxs) +# fft_rev = 14 #nFFT to plot +# hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot +# hei_rev = hei_rev[0] +# #print(hei_rev) +# +# #print numpy.absolute(rfunc[:,0,0,14]) +# +# gauss_fit, covariance = None, None +# for ih in range(self.minAltInd,self.maxAltInd): +# for ifreq in range(self.nFFTPoints): +# ''' +# ###ONLY FOR TEST: +# if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY +# fig, axs = plt.subplots(raxs, caxs) +# fig2, axs2 = plt.subplots(raxs, caxs) +# col_ax = 0 +# row_ax = 0 +# ''' +# #print(self.nPairs) +# for ii in range(self.nChan): #PARES DE CANALES SELF y CROSS +# # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS +# # continue +# # if not self.crosspairs[ii][0] in channelsProcssd: +# # channelA_ok = True +# #print("pair: ",self.crosspairs[ii]) +# ''' +# ###ONLY FOR TEST: +# if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1): +# 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] +# +# 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) +# ''' +# ###ONLY FOR TEST: +# if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY +# 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("CH "+str(self.channels[ii])) +# ''' +# except: +# mode = mean +# stdv = sigma +# #print("FIT FAIL") +# #continue +# +# +# #print(mode,stdv) +# #Removing echoes greater than mode + std_factor*stdv +# noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero() +# #noval tiene los indices que se van a remover +# #print("Chan ",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(" ",self.pairsArray[ii]) +# #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") +# try: +# pair = (self.channels[ii],self.channels[ii + 1]) +# except: +# pair = (99,99) +# #print("par ", pair) +# if ( pair in self.crosspairs): +# q = self.crosspairs.index(pair) +# #print("está aqui: ", q, (ii,ii + 1)) +# new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0]) +# cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra +# +# #if channelA_ok: +# #chA = self.channels.index(cross_pairs[0]) +# new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0]) +# spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A +# #channelA_ok = False +# +# # chB = self.channels.index(cross_pairs[1]) +# # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0]) +# # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B +# # +# # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A +# # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B +# ''' +# ###ONLY FOR TEST: +# if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY +# func2clean = 10*numpy.log10(numpy.absolute(spectra[:,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("CH "+str(self.channels[ii])) +# ''' +# ''' +# ###ONLY FOR TEST: +# col_ax += 1 #contador de ploteo columnas +# ##print(col_ax) +# ###ONLY FOR TEST: +# if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY +# 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() +# ''' +# ################################################################################################## +# +# #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 len(valid[0]) >0 : +# out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0]) +# +# 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]) +# +# 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 IntegrationFaradaySpectra(Operation): @@ -1313,7 +1314,7 @@ class IntegrationFaradaySpectra(Operation): self.navg = avg #self.ByLags = dataOut.ByLags ###REDEFINIR self.ByLags = False - self.maxProfilesInt = 1 + self.maxProfilesInt = 0 self.__nChannels = dataOut.nChannels if DPL != None: self.DPL=DPL @@ -1542,7 +1543,7 @@ class IntegrationFaradaySpectra(Operation): data_cspc = None data_dc = self.__buffer_dc #(CH, HEIGH) - self.maxProfilesInt = self.__profIndex + self.maxProfilesInt = self.__profIndex - 1 n = self.__profIndex - self.dataOutliers # n becomes a matrix self.__buffer_spc = [] @@ -1656,8 +1657,8 @@ class IntegrationFaradaySpectra(Operation): self.dataOut.nIncohInt *= self.n_ints - self.dataOut.max_nIncohInt = self.maxProfilesInt - #print(self.dataOut.max_nIncohInt) + #print("maxProfilesInt: ",self.maxProfilesInt) + self.dataOut.utctime = avgdatatime self.dataOut.flagNoData = False #print("Faraday Integration DONE...", self.dataOut.data_cspc) @@ -2130,7 +2131,6 @@ class IncohInt(Operation): dataOut.data_outlier = self.nOutliers dataOut.utctime = avgdatatime dataOut.flagNoData = False - dataOut.max_nIncohInt += self.__profIndex self.incohInt = 0 self.nOutliers = 0 self.__profIndex = 0 diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index fb99558..fa8d8a2 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -3250,7 +3250,7 @@ class RemoveProfileSats2(Operation): fnoise = scipy.signal.filtfilt(b, a, noise_ref) #noise_refdB = 10* numpy.log10(noise_ref) #print("Noise ",numpy.percentile(noise_ref,95)) - p85 = numpy.percentile(fnoise,85) + p85 = numpy.percentile(fnoise,83) mean_noise = fnoise.mean() if self.prev_pnoise != None: if mean_noise < (1.5 * self.prev_pnoise) : @@ -3280,7 +3280,7 @@ class RemoveProfileSats2(Operation): if index[0].size < int(self.navg*profiles): indexes = numpy.append(indexes, index[0]) - #print(indexes) + # from matplotlib import pyplot as plt # fig, ax = plt.subplots() # ax.plot(fpower) @@ -3288,7 +3288,8 @@ class RemoveProfileSats2(Operation): # ax.axhline(std, color='b') # plt.grid() # plt.show() - # print(indexes) + + #print(indexes) #outliers_IDs = outliers_IDs.astype(numpy.dtype('int64')) #outliers_IDs = numpy.unique(outliers_IDs) @@ -3301,15 +3302,19 @@ class RemoveProfileSats2(Operation): hist, bins = numpy.histogram(outs_lines,bins=my_bins) - hist_outliers_indexes = numpy.where(hist > self.thHistOutlier) #es outlier + hist_outliers_indexes = numpy.where(hist >= self.thHistOutlier) #es outlier hist_outliers_indexes = hist_outliers_indexes[0] if len(hist_outliers_indexes>0): hist_outliers_indexes = numpy.append(hist_outliers_indexes,hist_outliers_indexes[-1]+1) bins_outliers_indexes = [int(i) for i in (bins[hist_outliers_indexes])] # outlier_loc_index = [] - - outlier_loc_index = [e for n in range(len(bins_outliers_indexes)) for e in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n] + self.profileMargin) ] + #print("out indexes ", bins_outliers_indexes) + if len(bins_outliers_indexes) <= 3: + extprof = 0 + else: + extprof = self.profileMargin + outlier_loc_index = [e for n in range(len(bins_outliers_indexes)) for e in range(bins_outliers_indexes[n]-extprof,bins_outliers_indexes[n] + extprof) ] outlier_loc_index = numpy.asarray(outlier_loc_index) # if len(outlier_loc_index)>1: # ipmax = numpy.where(fpower==fpower.max())[0] @@ -3331,7 +3336,7 @@ class RemoveProfileSats2(Operation): # o = numpy.nanstd(dat) # #print(m, o, x.shape, y.shape) # #c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o)) - # c = ax[i][0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = 55, vmax = 65) + # c = ax[i][0].pcolormesh(x, y, dat.T, cmap ='jet', vmin = 60, vmax = 70) # ax[i][0].vlines(outs_lines,650,700, linestyles='dashed', label = 'outs', color='w') # #fig.colorbar(c) # ax[i][0].vlines(outlier_loc_index,700,750, linestyles='dashed', label = 'outs', color='r')