schain Commit - r1545:af752b8a14b4 · Jicamarca Repository

pc amisr

joabAM -

r1545:af752b8a14b4

parent child

Context file:

r1545:af752b8a14b4

Collapse all files

schainpy/model/proc/jroproc_spectra.py +1 -1

             # 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.model.data import _noise
             from schainpy.utils import log
             import matplotlib.pyplot as plt
             #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):
                     #print("run spc proc")
                     try:
                         type = self.dataIn.type.decode("utf-8")
                         self.dataIn.type = type
                     except:
                         pass
                     if self.dataIn.type == "Spectra":
                         try:
                             self.dataOut.copy(self.dataIn)
                         except Exception as e:
                             print("Error dataIn ",e)
                         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
                         #print(" volts ch,prof, h: ", self.dataIn.data.shape)
                         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
                     elif self.dataIn.type == "Parameters":
                         self.dataOut.data_spc = self.dataIn.data_spc
                         self.dataOut.data_cspc = self.dataIn.data_cspc
                         self.dataOut.data_outlier = self.dataIn.data_outlier
                         self.dataOut.nProfiles = self.dataIn.nProfiles
                         self.dataOut.nIncohInt = self.dataIn.nIncohInt
                         self.dataOut.nFFTPoints = self.dataIn.nFFTPoints
                         self.dataOut.ippFactor = self.dataIn.ippFactor
                         self.dataOut.max_nIncohInt = self.dataIn.max_nIncohInt
                         self.dataOut.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy()
                         self.dataOut.ipp = self.dataIn.ipp
                         #self.dataOut.abscissaList = self.dataIn.getVelRange(1)
                         #self.dataOut.spc_noise = self.dataIn.getNoise()
                         #self.dataOut.spc_range = (self.dataIn.getFreqRange(1) , self.dataIn.getAcfRange(1) , self.dataIn.getVelRange(1))
                         # self.dataOut.normFactor = self.dataIn.normFactor
                         if hasattr(self.dataIn, 'channelList'):
                             self.dataOut.channelList = self.dataIn.channelList
                         if hasattr(self.dataIn, 'pairsList'):
                             self.dataOut.pairsList = self.dataIn.pairsList
                             self.dataOut.groupList = self.dataIn.pairsList
                         self.dataOut.flagNoData = False
                         if hasattr(self.dataIn, 'ChanDist'): #Distances of receiver channels
                             self.dataOut.ChanDist = self.dataIn.ChanDist
                         else: self.dataOut.ChanDist = None
                         #if hasattr(self.dataIn, 'VelRange'): #Velocities range
                         #    self.dataOut.VelRange = self.dataIn.VelRange
                         #else: self.dataOut.VelRange = None
                     else:
                         raise ValueError("The type of input object {} 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
             class getNoiseB(Operation):
                 __slots__ =('offset','warnings', 'isConfig', 'minIndex','maxIndex','minIndexFFT','maxIndexFFT')
                 def __init__(self):
                     Operation.__init__(self)
                     self.isConfig = False
                 def setup(self, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
                     self.warnings = warnings
                     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):
                         if self.warnings:
                             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):
                         if self.warnings:
                             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]
                     #indices relativos a los puntos de fft, puede ser de acuerdo a velocidad o frecuencia
                     minIndexFFT = 0
                     maxIndexFFT = 0
                     # validacion de velocidades
                     indminPoint = None
                     indmaxPoint = None
                     if self.dataOut.type == 'Spectra':
                         if minVel == None and maxVel == None :
                             freqrange = self.dataOut.getFreqRange(1)
                             if minFreq == None:
                                 minFreq = freqrange[0]
                             if maxFreq == None:
                                 maxFreq = freqrange[-1]
                             if (minFreq < freqrange[0]) or (minFreq > maxFreq):
                                 if self.warnings:
                                     print('minFreq: %.2f is out of the frequency range' % (minFreq))
                                     print('minFreq is setting to %.2f' % (freqrange[0]))
                                 minFreq = freqrange[0]
                             if (maxFreq > freqrange[-1]) or (maxFreq < minFreq):
                                 if self.warnings:
                                     print('maxFreq: %.2f is out of the frequency range' % (maxFreq))
                                     print('maxFreq is setting to %.2f' % (freqrange[-1]))
                                 maxFreq = freqrange[-1]
                             indminPoint = numpy.where(freqrange >= minFreq)
                             indmaxPoint = numpy.where(freqrange <= maxFreq)
                         else:
                             velrange = self.dataOut.getVelRange(1)
                             if minVel == None:
                                 minVel = velrange[0]
                             if maxVel == None:
                                 maxVel = velrange[-1]
                             if (minVel < velrange[0]) or (minVel > maxVel):
                                 if self.warnings:
                                     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):
                                 if self.warnings:
                                     print('maxVel: %.2f is out of the velocity range' % (maxVel))
                                     print('maxVel is setting to %.2f' % (velrange[-1]))
                                 maxVel = velrange[-1]
                             indminPoint = numpy.where(velrange >= minVel)
                             indmaxPoint = numpy.where(velrange <= maxVel)
                     # 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
                     #############################################################3
                     # seleccion de indices para velocidades
                     if self.dataOut.type == 'Spectra':
                         try:
                             minIndexFFT = indminPoint[0][0]
                         except:
                             minIndexFFT = 0
                         try:
                             maxIndexFFT = indmaxPoint[0][-1]
                         except:
                             maxIndexFFT = len( self.dataOut.getFreqRange(1))
                     self.minIndex, self.maxIndex, self.minIndexFFT, self.maxIndexFFT = minIndex, maxIndex, minIndexFFT, maxIndexFFT
                     self.isConfig = True
                     if offset!=None:
                         self.offset = 10**(offset/10)
                     #print("config getNoiseB Done")
                 def run(self, dataOut, offset=None, minHei=None, maxHei=None,minVel=None, maxVel=None, minFreq= None, maxFreq=None, warnings=False):
                     self.dataOut = dataOut
                     if not self.isConfig:
                         self.setup(offset, minHei, maxHei,minVel, maxVel, minFreq, maxFreq, warnings)
                     self.dataOut.noise_estimation = None
                     noise = None
                     if self.dataOut.type == 'Voltage':
                         noise = self.dataOut.getNoise(ymin_index=self.minIndex, ymax_index=self.maxIndex)
                         #print(minIndex, maxIndex,minIndexVel, maxIndexVel)
                     elif self.dataOut.type == 'Spectra':
                         #print(self.minIndex, self.maxIndex,self.minIndexFFT, self.maxIndexFFT, self.dataOut.nIncohInt)
                         noise = numpy.zeros( self.dataOut.nChannels)
                         norm = 1
                         for channel in range( self.dataOut.nChannels):
                             if not hasattr(self.dataOut.nIncohInt,'__len__'):
                                 norm = 1
                             else:
                                 norm = self.dataOut.max_nIncohInt/self.dataOut.nIncohInt[channel, self.minIndex:self.maxIndex]
                             #print("norm nIncoh: ", norm ,self.dataOut.data_spc.shape)
                             daux =  self.dataOut.data_spc[channel,self.minIndexFFT:self.maxIndexFFT, self.minIndex:self.maxIndex]
                             daux = numpy.multiply(daux, norm)
                             #print("offset: ", self.offset, 10*numpy.log10(self.offset))
                             #noise[channel] = self.getNoiseByMean(daux)/self.offset
                             #print(daux.shape, daux)
                             noise[channel] = self.getNoiseByHS(daux, self.dataOut.max_nIncohInt)/self.offset
                             # data = numpy.mean(daux,axis=1)
                             # sortdata = numpy.sort(data, axis=None)
                             # noise[channel] = _noise.hildebrand_sekhon(sortdata, self.dataOut.max_nIncohInt)/self.offset
                         #noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
                     else:
                         noise = self.dataOut.getNoise(xmin_index=self.minIndexFFT, xmax_index=self.maxIndexFFT, ymin_index=self.minIndex, ymax_index=self.maxIndex)
                     self.dataOut.noise_estimation = noise.copy() # dataOut.noise
                     #print("2: ",10*numpy.log10(self.dataOut.noise_estimation/64))
                     #print(self.dataOut.flagNoData)
