schain Commit - r611:cdbd858cadba · Jicamarca Repository

Bug fixed: Padding decode data with zeros at the first heights was eliminated.

Miguel Valdez -

r611:cdbd858cadba

parent child

Context file:

r611:cdbd858cadba

Collapse all files

schainpy/model/data/jroamisr.py +1 -1

             import numpy
             import copy
             class Beam:
                 def __init__(self):
                     self.codeList = []
                     self.azimuthList = []
                     self.zenithList = []
             class AMISR:
                 def __init__(self):
                     self.flagNoData = True
                     self.data = None
                     self.utctime = None
                     self.type = "AMISR"
                     #propiedades para compatibilidad con Voltages
                     self.timeZone = 0#timezone like jroheader, difference in minutes between UTC and localtime
                     self.dstFlag = 0#self.dataIn.dstFlag
                     self.errorCount = 0#self.dataIn.errorCount
                     self.useLocalTime = True#self.dataIn.useLocalTime
                     self.radarControllerHeaderObj = None#self.dataIn.radarControllerHeaderObj.copy()
                     self.systemHeaderObj = None#self.dataIn.systemHeaderObj.copy()
                     self.channelList = [0]#self.dataIn.channelList esto solo aplica para el caso de AMISR
                     self.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
                     self.flagDiscontinuousBlock = None#self.dataIn.flagDiscontinuousBlock
                     #self.utctime = #self.firstdatatime
                     self.flagDecodeData = None#self.dataIn.flagDecodeData #asumo q la data esta decodificada
                     self.flagDeflipData = None#self.dataIn.flagDeflipData #asumo q la data esta sin flip
                     self.nCohInt = 1#self.dataIn.nCohInt
                     self.nIncohInt = 1
                     self.ippSeconds = None#self.dataIn.ippSeconds, segun el filename/Setup/Tufile
                     self.windowOfFilter = None#self.dataIn.windowOfFilter
                     self.timeInterval = None#self.dataIn.timeInterval*self.dataOut.nFFTPoints*self.dataOut.nIncohInt
                     self.frequency = None#self.dataIn.frequency
                     self.realtime = 0#self.dataIn.realtime
                     #actualizar en la lectura de datos
                     self.heightList = None#self.dataIn.heightList
                     self.nProfiles = None#Number of samples or nFFTPoints
                     self.nRecords = None
                     self.nBeams = None
                     self.nBaud = None#self.dataIn.nBaud
                     self.nCode = None#self.dataIn.nCode
                     self.code = None#self.dataIn.code
                     #consideracion para los Beams
                     self.beamCodeDict = None
                     self.beamRangeDict = None
                     self.beamcode = None
                     self.azimuth = None
                     self.zenith = None
                     self.gain = None
                     self.npulseByFrame = None
                     self.profileIndex = None
                     self.beam = Beam()
                 def copy(self, inputObj=None):
-                    if inputObj == None:
+                    if inputObj is None:
                         return copy.deepcopy(self)
                     for key in inputObj.__dict__.keys():
                         self.__dict__[key] = inputObj.__dict__[key]
                 def getNHeights(self):
                     return len(self.heightList)
                 def isEmpty(self):
                     return self.flagNoData
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt
                     return timeInterval
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
                 nHeights = property(getNHeights, "I'm the 'nHeights' property.")

schainpy/model/data/jrodata.py +2 -2

             '''
             $Author: murco $
             $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
             '''
             import copy
             import numpy
             import datetime
             from jroheaderIO import SystemHeader, RadarControllerHeader
             def getNumpyDtype(dataTypeCode):
                 if dataTypeCode == 0:
                     numpyDtype = numpy.dtype([('real','<i1'),('imag','<i1')])
                 elif dataTypeCode == 1:
                     numpyDtype = numpy.dtype([('real','<i2'),('imag','<i2')])
                 elif dataTypeCode == 2:
                     numpyDtype = numpy.dtype([('real','<i4'),('imag','<i4')])
                 elif dataTypeCode == 3:
                     numpyDtype = numpy.dtype([('real','<i8'),('imag','<i8')])
                 elif dataTypeCode == 4:
                     numpyDtype = numpy.dtype([('real','<f4'),('imag','<f4')])
                 elif dataTypeCode == 5:
                     numpyDtype = numpy.dtype([('real','<f8'),('imag','<f8')])
                 else:
                     raise ValueError, 'dataTypeCode was not defined'
                 return numpyDtype
             def getDataTypeCode(numpyDtype):
                 if numpyDtype == numpy.dtype([('real','<i1'),('imag','<i1')]):
                     datatype = 0
                 elif numpyDtype == numpy.dtype([('real','<i2'),('imag','<i2')]):
                     datatype = 1
                 elif numpyDtype == numpy.dtype([('real','<i4'),('imag','<i4')]):
                     datatype = 2
                 elif numpyDtype == numpy.dtype([('real','<i8'),('imag','<i8')]):
                     datatype = 3
                 elif numpyDtype == numpy.dtype([('real','<f4'),('imag','<f4')]):
                     datatype = 4
                 elif numpyDtype == numpy.dtype([('real','<f8'),('imag','<f8')]):
                     datatype = 5
                 else:
                     datatype = None
                 return datatype
             def hildebrand_sekhon(data, navg):
                 """
                 This method is for the objective determination of the noise level in Doppler spectra. This
                 implementation technique is based on the fact that the standard deviation of the spectral
                 densities is equal to the mean spectral density for white Gaussian noise
                 Inputs:
                     Data    :    heights
                     navg    :    numbers of averages
                 Return:
                     -1        :    any error
                     anoise    :    noise's level
                 """
                 sortdata = numpy.sort(data,axis=None)
                 lenOfData = len(sortdata)
                 nums_min = lenOfData/10
                 if (lenOfData/10) > 2:
                     nums_min = lenOfData/10
                 else:
                     nums_min = 2
                 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
                 stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1))
                 return lnoise
             class Beam:
                 def __init__(self):
                     self.codeList = []
                     self.azimuthList = []
                     self.zenithList = []
             class GenericData(object):
                 flagNoData = True
                 def __init__(self):
                     raise ValueError, "This class has not been implemented"
                 def copy(self, inputObj=None):
                     if inputObj == None:
                         return copy.deepcopy(self)
                     for key in inputObj.__dict__.keys():
                         self.__dict__[key] = inputObj.__dict__[key]
                 def deepcopy(self):
                     return copy.deepcopy(self)
                 def isEmpty(self):
                     return self.flagNoData
             class JROData(GenericData):
             #    m_BasicHeader = BasicHeader()
             #    m_ProcessingHeader = ProcessingHeader()
                 systemHeaderObj = SystemHeader()
                 radarControllerHeaderObj = RadarControllerHeader()
             #    data = None
                 type = None
                 datatype = None #dtype but in string
             #     dtype = None
             #    nChannels = None
             #    nHeights = None
                 nProfiles = None
                 heightList = None
                 channelList = None
                 flagDiscontinuousBlock = False
                 useLocalTime = False
                 utctime = None
                 timeZone = None
                 dstFlag = None
                 errorCount = None
                 blocksize = None
             #     nCode = None
             #
             #     nBaud = None
             #
             #     code = None
                 flagDecodeData = False #asumo q la data no esta decodificada
                 flagDeflipData = False #asumo q la data no esta sin flip
                 flagShiftFFT = False
             #     ippSeconds = None
             #     timeInterval = None
                 nCohInt = None
             #     noise = None
                 windowOfFilter = 1
                 #Speed of ligth
                 C = 3e8
                 frequency = 49.92e6
                 realtime = False
                 beacon_heiIndexList = None
                 last_block = None
                 blocknow = None
                 azimuth = None
                 zenith = None
                 beam = Beam()
                 profileIndex = None
                 def __init__(self):
                     raise ValueError, "This class has not been implemented"
                 def getNoise(self):
                     raise ValueError, "Not implemented"
                 def getNChannels(self):
                     return len(self.channelList)
                 def getChannelIndexList(self):
                     return range(self.nChannels)
                 def getNHeights(self):
                     return len(self.heightList)
                 def getHeiRange(self, extrapoints=0):
                     heis = self.heightList
             #        deltah = self.heightList[1] - self.heightList[0]
             #
             #        heis.append(self.heightList[-1])
                     return heis
                 def getltctime(self):
                     if self.useLocalTime:
                         return self.utctime - self.timeZone*60
                     return self.utctime
                 def getDatatime(self):
                     datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
                     return datatimeValue
                 def getTimeRange(self):
                     datatime = []
                     datatime.append(self.ltctime)
                     datatime.append(self.ltctime + self.timeInterval+60)
                     datatime = numpy.array(datatime)
                     return datatime
                 def getFmax(self):
                     PRF = 1./(self.ippSeconds * self.nCohInt)
                     fmax = PRF/2.
                     return fmax
                 def getVmax(self):
                     _lambda = self.C/self.frequency
                     vmax = self.getFmax() * _lambda
                     return vmax
                 def get_ippSeconds(self):
                     '''
                     '''
                     return self.radarControllerHeaderObj.ippSeconds
                 def set_ippSeconds(self, ippSeconds):
                     '''
                     '''
                     self.radarControllerHeaderObj.ippSeconds = ippSeconds
                     return
                 def get_dtype(self):
                     '''
                     '''
                     return getNumpyDtype(self.datatype)
                 def set_dtype(self, numpyDtype):
                     '''
                     '''
                     self.datatype = getDataTypeCode(numpyDtype)
                 def get_code(self):
                     '''
                     '''
                     return self.radarControllerHeaderObj.code
                 def set_code(self, code):
                     '''
                     '''
                     self.radarControllerHeaderObj.code = code
                     return
                 def get_ncode(self):
                     '''
                     '''
                     return self.radarControllerHeaderObj.nCode
                 def set_ncode(self, nCode):
                     '''
                     '''
                     self.radarControllerHeaderObj.nCode = nCode
                     return
                 def get_nbaud(self):
                     '''
                     '''
                     return self.radarControllerHeaderObj.nBaud
                 def set_nbaud(self, nBaud):
                     '''
                     '''
                     self.radarControllerHeaderObj.nBaud = nBaud
                     return
             #     def getTimeInterval(self):
             #
             #         raise IOError, "This method should be implemented inside each Class"
                 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
                 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
                 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
                 #noise = property(getNoise, "I'm the 'nHeights' property.")
                 datatime = property(getDatatime, "I'm the 'datatime' property")
                 ltctime = property(getltctime, "I'm the 'ltctime' property")
                 ippSeconds = property(get_ippSeconds, set_ippSeconds)
                 dtype = property(get_dtype, set_dtype)
             #     timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
                 code = property(get_code, set_code)
                 nCode = property(get_ncode, set_ncode)
                 nBaud = property(get_nbaud, set_nbaud)
             class Voltage(JROData):
                 #data es un numpy array de 2 dmensiones (canales, alturas)
                 data = None
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.useLocalTime = True
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Voltage"
                     self.data = None
             #         self.dtype = None
             #        self.nChannels = 0
             #        self.nHeights = 0
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #        self.channelIndexList = None
                     self.flagNoData = True
                     self.flagDiscontinuousBlock = False
                     self.utctime = None
                     self.timeZone = None
                     self.dstFlag = None
                     self.errorCount = None
                     self.nCohInt = None
                     self.blocksize = None
                     self.flagDecodeData = False #asumo q la data no esta decodificada
                     self.flagDeflipData = False #asumo q la data no esta sin flip
                     self.flagShiftFFT = False
                     self.flagDataAsBlock = False    #Asumo que la data es leida perfil a perfil
                     self.profileIndex = 0
                 def getNoisebyHildebrand(self, channel = None):
                     """
                     Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
                     Return:
                         noiselevel
                     """
                     if channel != None:
                         data = self.data[channel]
                         nChannels = 1
                     else:
                         data = self.data
                         nChannels = self.nChannels
                     noise = numpy.zeros(nChannels)
                     power = data * numpy.conjugate(data)
                     for thisChannel in range(nChannels):
                         if nChannels == 1:
                             daux = power[:].real
                         else:
                             daux = power[thisChannel,:].real
                         noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
                     return noise
                 def getNoise(self, type = 1, channel = None):
                     if type == 1:
                         noise = self.getNoisebyHildebrand(channel)
                     return 10*numpy.log10(noise)
                 def getPower(self, channel = None):
                     if channel != None:
                         data = self.data[channel]
                     else:
                         data = self.data
                     power = data * numpy.conjugate(data)
                     return 10*numpy.log10(power.real)
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt
                     return timeInterval
                 noise = property(getNoise, "I'm the 'nHeights' property.")
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
             class Spectra(JROData):
                 #data es un numpy array de 2 dmensiones (canales, perfiles, alturas)
                 data_spc = None
                 #data es un numpy array de 2 dmensiones (canales, pares, alturas)
                 data_cspc = None
                 #data es un numpy array de 2 dmensiones (canales, alturas)
                 data_dc = None
                 nFFTPoints = None
             #     nPairs = None
                 pairsList = None
                 nIncohInt = None
                 wavelength = None #Necesario para cacular el rango de velocidad desde la frecuencia
                 nCohInt = None #se requiere para determinar el valor de timeInterval
                 ippFactor = None
                 profileIndex = 0
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.useLocalTime = True
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Spectra"
             #        self.data = None
             #         self.dtype = None
             #        self.nChannels = 0
             #        self.nHeights = 0
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #        self.channelIndexList = None
                     self.pairsList = None
                     self.flagNoData = True
                     self.flagDiscontinuousBlock = False
                     self.utctime = None
                     self.nCohInt = None
                     self.nIncohInt = None
                     self.blocksize = None
                     self.nFFTPoints = None
                     self.wavelength = None
                     self.flagDecodeData = False #asumo q la data no esta decodificada
                     self.flagDeflipData = False #asumo q la data no esta sin flip
                     self.flagShiftFFT = False
                     self.ippFactor = 1
                     #self.noise = None
                     self.beacon_heiIndexList = []
                     self.noise_estimation = None
                 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
                     """
                     Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
                     Return:
                         noiselevel
                     """
                     noise = numpy.zeros(self.nChannels)
                     for channel in range(self.nChannels):
                         daux = self.data_spc[channel,xmin_index:xmax_index,ymin_index:ymax_index]
                         noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
                     return noise
                 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
                     if self.noise_estimation != None:
                         return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py
                     else:
                         noise = self.getNoisebyHildebrand(xmin_index, xmax_index, ymin_index, ymax_index)
                         return noise
                 def getFreqRange(self, extrapoints=0):
                     deltafreq = self.getFmax() / (self.nFFTPoints*self.ippFactor)
                     freqrange = deltafreq*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltafreq/2
                     return freqrange
                 def getVelRange(self, extrapoints=0):
                     deltav = self.getVmax() / (self.nFFTPoints*self.ippFactor)
                     velrange = deltav*(numpy.arange(self.nFFTPoints+extrapoints)-self.nFFTPoints/2.) - deltav/2
                     return velrange
                 def getNPairs(self):
                     return len(self.pairsList)
                 def getPairsIndexList(self):
                     return range(self.nPairs)
                 def getNormFactor(self):
                     pwcode = 1
                     if self.flagDecodeData:
                         pwcode = numpy.sum(self.code[0]**2)
                     #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
                     normFactor = self.nProfiles*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
                     return normFactor
                 def getFlagCspc(self):
-                    if self.data_cspc == None:
+                    if self.data_cspc is None:
                         return True
                     return False
                 def getFlagDc(self):
