From e987b2f0c1dad7db3362b9e9ce71ecdb696b828f 2021-10-27 13:34:46 From: joabAM Date: 2021-10-27 13:34:46 Subject: [PATCH] Operando clean Rayleigh con crossSpectra --- diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index 17b9192..f975e1e 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -116,6 +116,8 @@ class CrossSpectraPlot(Plot): zmax_coh = None zmin_phase = None zmax_phase = None + realChannels = None + crossPairs = None def setup(self): @@ -136,11 +138,17 @@ class CrossSpectraPlot(Plot): spc = dataOut.data_spc cspc = dataOut.data_cspc meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1)) - meta['pairs'] = dataOut.pairsList + rawPairs = list(combinations(list(range(dataOut.nChannels)), 2)) + #print(rawPairs) + meta['pairs'] = rawPairs + + if self.crossPairs == None: + self.crossPairs = dataOut.pairsList tmp = [] for n, pair in enumerate(meta['pairs']): + out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]]) coh = numpy.abs(out) phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi @@ -170,12 +178,15 @@ class CrossSpectraPlot(Plot): data = self.data[-1] cspc = data['cspc'] - + #print(self.crossPairs) for n in range(len(self.data.pairs)): - pair = self.data.pairs[n] + #pair = self.data.pairs[n] + pair = self.crossPairs[n] + coh = cspc[n*2] phase = cspc[n*2+1] ax = self.axes[2 * n] + if ax.firsttime: ax.plt = ax.pcolormesh(x, y, coh.T, vmin=0, @@ -238,10 +249,14 @@ class RTIPlot(Plot): self.y = self.data.yrange self.z = self.data[self.CODE] self.z = numpy.ma.masked_invalid(self.z) - if self.channelList != None: - self.titles = ['{} Channel {}'.format( - self.CODE.upper(), x) for x in self.channelList] - + try: + if self.channelList != None: + self.titles = ['{} Channel {}'.format( + self.CODE.upper(), x) for x in self.channelList] + except: + if self.channelList.any() != None: + self.titles = ['{} Channel {}'.format( + self.CODE.upper(), x) for x in self.channelList] if self.decimation is None: x, y, z = self.fill_gaps(self.x, self.y, self.z) else: diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index dcb4528..22f6622 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -487,6 +487,8 @@ 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): @@ -641,9 +643,10 @@ class CleanRayleigh(Operation): def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv): #print("OP cleanRayleigh") - #import matplotlib.pyplot as plt + import matplotlib.pyplot as plt #for k in range(149): - + channelsProcssd = [] + channelA_ok = False rfunc = cspectra.copy() #self.bloques #rfunc = cspectra #val_spc = spectra*0.0 #self.bloque0*0.0 @@ -651,24 +654,40 @@ class CleanRayleigh(Operation): #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) + ###ONLY FOR TEST: + raxs = math.ceil(math.sqrt(self.nPairs)) + 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): - # fig, axs = plt.subplots(raxs, caxs) - # fig2, axs2 = plt.subplots(raxs, caxs) - # col_ax = 0 - # row_ax = 0 - #print(len(self.nPairs)) - 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 + ###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.nPairs): #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]) + if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1): ###ONLY FOR TEST: + 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() @@ -692,22 +711,26 @@ class CleanRayleigh(Operation): sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist)) parg = [numpy.amax(y_dist),mean,sigma] - #newY = 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]))''' + + ###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("Pair "+str(self.crosspairs[ii])) + except: mode = mean stdv = sigma #print("FIT FAIL") + continue #print(mode,stdv) @@ -727,36 +750,45 @@ class CleanRayleigh(Operation): # 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 - chA = self.channels.index(cross_pairs[0]) - chB = self.channels.index(cross_pairs[1]) new_a = numpy.delete(cspectra[:,ii,ifreq,ih], noval[0]) cspectra[noval,ii,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra - new_b = numpy.delete(spectra[:,chA,ifreq,ih], noval[0]) - spectra[noval,chA,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A + + if channelA_ok: + chA = self.channels.index(cross_pairs[0]) + new_b = numpy.delete(spectra[:,chA,ifreq,ih], noval[0]) + spectra[noval,chA,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 - ''' - 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])) - ''' + ###ONLY FOR TEST: + if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY + 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 + ###ONLY FOR TEST: + 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 + ###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() + + + ''' + + channels = channels cross_pairs = cross_pairs #print("OUT NOVALL 2") @@ -786,14 +818,10 @@ class CleanRayleigh(Operation): 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()