-                    print("getNoise Done")
+                    print("getNoise Done", noise)
                     return self.dataOut
                 def getNoiseByMean(self,data):
                     #data debe estar ordenado
                     data = numpy.mean(data,axis=1)
                     sortdata = numpy.sort(data, axis=None)
                     #sortID=data.argsort()
                     #print(data.shape)
                     pnoise = None
                     j = 0
                     mean = numpy.mean(sortdata)
                     min = numpy.min(sortdata)
                     delta = mean - min
                     indexes = numpy.where(sortdata > (mean+delta))[0] #only array of indexes
                     #print(len(indexes))
                     if len(indexes)==0:
                         pnoise = numpy.mean(sortdata)
                     else:
                         j = indexes[0]
                         pnoise = numpy.mean(sortdata[0:j])
                     # from matplotlib import pyplot as plt
                     # plt.plot(sortdata)
                     # plt.vlines(j,(pnoise-delta),(pnoise+delta), color='r')
                     # plt.show()
                     #print("noise: ", 10*numpy.log10(pnoise))
                     return pnoise
                 def getNoiseByHS(self,data, navg):
                     #data debe estar ordenado
                     #data = numpy.mean(data,axis=1)
                     sortdata = numpy.sort(data, axis=None)
                     lenOfData = len(sortdata)
                     nums_min = lenOfData*0.05
                     if nums_min <= 5:
                         nums_min = 5
                     sump = 0.
                     sumq = 0.
                     j = 0
                     cont = 1
                     while((cont == 1)and(j < lenOfData)):
                         sump += sortdata[j]
                         sumq += sortdata[j]**2
                         #sumq -= sump**2
                         if j > nums_min:
                             rtest = float(j)/(j-1) + 1.0/0.1
                             #if ((sumq*j) > (sump**2)):
                             if ((sumq*j) > (rtest*sump**2)):
                                 j = j - 1
                                 sump = sump - sortdata[j]
                                 sumq = sumq - sortdata[j]**2
                                 cont = 0
                         j += 1
                     lnoise = sump / j
                     return lnoise
             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
                     if n != None :
                         self.byProfiles = True
                         self.nIntProfiles = n
                     else:
                         self.__integrationtime = timeInterval
                     self.__dataReady = False
                     self.isConfig = True
                 def run(self, dataOut,min_hei=None,max_hei=None, n=None, timeInterval=10,factor_stdv=2.5):
                     #print("runing cleanRayleigh")
                     if not self.isConfig :
                         self.setup(dataOut, min_hei,max_hei,n,timeInterval,factor_stdv)
                     tini=dataOut.utctime
                     if self.byProfiles:
                         if self.__profIndex ==  self.nIntProfiles:
                             self.__dataReady = True
                     else:
                         if (tini - self.__initime) >= self.__integrationtime:
                             self.__dataReady = True
                             self.__initime = tini
                     #if (tini.tm_min % 2) == 0 and (tini.tm_sec < 5 and self.fint==0):
                     if self.__dataReady:
                         self.__profIndex = 0
                         jspc = self.buffer
                         jcspc = self.buffer2
                         #jnoise = self.buffer3
                         self.buffer = dataOut.data_spc
                         self.buffer2 = dataOut.data_cspc
                         #self.buffer3 = dataOut.noise
                         self.currentTime = dataOut.utctime
                         if numpy.any(jspc) :
                             #print( jspc.shape, jcspc.shape)
                             jspc = numpy.reshape(jspc,(int(len(jspc)/self.nChannels),self.nChannels,self.nFFTPoints,self.nHeights))
                             try:
                                 jcspc= numpy.reshape(jcspc,(int(len(jcspc)/self.nPairs),self.nPairs,self.nFFTPoints,self.nHeights))
                             except:
                                 print("no cspc")
                             self.__dataReady = False
                             #print( jspc.shape, jcspc.shape)
                             dataOut.flagNoData = False
                         else:
                             dataOut.flagNoData = True
                             self.__dataReady = False
                             return dataOut
                     else:
                         #print( len(self.buffer))
                         if numpy.any(self.buffer):
                             self.buffer = numpy.concatenate((self.buffer,dataOut.data_spc), axis=0)
                             try:
                                 self.buffer2 = numpy.concatenate((self.buffer2,dataOut.data_cspc), axis=0)
                                 self.buffer3 += dataOut.data_dc
                             except:
                                 pass
                         else:
                             self.buffer = dataOut.data_spc
                             self.buffer2 = dataOut.data_cspc
                             self.buffer3 = dataOut.data_dc
                         #print self.index, self.fint
                         #print self.buffer2.shape
                         dataOut.flagNoData = True                   ## NOTE: ?? revisar LUEGO
                         self.__profIndex += 1
                         return dataOut                             ## NOTE: REV
                     #index = tini.tm_hour*12+tini.tm_min/5
                     '''
                     #REVISAR
                     '''
                     # jspc = jspc/self.nFFTPoints/self.normFactor
                     # jcspc = jcspc/self.nFFTPoints/self.normFactor
                     tmp_spectra,tmp_cspectra = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
                     dataOut.data_spc = tmp_spectra
                     dataOut.data_cspc = tmp_cspectra
                     #dataOut.data_spc,dataOut.data_cspc = self.cleanRayleigh(dataOut,jspc,jcspc,self.factor_stdv)
                     dataOut.data_dc = self.buffer3
                     dataOut.nIncohInt *= self.nIntProfiles
                     dataOut.max_nIncohInt = self.nIntProfiles
                     dataOut.utctime = self.currentTime             #tiempo promediado
                     #print("Time: ",time.localtime(dataOut.utctime))
                     # dataOut.data_spc = sat_spectra
                     # dataOut.data_cspc = sat_cspectra
                     self.buffer = 0
                     self.buffer2 = 0
                     self.buffer3 = 0
                     return dataOut
                 def cleanRayleigh(self,dataOut,spectra,cspectra,factor_stdv):
                     print("OP cleanRayleigh")
                     #import matplotlib.pyplot as plt
                     #for k in range(149):
                     #channelsProcssd = []
                     #channelA_ok = False
                     #rfunc = cspectra.copy() #self.bloques
                     rfunc = spectra.copy()
                     #rfunc = cspectra
                     #val_spc = spectra*0.0 #self.bloque0*0.0
                     #val_cspc = cspectra*0.0 #self.bloques*0.0
                     #in_sat_spectra = spectra.copy()  #self.bloque0
                     #in_sat_cspectra = cspectra.copy()  #self.bloques
                     ###ONLY FOR TEST:
                     raxs = math.ceil(math.sqrt(self.nPairs))
                     if raxs == 0:
                         raxs = 1
                     caxs = math.ceil(self.nPairs/raxs)
                     if self.nPairs <4:
                         raxs = 2
                         caxs = 2
                     #print(raxs, caxs)
                     fft_rev = 14 #nFFT to plot
                     hei_rev = ((self.heights >= 550) & (self.heights <= 551)).nonzero() #hei to plot
                     hei_rev = hei_rev[0]
                     #print(hei_rev)
                     #print numpy.absolute(rfunc[:,0,0,14])
                     gauss_fit, covariance = None, None
                     for ih in range(self.minAltInd,self.maxAltInd):
                         for ifreq in range(self.nFFTPoints):
                             '''
                             ###ONLY FOR TEST:
                             if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
                                 fig, axs = plt.subplots(raxs, caxs)
                                 fig2, axs2 = plt.subplots(raxs, caxs)
                             col_ax = 0
                             row_ax = 0
                             '''
                             #print(self.nPairs)
                             for ii in range(self.nChan):  #PARES DE CANALES SELF y CROSS
                                 # if self.crosspairs[ii][1]-self.crosspairs[ii][0] > 1: # APLICAR SOLO EN PARES CONTIGUOS
                                 #     continue
                                 # if not self.crosspairs[ii][0] in channelsProcssd:
                                 #     channelA_ok = True
                                 #print("pair: ",self.crosspairs[ii])
                                 '''
                                 ###ONLY FOR TEST:
                                 if (col_ax%caxs==0 and col_ax!=0 and self.nPairs !=1):
                                     col_ax = 0
                                     row_ax += 1
                                 '''
                                 func2clean = 10*numpy.log10(numpy.absolute(rfunc[:,ii,ifreq,ih]))           #Potencia?
                                 #print(func2clean.shape)
                                 val = (numpy.isfinite(func2clean)==True).nonzero()
                                 if len(val)>0:                                                              #limitador
                                     min_val = numpy.around(numpy.amin(func2clean)-2) #> (-40)
                                     if min_val <= -40 :
                                         min_val = -40
                                     max_val = numpy.around(numpy.amax(func2clean)+2) #< 200
                                     if max_val >= 200 :
                                         max_val = 200
                                     #print min_val, max_val
                                     step = 1
                                     #print("Getting bins and the histogram")
                                     x_dist = min_val + numpy.arange(1 + ((max_val-(min_val))/step))*step
                                     y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
                                     #print(len(y_dist),len(binstep[:-1]))
                                     #print(row_ax,col_ax, " ..")
                                     #print(self.pairsArray[ii][0],self.pairsArray[ii][1])
                                     mean = numpy.sum(x_dist * y_dist) / numpy.sum(y_dist)
                                     sigma = numpy.sqrt(numpy.sum(y_dist * (x_dist - mean)**2) / numpy.sum(y_dist))
                                     parg = [numpy.amax(y_dist),mean,sigma]
                                     newY = None
                                     try :
                                         gauss_fit, covariance = curve_fit(fit_func, x_dist, y_dist,p0=parg)
                                         mode = gauss_fit[1]
                                         stdv = gauss_fit[2]
                                         #print(" FIT OK",gauss_fit)
                                         '''
                                         ###ONLY FOR TEST:
                                         if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
                                             newY = fit_func(x_dist,gauss_fit[0],gauss_fit[1],gauss_fit[2])
                                             axs[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
                                             axs[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
                                             axs[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
                                         '''
                                     except:
                                         mode = mean
                                         stdv = sigma
                                         #print("FIT FAIL")
                                         #continue
                                     #print(mode,stdv)
                                     #Removing echoes greater than mode + std_factor*stdv
                                     noval = (abs(func2clean - mode)>=(factor_stdv*stdv)).nonzero()
                                     #noval tiene los indices que se van a remover
                                     #print("Chan ",ii," novals: ",len(noval[0]))
                                     if len(noval[0]) > 0: #forma de array (N,) es igual a longitud (N)
                                         novall = ((func2clean - mode) >= (factor_stdv*stdv)).nonzero()
                                         #print(novall)
                                         #print("   ",self.pairsArray[ii])
                                         #cross_pairs = self.pairsArray[ii]
                                             #Getting coherent echoes which are removed.
                                         # if len(novall[0]) > 0:
                                         #
                                         #     val_spc[novall[0],cross_pairs[0],ifreq,ih] = 1
                                         #     val_spc[novall[0],cross_pairs[1],ifreq,ih] = 1
                                         #     val_cspc[novall[0],ii,ifreq,ih] = 1
                                             #print("OUT NOVALL 1")
                                         try:
                                             pair = (self.channels[ii],self.channels[ii + 1])
                                         except:
                                             pair = (99,99)
                                         #print("par ", pair)
                                         if ( pair in self.crosspairs):
                                             q = self.crosspairs.index(pair)
                                             #print("está aqui: ", q, (ii,ii + 1))
                                             new_a = numpy.delete(cspectra[:,q,ifreq,ih], noval[0])
                                             cspectra[noval,q,ifreq,ih] = numpy.mean(new_a) #mean CrossSpectra
                                         #if channelA_ok:
                                         #chA = self.channels.index(cross_pairs[0])
                                         new_b = numpy.delete(spectra[:,ii,ifreq,ih], noval[0])
                                         spectra[noval,ii,ifreq,ih] = numpy.mean(new_b) #mean Spectra Pair A
                                         #channelA_ok = False
                                         # chB = self.channels.index(cross_pairs[1])
                                         # new_c = numpy.delete(spectra[:,chB,ifreq,ih], noval[0])
                                         # spectra[noval,chB,ifreq,ih] = numpy.mean(new_c) #mean Spectra Pair B
                                         #
                                         # channelsProcssd.append(self.crosspairs[ii][0]) # save channel A
                                         # channelsProcssd.append(self.crosspairs[ii][1]) # save channel B
                                     '''
                                     ###ONLY FOR TEST:
                                     if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
                                         func2clean = 10*numpy.log10(numpy.absolute(spectra[:,ii,ifreq,ih]))
                                         y_dist,binstep = numpy.histogram(func2clean,bins=range(int(min_val),int(max_val+2),step))
                                         axs2[row_ax,col_ax].plot(binstep[:-1],newY,color='red')
                                         axs2[row_ax,col_ax].plot(binstep[:-1],y_dist,color='green')
                                         axs2[row_ax,col_ax].set_title("CH "+str(self.channels[ii]))
                                     '''
                         '''
                             ###ONLY FOR TEST:
                                 col_ax += 1    #contador de ploteo columnas
                                 ##print(col_ax)
                             ###ONLY FOR TEST:
                             if ifreq ==fft_rev and ih==hei_rev: #TO VIEW A SIGNLE FREQUENCY
                                 title = str(dataOut.datatime)+" nFFT: "+str(ifreq)+"   Alt: "+str(self.heights[ih])+ " km"
                                 title2 = str(dataOut.datatime)+" nFFT: "+str(ifreq)+"   Alt: "+str(self.heights[ih])+ " km  CLEANED"
                                 fig.suptitle(title)
                                 fig2.suptitle(title2)
                                 plt.show()
                         '''
                         ##################################################################################################
                     #print("Getting average of the spectra and cross-spectra from incoherent echoes.")
                     out_spectra = numpy.zeros([self.nChan,self.nFFTPoints,self.nHeights], dtype=float) #+numpy.nan
                     out_cspectra = numpy.zeros([self.nPairs,self.nFFTPoints,self.nHeights], dtype=complex) #+numpy.nan
                     for ih in range(self.nHeights):
                         for ifreq in range(self.nFFTPoints):
                             for ich in range(self.nChan):
                                 tmp = spectra[:,ich,ifreq,ih]
                                 valid = (numpy.isfinite(tmp[:])==True).nonzero()
                                 if len(valid[0]) >0 :
                                     out_spectra[ich,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
                             for icr in range(self.nPairs):
                                 tmp = numpy.squeeze(cspectra[:,icr,ifreq,ih])
                                 valid = (numpy.isfinite(tmp)==True).nonzero()
                                 if len(valid[0]) > 0:
                                     out_cspectra[icr,ifreq,ih] = numpy.nansum(tmp)#/len(valid[0])
                     return out_spectra, out_cspectra
                 def REM_ISOLATED_POINTS(self,array,rth):
             #         import matplotlib.pyplot as plt
                     if rth == None :
                         rth = 4
                     #print("REM ISO")
                     num_prof = len(array[0,:,0])
                     num_hei = len(array[0,0,:])
                     n2d = len(array[:,0,0])
                     for ii in range(n2d) :
                         #print ii,n2d
                         tmp = array[ii,:,:]
                         #print tmp.shape, array[ii,101,:],array[ii,102,:]
                       # fig = plt.figure(figsize=(6,5))
                       # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
                       # ax = fig.add_axes([left, bottom, width, height])
                       # x = range(num_prof)
                       # y = range(num_hei)
                       # cp = ax.contour(y,x,tmp)
                       # ax.clabel(cp, inline=True,fontsize=10)
                       # plt.show()
                         #indxs = WHERE(FINITE(tmp) AND tmp GT 0,cindxs)
                         tmp = numpy.reshape(tmp,num_prof*num_hei)
                         indxs1 = (numpy.isfinite(tmp)==True).nonzero()
                         indxs2 = (tmp > 0).nonzero()
                         indxs1 = (indxs1[0])
                         indxs2 = indxs2[0]
                         #indxs1 = numpy.array(indxs1[0])
                         #indxs2 = numpy.array(indxs2[0])
                         indxs = None
                         #print indxs1 , indxs2
                         for iv in range(len(indxs2)):
                             indv = numpy.array((indxs1 == indxs2[iv]).nonzero())
                             #print len(indxs2), indv
                             if len(indv[0]) > 0  :
                                 indxs =  numpy.concatenate((indxs,indxs2[iv]), axis=None)
               #           print indxs
                         indxs = indxs[1:]
                         #print(indxs, len(indxs))
                         if len(indxs) < 4 :
                             array[ii,:,:] = 0.
                             return
                         xpos = numpy.mod(indxs ,num_hei)
                         ypos = (indxs / num_hei)
                         sx = numpy.argsort(xpos) # Ordering respect to "x" (time)
                         #print sx
                         xpos = xpos[sx]
                         ypos = ypos[sx]
                  # *********************************** Cleaning isolated points **********************************
                         ic = 0
                         while True :
                             r = numpy.sqrt(list(numpy.power((xpos[ic]-xpos),2)+ numpy.power((ypos[ic]-ypos),2)))
                             #no_coh = WHERE(FINITE(r) AND (r LE rth),cno_coh)
                             #plt.plot(r)
                             #plt.show()
                             no_coh1 = (numpy.isfinite(r)==True).nonzero()
                             no_coh2 = (r <= rth).nonzero()
                             #print r, no_coh1, no_coh2
                             no_coh1 = numpy.array(no_coh1[0])
                             no_coh2 = numpy.array(no_coh2[0])
                             no_coh = None
                             #print valid1 , valid2
                             for iv in range(len(no_coh2)):
                                 indv = numpy.array((no_coh1 == no_coh2[iv]).nonzero())
                                 if len(indv[0]) > 0  :
                                     no_coh =  numpy.concatenate((no_coh,no_coh2[iv]), axis=None)
                             no_coh = no_coh[1:]
                             #print len(no_coh), no_coh
                             if len(no_coh) < 4 :
                                 #print xpos[ic], ypos[ic], ic
                 #                plt.plot(r)
                 #                plt.show()
                                 xpos[ic] = numpy.nan
                                 ypos[ic] = numpy.nan
                             ic = ic + 1
                             if  (ic == len(indxs)) :
                                 break
                             #print( xpos, ypos)
                         indxs = (numpy.isfinite(list(xpos))==True).nonzero()
                         #print indxs[0]
                         if len(indxs[0]) < 4 :
                             array[ii,:,:] = 0.
                             return
                         xpos = xpos[indxs[0]]
                         ypos = ypos[indxs[0]]
                         for i in range(0,len(ypos)):
                             ypos[i]=int(ypos[i])
                         junk = tmp
                         tmp = junk*0.0
                         tmp[list(xpos + (ypos*num_hei))] = junk[list(xpos + (ypos*num_hei))]
                         array[ii,:,:] = numpy.reshape(tmp,(num_prof,num_hei))
                         #print array.shape
                         #tmp = numpy.reshape(tmp,(num_prof,num_hei))
                         #print tmp.shape
                         # fig = plt.figure(figsize=(6,5))
                         # left, bottom, width, height = 0.1, 0.1, 0.8, 0.8
                         # ax = fig.add_axes([left, bottom, width, height])
                         # x = range(num_prof)
                         # y = range(num_hei)
                         # cp = ax.contour(y,x,array[ii,:,:])
                         # ax.clabel(cp, inline=True,fontsize=10)
                         # plt.show()
                     return array
             class IntegrationFaradaySpectra(Operation):
                 __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_ints = None #matriz de numero de integracions (CH,HEI)
                 n = None
                 minHei_ind = None
                 maxHei_ind = None
                 navg = 1.0
                 factor = 0.0
                 dataoutliers = None # (CHANNELS, HEIGHTS)
                 def __init__(self):
                     Operation.__init__(self)
                 def setup(self, dataOut,n=None, timeInterval=None, overlapping=False, DPL=None, minHei=None, maxHei=None, avg=1,factor=0.75):
                     """
                     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 = []
                     self.__buffer_cspc = []
                     self.__buffer_dc = 0
                     self.__profIndex = 0
                     self.__dataReady = False
                     self.__byTime = False
                     self.factor = factor
                     self.navg = avg
                     #self.ByLags = dataOut.ByLags   ###REDEFINIR
                     self.ByLags = False
                     self.maxProfilesInt = 1
                     if DPL != None:
                         self.DPL=DPL
                     else:
                         #self.DPL=dataOut.DPL    ###REDEFINIR
                         self.DPL=0
                     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
                     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]
                     ind_list1 = numpy.where(self.dataOut.heightList >= minHei)
                     ind_list2 = numpy.where(self.dataOut.heightList <= maxHei)
                     self.minHei_ind = ind_list1[0][0]
                     self.maxHei_ind = ind_list2[0][-1]
                 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.append(data_spc)
                     if data_cspc is None:
                         self.__buffer_cspc = None
                     else:
                         self.__buffer_cspc.append(data_cspc)
                     if data_dc is None:
                         self.__buffer_dc = None
                     else:
                         self.__buffer_dc += data_dc
                     self.__profIndex += 1
                     return
                 def hildebrand_sekhon_Integration(self,sortdata,navg, factor):
                     #data debe estar ordenado
                     #sortdata = numpy.sort(data, axis=None)
                     #sortID=data.argsort()
                     lenOfData = len(sortdata)
                     nums_min = lenOfData*factor
                     if nums_min <= 5:
                         nums_min = 5
                     sump = 0.
                     sumq = 0.
                     j = 0
                     cont = 1
                     while((cont == 1)and(j < lenOfData)):
                         sump += sortdata[j]
                         sumq += sortdata[j]**2
                         if j > nums_min:
                             rtest = float(j)/(j-1) + 1.0/navg
                             if ((sumq*j) > (rtest*sump**2)):
                                 j = j - 1
                                 sump = sump - sortdata[j]
                                 sumq = sumq - sortdata[j]**2
                                 cont = 0
                         j += 1
                     #lnoise = sump / j
                     #print("H S done")
                     #return j,sortID
                     return j
                 def pushData(self):
                     """
                     Return the sum of the last profiles and the profiles used in the sum.
                     Affected:
                     self.__profileIndex
                     """
                     bufferH=None
                     buffer=None
                     buffer1=None
                     buffer_cspc=None
                     self.__buffer_spc=numpy.array(self.__buffer_spc)
                     try:
                         self.__buffer_cspc=numpy.array(self.__buffer_cspc)
                     except :
                         #print("No cpsc",e)
                         pass
                     #print("FREQ_DC",self.__buffer_spc.shape,self.__buffer_cspc.shape)
                     freq_dc = int(self.__buffer_spc.shape[2] / 2)
                     #print("FREQ_DC",freq_dc,self.__buffer_spc.shape,self.nHeights)
                     self.dataOutliers = numpy.zeros((self.nChannels,self.nHeights)) # --> almacen de outliers
                     for k in range(self.minHei_ind,self.maxHei_ind):
                         try:
                             buffer_cspc=numpy.copy(self.__buffer_cspc[:,:,:,k])
                         except:
                             #print("No cpsc",e)
                             pass
                         outliers_IDs_cspc=[]
                         cspc_outliers_exist=False
                         for i in range(self.nChannels):#dataOut.nChannels):
                             buffer1=numpy.copy(self.__buffer_spc[:,i,:,k])
                             indexes=[]
                             #sortIDs=[]
                             outliers_IDs=[]
                             for j in range(self.nProfiles): #frecuencias en el tiempo
                                 # if i==0 and j==freq_dc: #NOT CONSIDERING DC PROFILE AT CHANNEL 0
                                 #     continue
                                 # if i==1 and j==0: #NOT CONSIDERING DC PROFILE AT CHANNEL 1
                                 #     continue
                                 buffer=buffer1[:,j]
                                 sortdata = numpy.sort(buffer, axis=None)
                                 sortID=buffer.argsort()
                                 index = _noise.hildebrand_sekhon2(sortdata,self.navg)
                                 #index,sortID=self.hildebrand_sekhon_Integration(buffer,1,self.factor)
                                 # fig,ax = plt.subplots()
                                 # ax.set_title(str(k)+" "+str(j))
                                 # x=range(len(sortdata))
                                 # ax.scatter(x,sortdata)
                                 # ax.axvline(index)
                                 # plt.show()
                                 indexes.append(index)
                                 #sortIDs.append(sortID)
                                 outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
                             #print("Outliers: ",outliers_IDs)
                             outliers_IDs=numpy.array(outliers_IDs)
                             outliers_IDs=outliers_IDs.ravel()
                             outliers_IDs=numpy.unique(outliers_IDs)
                             outliers_IDs=outliers_IDs.astype(numpy.dtype('int64'))
                             indexes=numpy.array(indexes)
                             indexmin=numpy.min(indexes)
                             #print(indexmin,buffer1.shape[0], k)
                             # fig,ax = plt.subplots()
                             # ax.plot(sortdata)
                             # ax2 = ax.twinx()
                             # x=range(len(indexes))
                             # #plt.scatter(x,indexes)
                             # ax2.scatter(x,indexes)
                             # plt.show()
                             if indexmin != buffer1.shape[0]:
                                 if self.nChannels > 1:
                                     cspc_outliers_exist= True
                                 lt=outliers_IDs
                                 #avg=numpy.mean(buffer1[[t for t in range(buffer1.shape[0]) if t not in lt],:],axis=0)
                                 for p in list(outliers_IDs):
                                     #buffer1[p,:]=avg
                                     buffer1[p,:] = numpy.NaN
                                 self.dataOutliers[i,k] = len(outliers_IDs)
                             self.__buffer_spc[:,i,:,k]=numpy.copy(buffer1)
                             outliers_IDs_cspc=numpy.append(outliers_IDs_cspc,outliers_IDs)
                         outliers_IDs_cspc=outliers_IDs_cspc.astype(numpy.dtype('int64'))
                         if cspc_outliers_exist :
                             lt=outliers_IDs_cspc
                             #avg=numpy.mean(buffer_cspc[[t for t in range(buffer_cspc.shape[0]) if t not in lt],:],axis=0)
                             for p in list(outliers_IDs_cspc):
                                 #buffer_cspc[p,:]=avg
                                 buffer_cspc[p,:] = numpy.NaN
                         try:
                             self.__buffer_cspc[:,:,:,k]=numpy.copy(buffer_cspc)
                         except:
                             #print("No cpsc",e)
                             pass
                         #else:
                             #break
                     nOutliers = len(outliers_IDs)
                     #print("Outliers  n: ",self.dataOutliers,nOutliers)
                     buffer=None
                     bufferH=None
                     buffer1=None
                     buffer_cspc=None
                     buffer=None
                     #data_spc = numpy.sum(self.__buffer_spc,axis=0)
                     data_spc = numpy.nansum(self.__buffer_spc,axis=0)
                     try:
                         #data_cspc = numpy.sum(self.__buffer_cspc,axis=0)
                         data_cspc = numpy.nansum(self.__buffer_cspc,axis=0)
                     except:
                         #print("No cpsc",e)
                         pass
                     data_dc = self.__buffer_dc
                     #(CH, HEIGH)
                     self.maxProfilesInt = self.__profIndex
                     n = self.__profIndex - self.dataOutliers # n becomes a matrix
                     self.__buffer_spc = []
                     self.__buffer_cspc = []
                     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_ints = 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_ints = 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, DPL = None,timeInterval=None, overlapping=False, minHei=None, maxHei=None, avg=1, factor=0.75):
                     self.dataOut = dataOut
                     if n == 1:
                         return self.dataOut
                     if self.dataOut.nChannels == 1:
                         self.dataOut.data_cspc = None #si es un solo canal no vale la pena acumular DATOS
                     #print(self.dataOut.data_spc.shape, self.dataOut.data_cspc)
                     if not self.isConfig:
                         self.setup(self.dataOut, n, timeInterval, overlapping,DPL ,minHei, maxHei, avg, factor)
                         self.isConfig = True
                     if not self.ByLags:
                         self.nProfiles=self.dataOut.nProfiles
                         self.nChannels=self.dataOut.nChannels
                         self.nHeights=self.dataOut.nHeights
                         avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(self.dataOut.utctime,
                                                                                             self.dataOut.data_spc,
                                                                                             self.dataOut.data_cspc,
                                                                                             self.dataOut.data_dc)
                     else:
                         self.nProfiles=self.dataOut.nProfiles
                         self.nChannels=self.dataOut.nChannels
                         self.nHeights=self.dataOut.nHeights
                         avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc = self.integrate(self.dataOut.utctime,
                                                                                             self.dataOut.dataLag_spc,
                                                                                             self.dataOut.dataLag_cspc,
                                                                                             self.dataOut.dataLag_dc)
                     self.dataOut.flagNoData = True
                     if self.__dataReady:
                         if not self.ByLags:
                             if self.nChannels == 1:
                                 #print("f int", avgdata_spc.shape)
                                 self.dataOut.data_spc = avgdata_spc
                                 self.dataOut.data_cspc = avgdata_spc
                             else:
                                 self.dataOut.data_spc = numpy.squeeze(avgdata_spc)
                                 self.dataOut.data_cspc = numpy.squeeze(avgdata_cspc)
                             self.dataOut.data_dc = avgdata_dc
                             self.dataOut.data_outlier = self.dataOutliers
                         else:
                             self.dataOut.dataLag_spc = avgdata_spc
                             self.dataOut.dataLag_cspc = avgdata_cspc
                             self.dataOut.dataLag_dc = avgdata_dc
                             self.dataOut.data_spc=self.dataOut.dataLag_spc[:,:,:,self.dataOut.LagPlot]
                             self.dataOut.data_cspc=self.dataOut.dataLag_cspc[:,:,:,self.dataOut.LagPlot]
                             self.dataOut.data_dc=self.dataOut.dataLag_dc[:,:,self.dataOut.LagPlot]
                         self.dataOut.nIncohInt *= self.n_ints
                         self.dataOut.max_nIncohInt = self.maxProfilesInt
                         #print(self.dataOut.max_nIncohInt)
                         self.dataOut.utctime = avgdatatime
                         self.dataOut.flagNoData = False
                         #print("Faraday Integration DONE...")
                     #print(self.dataOut.flagNoData)
                     return self.dataOut
             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
                 incohInt = 0
                 nOutliers = 0
                 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
                     self.incohInt = 0
                     self.nOutliers = 0
                     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
                     """
                     if data_spc.all() == numpy.nan :
                         print("nan ")
                         return
                     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
                     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
                     if dataOut.flagNoData == True:
                         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)
                     self.incohInt += dataOut.nIncohInt
                     self.nOutliers += dataOut.data_outlier
                     if self.__dataReady:
                         #print("prof: ",dataOut.max_nIncohInt,self.__profIndex)
                         dataOut.data_spc = avgdata_spc
                         dataOut.data_cspc = avgdata_cspc
                         dataOut.data_dc = avgdata_dc
                         dataOut.nIncohInt = self.incohInt
                         dataOut.data_outlier = self.nOutliers
                         dataOut.utctime = avgdatatime
                         dataOut.flagNoData = False
                         dataOut.max_nIncohInt +=  self.__profIndex
                         self.incohInt = 0
                         self.nOutliers = 0
                         self.__profIndex = 0
                         #print("IncohInt Done")
                     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

