schain Commit - r1392:e987b2f0c1da · Jicamarca Repository

Operando clean Rayleigh con crossSpectra

joabAM -

r1392:e987b2f0c1da

parent child

Context file:

r1392:e987b2f0c1da

Collapse all files

schainpy/model/graphics/jroplot_spectra.py +22 -7

             # Copyright (c) 2012-2020 Jicamarca Radio Observatory
             # All rights reserved.
             #
             # Distributed under the terms of the BSD 3-clause license.
             """Classes to plot Spectra data
             """
             import os
             import numpy
             from schainpy.model.graphics.jroplot_base import Plot, plt, log
             class SpectraPlot(Plot):
                 '''
                 Plot for Spectra data
                 '''
                 CODE = 'spc'
                 colormap = 'jet'
                 plot_type = 'pcolor'
                 buffering = False
                 channelList = []
                 def setup(self):
                     self.nplots = len(self.data.channels)
                     self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
                     self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
                     self.height = 2.6 * self.nrows
                     self.cb_label = 'dB'
                     if self.showprofile:
                         self.width = 4 * self.ncols
                     else:
                         self.width = 3.5 * self.ncols
                     self.plots_adjust.update({'wspace': 0.4, 'hspace':0.4, 'left': 0.1, 'right': 0.9, 'bottom': 0.08})
                     self.ylabel = 'Range [km]'
                 def update(self, dataOut):
                     if self.channelList == None:
                         self.channelList = dataOut.channelList
                     data = {}
                     meta = {}
                     spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
                     data['spc'] = spc
                     data['rti'] = dataOut.getPower()
                     data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
                     meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
                     if self.CODE == 'spc_moments':
                         data['moments'] = dataOut.moments
                     return data, meta
                 def plot(self):
                     if self.xaxis == "frequency":
                         x = self.data.xrange[0]
                         self.xlabel = "Frequency (kHz)"
                     elif self.xaxis == "time":
                         x = self.data.xrange[1]
                         self.xlabel = "Time (ms)"
                     else:
                         x = self.data.xrange[2]
                         self.xlabel = "Velocity (m/s)"
                     if self.CODE == 'spc_moments':
                         x = self.data.xrange[2]
                         self.xlabel = "Velocity (m/s)"
                     self.titles = []
                     y = self.data.yrange
                     self.y = y
                     data = self.data[-1]
                     z = data['spc']
                     for n, ax in enumerate(self.axes):
                         noise = data['noise'][n]
                         if self.CODE == 'spc_moments':
                             mean = data['moments'][n, 1]
                         if ax.firsttime:
                             self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
                             self.xmin = self.xmin if self.xmin else -self.xmax
                             self.zmin = self.zmin if self.zmin else numpy.nanmin(z)
                             self.zmax = self.zmax if self.zmax else numpy.nanmax(z)
                             ax.plt = ax.pcolormesh(x, y, z[n].T,
                                                    vmin=self.zmin,
                                                    vmax=self.zmax,
                                                    cmap=plt.get_cmap(self.colormap)
                                                    )
                             if self.showprofile:
                                 ax.plt_profile = self.pf_axes[n].plot(
                                     data['rti'][n], y)[0]
                                 ax.plt_noise = self.pf_axes[n].plot(numpy.repeat(noise, len(y)), y,
                                                                     color="k", linestyle="dashed", lw=1)[0]
                             if self.CODE == 'spc_moments':
                                 ax.plt_mean = ax.plot(mean, y, color='k')[0]
                         else:
                             ax.plt.set_array(z[n].T.ravel())
                             if self.showprofile:
                                 ax.plt_profile.set_data(data['rti'][n], y)
                                 ax.plt_noise.set_data(numpy.repeat(noise, len(y)), y)
                             if self.CODE == 'spc_moments':
                                 ax.plt_mean.set_data(mean, y)
                         self.titles.append('CH {}: {:3.2f}dB'.format(self.channelList[n], noise))
             class CrossSpectraPlot(Plot):
                 CODE = 'cspc'
                 colormap = 'jet'
                 plot_type = 'pcolor'
                 zmin_coh = None
                 zmax_coh = None
                 zmin_phase = None
                 zmax_phase = None
+                realChannels = None
+                crossPairs = None
                 def setup(self):
                     self.ncols = 4
                     self.nplots = len(self.data.pairs) * 2
                     self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
                     self.width = 3.1 * self.ncols
                     self.height = 2.6 * self.nrows
                     self.ylabel = 'Range [km]'
                     self.showprofile = False
                     self.plots_adjust.update({'left': 0.08, 'right': 0.92, 'wspace': 0.5, 'hspace':0.4, 'top':0.95, 'bottom': 0.08})
                 def update(self, dataOut):
                     data = {}
                     meta = {}
                     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
                         tmp.append(coh)
                         tmp.append(phase)
                     data['cspc'] = numpy.array(tmp)
                     return data, meta
                 def plot(self):
                     if self.xaxis == "frequency":
                         x = self.data.xrange[0]
                         self.xlabel = "Frequency (kHz)"
                     elif self.xaxis == "time":
                         x = self.data.xrange[1]
                         self.xlabel = "Time (ms)"
                     else:
                         x = self.data.xrange[2]
                         self.xlabel = "Velocity (m/s)"
                     self.titles = []
                     y = self.data.yrange
                     self.y = y
                     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,
                                                    vmax=1,
                                                    cmap=plt.get_cmap(self.colormap_coh)
                                                    )
                         else:
                             ax.plt.set_array(coh.T.ravel())
                         self.titles.append(
                             'Coherence Ch{} * Ch{}'.format(pair[0], pair[1]))
                         ax = self.axes[2 * n + 1]
                         if ax.firsttime:
                             ax.plt = ax.pcolormesh(x, y, phase.T,
                                                    vmin=-180,
                                                    vmax=180,
                                                    cmap=plt.get_cmap(self.colormap_phase)
                                                    )
                         else:
                             ax.plt.set_array(phase.T.ravel())
                         self.titles.append('Phase CH{} * CH{}'.format(pair[0], pair[1]))
             class RTIPlot(Plot):
                 '''
                 Plot for RTI data
                 '''
                 CODE = 'rti'
                 colormap = 'jet'
                 plot_type = 'pcolorbuffer'
                 titles = None
                 channelList = []
                 def setup(self):
                     self.xaxis = 'time'
                     self.ncols = 1
                     print("dataChannels ",self.data.channels)
                     self.nrows = len(self.data.channels)
                     self.nplots = len(self.data.channels)
                     self.ylabel = 'Range [km]'
                     self.xlabel = 'Time'
                     self.cb_label = 'dB'
                     self.plots_adjust.update({'hspace':0.8, 'left': 0.1, 'bottom': 0.08, 'right':0.95})
                     self.titles = ['{} Channel {}'.format(
                         self.CODE.upper(), x) for x in range(self.nplots)]
                     print("SETUP")
                 def update(self, dataOut):
                     if len(self.channelList) == 0:
                         self.channelList = dataOut.channelList
                     data = {}
                     meta = {}
                     data['rti'] = dataOut.getPower()
                     data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
                     return data, meta
                 def plot(self):
                     self.x = self.data.times
                     self.y = self.data.yrange
                     self.z = self.data[self.CODE]
                     self.z = numpy.ma.masked_invalid(self.z)