-                    if self.data_dc == None:
+                    if self.data_dc is None:
                         return True
                     return False
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles
                     return timeInterval
                 nPairs = property(getNPairs, "I'm the 'nPairs' property.")
                 pairsIndexList = property(getPairsIndexList, "I'm the 'pairsIndexList' property.")
                 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
                 flag_cspc = property(getFlagCspc)
                 flag_dc = property(getFlagDc)
                 noise = property(getNoise, "I'm the 'nHeights' property.")
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
             class SpectraHeis(Spectra):
                 data_spc = None
                 data_cspc = None
                 data_dc = None
                 nFFTPoints = None
             #     nPairs = None
                 pairsList = None
                 nCohInt = None
                 nIncohInt = None
                 def __init__(self):
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "SpectraHeis"
             #         self.dtype = None
             #        self.nChannels = 0
             #        self.nHeights = 0
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #        self.channelIndexList = None
                     self.flagNoData = True
                     self.flagDiscontinuousBlock = False
             #         self.nPairs = 0
                     self.utctime = None
                     self.blocksize = None
                     self.profileIndex = 0
                     self.nCohInt = 1
                     self.nIncohInt = 1
                 def getNormFactor(self):
                     pwcode = 1
                     if self.flagDecodeData:
                         pwcode = numpy.sum(self.code[0]**2)
                     normFactor = self.nIncohInt*self.nCohInt*pwcode
                     return normFactor
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
                     return timeInterval
                 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
             class Fits(JROData):
                 heightList = None
                 channelList = None
                 flagNoData = True
                 flagDiscontinuousBlock = False
                 useLocalTime = False
                 utctime = None
                 timeZone = None
             #     ippSeconds = None
             #     timeInterval = None
                 nCohInt = None
                 nIncohInt = None
                 noise = None
                 windowOfFilter = 1
                 #Speed of ligth
                 C = 3e8
                 frequency = 49.92e6
                 realtime = False
                 def __init__(self):
                     self.type = "Fits"
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #         self.channelIndexList = None
                     self.flagNoData = True
                     self.utctime = None
                     self.nCohInt = 1
                     self.nIncohInt = 1
                     self.useLocalTime = True
                     self.profileIndex = 0
             #         self.utctime = None
             #         self.timeZone = None
             #         self.ltctime = None
             #         self.timeInterval = None
             #         self.header = None
             #         self.data_header = None
             #         self.data = None
             #         self.datatime = None
             #         self.flagNoData = False
             #         self.expName = ''
             #         self.nChannels = None
             #         self.nSamples = None
             #         self.dataBlocksPerFile = None
             #         self.comments = ''
             #
                 def getltctime(self):
                     if self.useLocalTime:
                         return self.utctime - self.timeZone*60
                     return self.utctime
                 def getDatatime(self):
                     datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
                     return datatime
                 def getTimeRange(self):
                     datatime = []
                     datatime.append(self.ltctime)
                     datatime.append(self.ltctime + self.timeInterval)
                     datatime = numpy.array(datatime)
                     return datatime
                 def getHeiRange(self):
                     heis = self.heightList
                     return heis
                 def isEmpty(self):
                     return self.flagNoData
                 def getNHeights(self):
                     return len(self.heightList)
                 def getNChannels(self):
                     return len(self.channelList)
                 def getChannelIndexList(self):
                     return range(self.nChannels)
                 def getNoise(self, type = 1):
                     #noise = numpy.zeros(self.nChannels)
                     if type == 1:
                         noise = self.getNoisebyHildebrand()
                     if type == 2:
                         noise = self.getNoisebySort()
                     if type == 3:
                         noise = self.getNoisebyWindow()
                     return noise
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
                     return timeInterval
                 datatime = property(getDatatime, "I'm the 'datatime' property")
                 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
                 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
                 channelIndexList = property(getChannelIndexList, "I'm the 'channelIndexList' property.")
                 noise = property(getNoise, "I'm the 'nHeights' property.")
                 datatime = property(getDatatime, "I'm the 'datatime' property")
                 ltctime = property(getltctime, "I'm the 'ltctime' property")
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
             class Correlation(JROData):
                 noise = None
                 SNR = None
                 pairsAutoCorr = None    #Pairs of Autocorrelation
                 #--------------------------------------------------
                 data_corr = None
                 data_volt = None
                 lagT = None # each element value is a profileIndex
                 lagR = None # each element value is in km
                 pairsList = None
                 calculateVelocity = None
                 nPoints = None
                 nAvg = None
                 bufferSize = None
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Correlation"
                     self.data = None
                     self.dtype = None
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
                     self.flagNoData = True
                     self.flagDiscontinuousBlock = False
                     self.utctime = None
                     self.timeZone = None
                     self.dstFlag = None
                     self.errorCount = None
                     self.blocksize = None
                     self.flagDecodeData = False #asumo q la data no esta decodificada
                     self.flagDeflipData = False #asumo q la data no esta sin flip
                     self.pairsList = None
                     self.nPoints = None
                 def getLagTRange(self, extrapoints=0):
                     lagTRange = self.lagT
                     diff = lagTRange[1] - lagTRange[0]
                     extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1]
                     lagTRange = numpy.hstack((lagTRange, extra))
                     return lagTRange
                 def getLagRRange(self, extrapoints=0):
                     return self.lagR
                 def getPairsList(self):
                     return self.pairsList
                 def getCalculateVelocity(self):
                     return self.calculateVelocity
                 def getNPoints(self):
                     return self.nPoints
                 def getNAvg(self):
                     return self.nAvg
                 def getBufferSize(self):
                     return self.bufferSize
                 def getPairsAutoCorr(self):
                     pairsList = self.pairsList
                     pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan
                     for l in range(len(pairsList)):
                         firstChannel = pairsList[l][0]
                         secondChannel = pairsList[l][1]
                         #Obteniendo pares de Autocorrelacion
                         if firstChannel == secondChannel:
                             pairsAutoCorr[firstChannel] = int(l)
                     pairsAutoCorr = pairsAutoCorr.astype(int)
                     return pairsAutoCorr
                 def getNoise(self, mode = 2):
                     indR = numpy.where(self.lagR == 0)[0][0]
                     indT = numpy.where(self.lagT == 0)[0][0]
                     jspectra0 = self.data_corr[:,:,indR,:]
                     jspectra = copy.copy(jspectra0)
                     num_chan = jspectra.shape[0]
                     num_hei = jspectra.shape[2]
                     freq_dc = jspectra.shape[1]/2
                     ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
                     if ind_vel[0]<0:
                         ind_vel[range(0,1)] = ind_vel[range(0,1)] + self.num_prof
                     if mode == 1:
                         jspectra[:,freq_dc,:] = (jspectra[:,ind_vel[1],:] + jspectra[:,ind_vel[2],:])/2 #CORRECCION
                     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(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
                     noise = jspectra0[:,freq_dc,:] - jspectra[:,freq_dc,:]
                     return noise
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt * self.nPoints
                     return timeInterval
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
             #     pairsList = property(getPairsList, "I'm the 'pairsList' property.")
             #     nPoints = property(getNPoints, "I'm the 'nPoints' property.")
                 calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.")
                 nAvg = property(getNAvg, "I'm the 'nAvg' property.")
                 bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.")
             class Parameters(JROData):
                 #Information from previous data
                 inputUnit = None        #Type of data to be processed
                 operation = None        #Type of operation to parametrize
                 normFactor = None       #Normalization Factor
                 groupList = None        #List of Pairs, Groups, etc
                 #Parameters
                 data_param = None       #Parameters obtained
                 data_pre = None         #Data Pre Parametrization
                 data_SNR = None         #Signal to Noise Ratio
             #     heightRange = None      #Heights
                 abscissaList = None    #Abscissa, can be velocities, lags or time
                 noise = None            #Noise Potency
                 utctimeInit = None      #Initial UTC time
                 paramInterval = None    #Time interval to calculate Parameters in seconds
                 #Fitting
                 data_error = None       #Error of the estimation
                 constants = None
                 library = None
                 #Output signal
                 outputInterval = None     #Time interval to calculate output signal in seconds
                 data_output = None       #Out signal
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Parameters"
                 def getTimeRange1(self):
                     datatime = []
                     if self.useLocalTime:
                         time1 = self.utctimeInit - self.timeZone*60
                     else:
                         time1 = utctimeInit
             #         datatime.append(self.utctimeInit)
             #         datatime.append(self.utctimeInit + self.outputInterval)
                     datatime.append(time1)
                     datatime.append(time1 + self.outputInterval)
                     datatime = numpy.array(datatime)
                     return datatime

schainpy/model/io/jroIO_amisr.py +1 -1

             '''
             @author: Daniel Suarez
             '''
             import os
             import sys
             import glob
             import fnmatch
             import datetime
             import time
             import re
             import h5py
             import numpy
             from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
             from schainpy.model.data.jroamisr import AMISR
             try:
                 from gevent import sleep
             except:
                 from time import sleep
             class RadacHeader():
                 def __init__(self, fp):
                     header = 'Raw11/Data/RadacHeader'
                     self.beamCodeByPulse = fp.get(header+'/BeamCode')
                     self.beamCode = fp.get('Raw11/Data/Beamcodes')
                     self.code = fp.get(header+'/Code')
                     self.frameCount = fp.get(header+'/FrameCount')
                     self.modeGroup = fp.get(header+'/ModeGroup')
                     self.nsamplesPulse = fp.get(header+'/NSamplesPulse')
                     self.pulseCount = fp.get(header+'/PulseCount')
                     self.radacTime = fp.get(header+'/RadacTime')
                     self.timeCount = fp.get(header+'/TimeCount')
                     self.timeStatus = fp.get(header+'/TimeStatus')
                     self.nrecords = self.pulseCount.shape[0] #nblocks
                     self.npulses = self.pulseCount.shape[1] #nprofile
                     self.nsamples = self.nsamplesPulse[0,0] #ngates
                     self.nbeams = self.beamCode.shape[1]
                 def getIndexRangeToPulse(self, idrecord=0):
                     #indexToZero = numpy.where(self.pulseCount.value[idrecord,:]==0)
                     #startPulseCountId = indexToZero[0][0]
                     #endPulseCountId = startPulseCountId - 1
                     #range1 = numpy.arange(startPulseCountId,self.npulses,1)
                     #range2 = numpy.arange(0,startPulseCountId,1)
                     #return range1, range2
                     zero = 0
                     npulse = max(self.pulseCount[0,:]+1)-1
                     looking_index = numpy.where(self.pulseCount.value[idrecord,:]==npulse)[0]
                     getLastIndex = looking_index[-1]
                     index_data = numpy.arange(0,getLastIndex+1,1)
                     index_buffer = numpy.arange(getLastIndex+1,self.npulses,1)
                     return index_data, index_buffer
             class AMISRReader(ProcessingUnit):
                 path = None
                 startDate = None
                 endDate = None
                 startTime = None
                 endTime = None
                 walk = None
                 isConfig = False
                 def __init__(self):
                     self.set = None
                     self.subset = None
                     self.extension_file = '.h5'
                     self.dtc_str = 'dtc'
                     self.dtc_id = 0
                     self.status = True
                     self.isConfig = False
                     self.dirnameList = []
                     self.filenameList = []
                     self.fileIndex = None
                     self.flagNoMoreFiles = False
                     self.flagIsNewFile = 0
                     self.filename = ''
                     self.amisrFilePointer = None
                     self.radacHeaderObj = None
                     self.dataOut = self.__createObjByDefault()
                     self.datablock = None
                     self.rest_datablock = None
                     self.range = None
                     self.idrecord_count = 0
                     self.profileIndex = 0
                     self.index_amisr_sample = None
                     self.index_amisr_buffer = None
                     self.beamCodeByFrame = None
                     self.radacTimeByFrame = None
                     #atributos originales tal y como esta en el archivo de datos
                     self.beamCodesFromFile = None
                     self.radacTimeFromFile = None
                     self.rangeFromFile = None
                     self.dataByFrame = None
                     self.dataset = None
                     self.beamCodeDict = {}
                     self.beamRangeDict = {}
                     #experiment cgf file
                     self.npulsesint_fromfile = None
                     self.recordsperfile_fromfile = None
                     self.nbeamcodes_fromfile = None
                     self.ngates_fromfile = None
                     self.ippSeconds_fromfile = None
                     self.frequency_h5file = None
                     self.__firstFile = True
                     self.buffer_radactime = None
                     self.index4_schain_datablock = None
                     self.index4_buffer = None
                     self.schain_datablock = None
                     self.buffer = None
                     self.linear_pulseCount = None
                     self.npulseByFrame = None
                     self.profileIndex_offset = None
                     self.timezone = 'ut'
                     self.__waitForNewFile = 20
                     self.__filename_online = None
                 def __createObjByDefault(self):
                     dataObj = AMISR()
                     return dataObj
                 def __setParameters(self,path='', startDate='',endDate='',startTime='', endTime='', walk=''):
                     self.path = path
                     self.startDate = startDate
                     self.endDate = endDate
                     self.startTime = startTime
                     self.endTime = endTime
                     self.walk = walk
                 def __checkPath(self):
                     if os.path.exists(self.path):
                         self.status = 1
                     else:
                         self.status = 0
                         print 'Path:%s does not exists'%self.path
                     return
                 def __selDates(self, amisr_dirname_format):
                     try:
                         year = int(amisr_dirname_format[0:4])
                         month = int(amisr_dirname_format[4:6])
                         dom = int(amisr_dirname_format[6:8])
                         thisDate = datetime.date(year,month,dom)
                         if (thisDate>=self.startDate and thisDate <= self.endDate):
                             return amisr_dirname_format
                     except:
                         return None
                 def __findDataForDates(self,online=False):
                     if not(self.status):
                         return None
                     pat = '\d+.\d+'
                     dirnameList = [re.search(pat,x) for x in os.listdir(self.path)]
                     dirnameList = filter(lambda x:x!=None,dirnameList)
                     dirnameList = [x.string for x in dirnameList]
                     if not(online):
                         dirnameList = [self.__selDates(x) for x in dirnameList]
                         dirnameList = filter(lambda x:x!=None,dirnameList)
                     if len(dirnameList)>0:
                         self.status = 1
                         self.dirnameList = dirnameList
                         self.dirnameList.sort()
                     else:
                         self.status = 0
                         return None
                 def __getTimeFromData(self):
                     startDateTime_Reader = datetime.datetime.combine(self.startDate,self.startTime)
                     endDateTime_Reader = datetime.datetime.combine(self.endDate,self.endTime)
                     print 'Filtering Files from %s to %s'%(startDateTime_Reader, endDateTime_Reader)
                     print '........................................'
                     filter_filenameList = []
                     self.filenameList.sort()
                     for i in range(len(self.filenameList)-1):
                         filename = self.filenameList[i]
                         fp = h5py.File(filename,'r')
                         time_str = fp.get('Time/RadacTimeString')
                         startDateTimeStr_File = time_str[0][0].split('.')[0]
                         junk = time.strptime(startDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
                         startDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
                         endDateTimeStr_File = time_str[-1][-1].split('.')[0]
                         junk = time.strptime(endDateTimeStr_File, '%Y-%m-%d %H:%M:%S')
                         endDateTime_File = datetime.datetime(junk.tm_year,junk.tm_mon,junk.tm_mday,junk.tm_hour, junk.tm_min, junk.tm_sec)
                         fp.close()
                         if self.timezone == 'lt':
                             startDateTime_File = startDateTime_File - datetime.timedelta(minutes = 300)
                             endDateTime_File = endDateTime_File - datetime.timedelta(minutes = 300)
                         if (endDateTime_File>=startDateTime_Reader and endDateTime_File<endDateTime_Reader):
                             #self.filenameList.remove(filename)
                             filter_filenameList.append(filename)
                     filter_filenameList.sort()
                     self.filenameList = filter_filenameList
                     return 1
                 def __filterByGlob1(self, dirName):
                     filter_files = glob.glob1(dirName, '*.*%s'%self.extension_file)
                     filterDict = {}
                     filterDict.setdefault(dirName)
                     filterDict[dirName] = filter_files
                     return filterDict
                 def __getFilenameList(self, fileListInKeys, dirList):
                     for value in fileListInKeys:
                         dirName = value.keys()[0]
                         for file in value[dirName]:
                             filename = os.path.join(dirName, file)
                             self.filenameList.append(filename)
                 def __selectDataForTimes(self, online=False):
                     #aun no esta implementado el filtro for tiempo
                     if not(self.status):
                         return None
                     dirList = [os.path.join(self.path,x) for x in self.dirnameList]
                     fileListInKeys = [self.__filterByGlob1(x) for x in dirList]
                     self.__getFilenameList(fileListInKeys, dirList)
                     if not(online):
                         #filtro por tiempo
                         if not(self.all):
                             self.__getTimeFromData()
                         if len(self.filenameList)>0:
                             self.status = 1
                             self.filenameList.sort()
                         else:
                             self.status = 0
                             return None
                     else:
                         #get the last file - 1
                         self.filenameList = [self.filenameList[-2]]
                     new_dirnameList = []
                     for dirname in self.dirnameList:
                         junk = numpy.array([dirname in x for x in self.filenameList])
                         junk_sum = junk.sum()
                         if junk_sum > 0:
                             new_dirnameList.append(dirname)
                     self.dirnameList = new_dirnameList
                     return 1
                 def __searchFilesOnline(self,
                                         path,
                                         walk=True):
                     startDate = datetime.datetime.utcnow().date()
                     endDate = datetime.datetime.utcnow().date()
                     self.__setParameters(path=path, startDate=startDate, endDate=endDate, walk=walk)
                     self.__checkPath()
                     self.__findDataForDates(online=True)
                     self.dirnameList = [self.dirnameList[-1]]
                     self.__selectDataForTimes(online=True)
                     return
                 def __searchFilesOffline(self,
                                         path,
                                         startDate,
                                         endDate,
                                         startTime=datetime.time(0,0,0),
                                         endTime=datetime.time(23,59,59),
                                         walk=True):
                     self.__setParameters(path, startDate, endDate, startTime, endTime, walk)
                     self.__checkPath()
                     self.__findDataForDates()
                     self.__selectDataForTimes()
                     for i in range(len(self.filenameList)):
                         print "%s" %(self.filenameList[i])
                     return
                 def __setNextFileOffline(self):
                     idFile = self.fileIndex
                     while (True):
                         idFile += 1
                         if not(idFile < len(self.filenameList)):
                             self.flagNoMoreFiles = 1
                             print "No more Files"
                             return 0
                         filename = self.filenameList[idFile]
                         amisrFilePointer = h5py.File(filename,'r')
                         break
                     self.flagIsNewFile = 1
                     self.fileIndex = idFile
                     self.filename = filename
                     self.amisrFilePointer = amisrFilePointer
                     print "Setting the file: %s"%self.filename
                     return 1
                 def __setNextFileOnline(self):
                     filename = self.filenameList[0]
                     if self.__filename_online != None:
                         self.__selectDataForTimes(online=True)
                         filename = self.filenameList[0]
                         while self.__filename_online == filename:
                             print 'waiting %d seconds to get a new file...'%(self.__waitForNewFile)
                             sleep(self.__waitForNewFile)
                             self.__selectDataForTimes(online=True)
                             filename = self.filenameList[0]
                     self.__filename_online = filename
                     self.amisrFilePointer = h5py.File(filename,'r')
                     self.flagIsNewFile = 1
                     self.filename = filename
                     print "Setting the file: %s"%self.filename
                     return 1
                 def __readHeader(self):
                     self.radacHeaderObj = RadacHeader(self.amisrFilePointer)
                     #update values from experiment cfg file
                     if self.radacHeaderObj.nrecords == self.recordsperfile_fromfile:
                         self.radacHeaderObj.nrecords = self.recordsperfile_fromfile
                     self.radacHeaderObj.nbeams = self.nbeamcodes_fromfile
                     self.radacHeaderObj.npulses = self.npulsesint_fromfile
                     self.radacHeaderObj.nsamples = self.ngates_fromfile
                     #looking index list for data
                     start_index = self.radacHeaderObj.pulseCount[0,:][0]
                     end_index = self.radacHeaderObj.npulses
                     range4data = range(start_index, end_index)
                     self.index4_schain_datablock = numpy.array(range4data)
                     buffer_start_index = 0
                     buffer_end_index = self.radacHeaderObj.pulseCount[0,:][0]
                     range4buffer = range(buffer_start_index, buffer_end_index)
                     self.index4_buffer = numpy.array(range4buffer)
                     self.linear_pulseCount = numpy.array(range4data + range4buffer)
                     self.npulseByFrame = max(self.radacHeaderObj.pulseCount[0,:]+1)
                     #get tuning frequency
                     frequency_h5file_dataset = self.amisrFilePointer.get('Rx'+'/TuningFrequency')
                     self.frequency_h5file = frequency_h5file_dataset[0,0]
                     self.flagIsNewFile = 1
                 def __getBeamCode(self):
                     self.beamCodeDict = {}
                     self.beamRangeDict = {}
                     beamCodeMap = self.amisrFilePointer.get('Setup/BeamcodeMap')
                     for i in range(len(self.radacHeaderObj.beamCode[0,:])):
                         self.beamCodeDict.setdefault(i)
                         self.beamRangeDict.setdefault(i)
                         beamcodeValue = self.radacHeaderObj.beamCode[0,i]
                         beamcodeIndex = numpy.where(beamCodeMap[:,0] == beamcodeValue)[0][0]
                         x = beamCodeMap[beamcodeIndex][1]
                         y = beamCodeMap[beamcodeIndex][2]
                         z = beamCodeMap[beamcodeIndex][3]
                         self.beamCodeDict[i] = [beamcodeValue, x, y, z]
                     just4record0 = self.radacHeaderObj.beamCodeByPulse[0,:]
                     for i in range(len(self.beamCodeDict.values())):
                         xx = numpy.where(just4record0==self.beamCodeDict.values()[i][0])
                         indexPulseByBeam = self.linear_pulseCount[xx[0]]
                         self.beamRangeDict[i] = indexPulseByBeam
                 def __getExpParameters(self):
                     if not(self.status):
                         return None
                     experimentCfgPath = os.path.join(self.path, self.dirnameList[0], 'Setup')
                     expFinder = glob.glob1(experimentCfgPath,'*.exp')
                     if len(expFinder)== 0:
                         self.status = 0
                         return None
                     experimentFilename = os.path.join(experimentCfgPath,expFinder[0])
                     f = open(experimentFilename)
                     lines = f.readlines()
                     f.close()
                     parmsList = ['npulsesint*','recordsperfile*','nbeamcodes*','ngates*']
                     filterList = [fnmatch.filter(lines, x) for x in parmsList]
                     values = [re.sub(r'\D',"",x[0]) for x in filterList]
                     self.npulsesint_fromfile = int(values[0])
                     self.recordsperfile_fromfile = int(values[1])
                     self.nbeamcodes_fromfile = int(values[2])
                     self.ngates_fromfile = int(values[3])
                     tufileFinder = fnmatch.filter(lines, 'tufile=*')
                     tufile = tufileFinder[0].split('=')[1].split('\n')[0]
                     tufile = tufile.split('\r')[0]
                     tufilename = os.path.join(experimentCfgPath,tufile)
                     f = open(tufilename)
                     lines = f.readlines()
                     f.close()
                     self.ippSeconds_fromfile = float(lines[1].split()[2])/1E6
                     self.status = 1
                 def __setIdsAndArrays(self):
                     self.dataByFrame = self.__setDataByFrame()
                     self.beamCodeByFrame = self.amisrFilePointer.get('Raw11/Data/RadacHeader/BeamCode').value[0, :]
                     self.readRanges()
                     self.index_amisr_sample, self.index_amisr_buffer = self.radacHeaderObj.getIndexRangeToPulse(0)
                     self.radacTimeByFrame = numpy.zeros(self.radacHeaderObj.npulses)
                     if len(self.index_amisr_buffer) > 0:
                         self.buffer_radactime = numpy.zeros_like(self.radacTimeByFrame)
                 def __setNextFile(self,online=False):
                     if not(online):
                         newFile = self.__setNextFileOffline()
                     else:
                         newFile = self.__setNextFileOnline()
                     if not(newFile):
                         return 0
                     self.__readHeader()
                     if self.__firstFile:
                         self.__setIdsAndArrays()
                         self.__firstFile = False
                     self.__getBeamCode()
                     self.readDataBlock()
                 def setup(self,path=None,
                                 startDate=None,
                                 endDate=None,
                                 startTime=datetime.time(0,0,0),
                                 endTime=datetime.time(23,59,59),
                                 walk=True,
                                 timezone='ut',
                                 all=0,
                                 online=False):
                     self.timezone = timezone
                     self.all = all
                     self.online = online
                     if not(online):
                         #Busqueda de archivos offline
                         self.__searchFilesOffline(path, startDate, endDate, startTime, endTime, walk)
                     else:
                         self.__searchFilesOnline(path, walk)
                     if not(self.filenameList):
                         print "There is no files into the folder: %s"%(path)
                         sys.exit(-1)
                     self.__getExpParameters()
                     self.fileIndex = -1
                     self.__setNextFile(online)
             #         first_beamcode = self.radacHeaderObj.beamCodeByPulse[0,0]
             #         index = numpy.where(self.radacHeaderObj.beamCodeByPulse[0,:]!=first_beamcode)[0][0]
                     self.profileIndex_offset = self.radacHeaderObj.pulseCount[0,:][0]
                     self.profileIndex = self.profileIndex_offset
                 def readRanges(self):
                     dataset = self.amisrFilePointer.get('Raw11/Data/Samples/Range')
                     self.rangeFromFile = numpy.reshape(dataset.value,(-1))
                     return self.rangeFromFile
                 def readRadacTime(self,idrecord, range1, range2):
                     self.radacTimeFromFile = self.radacHeaderObj.radacTime.value
                     radacTimeByFrame = numpy.zeros((self.radacHeaderObj.npulses))
                     #radacTimeByFrame = dataset[idrecord - 1,range1]
                     #radacTimeByFrame = dataset[idrecord,range2]
                     return radacTimeByFrame
                 def readBeamCode(self, idrecord, range1, range2):
                     dataset = self.amisrFilePointer.get('Raw11/Data/RadacHeader/BeamCode')
                     beamcodeByFrame = numpy.zeros((self.radacHeaderObj.npulses))
                     self.beamCodesFromFile = dataset.value
                     #beamcodeByFrame[range1] = dataset[idrecord - 1, range1]
                     #beamcodeByFrame[range2] = dataset[idrecord, range2]
                     beamcodeByFrame[range1] = dataset[idrecord, range1]
                     beamcodeByFrame[range2] = dataset[idrecord, range2]
                     return beamcodeByFrame
                 def __setDataByFrame(self):
                     ndata = 2 # porque es complejo
                     dataByFrame = numpy.zeros((self.radacHeaderObj.npulses, self.radacHeaderObj.nsamples, ndata))
                     return dataByFrame
                 def __readDataSet(self):
                     dataset = self.amisrFilePointer.get('Raw11/Data/Samples/Data')
                     return dataset
                 def __setDataBlock(self,):
                     real = self.dataByFrame[:,:,0] #asumo que 0 es real
                     imag = self.dataByFrame[:,:,1] #asumo que 1 es imaginario
                     datablock = real + imag*1j #armo el complejo
                     return datablock
                 def readSamples_version1(self,idrecord):
                     #estas tres primeras lineas solo se deben ejecutar una vez
                     if self.flagIsNewFile:
                         #reading dataset
                         self.dataset = self.__readDataSet()
                         self.flagIsNewFile = 0
                     if idrecord == 0:
                         self.dataByFrame[self.index4_schain_datablock, : ,:] = self.dataset[0, self.index_amisr_sample,:,:]
                         self.radacTimeByFrame[self.index4_schain_datablock] = self.radacHeaderObj.radacTime[0, self.index_amisr_sample]
                         datablock = self.__setDataBlock()
                         if len(self.index_amisr_buffer) > 0:
                             self.buffer = self.dataset[0, self.index_amisr_buffer,:,:]
                             self.buffer_radactime = self.radacHeaderObj.radacTime[0, self.index_amisr_buffer]
                         return datablock
                     if len(self.index_amisr_buffer) > 0:
                         self.dataByFrame[self.index4_buffer,:,:] = self.buffer.copy()
                         self.radacTimeByFrame[self.index4_buffer] = self.buffer_radactime.copy()
                     self.dataByFrame[self.index4_schain_datablock,:,:] = self.dataset[idrecord, self.index_amisr_sample,:,:]
                     self.radacTimeByFrame[self.index4_schain_datablock] = self.radacHeaderObj.radacTime[idrecord, self.index_amisr_sample]
                     datablock = self.__setDataBlock()
                     if len(self.index_amisr_buffer) > 0:
                         self.buffer = self.dataset[idrecord, self.index_amisr_buffer, :, :]
                         self.buffer_radactime = self.radacHeaderObj.radacTime[idrecord, self.index_amisr_buffer]
                     return datablock
                 def readSamples(self,idrecord):
                     if self.flagIsNewFile:
                         self.dataByFrame = self.__setDataByFrame()
                         self.beamCodeByFrame = self.amisrFilePointer.get('Raw11/Data/RadacHeader/BeamCode').value[idrecord, :]
                         #reading ranges
                         self.readRanges()
                         #reading dataset
                         self.dataset = self.__readDataSet()
                     self.flagIsNewFile = 0
                     self.radacTimeByFrame = self.radacHeaderObj.radacTime.value[idrecord, :]
                     self.dataByFrame = self.dataset[idrecord, :, :, :]
                     datablock = self.__setDataBlock()
                     return datablock
                 def readDataBlock(self):
                     self.datablock = self.readSamples_version1(self.idrecord_count)
                     #self.datablock = self.readSamples(self.idrecord_count)
                     #print 'record:', self.idrecord_count
                     self.idrecord_count += 1
                     self.profileIndex = 0
                     if self.idrecord_count >= self.radacHeaderObj.nrecords:
                         self.idrecord_count = 0
                         self.flagIsNewFile = 1
                 def readNextBlock(self):
                     self.readDataBlock()
                     if self.flagIsNewFile:
                         self.__setNextFile(self.online)
                     pass
                 def __hasNotDataInBuffer(self):
                     #self.radacHeaderObj.npulses debe ser otra variable para considerar el numero de pulsos a tomar en el primer y ultimo record
                     if self.profileIndex >= self.radacHeaderObj.npulses:
                         return 1
                     return 0
                 def printUTC(self):
                     print self.dataOut.utctime
                     print ''
                 def setObjProperties(self):
                     self.dataOut.heightList = self.rangeFromFile/1000.0 #km
                     self.dataOut.nProfiles = self.radacHeaderObj.npulses
                     self.dataOut.nRecords = self.radacHeaderObj.nrecords
                     self.dataOut.nBeams = self.radacHeaderObj.nbeams
                     self.dataOut.ippSeconds = self.ippSeconds_fromfile
             #         self.dataOut.timeInterval = self.dataOut.ippSeconds * self.dataOut.nCohInt
                     self.dataOut.frequency = self.frequency_h5file
                     self.dataOut.npulseByFrame = self.npulseByFrame
                     self.dataOut.nBaud = None
                     self.dataOut.nCode = None
                     self.dataOut.code = None
                     self.dataOut.beamCodeDict = self.beamCodeDict
                     self.dataOut.beamRangeDict = self.beamRangeDict
                     if self.timezone == 'lt':
                         self.dataOut.timeZone = time.timezone / 60. #get the timezone in minutes
                     else:
                         self.dataOut.timeZone = 0 #by default time is UTC
                 def getData(self):
                     if self.flagNoMoreFiles:
                         self.dataOut.flagNoData = True
                         print 'Process finished'
                         return 0
                     if self.__hasNotDataInBuffer():
                         self.readNextBlock()
-                    if self.datablock == None: # setear esta condicion cuando no hayan datos por leers
+                    if self.datablock is None: # setear esta condicion cuando no hayan datos por leers
                         self.dataOut.flagNoData = True
                         return 0
                     self.dataOut.data = numpy.reshape(self.datablock[self.profileIndex,:],(1,-1))
                     self.dataOut.utctime = self.radacTimeByFrame[self.profileIndex]
                     self.dataOut.profileIndex = self.profileIndex
                     self.dataOut.flagNoData = False
                     self.profileIndex += 1
                     return self.dataOut.data
                 def run(self, **kwargs):
                     if not(self.isConfig):
                         self.setup(**kwargs)
                         self.setObjProperties()
                         self.isConfig = True
                     self.getData()

schainpy/model/io/jroIO_voltage.py +1 -1

             '''
             Created on Jul 2, 2014
             @author: roj-idl71
             '''
             import numpy
             from jroIO_base import LOCALTIME, JRODataReader, JRODataWriter
             from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation
             from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader, SystemHeader, RadarControllerHeader, ProcessingHeader
             from schainpy.model.data.jrodata import Voltage
             class VoltageReader(JRODataReader, ProcessingUnit):
                 """
                 Esta clase permite leer datos de voltage desde archivos en formato rawdata (.r). La lectura
                 de los datos siempre se realiza por bloques. Los datos leidos (array de 3 dimensiones:
                 perfiles*alturas*canales) son almacenados en la variable "buffer".
                                          perfiles * alturas * canales
                 Esta clase contiene instancias (objetos) de las clases BasicHeader, SystemHeader,
                 RadarControllerHeader y Voltage. Los tres primeros se usan para almacenar informacion de la
                 cabecera de datos (metadata), y el cuarto (Voltage) para obtener y almacenar un perfil de
                 datos desde el "buffer" cada vez que se ejecute el metodo "getData".
                 Example:
                     dpath = "/home/myuser/data"
                     startTime = datetime.datetime(2010,1,20,0,0,0,0,0,0)
                     endTime = datetime.datetime(2010,1,21,23,59,59,0,0,0)
                     readerObj = VoltageReader()
                     readerObj.setup(dpath, startTime, endTime)
                     while(True):
                         #to get one profile
                         profile =  readerObj.getData()
                         #print the profile
                         print profile
                         #If you want to see all datablock
                         print readerObj.datablock
                         if readerObj.flagNoMoreFiles:
                             break
                 """
                 ext = ".r"
                 optchar = "D"
                 dataOut = None
                 def __init__(self):
                     """
                     Inicializador de la clase VoltageReader para la lectura de datos de voltage.
                     Input:
                         dataOut    :    Objeto de la clase Voltage. Este objeto sera utilizado para
                                           almacenar un perfil de datos cada vez que se haga un requerimiento
                                           (getData). El perfil sera obtenido a partir del buffer de datos,
                                           si el buffer esta vacio se hara un nuevo proceso de lectura de un
                                           bloque de datos.
                                           Si este parametro no es pasado se creara uno internamente.
                     Variables afectadas:
                         self.dataOut
                     Return:
                         None
                     """
                     ProcessingUnit.__init__(self)
                     self.isConfig = False
                     self.datablock = None
                     self.utc = 0
                     self.ext = ".r"
                     self.optchar = "D"
                     self.basicHeaderObj = BasicHeader(LOCALTIME)
                     self.systemHeaderObj = SystemHeader()
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.processingHeaderObj = ProcessingHeader()
                     self.online = 0
                     self.fp = None
                     self.idFile = None
                     self.dtype = None
                     self.fileSizeByHeader = None
                     self.filenameList = []
                     self.filename = None
                     self.fileSize = None
                     self.firstHeaderSize = 0
                     self.basicHeaderSize = 24
                     self.pathList = []
                     self.filenameList = []
                     self.lastUTTime = 0
                     self.maxTimeStep = 30
                     self.flagNoMoreFiles = 0
                     self.set = 0
                     self.path = None
                     self.profileIndex = 2**32-1
                     self.delay  = 3   #seconds
                     self.nTries  = 3  #quantity tries
                     self.nFiles = 3   #number of files for searching
                     self.nReadBlocks = 0
                     self.flagIsNewFile = 1
                     self.__isFirstTimeOnline = 1
             #         self.ippSeconds = 0
                     self.flagDiscontinuousBlock = 0
                     self.flagIsNewBlock = 0
                     self.nTotalBlocks = 0
                     self.blocksize = 0
                     self.dataOut = self.createObjByDefault()
                     self.nTxs = 1
                     self.txIndex = 0
                 def createObjByDefault(self):
                     dataObj = Voltage()
                     return dataObj
                 def __hasNotDataInBuffer(self):
                     if self.profileIndex >= self.processingHeaderObj.profilesPerBlock:
                         return 1
                     return 0
                 def getBlockDimension(self):
                     """
                     Obtiene la cantidad de puntos a leer por cada bloque de datos
                     Affected:
                         self.blocksize
                     Return:
                         None
                     """
                     pts2read = self.processingHeaderObj.profilesPerBlock * self.processingHeaderObj.nHeights * self.systemHeaderObj.nChannels
                     self.blocksize = pts2read
                 def readBlock(self):
                     """
                     readBlock lee el bloque de datos desde la posicion actual del puntero del archivo
                     (self.fp) y actualiza todos los parametros relacionados al bloque de datos
                     (metadata + data). La data leida es almacenada en el buffer y el contador del buffer
                     es seteado a 0
                     Inputs:
                         None
                     Return:
                         None
                     Affected:
                         self.profileIndex
                         self.datablock
                         self.flagIsNewFile
                         self.flagIsNewBlock
                         self.nTotalBlocks
                     Exceptions:
                         Si un bloque leido no es un bloque valido
                     """
                     current_pointer_location = self.fp.tell()
                     junk = numpy.fromfile( self.fp, self.dtype, self.blocksize )
                     try:
                         junk = junk.reshape( (self.processingHeaderObj.profilesPerBlock, self.processingHeaderObj.nHeights, self.systemHeaderObj.nChannels) )
                     except:
                         #print "The read block (%3d) has not enough data" %self.nReadBlocks
                         if self.waitDataBlock(pointer_location=current_pointer_location):
                             junk = numpy.fromfile( self.fp, self.dtype, self.blocksize )
                             junk = junk.reshape( (self.processingHeaderObj.profilesPerBlock, self.processingHeaderObj.nHeights, self.systemHeaderObj.nChannels) )
             #             return 0
                     junk = numpy.transpose(junk, (2,0,1))
                     self.datablock = junk['real'] + junk['imag']*1j
                     self.profileIndex = 0
                     self.flagIsNewFile = 0
                     self.flagIsNewBlock = 1
                     self.nTotalBlocks += 1
                     self.nReadBlocks += 1
                     return 1
                 def getFirstHeader(self):
                     self.dataOut.systemHeaderObj = self.systemHeaderObj.copy()
                     self.dataOut.radarControllerHeaderObj = self.radarControllerHeaderObj.copy()
                     if self.nTxs > 1:
                         self.dataOut.radarControllerHeaderObj.ippSeconds = self.radarControllerHeaderObj.ippSeconds/self.nTxs
             #         self.dataOut.timeInterval = self.radarControllerHeaderObj.ippSeconds * self.processingHeaderObj.nCohInt
             #
             #         if self.radarControllerHeaderObj.code != None:
             #
             #             self.dataOut.nCode = self.radarControllerHeaderObj.nCode
             #
             #             self.dataOut.nBaud = self.radarControllerHeaderObj.nBaud
             #
             #             self.dataOut.code = self.radarControllerHeaderObj.code
                     self.dataOut.dtype = self.dtype
                     self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock*self.nTxs
                     if self.processingHeaderObj.nHeights % self.nTxs != 0:
                         raise ValueError, "nTxs (%d) should be a multiple of nHeights (%d)" %(self.nTxs, self.processingHeaderObj.nHeights)
                     xf = self.processingHeaderObj.firstHeight + int(self.processingHeaderObj.nHeights/self.nTxs)*self.processingHeaderObj.deltaHeight
                     self.dataOut.heightList = numpy.arange(self.processingHeaderObj.firstHeight, xf, self.processingHeaderObj.deltaHeight)
                     self.dataOut.channelList = range(self.systemHeaderObj.nChannels)
                     self.dataOut.nCohInt = self.processingHeaderObj.nCohInt
                     self.dataOut.flagShiftFFT = False
                     self.dataOut.flagDecodeData = False #asumo q la data no esta decodificada
                     self.dataOut.flagDeflipData = False #asumo q la data no esta sin flip
                     self.dataOut.flagShiftFFT = False
                 def getData(self):
                     """
                     getData obtiene una unidad de datos del buffer de lectura, un perfil,  y la copia al objeto self.dataOut
                     del tipo "Voltage" con todos los parametros asociados a este (metadata). cuando no hay datos
                     en el buffer de lectura es necesario hacer una nueva lectura de los bloques de datos usando
                     "readNextBlock"
                     Ademas incrementa el contador del buffer "self.profileIndex" en 1.
                     Return:
                         Si el flag self.getByBlock ha sido seteado el bloque completo es copiado a self.dataOut y el self.profileIndex
                         es igual al total de perfiles leidos desde el archivo.
                         Si self.getByBlock == False:
                             self.dataOut.data = buffer[:, thisProfile, :]
                             shape = [nChannels, nHeis]
                         Si self.getByBlock == True:
                             self.dataOut.data = buffer[:, :, :]
                             shape = [nChannels, nProfiles, nHeis]
                     Variables afectadas:
                         self.dataOut
                         self.profileIndex
                     Affected:
                         self.dataOut
                         self.profileIndex
                         self.flagDiscontinuousBlock
                         self.flagIsNewBlock
                     """
                     if self.flagNoMoreFiles:
                         self.dataOut.flagNoData = True
                         print 'Process finished'
                         return 0
                     self.flagDiscontinuousBlock = 0
                     self.flagIsNewBlock = 0
                     if self.__hasNotDataInBuffer():
                         if not( self.readNextBlock() ):
                             return 0
                         self.getFirstHeader()
-                    if self.datablock == None:
+                    if self.datablock is None:
                         self.dataOut.flagNoData = True
                         return 0
                     if not self.getByBlock:
                         """
                         Return profile by profile
                         If nTxs > 1 then one profile is divided by nTxs and number of total
                         blocks is increased by nTxs (nProfiles *= nTxs)
                         """
                         self.dataOut.flagDataAsBlock = False
                         if self.nTxs == 1:
                             self.dataOut.data = self.datablock[:,self.profileIndex,:]
                             self.dataOut.profileIndex = self.profileIndex
                             self.profileIndex += 1
                         else:
                             iniHei_ForThisTx = (self.txIndex)*int(self.processingHeaderObj.nHeights/self.nTxs)
                             endHei_ForThisTx = (self.txIndex+1)*int(self.processingHeaderObj.nHeights/self.nTxs)
             #                 print iniHei_ForThisTx, endHei_ForThisTx
                             self.dataOut.data = self.datablock[:, self.profileIndex, iniHei_ForThisTx:endHei_ForThisTx]
                             self.dataOut.profileIndex = self.profileIndex*self.nTxs + self.txIndex
                             self.txIndex += 1
                             if self.txIndex == self.nTxs:
                                 self.txIndex = 0
                                 self.profileIndex += 1
                     else:
                         """
                         Return all block
                         """
                         self.dataOut.flagDataAsBlock = True
                         self.dataOut.data = self.datablock
                         self.dataOut.profileIndex = self.processingHeaderObj.profilesPerBlock
                         self.profileIndex = self.processingHeaderObj.profilesPerBlock
                     self.dataOut.flagNoData = False
                     self.getBasicHeader()
                     self.dataOut.realtime = self.online
                     return self.dataOut.data
             class VoltageWriter(JRODataWriter, Operation):
                 """
                 Esta clase permite escribir datos de voltajes a archivos procesados (.r). La escritura
                 de los datos siempre se realiza por bloques.
                 """
                 ext = ".r"
                 optchar = "D"
                 shapeBuffer = None
                 def __init__(self):
                     """
                     Inicializador de la clase VoltageWriter para la escritura de datos de espectros.
                     Affected:
                         self.dataOut
                     Return: None
                     """
                     Operation.__init__(self)
                     self.nTotalBlocks = 0
                     self.profileIndex = 0
                     self.isConfig = False
                     self.fp = None
                     self.flagIsNewFile = 1
                     self.nTotalBlocks = 0
                     self.flagIsNewBlock = 0
                     self.setFile = None
                     self.dtype = None
                     self.path = None
                     self.filename = None
                     self.basicHeaderObj = BasicHeader(LOCALTIME)
                     self.systemHeaderObj = SystemHeader()
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.processingHeaderObj = ProcessingHeader()
                 def hasAllDataInBuffer(self):
                     if self.profileIndex >= self.processingHeaderObj.profilesPerBlock:
                         return 1
                     return 0
                 def setBlockDimension(self):
                     """
                     Obtiene las formas dimensionales del los subbloques de datos que componen un bloque
                     Affected:
                         self.shape_spc_Buffer
                         self.shape_cspc_Buffer
                         self.shape_dc_Buffer
                     Return: None
                     """
                     self.shapeBuffer = (self.processingHeaderObj.profilesPerBlock,
                                         self.processingHeaderObj.nHeights,
                                         self.systemHeaderObj.nChannels)
                     self.datablock = numpy.zeros((self.systemHeaderObj.nChannels,
                                                  self.processingHeaderObj.profilesPerBlock,
                                                  self.processingHeaderObj.nHeights),
                                                  dtype=numpy.dtype('complex64'))
                 def writeBlock(self):
                     """
                     Escribe el buffer en el file designado
                     Affected:
                         self.profileIndex
                         self.flagIsNewFile
                         self.flagIsNewBlock
                         self.nTotalBlocks
                         self.blockIndex
                     Return: None
                     """
                     data = numpy.zeros( self.shapeBuffer, self.dtype )
                     junk = numpy.transpose(self.datablock, (1,2,0))
                     data['real'] = junk.real
                     data['imag'] = junk.imag
                     data = data.reshape( (-1) )
                     data.tofile( self.fp )
                     self.datablock.fill(0)
                     self.profileIndex = 0
                     self.flagIsNewFile = 0
                     self.flagIsNewBlock = 1
                     self.blockIndex += 1
                     self.nTotalBlocks += 1
             #         print "[Writing] Block = %04d" %self.blockIndex
                 def putData(self):
                     """
                     Setea un bloque de datos y luego los escribe en un file
                     Affected:
                         self.flagIsNewBlock
                         self.profileIndex
                     Return:
 :    Si no hay data o no hay mas files que puedan escribirse
 :    Si se escribio la data de un bloque en un file
                     """
                     if self.dataOut.flagNoData:
                         return 0
                     self.flagIsNewBlock = 0
                     if self.dataOut.flagDiscontinuousBlock:
                         self.datablock.fill(0)
                         self.profileIndex = 0
                         self.setNextFile()
                     if self.profileIndex == 0:
                         self.setBasicHeader()
                     self.datablock[:,self.profileIndex,:] = self.dataOut.data
                     self.profileIndex += 1
                     if self.hasAllDataInBuffer():
                         #if self.flagIsNewFile:
                         self.writeNextBlock()
             #            self.setFirstHeader()
                     return 1
                 def __getProcessFlags(self):
                     processFlags = 0
                     dtype0 = numpy.dtype([('real','<i1'),('imag','<i1')])
                     dtype1 = numpy.dtype([('real','<i2'),('imag','<i2')])
                     dtype2 = numpy.dtype([('real','<i4'),('imag','<i4')])
                     dtype3 = numpy.dtype([('real','<i8'),('imag','<i8')])
                     dtype4 = numpy.dtype([('real','<f4'),('imag','<f4')])
                     dtype5 = numpy.dtype([('real','<f8'),('imag','<f8')])
                     dtypeList = [dtype0, dtype1, dtype2, dtype3, dtype4, dtype5]
                     datatypeValueList =  [PROCFLAG.DATATYPE_CHAR,
                                        PROCFLAG.DATATYPE_SHORT,
                                        PROCFLAG.DATATYPE_LONG,
                                        PROCFLAG.DATATYPE_INT64,
                                        PROCFLAG.DATATYPE_FLOAT,
                                        PROCFLAG.DATATYPE_DOUBLE]
                     for index in range(len(dtypeList)):
                         if self.dataOut.dtype == dtypeList[index]:
                             dtypeValue = datatypeValueList[index]
                             break
                     processFlags += dtypeValue
                     if self.dataOut.flagDecodeData:
                         processFlags += PROCFLAG.DECODE_DATA
                     if self.dataOut.flagDeflipData:
                         processFlags += PROCFLAG.DEFLIP_DATA
                     if self.dataOut.code != None:
                         processFlags += PROCFLAG.DEFINE_PROCESS_CODE
                     if self.dataOut.nCohInt > 1:
                         processFlags += PROCFLAG.COHERENT_INTEGRATION
                     return processFlags
                 def __getBlockSize(self):
                     '''
                     Este metodos determina el cantidad de bytes para un bloque de datos de tipo Voltage
                     '''
                     dtype0 = numpy.dtype([('real','<i1'),('imag','<i1')])
                     dtype1 = numpy.dtype([('real','<i2'),('imag','<i2')])
                     dtype2 = numpy.dtype([('real','<i4'),('imag','<i4')])
                     dtype3 = numpy.dtype([('real','<i8'),('imag','<i8')])
                     dtype4 = numpy.dtype([('real','<f4'),('imag','<f4')])
                     dtype5 = numpy.dtype([('real','<f8'),('imag','<f8')])
                     dtypeList = [dtype0, dtype1, dtype2, dtype3, dtype4, dtype5]
                     datatypeValueList = [1,2,4,8,4,8]
                     for index in range(len(dtypeList)):
                         if self.dataOut.dtype == dtypeList[index]:
                             datatypeValue = datatypeValueList[index]
                             break
                     blocksize = int(self.dataOut.nHeights * self.dataOut.nChannels * self.profilesPerBlock * datatypeValue * 2)
                     return blocksize
                 def setFirstHeader(self):
                     """
                     Obtiene una copia del First Header
                     Affected:
                         self.systemHeaderObj
                         self.radarControllerHeaderObj
                         self.dtype
                     Return:
                         None
                     """
                     self.systemHeaderObj = self.dataOut.systemHeaderObj.copy()
                     self.systemHeaderObj.nChannels = self.dataOut.nChannels
                     self.radarControllerHeaderObj = self.dataOut.radarControllerHeaderObj.copy()
                     self.setBasicHeader()
                     processingHeaderSize = 40 # bytes
                     self.processingHeaderObj.dtype = 0 # Voltage
                     self.processingHeaderObj.blockSize = self.__getBlockSize()
                     self.processingHeaderObj.profilesPerBlock = self.profilesPerBlock
                     self.processingHeaderObj.dataBlocksPerFile = self.blocksPerFile
                     self.processingHeaderObj.nWindows = 1 #podria ser 1 o self.dataOut.processingHeaderObj.nWindows
                     self.processingHeaderObj.processFlags = self.__getProcessFlags()
                     self.processingHeaderObj.nCohInt = self.dataOut.nCohInt
                     self.processingHeaderObj.nIncohInt = 1 # Cuando la data de origen es de tipo Voltage
                     self.processingHeaderObj.totalSpectra = 0 # Cuando la data de origen es de tipo Voltage
             #         if self.dataOut.code != None:
             #             self.processingHeaderObj.code = self.dataOut.code
             #             self.processingHeaderObj.nCode = self.dataOut.nCode
             #             self.processingHeaderObj.nBaud = self.dataOut.nBaud
             #             codesize = int(8 + 4 * self.dataOut.nCode * self.dataOut.nBaud)
             #             processingHeaderSize += codesize
                     if self.processingHeaderObj.nWindows != 0:
                         self.processingHeaderObj.firstHeight = self.dataOut.heightList[0]
                         self.processingHeaderObj.deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
                         self.processingHeaderObj.nHeights = self.dataOut.nHeights
                         self.processingHeaderObj.samplesWin = self.dataOut.nHeights
                         processingHeaderSize += 12
                     self.processingHeaderObj.size = processingHeaderSize

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

             '''
             $Author: murco $
             $Id: jroproc_base.py 1 2012-11-12 18:56:07Z murco $
             '''
             class ProcessingUnit(object):
                 """
                 Esta es la clase base para el procesamiento de datos.
                 Contiene el metodo "call" para llamar operaciones. Las operaciones pueden ser:
                     - Metodos internos (callMethod)
                     - Objetos del tipo Operation (callObject). Antes de ser llamados, estos objetos
                       tienen que ser agreagados con el metodo "add".
                 """
                 # objeto de datos de entrada (Voltage, Spectra o Correlation)
                 dataIn = None
                 dataInList = []
                 # objeto de datos de entrada (Voltage, Spectra o Correlation)
                 dataOut = None
                 operations2RunDict = None
                 isConfig = False
                 def __init__(self):
                     self.dataIn = None
                     self.dataInList = []
                     self.dataOut = None
                     self.operations2RunDict = {}
                     self.isConfig = False
                 def addOperation(self, opObj, objId):
                     """
                     Agrega un objeto del tipo "Operation" (opObj) a la lista de objetos "self.objectList" y retorna el
                     identificador asociado a este objeto.
                     Input:
                         object    :    objeto de la clase "Operation"
                     Return:
                         objId    :    identificador del objeto, necesario para ejecutar la operacion
                     """
                     self.operations2RunDict[objId] = opObj
                     return objId
                 def getOperationObj(self, objId):
                     if objId not in self.operations2RunDict.keys():
                         return None
                     return self.operations2RunDict[objId]
                 def operation(self, **kwargs):
                     """
                     Operacion directa sobre la data (dataOut.data). Es necesario actualizar los valores de los
                     atributos del objeto dataOut
                     Input:
                         **kwargs    :    Diccionario de argumentos de la funcion a ejecutar
                     """
                     raise ValueError, "ImplementedError"
                 def callMethod(self, name, **kwargs):
                     """
                     Ejecuta el metodo con el nombre "name" y con argumentos **kwargs de la propia clase.
                     Input:
                         name        :    nombre del metodo a ejecutar
                         **kwargs     :    diccionario con los nombres y valores de la funcion a ejecutar.
                     """
                     #Checking the inputs
                     if name == 'run':
                         if not self.checkInputs():
                             self.dataOut.flagNoData = True
                             return False
                     else:
                         #Si no es un metodo RUN la entrada es la misma dataOut (interna)
                         if self.dataOut.isEmpty():
                             return False
                     #Getting the pointer to method
                     methodToCall = getattr(self, name)
                     #Executing the self method
                     methodToCall(**kwargs)
                     #Checkin the outputs
             #         if name == 'run':
             #             pass
             #         else:
             #             pass
             #
             #         if name != 'run':
             #             return True
-                    if self.dataOut == None:
+                    if self.dataOut is None:
                         return False
                     if self.dataOut.isEmpty():
                         return False
                     return True
                 def callObject(self, objId, **kwargs):
                     """
                     Ejecuta la operacion asociada al identificador del objeto "objId"
                     Input:
                         objId        :    identificador del objeto a ejecutar
                         **kwargs    :    diccionario con los nombres y valores de la funcion a ejecutar.
                     Return:
                         None
                     """
                     if self.dataOut.isEmpty():
                         return False
                     externalProcObj = self.operations2RunDict[objId]
                     externalProcObj.run(self.dataOut, **kwargs)
                     return True
                 def call(self, opType, opName=None, opId=None, **kwargs):
                     """
                     Return True si ejecuta la operacion interna nombrada "opName" o la operacion externa
                     identificada con el id "opId"; con los argumentos "**kwargs".
                     False si la operacion no se ha ejecutado.
                     Input:
                         opType    :    Puede ser "self" o "external"
                             La operacion puede ser de dos tipos (callMethod or callObject):
 . Un metodo propio de esta clase:
                                 opType = "self"
 . El metodo "run" de un objeto del tipo Operation o de un derivado de ella:
                                 opType = "other" or "external".
                         opName    : Si la operacion es interna (opType = 'self'), entonces el "opName" sera
                                     usada para llamar a un metodo interno de la clase Processing
                         opId    :    Si la operacion es externa (opType = 'other'), entonces el "opId" sera
                                     usada para llamar al metodo "run" de la clase Operation registrada con ese Id
                     Exception:
                            Este objeto de tipo Operation debe de haber sido agregado antes con el metodo:
                            "addOperation" e identificado con el valor "opId"  = el id de la operacion.
                            De lo contrario retornara un error del tipo IOError
                     """
                     if opType == 'self':
                         if not opName:
                             raise IOError, "opName parameter should be defined"
                         sts = self.callMethod(opName, **kwargs)
                     if opType == 'other' or opType == 'external':
                         if not opId:
                             raise IOError, "opId parameter should be defined"
                         if opId not in self.operations2RunDict.keys():
                             raise IOError, "This id operation have not been registered"
                         sts = self.callObject(opId, **kwargs)
                     return sts
                 def setInput(self, dataIn):
                     self.dataIn = dataIn
                     self.dataInList.append(dataIn)
                 def getOutputObj(self):
                     return self.dataOut
                 def checkInputs(self):
                     for thisDataIn in self.dataInList:
                         if thisDataIn.isEmpty():
                             return False
                     return True
                 def setup(self):
                     raise ValueError, "Not implemented"
                 def run(self):
                     raise ValueError, "Not implemented"
                 def close(self):
                     #Close every thread, queue or any other object here is it is neccesary.
                     return
             class Operation(object):
                 """
                 Clase base para definir las operaciones adicionales que se pueden agregar a la clase ProcessingUnit
                 y necesiten acumular informacion previa de los datos a procesar. De preferencia usar un buffer de
                 acumulacion dentro de esta clase
                 Ejemplo: Integraciones coherentes, necesita la informacion previa de los n perfiles anteriores (bufffer)
                 """
                 __buffer = None
                 isConfig = False
                 def __init__(self):
                     self.__buffer = None
                     self.isConfig = False
                 def setup(self):
                     self.isConfig = True
                     raise ValueError, "Not implemented"
                 def run(self, dataIn, **kwargs):
                     """
                     Realiza las operaciones necesarias sobre la dataIn.data y actualiza los
                     atributos del objeto dataIn.
                     Input:
                         dataIn    :    objeto del tipo JROData
                     Return:
                         None
                     Affected:
                         __buffer    :    buffer de recepcion de datos.
                     """
                     if not self.isConfig:
                         self.setup(**kwargs)
                     raise ValueError, "ImplementedError"
                 def close(self):
                     pass

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

             import numpy
             from jroproc_base import ProcessingUnit, Operation
             from schainpy.model.data.jrodata import Correlation
             class CorrelationProc(ProcessingUnit):
                 def __init__(self):
                     ProcessingUnit.__init__(self)
                     self.objectDict = {}
                     self.buffer = None
                     self.firstdatatime = None
                     self.profIndex = 0
                     self.dataOut = Correlation()
                 def __updateObjFromVoltage(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.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.nHeights = self.dataIn.nHeights
             #         self.dataOut.nChannels = self.dataIn.nChannels
                     self.dataOut.nBaud = self.dataIn.nBaud
                     self.dataOut.nCode = self.dataIn.nCode
                     self.dataOut.code = self.dataIn.code
             #        self.dataOut.nProfiles = self.dataOut.nFFTPoints
                     self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
                     self.dataOut.utctime = self.firstdatatime
                     self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
                     self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
             #        self.dataOut.nCohInt = self.dataIn.nCohInt
             #        self.dataOut.nIncohInt = 1
                     self.dataOut.ippSeconds = self.dataIn.ippSeconds
             #        self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
             #         self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nPoints
                 def removeDC(self, jspectra):
                     nChannel = jspectra.shape[0]
                     for i in range(nChannel):
                         jspectra_tmp = jspectra[i,:,:]
                         jspectra_DC = numpy.mean(jspectra_tmp,axis = 0)
                         jspectra_tmp = jspectra_tmp - jspectra_DC
                         jspectra[i,:,:] = jspectra_tmp
                     return jspectra
                 def removeNoise(self, mode = 2):
                     indR = numpy.where(self.dataOut.lagR == 0)[0][0]
                     indT = numpy.where(self.dataOut.lagT == 0)[0][0]
                     jspectra = self.dataOut.data_corr[:,:,indR,:]
                     num_chan = jspectra.shape[0]
                     num_hei = jspectra.shape[2]
                     freq_dc = indT
                     ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
                     NPot = self.dataOut.getNoise(mode)
                     jspectra[:,freq_dc,:] = jspectra[:,freq_dc,:] - NPot
                     SPot = jspectra[:,freq_dc,:]
                     pairsAutoCorr = self.dataOut.getPairsAutoCorr()
             #         self.dataOut.signalPotency = SPot
                     self.dataOut.noise = NPot
                     self.dataOut.SNR = (SPot/NPot)[pairsAutoCorr]
                     self.dataOut.data_corr[:,:,indR,:] = jspectra
                     return 1
                 def calculateNormFactor(self):
                     pairsList = self.dataOut.pairsList
                     pairsAutoCorr = self.dataOut.pairsAutoCorr
                     nHeights = self.dataOut.nHeights
                     nPairs = len(pairsList)
                     normFactor = numpy.zeros((nPairs,nHeights))
                     indR = numpy.where(self.dataOut.lagR == 0)[0][0]
                     indT = numpy.where(self.dataOut.lagT == 0)[0][0]
                     for l in range(len(pairsList)):
                         firstChannel = pairsList[l][0]
                         secondChannel = pairsList[l][1]
                         AC1 = pairsAutoCorr[firstChannel]
                         AC2 = pairsAutoCorr[secondChannel]
                         if (AC1 >= 0 and AC2 >= 0):
                             data1 = numpy.abs(self.dataOut.data_corr[AC1,:,indR,:])
                             data2 = numpy.abs(self.dataOut.data_corr[AC2,:,indR,:])
                             maxim1 = data1.max(axis = 0)
                             maxim2 = data1.max(axis = 0)
                             maxim = numpy.sqrt(maxim1*maxim2)
                         else:
                             #In case there is no autocorrelation for the pair
                             data = numpy.abs(self.dataOut.data_corr[l,:,indR,:])
                             maxim = numpy.max(data, axis = 0)
                         normFactor[l,:] = maxim
                     self.dataOut.normFactor = normFactor
                     return 1
                 def run(self, lagT=None, lagR=None, pairsList=None,
                          nPoints=None, nAvg=None, bufferSize=None,
                          fullT = False, fullR = False, removeDC = False):
                     self.dataOut.flagNoData = True
                     if self.dataIn.type == "Correlation":
                         self.dataOut.copy(self.dataIn)
                         return
                     if self.dataIn.type == "Voltage":
                         if pairsList == None:
                             pairsList = [numpy.array([0,0])]
                         if nPoints == None:
                             nPoints = 128
                         #------------------------------------------------------------
                         #Condicionales para calcular Correlaciones en Tiempo y Rango
                         if fullT:
                             lagT = numpy.arange(nPoints*2 - 1) - nPoints + 1
                         elif lagT == None:
                             lagT = numpy.array([0])
                         else:
                             lagT = numpy.array(lagT)
                         if fullR:
                             lagR = numpy.arange(self.dataOut.nHeights)
                         elif lagR == None:
                             lagR = numpy.array([0])
                         #-------------------------------------------------------------
                         if nAvg == None:
                             nAvg = 1
                         if bufferSize == None:
                             bufferSize = 0
                         deltaH = self.dataIn.heightList[1] - self.dataIn.heightList[0]
                         self.dataOut.lagR = numpy.round(numpy.array(lagR)/deltaH)
                         self.dataOut.pairsList = pairsList
                         self.dataOut.nPoints = nPoints
             #             channels = numpy.sort(list(set(list(itertools.chain.from_iterable(pairsList)))))
-                        if self.buffer == None:
+                        if self.buffer is None:
                             self.buffer = numpy.zeros((self.dataIn.nChannels,self.dataIn.nProfiles,self.dataIn.nHeights),dtype='complex')
                         self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()[:,:]
                         self.profIndex += 1
                         if self.firstdatatime == None:
                             self.firstdatatime = self.dataIn.utctime
                         if self.profIndex == nPoints:
                             tmp = self.buffer[:,0:nPoints,:]
                             self.buffer = None
                             self.buffer = tmp
                             #---------------    Remover DC    ------------
                             if removeDC:
                                 self.buffer = self.removeDC(self.buffer)
                             #---------------------------------------------
                             self.dataOut.data_volts = self.buffer
                             self.__updateObjFromVoltage()
                             self.dataOut.data_corr = numpy.zeros((len(pairsList),
                                                              len(lagT),len(lagR),
                                                              self.dataIn.nHeights),
                                                             dtype='complex')
                             for l in range(len(pairsList)):
                                 firstChannel = pairsList[l][0]
                                 secondChannel = pairsList[l][1]
                                 tmp = None
                                 tmp = numpy.zeros((len(lagT),len(lagR),self.dataIn.nHeights),dtype='complex')
                                 for t in range(len(lagT)):
                                     for r in range(len(lagR)):
                                         idxT = lagT[t]
                                         idxR = lagR[r]
                                         if idxT >= 0:
                                             vStacked = numpy.vstack((self.buffer[secondChannel,idxT:,:],
                                                                      numpy.zeros((idxT,self.dataIn.nHeights),dtype='complex')))
                                         else:
                                             vStacked = numpy.vstack((numpy.zeros((-idxT,self.dataIn.nHeights),dtype='complex'),
                                                                      self.buffer[secondChannel,:(nPoints + idxT),:]))
                                         if idxR >= 0:
                                             hStacked = numpy.hstack((vStacked[:,idxR:],numpy.zeros((nPoints,idxR),dtype='complex')))
                                         else:
                                             hStacked = numpy.hstack((numpy.zeros((nPoints,-idxR),dtype='complex'),vStacked[:,(self.dataOut.nHeights + idxR)]))
                                         tmp[t,r,:] = numpy.sum((numpy.conjugate(self.buffer[firstChannel,:,:])*hStacked),axis=0)
                                         hStacked = None
                                         vStacked = None
                                 self.dataOut.data_corr[l,:,:,:] = tmp[:,:,:]
                             #Se Calcula los factores de Normalizacion
                             self.dataOut.pairsAutoCorr = self.dataOut.getPairsAutoCorr()
                             self.dataOut.lagT = lagT*self.dataIn.ippSeconds*self.dataIn.nCohInt
                             self.dataOut.lagR = lagR
                             self.calculateNormFactor()
                             self.dataOut.flagNoData = False
                             self.buffer = None
                             self.firstdatatime = None
                             self.profIndex = 0
                             return

schainpy/model/proc/jroproc_heispectra.py +2 -2

             import numpy
             from jroproc_base import ProcessingUnit, Operation
             from schainpy.model.data.jrodata import SpectraHeis
             class SpectraHeisProc(ProcessingUnit):
                 def __init__(self):
                     ProcessingUnit.__init__(self)
             #        self.buffer = None
             #        self.firstdatatime = None
             #        self.profIndex = 0
                     self.dataOut = SpectraHeis()
                 def __updateObjFromVoltage(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.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 = self.dataIn.dtype
                     self.dataOut.dtype = numpy.dtype([('real','<f4'),('imag','<f4')])
             #        self.dataOut.nHeights = self.dataIn.nHeights
             #        self.dataOut.nChannels = self.dataIn.nChannels
                     self.dataOut.nBaud = self.dataIn.nBaud
                     self.dataOut.nCode = self.dataIn.nCode
                     self.dataOut.code = self.dataIn.code
             #        self.dataOut.nProfiles = 1
             #        self.dataOut.nProfiles = self.dataOut.nFFTPoints
                     self.dataOut.nFFTPoints = self.dataIn.nHeights
             #        self.dataOut.channelIndexList = self.dataIn.channelIndexList
             #        self.dataOut.flagNoData = self.dataIn.flagNoData
                     self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
                     self.dataOut.utctime = self.dataIn.utctime
             #        self.dataOut.utctime = self.firstdatatime
                     self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
                     self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
             #        self.dataOut.flagShiftFFT = self.dataIn.flagShiftFFT
                     self.dataOut.nCohInt = self.dataIn.nCohInt
                     self.dataOut.nIncohInt = 1
             #         self.dataOut.ippSeconds= self.dataIn.ippSeconds
                     self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
             #         self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nIncohInt
             #        self.dataOut.set=self.dataIn.set
             #        self.dataOut.deltaHeight=self.dataIn.deltaHeight
                 def __updateObjFromFits(self):
                     self.dataOut.utctime = self.dataIn.utctime
                     self.dataOut.channelIndexList = self.dataIn.channelIndexList
                     self.dataOut.channelList = self.dataIn.channelList
                     self.dataOut.heightList = self.dataIn.heightList
                     self.dataOut.data_spc = self.dataIn.data
                     self.dataOut.ippSeconds = self.dataIn.ippSeconds
                     self.dataOut.nCohInt = self.dataIn.nCohInt
                     self.dataOut.nIncohInt = self.dataIn.nIncohInt
             #         self.dataOut.timeInterval = self.dataIn.timeInterval
                     self.dataOut.timeZone = self.dataIn.timeZone
                     self.dataOut.useLocalTime = True
             #         self.dataOut.
             #         self.dataOut.
                 def __getFft(self):
                     fft_volt = numpy.fft.fft(self.dataIn.data, axis=1)
                     fft_volt = numpy.fft.fftshift(fft_volt,axes=(1,))
                     spc = numpy.abs(fft_volt * numpy.conjugate(fft_volt))/(self.dataOut.nFFTPoints)
                     self.dataOut.data_spc = spc
                 def run(self):
                     self.dataOut.flagNoData = True
                     if self.dataIn.type == "Fits":
                         self.__updateObjFromFits()
                         self.dataOut.flagNoData = False
                         return
                     if self.dataIn.type == "SpectraHeis":
                         self.dataOut.copy(self.dataIn)
                         return
                     if self.dataIn.type == "Voltage":
                         self.__updateObjFromVoltage()
                         self.__getFft()
                         self.dataOut.flagNoData = False
                         return
                     raise ValueError, "The type object %s is not valid"%(self.dataIn.type)
                 def selectChannels(self, channelList):
                     channelIndexList = []
                     for channel in channelList:
                         index = self.dataOut.channelList.index(channel)
                         channelIndexList.append(index)
                     self.selectChannelsByIndex(channelIndexList)
                 def selectChannelsByIndex(self, 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:
                         self.dataOut.data
                         self.dataOut.channelIndexList
                         self.dataOut.nChannels
                         self.dataOut.m_ProcessingHeader.totalSpectra
                         self.dataOut.systemHeaderObj.numChannels
                         self.dataOut.m_ProcessingHeader.blockSize
                     Return:
                         None
                     """
                     for channelIndex in channelIndexList:
                         if channelIndex not in self.dataOut.channelIndexList:
                             print channelIndexList
                             raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
             #         nChannels = len(channelIndexList)
                     data_spc = self.dataOut.data_spc[channelIndexList,:]
                     self.dataOut.data_spc = data_spc
                     self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
                     return 1
             class IncohInt4SpectraHeis(Operation):
                 isConfig = False
                 __profIndex = 0
                 __withOverapping  = False
                 __byTime = False
                 __initime = None
                 __lastdatatime = None
                 __integrationtime = None
                 __buffer = None
                 __dataReady = False
                 n = None
                 def __init__(self):
                     Operation.__init__(self)
             #         self.isConfig = False
                 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 = None
                     self.__dataReady = False
                     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.__withOverapping = True
                         self.__buffer = None
                     else:
                         self.__withOverapping = 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.__withOverapping:
                         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 == None:
+                    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.__withOverapping:
                         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)
                     n = self.__profIndex
                     return data, n
                 def byProfiles(self, data):
                     self.__dataReady = False
                     avgdata = None
             #         n = None
                     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 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 == None:
+                    if avgdata is None:
                         return None, None
                     avgdatatime = self.__initime
                     deltatime = datatime -self.__lastdatatime
                     if not self.__withOverapping:
                         self.__initime = datatime
                     else:
                         self.__initime += deltatime
                     return avgdata, avgdatatime
                 def run(self, dataOut, **kwargs):
                     if not self.isConfig:
                         self.setup(**kwargs)
                         self.isConfig = True
                     avgdata, avgdatatime = self.integrate(dataOut.data_spc, dataOut.utctime)
             #        dataOut.timeInterval *= n
                     dataOut.flagNoData = True
                     if self.__dataReady:
                         dataOut.data_spc = avgdata
                         dataOut.nIncohInt *= self.n
             #            dataOut.nCohInt *= self.n
                         dataOut.utctime = avgdatatime
             #             dataOut.timeInterval = dataOut.ippSeconds * dataOut.nIncohInt
             #            dataOut.timeInterval = self.__timeInterval*self.n
                         dataOut.flagNoData = False

schainpy/model/proc/jroproc_parameters.py +5 -5

             import numpy
             import math
             from scipy import optimize
             from scipy import interpolate
             from scipy import signal
             from scipy import stats
             import re
             import datetime
             import copy
             import sys
             import importlib
             import itertools
             from jroproc_base import ProcessingUnit, Operation
             from schainpy.model.data.jrodata import Parameters
             class ParametersProc(ProcessingUnit):
                 nSeconds = None
                 def __init__(self):
                     ProcessingUnit.__init__(self)
             #         self.objectDict = {}
                     self.buffer = None
                     self.firstdatatime = None
                     self.profIndex = 0
                     self.dataOut = Parameters()
                 def __updateObjFromInput(self):
                     self.dataOut.inputUnit = self.dataIn.type
                     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.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.nHeights = self.dataIn.nHeights
             #         self.dataOut.nChannels = self.dataIn.nChannels
                     self.dataOut.nBaud = self.dataIn.nBaud
                     self.dataOut.nCode = self.dataIn.nCode
                     self.dataOut.code = self.dataIn.code
             #        self.dataOut.nProfiles = self.dataOut.nFFTPoints
                     self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
                     self.dataOut.utctime = self.firstdatatime
                     self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
                     self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
             #        self.dataOut.nCohInt = self.dataIn.nCohInt
             #        self.dataOut.nIncohInt = 1
                     self.dataOut.ippSeconds = self.dataIn.ippSeconds
             #        self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
                     self.dataOut.timeInterval = self.dataIn.timeInterval
                     self.dataOut.heightList = self.dataIn.getHeiRange()
                     self.dataOut.frequency = self.dataIn.frequency
                 def run(self, nSeconds = None, nProfiles = None):
                     if self.firstdatatime == None:
                         self.firstdatatime = self.dataIn.utctime
                     #----------------------    Voltage Data    ---------------------------
                     if self.dataIn.type == "Voltage":
                         self.dataOut.flagNoData = True
                         if nSeconds != None:
                             self.nSeconds = nSeconds
                             self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt)))
-                        if self.buffer == None:
+                        if self.buffer is None:
                                 self.buffer = numpy.zeros((self.dataIn.nChannels,
                                                            self.nProfiles,
                                                            self.dataIn.nHeights),
                                                            dtype='complex')
                         self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
                         self.profIndex += 1
                         if self.profIndex == self.nProfiles:
                             self.__updateObjFromInput()
                             self.dataOut.data_pre = self.buffer.copy()
                             self.dataOut.paramInterval = nSeconds
                             self.dataOut.flagNoData = False
                             self.buffer = None
                             self.firstdatatime = None
                             self.profIndex = 0
                             return
                     #----------------------    Spectra Data    ---------------------------
                     if self.dataIn.type == "Spectra":
                         self.dataOut.data_pre = self.dataIn.data_spc.copy()
                         self.dataOut.abscissaList = self.dataIn.getVelRange(1)
                         self.dataOut.noise = self.dataIn.getNoise()
                         self.dataOut.normFactor = self.dataIn.normFactor
                         self.dataOut.groupList = self.dataIn.pairsList
                         self.dataOut.flagNoData = False
                     #----------------------    Correlation Data    ---------------------------
                     if self.dataIn.type == "Correlation":
                         lagRRange = self.dataIn.lagR
                         indR = numpy.where(lagRRange == 0)[0][0]
                         self.dataOut.data_pre = self.dataIn.data_corr.copy()[:,:,indR,:]
                         self.dataOut.abscissaList = self.dataIn.getLagTRange(1)
                         self.dataOut.noise = self.dataIn.noise
                         self.dataOut.normFactor = self.dataIn.normFactor
                         self.dataOut.data_SNR = self.dataIn.SNR
                         self.dataOut.groupList = self.dataIn.pairsList
                         self.dataOut.flagNoData = False
                             #----------------------    Correlation Data    ---------------------------
                     if self.dataIn.type == "Parameters":
                         self.dataOut.copy(self.dataIn)
                         self.dataOut.flagNoData = False
                         return True
                     self.__updateObjFromInput()
                     self.firstdatatime = None
                     self.dataOut.utctimeInit = self.dataIn.utctime
                     self.dataOut.outputInterval = self.dataIn.timeInterval
                 #-------------------    Get Moments    ----------------------------------
                 def GetMoments(self, channelList = None):
                     '''
                     Function GetMoments()
                     Input:
                         channelList    :    simple channel list to select e.g. [2,3,7]
                         self.dataOut.data_pre
                         self.dataOut.abscissaList
                         self.dataOut.noise
                     Affected:
                         self.dataOut.data_param
                         self.dataOut.data_SNR
                     '''
                     data = self.dataOut.data_pre
                     absc = self.dataOut.abscissaList[:-1]
                     noise = self.dataOut.noise
                     data_param = numpy.zeros((data.shape[0], 4, data.shape[2]))
                     if channelList== None:
                         channelList = self.dataIn.channelList
                     self.dataOut.channelList = channelList
                     for ind in channelList:
                         data_param[ind,:,:] = self.__calculateMoments(data[ind,:,:], absc, noise[ind])
                     self.dataOut.data_param = data_param[:,1:,:]
                     self.dataOut.data_SNR = data_param[:,0]
                     return
                 def __calculateMoments(self, oldspec, oldfreq, n0, nicoh = None, graph = None, smooth = None, type1 = None, fwindow = None, snrth = None, dc = None, aliasing = None, oldfd = None, wwauto = None):
                     if (nicoh == None): nicoh = 1
                     if (graph == None): graph = 0
                     if (smooth == None): smooth = 0
                     elif (self.smooth < 3): smooth = 0
                     if (type1 == None): type1 = 0
-                    if (fwindow == None): fwindow = numpy.zeros(oldfreq.size) + 1
+                    if (fwindow is None): fwindow = numpy.zeros(oldfreq.size) + 1
                     if (snrth == None): snrth = -3
                     if (dc == None): dc = 0
                     if (aliasing == None): aliasing = 0
                     if (oldfd == None): oldfd = 0
                     if (wwauto == None): wwauto = 0
                     if (n0 < 1.e-20):   n0 = 1.e-20
                     freq = oldfreq
                     vec_power = numpy.zeros(oldspec.shape[1])
                     vec_fd = numpy.zeros(oldspec.shape[1])
                     vec_w = numpy.zeros(oldspec.shape[1])
                     vec_snr = numpy.zeros(oldspec.shape[1])
                     for ind in range(oldspec.shape[1]):
                         spec = oldspec[:,ind]
                         aux = spec*fwindow
                         max_spec = aux.max()
                         m = list(aux).index(max_spec)
                         #Smooth
                         if (smooth == 0):   spec2 = spec
                         else:   spec2 = scipy.ndimage.filters.uniform_filter1d(spec,size=smooth)
                         #    Calculo de Momentos
                         bb = spec2[range(m,spec2.size)]
                         bb = (bb<n0).nonzero()
                         bb = bb[0]
                         ss = spec2[range(0,m + 1)]
                         ss = (ss<n0).nonzero()
                         ss = ss[0]
                         if (bb.size == 0):
                             bb0 = spec.size - 1 - m
                         else:
                             bb0 = bb[0] - 1
                             if (bb0 < 0):
                                 bb0 = 0
                         if (ss.size == 0):   ss1 = 1
                         else: ss1 = max(ss) + 1
                         if (ss1 > m):   ss1 = m
                         valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1
                         power = ((spec2[valid] - n0)*fwindow[valid]).sum()
                         fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power
                         w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power)
                         snr = (spec2.mean()-n0)/n0
                         if (snr < 1.e-20) :
                             snr = 1.e-20
                         vec_power[ind] = power
                         vec_fd[ind] = fd
                         vec_w[ind] = w
                         vec_snr[ind] = snr
                     moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w))
                     return moments
                 #------------------    Get SA Parameters    --------------------------
                 def GetSAParameters(self):
                     pairslist = self.dataOut.groupList
                     num_pairs = len(pairslist)
                     vel = self.dataOut.abscissaList
                     spectra = self.dataOut.data_pre
                     cspectra = self.dataIn.data_cspc
                     delta_v = vel[1] - vel[0]
                     #Calculating the power spectrum
                     spc_pow = numpy.sum(spectra, 3)*delta_v
                     #Normalizing Spectra
                     norm_spectra = spectra/spc_pow
                     #Calculating the norm_spectra at peak
                     max_spectra = numpy.max(norm_spectra, 3)
                     #Normalizing Cross Spectra
                     norm_cspectra = numpy.zeros(cspectra.shape)
                     for i in range(num_chan):
                         norm_cspectra[i,:,:] = cspectra[i,:,:]/numpy.sqrt(spc_pow[pairslist[i][0],:]*spc_pow[pairslist[i][1],:])
                     max_cspectra = numpy.max(norm_cspectra,2)
                     max_cspectra_index = numpy.argmax(norm_cspectra, 2)
                     for i in range(num_pairs):
                         cspc_par[i,:,:] = __calculateMoments(norm_cspectra)
                 #-------------------    Get Lags    ----------------------------------
                 def GetLags(self):
                     '''
                     Function GetMoments()
                     Input:
                         self.dataOut.data_pre
                         self.dataOut.abscissaList
                         self.dataOut.noise
                         self.dataOut.normFactor
                         self.dataOut.data_SNR
                         self.dataOut.groupList
                         self.dataOut.nChannels
                     Affected:
                         self.dataOut.data_param
                     '''
                     data = self.dataOut.data_pre
                     normFactor = self.dataOut.normFactor
                     nHeights = self.dataOut.nHeights
                     absc = self.dataOut.abscissaList[:-1]
                     noise = self.dataOut.noise
                     SNR = self.dataOut.data_SNR
                     pairsList = self.dataOut.groupList
                     nChannels = self.dataOut.nChannels
                     pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
                     self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights))
                     dataNorm = numpy.abs(data)
                     for l in range(len(pairsList)):
                         dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:]
                     self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc)
                     self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, pairsAutoCorr, absc)
                     return
                 def __getPairsAutoCorr(self, pairsList, nChannels):
                     pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
                     for l in range(len(pairsList)):
                         firstChannel = pairsList[l][0]
                         secondChannel = pairsList[l][1]
                         #Obteniendo pares de Autocorrelacion
                         if firstChannel == secondChannel:
                             pairsAutoCorr[firstChannel] = int(l)
                     pairsAutoCorr = pairsAutoCorr.astype(int)
                     pairsCrossCorr = range(len(pairsList))
                     pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
                     return pairsAutoCorr, pairsCrossCorr
                 def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange):
                     Pt0 = data.shape[1]/2
                     #Funcion de Autocorrelacion
                     dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0)
                     #Obtencion Indice de TauCross
                     indCross = data[pairsCrossCorr,:,:].argmax(axis = 1)
                     #Obtencion Indice de TauAuto
                     indAuto = numpy.zeros(indCross.shape,dtype = 'int')
                     CCValue = data[pairsCrossCorr,Pt0,:]
                     for i in range(pairsCrossCorr.size):
                         indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0)
                     #Obtencion de TauCross y TauAuto
                     tauCross = lagTRange[indCross]
                     tauAuto  = lagTRange[indAuto]
                     Nan1, Nan2 = numpy.where(tauCross == lagTRange[0])
                     tauCross[Nan1,Nan2] = numpy.nan
                     tauAuto[Nan1,Nan2] = numpy.nan
                     tau = numpy.vstack((tauCross,tauAuto))
                     return tau
                 def __calculateLag1Phase(self, data, pairs, lagTRange):
                     data1 = stats.nanmean(data[pairs,:,:], axis = 0)
                     lag1 = numpy.where(lagTRange == 0)[0][0] + 1
                     phase = numpy.angle(data1[lag1,:])
                     return phase
                     #-------------------    Detect Meteors  ------------------------------
                 def MeteorDetection(self, hei_ref = None, tauindex = 0,
                                   predefinedPhaseShifts = None, centerReceiverIndex = 2, saveAll = False,
                                   cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25,
                                   noise_timeStep = 4, noise_multiple = 4,
                                   multDet_timeLimit = 1, multDet_rangeLimit = 3,
                                   phaseThresh = 20, SNRThresh = 8,
                                   hmin = 70, hmax=110, azimuth = 0) :
                     '''
                     Function DetectMeteors()
                         Project developed with paper:
                         HOLDSWORTH ET AL. 2004
                     Input:
                         self.dataOut.data_pre
                         centerReceiverIndex:      From the channels, which is the center receiver
                         hei_ref:                  Height reference for the Beacon signal extraction
                         tauindex:
                         predefinedPhaseShifts:    Predefined phase offset for the voltge signals
                         cohDetection:             Whether to user Coherent detection or not
                         cohDet_timeStep:          Coherent Detection calculation time step
                         cohDet_thresh:            Coherent Detection phase threshold to correct phases
                         noise_timeStep:           Noise calculation time step
                         noise_multiple:           Noise multiple to define signal threshold
                         multDet_timeLimit:        Multiple Detection Removal time limit in seconds
                         multDet_rangeLimit:       Multiple Detection Removal range limit in km
                         phaseThresh:              Maximum phase difference between receiver to be consider a meteor
                         SNRThresh:                Minimum SNR threshold of the meteor signal to be consider a meteor
                         hmin:                     Minimum Height of the meteor to use it in the further wind estimations
                         hmax:                     Maximum Height of the meteor to use it in the further wind estimations
                         azimuth:                  Azimuth angle correction
                     Affected:
                         self.dataOut.data_param
                     Rejection Criteria (Errors):
 : No error; analysis OK
 : SNR < SNR threshold
 : angle of arrival (AOA) ambiguously determined
 : AOA estimate not feasible
 : Large difference in AOAs obtained from different antenna baselines
 : echo at start or end of time series
 : echo less than 5 examples long; too short for analysis
 : echo rise exceeds 0.3s
 : echo decay time less than twice rise time
 : large power level before echo
 : large power level after echo
 : poor fit to amplitude for estimation of decay time
 : poor fit to CCF phase variation for estimation of radial drift velocity
 : height unresolvable echo: not valid height within 70 to 110 km
 : height ambiguous echo: more then one possible height within 70 to 110 km
 : radial drift velocity or projected horizontal velocity exceeds 200 m/s
 : oscilatory echo, indicating event most likely not an underdense echo
 : phase difference in meteor Reestimation
                     Data Storage:
                         Meteors for Wind Estimation   (8):
                         Day    Hour    |    Range    Height
                         Azimuth    Zenith    errorCosDir
                         VelRad    errorVelRad
                         TypeError
                      '''
                     #Get Beacon signal
                     newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex])
                     if hei_ref != None:
                         newheis = numpy.where(self.dataOut.heightList>hei_ref)
                     heiRang = self.dataOut.getHeiRange()
                     #Pairs List
                     pairslist = []
                     nChannel = self.dataOut.nChannels
                     for i in range(nChannel):
                         if i != centerReceiverIndex:
                             pairslist.append((centerReceiverIndex,i))
                     #****************REMOVING HARDWARE PHASE DIFFERENCES***************
                     # see if the user put in pre defined phase shifts
                     voltsPShift = self.dataOut.data_pre.copy()
                     if predefinedPhaseShifts != None:
                         hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180
                     elif beaconPhaseShifts:
                         #get hardware phase shifts using beacon signal
                         hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10)
                         hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0)
                     else:
                         hardwarePhaseShifts = numpy.zeros(5)
                     voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex')
                     for i in range(self.dataOut.data_pre.shape[0]):
                         voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i])
                     #******************END OF REMOVING HARDWARE PHASE DIFFERENCES*********
                     #Remove DC
                     voltsDC = numpy.mean(voltsPShift,1)
                     voltsDC = numpy.mean(voltsDC,1)
                     for i in range(voltsDC.shape[0]):
                         voltsPShift[i] = voltsPShift[i] - voltsDC[i]
                     #Don't considerate last heights, theyre used to calculate Hardware Phase Shift
                     voltsPShift = voltsPShift[:,:,:newheis[0][0]]
                     #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) **********
                     #Coherent Detection
                     if cohDetection:
                         #use coherent detection to get the net power
                         cohDet_thresh = cohDet_thresh*numpy.pi/180
                         voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist, cohDet_thresh)
                     #Non-coherent detection!
                     powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0)
                     #********** END OF COH/NON-COH POWER CALCULATION**********************
                     #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS ****************
                     #Get noise
                     noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval)
             #         noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval)
                     #Get signal threshold
                     signalThresh = noise_multiple*noise
                     #Meteor echoes detection
                     listMeteors = self.__findMeteors(powerNet, signalThresh)
                     #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION **********
                     #************** REMOVE MULTIPLE DETECTIONS (3.5) ***************************
                     #Parameters
                     heiRange = self.dataOut.getHeiRange()
                     rangeInterval = heiRange[1] - heiRange[0]
                     rangeLimit = multDet_rangeLimit/rangeInterval
                     timeLimit = multDet_timeLimit/self.dataOut.timeInterval
                     #Multiple detection removals
                     listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit)
                     #************ END OF REMOVE MULTIPLE DETECTIONS **********************
                     #*********************     METEOR REESTIMATION  (3.7, 3.8, 3.9, 3.10)   ********************
                     #Parameters
                     phaseThresh = phaseThresh*numpy.pi/180
                     thresh = [phaseThresh, noise_multiple, SNRThresh]
                     #Meteor reestimation  (Errors N 1, 6, 12, 17)
                     listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency)
             #         listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise)
                     #Estimation of decay times (Errors N 7, 8, 11)
                     listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency)
                     #*******************     END OF METEOR REESTIMATION    *******************
                     #*********************    METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13)    **************************
                     #Calculating Radial Velocity (Error N 15)
                     radialStdThresh = 10
                     listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval)
                     if len(listMeteors4) > 0:
                         pairsList = []
                         pairx = (0,3)
                         pairy = (1,2)
                         pairsList.append(pairx)
                         pairsList.append(pairy)
                         #Setting New Array
                         date = repr(self.dataOut.datatime)
                         arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang)
                         meteorOps = MeteorOperations()
                         jph = numpy.array([0,0,0,0])
                         h = (hmin,hmax)
                         arrayParameters = meteorOps.getMeteorParams(arrayParameters, azimuth, h, pairsList, jph)
             #             #Calculate AOA (Error N 3, 4)
             #             #JONES ET AL. 1998
             #             error = arrayParameters[:,-1]
             #             AOAthresh = numpy.pi/8
             #             phases = -arrayParameters[:,9:13]
             #             arrayParameters[:,4:7], arrayParameters[:,-1] = meteorOps.getAOA(phases, pairsList, error, AOAthresh, azimuth)
             #
             #             #Calculate Heights (Error N 13 and 14)
             #             error = arrayParameters[:,-1]
             #             Ranges = arrayParameters[:,2]
             #             zenith = arrayParameters[:,5]
             #             arrayParameters[:,3], arrayParameters[:,-1] = meteorOps.getHeights(Ranges, zenith, error, hmin, hmax)
             #             error = arrayParameters[:,-1]
                     #*********************    END OF PARAMETERS CALCULATION    **************************
                     #***************************+     PASS DATA TO NEXT STEP    **********************
                         arrayFinal = arrayParameters.reshape((1,arrayParameters.shape[0],arrayParameters.shape[1]))
                         self.dataOut.data_param = arrayFinal
                     return
                 def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n):
                     minIndex = min(newheis[0])
                     maxIndex = max(newheis[0])
                     voltage = voltage0[:,:,minIndex:maxIndex+1]
                     nLength = voltage.shape[1]/n
                     nMin = 0
                     nMax = 0
                     phaseOffset = numpy.zeros((len(pairslist),n))
                     for i in range(n):
                         nMax += nLength
                         phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0]))
                         phaseCCF = numpy.mean(phaseCCF, axis = 2)
                         phaseOffset[:,i] = phaseCCF.transpose()
                         nMin = nMax
             #         phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist)
                     #Remove Outliers
                     factor = 2
                     wt = phaseOffset - signal.medfilt(phaseOffset,(1,5))
                     dw = numpy.std(wt,axis = 1)
                     dw = dw.reshape((dw.size,1))
                     ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor))
                     phaseOffset[ind] = numpy.nan
                     phaseOffset = stats.nanmean(phaseOffset, axis=1)
                     return phaseOffset
                 def __shiftPhase(self, data, phaseShift):
                     #this will shift the phase of a complex number
                     dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j)
                     return dataShifted
                 def __estimatePhaseDifference(self, array, pairslist):
                     nChannel = array.shape[0]
                     nHeights = array.shape[2]
                     numPairs = len(pairslist)
             #         phaseCCF = numpy.zeros((nChannel, 5, nHeights))
                     phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2]))
                     #Correct phases
                     derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:]
                     indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
                     if indDer[0].shape[0] > 0:
                         for i in range(indDer[0].shape[0]):
                             signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]])
                             phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi
             #         for j in range(numSides):
             #             phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2])
             #             phaseCCF[j,:,:] = numpy.angle(phaseCCFAux)
             #
                     #Linear
                     phaseInt = numpy.zeros((numPairs,1))
                     angAllCCF = phaseCCF[:,[0,1,3,4],0]
                     for j in range(numPairs):
                         fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:])
                         phaseInt[j] = fit[1]
                     #Phase Differences
                     phaseDiff = phaseInt - phaseCCF[:,2,:]
                     phaseArrival = phaseInt.reshape(phaseInt.size)
                     #Dealias
                     indAlias = numpy.where(phaseArrival > numpy.pi)
                     phaseArrival[indAlias] -= 2*numpy.pi
                     indAlias = numpy.where(phaseArrival < -numpy.pi)
                     phaseArrival[indAlias] += 2*numpy.pi
                     return phaseDiff, phaseArrival
                 def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh):
                     #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power
                     #find the phase shifts of each channel over 1 second intervals
                     #only look at ranges below the beacon signal
                     numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
                     numBlocks = int(volts.shape[1]/numProfPerBlock)
                     numHeights = volts.shape[2]
                     nChannel = volts.shape[0]
                     voltsCohDet = volts.copy()
                     pairsarray = numpy.array(pairslist)
                     indSides = pairsarray[:,1]
             #         indSides = numpy.array(range(nChannel))
             #         indSides = numpy.delete(indSides, indCenter)
             #
             #         listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0)
                     listBlocks = numpy.array_split(volts, numBlocks, 1)
                     startInd = 0
                     endInd = 0
                     for i in range(numBlocks):
                         startInd = endInd
                         endInd = endInd + listBlocks[i].shape[1]
                         arrayBlock = listBlocks[i]
             #             arrayBlockCenter = listCenter[i]
                         #Estimate the Phase Difference
                         phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist)
                         #Phase Difference RMS
                         arrayPhaseRMS = numpy.abs(phaseDiff)
                         phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0)
                         indPhase = numpy.where(phaseRMSaux==4)
                         #Shifting
                         if indPhase[0].shape[0] > 0:
                             for j in range(indSides.size):
                                 arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose())
                             voltsCohDet[:,startInd:endInd,:] = arrayBlock
                     return voltsCohDet
                 def __calculateCCF(self, volts, pairslist ,laglist):
                     nHeights = volts.shape[2]
                     nPoints = volts.shape[1]
                     voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex')
                     for i in range(len(pairslist)):
                         volts1 = volts[pairslist[i][0]]
                         volts2 = volts[pairslist[i][1]]
                         for t in range(len(laglist)):
                             idxT = laglist[t]
                             if idxT >= 0:
                                 vStacked = numpy.vstack((volts2[idxT:,:],
                                                        numpy.zeros((idxT, nHeights),dtype='complex')))
                             else:
                                 vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'),
                                                        volts2[:(nPoints + idxT),:]))
                             voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0)
                             vStacked = None
                     return voltsCCF
                 def __getNoise(self, power, timeSegment, timeInterval):
                     numProfPerBlock = numpy.ceil(timeSegment/timeInterval)
                     numBlocks = int(power.shape[0]/numProfPerBlock)
                     numHeights = power.shape[1]
                     listPower = numpy.array_split(power, numBlocks, 0)
                     noise = numpy.zeros((power.shape[0], power.shape[1]))
                     noise1 = numpy.zeros((power.shape[0], power.shape[1]))
                     startInd = 0
                     endInd = 0
                     for i in range(numBlocks):             #split por canal
                         startInd = endInd
                         endInd = endInd + listPower[i].shape[0]
                         arrayBlock = listPower[i]
                         noiseAux = numpy.mean(arrayBlock, 0)
             #             noiseAux = numpy.median(noiseAux)
             #             noiseAux = numpy.mean(arrayBlock)
                         noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux
                         noiseAux1 = numpy.mean(arrayBlock)
                         noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1
                     return noise, noise1
                 def __findMeteors(self, power, thresh):
                     nProf = power.shape[0]
                     nHeights = power.shape[1]
                     listMeteors = []
                     for i in range(nHeights):
                         powerAux = power[:,i]
                         threshAux = thresh[:,i]
                         indUPthresh = numpy.where(powerAux > threshAux)[0]
                         indDNthresh = numpy.where(powerAux <= threshAux)[0]
                         j = 0
                         while (j < indUPthresh.size - 2):
                             if (indUPthresh[j + 2] == indUPthresh[j] + 2):
                                 indDNAux = numpy.where(indDNthresh > indUPthresh[j])
                                 indDNthresh = indDNthresh[indDNAux]
                                 if (indDNthresh.size > 0):
                                     indEnd = indDNthresh[0] - 1
                                     indInit = indUPthresh[j]
                                     meteor = powerAux[indInit:indEnd + 1]
                                     indPeak = meteor.argmax() + indInit
                                     FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0)))
                                     listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!!
                                     j = numpy.where(indUPthresh == indEnd)[0] + 1
                                 else: j+=1
                             else: j+=1
                     return listMeteors
                 def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit):
                     arrayMeteors = numpy.asarray(listMeteors)
                     listMeteors1 = []
                     while arrayMeteors.shape[0] > 0:
                         FLAs = arrayMeteors[:,4]
                         maxFLA = FLAs.argmax()
                         listMeteors1.append(arrayMeteors[maxFLA,:])
                         MeteorInitTime = arrayMeteors[maxFLA,1]
                         MeteorEndTime = arrayMeteors[maxFLA,3]
                         MeteorHeight = arrayMeteors[maxFLA,0]
                         #Check neighborhood
                         maxHeightIndex = MeteorHeight + rangeLimit
                         minHeightIndex = MeteorHeight - rangeLimit
                         minTimeIndex = MeteorInitTime - timeLimit
                         maxTimeIndex = MeteorEndTime + timeLimit
                         #Check Heights
                         indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex)
                         indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex)
                         indBoth = numpy.where(numpy.logical_and(indTime,indHeight))
                         arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0)
                     return listMeteors1
                 def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency):
                     numHeights = volts.shape[2]
                     nChannel = volts.shape[0]
                     thresholdPhase = thresh[0]
                     thresholdNoise = thresh[1]
                     thresholdDB = float(thresh[2])
                     thresholdDB1 = 10**(thresholdDB/10)
                     pairsarray = numpy.array(pairslist)
                     indSides = pairsarray[:,1]
                     pairslist1 = list(pairslist)
                     pairslist1.append((0,1))
                     pairslist1.append((3,4))
                     listMeteors1 = []
                     listPowerSeries = []
                     listVoltageSeries = []
                     #volts has the war data
                     if frequency == 30e6:
                         timeLag = 45*10**-3
                     else:
                         timeLag = 15*10**-3
                     lag = numpy.ceil(timeLag/timeInterval)
                     for i in range(len(listMeteors)):
                         ######################   3.6 - 3.7 PARAMETERS REESTIMATION    #########################
                         meteorAux = numpy.zeros(16)
                         #Loading meteor Data (mHeight, mStart, mPeak, mEnd)
                         mHeight = listMeteors[i][0]
                         mStart = listMeteors[i][1]
                         mPeak = listMeteors[i][2]
                         mEnd = listMeteors[i][3]
                         #get the volt data between the start and end times of the meteor
                         meteorVolts = volts[:,mStart:mEnd+1,mHeight]
                         meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
                         #3.6. Phase Difference estimation
                         phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist)
                         #3.7. Phase difference removal & meteor start, peak and end times reestimated
                         #meteorVolts0.- all Channels, all Profiles
                         meteorVolts0 = volts[:,:,mHeight]
                         meteorThresh = noise[:,mHeight]*thresholdNoise
                         meteorNoise = noise[:,mHeight]
                         meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting
                         powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0)  #Power
                         #Times reestimation
                         mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0]
                         if mStart1.size > 0:
                             mStart1 = mStart1[-1] + 1
                         else:
                             mStart1 = mPeak
                         mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1
                         mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0]
                         if mEndDecayTime1.size == 0:
                             mEndDecayTime1 = powerNet0.size
                         else:
                             mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1
             #             mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax()
                         #meteorVolts1.- all Channels, from start to end
                         meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1]
                         meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1]
                         if meteorVolts2.shape[1] == 0:
                             meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1]
                         meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1)
                         meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1)
                         #####################    END PARAMETERS REESTIMATION    #########################
                         #####################   3.8 PHASE DIFFERENCE REESTIMATION  ########################
             #             if mEnd1 - mStart1 > 4:       #Error Number 6: echo less than 5 samples long; too short for analysis
                         if meteorVolts2.shape[1] > 0:
                             #Phase Difference re-estimation
                             phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1)   #Phase Difference Estimation
             #                 phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist)
                             meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1])
                             phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1))
                             meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4])     #Phase Shifting
                             #Phase Difference RMS
                             phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1)))
                             powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0)
                             #Data from Meteor
                             mPeak1 = powerNet1.argmax() + mStart1
                             mPeakPower1 = powerNet1.max()
                             noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight])
                             mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux
                             Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1])
                             Meteor1 = numpy.hstack((Meteor1,phaseDiffint))
                             PowerSeries  = powerNet0[mStart1:mEndDecayTime1 + 1]
                             #Vectorize
                             meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]
                             meteorAux[7:11] = phaseDiffint[0:4]
                             #Rejection Criterions
                             if phaseRMS1 > thresholdPhase:  #Error Number 17: Phase variation
                                 meteorAux[-1] = 17
                             elif mSNR1 < thresholdDB1:      #Error Number 1: SNR < threshold dB
                                 meteorAux[-1] = 1
                         else:
                             meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd]
                             meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis
                             PowerSeries = 0
                         listMeteors1.append(meteorAux)
                         listPowerSeries.append(PowerSeries)
                         listVoltageSeries.append(meteorVolts1)
                     return listMeteors1, listPowerSeries, listVoltageSeries
                 def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency):
                     threshError = 10
                     #Depending if it is 30 or 50 MHz
                     if frequency == 30e6:
                         timeLag = 45*10**-3
                     else:
                         timeLag = 15*10**-3
                     lag = numpy.ceil(timeLag/timeInterval)
                     listMeteors1 = []
                     for i in range(len(listMeteors)):
                         meteorPower = listPower[i]
                         meteorAux = listMeteors[i]
                         if meteorAux[-1] == 0:
                             try:
                                 indmax = meteorPower.argmax()
                                 indlag = indmax + lag
                                 y = meteorPower[indlag:]
                                 x = numpy.arange(0, y.size)*timeLag
                                 #first guess
                                 a = y[0]
                                 tau = timeLag
                                 #exponential fit
                                 popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau])
                                 y1 = self.__exponential_function(x, *popt)
                                 #error estimation
                                 error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size))
                                 decayTime = popt[1]
                                 riseTime = indmax*timeInterval
                                 meteorAux[11:13] = [decayTime, error]
                                 #Table items 7, 8 and 11
                                 if (riseTime > 0.3):            #Number 7: Echo rise exceeds 0.3s
                                     meteorAux[-1] = 7
                                 elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time
                                     meteorAux[-1] = 8
                                 if (error > threshError):       #Number 11: Poor fit to amplitude for estimation of decay time
                                     meteorAux[-1] = 11
                             except:
                                 meteorAux[-1] = 11
                         listMeteors1.append(meteorAux)
                     return listMeteors1
                 #Exponential Function
                 def __exponential_function(self, x, a, tau):
                     y = a*numpy.exp(-x/tau)
                     return y
                 def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist,  timeInterval):
                     pairslist1 = list(pairslist)
                     pairslist1.append((0,1))
                     pairslist1.append((3,4))
                     numPairs = len(pairslist1)
                     #Time Lag
                     timeLag = 45*10**-3
                     c = 3e8
                     lag = numpy.ceil(timeLag/timeInterval)
                     freq = 30e6
                     listMeteors1 = []
                     for i in range(len(listMeteors)):
                         meteorAux = listMeteors[i]
                         if meteorAux[-1] == 0:
                             mStart = listMeteors[i][1]
                             mPeak = listMeteors[i][2]
                             mLag = mPeak - mStart + lag
                             #get the volt data between the start and end times of the meteor
                             meteorVolts = listVolts[i]
                             meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1)
                             #Get CCF
                             allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2])
                             #Method 2
                             slopes = numpy.zeros(numPairs)
                             time = numpy.array([-2,-1,1,2])*timeInterval
                             angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0])
                             #Correct phases
                             derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1]
                             indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi)
                             if indDer[0].shape[0] > 0:
                                 for i in range(indDer[0].shape[0]):
                                     signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]])
                                     angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi
             #                     fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]]))
                             for j in range(numPairs):
                                 fit = stats.linregress(time, angAllCCF[j,:])
                                 slopes[j] = fit[0]
                             #Remove Outlier
             #                 indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
             #                 slopes = numpy.delete(slopes,indOut)
             #                 indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes)))
             #                 slopes = numpy.delete(slopes,indOut)
                             radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq)
                             radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq)
                             meteorAux[-2] = radialError
                             meteorAux[-3] = radialVelocity
                             #Setting Error
                             #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s
                             if numpy.abs(radialVelocity) > 200:
                                 meteorAux[-1] = 15
                             #Number 12: Poor fit to CCF variation for estimation of radial drift velocity
                             elif radialError > radialStdThresh:
                                 meteorAux[-1] = 12
                         listMeteors1.append(meteorAux)
                     return listMeteors1
                 def __setNewArrays(self, listMeteors, date, heiRang):
                     #New arrays
                     arrayMeteors = numpy.array(listMeteors)
                     arrayParameters = numpy.zeros((len(listMeteors), 14))
                     #Date inclusion
                     date = re.findall(r'\((.*?)\)', date)
                     date = date[0].split(',')
                     date = map(int, date)
                     date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]]
                     arrayDate = numpy.tile(date, (len(listMeteors), 1))
                     #Meteor array
             #         arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)]
             #         arrayMeteors = numpy.hstack((arrayDate, arrayMeteors))
                     #Parameters Array
                     arrayParameters[:,:2] = arrayDate #Date
                     arrayParameters[:,2] = heiRang[arrayMeteors[:,0].astype(int)] #Range
                     arrayParameters[:,7:9] = arrayMeteors[:,-3:-1] #Radial velocity and its error
                     arrayParameters[:,9:13] = arrayMeteors[:,7:11]  #Phases
                     arrayParameters[:,-1] = arrayMeteors[:,-1]  #Error
                     return arrayParameters
                 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
                     arrayAOA = numpy.zeros((phases.shape[0],3))
                     cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
                     arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
                     cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
                     arrayAOA[:,2] = cosDirError
                     azimuthAngle = arrayAOA[:,0]
                     zenithAngle = arrayAOA[:,1]
                     #Setting Error
                     #Number 3: AOA not fesible
                     indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
                     error[indInvalid] = 3
                     #Number 4: Large difference in AOAs obtained from different antenna baselines
                     indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
                     error[indInvalid] = 4
                     return arrayAOA, error
                 def __getDirectionCosines(self, arrayPhase, pairsList):
                     #Initializing some variables
                     ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
                     ang_aux = ang_aux.reshape(1,ang_aux.size)
                     cosdir = numpy.zeros((arrayPhase.shape[0],2))
                     cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
                     for i in range(2):
                         #First Estimation
                         phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
                         #Dealias
                         indcsi = numpy.where(phi0_aux > numpy.pi)
                         phi0_aux[indcsi] -= 2*numpy.pi
                         indcsi = numpy.where(phi0_aux < -numpy.pi)
                         phi0_aux[indcsi] += 2*numpy.pi
                         #Direction Cosine 0
                         cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
                         #Most-Accurate Second Estimation
                         phi1_aux =  arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
                         phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
                         #Direction Cosine 1
                         cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
                         #Searching the correct Direction Cosine
                         cosdir0_aux = cosdir0[:,i]
                         cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
                         #Minimum Distance
                         cosDiff = (cosdir1 - cosdir0_aux)**2
                         indcos = cosDiff.argmin(axis = 1)
                         #Saving Value obtained
                         cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
                     return cosdir0, cosdir
                 def __calculateAOA(self, cosdir, azimuth):
                     cosdirX = cosdir[:,0]
                     cosdirY = cosdir[:,1]
                     zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
                     azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
                     angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
                     return angles
                 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
                     Ramb = 375  #Ramb = c/(2*PRF)
                     Re = 6371   #Earth Radius
                     heights = numpy.zeros(Ranges.shape)
                     R_aux = numpy.array([0,1,2])*Ramb
                     R_aux = R_aux.reshape(1,R_aux.size)
                     Ranges = Ranges.reshape(Ranges.size,1)
                     Ri = Ranges + R_aux
                     hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
                     #Check if there is a height between 70 and 110 km
                     h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
                     ind_h = numpy.where(h_bool == 1)[0]
                     hCorr = hi[ind_h, :]
                     ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
                     hCorr = hi[ind_hCorr]
                     heights[ind_h] = hCorr
                     #Setting Error
                     #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
                     #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
                     indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
                     error[indInvalid2] = 14
                     indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
                     error[indInvalid1] = 13
                     return heights, error
                 def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None):
                     '''
                     Function GetMoments()
                     Input:
                     Output:
                     Variables modified:
                     '''
                     if path != None:
                         sys.path.append(path)
                     self.dataOut.library = importlib.import_module(file)
                     #To be inserted as a parameter
                     groupArray = numpy.array(groupList)
             #         groupArray = numpy.array([[0,1],[2,3]])
                     self.dataOut.groupList = groupArray
                     nGroups = groupArray.shape[0]
                     nChannels = self.dataIn.nChannels
                     nHeights=self.dataIn.heightList.size
                     #Parameters Array
                     self.dataOut.data_param = None
                     #Set constants
                     constants = self.dataOut.library.setConstants(self.dataIn)
                     self.dataOut.constants = constants
                     M = self.dataIn.normFactor
                     N = self.dataIn.nFFTPoints
                     ippSeconds = self.dataIn.ippSeconds
                     K = self.dataIn.nIncohInt
                     pairsArray = numpy.array(self.dataIn.pairsList)
                     #List of possible combinations
                     listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2)
                     indCross = numpy.zeros(len(list(listComb)), dtype = 'int')
                     if getSNR:
                         listChannels = groupArray.reshape((groupArray.size))
                         listChannels.sort()
                         noise = self.dataIn.getNoise()
                         self.dataOut.data_SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels])
                     for i in range(nGroups):
                         coord = groupArray[i,:]
                         #Input data array
                         data = self.dataIn.data_spc[coord,:,:]/(M*N)
                         data = data.reshape((data.shape[0]*data.shape[1],data.shape[2]))
                         #Cross Spectra data array for Covariance Matrixes
                         ind = 0
                         for pairs in listComb:
                             pairsSel = numpy.array([coord[x],coord[y]])
                             indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0])
                             ind += 1
                         dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N)
                         dataCross = dataCross**2/K
                         for h in range(nHeights):
             #                 print self.dataOut.heightList[h]
                             #Input
                             d = data[:,h]
                             #Covariance Matrix
                             D = numpy.diag(d**2/K)
                             ind = 0
                             for pairs in listComb:
                                 #Coordinates in Covariance Matrix
                                 x = pairs[0]
                                 y = pairs[1]
                                 #Channel Index
                                 S12 = dataCross[ind,:,h]
                                 D12 = numpy.diag(S12)
                                 #Completing Covariance Matrix with Cross Spectras
                                 D[x*N:(x+1)*N,y*N:(y+1)*N] = D12
                                 D[y*N:(y+1)*N,x*N:(x+1)*N] = D12
                                 ind += 1
                             Dinv=numpy.linalg.inv(D)
                             L=numpy.linalg.cholesky(Dinv)
                             LT=L.T
                             dp = numpy.dot(LT,d)
                             #Initial values
                             data_spc = self.dataIn.data_spc[coord,:,h]
                             if (h>0)and(error1[3]<5):
                                 p0 = self.dataOut.data_param[i,:,h-1]
                             else:
                                 p0 = numpy.array(self.dataOut.library.initialValuesFunction(data_spc, constants, i))
                             try:
                                 #Least Squares
                                 minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True)
             #                   minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants))
                                 #Chi square error
                                 error0 = numpy.sum(infodict['fvec']**2)/(2*N)
                                 #Error with Jacobian
                                 error1 = self.dataOut.library.errorFunction(minp,constants,LT)
                             except:
                                 minp = p0*numpy.nan
                                 error0 = numpy.nan
                                 error1 = p0*numpy.nan
                             #Save
-                            if self.dataOut.data_param == None:
+                            if self.dataOut.data_param is None:
                                 self.dataOut.data_param = numpy.zeros((nGroups, p0.size, nHeights))*numpy.nan
                                 self.dataOut.data_error = numpy.zeros((nGroups, p0.size + 1, nHeights))*numpy.nan
                             self.dataOut.data_error[i,:,h] = numpy.hstack((error0,error1))
                             self.dataOut.data_param[i,:,h] = minp
                     return
                 def __residFunction(self, p, dp, LT, constants):
                     fm = self.dataOut.library.modelFunction(p, constants)
                     fmp=numpy.dot(LT,fm)
                     return  dp-fmp
                 def __getSNR(self, z, noise):
                     avg = numpy.average(z, axis=1)
                     SNR = (avg.T-noise)/noise
                     SNR = SNR.T
                     return SNR
                 def __chisq(p,chindex,hindex):
                     #similar to Resid but calculates CHI**2
                     [LT,d,fm]=setupLTdfm(p,chindex,hindex)
                     dp=numpy.dot(LT,d)
                     fmp=numpy.dot(LT,fm)
                     chisq=numpy.dot((dp-fmp).T,(dp-fmp))
                     return chisq
             class WindProfiler(Operation):
                 __isConfig = False
                 __initime = None
                 __lastdatatime = None
                 __integrationtime = None
                 __buffer = None
                 __dataReady = False
                 __firstdata = None
                 n = None
                 def __init__(self):
                     Operation.__init__(self)
                 def __calculateCosDir(self, elev, azim):
                     zen = (90 - elev)*numpy.pi/180
                     azim = azim*numpy.pi/180
                     cosDirX = numpy.sqrt((1-numpy.cos(zen)**2)/((1+numpy.tan(azim)**2)))
                     cosDirY = numpy.sqrt(1-numpy.cos(zen)**2-cosDirX**2)
                     signX = numpy.sign(numpy.cos(azim))
                     signY = numpy.sign(numpy.sin(azim))
                     cosDirX = numpy.copysign(cosDirX, signX)
                     cosDirY = numpy.copysign(cosDirY, signY)
                     return cosDirX, cosDirY
                 def __calculateAngles(self, theta_x, theta_y, azimuth):
                     dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2)
                     zenith_arr = numpy.arccos(dir_cosw)
                     azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180
                     dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr)
                     dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr)
                     return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw
                 def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly):
             #
                     if horOnly:
                         A = numpy.c_[dir_cosu,dir_cosv]
                     else:
                         A = numpy.c_[dir_cosu,dir_cosv,dir_cosw]
                     A = numpy.asmatrix(A)
                     A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose()
                     return A1
                 def __correctValues(self, heiRang, phi, velRadial, SNR):
                     listPhi = phi.tolist()
                     maxid = listPhi.index(max(listPhi))
                     minid = listPhi.index(min(listPhi))
                     rango = range(len(phi))
                #     rango = numpy.delete(rango,maxid)
                     heiRang1 = heiRang*math.cos(phi[maxid])
                     heiRangAux = heiRang*math.cos(phi[minid])
                     indOut = (heiRang1 < heiRangAux[0]).nonzero()
                     heiRang1 = numpy.delete(heiRang1,indOut)
                     velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
                     SNR1 = numpy.zeros([len(phi),len(heiRang1)])
                     for i in rango:
                         x = heiRang*math.cos(phi[i])
                         y1 = velRadial[i,:]
                         f1 = interpolate.interp1d(x,y1,kind = 'cubic')
                         x1 = heiRang1
                         y11 = f1(x1)
                         y2 = SNR[i,:]
                         f2 = interpolate.interp1d(x,y2,kind = 'cubic')
                         y21 = f2(x1)
                         velRadial1[i,:] = y11
                         SNR1[i,:] = y21
                     return heiRang1, velRadial1, SNR1
                 def __calculateVelUVW(self, A, velRadial):
                     #Operacion Matricial
             #         velUVW = numpy.zeros((velRadial.shape[1],3))
             #         for ind in range(velRadial.shape[1]):
             #             velUVW[ind,:] = numpy.dot(A,velRadial[:,ind])
             #         velUVW = velUVW.transpose()
                     velUVW = numpy.zeros((A.shape[0],velRadial.shape[1]))
                     velUVW[:,:] = numpy.dot(A,velRadial)
                     return velUVW
                 def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0):
                     """
                     Function that implements Doppler Beam Swinging (DBS) technique.
                     Input:    Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
                                 Direction correction (if necessary), Ranges and SNR
                     Output:    Winds estimation (Zonal, Meridional and Vertical)
                     Parameters affected:    Winds, height range, SNR
                     """
                     azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth)
                     heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0)
                     A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly)
                     #Calculo de Componentes de la velocidad con DBS
                     winds = self.__calculateVelUVW(A,velRadial1)
                     return winds, heiRang1, SNR1
                 def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None):
                     posx = numpy.asarray(posx)
                     posy = numpy.asarray(posy)
                     #Rotacion Inversa para alinear con el azimuth
                     if azimuth!= None:
                         azimuth = azimuth*math.pi/180
                         posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth)
                         posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth)
                     else:
                         posx1 = posx
                         posy1 = posy
                     #Calculo de Distancias
                     distx = numpy.zeros(pairsCrossCorr.size)
                     disty = numpy.zeros(pairsCrossCorr.size)
                     dist = numpy.zeros(pairsCrossCorr.size)
                     ang = numpy.zeros(pairsCrossCorr.size)
                     for i in range(pairsCrossCorr.size):
                         distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]]
                         disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]]
                         dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2)
                         ang[i] = numpy.arctan2(disty[i],distx[i])
                     #Calculo de Matrices
                     nPairs = len(pairs)
                     ang1 = numpy.zeros((nPairs, 2, 1))
                     dist1 = numpy.zeros((nPairs, 2, 1))
                     for j in range(nPairs):
                         dist1[j,0,0] = dist[pairs[j][0]]
                         dist1[j,1,0] = dist[pairs[j][1]]
                         ang1[j,0,0] = ang[pairs[j][0]]
                         ang1[j,1,0] = ang[pairs[j][1]]
                     return distx,disty, dist1,ang1
                 def __calculateVelVer(self, phase, lagTRange, _lambda):
                     Ts = lagTRange[1] - lagTRange[0]
                     velW = -_lambda*phase/(4*math.pi*Ts)
                     return velW
                 def __calculateVelHorDir(self, dist, tau1, tau2, ang):
                     nPairs = tau1.shape[0]
                     vel = numpy.zeros((nPairs,3,tau1.shape[2]))
                     angCos = numpy.cos(ang)
                     angSin = numpy.sin(ang)
                     vel0 = dist*tau1/(2*tau2**2)
                     vel[:,0,:] = (vel0*angCos).sum(axis = 1)
                     vel[:,1,:] = (vel0*angSin).sum(axis = 1)
                     ind = numpy.where(numpy.isinf(vel))
                     vel[ind] = numpy.nan
                     return vel
                 def __getPairsAutoCorr(self, pairsList, nChannels):
                     pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan
                     for l in range(len(pairsList)):
                         firstChannel = pairsList[l][0]
                         secondChannel = pairsList[l][1]
                         #Obteniendo pares de Autocorrelacion
                         if firstChannel == secondChannel:
                             pairsAutoCorr[firstChannel] = int(l)
                     pairsAutoCorr = pairsAutoCorr.astype(int)
                     pairsCrossCorr = range(len(pairsList))
                     pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr)
                     return pairsAutoCorr, pairsCrossCorr
                 def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor):
                     """
                     Function that implements Spaced Antenna (SA) technique.
                     Input:    Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth,
                                 Direction correction (if necessary), Ranges and SNR
                     Output:    Winds estimation (Zonal, Meridional and Vertical)
                     Parameters affected:    Winds
                     """
                     #Cross Correlation pairs obtained
                     pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels)
                     pairsArray = numpy.array(pairsList)[pairsCrossCorr]
                     pairsSelArray = numpy.array(pairsSelected)
                     pairs = []
                     #Wind estimation pairs obtained
                     for i in range(pairsSelArray.shape[0]/2):
                         ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0]
                         ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0]
                         pairs.append((ind1,ind2))
                     indtau = tau.shape[0]/2
                     tau1 = tau[:indtau,:]
                     tau2 = tau[indtau:-1,:]
                     tau1 = tau1[pairs,:]
                     tau2 = tau2[pairs,:]
                     phase1 = tau[-1,:]
                     #---------------------------------------------------------------------
                     #Metodo Directo
                     distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth)
                     winds = self.__calculateVelHorDir(dist, tau1, tau2, ang)
                     winds = stats.nanmean(winds, axis=0)
                     #---------------------------------------------------------------------
                     #Metodo General
             #         distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth)
             #         #Calculo Coeficientes de Funcion de Correlacion
             #         F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n)
             #         #Calculo de Velocidades
             #         winds = self.calculateVelUV(F,G,A,B,H)
                     #---------------------------------------------------------------------
                     winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda)
                     winds = correctFactor*winds
                     return winds
                 def __checkTime(self, currentTime, paramInterval, outputInterval):
                     dataTime = currentTime + paramInterval
                     deltaTime = dataTime - self.__initime
                     if deltaTime >= outputInterval or deltaTime < 0:
                         self.__dataReady = True
                     return
                 def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax):
                     '''
                     Function that implements winds estimation technique with detected meteors.
                     Input:    Detected meteors, Minimum meteor quantity to wind estimation
                     Output:    Winds estimation (Zonal and Meridional)
                     Parameters affected:    Winds
                     '''
                     print arrayMeteor.shape
                     #Settings
                     nInt = (heightMax - heightMin)/2
                     winds = numpy.zeros((2,nInt))*numpy.nan
                     #Filter errors
                     error = numpy.where(arrayMeteor[0,:,-1] == 0)[0]
                     finalMeteor = arrayMeteor[0,error,:]
                     #Meteor Histogram
                     finalHeights = finalMeteor[:,3]
                     hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax))
                     nMeteorsPerI = hist[0]
                     heightPerI = hist[1]
                     #Sort of meteors
                     indSort = finalHeights.argsort()
                     finalMeteor2 = finalMeteor[indSort,:]
                     #    Calculating winds
                     ind1 = 0
                     ind2 = 0
                     for i in range(nInt):
                         nMet = nMeteorsPerI[i]
                         ind1 = ind2
                         ind2 = ind1 + nMet
                         meteorAux = finalMeteor2[ind1:ind2,:]
                         if meteorAux.shape[0] >= meteorThresh:
                             vel = meteorAux[:, 7]
                             zen = meteorAux[:, 5]*numpy.pi/180
                             azim = meteorAux[:, 4]*numpy.pi/180
                             n = numpy.cos(zen)
                     #         m = (1 - n**2)/(1 - numpy.tan(azim)**2)
                     #         l = m*numpy.tan(azim)
                             l = numpy.sin(zen)*numpy.sin(azim)
                             m = numpy.sin(zen)*numpy.cos(azim)
                             A = numpy.vstack((l, m)).transpose()
                             A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose())
                             windsAux = numpy.dot(A1, vel)
                             winds[0,i] = windsAux[0]
                             winds[1,i] = windsAux[1]
                     return winds, heightPerI[:-1]
                 def run(self, dataOut, technique, **kwargs):
                     param = dataOut.data_param
             #         if dataOut.abscissaList != None:
             #             absc = dataOut.abscissaList[:-1]
                     noise = dataOut.noise
                     heightList = dataOut.heightList
                     SNR = dataOut.data_SNR
                     if technique == 'DBS':
                         if kwargs.has_key('dirCosx') and kwargs.has_key('dirCosy'):
                             theta_x = numpy.array(kwargs['dirCosx'])
                             theta_y = numpy.array(kwargs['dirCosy'])
                         else:
                             elev = numpy.array(kwargs['elevation'])
                             azim = numpy.array(kwargs['azimuth'])
                             theta_x, theta_y = self.__calculateCosDir(elev, azim)
                         azimuth = kwargs['correctAzimuth']
                         if kwargs.has_key('horizontalOnly'):
                             horizontalOnly = kwargs['horizontalOnly']
                         else:   horizontalOnly = False
                         if kwargs.has_key('correctFactor'):
                             correctFactor = kwargs['correctFactor']
                         else:   correctFactor = 1
                         if kwargs.has_key('channelList'):
                             channelList = kwargs['channelList']
                             if len(channelList) == 2:
                                 horizontalOnly = True
                             arrayChannel = numpy.array(channelList)
                             param = param[arrayChannel,:,:]
                             theta_x = theta_x[arrayChannel]
                             theta_y = theta_y[arrayChannel]
                         velRadial0 = param[:,1,:] #Radial velocity
                         dataOut.data_output, dataOut.heightList, dataOut.data_SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightList, SNR) #DBS Function
                         dataOut.utctimeInit = dataOut.utctime
                         dataOut.outputInterval = dataOut.timeInterval
                     elif technique == 'SA':
                         #Parameters
                         position_x = kwargs['positionX']
                         position_y = kwargs['positionY']
                         azimuth = kwargs['azimuth']
                         if kwargs.has_key('crosspairsList'):
                             pairs = kwargs['crosspairsList']
                         else:
                             pairs = None
                         if kwargs.has_key('correctFactor'):
                             correctFactor = kwargs['correctFactor']
                         else:
                             correctFactor = 1
                         tau = dataOut.data_param
                         _lambda = dataOut.C/dataOut.frequency
                         pairsList = dataOut.groupList
                         nChannels = dataOut.nChannels
                         dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor)
                         dataOut.utctimeInit = dataOut.utctime
                         dataOut.outputInterval = dataOut.timeInterval
                     elif technique == 'Meteors':
                         dataOut.flagNoData = True
                         self.__dataReady = False
                         if kwargs.has_key('nHours'):
                             nHours = kwargs['nHours']
                         else:
                             nHours = 1
                         if kwargs.has_key('meteorsPerBin'):
                             meteorThresh = kwargs['meteorsPerBin']
                         else:
                             meteorThresh = 6
                         if kwargs.has_key('hmin'):
                             hmin = kwargs['hmin']
                         else:   hmin = 70
                         if kwargs.has_key('hmax'):
                             hmax = kwargs['hmax']
                         else:   hmax = 110
                         dataOut.outputInterval = nHours*3600
                         if self.__isConfig == False:
             #                 self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
                             #Get Initial LTC time
                             self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
                             self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
                             self.__isConfig = True
-                        if self.__buffer == None:
+                        if self.__buffer is None:
                             self.__buffer = dataOut.data_param
                             self.__firstdata = copy.copy(dataOut)
                         else:
                             self.__buffer = numpy.hstack((self.__buffer, dataOut.data_param))
                         self.__checkTime(dataOut.utctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
                         if self.__dataReady:
                             dataOut.utctimeInit = self.__initime
                             self.__initime += dataOut.outputInterval #to erase time offset
                             dataOut.data_output, dataOut.heightList = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax)
                             dataOut.flagNoData = False
                             self.__buffer = None
                     return
             class EWDriftsEstimation(Operation):
                 def __init__(self):
                     Operation.__init__(self)
                 def __correctValues(self, heiRang, phi, velRadial, SNR):
                     listPhi = phi.tolist()
                     maxid = listPhi.index(max(listPhi))
                     minid = listPhi.index(min(listPhi))
                     rango = range(len(phi))
                #     rango = numpy.delete(rango,maxid)
                     heiRang1 = heiRang*math.cos(phi[maxid])
                     heiRangAux = heiRang*math.cos(phi[minid])
                     indOut = (heiRang1 < heiRangAux[0]).nonzero()
                     heiRang1 = numpy.delete(heiRang1,indOut)
                     velRadial1 = numpy.zeros([len(phi),len(heiRang1)])
                     SNR1 = numpy.zeros([len(phi),len(heiRang1)])
                     for i in rango:
                         x = heiRang*math.cos(phi[i])
                         y1 = velRadial[i,:]
                         f1 = interpolate.interp1d(x,y1,kind = 'cubic')
                         x1 = heiRang1
                         y11 = f1(x1)
                         y2 = SNR[i,:]
                         f2 = interpolate.interp1d(x,y2,kind = 'cubic')
                         y21 = f2(x1)
                         velRadial1[i,:] = y11
                         SNR1[i,:] = y21
                     return heiRang1, velRadial1, SNR1
                 def run(self, dataOut, zenith, zenithCorrection):
                     heiRang = dataOut.heightList
                     velRadial = dataOut.data_param[:,3,:]
                     SNR = dataOut.data_SNR
                     zenith = numpy.array(zenith)
                     zenith -= zenithCorrection
                     zenith *= numpy.pi/180
                     heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR)
                     alp = zenith[0]
                     bet = zenith[1]
                     w_w = velRadial1[0,:]
                     w_e = velRadial1[1,:]
                     w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp))
                     u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp))
                     winds = numpy.vstack((u,w))
                     dataOut.heightList = heiRang1
                     dataOut.data_output = winds
                     dataOut.data_SNR = SNR1
                     dataOut.utctimeInit = dataOut.utctime
                     dataOut.outputInterval = dataOut.timeInterval
                     return
             class PhaseCalibration(Operation):
                 __buffer = None
                 __initime = None
                 __dataReady = False
                 __isConfig = False
                 def __checkTime(self, currentTime, initTime, paramInterval, outputInterval):
                     dataTime = currentTime + paramInterval
                     deltaTime = dataTime - initTime
                     if deltaTime >= outputInterval or deltaTime < 0:
                         return True
                     return False
                 def __getGammas(self, pairs, k, d, phases):
                     gammas = numpy.zeros(2)
                     for i in range(len(pairs)):
                         pairi = pairs[i]
                         #Calculating gamma
                         jdcos = phases[:,pairi[1]]/(k*d[pairi[1]])
                         jgamma = numpy.angle(numpy.exp(1j*(k*d[pairi[0]]*jdcos - phases[:,pairi[0]])))
                         #Revised distribution
                         jgammaArray = numpy.hstack((jgamma,jgamma+0.5*numpy.pi,jgamma-0.5*numpy.pi))
                         #Histogram
                         nBins = 64.0
                         rmin = -0.5*numpy.pi
                         rmax = 0.5*numpy.pi
                         phaseHisto = numpy.histogram(jgammaArray, bins=nBins, range=(rmin,rmax))
                         meteorsY = phaseHisto[0]
                         phasesX = phaseHisto[1][:-1]
                         width = phasesX[1] - phasesX[0]
                         phasesX += width/2
                         #Gaussian aproximation
                         bpeak = meteorsY.argmax()
                         peak = meteorsY.max()
                         jmin = bpeak - 5
                         jmax = bpeak + 5 + 1
                         if jmin<0:
                             jmin = 0
                             jmax = 6
                         elif jmax > meteorsY.size:
                             jmin = meteorsY.size - 6
                             jmax = meteorsY.size
                         x0 = numpy.array([peak,bpeak,50])
                         coeff = optimize.leastsq(self.__residualFunction, x0, args=(meteorsY[jmin:jmax], phasesX[jmin:jmax]))
                         #Gammas
                         gammas[i] = coeff[0][1]
             #             gammas[i] = bpeak
                     return gammas
                 def __residualFunction(self, coeffs, y, t):
                     return y - self.__gauss_function(t, coeffs)
                 def __gauss_function(self, t, coeffs):
                     return coeffs[0]*numpy.exp(-0.5*((t - coeffs[1]) / coeffs[2])**2)
                 def __getPhases(self, azimuth, h, pairsList, d, gammas, meteorsArray):
                     meteorOps = MeteorOperations()
                     nchan = 4
                     pairx = pairsList[0]
                     pairy = pairsList[1]
                     center_xangle = 0
                     center_yangle = 0
                     range_angle = numpy.array([8*numpy.pi,numpy.pi,numpy.pi/2,numpy.pi/4])
                     ntimes = len(range_angle)
                     nstepsx = 20.0
                     nstepsy = 20.0
                     for iz in range(ntimes):
                         min_xangle = -range_angle[iz]/2 + center_xangle
                         max_xangle = range_angle[iz]/2 + center_xangle
                         min_yangle = -range_angle[iz]/2 + center_yangle
                         max_yangle = range_angle[iz]/2 + center_yangle
                         inc_x = (max_xangle-min_xangle)/nstepsx
                         inc_y = (max_yangle-min_yangle)/nstepsy
                         alpha_y = numpy.arange(nstepsy)*inc_y + min_yangle
                         alpha_x = numpy.arange(nstepsx)*inc_x + min_xangle
                         penalty = numpy.zeros((nstepsx,nstepsy))
                         jph_array = numpy.zeros((nchan,nstepsx,nstepsy))
                         jph = numpy.zeros(nchan)
                         # Iterations looking for the offset
                         for iy in range(int(nstepsy)):
                             for ix in range(int(nstepsx)):
                                 jph[pairy[1]] = alpha_y[iy]
                                 jph[pairy[0]] = -gammas[1] + alpha_y[iy]*d[pairy[0]]/d[pairy[1]]
                                 jph[pairx[1]] = alpha_x[ix]
                                 jph[pairx[0]] = -gammas[0] + alpha_x[ix]*d[pairx[0]]/d[pairx[1]]
                                 jph_array[:,ix,iy] = jph
                                 meteorsArray1 = meteorOps.getMeteorParams(meteorsArray, azimuth, h, pairsList, jph)
                                 error = meteorsArray1[:,-1]
                                 ind1 = numpy.where(error==0)[0]
                                 penalty[ix,iy] = ind1.size
                         i,j = numpy.unravel_index(penalty.argmax(), penalty.shape)
                         phOffset = jph_array[:,i,j]
                         center_xangle = phOffset[pairx[1]]
                         center_yangle = phOffset[pairy[1]]
                     phOffset = numpy.angle(numpy.exp(1j*jph_array[:,i,j]))
                     phOffset = phOffset*180/numpy.pi
                     return phOffset
                 def run(self, dataOut, pairs, distances, hmin, hmax, nHours = None):
                     dataOut.flagNoData = True
                     self.__dataReady = False
                     if nHours == None:
                         nHours = 1
                     dataOut.outputInterval = nHours*3600
                     if self.__isConfig == False:
             #                 self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03)
                         #Get Initial LTC time
                         self.__initime = datetime.datetime.utcfromtimestamp(dataOut.utctime)
                         self.__initime = (self.__initime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds()
                         self.__isConfig = True
-                    if self.__buffer == None:
+                    if self.__buffer is None:
                         self.__buffer = dataOut.data_param.copy()
                     else:
                         self.__buffer = numpy.hstack((self.__buffer, dataOut.data_param))
                     self.__dataReady = self.__checkTime(dataOut.utctime, self.__initime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready
                     if self.__dataReady:
                         dataOut.utctimeInit = self.__initime
                         self.__initime += dataOut.outputInterval #to erase time offset
                         freq = dataOut.frequency
                         c = dataOut.C #m/s
                         lamb = c/freq
                         k = 2*numpy.pi/lamb
                         azimuth = 0
                         h = (hmin, hmax)
                         pairsList = ((0,3),(1,2))
                         meteorsArray = self.__buffer[0,:,:]
                         error = meteorsArray[:,-1]
                         boolError = (error==0)|(error==3)|(error==4)|(error==13)|(error==14)
                         ind1 = numpy.where(boolError)[0]
                         meteorsArray = meteorsArray[ind1,:]
                         meteorsArray[:,-1] = 0
                         phases = meteorsArray[:,9:13]
                         #Calculate Gammas
                         gammas = self.__getGammas(pairs, k, distances, phases)
             #             gammas = numpy.array([-21.70409463,45.76935864])*numpy.pi/180
                         #Calculate Phases
                         phasesOff = self.__getPhases(azimuth, h, pairsList, distances, gammas, meteorsArray)
                         phasesOff = phasesOff.reshape((1,phasesOff.size))
                         dataOut.data_output = -phasesOff
                         dataOut.flagNoData = False
                         self.__buffer = None
                     return
             class MeteorOperations():
                 def __init__(self):
                     return
                 def getMeteorParams(self, arrayParameters0, azimuth, h, pairsList, jph):
                     arrayParameters = arrayParameters0.copy()
                     hmin = h[0]
                     hmax = h[1]
                     #Calculate AOA (Error N 3, 4)
                     #JONES ET AL. 1998
                     AOAthresh = numpy.pi/8
                     error = arrayParameters[:,-1]
                     phases = -arrayParameters[:,9:13] + jph
                     arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth)
                     #Calculate Heights (Error N 13 and 14)
                     error = arrayParameters[:,-1]
                     Ranges = arrayParameters[:,2]
                     zenith = arrayParameters[:,5]
                     arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax)
                     #----------------------- Get Final data    ------------------------------------
             #         error = arrayParameters[:,-1]
             #         ind1 = numpy.where(error==0)[0]
             #         arrayParameters = arrayParameters[ind1,:]
                     return arrayParameters
                 def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth):
                     arrayAOA = numpy.zeros((phases.shape[0],3))
                     cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList)
                     arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth)
                     cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1)
                     arrayAOA[:,2] = cosDirError
                     azimuthAngle = arrayAOA[:,0]
                     zenithAngle = arrayAOA[:,1]
                     #Setting Error
                     #Number 3: AOA not fesible
                     indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0]
                     error[indInvalid] = 3
                     #Number 4: Large difference in AOAs obtained from different antenna baselines
                     indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0]
                     error[indInvalid] = 4
                     return arrayAOA, error
                 def __getDirectionCosines(self, arrayPhase, pairsList):
                     #Initializing some variables
                     ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi
                     ang_aux = ang_aux.reshape(1,ang_aux.size)
                     cosdir = numpy.zeros((arrayPhase.shape[0],2))
                     cosdir0 = numpy.zeros((arrayPhase.shape[0],2))
                     for i in range(2):
                         #First Estimation
                         phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]]
                         #Dealias
                         indcsi = numpy.where(phi0_aux > numpy.pi)
                         phi0_aux[indcsi] -= 2*numpy.pi
                         indcsi = numpy.where(phi0_aux < -numpy.pi)
                         phi0_aux[indcsi] += 2*numpy.pi
                         #Direction Cosine 0
                         cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5)
                         #Most-Accurate Second Estimation
                         phi1_aux =  arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]]
                         phi1_aux = phi1_aux.reshape(phi1_aux.size,1)
                         #Direction Cosine 1
                         cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5)
                         #Searching the correct Direction Cosine
                         cosdir0_aux = cosdir0[:,i]
                         cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1)
                         #Minimum Distance
                         cosDiff = (cosdir1 - cosdir0_aux)**2
                         indcos = cosDiff.argmin(axis = 1)
                         #Saving Value obtained
                         cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos]
                     return cosdir0, cosdir
                 def __calculateAOA(self, cosdir, azimuth):
                     cosdirX = cosdir[:,0]
                     cosdirY = cosdir[:,1]
                     zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi
                     azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east
                     angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose()
                     return angles
                 def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight):
                     Ramb = 375  #Ramb = c/(2*PRF)
                     Re = 6371   #Earth Radius
                     heights = numpy.zeros(Ranges.shape)
                     R_aux = numpy.array([0,1,2])*Ramb
                     R_aux = R_aux.reshape(1,R_aux.size)
                     Ranges = Ranges.reshape(Ranges.size,1)
                     Ri = Ranges + R_aux
                     hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re
                     #Check if there is a height between 70 and 110 km
                     h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1)
                     ind_h = numpy.where(h_bool == 1)[0]
                     hCorr = hi[ind_h, :]
                     ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight))
                     hCorr = hi[ind_hCorr]
                     heights[ind_h] = hCorr
                     #Setting Error
                     #Number 13: Height unresolvable echo: not valid height within 70 to 110 km
                     #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km
                     indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0]
                     error[indInvalid2] = 14
                     indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0]
                     error[indInvalid1] = 13
                     return heights, error
   No newline at end of file

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

             import numpy
             import math
             from jroproc_base import ProcessingUnit, Operation
             from schainpy.model.data.jrodata import Spectra
             from schainpy.model.data.jrodata import hildebrand_sekhon
             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
                 def __updateObjFromInput(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.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.nHeights = self.dataIn.nHeights
             #        self.dataOut.nChannels = self.dataIn.nChannels
                     self.dataOut.nBaud = self.dataIn.nBaud
                     self.dataOut.nCode = self.dataIn.nCode
                     self.dataOut.code = self.dataIn.code
                     self.dataOut.nProfiles = self.dataOut.nFFTPoints
             #        self.dataOut.channelIndexList = self.dataIn.channelIndexList
                     self.dataOut.flagDiscontinuousBlock = self.dataIn.flagDiscontinuousBlock
                     self.dataOut.utctime = self.firstdatatime
                     self.dataOut.flagDecodeData = self.dataIn.flagDecodeData #asumo q la data esta decodificada
                     self.dataOut.flagDeflipData = self.dataIn.flagDeflipData #asumo q la data esta sin flip
             #        self.dataOut.flagShiftFFT = self.dataIn.flagShiftFFT
                     self.dataOut.nCohInt = self.dataIn.nCohInt
                     self.dataOut.nIncohInt = 1
             #         self.dataOut.ippSeconds = self.dataIn.ippSeconds
                     self.dataOut.windowOfFilter = self.dataIn.windowOfFilter
             #         self.dataOut.timeInterval = self.dataIn.timeInterval*self.dataOut.nFFTPoints*self.dataOut.nIncohInt
                     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.beam.codeList = self.dataIn.beam.codeList
                     self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
                     self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
                 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=[], ippFactor=None):
                     self.dataOut.flagNoData = True
                     if self.dataIn.type == "Spectra":
                         self.dataOut.copy(self.dataIn)
                         return True
                     if self.dataIn.type == "Voltage":
                         if nFFTPoints == None:
                             raise ValueError, "This SpectraProc.run() need nFFTPoints input variable"
                         if nProfiles == None:
                             nProfiles = nFFTPoints
             #                 raise ValueError, "This SpectraProc.run() need nProfiles input variable"
                         if ippFactor == None:
                             ippFactor = 1
                         self.dataOut.ippFactor = ippFactor
                         self.dataOut.nFFTPoints = nFFTPoints
                         self.dataOut.pairsList = pairsList
-                        if self.buffer == None:
+                        if self.buffer is None:
                             self.buffer = numpy.zeros((self.dataIn.nChannels,
                                                        nProfiles,
                                                        self.dataIn.nHeights),
                                                     dtype='complex')
                             self.id_min = 0
                             self.id_max = self.dataIn.data.shape[1]
                         if len(self.dataIn.data.shape) == 2:
                             self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()
                             self.profIndex += 1
                         else:
                             if self.dataIn.data.shape[1] == nProfiles:
                                 self.buffer = self.dataIn.data.copy()
                                 self.profIndex = nProfiles
                             elif self.dataIn.data.shape[1] < nProfiles:
                                 self.buffer[:,self.id_min:self.id_max,:] = self.dataIn.data
                                 self.profIndex += self.dataIn.data.shape[1]
                                 self.id_min += self.dataIn.data.shape[1]
                                 self.id_max += self.dataIn.data.shape[1]
                             else:
                                 raise ValueError, "The type object %s has %d profiles, it should be equal to %d profiles"%(self.dataIn.type,self.dataIn.data.shape[1],nProfiles)
                                 self.dataOut.flagNoData = True
                                 return 0
                         if self.firstdatatime == None:
                             self.firstdatatime = self.dataIn.utctime
                         if self.profIndex == nProfiles:
                             self.__updateObjFromInput()
                             self.__getFft()
                             self.dataOut.flagNoData = False
                             self.buffer = None
                             self.firstdatatime = None
                             self.profIndex = 0
                         return True
                     raise ValueError, "The type of input object '%s' is not valid"%(self.dataIn.type)
                 def __selectPairs(self, channelList=None):
                     if channelList == None:
                         return
                     pairsIndexListSelected = []
                     for pairIndex in self.dataOut.pairsIndexList:
                         #First pair
                         if self.dataOut.pairsList[pairIndex][0] not in channelList:
                             continue
                         #Second pair
                         if self.dataOut.pairsList[pairIndex][1] not in channelList:
                             continue
                         pairsIndexListSelected.append(pairIndex)
                     if not pairsIndexListSelected:
                         self.dataOut.data_cspc = None
                         self.dataOut.pairsList = []
                         return
                     self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
                     self.dataOut.pairsList = self.dataOut.pairsList[pairsIndexListSelected]
                     return
                 def selectChannels(self, channelList):
                     channelIndexList = []
                     for channel in channelList:
                         if channel not in self.dataOut.channelList:
                             raise ValueError, "Error selecting channels: The value %d in channelList is not valid.\nAvailable channels = %s" %(channel, str(self.dataOut.channelList))
                         index = self.dataOut.channelList.index(channel)
                         channelIndexList.append(index)
                     self.selectChannelsByIndex(channelIndexList)
                 def selectChannelsByIndex(self, 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:
                         self.dataOut.data_spc
                         self.dataOut.channelIndexList
                         self.dataOut.nChannels
                     Return:
                         None
                     """
                     for channelIndex in channelIndexList:
                         if channelIndex not in self.dataOut.channelIndexList:
                             raise ValueError, "Error selecting channels: The value %d in channelIndexList is not valid.\nAvailable channel indexes = " %(channelIndex, self.dataOut.channelIndexList)
             #         nChannels = len(channelIndexList)
                     data_spc = self.dataOut.data_spc[channelIndexList,:]
                     data_dc = self.dataOut.data_dc[channelIndexList,:]
                     self.dataOut.data_spc = data_spc
                     self.dataOut.data_dc = data_dc
                     self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
             #        self.dataOut.nChannels = nChannels
                     self.__selectPairs(self.dataOut.channelList)
                     return 1
                 def selectHeights(self, minHei, maxHei):
                     """
                     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
                     """
                     if (minHei > maxHei):
                         raise ValueError, "Error selecting heights: Height range (%d,%d) is not valid" % (minHei, maxHei)
                     if (minHei < self.dataOut.heightList[0]):
                         minHei = self.dataOut.heightList[0]
                     if (maxHei > self.dataOut.heightList[-1]):
                         maxHei = self.dataOut.heightList[-1]
             #            raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
                     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)
                     self.selectHeightsByIndex(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 != None:
                         data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
                         #data_cspc = data_cspc[:,:,beacon_heiIndexList]
                     data_dc = None
                     if self.dataOut.data_dc != 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 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_spc
                         self.dataOut.data_cspc
                         self.dataOut.data_dc
                         self.dataOut.heightList
                     Return:
 si el metodo se ejecuto con exito caso contrario devuelve 0
                     """
                     if (minIndex < 0) or (minIndex > maxIndex):
                         raise ValueError, "Error selecting heights by index: Index range in (%d,%d) is not valid" % (minIndex, maxIndex)
                     if (maxIndex >= self.dataOut.nHeights):
                         maxIndex = self.dataOut.nHeights-1
             #            raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
             #         nHeights = maxIndex - minIndex + 1
                     #Spectra
                     data_spc = self.dataOut.data_spc[:,:,minIndex:maxIndex+1]
                     data_cspc = None
                     if self.dataOut.data_cspc != None:
                         data_cspc = self.dataOut.data_cspc[:,:,minIndex:maxIndex+1]
                     data_dc = None
                     if self.dataOut.data_dc != 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
                 def removeDC(self, mode = 2):
                     jspectra = self.dataOut.data_spc
                     jcspectra = self.dataOut.data_cspc
                     num_chan = jspectra.shape[0]
                     num_hei = jspectra.shape[2]
                     if jcspectra != None:
                         jcspectraExist = True
                         num_pairs = jcspectra.shape[0]
                     else:   jcspectraExist = False
                     freq_dc = jspectra.shape[1]/2
                     ind_vel = numpy.array([-2,-1,1,2]) + freq_dc
                     if ind_vel[0]<0:
                         ind_vel[range(0,1)] = ind_vel[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(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 1
                 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 == None:
                         count_hei = num_hei/2   #Como es entero no importa
                         hei_interf = numpy.asmatrix(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 = 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(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[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/ math.sqrt(num_incoh))
                         noiseid = noiseid[0]
                         cnoiseid = noiseid.size
                         interfid = numpy.where(jspc_interf > tmp_noise/ math.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/math.sqrt(num_incoh))).nonzero()
                                 cind = len(ind)
                                 if (cind > 0):
                                     jspectra[ich,interfid[ip],ind] = tmp_noise*(1 + (numpy.random.uniform(cind) - 0.5)/math.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(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/math.sqrt(num_incoh))).nonzero()
                         jspectra[ich,indAux[0],indAux[1]] = tmp_noise * (1 - 1/math.sqrt(num_incoh))
                     #Remocion de Interferencia en el Cross Spectra
                     if jcspectra == None: return jspectra, jcspectra
                     num_pairs = 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[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 = numpy.median(numpy.real(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
                         median_imag = numpy.median(numpy.imag(junkcspc_interf[mask_prof[ind[range(3*num_prof/4)]],:]))
                         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(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 setRadarFrequency(self, frequency=None):
                     if frequency != None:
                         self.dataOut.frequency = frequency
                     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,:,:]
                         noise[channel] = hildebrand_sekhon(daux, self.dataOut.nIncohInt)
                     self.dataOut.noise_estimation = noise.copy()
                     return 1
             class IncohInt(Operation):
                 __profIndex = 0
                 __withOverapping  = False
                 __byTime = False
                 __initime = None
                 __lastdatatime = None
                 __integrationtime = None
                 __buffer_spc = None
                 __buffer_cspc = None
                 __buffer_dc = None
                 __dataReady = False
                 __timeInterval = None
                 n = None
                 def __init__(self):
                     Operation.__init__(self)
             #         self.isConfig = False
                 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 = None
                     self.__buffer_cspc = None
                     self.__buffer_dc = None
                     self.__dataReady = False
                     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 #if (type(timeInterval)!=integer) -> change this line
                         self.n = 9999
                         self.__byTime = True
                     if overlapping:
                         self.__withOverapping = True
                     else:
                         self.__withOverapping = False
                         self.__buffer_spc = 0
                         self.__buffer_cspc = 0
                         self.__buffer_dc = 0
                     self.__profIndex = 0
                 def putData(self, data_spc, data_cspc, data_dc):
                     """
                     Add a profile to the __buffer_spc and increase in one the __profileIndex
                     """
                     if not self.__withOverapping:
                         self.__buffer_spc += data_spc
-                        if data_cspc == None:
+                        if data_cspc is None:
                             self.__buffer_cspc = None
                         else:
                             self.__buffer_cspc += data_cspc
-                        if data_dc == None:
+                        if data_dc is None:
                             self.__buffer_dc = None
                         else:
                             self.__buffer_dc += data_dc
                         self.__profIndex += 1
                         return
                     #Overlapping data
                     nChannels, nFFTPoints, nHeis = data_spc.shape
                     data_spc = numpy.reshape(data_spc, (1, nChannels, nFFTPoints, nHeis))
                     if data_cspc != None:
                         data_cspc = numpy.reshape(data_cspc, (1, -1, nFFTPoints, nHeis))
                     if data_dc != None:
                         data_dc = numpy.reshape(data_dc, (1, -1, nHeis))
                     #If the buffer is empty then it takes the data value
-                    if self.__buffer_spc == None:
+                    if self.__buffer_spc is None:
                         self.__buffer_spc = data_spc
-                        if data_cspc == None:
+                        if data_cspc is None:
                             self.__buffer_cspc = None
                         else:
                             self.__buffer_cspc += data_cspc
-                        if data_dc == None:
+                        if data_dc is None:
                             self.__buffer_dc = None
                         else:
                             self.__buffer_dc += data_dc
                         self.__profIndex += 1
                         return
                     #If the buffer length is lower than n then stakcing the data value
                     if self.__profIndex < self.n:
                         self.__buffer_spc = numpy.vstack((self.__buffer_spc, data_spc))
                         if data_cspc != None:
                             self.__buffer_cspc = numpy.vstack((self.__buffer_cspc, data_cspc))
                         if data_dc != None:
                             self.__buffer_dc = numpy.vstack((self.__buffer_dc, data_dc))
                         self.__profIndex += 1
                         return
                     #If the buffer length is equal to n then replacing the last buffer value with the data value
                     self.__buffer_spc = numpy.roll(self.__buffer_spc, -1, axis=0)
                     self.__buffer_spc[self.n-1] = data_spc
                     if data_cspc != None:
                         self.__buffer_cspc = numpy.roll(self.__buffer_cspc, -1, axis=0)
                         self.__buffer_cspc[self.n-1] = data_cspc
                     if data_dc != None:
                         self.__buffer_dc = numpy.roll(self.__buffer_dc, -1, axis=0)
                         self.__buffer_dc[self.n-1] = data_dc
                     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
                     """
                     data_spc = None
                     data_cspc  = None
                     data_dc = None
                     if not self.__withOverapping:
                         data_spc = self.__buffer_spc
                         data_cspc = self.__buffer_cspc
                         data_dc = self.__buffer_dc
                         n = self.__profIndex
                         self.__buffer_spc = 0
                         self.__buffer_cspc = 0
                         self.__buffer_dc = 0
                         self.__profIndex = 0
                         return data_spc, data_cspc, data_dc, n
                     #Integration with Overlapping
                     data_spc = numpy.sum(self.__buffer_spc, axis=0)
                     if self.__buffer_cspc != None:
                         data_cspc = numpy.sum(self.__buffer_cspc, axis=0)
                     if self.__buffer_dc != None:
                         data_dc = numpy.sum(self.__buffer_dc, axis=0)
                     n = self.__profIndex
                     return data_spc, data_cspc, data_dc, n
                 def byProfiles(self, *args):
                     self.__dataReady = False
                     avgdata_spc = None
                     avgdata_cspc = None
                     avgdata_dc = None
             #         n = None
                     self.putData(*args)
                     if self.__profIndex == self.n:
                         avgdata_spc, avgdata_cspc, avgdata_dc, n = self.pushData()
                         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
                     n = 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.__initime == None:
                         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)
                     self.__lastdatatime = datatime
-                    if avgdata_spc == None:
+                    if avgdata_spc is None:
                         return None, None, None, None
                     avgdatatime = self.__initime
                     try:
                         self.__timeInterval = (self.__lastdatatime - self.__initime)/(self.n - 1)
                     except:
                         self.__timeInterval = self.__lastdatatime - self.__initime
                     deltatime = datatime -self.__lastdatatime
                     if not self.__withOverapping:
                         self.__initime = datatime
                     else:
                         self.__initime += deltatime
                     return avgdatatime, avgdata_spc, avgdata_cspc, avgdata_dc
                 def run(self, dataOut, n=None, timeInterval=None, overlapping=False):
                     if n==1:
                         dataOut.flagNoData = False
                         return
                     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)
             #        dataOut.timeInterval *= n
                     dataOut.flagNoData = True
                     if self.__dataReady:
                         dataOut.data_spc = avgdata_spc
                         dataOut.data_cspc = avgdata_cspc
                         dataOut.data_dc = avgdata_dc
                         dataOut.nIncohInt *= self.n
                         dataOut.utctime = avgdatatime
                         #dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt * dataOut.nIncohInt * dataOut.nFFTPoints
             #             dataOut.timeInterval = self.__timeInterval*self.n
                         dataOut.flagNoData = False

schainpy/model/proc/jroproc_voltage.py +4 -4

             import numpy
             from jroproc_base import ProcessingUnit, Operation
             from schainpy.model.data.jrodata import Voltage
             class VoltageProc(ProcessingUnit):
                 def __init__(self):
                     ProcessingUnit.__init__(self)
             #         self.objectDict = {}
                     self.dataOut = Voltage()
                     self.flip = 1
                 def run(self):
                     if self.dataIn.type == 'AMISR':
                         self.__updateObjFromAmisrInput()
                     if self.dataIn.type == 'Voltage':
                         self.dataOut.copy(self.dataIn)
             #         self.dataOut.copy(self.dataIn)
                 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
             #
             #        pass#
             #
             #    def init(self):
             #
             #
             #        if self.dataIn.type == 'AMISR':
             #            self.__updateObjFromAmisrInput()
             #
             #        if self.dataIn.type == 'Voltage':
             #            self.dataOut.copy(self.dataIn)
             #        # No necesita copiar en cada init() los atributos de dataIn
             #        # la copia deberia hacerse por cada nuevo bloque de datos
                 def selectChannels(self, channelList):
                     channelIndexList = []
                     for channel in channelList:
                         if channel not in self.dataOut.channelList:
                             raise ValueError, "Channel %d is not in %s" %(channel, str(self.dataOut.channelList))
                         index = self.dataOut.channelList.index(channel)
                         channelIndexList.append(index)
                     self.selectChannelsByIndex(channelIndexList)
                 def selectChannelsByIndex(self, 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:
                         self.dataOut.data
                         self.dataOut.channelIndexList
                         self.dataOut.nChannels
                         self.dataOut.m_ProcessingHeader.totalSpectra
                         self.dataOut.systemHeaderObj.numChannels
                         self.dataOut.m_ProcessingHeader.blockSize
                     Return:
                         None
                     """
                     for channelIndex in channelIndexList:
                         if channelIndex not in self.dataOut.channelIndexList:
                             print channelIndexList
                             raise ValueError, "The value %d in channelIndexList is not valid" %channelIndex
             #         nChannels = len(channelIndexList)
                     if self.dataOut.flagDataAsBlock:
                         """
                         Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
                         """
                         data = self.dataOut.data[channelIndexList,:,:]
                     else:
                         data = self.dataOut.data[channelIndexList,:]
                     self.dataOut.data = data
                     self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
             #        self.dataOut.nChannels = nChannels
                     return 1
                 def selectHeights(self, minHei=None, maxHei=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
                     """
                     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):
                         raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
                     if (maxHei > self.dataOut.heightList[-1]):
                         maxHei = self.dataOut.heightList[-1]
             #            raise ValueError, "some value in (%d,%d) is not valid" % (minHei, maxHei)
                     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)
                     self.selectHeightsByIndex(minIndex, maxIndex)
                     return 1
                 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 (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
             #            raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex)
             #         nHeights = maxIndex - minIndex + 1
                     #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)
                     return 1
                 def filterByHeights(self, window):
                     deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
                     if window == None:
                         window = (self.dataOut.radarControllerHeaderObj.txA/self.dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
                     newdelta = deltaHeight * window
                     r = self.dataOut.nHeights % window
                     newheights = (self.dataOut.nHeights-r)/window
                     if newheights <= 1:
                         raise ValueError, "filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(self.dataOut.nHeights, window)
                     if self.dataOut.flagDataAsBlock:
                         """
                         Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
                         """
                         buffer = self.dataOut.data[:, :, 0:self.dataOut.nHeights-r]
                         buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nProfiles,self.dataOut.nHeights/window,window)
                         buffer = numpy.sum(buffer,3)
                     else:
                         buffer = self.dataOut.data[:,0:self.dataOut.nHeights-r]
                         buffer = buffer.reshape(self.dataOut.nChannels,self.dataOut.nHeights/window,window)
                         buffer = numpy.sum(buffer,2)
                     self.dataOut.data = buffer
                     self.dataOut.heightList = self.dataOut.heightList[0] + numpy.arange( newheights )*newdelta
                     self.dataOut.windowOfFilter = window
                 def setH0(self, h0, deltaHeight = None):
                     if not deltaHeight:
                         deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0]
                     nHeights = self.dataOut.nHeights
                     newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
                     self.dataOut.heightList = newHeiRange
                 def deFlip(self, channelList = []):
                     data = self.dataOut.data.copy()
                     if self.dataOut.flagDataAsBlock:
                         flip = self.flip
                         profileList = range(self.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 self.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 self.dataOut.channelList:
                                     continue
                                 data[thisChannel,:] = data[thisChannel,:]*self.flip
                         self.flip *= -1.
                     self.dataOut.data = data
                 def setRadarFrequency(self, frequency=None):
                     if frequency != None:
                         self.dataOut.frequency = frequency
                     return 1
             class CohInt(Operation):
                 isConfig = False
                 __profIndex = 0
                 __withOverapping  = False
                 __byTime = False
                 __initime = None
                 __lastdatatime = None
                 __integrationtime = None
                 __buffer = None
                 __dataReady = False
                 n = None
                 def __init__(self):
                     Operation.__init__(self)
             #         self.isConfig = False
                 def setup(self, n=None, timeInterval=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
                     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.__withOverapping = True
                         self.__buffer = None
                     else:
                         self.__withOverapping = 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.__withOverapping:
                         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 == None:
+                    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.__withOverapping:
                         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)
                     n = self.__profIndex
                     return data, n
                 def byProfiles(self, data):
                     self.__dataReady = False
                     avgdata = None
             #         n = None
                     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 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 == None:
+                    if avgdata is None:
                         return None, None
                     avgdatatime = self.__initime
                     deltatime = datatime -self.__lastdatatime
                     if not self.__withOverapping:
                         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, **kwargs):
                     if not self.isConfig:
                         self.setup(**kwargs)
                         self.isConfig = True
                     if dataOut.flagDataAsBlock:
                         """
                         Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
                         """
                         avgdata, avgdatatime = self.integrateByBlock(dataOut)
                     else:
                         avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime)
             #        dataOut.timeInterval *= n
                     dataOut.flagNoData = True
                     if self.__dataReady:
                         dataOut.data = avgdata
                         dataOut.nCohInt *= self.n
                         dataOut.utctime = avgdatatime
             #             dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
                         dataOut.flagNoData = False
             class Decoder(Operation):
                 isConfig = False
                 __profIndex = 0
                 code = None
                 nCode = None
                 nBaud = None
                 def __init__(self):
                     Operation.__init__(self)
                     self.times = None
                     self.osamp = None
             #         self.__setValues = False
                     self.isConfig = 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 dataOut.flagDataAsBlock:
                         self.ndatadec = self.__nHeis #- self.nBaud + 1
                         self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
                     else:
                         __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))
                         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)
                     datadec = data#[:,:]
                     return datadec
                 def convolutionInFreqOpt(self, data):
                     raise NotImplementedError
             #         fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
             #
             #         data = cfunctions.decoder(fft_code, data)
             #
             #         datadec = data#[:,:]
             #
             #         return datadec
                 def convolutionInTime(self, data):
                     code = self.code[self.__profIndex]
                     for i in range(self.__nChannels):