schainpy/model/proc/jroproc_voltage.py +38 -32

             import sys
             import numpy,math
             from scipy import interpolate
             from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
             from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
             from schainpy.utils import log
             from schainpy.model.io.utils import getHei_index
             from time import time
-            #import datetime
             import numpy
             #import copy
             from schainpy.model.data import _noise
             class VoltageProc(ProcessingUnit):
                 def __init__(self):
                     ProcessingUnit.__init__(self)
                     self.dataOut = Voltage()
                     self.flip = 1
                     self.setupReq = False
                 def run(self):
                     #print("running volt proc")
                     if self.dataIn.type == 'AMISR':
                         self.__updateObjFromAmisrInput()
                     if  self.dataOut.buffer_empty:
                         if self.dataIn.type == 'Voltage':
                             self.dataOut.copy(self.dataIn)
                             #print("new volts reading")
                 def __updateObjFromAmisrInput(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
                     self.dataOut.flagNoData = self.dataIn.flagNoData
                     self.dataOut.data = self.dataIn.data
                     self.dataOut.utctime = self.dataIn.utctime
                     self.dataOut.channelList = self.dataIn.channelList
                     #self.dataOut.timeInterval = self.dataIn.timeInterval
                     self.dataOut.heightList = self.dataIn.heightList
                     self.dataOut.nProfiles = self.dataIn.nProfiles
                     self.dataOut.nCohInt = self.dataIn.nCohInt
                     self.dataOut.ippSeconds = self.dataIn.ippSeconds
                     self.dataOut.frequency = self.dataIn.frequency
                     self.dataOut.azimuth = self.dataIn.azimuth
                     self.dataOut.zenith = self.dataIn.zenith
                     self.dataOut.beam.codeList = self.dataIn.beam.codeList
                     self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
                     self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
             class selectChannels(Operation):
                 def run(self, dataOut, channelList=None):
                     self.channelList = channelList
                     if self.channelList == None:
                         print("Missing channelList")
                         return dataOut
                     channelIndexList = []
                     if  type(dataOut.channelList) is not list:    #leer array desde HDF5
                         try:
                             dataOut.channelList = dataOut.channelList.tolist()
                         except Exception as e:
                             print("Select Channels: ",e)
                     for channel in self.channelList:
                         if channel not in dataOut.channelList:
                             raise ValueError("Channel %d is not in %s" %(channel, str(dataOut.channelList)))
                         index = dataOut.channelList.index(channel)
                         channelIndexList.append(index)
                     dataOut = self.selectChannelsByIndex(dataOut,channelIndexList)
                     return dataOut
                 def selectChannelsByIndex(self, dataOut, channelIndexList):
                     """
                     Selecciona un bloque de datos en base a canales segun el channelIndexList
                     Input:
                         channelIndexList    :    lista sencilla de canales a seleccionar por ej. [2,3,7]
                     Affected:
                         dataOut.data
                         dataOut.channelIndexList
                         dataOut.nChannels
                         dataOut.m_ProcessingHeader.totalSpectra
                         dataOut.systemHeaderObj.numChannels
                         dataOut.m_ProcessingHeader.blockSize
                     Return:
                         None
                     """
                     #print("selectChannelsByIndex")
                     # for channelIndex in channelIndexList:
                     #     if channelIndex not in dataOut.channelIndexList:
                     #         raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
                     if dataOut.type == 'Voltage':
                         if dataOut.flagDataAsBlock:
                             """
                             Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
                             """
                             data = dataOut.data[channelIndexList,:,:]
                         else:
                             data = dataOut.data[channelIndexList,:]
                         dataOut.data = data
                         # dataOut.channelList = [dataOut.channelList[i] for i in channelIndexList]
                         dataOut.channelList = range(len(channelIndexList))
                     elif dataOut.type == 'Spectra':
                         if hasattr(dataOut, 'data_spc'):
                             if dataOut.data_spc is None:
                                 raise ValueError("data_spc is None")
                                 return dataOut
                             else:
                                 data_spc = dataOut.data_spc[channelIndexList, :]
                                 dataOut.data_spc = data_spc
                         # if hasattr(dataOut, 'data_dc') :# and
                         #     if dataOut.data_dc is None:
                         #         raise ValueError("data_dc is None")
                         #         return dataOut
                         #     else:
                         #         data_dc = dataOut.data_dc[channelIndexList, :]
                         #         dataOut.data_dc = data_dc
                         # dataOut.channelList = [dataOut.channelList[i] for i in channelIndexList]
                         dataOut.channelList = channelIndexList
                         dataOut = self.__selectPairsByChannel(dataOut,channelIndexList)
                     return dataOut
                 def __selectPairsByChannel(self, dataOut, channelList=None):
                     #print("__selectPairsByChannel")
                     if channelList == None:
                         return
                     pairsIndexListSelected = []
                     for pairIndex in dataOut.pairsIndexList:
                         # First pair
                         if dataOut.pairsList[pairIndex][0] not in channelList:
                             continue
                         # Second pair
                         if dataOut.pairsList[pairIndex][1] not in channelList:
                             continue
                         pairsIndexListSelected.append(pairIndex)
                     if not pairsIndexListSelected:
                         dataOut.data_cspc = None
                         dataOut.pairsList = []
                         return
                     dataOut.data_cspc = dataOut.data_cspc[pairsIndexListSelected]
                     dataOut.pairsList = [dataOut.pairsList[i]
                                               for i in pairsIndexListSelected]
                     return dataOut
             class selectHeights(Operation):
                 def run(self, dataOut, minHei=None, maxHei=None, minIndex=None, maxIndex=None):
                     """
                     Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango
                     minHei <= height <= maxHei
                     Input:
                         minHei    :    valor minimo de altura a considerar
                         maxHei    :    valor maximo de altura a considerar
                     Affected:
                         Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
                     Return:
 si el metodo se ejecuto con exito caso contrario devuelve 0
                     """
                     self.dataOut = dataOut
                     if minHei and maxHei:
                         if (minHei < dataOut.heightList[0]):
                             minHei = dataOut.heightList[0]
                         if (maxHei > dataOut.heightList[-1]):
                             maxHei = dataOut.heightList[-1]
                         minIndex = 0
                         maxIndex = 0
                         heights = 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)
                     self.selectHeightsByIndex(minIndex, maxIndex)
                     return dataOut
                 def selectHeightsByIndex(self, minIndex, maxIndex):
                     """
                     Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
                     minIndex <= index <= maxIndex
                     Input:
                         minIndex    :    valor de indice minimo de altura a considerar
                         maxIndex    :    valor de indice maximo de altura a considerar
                     Affected:
                         self.dataOut.data
                         self.dataOut.heightList
                     Return:
 si el metodo se ejecuto con exito caso contrario devuelve 0
                     """
                     if self.dataOut.type == 'Voltage':
                         if (minIndex < 0) or (minIndex > maxIndex):
                             raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
                         if (maxIndex >= self.dataOut.nHeights):
                             maxIndex = self.dataOut.nHeights
                         #voltage
                         if self.dataOut.flagDataAsBlock:
                             """
                             Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
                             """
                             data = self.dataOut.data[:,:, minIndex:maxIndex]
                         else:
                             data = self.dataOut.data[:, minIndex:maxIndex]
                         #         firstHeight = self.dataOut.heightList[minIndex]
                         self.dataOut.data = data
                         self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
                         if self.dataOut.nHeights <= 1:
                             raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
                     elif self.dataOut.type == 'Spectra':
                         if (minIndex < 0) or (minIndex > maxIndex):
                             raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
                                 minIndex, maxIndex))
                         if (maxIndex >= self.dataOut.nHeights):
                             maxIndex = self.dataOut.nHeights - 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.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
                     return 1
             class filterByHeights(Operation):
                 def run(self, dataOut, window):
                     deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                     if window == None:
                         window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
                     newdelta = deltaHeight * window
                     r = dataOut.nHeights % window
                     newheights = (dataOut.nHeights-r)/window
                     if newheights <= 1:
                         raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
                     if dataOut.flagDataAsBlock:
                         """
                         Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
                         """
                         buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)]
                         buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window)
                         buffer = numpy.sum(buffer,3)
                     else:
                         buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
                         buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
                         buffer = numpy.sum(buffer,2)
                     dataOut.data = buffer
                     dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
                     dataOut.windowOfFilter = window
                     return dataOut
             class setH0(Operation):
                 def run(self, dataOut, h0, deltaHeight = None):
                     if not deltaHeight:
                         deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                     nHeights = dataOut.nHeights
                     newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
                     dataOut.heightList = newHeiRange
                     return dataOut
             class deFlip(Operation):
                 def run(self, dataOut, channelList = []):
                     data = dataOut.data.copy()
                     if dataOut.flagDataAsBlock:
                         flip = self.flip
                         profileList = list(range(dataOut.nProfiles))
                         if not channelList:
                             for thisProfile in profileList:
                                 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
                                 flip *= -1.0
                         else:
                             for thisChannel in channelList:
                                 if thisChannel not in dataOut.channelList:
                                     continue
                                 for thisProfile in profileList:
                                     data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
                                     flip *= -1.0
                         self.flip = flip
                     else:
                         if not channelList:
                             data[:,:] = data[:,:]*self.flip
                         else:
                             for thisChannel in channelList:
                                 if thisChannel not in dataOut.channelList:
                                     continue
                                 data[thisChannel,:] = data[thisChannel,:]*self.flip
                         self.flip *= -1.
                     dataOut.data = data
                     return dataOut
             class setAttribute(Operation):
                 '''
                 Set an arbitrary attribute(s) to dataOut
                 '''
                 def __init__(self):
                     Operation.__init__(self)
                     self._ready = False
                 def run(self, dataOut, **kwargs):
                     for key, value in kwargs.items():
                         setattr(dataOut, key, value)
                     return dataOut
             @MPDecorator
             class printAttribute(Operation):
                 '''
                 Print an arbitrary attribute of dataOut
                 '''
                 def __init__(self):
                     Operation.__init__(self)
                 def run(self, dataOut, attributes):
                     if isinstance(attributes, str):
                         attributes = [attributes]
                     for attr in attributes:
                         if hasattr(dataOut, attr):
                             log.log(getattr(dataOut, attr), attr)
             class interpolateHeights(Operation):
                 def run(self, dataOut, topLim, botLim):
                     #69 al 72 para julia
                     #82-84 para meteoros
                     if len(numpy.shape(dataOut.data))==2:
                         sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
                         sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
                         #dataOut.data[:,botLim:limSup+1] = sampInterp
                         dataOut.data[:,botLim:topLim+1] = sampInterp
                     else:
                         nHeights = dataOut.data.shape[2]
                         x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
                         y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
                         f = interpolate.interp1d(x, y, axis = 2)
                         xnew = numpy.arange(botLim,topLim+1)
                         ynew = f(xnew)
                         dataOut.data[:,:,botLim:topLim+1]  = ynew
                     return dataOut
             class CohInt(Operation):
                 isConfig = False
                 __profIndex = 0
                 __byTime = False
                 __initime = None
                 __lastdatatime = None
                 __integrationtime = None
                 __buffer = None
                 __bufferStride = []
                 __dataReady = False
                 __profIndexStride = 0
                 __dataToPutStride = False
                 n = None
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=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 = None
                     self.__dataReady = False
                     self.byblock = byblock
                     self.stride = stride
                     if n == None and timeInterval == None:
                         raise ValueError("n or timeInterval should be specified ...")
                     if n != None:
                         self.n = n
                         self.__byTime = False
                     else:
                         self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
                         self.n = 9999
                         self.__byTime = True
                     if overlapping:
                         self.__withOverlapping = True
                         self.__buffer = None
                     else:
                         self.__withOverlapping = False
                         self.__buffer = 0
                     self.__profIndex = 0
                 def putData(self, data):
                     """
                     Add a profile to the __buffer and increase in one the __profileIndex
                     """
                     if not self.__withOverlapping:
                         self.__buffer += data.copy()
                         self.__profIndex += 1
                         return
                     #Overlapping data
                     nChannels, nHeis = data.shape
                     data = numpy.reshape(data, (1, nChannels, nHeis))
                     #If the buffer is empty then it takes the data value
                     if self.__buffer is None:
                         self.__buffer = data
                         self.__profIndex += 1
                         return
                     #If the buffer length is lower than n then stakcing the data value
                     if self.__profIndex < self.n:
                         self.__buffer = numpy.vstack((self.__buffer, data))
                         self.__profIndex += 1
                         return
                     #If the buffer length is equal to n then replacing the last buffer value with the data value
                     self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
                     self.__buffer[self.n-1] = data
                     self.__profIndex = self.n
                     return
                 def pushData(self):
                     """
                     Return the sum of the last profiles and the profiles used in the sum.
                     Affected:
                     self.__profileIndex
                     """
                     if not self.__withOverlapping:
                         data = self.__buffer
                         n = self.__profIndex
                         self.__buffer = 0
                         self.__profIndex = 0
                         return data, n
                     #Integration with Overlapping
                     data = numpy.sum(self.__buffer, axis=0)
                     # print data
                     # raise
                     n = self.__profIndex
                     return data, n
                 def byProfiles(self, data):
                     self.__dataReady = False
                     avgdata = None
                     #         n = None
                     # print data
                     # raise
                     self.putData(data)
                     if self.__profIndex == self.n:
                         avgdata, n = self.pushData()
                         self.__dataReady = True
                     return avgdata
                 def byTime(self, data, datatime):
                     self.__dataReady = False
                     avgdata = None
                     n = None
                     self.putData(data)
                     if (datatime - self.__initime) >= self.__integrationtime:
                         avgdata, n = self.pushData()
                         self.n = n
                         self.__dataReady = True
                     return avgdata
                 def integrateByStride(self, data, datatime):
                     # print data
                     if self.__profIndex == 0:
                         self.__buffer = [[data.copy(), datatime]]
                     else:
                         self.__buffer.append([data.copy(),datatime])
                     self.__profIndex += 1
                     self.__dataReady = False
                     if self.__profIndex == self.n * self.stride :
                         self.__dataToPutStride = True
                         self.__profIndexStride = 0
                         self.__profIndex = 0
                         self.__bufferStride = []
                         for i in range(self.stride):
                             current = self.__buffer[i::self.stride]
                             data = numpy.sum([t[0] for t in current], axis=0)
                             avgdatatime = numpy.average([t[1] for t in current])
                             # print data
                             self.__bufferStride.append((data, avgdatatime))
                     if self.__dataToPutStride:
                         self.__dataReady = True
                         self.__profIndexStride += 1
                         if self.__profIndexStride == self.stride:
                             self.__dataToPutStride = False
                         # print self.__bufferStride[self.__profIndexStride - 1]
                         # raise
                         return self.__bufferStride[self.__profIndexStride - 1]
                     return None, None
                 def integrate(self, data, datatime=None):
                     if self.__initime == None:
                         self.__initime = datatime
                     if self.__byTime:
                         avgdata = self.byTime(data, datatime)
                     else:
                         avgdata = self.byProfiles(data)
                     self.__lastdatatime = datatime
                     if avgdata is None:
                         return None, None
                     avgdatatime = self.__initime
                     deltatime = datatime - self.__lastdatatime
                     if not self.__withOverlapping:
                         self.__initime = datatime
                     else:
                         self.__initime += deltatime
                     return avgdata, avgdatatime
                 def integrateByBlock(self, dataOut):
                     times = int(dataOut.data.shape[1]/self.n)
                     avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
                     id_min = 0
                     id_max = self.n
                     for i in range(times):
                         junk = dataOut.data[:,id_min:id_max,:]
                         avgdata[:,i,:] = junk.sum(axis=1)
                         id_min += self.n
                         id_max += self.n
                     timeInterval = dataOut.ippSeconds*self.n
                     avgdatatime = (times - 1) * timeInterval + dataOut.utctime
                     self.__dataReady = True
                     return avgdata, avgdatatime
                 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
                     if not self.isConfig:
                         self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
                         self.isConfig = True
                     if dataOut.flagDataAsBlock:
                         """
                         Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
                         """
                         avgdata, avgdatatime = self.integrateByBlock(dataOut)
                         dataOut.nProfiles /= self.n
                     else:
                         if stride is None:
                             avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
                         else:
                             avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime)
                     #   dataOut.timeInterval *= n
                     dataOut.flagNoData = True
                     if self.__dataReady:
                         dataOut.data = avgdata
                         if not dataOut.flagCohInt:
                             dataOut.nCohInt *= self.n
                             dataOut.flagCohInt = True
                         dataOut.utctime = avgdatatime
                         # print avgdata, avgdatatime
                         # raise
                         #   dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
                         dataOut.flagNoData = False
                     return dataOut
             class Decoder(Operation):
                 isConfig = False
                 __profIndex = 0
                 code = None
                 nCode = None
                 nBaud = None
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.times = None
                     self.osamp = None
                 #         self.__setValues = False
                     self.isConfig = False
                     self.setupReq = False
                 def setup(self, code, osamp, dataOut):
                     self.__profIndex = 0
                     self.code = code
                     self.nCode = len(code)
                     self.nBaud = len(code[0])
                     if (osamp != None) and (osamp >1):
                         self.osamp = osamp
                         self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
                         self.nBaud = self.nBaud*self.osamp
                     self.__nChannels = dataOut.nChannels
                     self.__nProfiles = dataOut.nProfiles
                     self.__nHeis = dataOut.nHeights
                     if self.__nHeis < self.nBaud:
                         raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
                     #Frequency
                     __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
                     __codeBuffer[:,0:self.nBaud] = self.code
                     self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
                     if dataOut.flagDataAsBlock:
                         self.ndatadec = self.__nHeis #- self.nBaud + 1
                         self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
                     else:
                         #Time
                         self.ndatadec = self.__nHeis #- self.nBaud + 1
                         self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
                 def __convolutionInFreq(self, data):
                     fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
                     fft_data = numpy.fft.fft(data, axis=1)
                     conv = fft_data*fft_code
                     data = numpy.fft.ifft(conv,axis=1)
                     return data
                 def __convolutionInFreqOpt(self, data):
                     raise NotImplementedError
                 def __convolutionInTime(self, data):
                     code = self.code[self.__profIndex]
                     for i in range(self.__nChannels):
                         self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
                     return self.datadecTime
                 def __convolutionByBlockInTime(self, data):
                     repetitions = int(self.__nProfiles / self.nCode)
                     junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
                     junk = junk.flatten()
                     code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
                     profilesList = range(self.__nProfiles)
                     for i in range(self.__nChannels):
                         for j in profilesList:
                             self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:]
                     return self.datadecTime
                 def __convolutionByBlockInFreq(self, data):
                     raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
                     fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
                     fft_data = numpy.fft.fft(data, axis=2)
                     conv = fft_data*fft_code
                     data = numpy.fft.ifft(conv,axis=2)
                     return data
                 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
                     if dataOut.flagDecodeData:
                         print("This data is already decoded, recoding again ...")
                     if not self.isConfig:
                         if code is None:
                             if dataOut.code is None:
                                 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
                             code = dataOut.code
                         else:
                             code = numpy.array(code).reshape(nCode,nBaud)
                         self.setup(code, osamp, dataOut)
                         self.isConfig = True
                         if mode == 3:
                             sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
                         if times != None:
                             sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
                     if self.code is None:
                         print("Fail decoding: Code is not defined.")
                         return
                     self.__nProfiles = dataOut.nProfiles
                     datadec = None
                     if mode == 3:
                         mode = 0
                     if dataOut.flagDataAsBlock:
                         """
                         Decoding when data have been read as block,
                         """
                         if mode == 0:
                             datadec = self.__convolutionByBlockInTime(dataOut.data)
                         if mode == 1:
                             datadec = self.__convolutionByBlockInFreq(dataOut.data)
                     else:
                         """
                         Decoding when data have been read profile by profile
                         """
                         if mode == 0:
                             datadec = self.__convolutionInTime(dataOut.data)
                         if mode == 1:
                             datadec = self.__convolutionInFreq(dataOut.data)
                         if mode == 2:
                             datadec = self.__convolutionInFreqOpt(dataOut.data)
                     if datadec is None:
                         raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
                     dataOut.code = self.code
                     dataOut.nCode = self.nCode
                     dataOut.nBaud = self.nBaud
                     dataOut.data = datadec
                     dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
                     dataOut.flagDecodeData = True #asumo q la data esta decodificada
                     if self.__profIndex == self.nCode-1:
                         self.__profIndex = 0
                         return dataOut
                     self.__profIndex += 1
                     return dataOut
                 #        dataOut.flagDeflipData = True #asumo q la data no esta sin flip
             class ProfileConcat(Operation):
                 isConfig = False
                 buffer = None
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.profileIndex = 0
                 def reset(self):
                     self.buffer = numpy.zeros_like(self.buffer)
                     self.start_index = 0
                     self.times = 1
                 def setup(self, data, m, n=1):
                     self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
                     self.nHeights = data.shape[1]#.nHeights
                     self.start_index = 0
                     self.times = 1
                 def concat(self, data):
                     self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
                     self.start_index = self.start_index + self.nHeights
                 def run(self, dataOut, m):
                     dataOut.flagNoData = True
                     if not self.isConfig:
                         self.setup(dataOut.data, m, 1)
                         self.isConfig = True
                     if dataOut.flagDataAsBlock:
                         raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
                     else:
                         self.concat(dataOut.data)
                         self.times += 1
                         if self.times > m:
                             dataOut.data = self.buffer
                             self.reset()
                             dataOut.flagNoData = False
                             # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
                             deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                             xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
                             dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
                             dataOut.ippSeconds *= m
                     return dataOut
             class ProfileSelector(Operation):
                 profileIndex = None
                 # Tamanho total de los perfiles
                 nProfiles = None
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.profileIndex = 0
                 def incProfileIndex(self):
                     self.profileIndex += 1
                     if self.profileIndex >= self.nProfiles:
                         self.profileIndex = 0
                 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
                     if profileIndex < minIndex:
                         return False
                     if profileIndex > maxIndex:
                         return False
                     return True
                 def isThisProfileInList(self, profileIndex, profileList):
                     if profileIndex not in profileList:
                         return False
                     return True
                 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
                     """
                     ProfileSelector:
                     Inputs:
                         profileList        :    Index of profiles selected. Example: profileList = (0,1,2,7,8)
                         profileRangeList    :    Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
                         rangeList            :    List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
                     """
                     if rangeList is not None:
                         if type(rangeList[0]) not in (tuple, list):
                             rangeList = [rangeList]
                     dataOut.flagNoData = True
                     if dataOut.flagDataAsBlock:
                         """
                         data dimension  = [nChannels, nProfiles, nHeis]
                         """
                         if profileList != None:
                             dataOut.data = dataOut.data[:,profileList,:]
                         if profileRangeList != None:
                             minIndex = profileRangeList[0]
                             maxIndex = profileRangeList[1]
                             profileList = list(range(minIndex, maxIndex+1))
                             dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
                         if rangeList != None:
                             profileList = []
                             for thisRange in rangeList:
                                 minIndex = thisRange[0]
                                 maxIndex = thisRange[1]
                                 profileList.extend(list(range(minIndex, maxIndex+1)))
                             dataOut.data = dataOut.data[:,profileList,:]
                         dataOut.nProfiles = len(profileList)
                         dataOut.profileIndex = dataOut.nProfiles - 1
                         dataOut.flagNoData = False
                         return dataOut
                     """
                     data dimension  = [nChannels, nHeis]
                     """
                     if profileList != None:
                         if self.isThisProfileInList(dataOut.profileIndex, profileList):
                             self.nProfiles = len(profileList)
                             dataOut.nProfiles = self.nProfiles
                             dataOut.profileIndex = self.profileIndex
                             dataOut.flagNoData = False
                             self.incProfileIndex()
                         return dataOut
                     if profileRangeList != None:
                         minIndex = profileRangeList[0]
                         maxIndex = profileRangeList[1]
                         if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
                             self.nProfiles = maxIndex - minIndex + 1
                             dataOut.nProfiles = self.nProfiles
                             dataOut.profileIndex = self.profileIndex
                             dataOut.flagNoData = False
                             self.incProfileIndex()
                         return dataOut
                     if rangeList != None:
                         nProfiles = 0
                         for thisRange in rangeList:
                             minIndex = thisRange[0]
                             maxIndex = thisRange[1]
                             nProfiles += maxIndex - minIndex + 1
                         for thisRange in rangeList:
                             minIndex = thisRange[0]
                             maxIndex = thisRange[1]
                             if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
                                 self.nProfiles = nProfiles
                                 dataOut.nProfiles = self.nProfiles
                                 dataOut.profileIndex = self.profileIndex
                                 dataOut.flagNoData = False
                                 self.incProfileIndex()
                                 break
                         return dataOut
                     if beam != None: #beam is only for AMISR data
                         if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
                             dataOut.flagNoData = False
                             dataOut.profileIndex = self.profileIndex
                             self.incProfileIndex()
                         return dataOut
                     raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
             class Reshaper(Operation):
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.__buffer = None
                     self.__nitems = 0
                 def __appendProfile(self, dataOut, nTxs):
                     if self.__buffer is None:
                         shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
                         self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
                     ini = dataOut.nHeights * self.__nitems
                     end = ini + dataOut.nHeights
                     self.__buffer[:, ini:end] = dataOut.data
                     self.__nitems += 1
                     return int(self.__nitems*nTxs)
                 def __getBuffer(self):
                     if self.__nitems == int(1./self.__nTxs):
                         self.__nitems = 0
                         return self.__buffer.copy()
                     return None
                 def __checkInputs(self, dataOut, shape, nTxs):
                     if shape is None and nTxs is None:
                         raise ValueError("Reshaper: shape of factor should be defined")
                     if nTxs:
                         if nTxs < 0:
                             raise ValueError("nTxs should be greater than 0")
                         if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
                             raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
                         shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
                         return shape, nTxs
                     if len(shape) != 2 and len(shape) !=  3:
                         raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
                     if len(shape) == 2:
                         shape_tuple = [dataOut.nChannels]
                         shape_tuple.extend(shape)
                     else:
                         shape_tuple = list(shape)
                     nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
                     return shape_tuple, nTxs
                 def run(self, dataOut, shape=None, nTxs=None):
                     shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
                     dataOut.flagNoData = True
                     profileIndex = None
                     if dataOut.flagDataAsBlock:
                         dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
                         dataOut.flagNoData = False
                         profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
                     else:
                         if self.__nTxs < 1:
                             self.__appendProfile(dataOut, self.__nTxs)
                             new_data = self.__getBuffer()
                             if new_data is not None:
                                 dataOut.data = new_data
                                 dataOut.flagNoData = False
                                 profileIndex = dataOut.profileIndex*nTxs
                         else:
                             raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
                     deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                     dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
                     dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
                     dataOut.profileIndex = profileIndex
                     dataOut.ippSeconds /= self.__nTxs
                     return dataOut
             class SplitProfiles(Operation):
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                 def run(self, dataOut, n):
                     dataOut.flagNoData = True
                     profileIndex = None
                     if dataOut.flagDataAsBlock:
                         #nchannels, nprofiles, nsamples
                         shape = dataOut.data.shape
                         if shape[2] % n != 0:
                             raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2]))
                         new_shape = shape[0], shape[1]*n, int(shape[2]/n)
                         dataOut.data = numpy.reshape(dataOut.data, new_shape)
                         dataOut.flagNoData = False
                         profileIndex = int(dataOut.nProfiles/n) - 1
                     else:
                         raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
                     deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                     dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
                     dataOut.nProfiles = int(dataOut.nProfiles*n)
                     dataOut.profileIndex = profileIndex
                     dataOut.ippSeconds /= n
                     return dataOut
             class CombineProfiles(Operation):
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.__remData = None
                     self.__profileIndex = 0
                 def run(self, dataOut, n):
                     dataOut.flagNoData = True
                     profileIndex = None
                     if dataOut.flagDataAsBlock:
                         #nchannels, nprofiles, nsamples
                         shape = dataOut.data.shape
                         new_shape = shape[0], shape[1]/n, shape[2]*n
                         if shape[1] % n != 0:
                             raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
                         dataOut.data = numpy.reshape(dataOut.data, new_shape)
                         dataOut.flagNoData = False
                         profileIndex = int(dataOut.nProfiles*n) - 1
                     else:
                         #nchannels, nsamples
                         if self.__remData is None:
                             newData = dataOut.data
                         else:
                             newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
                         self.__profileIndex += 1
                         if self.__profileIndex < n:
                             self.__remData = newData
                             #continue
                             return
                         self.__profileIndex = 0
                         self.__remData = None
                         dataOut.data = newData
                         dataOut.flagNoData = False
                         profileIndex = dataOut.profileIndex/n
                     deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                     dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
                     dataOut.nProfiles = int(dataOut.nProfiles/n)
                     dataOut.profileIndex = profileIndex
                     dataOut.ippSeconds *= n
                     return dataOut
             class PulsePairVoltage(Operation):
                 '''
                 Function PulsePair(Signal Power, Velocity)
                 The real component of Lag[0] provides Intensity Information
                 The imag component of Lag[1] Phase provides Velocity Information
                 Configuration Parameters:
                 nPRF = Number of Several PRF
                 theta = Degree Azimuth angel Boundaries
                 Input:
                       self.dataOut
                       lag[N]
                 Affected:
                       self.dataOut.spc
                 '''
                 isConfig       = False
                 __profIndex    = 0
                 __initime      = None
                 __lastdatatime = None
                 __buffer       = None
                 noise          = None
                 __dataReady    = False
                 n              = None
                 __nch          = 0
                 __nHeis        = 0
                 removeDC       = False
                 ipp            = None
                 lambda_        = 0
                 def __init__(self,**kwargs):
                     Operation.__init__(self,**kwargs)
                 def setup(self, dataOut, n = None, removeDC=False):
                     '''
                     n= Numero de PRF's de entrada
                     '''
                     self.__initime        = None
                     self.__lastdatatime   = 0
                     self.__dataReady      = False
                     self.__buffer         = 0
                     self.__profIndex      = 0
                     self.noise            = None
                     self.__nch            = dataOut.nChannels
                     self.__nHeis          = dataOut.nHeights
                     self.removeDC         = removeDC
                     self.lambda_          = 3.0e8/(9345.0e6)
                     self.ippSec           = dataOut.ippSeconds
                     self.nCohInt          = dataOut.nCohInt
                     if n == None:
                         raise ValueError("n should be specified.")
                     if n != None:
                         if n<2:
                             raise ValueError("n should be greater than 2")
                     self.n       = n
                     self.__nProf = n
                     self.__buffer = numpy.zeros((dataOut.nChannels,
                                                        n,
                                                        dataOut.nHeights),
                                                       dtype='complex')
                 def putData(self,data):
                     '''
                     Add a profile to he __buffer and increase in one the __profiel Index
                     '''
                     self.__buffer[:,self.__profIndex,:]= data
                     self.__profIndex      += 1
                     return
                 def pushData(self,dataOut):
                     '''
                     Return the PULSEPAIR and the profiles used in the operation
                     Affected :  self.__profileIndex
                     '''
                     #----------------- Remove DC-----------------------------------
                     if self.removeDC==True:
                         mean    = numpy.mean(self.__buffer,1)
                         tmp     = mean.reshape(self.__nch,1,self.__nHeis)
                         dc= numpy.tile(tmp,[1,self.__nProf,1])
                         self.__buffer = self.__buffer -  dc
                     #------------------Calculo de Potencia ------------------------
                     pair0       = self.__buffer*numpy.conj(self.__buffer)
                     pair0       = pair0.real
                     lag_0       = numpy.sum(pair0,1)
                     #------------------Calculo de Ruido x canal--------------------
                     self.noise  = numpy.zeros(self.__nch)
                     for i in range(self.__nch):
                         daux         = numpy.sort(pair0[i,:,:],axis= None)
                         self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
                     self.noise       = self.noise.reshape(self.__nch,1)
                     self.noise       = numpy.tile(self.noise,[1,self.__nHeis])
                     noise_buffer     = self.noise.reshape(self.__nch,1,self.__nHeis)
                     noise_buffer     = numpy.tile(noise_buffer,[1,self.__nProf,1])
                     #------------------ Potencia recibida= P , Potencia senal = S , Ruido= N--
                     #------------------   P= S+N  ,P=lag_0/N ---------------------------------
                     #-------------------- Power --------------------------------------------------
                     data_power       = lag_0/(self.n*self.nCohInt)
                     #------------------  Senal  ---------------------------------------------------
                     data_intensity   = pair0 - noise_buffer
                     data_intensity   = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
                     #data_intensity   = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
                     for i in range(self.__nch):
                         for j in range(self.__nHeis):
                             if data_intensity[i][j]  < 0:
                                 data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
                     #----------------- Calculo de Frecuencia y Velocidad doppler--------
                     pair1            = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
                     lag_1            = numpy.sum(pair1,1)
                     data_freq        = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
                     data_velocity    = (self.lambda_/2.0)*data_freq
                     #---------------- Potencia promedio estimada de la Senal-----------
                     lag_0            = lag_0/self.n
                     S                = lag_0-self.noise
                     #---------------- Frecuencia Doppler promedio ---------------------
                     lag_1            = lag_1/(self.n-1)
                     R1               = numpy.abs(lag_1)
                     #---------------- Calculo del SNR----------------------------------
                     data_snrPP       = S/self.noise
                     for i in range(self.__nch):
                         for j in range(self.__nHeis):
                             if data_snrPP[i][j]  < 1.e-20:
                                 data_snrPP[i][j] = 1.e-20
                     #----------------- Calculo del ancho espectral ----------------------
                     L                = S/R1
                     L                = numpy.where(L<0,1,L)
                     L                = numpy.log(L)
                     tmp              = numpy.sqrt(numpy.absolute(L))
                     data_specwidth   = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
                     n                = self.__profIndex
                     self.__buffer    = numpy.zeros((self.__nch, self.__nProf,self.__nHeis),  dtype='complex')
                     self.__profIndex = 0
                     return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n
                 def pulsePairbyProfiles(self,dataOut):
                     self.__dataReady     =  False
                     data_power           =  None
                     data_intensity       =  None
                     data_velocity        =  None
                     data_specwidth       =  None
                     data_snrPP           =  None
                     self.putData(data=dataOut.data)
                     if self.__profIndex  == self.n:
                         data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n   = self.pushData(dataOut=dataOut)
                         self.__dataReady                   = True
                     return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth
                 def pulsePairOp(self, dataOut, datatime= None):
                     if self.__initime == None:
                         self.__initime = datatime
                     data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
                     self.__lastdatatime           = datatime
                     if data_power is None:
                         return None, None, None,None,None,None
                     avgdatatime    = self.__initime
                     deltatime      = datatime - self.__lastdatatime
                     self.__initime = datatime
                     return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime
                 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
                     if not self.isConfig:
                         self.setup(dataOut = dataOut, n    = n , removeDC=removeDC , **kwargs)
                         self.isConfig   = True
                     data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
                     dataOut.flagNoData                         = True
                     if self.__dataReady:
                         dataOut.nCohInt        *= self.n
                         dataOut.dataPP_POW      = data_intensity # S
                         dataOut.dataPP_POWER    = data_power     # P
                         dataOut.dataPP_DOP      = data_velocity
                         dataOut.dataPP_SNR      = data_snrPP
                         dataOut.dataPP_WIDTH    = data_specwidth
                         dataOut.PRFbyAngle      = self.n         #numero de PRF*cada angulo rotado que equivale a un tiempo.
                         dataOut.utctime         = avgdatatime
                         dataOut.flagNoData      = False
                     return dataOut
             # import collections
             # from scipy.stats import mode
             #
             # class Synchronize(Operation):
             #
             #     isConfig = False
             #     __profIndex = 0
             #
             #     def __init__(self, **kwargs):
             #
             #         Operation.__init__(self, **kwargs)
             # #         self.isConfig = False
             #         self.__powBuffer = None
             #         self.__startIndex = 0
             #         self.__pulseFound = False
             #
             #     def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
             #
             #         #Read data
             #
             #         powerdB = dataOut.getPower(channel = channel)
             #         noisedB = dataOut.getNoise(channel = channel)[0]
             #
             #         self.__powBuffer.extend(powerdB.flatten())
             #
             #         dataArray = numpy.array(self.__powBuffer)
             #
             #         filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
             #
             #         maxValue = numpy.nanmax(filteredPower)
             #
             #         if maxValue < noisedB + 10:
             #             #No se encuentra ningun pulso de transmision
             #             return None
             #
             #         maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
             #
             #         if len(maxValuesIndex) < 2:
             #             #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
             #             return None
             #
             #         phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
             #
             #         #Seleccionar solo valores con un espaciamiento de nSamples
             #         pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
             #
             #         if len(pulseIndex) < 2:
             #             #Solo se encontro un pulso de transmision con ancho mayor a 1
             #             return None
             #
             #         spacing = pulseIndex[1:] - pulseIndex[:-1]
             #
             #         #remover senales que se distancien menos de 10 unidades o muestras
             #         #(No deberian existir IPP menor a 10 unidades)
             #
             #         realIndex = numpy.where(spacing > 10 )[0]
             #
             #         if len(realIndex) < 2:
             #             #Solo se encontro un pulso de transmision con ancho mayor a 1
             #             return None
             #
             #         #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
             #         realPulseIndex = pulseIndex[realIndex]
             #
             #         period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
             #
             #         print "IPP = %d samples" %period
             #
             #         self.__newNSamples = dataOut.nHeights #int(period)
             #         self.__startIndex = int(realPulseIndex[0])
             #
             #         return 1
             #
             #
             #     def setup(self, nSamples, nChannels, buffer_size = 4):
             #
             #         self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
             #                                           maxlen = buffer_size*nSamples)
             #
             #         bufferList = []
             #
             #         for i in range(nChannels):
             #             bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) +  numpy.NAN,
             #                                           maxlen = buffer_size*nSamples)
             #
             #             bufferList.append(bufferByChannel)
             #
             #         self.__nSamples = nSamples
             #         self.__nChannels = nChannels
             #         self.__bufferList = bufferList
             #
             #     def run(self, dataOut, channel = 0):
             #
             #         if not self.isConfig:
             #             nSamples = dataOut.nHeights
             #             nChannels = dataOut.nChannels
             #             self.setup(nSamples, nChannels)
             #             self.isConfig = True
             #
             #         #Append new data to internal buffer
             #         for thisChannel in range(self.__nChannels):
             #             bufferByChannel = self.__bufferList[thisChannel]
             #             bufferByChannel.extend(dataOut.data[thisChannel])
             #
             #         if self.__pulseFound:
             #             self.__startIndex -= self.__nSamples
             #
             #         #Finding Tx Pulse
             #         if not self.__pulseFound:
             #             indexFound = self.__findTxPulse(dataOut, channel)
             #
             #             if indexFound == None:
             #                 dataOut.flagNoData = True
             #                 return
             #
             #             self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
             #             self.__pulseFound = True
             #             self.__startIndex = indexFound
             #
             #         #If pulse was found ...
             #         for thisChannel in range(self.__nChannels):
             #             bufferByChannel = self.__bufferList[thisChannel]
             #             #print self.__startIndex
             #             x = numpy.array(bufferByChannel)
             #             self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
             #
             #         deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
             #         dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
             # #             dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
             #
             #         dataOut.data = self.__arrayBuffer
             #
             #         self.__startIndex += self.__newNSamples
             #
             #         return
             class SSheightProfiles(Operation):
                 step          = None
                 nsamples      = None
                 bufferShape   = None
                 profileShape  = None
                 sshProfiles   = None
                 profileIndex  = None
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.isConfig = False
                 def setup(self,dataOut ,step = None , nsamples = None):
                     if step == None and nsamples == None:
                         raise ValueError("step or nheights should be specified ...")
                     self.step         = step
                     self.nsamples     = nsamples
                     self.__nChannels  = dataOut.nChannels
                     self.__nProfiles  = dataOut.nProfiles
                     self.__nHeis      = dataOut.nHeights
                     shape             = dataOut.data.shape #nchannels, nprofiles, nsamples
                     residue     =  (shape[1] - self.nsamples) % self.step
                     if residue != 0:
                         print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,shape[1] - self.nsamples,residue))
                     deltaHeight      =  dataOut.heightList[1] - dataOut.heightList[0]
                     numberProfile    =  self.nsamples
                     numberSamples    =  (shape[1] - self.nsamples)/self.step
                     self.bufferShape  = int(shape[0]), int(numberSamples), int(numberProfile)  # nchannels, nsamples , nprofiles
                     self.profileShape = int(shape[0]), int(numberProfile), int(numberSamples)  # nchannels, nprofiles, nsamples
                     self.buffer       = numpy.zeros(self.bufferShape , dtype=numpy.complex)
                     self.sshProfiles  = numpy.zeros(self.profileShape, dtype=numpy.complex)
                 def run(self, dataOut, step, nsamples, code = None, repeat = None):
                     dataOut.flagNoData      = True
                     profileIndex            = None
                     #print("nProfiles, nHeights ",dataOut.nProfiles, dataOut.nHeights)
                     #print(dataOut.getFreqRange(1)/1000.)
                     #exit(1)
                     if dataOut.flagDataAsBlock:
                         dataOut.data = numpy.average(dataOut.data,axis=1)
                         #print("jee")
                     dataOut.flagDataAsBlock = False
                     if not self.isConfig:
                         self.setup(dataOut, step=step , nsamples=nsamples)
                         #print("Setup done")
                         self.isConfig = True
                     if code is not None:
                         code = numpy.array(code)
                         code_block = code
                         if repeat is not None:
                             code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
                             #print(code_block.shape)
                     for i in range(self.buffer.shape[1]):
                         if code is not None:
                             self.buffer[:,i] = dataOut.data[:,i*self.step:i*self.step + self.nsamples]*code_block
                         else:
                             self.buffer[:,i]    = dataOut.data[:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
                         #self.buffer[:,j,self.__nHeis-j*self.step - self.nheights:self.__nHeis-j*self.step] = numpy.flip(dataOut.data[:,j*self.step:j*self.step + self.nheights])
                     for j in range(self.buffer.shape[0]):
                         self.sshProfiles[j] = numpy.transpose(self.buffer[j])
                     profileIndex  =  self.nsamples
                     deltaHeight   =  dataOut.heightList[1] - dataOut.heightList[0]
                     ippSeconds    =  (deltaHeight*1.0e-6)/(0.15)
                     #print("ippSeconds, dH: ",ippSeconds,deltaHeight)
                     try:
                         if dataOut.concat_m  is not None:
                             ippSeconds= ippSeconds/float(dataOut.concat_m)
                             #print "Profile concat %d"%dataOut.concat_m
                     except:
                         pass
                     dataOut.data            = self.sshProfiles
                     dataOut.flagNoData      = False
                     dataOut.heightList      = numpy.arange(self.buffer.shape[1]) *self.step*deltaHeight + dataOut.heightList[0]
                     dataOut.nProfiles       = int(dataOut.nProfiles*self.nsamples)
                     dataOut.profileIndex    = profileIndex
                     dataOut.flagDataAsBlock = True
                     dataOut.ippSeconds      = ippSeconds
                     dataOut.step            = self.step
                     #print(numpy.shape(dataOut.data))
                     #exit(1)
                     #print("new data shape and time:", dataOut.data.shape, dataOut.utctime)
                     return dataOut
             ################################################################################3############################3
             ################################################################################3############################3
             ################################################################################3############################3
             ################################################################################3############################3
             class SSheightProfiles2(Operation):
                 '''
                 Procesa por perfiles y por bloques
                 '''
                 bufferShape   = None
                 profileShape  = None
                 sshProfiles   = None
                 profileIndex  = None
                 #nsamples      = None
                 #step          = None
                 #deltaHeight = None
                 #init_range  = None
                 __slots__ = ('step', 'nsamples', 'deltaHeight', 'init_range', 'isConfig', '__nChannels',
                             '__nProfiles', '__nHeis', 'deltaHeight', 'new_nHeights')
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.isConfig = False
                 def setup(self,dataOut ,step = None , nsamples = None):
                     if step == None and nsamples == None:
                         raise ValueError("step or nheights should be specified ...")
                     self.step         = step
                     self.nsamples     = nsamples
                     self.__nChannels  = int(dataOut.nChannels)
                     self.__nProfiles  = int(dataOut.nProfiles)
                     self.__nHeis      = int(dataOut.nHeights)
                     residue     =  (self.__nHeis - self.nsamples) % self.step
                     if residue != 0:
                         print("The residue is %d, step=%d should be multiple of %d to avoid loss of %d samples"%(residue,step,shape[1] - self.nsamples,residue))
                     self.deltaHeight      =  dataOut.heightList[1] - dataOut.heightList[0]
                     self.init_range = dataOut.heightList[0]
                     #numberProfile    =  self.nsamples
                     numberSamples    =  (self.__nHeis - self.nsamples)/self.step
                     self.new_nHeights    =  numberSamples
                     self.bufferShape  = int(self.__nChannels), int(numberSamples), int(self.nsamples)  # nchannels, nsamples , nprofiles
                     self.profileShape = int(self.__nChannels), int(self.nsamples), int(numberSamples)  # nchannels, nprofiles, nsamples
                     self.buffer       = numpy.zeros(self.bufferShape , dtype=numpy.complex)
                     self.sshProfiles  = numpy.zeros(self.profileShape, dtype=numpy.complex)
                 def getNewProfiles(self, data, code=None, repeat=None):
                     if code is not None:
                         code = numpy.array(code)
                         code_block = code
                         if repeat is not None:
                             code_block = numpy.repeat(code_block, repeats=repeat, axis=1)
                     if data.ndim == 2:
                         data = data.reshape(1,1,self.__nHeis )
                     #print("buff, data, :",self.buffer.shape, data.shape,self.sshProfiles.shape)
                     for i in range(int(self.new_nHeights)): #nuevas alturas
                         if code is not None:
                             self.buffer[:,i,:] = data[:,:,i*self.step:i*self.step + self.nsamples]*code_block
                         else:
                             self.buffer[:,i,:] = data[:,:,i*self.step:i*self.step + self.nsamples]#*code[dataOut.profileIndex,:]
                     for j in range(self.__nChannels): #en los cananles
                         self.sshProfiles[j,:,:] = numpy.transpose(self.buffer[j,:,:])
                     #print("new profs Done")
                 def run(self, dataOut, step, nsamples, code = None, repeat = None):
                     if dataOut.flagNoData == True:
                         return dataOut
                     dataOut.flagNoData = True
                     #print("init  data shape:", dataOut.data.shape)
                     #print("ch: {}  prof: {}  hs: {}".format(int(dataOut.nChannels),
                     #                int(dataOut.nProfiles),int(dataOut.nHeights)))
                     profileIndex            = None
                     # if not dataOut.flagDataAsBlock:
                     #     dataOut.nProfiles = 1
                     if not self.isConfig:
                         self.setup(dataOut, step=step , nsamples=nsamples)
                         #print("Setup done")
                         self.isConfig = True
                     dataBlock = None
                     nprof = 1
                     if dataOut.flagDataAsBlock:
                         nprof = int(dataOut.nProfiles)
                     #print("dataOut nProfiles:", dataOut.nProfiles)
                     for profile in range(nprof):
                         if dataOut.flagDataAsBlock:
                             #print("read blocks")
                             self.getNewProfiles(dataOut.data[:,profile,:], code=code, repeat=repeat)
                         else:
                             #print("read profiles")
                             self.getNewProfiles(dataOut.data, code=code, repeat=repeat) #only one channe
                         if profile == 0:
                             dataBlock = self.sshProfiles.copy()
                         else:   #by blocks
                             dataBlock = numpy.concatenate((dataBlock,self.sshProfiles), axis=1) #profile axis
                             #print("by blocks: ",dataBlock.shape, self.sshProfiles.shape)
                     profileIndex  =  self.nsamples
                     #deltaHeight   =  dataOut.heightList[1] - dataOut.heightList[0]
                     ippSeconds    =  (self.deltaHeight*1.0e-6)/(0.15)
                     dataOut.data            = dataBlock
                     #print("show me: ",self.step,self.deltaHeight, dataOut.heightList, self.new_nHeights)
                     dataOut.heightList      = numpy.arange(int(self.new_nHeights)) *self.step*self.deltaHeight + self.init_range
                     dataOut.ippSeconds      = ippSeconds
                     dataOut.step            = self.step
                     dataOut.flagNoData      = False
                     if dataOut.flagDataAsBlock:
                         dataOut.nProfiles       = int(dataOut.nProfiles*self.nsamples)
                     else:
                         dataOut.nProfiles       = int(self.nsamples)
                     dataOut.profileIndex    = dataOut.nProfiles
                     dataOut.flagDataAsBlock = True
                     dataBlock = None
                     #print("new data shape:", dataOut.data.shape, dataOut.utctime)
                     return dataOut
-            #import skimage.color
-            #import skimage.io
-            #import matplotlib.pyplot as plt
             class removeProfileByFaradayHS(Operation):
                 '''
                 '''
                 __buffer_data = []
                 __buffer_times = []
                 buffer = None
                 outliers_IDs_list = []
                 __slots__ = ('n','navg','profileMargin','thHistOutlier','minHei_idx','maxHei_idx','nHeights',
                     '__dh','first_utcBlock','__profIndex','init_prof','end_prof','lenProfileOut','nChannels',
                     '__count_exec','__initime','__dataReady','__ipp')
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.isConfig = False
                 def setup(self,dataOut, n=None , navg=0.8, profileMargin=50,thHistOutlier=3, minHei=None, maxHei=None):
                     if n == None and timeInterval == None:
                         raise ValueError("nprofiles or timeInterval should be specified ...")
                     if n != None:
                         self.n = n
                     self.navg = navg
                     self.profileMargin = profileMargin
                     self.thHistOutlier = thHistOutlier
                     self.__profIndex = 0
                     self.buffer = None
                     self._ipp = dataOut.ippSeconds
                     self.n_prof_released = 0
                     self.heightList = dataOut.heightList
                     self.init_prof = 0
                     self.end_prof = 0
                     self.__count_exec = 0
                     self.__profIndex = 0
                     self.first_utcBlock = None
                     self.__dh = dataOut.heightList[1] - dataOut.heightList[0]
                     minHei = minHei
                     maxHei = maxHei
                     if minHei==None :
                         minHei = dataOut.heightList[0]
                     if maxHei==None :
                         maxHei = dataOut.heightList[-1]
                     self.minHei_idx,self.maxHei_idx =  getHei_index(minHei, maxHei, dataOut.heightList)
                     self.nChannels = dataOut.nChannels
                     self.nHeights = dataOut.nHeights
                 def filterSatsProfiles(self):
                     data = self.__buffer_data
                     #print(data.shape)
                     nChannels, profiles, heights = data.shape
                     indexes=[]
                     outliers_IDs=[]
                     for c in range(nChannels):
                         for h in range(self.minHei_idx, self.maxHei_idx):
                             power = data[c,:,h] * numpy.conjugate(data[c,:,h])
                             power = power.real
                             #power = (numpy.abs(data[c,:,h].real))
                             sortdata = numpy.sort(power, axis=None)
                             sortID=power.argsort()
                             index = _noise.hildebrand_sekhon2(sortdata,self.navg) #0.75-> buen valor
                             indexes.append(index)
                             outliers_IDs=numpy.append(outliers_IDs,sortID[index:])
                             # print(outliers_IDs)
                             # fig,ax = plt.subplots()
                             # #ax.set_title(str(k)+" "+str(j))
                             # x=range(len(sortdata))
                             # ax.scatter(x,sortdata)
                             # ax.axvline(index)
                             # plt.grid()
                             # plt.show()
                     outliers_IDs = outliers_IDs.astype(numpy.dtype('int64'))
                     outliers_IDs = numpy.unique(outliers_IDs)
                     outs_lines = numpy.sort(outliers_IDs)
                     # #print("outliers Ids: ", outs_lines, outs_lines.shape)
                     #hist, bin_edges = numpy.histogram(outs_lines, bins=10, density=True)
                     #Agrupando el histograma de outliers,
                     #my_bins = numpy.linspace(0,int(profiles), int(profiles/100), endpoint=False)
                     my_bins = numpy.linspace(0,9600, 96, endpoint=False)
                     hist, bins =  numpy.histogram(outs_lines,bins=my_bins)
                     hist_outliers_indexes = numpy.where(hist > self.thHistOutlier)       #es outlier
                     #print(hist_outliers_indexes[0])
                     bins_outliers_indexes = [int(i) for i in bins[hist_outliers_indexes]] #
                     #print(bins_outliers_indexes)
                     outlier_loc_index = []
                     # for n in range(len(bins_outliers_indexes)-1):
                     #     for k in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n+1]+self.profileMargin):
                     #         outlier_loc_index.append(k)
                     outlier_loc_index = [e for n in range(len(bins_outliers_indexes)-1) for e in range(bins_outliers_indexes[n]-self.profileMargin,bins_outliers_indexes[n+1]+self.profileMargin) ]
                     outlier_loc_index = numpy.asarray(outlier_loc_index)
                     #print(len(numpy.unique(outlier_loc_index)), numpy.unique(outlier_loc_index))
                     # x, y = numpy.meshgrid(numpy.arange(profiles), self.heightList)
                     # fig, ax = plt.subplots(1,2,figsize=(8, 6))
                     #
                     # dat = data[0,:,:].real
                     # m = numpy.nanmean(dat)
                     # o = numpy.nanstd(dat)
                     # #print(m, o, x.shape, y.shape)
                     # c = ax[0].pcolormesh(x, y, dat.T, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
                     # ax[0].vlines(outs_lines,200,600, linestyles='dashed', label = 'outs', color='w')
                     # fig.colorbar(c)
                     # ax[0].vlines(outlier_loc_index,650,750, linestyles='dashed', label = 'outs', color='r')
                     # ax[1].hist(outs_lines,bins=my_bins)
                     # plt.show()
                     self.outliers_IDs_list = numpy.unique(outlier_loc_index)
                     return data
                 def cleanOutliersByBlock(self):
+                    import matplotlib.pyplot as plt
+                    import datetime
                     #print(self.__buffer_data[0].shape)
                     data = self.__buffer_data#.copy()
-                    #print("cleaning shape inpt: ",data.shape)
+                    print("cleaning shape inpt: ",data.shape)
-                    '''
                     self.__buffer_data = []
-                    spectrum = numpy.fft.fft2(data, axes=(0,2))
+                    spectrum = numpy.fft.fft2(data[:,:,self.minHei_idx:], axes=(0,2))
-                    #print("spc : ",spectrum.shape)
+                    print("spc : ",spectrum.shape)
                     (nch,nsamples, nh) = spectrum.shape
                     data2 = None
                     #print(data.shape)
                     spectrum2 = spectrum.copy()
                     for ch in range(nch):
                         dh = self.__dh
                         dt1 = (dh*1.0e-6)/(0.15)
                         dt2 = self.__buffer_times[1]-self.__buffer_times[0]
                         freqv = numpy.fft.fftfreq(nh, d=dt1)
                         freqh = numpy.fft.fftfreq(self.n, d=dt2)
                         #print("spc loop: ")
                         x, y = numpy.meshgrid(numpy.sort(freqh),numpy.sort(freqv))
                         z = numpy.abs(spectrum[ch,:,:])
+                        phase = numpy.angle(spectrum[ch,:,:])
                         # Find all peaks higher than the 98th percentile
-                        peaks = z < numpy.percentile(z, 98)
+                        peaks = z < numpy.percentile(z, 99)
                         #print(peaks)
                         # Set those peak coefficients to zero
                         spectrum2 = spectrum2 * peaks.astype(int)
-                        data2 = numpy.fft.ifft2(spectrum2)
+                        data2 = numpy.fft.ifft2(spectrum2,axes=(0,2))
                         dat = numpy.log10(z.T)
-                        dat2 = numpy.log10(spectrum2.T)
+                        pdat = numpy.log10(phase.T)
+                        dat2 = numpy.log10(numpy.abs(spectrum2.T))
-                        # m = numpy.mean(dat)
-                        # o = numpy.std(dat)
+                        m = numpy.mean(dat)
-                        # fig, ax = plt.subplots(2,1,figsize=(8, 6))
+                        o = numpy.std(dat)
+                        fig, ax = plt.subplots(1,2,figsize=(12, 6))
-                        # c = ax[0].pcolormesh(x, y, dat, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
-                        # #c = ax.pcolor( z.T , cmap ='gray', vmin = (m-2*o), vmax = (m+2*o))
+                        colormap = 'jet'
-                        # date_time = datetime.datetime.fromtimestamp(self.__buffer_times[0]).strftime('%Y-%m-%d %H:%M:%S.%f')
+                        #c = ax[0].pcolormesh(x, y, dat, cmap =colormap, vmin = (m-2*o)/2, vmax = (m+2*o))
-                        # #strftime('%Y-%m-%d %H:%M:%S')
+                        c = ax[0].pcolormesh(x, y, dat, cmap =colormap, vmin = 4.2, vmax = 5.0)
-                        # ax[0].set_title('Spectrum magnitude '+date_time)
+                        fig.colorbar(c, ax=ax[0])
-                        # fig.canvas.set_window_title('Spectrum magnitude {} '.format(self.n)+date_time)
+                        #print("aqui estoy", dat.shape)
+                        #c = ax.pcolor( z.T , cmap ='gray', vmin = (m-2*o), vmax = (m+2*o))
+                        date_time = datetime.datetime.fromtimestamp(self.__buffer_times[0]).strftime('%Y-%m-%d %H:%M:%S.%f')
-                        # c = ax[1].pcolormesh(x, y, dat, cmap ='YlGnBu', vmin = (m-2*o), vmax = (m+2*o))
+                        #strftime('%Y-%m-%d %H:%M:%S')
-                        # fig.colorbar(c)
+                        #ax[0].set_title('Spectrum magnitude '+date_time)
-                        # plt.show()
+                        fig.canvas.set_window_title('Spectrum magnitude {} '.format(self.n)+date_time)
+                        #print("aqui estoy2",dat2[:,:,0].shape)
+                        c = ax[1].pcolormesh(x, y, pdat, cmap =colormap, vmin = -0.0, vmax = 0.5)
+                        #c = ax[1].pcolormesh(x, y, dat2[:,:,0], cmap =colormap, vmin = (m-2*o)/2, vmax = (m+2*o)-1)
+                        #print("aqui estoy3")
+                        fig.colorbar(c, ax=ax[1])
+                        plt.show()
                     #print(data2.shape)
-                    data = data2
+                    #data = data2
+                    #cleanBlock = numpy.fft.ifft2(data, axes=(0,2)).reshape()
-                    #cleanBlock = numpy.fft.ifft2(spectrum, axes=(0,2)).reshape()
-                    '''
                     #print("cleanOutliersByBlock Done")
-                    return self.filterSatsProfiles()
+                    return data
                 def fillBuffer(self, data, datatime):
                     if self.__profIndex == 0:
                         self.__buffer_data = data.copy()
                     else:
                         self.__buffer_data = numpy.concatenate((self.__buffer_data,data), axis=1)#en perfiles
                     self.__profIndex += 1
-                    #self.__buffer_times.append(datatime)
+                    self.__buffer_times.append(datatime)
                 def getData(self, data, datatime=None):
                     if self.__profIndex  == 0:
                         self.__initime = datatime
                     self.__dataReady = False
                     self.fillBuffer(data, datatime)
                     dataBlock = None
                     if self.__profIndex == self.n:
                         #print("apnd : ",data)
                         #dataBlock = self.cleanOutliersByBlock()
                         dataBlock = self.filterSatsProfiles()
                         self.__dataReady = True
                     return dataBlock
                     if dataBlock is None:
                         return None, None
                     return dataBlock
                 def releaseBlock(self):
                     if self.n % self.lenProfileOut != 0:
                         raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
                         return None
                     data = self.buffer[:,self.init_prof:self.end_prof:,:]  #ch, prof, alt
                     self.init_prof = self.end_prof
                     self.end_prof += self.lenProfileOut
                     #print("data release shape: ",dataOut.data.shape, self.end_prof)
                     self.n_prof_released += 1
                     #print("f_no_data ", dataOut.flagNoData)
                     return data
                 def run(self, dataOut, n=None, navg=0.8, nProfilesOut=1, profile_margin=50,th_hist_outlier=3,minHei=None, maxHei=None):
                     #print("run op buffer 2D",dataOut.ippSeconds)
                     # self.nChannels = dataOut.nChannels
                     # self.nHeights = dataOut.nHeights
                     if not self.isConfig:
                         #print("init p idx: ", dataOut.profileIndex )
                         self.setup(dataOut,n=n, navg=navg,profileMargin=profile_margin,
                                 thHistOutlier=th_hist_outlier,minHei=minHei, maxHei=maxHei)
                         self.isConfig = True
                     dataBlock = None
                     if not dataOut.buffer_empty: #hay datos acumulados
                         if self.init_prof == 0:
                             self.n_prof_released = 0
                             self.lenProfileOut = nProfilesOut
                             dataOut.flagNoData = False
                             #print("tp 2 ",dataOut.data.shape)
                             self.init_prof = 0
                             self.end_prof = self.lenProfileOut
                             dataOut.nProfiles = self.lenProfileOut
                         if nProfilesOut == 1:
                             dataOut.flagDataAsBlock = False
                         else:
                             dataOut.flagDataAsBlock = True
                         #print("prof: ",self.init_prof)
                         dataOut.flagNoData = False
                         if numpy.isin(self.n_prof_released, self.outliers_IDs_list):
                             #print("omitting: ", self.n_prof_released)
                             dataOut.flagNoData = True
                         dataOut.ippSeconds = self._ipp
                         dataOut.utctime = self.first_utcBlock + self.init_prof*self._ipp
                         # print("time: ", dataOut.utctime, self.first_utcBlock, self.init_prof,self._ipp,dataOut.ippSeconds)
                         #dataOut.data = self.releaseBlock()
                         #########################################################3
                         if self.n % self.lenProfileOut != 0:
                             raise ValueError("lenProfileOut %d must be submultiple of nProfiles %d" %(self.lenProfileOut, self.n))
                             return None
                         dataOut.data = self.buffer[:,self.init_prof:self.end_prof:,:]  #ch, prof, alt
                         self.init_prof = self.end_prof
                         self.end_prof += self.lenProfileOut
                         #print("data release shape: ",dataOut.data.shape, self.end_prof)
                         self.n_prof_released += 1
                         if self.end_prof >= (self.n +self.lenProfileOut):
                             self.init_prof = 0
                             self.__profIndex = 0
                             self.buffer = None
                             dataOut.buffer_empty = True
                             self.outliers_IDs_list = []
                             self.n_prof_released = 0
                             dataOut.flagNoData = False #enviar ultimo aunque sea outlier :(
                             #print("cleaning...", dataOut.buffer_empty)
                         dataOut.profileIndex    = 0 #self.lenProfileOut
                         ####################################################################
                         return dataOut
                     #print("tp 223 ",dataOut.data.shape)
                     dataOut.flagNoData = True
                     try:
                         #dataBlock = self.getData(dataOut.data.reshape(self.nChannels,1,self.nHeights), dataOut.utctime)
                         dataBlock = self.getData(numpy.reshape(dataOut.data,(self.nChannels,1,self.nHeights)), dataOut.utctime)
                         self.__count_exec +=1
                     except Exception as e:
                         print("Error getting profiles data",self.__count_exec )
                         print(e)
                         sys.exit()
                     if self.__dataReady:
                         #print("omitting: ", len(self.outliers_IDs_list))
                         self.__count_exec = 0
                         #dataOut.data =
                         #self.buffer = numpy.flip(dataBlock, axis=1)
                         self.buffer = dataBlock
                         self.first_utcBlock = self.__initime
                         dataOut.utctime = self.__initime
                         dataOut.nProfiles  = self.__profIndex
                         #dataOut.flagNoData = False
                         self.init_prof = 0
                         self.__profIndex = 0
                         self.__initime = None
                         dataBlock = None
                         self.__buffer_times = []
                         dataOut.error = False
                         dataOut.useInputBuffer = True
                         dataOut.buffer_empty = False
                         #print("1 ch: {}  prof: {}  hs: {}".format(int(dataOut.nChannels),int(dataOut.nProfiles),int(dataOut.nHeights)))
                     #print(self.__count_exec)
                     return dataOut
             class RemoveProfileSats(Operation):
                 '''
                 Omite los perfiles contaminados con señal de satelites,
                 In: minHei = min_sat_range
                     max_sat_range
                     min_hei_ref
                     max_hei_ref
                     th = diference between profiles mean, ref and sats
                 Out:
                     profile clean
                 '''
                 isConfig = False
                 min_sats = 0
                 max_sats = 999999999
                 min_ref= 0
                 max_ref= 9999999999
                 needReshape = False
                 count = 0
                 thdB = 0
                 byRanges = False
                 min_sats = None
                 max_sats = None
                 noise = 0
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.isConfig = False
                 def setup(self, dataOut, minHei, maxHei, minRef, maxRef, th, thdB, rangeHeiList):
                     if rangeHeiList!=None:
                         self.byRanges = True
                     else:
                         if minHei==None or  maxHei==None :
                             raise ValueError("Parameters heights are required")
                     if minRef==None or  maxRef==None:
                         raise ValueError("Parameters heights are required")
                     if self.byRanges:
                         self.min_sats = []
                         self.max_sats = []
                         for min,max in rangeHeiList:
                             a,b = getHei_index(min, max, dataOut.heightList)
                             self.min_sats.append(a)
                             self.max_sats.append(b)
                     else:
                         self.min_sats, self.max_sats = getHei_index(minHei, maxHei, dataOut.heightList)
                     self.min_ref, self.max_ref = getHei_index(minRef, maxRef, dataOut.heightList)
                     self.th = th
                     self.thdB = thdB
                     self.isConfig = True
                 def compareRanges(self,data, minHei,maxHei):
                     # ref =  data[0,self.min_ref:self.max_ref] * numpy.conjugate(data[0,self.min_ref:self.max_ref])
                     # p_ref = 10*numpy.log10(ref.real)
                     # m_ref =  numpy.mean(p_ref)
                     m_ref = self.noise
                     sats = data[0,minHei:maxHei] * numpy.conjugate(data[0,minHei:maxHei])
                     p_sats = 10*numpy.log10(sats.real)
                     m_sats = numpy.mean(p_sats)
                     if m_sats > (m_ref + self.th): #and (m_sats > self.thdB):
                         #print("msats: ",m_sats," \tmRef: ", m_ref, "\t",(m_sats - m_ref))
                         #print("Removing profiles...")
                         return False
                     return True
                 def isProfileClean(self, data):
                     '''
                     Analiza solo 1 canal, y descarta todos...
                     '''
                     clean = True
                     if self.byRanges:
                         for n in range(len(self.min_sats)):
                             c = self.compareRanges(data,self.min_sats[n],self.max_sats[n])
                             clean = clean and c
                     else:
                         clean =  (self.compareRanges(data, self.min_sats,self.max_sats))
                     return clean
                 def run(self, dataOut, minHei=None, maxHei=None, minRef=None, maxRef=None, th=5, thdB=65, rangeHeiList=None):
                     dataOut.flagNoData = True
                     if not self.isConfig:
                         self.setup(dataOut, minHei, maxHei, minRef, maxRef, th, thdB, rangeHeiList)
                         self.isConfig = True
                     #print(self.min_sats,self.max_sats)
                     if dataOut.flagDataAsBlock:
                         raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
                     else:
                         self.noise =10*numpy.log10(dataOut.getNoisebyHildebrand(ymin_index=self.min_ref, ymax_index=self.max_ref))
                         if not self.isProfileClean(dataOut.data):
                             return dataOut
                             #dataOut.data = numpy.full((dataOut.nChannels,dataOut.nHeights),numpy.NAN)
                     #self.count += 1
                     dataOut.flagNoData = False
                     return 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