-                    if self.channelList != None:
+                    try:
-                        self.titles = ['{} Channel {}'.format(
+                        if self.channelList != None:
-                            self.CODE.upper(), x) for x in self.channelList]
+                            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:
                         x, y, z = self.fill_gaps(*self.decimate())
                     for n, ax in enumerate(self.axes):
                         self.zmin = self.zmin if self.zmin else numpy.min(self.z)
                         self.zmax = self.zmax if self.zmax else numpy.max(self.z)
                         data = self.data[-1]
                         if ax.firsttime:
                             ax.plt = ax.pcolormesh(x, y, z[n].T,
                                                    vmin=self.zmin,
                                                    vmax=self.zmax,
                                                    cmap=plt.get_cmap(self.colormap)
                                                    )
                             if self.showprofile:
                                 ax.plot_profile = self.pf_axes[n].plot(
                                     data['rti'][n], self.y)[0]
                                 ax.plot_noise = self.pf_axes[n].plot(numpy.repeat(data['noise'][n], len(self.y)), self.y,
                                                                      color="k", linestyle="dashed", lw=1)[0]
                         else:
                             ax.collections.remove(ax.collections[0])
                             ax.plt = ax.pcolormesh(x, y, z[n].T,
                                                    vmin=self.zmin,
                                                    vmax=self.zmax,
                                                    cmap=plt.get_cmap(self.colormap)
                                                    )
                             if self.showprofile:
                                 ax.plot_profile.set_data(data['rti'][n], self.y)
                                 ax.plot_noise.set_data(numpy.repeat(
                                     data['noise'][n], len(self.y)), self.y)
             class CoherencePlot(RTIPlot):
                 '''
                 Plot for Coherence data
                 '''
                 CODE = 'coh'
                 def setup(self):
                     self.xaxis = 'time'
                     self.ncols = 1
                     self.nrows = len(self.data.pairs)
                     self.nplots = len(self.data.pairs)
                     self.ylabel = 'Range [km]'
                     self.xlabel = 'Time'
                     self.plots_adjust.update({'hspace':0.6, 'left': 0.1, 'bottom': 0.1,'right':0.95})
                     if self.CODE == 'coh':
                         self.cb_label = ''
                         self.titles = [
                             'Coherence Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
                     else:
                         self.cb_label = 'Degrees'
                         self.titles = [
                             'Phase Map Ch{} * Ch{}'.format(x[0], x[1]) for x in self.data.pairs]
                 def update(self, dataOut):
                     data = {}
                     meta = {}
                     data['coh'] = dataOut.getCoherence()
                     meta['pairs'] = dataOut.pairsList
                     return data, meta
             class PhasePlot(CoherencePlot):
                 '''
                 Plot for Phase map data
                 '''
                 CODE = 'phase'
                 colormap = 'seismic'
                 def update(self, dataOut):
                     data = {}
                     meta = {}
                     data['phase'] = dataOut.getCoherence(phase=True)
                     meta['pairs'] = dataOut.pairsList
                     return data, meta
             class NoisePlot(Plot):
                 '''
                 Plot for noise
                 '''
                 CODE = 'noise'
                 plot_type = 'scatterbuffer'
                 def setup(self):
                     self.xaxis = 'time'
                     self.ncols = 1
                     self.nrows = 1
                     self.nplots = 1
                     self.ylabel = 'Intensity [dB]'
                     self.xlabel = 'Time'
                     self.titles = ['Noise']
                     self.colorbar = False
                     self.plots_adjust.update({'right': 0.85 })
                 def update(self, dataOut):
                     data = {}
                     meta = {}
                     data['noise'] = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor).reshape(dataOut.nChannels, 1)
                     meta['yrange'] = numpy.array([])
                     return data, meta
                 def plot(self):
                     x = self.data.times
                     xmin = self.data.min_time
                     xmax = xmin + self.xrange * 60 * 60
                     Y = self.data['noise']
                     if self.axes[0].firsttime:
                         self.ymin = numpy.nanmin(Y) - 5
                         self.ymax = numpy.nanmax(Y) + 5
                         for ch in self.data.channels:
                             y = Y[ch]
                             self.axes[0].plot(x, y, lw=1, label='Ch{}'.format(ch))
                         plt.legend(bbox_to_anchor=(1.18, 1.0))
                     else:
                         for ch in self.data.channels:
                             y = Y[ch]
                             self.axes[0].lines[ch].set_data(x, y)
             class PowerProfilePlot(Plot):
                 CODE = 'pow_profile'
                 plot_type = 'scatter'
                 def setup(self):
                     self.ncols = 1
                     self.nrows = 1
                     self.nplots = 1
                     self.height = 4
                     self.width = 3
                     self.ylabel = 'Range [km]'
                     self.xlabel = 'Intensity [dB]'
                     self.titles = ['Power Profile']
                     self.colorbar = False
                 def update(self, dataOut):
                     data = {}
                     meta = {}
                     data[self.CODE] = dataOut.getPower()
                     return data, meta
                 def plot(self):
                     y = self.data.yrange
                     self.y = y
                     x = self.data[-1][self.CODE]
                     if self.xmin is None: self.xmin = numpy.nanmin(x)*0.9
                     if self.xmax is None: self.xmax = numpy.nanmax(x)*1.1
                     if self.axes[0].firsttime:
                         for ch in self.data.channels:
                             self.axes[0].plot(x[ch], y, lw=1, label='Ch{}'.format(ch))
                         plt.legend()
                     else:
                         for ch in self.data.channels:
                             self.axes[0].lines[ch].set_data(x[ch], y)
             class SpectraCutPlot(Plot):
                 CODE = 'spc_cut'
                 plot_type = 'scatter'
                 buffering = False
                 def setup(self):
                     self.nplots = len(self.data.channels)
                     self.ncols = int(numpy.sqrt(self.nplots) + 0.9)
                     self.nrows = int((1.0 * self.nplots / self.ncols) + 0.9)
                     self.width = 3.4 * self.ncols + 1.5
                     self.height = 3 * self.nrows
                     self.ylabel = 'Power [dB]'
                     self.colorbar = False
                     self.plots_adjust.update({'left':0.1, 'hspace':0.3, 'right': 0.75, 'bottom':0.08})
                 def update(self, dataOut):
                     data = {}
                     meta = {}
                     spc = 10*numpy.log10(dataOut.data_spc/dataOut.normFactor)
                     data['spc'] = spc
                     meta['xrange'] = (dataOut.getFreqRange(1)/1000., dataOut.getAcfRange(1), dataOut.getVelRange(1))
                     return data, meta
                 def plot(self):
                     if self.xaxis == "frequency":
                         x = self.data.xrange[0][1:]
                         self.xlabel = "Frequency (kHz)"
                     elif self.xaxis == "time":
                         x = self.data.xrange[1]
                         self.xlabel = "Time (ms)"
                     else:
                         x = self.data.xrange[2]
                         self.xlabel = "Velocity (m/s)"
                     self.titles = []
                     y = self.data.yrange
                     z = self.data[-1]['spc']
                     if self.height_index:
                         index = numpy.array(self.height_index)
                     else:
                         index = numpy.arange(0, len(y), int((len(y))/9))
                     for n, ax in enumerate(self.axes):
                         if ax.firsttime:
                             self.xmax = self.xmax if self.xmax else numpy.nanmax(x)
                             self.xmin = self.xmin if self.xmin else -self.xmax
                             self.ymin = self.ymin if self.ymin else numpy.nanmin(z)
                             self.ymax = self.ymax if self.ymax else numpy.nanmax(z)
                             ax.plt = ax.plot(x, z[n, :, index].T)
                             labels = ['Range = {:2.1f}km'.format(y[i]) for i in index]
                             self.figures[0].legend(ax.plt, labels, loc='center right')
                         else:
                             for i, line in enumerate(ax.plt):
                                 line.set_data(x, z[n, :, index[i]])
                         self.titles.append('CH {}'.format(n))
             class BeaconPhase(Plot):
                 __isConfig = None
                 __nsubplots = None
                 PREFIX = 'beacon_phase'
                 def __init__(self):
                     Plot.__init__(self)
                     self.timerange = 24*60*60
                     self.isConfig = False
                     self.__nsubplots = 1
                     self.counter_imagwr = 0
                     self.WIDTH = 800
                     self.HEIGHT = 400
                     self.WIDTHPROF = 120
                     self.HEIGHTPROF = 0
                     self.xdata = None
                     self.ydata = None
                     self.PLOT_CODE = BEACON_CODE
                     self.FTP_WEI = None
                     self.EXP_CODE = None
                     self.SUB_EXP_CODE = None
                     self.PLOT_POS = None
                     self.filename_phase = None
                     self.figfile = None
                     self.xmin = None
                     self.xmax = None
                 def getSubplots(self):
                     ncol = 1
                     nrow = 1
                     return nrow, ncol
                 def setup(self, id, nplots, wintitle, showprofile=True, show=True):
                     self.__showprofile = showprofile
                     self.nplots = nplots
                     ncolspan = 7
                     colspan = 6
                     self.__nsubplots = 2
                     self.createFigure(id = id,
                                       wintitle = wintitle,
                                       widthplot = self.WIDTH+self.WIDTHPROF,
                                       heightplot = self.HEIGHT+self.HEIGHTPROF,
                                       show=show)
                     nrow, ncol = self.getSubplots()
                     self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1)
                 def save_phase(self, filename_phase):
                     f = open(filename_phase,'w+')
                     f.write('\n\n')
                     f.write('JICAMARCA RADIO OBSERVATORY - Beacon Phase \n')
                     f.write('DD MM YYYY  HH MM SS   pair(2,0) pair(2,1) pair(2,3) pair(2,4)\n\n' )
                     f.close()
                 def save_data(self, filename_phase, data, data_datetime):
                     f=open(filename_phase,'a')
                     timetuple_data = data_datetime.timetuple()
                     day = str(timetuple_data.tm_mday)
                     month = str(timetuple_data.tm_mon)
                     year = str(timetuple_data.tm_year)
                     hour = str(timetuple_data.tm_hour)
                     minute = str(timetuple_data.tm_min)
                     second = str(timetuple_data.tm_sec)
                     f.write(day+' '+month+' '+year+'  '+hour+' '+minute+' '+second+'   '+str(data[0])+'   '+str(data[1])+'   '+str(data[2])+'   '+str(data[3])+'\n')
                     f.close()
                 def plot(self):
                     log.warning('TODO: Not yet implemented...')
                 def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True',
                         xmin=None, xmax=None, ymin=None, ymax=None, hmin=None, hmax=None,
                         timerange=None,
                         save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1,
                         server=None, folder=None, username=None, password=None,
                         ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0):
                     if dataOut.flagNoData:
                         return dataOut
                     if not isTimeInHourRange(dataOut.datatime, xmin, xmax):
                         return
                     if pairsList == None:
                         pairsIndexList = dataOut.pairsIndexList[:10]
                     else:
                         pairsIndexList = []
                         for pair in pairsList:
                             if pair not in dataOut.pairsList:
                                 raise ValueError("Pair %s is not in dataOut.pairsList" %(pair))
                             pairsIndexList.append(dataOut.pairsList.index(pair))
                     if pairsIndexList == []:
                         return
              #         if len(pairsIndexList) > 4:
              #             pairsIndexList = pairsIndexList[0:4]
                     hmin_index = None
                     hmax_index = None
                     if hmin != None and hmax != None:
                         indexes = numpy.arange(dataOut.nHeights)
                         hmin_list = indexes[dataOut.heightList >= hmin]
                         hmax_list = indexes[dataOut.heightList <= hmax]
                         if hmin_list.any():
                             hmin_index = hmin_list[0]
                         if hmax_list.any():
                             hmax_index = hmax_list[-1]+1
                     x = dataOut.getTimeRange()
                     thisDatetime = dataOut.datatime
                     title = wintitle + " Signal Phase" # : %s" %(thisDatetime.strftime("%d-%b-%Y"))
                     xlabel = "Local Time"
                     ylabel = "Phase (degrees)"
                     update_figfile = False
                     nplots = len(pairsIndexList)
                     #phase = numpy.zeros((len(pairsIndexList),len(dataOut.beacon_heiIndexList)))
                     phase_beacon = numpy.zeros(len(pairsIndexList))
                     for i in range(nplots):
                         pair = dataOut.pairsList[pairsIndexList[i]]
                         ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i], :, hmin_index:hmax_index], axis=0)
                         powa = numpy.average(dataOut.data_spc[pair[0], :, hmin_index:hmax_index], axis=0)
                         powb = numpy.average(dataOut.data_spc[pair[1], :, hmin_index:hmax_index], axis=0)
                         avgcoherenceComplex = ccf/numpy.sqrt(powa*powb)
                         phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi
                         if dataOut.beacon_heiIndexList:
                             phase_beacon[i] = numpy.average(phase[dataOut.beacon_heiIndexList])
                         else:
                             phase_beacon[i] = numpy.average(phase)
                     if not self.isConfig:
                         nplots = len(pairsIndexList)
                         self.setup(id=id,
                                    nplots=nplots,
                                    wintitle=wintitle,
                                    showprofile=showprofile,
                                    show=show)
                         if timerange != None:
                             self.timerange = timerange
                         self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange)
                         if ymin == None: ymin = 0
                         if ymax == None: ymax = 360
                         self.FTP_WEI = ftp_wei
                         self.EXP_CODE = exp_code
                         self.SUB_EXP_CODE = sub_exp_code
                         self.PLOT_POS = plot_pos
                         self.name = thisDatetime.strftime("%Y%m%d_%H%M%S")
                         self.isConfig = True
                         self.figfile = figfile
                         self.xdata = numpy.array([])
                         self.ydata = numpy.array([])
                         update_figfile = True
                         #open file beacon phase
                         path = '%s%03d' %(self.PREFIX, self.id)
                         beacon_file = os.path.join(path,'%s.txt'%self.name)
                         self.filename_phase = os.path.join(figpath,beacon_file)
                         #self.save_phase(self.filename_phase)
                     #store data beacon phase
                     #self.save_data(self.filename_phase, phase_beacon, thisDatetime)
                     self.setWinTitle(title)
                     title = "Phase Plot %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S"))
                     legendlabels = ["Pair (%d,%d)"%(pair[0], pair[1]) for pair in dataOut.pairsList]
                     axes = self.axesList[0]
                     self.xdata = numpy.hstack((self.xdata, x[0:1]))
                     if len(self.ydata)==0:
                         self.ydata = phase_beacon.reshape(-1,1)
                     else:
                         self.ydata = numpy.hstack((self.ydata, phase_beacon.reshape(-1,1)))
                     axes.pmultilineyaxis(x=self.xdata, y=self.ydata,
                                 xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax,
                                 xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid",
                                 XAxisAsTime=True, grid='both'
                                 )
                     self.draw()
                     if dataOut.ltctime >= self.xmax:
                         self.counter_imagwr = wr_period
                         self.isConfig = False
                         update_figfile = True
                     self.save(figpath=figpath,
                               figfile=figfile,
                               save=save,
                               ftp=ftp,
                               wr_period=wr_period,
                               thisDatetime=thisDatetime,
                               update_figfile=update_figfile)
                     return dataOut