-                        self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='same')
+                        self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
                     return self.datadecTime
                 def convolutionByBlockInTime(self, data):
                     repetitions = 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))
                     for i in range(self.__nChannels):
                         for j in range(self.__nProfiles):
                             self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='same')
                     return self.datadecTime
                 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
                     if not self.isConfig:
-                        if code == None:
+                        if code is None:
                             code = dataOut.code
                         else:
                             code = numpy.array(code).reshape(nCode,nBaud)
                         self.setup(code, osamp, dataOut)
                         self.isConfig = True
                     if dataOut.flagDataAsBlock:
                         """
                         Decoding when data have been read as block,
                         """
                         datadec = self.convolutionByBlockInTime(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)
                     dataOut.code = self.code
                     dataOut.nCode = self.nCode
                     dataOut.nBaud = self.nBaud
                     dataOut.data = datadec
                     dataOut.heightList = dataOut.heightList[0:self.ndatadec]
                     dataOut.flagDecodeData = True #asumo q la data esta decodificada
                     if self.__profIndex == self.nCode-1:
                         self.__profIndex = 0
                         return 1
                     self.__profIndex += 1
                     return 1
             #        dataOut.flagDeflipData = True #asumo q la data no esta sin flip
             class ProfileConcat(Operation):
                 isConfig = False
                 buffer = None
                 def __init__(self):
                     Operation.__init__(self)
                     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.nHeights
                     self.start_index = 0
                     self.times = 1
                 def concat(self, data):
                     self.buffer[:,self.start_index:self.profiles*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
             class ProfileSelector(Operation):
                 profileIndex = None
                 # Tamanho total de los perfiles
                 nProfiles = None
                 def __init__(self):
                     Operation.__init__(self)
                     self.profileIndex = 0
                 def incIndex(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))
                     """
                     dataOut.flagNoData = True
                     if nProfiles:
                         self.nProfiles = dataOut.nProfiles
                     else:
                         self.nProfiles = nProfiles
                     if dataOut.flagDataAsBlock:
                         """
                         data dimension  = [nChannels, nProfiles, nHeis]
                         """
                         if profileList != None:
                             dataOut.data = dataOut.data[:,profileList,:]
                             dataOut.nProfiles = len(profileList)
                             dataOut.profileIndex = dataOut.nProfiles - 1
                         if profileRangeList != None:
                             minIndex = profileRangeList[0]
                             maxIndex = profileRangeList[1]
                             dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
                             dataOut.nProfiles = maxIndex - minIndex + 1
                             dataOut.profileIndex = dataOut.nProfiles - 1
                         if rangeList != None:
                             raise ValueError, "Profile Selector: Not implemented for rangeList yet"
                         dataOut.flagNoData = False
                         return True
                     else:
                         """
                         data dimension  = [nChannels, nHeis]
                         """
                         if profileList != None:
                             dataOut.nProfiles = len(profileList)
                             if self.isThisProfileInList(dataOut.profileIndex, profileList):
                                 dataOut.flagNoData = False
                                 dataOut.profileIndex = self.profileIndex
                                 self.incIndex()
                             return True
                         if profileRangeList != None:
                             minIndex = profileRangeList[0]
                             maxIndex = profileRangeList[1]
                             dataOut.nProfiles = maxIndex - minIndex + 1
                             if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
                                 dataOut.flagNoData = False
                                 dataOut.profileIndex = self.profileIndex
                                 self.incIndex()
                             return True
                         if rangeList != None:
                             nProfiles = 0
                             for thisRange in rangeList:
                                 minIndex = thisRange[0]
                                 maxIndex = thisRange[1]
                                 nProfiles += maxIndex - minIndex + 1
                             dataOut.nProfiles = nProfiles
                             for thisRange in rangeList:
                                 minIndex = thisRange[0]
                                 maxIndex = thisRange[1]
                                 if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
             #                         print "profileIndex = ", dataOut.profileIndex
                                     dataOut.flagNoData = False
                                     dataOut.profileIndex = self.profileIndex
                                     self.incIndex()
                                     break
                             return True
                         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.incIndex()
                             return 1
                     raise ValueError, "ProfileSelector needs profileList, profileRangeList or rangeList parameter"
                     return 0
             class Reshaper(Operation):
                 def __init__(self):
                     Operation.__init__(self)
                     self.updateNewHeights = True
                 def run(self, dataOut, shape):
                     if not dataOut.flagDataAsBlock:
                         raise ValueError, "Reshaper can only be used when voltage have been read as Block, getBlock = True"
                     if len(shape) != 3:
                         raise ValueError, "shape len should be equal to 3, (nChannels, nProfiles, nHeis)"
                     shape_tuple = tuple(shape)
                     dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
                     dataOut.flagNoData = False
                     if self.updateNewHeights:
                         old_nheights = dataOut.nHeights
                         new_nheights = dataOut.data.shape[2]
                         factor = 1.0*new_nheights / old_nheights
                         deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                         xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * factor
                         dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
                         dataOut.nProfiles = dataOut.data.shape[1]
                         dataOut.ippSeconds *= factor
             import collections
             from scipy.stats import mode
             class Synchronize(Operation):
                 isConfig = False
                 __profIndex = 0
                 def __init__(self):
                     Operation.__init__(self)
             #         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

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