##// END OF EJS Templates
schain
RSS
Operando clean Rayleigh con crossSpectra
joabAM -
r1392:e987b2f0c1da
parent child
Show More
Side by Side Unified
Context file:
r1392:e987b2f0c1da Loading diff...:
Show file before | Show file after | Hide comments
@@ -116,6 +116,8 class CrossSpectraPlot(Plot):
116 zmax_coh = None
116 zmax_coh = None
117 zmin_phase = None
117 zmin_phase = None
118 zmax_phase = None
118 zmax_phase = None
119 realChannels = None
120 crossPairs = None
119
121
120 def setup(self):
122 def setup(self):
121
123
@@ -136,11 +138,17 class CrossSpectraPlot(Plot):
136 spc = dataOut.data_spc
138 spc = dataOut.data_spc
137 cspc = dataOut.data_cspc
139 cspc = dataOut.data_cspc
138 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
140 meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
139 meta['pairs'] = dataOut.pairsList
141 rawPairs = list(combinations(list(range(dataOut.nChannels)), 2))
142 #print(rawPairs)
143 meta['pairs'] = rawPairs
144
145 if self.crossPairs == None:
146 self.crossPairs = dataOut.pairsList
140
147
141 tmp = []
148 tmp = []
142
149
143 for n, pair in enumerate(meta['pairs']):
150 for n, pair in enumerate(meta['pairs']):
151
144 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
152 out = cspc[n] / numpy.sqrt(spc[pair[0]] * spc[pair[1]])
145 coh = numpy.abs(out)
153 coh = numpy.abs(out)
146 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
154 phase = numpy.arctan2(out.imag, out.real) * 180 / numpy.pi
@@ -170,12 +178,15 class CrossSpectraPlot(Plot):
170
178
171 data = self.data[-1]
179 data = self.data[-1]
172 cspc = data['cspc']
180 cspc = data['cspc']
173
181 #print(self.crossPairs)
174 for n in range(len(self.data.pairs)):
182 for n in range(len(self.data.pairs)):
175 pair = self.data.pairs[n]
183 #pair = self.data.pairs[n]
184 pair = self.crossPairs[n]
185
176 coh = cspc[n*2]
186 coh = cspc[n*2]
177 phase = cspc[n*2+1]
187 phase = cspc[n*2+1]
178 ax = self.axes[2 * n]
188 ax = self.axes[2 * n]
189
179 if ax.firsttime:
190 if ax.firsttime:
180 ax.plt = ax.pcolormesh(x, y, coh.T,
191 ax.plt = ax.pcolormesh(x, y, coh.T,
181 vmin=0,
192 vmin=0,
@@ -238,10 +249,14 class RTIPlot(Plot):
238 self.y = self.data.yrange
249 self.y = self.data.yrange
239 self.z = self.data[self.CODE]
250 self.z = self.data[self.CODE]
240 self.z = numpy.ma.masked_invalid(self.z)
251 self.z = numpy.ma.masked_invalid(self.z)
241 if self.channelList != None:
252 try:
242 self.titles = ['{} Channel {}'.format(
253 if self.channelList != None:
243 self.CODE.upper(), x) for x in self.channelList]
254 self.titles = ['{} Channel {}'.format(
244
255 self.CODE.upper(), x) for x in self.channelList]
256 except:
257 if self.channelList.any() != None:
258 self.titles = ['{} Channel {}'.format(
259 self.CODE.upper(), x) for x in self.channelList]
245 if self.decimation is None:
260 if self.decimation is None:
246 x, y, z = self.fill_gaps(self.x, self.y, self.z)
261 x, y, z = self.fill_gaps(self.x, self.y, self.z)
247 else:
262 else:
Show file before | Show file after | Hide comments
@@ -487,6 +487,8 def fit_func( x, a0, a1, a2): #, a3, a4, a5):
487 z = (x - a1) / a2
487 z = (x - a1) / a2
488 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
488 y = a0 * numpy.exp(-z**2 / a2) #+ a3 + a4 * x + a5 * x**2
489 return y
489 return y
490
491
490 class CleanRayleigh(Operation):
492 class CleanRayleigh(Operation):
491
493
492 def __init__(self):
494 def __init__(self):
@@ -641,9 +643,10 class CleanRayleigh(Operation):
641
643
642 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
644 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
643 #print("OP cleanRayleigh")
645 #print("OP cleanRayleigh")
644 #import matplotlib.pyplot as plt
646 import matplotlib.pyplot as plt
645 #for k in range(149):
647 #for k in range(149):
646
648 channelsProcssd = []
649 channelA_ok = False
647 rfunc = cspectra.copy() #self.bloques
650 rfunc = cspectra.copy() #self.bloques
648 #rfunc = cspectra
651 #rfunc = cspectra
649 #val_spc = spectra*0.0 #self.bloque0*0.0
652 #val_spc = spectra*0.0 #self.bloque0*0.0
@@ -651,24 +654,40 class CleanRayleigh(Operation):
651 #in_sat_spectra = spectra.copy() #self.bloque0
654 #in_sat_spectra = spectra.copy() #self.bloque0
652 #in_sat_cspectra = cspectra.copy() #self.bloques
655 #in_sat_cspectra = cspectra.copy() #self.bloques
653
656
654 #raxs = math.ceil(math.sqrt(self.nPairs))
655 #caxs = math.ceil(self.nPairs/raxs)
656
657
657 #print(self.hval)
658 ###ONLY FOR TEST:
659 raxs = math.ceil(math.sqrt(self.nPairs))
660 caxs = math.ceil(self.nPairs/raxs)
661 if self.nPairs <4:
662 raxs = 2
663 caxs = 2
664 #print(raxs, caxs)
665 fft_rev = 14 #nFFT to plot
666 hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot
667 hei_rev = hei_rev[0]
668 #print(hei_rev)
669
658 #print numpy.absolute(rfunc[:,0,0,14])
670 #print numpy.absolute(rfunc[:,0,0,14])
671
659 gauss_fit, covariance = None, None
672 gauss_fit, covariance = None, None
660 for ih in range(self.minAltInd,self.maxAltInd):
673 for ih in range(self.minAltInd,self.maxAltInd):
661 for ifreq in range(self.nFFTPoints):
674 for ifreq in range(self.nFFTPoints):
662 # fig, axs = plt.subplots(raxs, caxs)
675 ###ONLY FOR TEST:
663 # fig2, axs2 = plt.subplots(raxs, caxs)
676 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
664 # col_ax = 0
677 fig, axs = plt.subplots(raxs, caxs)
665 # row_ax = 0
678 fig2, axs2 = plt.subplots(raxs, caxs)
666 #print(len(self.nPairs))
679 col_ax = 0
667 for ii in range(self.nPairs): #PARES DE CANALES SELF y CROSS
680 row_ax = 0
668 #print("ii: ",ii)
681 #print(self.nPairs)
669 # if (col_ax%caxs==0 and col_ax!=0):
682 for ii in range(self.nPairs): #PARES DE CANALES SELF y CROSS
670 # col_ax = 0
683 if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
671 # row_ax += 1
684 continue
685 if not self.crosspairs[ii][0] in channelsProcssd:
686 channelA_ok = True
687 #print("pair: ",self.crosspairs[ii])
688 if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1): ###ONLY FOR TEST:
689 col_ax = 0
690 row_ax += 1
672 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
691 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih])) #Potencia?
673 #print(func2clean.shape)
692 #print(func2clean.shape)
674 val = (numpy.isfinite(func2clean)==True).nonzero()
693 val = (numpy.isfinite(func2clean)==True).nonzero()
@@ -692,22 +711,26 class CleanRayleigh(Operation):
692 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
711 sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
693 parg = [numpy.amax(y_dist),mean,sigma]
712 parg = [numpy.amax(y_dist),mean,sigma]
694
713
695 #newY = None
714 newY = None
696
715
697 try :
716 try :
698 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
717 gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
699 mode = gauss_fit[1]
718 mode = gauss_fit[1]
700 stdv = gauss_fit[2]
719 stdv = gauss_fit[2]
701 #print(" FIT OK",gauss_fit)
720 #print(" FIT OK",gauss_fit)
702 '''
721
703 newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
722 ###ONLY FOR TEST:
704 axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
723 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
705 axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
724 newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
706 axs[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))'''
725 axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
726 axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
727 axs[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))
728
707 except:
729 except:
708 mode = mean
730 mode = mean
709 stdv = sigma
731 stdv = sigma
710 #print("FIT FAIL")
732 #print("FIT FAIL")
733 continue
711
734
712
735
713 #print(mode,stdv)
736 #print(mode,stdv)
@@ -727,36 +750,45 class CleanRayleigh(Operation):
727 # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
750 # val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
728 # val_cspc[novall[0],ii,ifreq,ih] = 1
751 # val_cspc[novall[0],ii,ifreq,ih] = 1
729 #print("OUT NOVALL 1")
752 #print("OUT NOVALL 1")
730 #Removing coherent from ISR data
731 chA = self.channels.index(cross_pairs[0])
732 chB = self.channels.index(cross_pairs[1])
733
753
734 new_a = numpy.delete(cspectra[:,ii,ifreq,ih], noval[0])
754 new_a = numpy.delete(cspectra[:,ii,ifreq,ih], noval[0])
735 cspectra[noval,ii,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
755 cspectra[noval,ii,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
736 new_b = numpy.delete(spectra[:,chA,ifreq,ih], noval[0])
756
737 spectra[noval,chA,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
757 if channelA_ok:
758 chA = self.channels.index(cross_pairs[0])
759 new_b = numpy.delete(spectra[:,chA,ifreq,ih], noval[0])
760 spectra[noval,chA,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
761 channelA_ok = False
762 chB = self.channels.index(cross_pairs[1])
738 new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
763 new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
739 spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
764 spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
740
765
766 channelsProcssd.append(self.crosspairs[ii][0]) # save channel A
767 channelsProcssd.append(self.crosspairs[ii][1]) # save channel B
741
768
742 '''
769 ###ONLY FOR TEST:
743 func2clean = 10*numpy.log10(numpy.absolute(cspectra[:,ii,ifreq,ih]))
770 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
744 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
771 func2clean = 10*numpy.log10(numpy.absolute(cspectra[:,ii,ifreq,ih]))
745 axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
772 y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
746 axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
773 axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
747 axs2[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))
774 axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
748 '''
775 axs2[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))
749
776
750 #col_ax += 1 #contador de ploteo columnas
777 ###ONLY FOR TEST:
778 col_ax += 1 #contador de ploteo columnas
751 ##print(col_ax)
779 ##print(col_ax)
752 '''
780 ###ONLY FOR TEST:
753 title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
781 if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
754 title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
782 title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km"
755 fig.suptitle(title)
783 title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+" Alt: "+str(self.heights[ih])+ " km CLEANED"
756 fig2.suptitle(title2)
784 fig.suptitle(title)
757 plt.show()'''
785 fig2.suptitle(title2)
758
786 plt.show()
759 ''' channels = channels
787
788
789 '''
790
791 channels = channels
760 cross_pairs = cross_pairs
792 cross_pairs = cross_pairs
761 #print("OUT NOVALL 2")
793 #print("OUT NOVALL 2")
762
794
@@ -786,14 +818,10 class CleanRayleigh(Operation):
786 for ich in range(self.nChan):
818 for ich in range(self.nChan):
787 tmp = spectra[:,ich,ifreq,ih]
819 tmp = spectra[:,ich,ifreq,ih]
788 valid = (numpy.isfinite(tmp[:])==True).nonzero()
820 valid = (numpy.isfinite(tmp[:])==True).nonzero()
789 # if ich == 0 and ifreq == 0 and ih == 17 :
821
790 # print tmp
791 # print valid
792 # print len(valid[0])
793 #print('TMP',tmp)
794 if len(valid[0]) >0 :
822 if len(valid[0]) >0 :
795 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
823 out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
796 #for icr in range(nPairs):
824
797 for icr in range(self.nPairs):
825 for icr in range(self.nPairs):
798 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
826 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
799 valid = (numpy.isfinite(tmp)==True).nonzero()
827 valid = (numpy.isfinite(tmp)==True).nonzero()
General Comments 0
You need to be logged in to leave comments. Login now