schainpy/model/proc/jroproc_spectra.py +76 -48

             # Copyright (c) 2012-2020 Jicamarca Radio Observatory
             # All rights reserved.
             #
             # Distributed under the terms of the BSD 3-clause license.
             """Spectra processing Unit and operations
             Here you will find the processing unit `SpectraProc` and several operations
             to work with Spectra data type
             """
             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):
                 def __init__(self):
                     ProcessingUnit.__init__(self)
                     self.buffer = None
                     self.firstdatatime = None
                     self.profIndex = 0
                     self.dataOut = Spectra()
                     self.id_min = None
                     self.id_max = None
                     self.setupReq = False #Agregar a todas las unidades de proc
                 def __updateSpecFromVoltage(self):
                     self.dataOut.timeZone = self.dataIn.timeZone
                     self.dataOut.dstFlag = self.dataIn.dstFlag
                     self.dataOut.errorCount = self.dataIn.errorCount
                     self.dataOut.useLocalTime = self.dataIn.useLocalTime
                     try:
                         self.dataOut.processingHeaderObj = self.dataIn.processingHeaderObj.copy()
                     except:
                         pass
                     self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
                     self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy()
                     self.dataOut.channelList = self.dataIn.channelList
                     self.dataOut.heightList = self.dataIn.heightList
                     self.dataOut.dtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
                     self.dataOut.nProfiles = self.dataOut.nFFTPoints
                     self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
                     self.dataOut.utctime = self.firstdatatime
                     self.dataOut.flagDecodeData = self.dataIn.flagDecodeData
                     self.dataOut.flagDeflipData = self.dataIn.flagDeflipData
                     self.dataOut.flagShiftFFT = False
                     self.dataOut.nCohInt = self.dataIn.nCohInt
                     self.dataOut.nIncohInt = 1
                     self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
                     self.dataOut.frequency = self.dataIn.frequency
                     self.dataOut.realtime = self.dataIn.realtime
                     self.dataOut.azimuth = self.dataIn.azimuth
                     self.dataOut.zenith = self.dataIn.zenith
                     self.dataOut.codeList = self.dataIn.codeList
                     self.dataOut.azimuthList = self.dataIn.azimuthList
                     self.dataOut.elevationList = self.dataIn.elevationList
                 def __getFft(self):
                     """
                     Convierte valores de Voltaje a Spectra
                     Affected:
                         self.dataOut.data_spc
                         self.dataOut.data_cspc
                         self.dataOut.data_dc
                         self.dataOut.heightList
                         self.profIndex
                         self.buffer
                         self.dataOut.flagNoData
                     """
                     fft_volt = numpy.fft.fft(
                         self.buffer, n=self.dataOut.nFFTPoints, axis=1)
                     fft_volt = fft_volt.astype(numpy.dtype('complex'))
                     dc = fft_volt[:, 0, :]
                     # calculo de self-spectra
                     fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,))
                     spc = fft_volt * numpy.conjugate(fft_volt)
                     spc = spc.real
                     blocksize = 0
                     blocksize += dc.size
                     blocksize += spc.size
                     cspc = None
                     pairIndex = 0
                     if self.dataOut.pairsList != None:
                         # calculo de cross-spectra
                         cspc = numpy.zeros(
                             (self.dataOut.nPairs, self.dataOut.nFFTPoints, self.dataOut.nHeights), dtype='complex')
                         for pair in self.dataOut.pairsList:
                             if pair[0] not in self.dataOut.channelList:
                                 raise ValueError("Error getting CrossSpectra: pair 0 of %s is not in channelList = %s" % (
                                     str(pair), str(self.dataOut.channelList)))
                             if pair[1] not in self.dataOut.channelList:
                                 raise ValueError("Error getting CrossSpectra: pair 1 of %s is not in channelList = %s" % (
                                     str(pair), str(self.dataOut.channelList)))
                             cspc[pairIndex, :, :] = fft_volt[pair[0], :, :] * \
                                 numpy.conjugate(fft_volt[pair[1], :, :])
                             pairIndex += 1
                         blocksize += cspc.size
                     self.dataOut.data_spc = spc
                     self.dataOut.data_cspc = cspc
                     self.dataOut.data_dc = dc
                     self.dataOut.blockSize = blocksize
                     self.dataOut.flagShiftFFT = False
                 def run(self, nProfiles=None, nFFTPoints=None, pairsList=None, ippFactor=None, shift_fft=False):
                     if self.dataIn.type == "Spectra":
                         self.dataOut.copy(self.dataIn)
                         if shift_fft:
                             #desplaza a la derecha en el eje 2 determinadas posiciones
                             shift = int(self.dataOut.nFFTPoints/2)
                             self.dataOut.data_spc = numpy.roll(self.dataOut.data_spc, shift , axis=1)
                             if self.dataOut.data_cspc is not None:
                                 #desplaza a la derecha en el eje 2 determinadas posiciones
                                 self.dataOut.data_cspc = numpy.roll(self.dataOut.data_cspc, shift, axis=1)
                         if pairsList:
                             self.__selectPairs(pairsList)
                     elif self.dataIn.type == "Voltage":
                         self.dataOut.flagNoData = True
                         if nFFTPoints == None:
                             raise ValueError("This SpectraProc.run() need nFFTPoints input variable")
                         if nProfiles == None:
                             nProfiles = nFFTPoints
                         if ippFactor == None:
                             self.dataOut.ippFactor = 1
                         self.dataOut.nFFTPoints = nFFTPoints
                         if self.buffer is None:
                             self.buffer = numpy.zeros((self.dataIn.nChannels,
                                                        nProfiles,
                                                        self.dataIn.nHeights),
                                                       dtype='complex')
                         if self.dataIn.flagDataAsBlock:
                             nVoltProfiles = self.dataIn.data.shape[1]
                             if nVoltProfiles == nProfiles:
                                 self.buffer = self.dataIn.data.copy()
                                 self.profIndex = nVoltProfiles
                             elif nVoltProfiles < nProfiles:
                                 if self.profIndex == 0:
                                     self.id_min = 0
                                     self.id_max = nVoltProfiles
                                 self.buffer[:, self.id_min:self.id_max,
                                             :] = self.dataIn.data
                                 self.profIndex += nVoltProfiles
                                 self.id_min += nVoltProfiles
                                 self.id_max += nVoltProfiles
                             else:
                                 raise ValueError("The type object %s has %d profiles, it should just has %d profiles" % (
                                     self.dataIn.type, self.dataIn.data.shape[1], nProfiles))
                                 self.dataOut.flagNoData = True
                         else:
                             self.buffer[:, self.profIndex, :] = self.dataIn.data.copy()
                             self.profIndex += 1
                         if self.firstdatatime == None:
                             self.firstdatatime = self.dataIn.utctime
                         if self.profIndex == nProfiles:
                             self.__updateSpecFromVoltage()
                             if pairsList == None:
                                 self.dataOut.pairsList = [pair for pair in itertools.combinations(self.dataOut.channelList, 2)]
                             else:
                                 self.dataOut.pairsList = pairsList
                             self.__getFft()
                             self.dataOut.flagNoData = False
                             self.firstdatatime = None
                             self.profIndex = 0
                     else:
                         raise ValueError("The type of input object '%s' is not valid".format(
                             self.dataIn.type))
                 def __selectPairs(self, pairsList):
                     if not pairsList:
                         return
                     pairs = []
                     pairsIndex = []
                     for pair in pairsList:
                         if pair[0] not in self.dataOut.channelList or pair[1] not in self.dataOut.channelList:
                             continue
                         pairs.append(pair)
                         pairsIndex.append(pairs.index(pair))
                     self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndex]
                     self.dataOut.pairsList = pairs
                     return
                 def selectFFTs(self, minFFT, maxFFT ):
                     """
                     Selecciona un bloque de datos en base a un grupo de valores de puntos FFTs segun el rango
                     minFFT<= FFT <= maxFFT
                     """
                     if (minFFT > maxFFT):
                         raise ValueError("Error selecting heights: Height range (%d,%d) is not valid" % (minFFT, maxFFT))
                     if (minFFT < self.dataOut.getFreqRange()[0]):
                         minFFT = self.dataOut.getFreqRange()[0]
                     if (maxFFT > self.dataOut.getFreqRange()[-1]):
                         maxFFT = self.dataOut.getFreqRange()[-1]
                     minIndex = 0
                     maxIndex = 0
                     FFTs = self.dataOut.getFreqRange()
                     inda = numpy.where(FFTs >= minFFT)
                     indb = numpy.where(FFTs <= maxFFT)
                     try:
                         minIndex = inda[0][0]
                     except:
                         minIndex = 0
                     try:
                         maxIndex = indb[0][-1]
                     except:
                         maxIndex = len(FFTs)
                     self.selectFFTsByIndex(minIndex, maxIndex)
                     return 1
                 def getBeaconSignal(self, tauindex=0, channelindex=0, hei_ref=None):
                     newheis = numpy.where(
                         self.dataOut.heightList > self.dataOut.radarControllerHeaderObj.Taus[tauindex])
                     if hei_ref != None:
                         newheis = numpy.where(self.dataOut.heightList > hei_ref)
                     minIndex = min(newheis[0])
                     maxIndex = max(newheis[0])
                     data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
                     heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
                     # determina indices
                     nheis = int(self.dataOut.radarControllerHeaderObj.txB /
                                 (self.dataOut.heightList[1] - self.dataOut.heightList[0]))
                     avg_dB = 10 * \
                         numpy.log10(numpy.sum(data_spc[channelindex, :, :], axis=0))
                     beacon_dB = numpy.sort(avg_dB)[-nheis:]
                     beacon_heiIndexList = []
                     for val in avg_dB.tolist():
                         if val >= beacon_dB[0]:
                             beacon_heiIndexList.append(avg_dB.tolist().index(val))
                     #data_spc = data_spc[:,:,beacon_heiIndexList]
                     data_cspc = None
                     if self.dataOut.data_cspc is not None:
                         data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
                         #data_cspc = data_cspc[:,:,beacon_heiIndexList]
                     data_dc = None
                     if self.dataOut.data_dc is not None:
                         data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
                         #data_dc = data_dc[:,beacon_heiIndexList]
                     self.dataOut.data_spc = data_spc
                     self.dataOut.data_cspc = data_cspc
                     self.dataOut.data_dc = data_dc
                     self.dataOut.heightList = heightList
                     self.dataOut.beacon_heiIndexList = beacon_heiIndexList
                     return 1
                 def selectFFTsByIndex(self, minIndex, maxIndex):
                     """
                     """
                     if (minIndex < 0) or (minIndex > maxIndex):
                         raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (minIndex, maxIndex))
                     if (maxIndex >= self.dataOut.nProfiles):
                         maxIndex = self.dataOut.nProfiles-1
                     #Spectra
                     data_spc = self.dataOut.data_spc[:,minIndex:maxIndex+1,:]
                     data_cspc = None
                     if self.dataOut.data_cspc is not None:
                         data_cspc = self.dataOut.data_cspc[:,minIndex:maxIndex+1,:]
                     data_dc = None
                     if self.dataOut.data_dc is not None:
                         data_dc = self.dataOut.data_dc[minIndex:maxIndex+1,:]
                     self.dataOut.data_spc = data_spc
                     self.dataOut.data_cspc = data_cspc
                     self.dataOut.data_dc = data_dc
                     self.dataOut.ippSeconds = self.dataOut.ippSeconds*(self.dataOut.nFFTPoints / numpy.shape(data_cspc)[1])
                     self.dataOut.nFFTPoints = numpy.shape(data_cspc)[1]
                     self.dataOut.profilesPerBlock = numpy.shape(data_cspc)[1]
                     return 1
                 def getNoise(self, minHei=None, maxHei=None, minVel=None, maxVel=None):
                     # validacion de rango
                     if minHei == None:
                         minHei = self.dataOut.heightList[0]
                     if maxHei == None:
                         maxHei = self.dataOut.heightList[-1]
                     if (minHei < self.dataOut.heightList[0]) or (minHei > maxHei):
                         print('minHei: %.2f is out of the heights range' % (minHei))
                         print('minHei is setting to %.2f' % (self.dataOut.heightList[0]))
                         minHei = self.dataOut.heightList[0]
                     if (maxHei > self.dataOut.heightList[-1]) or (maxHei < minHei):
                         print('maxHei: %.2f is out of the heights range' % (maxHei))
                         print('maxHei is setting to %.2f' % (self.dataOut.heightList[-1]))
                         maxHei = self.dataOut.heightList[-1]
                     # validacion de velocidades
                     velrange = self.dataOut.getVelRange(1)
                     if minVel == None:
                         minVel = velrange[0]
                     if maxVel == None:
                         maxVel = velrange[-1]
                     if (minVel < velrange[0]) or (minVel > maxVel):
                         print('minVel: %.2f is out of the velocity range' % (minVel))
                         print('minVel is setting to %.2f' % (velrange[0]))
                         minVel = velrange[0]
                     if (maxVel > velrange[-1]) or (maxVel < minVel):
                         print('maxVel: %.2f is out of the velocity range' % (maxVel))
                         print('maxVel is setting to %.2f' % (velrange[-1]))
                         maxVel = velrange[-1]
                     # seleccion de indices para rango
                     minIndex = 0
                     maxIndex = 0
                     heights = self.dataOut.heightList
                     inda = numpy.where(heights >= minHei)
                     indb = numpy.where(heights <= maxHei)
                     try:
                         minIndex = inda[0][0]
                     except:
                         minIndex = 0
                     try:
                         maxIndex = indb[0][-1]
                     except:
                         maxIndex = len(heights)
                     if (minIndex < 0) or (minIndex > maxIndex):
                         raise ValueError("some value in (%d,%d) is not valid" % (
                             minIndex, maxIndex))
                     if (maxIndex >= self.dataOut.nHeights):
                         maxIndex = self.dataOut.nHeights - 1
                     # seleccion de indices para velocidades
                     indminvel = numpy.where(velrange >= minVel)
                     indmaxvel = numpy.where(velrange <= maxVel)
                     try:
                         minIndexVel = indminvel[0][0]
                     except:
                         minIndexVel = 0
                     try:
                         maxIndexVel = indmaxvel[0][-1]
                     except:
                         maxIndexVel = len(velrange)
                     # seleccion del espectro
                     data_spc = self.dataOut.data_spc[:,
                                                      minIndexVel:maxIndexVel + 1, minIndex:maxIndex + 1]
                     # estimacion de ruido
                     noise = numpy.zeros(self.dataOut.nChannels)
                     for channel in range(self.dataOut.nChannels):
                         daux = data_spc[channel, :, :]
                         sortdata = numpy.sort(daux, axis=None)
                         noise[channel] = hildebrand_sekhon(sortdata, self.dataOut.nIncohInt)
                     self.dataOut.noise_estimation = noise.copy()
                     return 1
             class removeDC(Operation):
                 def run(self, dataOut, mode=2):
                     self.dataOut = dataOut
                     jspectra = self.dataOut.data_spc
                     jcspectra = self.dataOut.data_cspc
                     num_chan = jspectra.shape[0]
                     num_hei = jspectra.shape[2]
                     if jcspectra is not None:
                         jcspectraExist = True
                         num_pairs = jcspectra.shape[0]
                     else:
                         jcspectraExist = False
                     freq_dc = int(jspectra.shape[1] / 2)
                     ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
                     ind_vel = ind_vel.astype(int)
                     if ind_vel[0] < 0:
                         ind_vel[list(range(0, 1))] = ind_vel[list(range(0, 1))] + self.num_prof
                     if mode == 1:
                         jspectra[:, freq_dc, :] = (
                             jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2  # CORRECCION
                         if jcspectraExist:
                             jcspectra[:, freq_dc, :] = (
                                 jcspectra[:, ind_vel[1], :] + jcspectra[:, ind_vel[2], :]) / 2
                     if mode == 2:
                         vel = numpy.array([-2, -1, 1, 2])
                         xx = numpy.zeros([4, 4])
                         for fil in range(4):
                             xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
                         xx_inv = numpy.linalg.inv(xx)
                         xx_aux = xx_inv[0, :]
                         for ich in range(num_chan):
                             yy = jspectra[ich, ind_vel, :]
                             jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
                             junkid = jspectra[ich, freq_dc, :] <= 0
                             cjunkid = sum(junkid)
                             if cjunkid.any():
                                 jspectra[ich, freq_dc, junkid.nonzero()] = (
                                     jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
                         if jcspectraExist:
                             for ip in range(num_pairs):
                                 yy = jcspectra[ip, ind_vel, :]
                                 jcspectra[ip, freq_dc, :] = numpy.dot(xx_aux, yy)
                     self.dataOut.data_spc = jspectra
                     self.dataOut.data_cspc = jcspectra
                     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
                 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
                     #print("CHANNELS: ",[x for x in self.channels])
                     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
                     #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.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
+                    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
                     #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)
+                    ###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)
+                            ###ONLY FOR TEST:
-                            # fig2, axs2 = plt.subplots(raxs, caxs)
+                            if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
-                            # col_ax = 0
+                                fig, axs = plt.subplots(raxs, caxs)
-                            # row_ax = 0
+                                fig2, axs2 = plt.subplots(raxs, caxs)
-                            #print(len(self.nPairs))
+                            col_ax = 0
-                            for ii in range(self.nPairs):                       #PARES DE CANALES SELF y CROSS
+                            row_ax = 0
-                                #print("ii: ",ii)
+                            #print(self.nPairs)
-                                # if (col_ax%caxs==0 and col_ax!=0):
+                            for ii in range(self.nPairs):  #PARES DE CANALES SELF y CROSS
-                                #     col_ax = 0
+                                if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
-                                #     row_ax += 1
+                                    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()
                                 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
+                                    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])
+                                        ###ONLY FOR TEST:
-                                        axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
+                                        if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
-                                        axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
+                                            newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
-                                        axs[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))'''
+                                            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)
                                     #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("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("   ",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")
-                                            #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
-                                    '''
+                                    ###ONLY FOR TEST:
-                                    func2clean = 10*numpy.log10(numpy.absolute(cspectra[:,ii,ifreq,ih]))
+                                    if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
-                                    y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
+                                        func2clean = 10*numpy.log10(numpy.absolute(cspectra[:,ii,ifreq,ih]))
-                                    axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
+                                        y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
-                                    axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
+                                        axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
-                                    axs2[row_ax,col_ax].set_title("Pair "+str(self.crosspairs[ii]))
+                                        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)
-                                '''
+                            ###ONLY FOR TEST:
-                            title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+"   Alt: "+str(self.heights[ih])+ " km"
+                            if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
-                            title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+"   Alt: "+str(self.heights[ih])+ " km  CLEANED"
+                                title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+"   Alt: "+str(self.heights[ih])+ " km"
-                            fig.suptitle(title)
+                                title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+"   Alt: "+str(self.heights[ih])+ " km  CLEANED"
-                            fig2.suptitle(title2)
+                                fig.suptitle(title)
-                            plt.show()'''
+                                fig2.suptitle(title2)
+                                plt.show()
-                        '''                channels = channels
+                        '''
+                                       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):
                     cspc = self.dataOut.data_cspc
                     spc = self.dataOut.data_spc
                     Heights = numpy.arange(cspc.shape[2])
                     realCspc = numpy.abs(cspc)
                     for i in range(cspc.shape[0]):
                         LinePower= numpy.sum(realCspc[i], axis=0)
                         Threshold = numpy.amax(LinePower)-numpy.sort(LinePower)[len(Heights)-int(len(Heights)*0.1)]
                         SelectedHeights = Heights[ numpy.where( LinePower < Threshold ) ]
                         InterferenceSum = numpy.sum( realCspc[i,:,SelectedHeights], axis=0 )
                         InterferenceThresholdMin = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.98)]
                         InterferenceThresholdMax = numpy.sort(InterferenceSum)[int(len(InterferenceSum)*0.99)]
                         InterferenceRange = numpy.where( ([InterferenceSum > InterferenceThresholdMin]))# , InterferenceSum < InterferenceThresholdMax]) )
                         #InterferenceRange = numpy.where( ([InterferenceRange < InterferenceThresholdMax]))
                         if len(InterferenceRange)<int(cspc.shape[1]*0.3):
                             cspc[i,InterferenceRange,:] = numpy.NaN
                     self.dataOut.data_cspc = cspc
                 def removeInterference(self, interf = 2, hei_interf = None, nhei_interf = None, offhei_interf = None):
                     jspectra = self.dataOut.data_spc
                     jcspectra = self.dataOut.data_cspc
                     jnoise = self.dataOut.getNoise()
                     num_incoh = self.dataOut.nIncohInt
                     num_channel = jspectra.shape[0]
                     num_prof = jspectra.shape[1]
                     num_hei = jspectra.shape[2]
                     # hei_interf
                     if hei_interf is None:
                         count_hei = int(num_hei / 2)
                         hei_interf = numpy.asmatrix(list(range(count_hei))) + num_hei - count_hei
                         hei_interf = numpy.asarray(hei_interf)[0]
                     # nhei_interf
                     if (nhei_interf == None):
                         nhei_interf = 5
                     if (nhei_interf < 1):
                         nhei_interf = 1
                     if (nhei_interf > count_hei):
                         nhei_interf = count_hei
                     if (offhei_interf == None):
                         offhei_interf = 0
                     ind_hei = list(range(num_hei))
             #         mask_prof = numpy.asarray(range(num_prof - 2)) + 1
             #         mask_prof[range(num_prof/2 - 1,len(mask_prof))] += 1
                     mask_prof = numpy.asarray(list(range(num_prof)))
                     num_mask_prof = mask_prof.size
                     comp_mask_prof = [0, num_prof / 2]
                     # noise_exist:    Determina si la variable jnoise ha sido definida y contiene la informacion del ruido de cada canal
                     if (jnoise.size < num_channel or numpy.isnan(jnoise).any()):
                         jnoise = numpy.nan
                     noise_exist = jnoise[0] < numpy.Inf
                     # Subrutina de Remocion de la Interferencia
                     for ich in range(num_channel):
                         # Se ordena los espectros segun su potencia (menor a mayor)
                         power = jspectra[ich, mask_prof, :]
                         power = power[:, hei_interf]
                         power = power.sum(axis=0)
                         psort = power.ravel().argsort()
                         # Se estima la interferencia promedio en los Espectros de Potencia empleando
                         junkspc_interf = jspectra[ich, :, hei_interf[psort[list(range(
                             offhei_interf, nhei_interf + offhei_interf))]]]
                         if noise_exist:
                             #    tmp_noise = jnoise[ich] / num_prof
                             tmp_noise = jnoise[ich]
                         junkspc_interf = junkspc_interf - tmp_noise
                         #junkspc_interf[:,comp_mask_prof] = 0
                         jspc_interf = junkspc_interf.sum(axis=0) / nhei_interf
                         jspc_interf = jspc_interf.transpose()
                         # Calculando el espectro de interferencia promedio
                         noiseid = numpy.where(
                             jspc_interf <= tmp_noise / numpy.sqrt(num_incoh))
                         noiseid = noiseid[0]
                         cnoiseid = noiseid.size
                         interfid = numpy.where(
                             jspc_interf > tmp_noise / numpy.sqrt(num_incoh))
                         interfid = interfid[0]
                         cinterfid = interfid.size
                         if (cnoiseid > 0):
                             jspc_interf[noiseid] = 0
                         # Expandiendo los perfiles a limpiar
                         if (cinterfid > 0):
                             new_interfid = (
                                 numpy.r_[interfid - 1, interfid, interfid + 1] + num_prof) % num_prof
                             new_interfid = numpy.asarray(new_interfid)
                             new_interfid = {x for x in new_interfid}
                             new_interfid = numpy.array(list(new_interfid))
                             new_cinterfid = new_interfid.size
                         else:
                             new_cinterfid = 0
                         for ip in range(new_cinterfid):
                             ind = junkspc_interf[:, new_interfid[ip]].ravel().argsort()
                             jspc_interf[new_interfid[ip]
                                         ] = junkspc_interf[ind[nhei_interf // 2], new_interfid[ip]]
                         jspectra[ich, :, ind_hei] = jspectra[ich, :,
                                                              ind_hei] - jspc_interf  # Corregir indices
                         # Removiendo la interferencia del punto de mayor interferencia
                         ListAux = jspc_interf[mask_prof].tolist()
                         maxid = ListAux.index(max(ListAux))
                         if cinterfid > 0:
                             for ip in range(cinterfid * (interf == 2) - 1):
                                 ind = (jspectra[ich, interfid[ip], :] < tmp_noise *
                                        (1 + 1 / numpy.sqrt(num_incoh))).nonzero()
                                 cind = len(ind)
                                 if (cind > 0):
                                     jspectra[ich, interfid[ip], ind] = tmp_noise * \
                                         (1 + (numpy.random.uniform(cind) - 0.5) /
                                          numpy.sqrt(num_incoh))
                             ind = numpy.array([-2, -1, 1, 2])
                             xx = numpy.zeros([4, 4])
                             for id1 in range(4):
                                 xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
                             xx_inv = numpy.linalg.inv(xx)
                             xx = xx_inv[:, 0]
                             ind = (ind + maxid + num_mask_prof) % num_mask_prof
                             yy = jspectra[ich, mask_prof[ind], :]
                             jspectra[ich, mask_prof[maxid], :] = numpy.dot(
                                 yy.transpose(), xx)
                         indAux = (jspectra[ich, :, :] < tmp_noise *
                                   (1 - 1 / numpy.sqrt(num_incoh))).nonzero()
                         jspectra[ich, indAux[0], indAux[1]] = tmp_noise * \
                             (1 - 1 / numpy.sqrt(num_incoh))
                     # Remocion de Interferencia en el Cross Spectra
                     if jcspectra is None:
                         return jspectra, jcspectra
                     num_pairs = int(jcspectra.size / (num_prof * num_hei))
                     jcspectra = jcspectra.reshape(num_pairs, num_prof, num_hei)
                     for ip in range(num_pairs):
                         #-------------------------------------------
                         cspower = numpy.abs(jcspectra[ip, mask_prof, :])
                         cspower = cspower[:, hei_interf]
                         cspower = cspower.sum(axis=0)
                         cspsort = cspower.ravel().argsort()
                         junkcspc_interf = jcspectra[ip, :, hei_interf[cspsort[list(range(
                             offhei_interf, nhei_interf + offhei_interf))]]]
                         junkcspc_interf = junkcspc_interf.transpose()
                         jcspc_interf = junkcspc_interf.sum(axis=1) / nhei_interf
                         ind = numpy.abs(jcspc_interf[mask_prof]).ravel().argsort()
                         median_real = int(numpy.median(numpy.real(
                             junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
                         median_imag = int(numpy.median(numpy.imag(
                             junkcspc_interf[mask_prof[ind[list(range(3 * num_prof // 4))]], :])))
                         comp_mask_prof = [int(e) for e in comp_mask_prof]
                         junkcspc_interf[comp_mask_prof, :] = numpy.complex(
                             median_real, median_imag)
                         for iprof in range(num_prof):
                             ind = numpy.abs(junkcspc_interf[iprof, :]).ravel().argsort()
                             jcspc_interf[iprof] = junkcspc_interf[iprof, ind[nhei_interf // 2]]
                         # Removiendo la Interferencia
                         jcspectra[ip, :, ind_hei] = jcspectra[ip,
                                                               :, ind_hei] - jcspc_interf
                         ListAux = numpy.abs(jcspc_interf[mask_prof]).tolist()
                         maxid = ListAux.index(max(ListAux))
                         ind = numpy.array([-2, -1, 1, 2])
                         xx = numpy.zeros([4, 4])
                         for id1 in range(4):
                             xx[:, id1] = ind[id1]**numpy.asarray(list(range(4)))
                         xx_inv = numpy.linalg.inv(xx)
                         xx = xx_inv[:, 0]
                         ind = (ind + maxid + num_mask_prof) % num_mask_prof
                         yy = jcspectra[ip, mask_prof[ind], :]
                         jcspectra[ip, mask_prof[maxid], :] = numpy.dot(yy.transpose(), xx)
                     # Guardar Resultados
                     self.dataOut.data_spc = jspectra
                     self.dataOut.data_cspc = jcspectra
                     return 1
                 def run(self, dataOut, interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None, mode=1):
                     self.dataOut = dataOut
                     if mode == 1:
                         self.removeInterference(interf = 2,hei_interf = None, nhei_interf = None, offhei_interf = None)
                     elif mode == 2:
                         self.removeInterference2()
                     return self.dataOut
             class IncohInt(Operation):
                 __profIndex = 0
                 __withOverapping = False
                 __byTime = False
                 __initime = None
                 __lastdatatime = None
                 __integrationtime = None
                 __buffer_spc = None
                 __buffer_cspc = None
                 __buffer_dc = None
                 __dataReady = False
                 __timeInterval = None
                 n = None
                 def __init__(self):
                     Operation.__init__(self)
                 def setup(self, n=None, timeInterval=None, overlapping=False):
                     """
                     Set the parameters of the integration class.
                     Inputs:
                         n        :    Number of coherent integrations
                         timeInterval   :    Time of integration. If the parameter "n" is selected this one does not work
                         overlapping    :
                     """
                     self.__initime = None
                     self.__lastdatatime = 0
                     self.__buffer_spc = 0
                     self.__buffer_cspc = 0
                     self.__buffer_dc = 0
                     self.__profIndex = 0
                     self.__dataReady = False
                     self.__byTime = False
                     if n is None and timeInterval is None:
                         raise ValueError("n or timeInterval should be specified ...")
                     if n is not None:
                         self.n = int(n)
                     else:
                         self.__integrationtime = int(timeInterval)
                         self.n = None
                         self.__byTime = True
                 def putData(self, data_spc, data_cspc, data_dc):
                     """
                     Add a profile to the __buffer_spc and increase in one the __profileIndex
                     """
                     self.__buffer_spc += data_spc
                     if data_cspc is None:
                         self.__buffer_cspc = None
                     else:
                         self.__buffer_cspc += data_cspc
                     if data_dc is None:
                         self.__buffer_dc = None
                     else:
                         self.__buffer_dc += data_dc
                     self.__profIndex += 1
                     return
                 def pushData(self):
                     """
                     Return the sum of the last profiles and the profiles used in the sum.
                     Affected:
                     self.__profileIndex
                     """
                     data_spc = self.__buffer_spc
                     data_cspc = self.__buffer_cspc
                     data_dc = self.__buffer_dc
                     n = self.__profIndex
                     self.__buffer_spc = 0
                     self.__buffer_cspc = 0
                     self.__buffer_dc = 0
                     self.__profIndex = 0
                     return data_spc, data_cspc, data_dc, n
                 def byProfiles(self, *args):
                     self.__dataReady = False
                     avgdata_spc = None
                     avgdata_cspc = None
                     avgdata_dc = None
                     self.putData(*args)
                     if self.__profIndex == self.n:
                         avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
                         self.n = n
                         self.__dataReady = True
                     return avgdata_spc, avgdata_cspc, avgdata_dc
                 def byTime(self, datatime, *args):
                     self.__dataReady = False
                     avgdata_spc = None
                     avgdata_cspc = None
                     avgdata_dc = None
                     self.putData(*args)
                     if (datatime - self.__initime) >= self.__integrationtime:
                         avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
                         self.n = n
                         self.__dataReady = True
                     return avgdata_spc, avgdata_cspc, avgdata_dc
                 def integrate(self, datatime, *args):
                     if self.__profIndex == 0:
                         self.__initime = datatime
                     if self.__byTime:
                         avgdata_spc, avgdata_cspc, avgdata_dc = self.byTime(
                             datatime, *args)
                     else:
                         avgdata_spc, avgdata_cspc, avgdata_dc = self.byProfiles(*args)
                     if not self.__dataReady:
                         return None, None, None, None
                     return self.__initime, avgdata_spc, avgdata_cspc, avgdata_dc
                 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
                     if n == 1:
                         return dataOut
                     dataOut.flagNoData = True
                     if not self.isConfig:
                         self.setup(n, timeInterval, overlapping)
                         self.isConfig = True
                     avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(dataOut.utctime,
                                                                                         dataOut.data_spc,
                                                                                         dataOut.data_cspc,
                                                                                         dataOut.data_dc)
                     if self.__dataReady:
                         dataOut.data_spc = avgdata_spc
                         dataOut.data_cspc = avgdata_cspc
                         dataOut.data_dc = avgdata_dc
                         dataOut.nIncohInt *= self.n
                         dataOut.utctime = avgdatatime
                         dataOut.flagNoData = False
                     return dataOut
             class dopplerFlip(Operation):
                 def run(self, dataOut):
                     # arreglo 1: (num_chan, num_profiles, num_heights)
                     self.dataOut = dataOut
                     # JULIA-oblicua, indice 2
                     # arreglo 2: (num_profiles, num_heights)
                     jspectra = self.dataOut.data_spc[2]
                     jspectra_tmp = numpy.zeros(jspectra.shape)
                     num_profiles = jspectra.shape[0]
                     freq_dc = int(num_profiles / 2)
                     # Flip con for
                     for j in range(num_profiles):
                         jspectra_tmp[num_profiles-j-1]= jspectra[j]
                     # Intercambio perfil de DC con perfil inmediato anterior
                     jspectra_tmp[freq_dc-1]= jspectra[freq_dc-1]
                     jspectra_tmp[freq_dc]= jspectra[freq_dc]
                     # canal modificado es re-escrito en el arreglo de canales
                     self.dataOut.data_spc[2] = jspectra_tmp
                     return self.dataOut

General Comments 0

Write
Preview

You need to be logged in to leave comments. Login now

No TODOs yet

	Site-wide shortcuts
/	Use quick search box
g h	Goto home page
g g	Goto my private gists page
g G	Goto my public gists page
g 0-9	Goto bookmarked items from 0-9
n r	New repository page
n g	New gist page

	Repositories
g s	Goto summary page
g c	Goto changelog page
g f	Goto files page
g F	Goto files page with file search activated
g p	Goto pull requests page
g o	Goto repository settings
g O	Goto repository access permissions settings
t s	Toggle sidebar on some pages