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