diff --git a/schainpy/controller.py b/schainpy/controller.py index 2fe8625..cf357aa 100644 --- a/schainpy/controller.py +++ b/schainpy/controller.py @@ -6,11 +6,8 @@ from xml.etree.ElementTree import Element, SubElement, ElementTree from xml.etree import ElementTree as ET from xml.dom import minidom -import sys import datetime -from model.jrodataIO import * -from model.jroprocessing import * -from model.jroplot import * +from model import * def prettify(elem): """Return a pretty-printed XML string for the Element. @@ -218,7 +215,7 @@ class OperationConf(): if self.type == 'self': raise ValueError, "This operation type cannot be created" - if self.type == 'other': + if self.type == 'external' or self.type == 'other': className = eval(self.name) opObj = className() @@ -287,7 +284,7 @@ class ProcUnitConf(): self.opConfObjList = [] - self.addOperation(name='init', optype='self') + self.addOperation(name='run', optype='self') def addParameter(self, **kwargs): @@ -377,7 +374,10 @@ class ProcUnitConf(): kwargs[parmConfObj.name] = parmConfObj.getValue() #print "\tRunning the '%s' operation with %s" %(opConfObj.name, opConfObj.id) - sts = self.procUnitObj.call(opConfObj, **kwargs) + sts = self.procUnitObj.call(opType = opConfObj.type, + opName = opConfObj.name, + opId = opConfObj.id, + **kwargs) finalSts = finalSts or sts return finalSts @@ -406,7 +406,7 @@ class ReadUnitConf(ProcUnitConf): return self.ELEMENTNAME - def setup(self, id, name, datatype, path, startDate, endDate, startTime, endTime, **kwargs): + def setup(self, id, name, datatype, path="", startDate="", endDate="", startTime="", endTime="", **kwargs): self.id = id self.name = name @@ -471,13 +471,13 @@ class Project(): self.name = name self.description = description - def addReadUnit(self, datatype, path, startDate='', endDate='', startTime='', endTime='', **kwargs): + def addReadUnit(self, datatype, **kwargs): id = self.__getNewId() - name = '%sReader' %(datatype) + name = '%s' %(datatype) readUnitConfObj = ReadUnitConf() - readUnitConfObj.setup(id, name, datatype, path, startDate, endDate, startTime, endTime, **kwargs) + readUnitConfObj.setup(id, name, datatype, **kwargs) self.procUnitConfObjDict[readUnitConfObj.getId()] = readUnitConfObj @@ -486,7 +486,7 @@ class Project(): def addProcUnit(self, datatype, inputId): id = self.__getNewId() - name = '%sProc' %(datatype) + name = '%s' %(datatype) procUnitConfObj = ProcUnitConf() procUnitConfObj.setup(id, name, datatype, inputId) @@ -569,25 +569,27 @@ class Project(): for procUnitConfObj in self.procUnitConfObjDict.values(): procUnitConfObj.createObjects() - def __connect(self, objIN, obj): + def __connect(self, objIN, thisObj): - obj.setInput(objIN.getOutput()) + thisObj.setInput(objIN.getOutputObj()) def connectObjects(self): - for puConfObj in self.procUnitConfObjDict.values(): + for thisPUConfObj in self.procUnitConfObjDict.values(): - inputId = puConfObj.getInputId() + inputId = thisPUConfObj.getInputId() if int(inputId) == 0: continue + #Get input object puConfINObj = self.procUnitConfObjDict[inputId] + puObjIN = puConfINObj.getProcUnitObj() - puObj = puConfObj.getProcUnitObj() - puINObj = puConfINObj.getProcUnitObj() + #Get current object + thisPUObj = thisPUConfObj.getProcUnitObj() - self.__connect(puINObj, puObj) + self.__connect(puObjIN, thisPUObj) def run(self): @@ -637,7 +639,7 @@ if __name__ == '__main__': opObj10.addParameter(name='minHei', value='90', format='float') opObj10.addParameter(name='maxHei', value='180', format='float') - opObj12 = procUnitConfObj0.addOperation(name='CohInt', optype='other') + opObj12 = procUnitConfObj0.addOperation(name='CohInt', optype='external') opObj12.addParameter(name='n', value='10', format='int') procUnitConfObj1 = controllerObj.addProcUnit(datatype='Spectra', inputId=procUnitConfObj0.getId()) @@ -645,14 +647,14 @@ if __name__ == '__main__': # procUnitConfObj1.addParameter(name='pairList', value='(0,1),(0,2),(1,2)', format='') - opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other') + opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='external') opObj11.addParameter(name='idfigure', value='1', format='int') opObj11.addParameter(name='wintitle', value='SpectraPlot0', format='str') opObj11.addParameter(name='zmin', value='40', format='int') opObj11.addParameter(name='zmax', value='90', format='int') opObj11.addParameter(name='showprofile', value='1', format='int') -# opObj11 = procUnitConfObj1.addOperation(name='CrossSpectraPlot', optype='other') +# opObj11 = procUnitConfObj1.addOperation(name='CrossSpectraPlot', optype='external') # opObj11.addParameter(name='idfigure', value='2', format='int') # opObj11.addParameter(name='wintitle', value='CrossSpectraPlot', format='str') # opObj11.addParameter(name='zmin', value='40', format='int') @@ -661,21 +663,21 @@ if __name__ == '__main__': # procUnitConfObj2 = controllerObj.addProcUnit(datatype='Voltage', inputId=procUnitConfObj0.getId()) # -# opObj12 = procUnitConfObj2.addOperation(name='CohInt', optype='other') +# opObj12 = procUnitConfObj2.addOperation(name='CohInt', optype='external') # opObj12.addParameter(name='n', value='2', format='int') # opObj12.addParameter(name='overlapping', value='1', format='int') # # procUnitConfObj3 = controllerObj.addProcUnit(datatype='Spectra', inputId=procUnitConfObj2.getId()) # procUnitConfObj3.addParameter(name='nFFTPoints', value='32', format='int') # -# opObj11 = procUnitConfObj3.addOperation(name='SpectraPlot', optype='other') +# opObj11 = procUnitConfObj3.addOperation(name='SpectraPlot', optype='external') # opObj11.addParameter(name='idfigure', value='2', format='int') # opObj11.addParameter(name='wintitle', value='SpectraPlot1', format='str') # opObj11.addParameter(name='zmin', value='40', format='int') # opObj11.addParameter(name='zmax', value='90', format='int') # opObj11.addParameter(name='showprofile', value='1', format='int') -# opObj11 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other') +# opObj11 = procUnitConfObj1.addOperation(name='RTIPlot', optype='external') # opObj11.addParameter(name='idfigure', value='10', format='int') # opObj11.addParameter(name='wintitle', value='RTI', format='str') ## opObj11.addParameter(name='xmin', value='21', format='float') @@ -688,10 +690,10 @@ if __name__ == '__main__': # opObj10 = procUnitConfObj1.addOperation(name='selectChannels') # opObj10.addParameter(name='channelList', value='0,2,4,6', format='intlist') # -# opObj12 = procUnitConfObj1.addOperation(name='IncohInt', optype='other') +# opObj12 = procUnitConfObj1.addOperation(name='IncohInt', optype='external') # opObj12.addParameter(name='n', value='2', format='int') # -# opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other') +# opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='external') # opObj11.addParameter(name='idfigure', value='2', format='int') # opObj11.addParameter(name='wintitle', value='SpectraPlot10', format='str') # opObj11.addParameter(name='zmin', value='70', format='int') @@ -700,10 +702,10 @@ if __name__ == '__main__': # opObj10 = procUnitConfObj1.addOperation(name='selectChannels') # opObj10.addParameter(name='channelList', value='2,6', format='intlist') # -# opObj12 = procUnitConfObj1.addOperation(name='IncohInt', optype='other') +# opObj12 = procUnitConfObj1.addOperation(name='IncohInt', optype='external') # opObj12.addParameter(name='n', value='2', format='int') # -# opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other') +# opObj11 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='external') # opObj11.addParameter(name='idfigure', value='3', format='int') # opObj11.addParameter(name='wintitle', value='SpectraPlot10', format='str') # opObj11.addParameter(name='zmin', value='70', format='int') @@ -720,10 +722,10 @@ if __name__ == '__main__': # procUnitConfObj2 = controllerObj.addProcUnit(datatype='Spectra', inputId=procUnitConfObj1.getId()) # -# opObj21 = procUnitConfObj2.addOperation(name='IncohInt', optype='other') +# opObj21 = procUnitConfObj2.addOperation(name='IncohInt', optype='external') # opObj21.addParameter(name='n', value='2', format='int') # -# opObj11 = procUnitConfObj2.addOperation(name='SpectraPlot', optype='other') +# opObj11 = procUnitConfObj2.addOperation(name='SpectraPlot', optype='external') # opObj11.addParameter(name='idfigure', value='4', format='int') # opObj11.addParameter(name='wintitle', value='SpectraPlot OBJ 2', format='str') # opObj11.addParameter(name='zmin', value='70', format='int') diff --git a/schainpy/model/__init__.py b/schainpy/model/__init__.py index 2041e82..576cc69 100644 --- a/schainpy/model/__init__.py +++ b/schainpy/model/__init__.py @@ -1,3 +1,4 @@ -import jrodata -import jrodataIO -import jroprocessing +from model.data.jrodata import * +from model.io.jrodataIO import * +from model.proc.jroprocessing import * +from model.graphics.jroplot import * diff --git a/schainpy/model/data/__init__.py b/schainpy/model/data/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/schainpy/model/data/__init__.py diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py new file mode 100644 index 0000000..feba054 --- /dev/null +++ b/schainpy/model/data/jrodata.py @@ -0,0 +1,782 @@ +''' + +$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',' 2: + nums_min = lenOfData/10 + else: + nums_min = 2 + + sump = 0. + + sumq = 0. + + j = 0 + + cont = 1 + + while((cont==1)and(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 + + lnoise = sump /j + stdv = numpy.sqrt((sumq - lnoise**2)/(j - 1)) + return lnoise + +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 + + flagTimeBlock = 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 + + 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) + + 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) + + 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) + +class Voltage(JROData): + + #data es un numpy array de 2 dmensiones (canales, alturas) + data = None + + def __init__(self): + ''' + Constructor + ''' + + 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.flagTimeBlock = 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 + + + def getNoisebyHildebrand(self): + """ + Determino el nivel de ruido usando el metodo Hildebrand-Sekhon + + Return: + noiselevel + """ + + for channel in range(self.nChannels): + daux = self.data_spc[channel,:,:] + self.noise[channel] = hildebrand_sekhon(daux, self.nCohInt) + + return self.noise + + def getNoise(self, type = 1): + + self.noise = numpy.zeros(self.nChannels) + + if type == 1: + noise = self.getNoisebyHildebrand() + + return 10*numpy.log10(noise) + + noise = property(getNoise, "I'm the 'nHeights' 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 + + def __init__(self): + ''' + Constructor + ''' + + 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.flagTimeBlock = 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): + """ + 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,:,:] + noise[channel] = hildebrand_sekhon(daux, self.nIncohInt) + + return noise + + def getNoise(self): + if self.noise_estimation != None: + return self.noise_estimation #this was estimated by getNoise Operation defined in jroproc_spectra.py + else: + noise = self.getNoisebyHildebrand() + 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: + return True + + return False + + def getFlagDc(self): + + if self.data_dc == None: + return True + + return False + + 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.") + +class SpectraHeis(Spectra): + + data_spc = None + + data_cspc = None + + data_dc = None + + nFFTPoints = None + +# nPairs = None + + pairsList = 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.flagTimeBlock = False + +# self.nPairs = 0 + + self.utctime = None + + self.blocksize = None + +class Fits: + + heightList = None + + channelList = None + + flagNoData = True + + flagTimeBlock = 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 = None + + self.nIncohInt = None + + self.useLocalTime = True + +# 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): + + self.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 + + 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") + + ltctime = property(getltctime, "I'm the 'ltctime' property") + +class AMISR: + def __init__(self): + self.flagNoData = True + self.data = None + self.utctime = None + self.type = "AMISR" + + #propiedades para compatibilidad con Voltages + self.timeZone = 300#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',' 0: +# + fp.seek(endFp) + + except Exception, e: + print "RadarControllerHeader: " + e + return 0 + + return 1 + + def write(self, fp): + headerTuple = (self.size, + self.expType, + self.nTx, + self.ipp, + self.txA, + self.txB, + self.nWindows, + self.numTaus, + self.codeType, + self.line6Function, + self.line5Function, + self.fClock, + self.prePulseBefore, + self.prePulserAfter, + self.rangeIpp, + self.rangeTxA, + self.rangeTxB) + + header = numpy.array(headerTuple,RADAR_STRUCTURE) + header.tofile(fp) + + #dynamic = self.dynamic + #dynamic.tofile(fp) + + samplingWindow = self.samplingWindow + samplingWindow.tofile(fp) + + if self.numTaus > 0: + self.Taus.tofile(fp) + + nCode = numpy.array(self.nCode, ' 0: + self.flag_cspc = True + +# except Exception, e: +# print "Error ProcessingHeader: " +# return 0 + + return 1 + + def write(self, fp): + headerTuple = (self.size, + self.dtype, + self.blockSize, + self.profilesPerBlock, + self.dataBlocksPerFile, + self.nWindows, + self.processFlags, + self.nCohInt, + self.nIncohInt, + self.totalSpectra) + + header = numpy.array(headerTuple,PROCESSING_STRUCTURE) + header.tofile(fp) + + if self.nWindows != 0: + sampleWindowTuple = (self.firstHeight,self.deltaHeight,self.samplesWin) + samplingWindow = numpy.array(sampleWindowTuple,SAMPLING_STRUCTURE) + samplingWindow.tofile(fp) + + + if self.totalSpectra != 0: + spectraComb = numpy.array([],numpy.dtype('u1')) + spectraComb = self.spectraComb + spectraComb.tofile(fp) + +# if self.processFlags & PROCFLAG.DEFINE_PROCESS_CODE == PROCFLAG.DEFINE_PROCESS_CODE: +# nCode = numpy.array([self.nCode], numpy.dtype('u4')) #Probar con un dato que almacene codigo, hasta el momento no se hizo la prueba +# nCode.tofile(fp) +# +# nBaud = numpy.array([self.nBaud], numpy.dtype('u4')) +# nBaud.tofile(fp) +# +# code = self.code.reshape(self.nCode*self.nBaud) +# code = code.astype(numpy.dtype('= 30.: + return False + return True + + class FTP_Thread (threading.Thread): def __init__(self): threading.Thread.__init__(self) @@ -288,6 +296,10 @@ class Figure: raise ValueError, "This method is not implemented" + def close(self): + + self.__driver.show(True) + axesList = property(getAxesObjList) diff --git a/schainpy/model/graphics/jroplot.py b/schainpy/model/graphics/jroplot.py new file mode 100644 index 0000000..7a53039 --- /dev/null +++ b/schainpy/model/graphics/jroplot.py @@ -0,0 +1,3 @@ +from jroplot_voltage import * +from jroplot_spectra import * +from jroplot_heispectra import * diff --git a/schainpy/model/graphics/jroplot_heispectra.py b/schainpy/model/graphics/jroplot_heispectra.py new file mode 100644 index 0000000..67bd648 --- /dev/null +++ b/schainpy/model/graphics/jroplot_heispectra.py @@ -0,0 +1,318 @@ +''' + +@author: Daniel Suarez +''' + +import os +import datetime +import numpy + +from figure import Figure, isRealtime + +class SpectraHeisScope(Figure): + + + isConfig = None + __nsubplots = None + + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'spc' + + def __init__(self): + + self.isConfig = False + self.__nsubplots = 1 + + self.WIDTH = 230 + self.HEIGHT = 250 + self.WIDTHPROF = 120 + self.HEIGHTPROF = 0 + self.counter_imagwr = 0 + + def getSubplots(self): + + ncol = int(numpy.sqrt(self.nplots)+0.9) + nrow = int(self.nplots*1./ncol + 0.9) + + return nrow, ncol + + def setup(self, id, nplots, wintitle, show): + + showprofile = False + self.__showprofile = showprofile + self.nplots = nplots + + ncolspan = 1 + colspan = 1 + if showprofile: + ncolspan = 3 + colspan = 2 + self.__nsubplots = 2 + + self.createFigure(id = id, + wintitle = wintitle, + widthplot = self.WIDTH + self.WIDTHPROF, + heightplot = self.HEIGHT + self.HEIGHTPROF, + show = show) + + nrow, ncol = self.getSubplots() + + counter = 0 + for y in range(nrow): + for x in range(ncol): + + if counter >= self.nplots: + break + + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1) + + if showprofile: + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1) + + counter += 1 + + + def run(self, dataOut, id, wintitle="", channelList=None, + xmin=None, xmax=None, ymin=None, ymax=None, save=False, + figpath='./', figfile=None, ftp=False, wr_period=1, show=True, + server=None, folder=None, username=None, password=None): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + """ + + if dataOut.realtime: + if not(isRealtime(utcdatatime = dataOut.utctime)): + print 'Skipping this plot function' + return + + if channelList == None: + channelIndexList = dataOut.channelIndexList + else: + channelIndexList = [] + for channel in channelList: + if channel not in dataOut.channelList: + raise ValueError, "Channel %d is not in dataOut.channelList" + channelIndexList.append(dataOut.channelList.index(channel)) + +# x = dataOut.heightList + c = 3E8 + deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] + #deberia cambiar para el caso de 1Mhz y 100KHz + x = numpy.arange(-1*dataOut.nHeights/2.,dataOut.nHeights/2.)*(c/(2*deltaHeight*dataOut.nHeights*1000)) + #para 1Mhz descomentar la siguiente linea + #x= x/(10000.0) +# y = dataOut.data[channelIndexList,:] * numpy.conjugate(dataOut.data[channelIndexList,:]) +# y = y.real + datadB = 10.*numpy.log10(dataOut.data_spc) + y = datadB + + #thisDatetime = dataOut.datatime + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "" + #para 1Mhz descomentar la siguiente linea + #xlabel = "Frequency x 10000" + ylabel = "Intensity (dB)" + + if not self.isConfig: + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle=wintitle, + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = numpy.nanmin(y) + if ymax == None: ymax = numpy.nanmax(y) + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(len(self.axesList)): + ychannel = y[i,:] + str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S")) + title = "Channel %d: %4.2fdB: %s" %(i, numpy.max(ychannel), str_datetime) + axes = self.axesList[i] + axes.pline(x, ychannel, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title, grid='both') + + + self.draw() + + if save: + date = thisDatetime.strftime("%Y%m%d_%H%M%S") + if figfile == None: + figfile = self.getFilename(name = date) + + self.saveFigure(figpath, figfile) + + self.counter_imagwr += 1 + if (ftp and (self.counter_imagwr==wr_period)): + ftp_filename = os.path.join(figpath,figfile) + self.sendByFTP_Thread(ftp_filename, server, folder, username, password) + self.counter_imagwr = 0 + +class RTIfromSpectraHeis(Figure): + + isConfig = None + __nsubplots = None + + PREFIX = 'rtinoise' + + def __init__(self): + + self.timerange = 24*60*60 + self.isConfig = False + self.__nsubplots = 1 + + self.WIDTH = 820 + self.HEIGHT = 200 + self.WIDTHPROF = 120 + self.HEIGHTPROF = 0 + self.counter_imagwr = 0 + self.xdata = None + self.ydata = None + + def getSubplots(self): + + ncol = 1 + nrow = 1 + + return nrow, ncol + + def setup(self, id, nplots, wintitle, showprofile=True, show=True): + + self.__showprofile = showprofile + self.nplots = nplots + + ncolspan = 7 + colspan = 6 + self.__nsubplots = 2 + + self.createFigure(id = id, + wintitle = wintitle, + widthplot = self.WIDTH+self.WIDTHPROF, + heightplot = self.HEIGHT+self.HEIGHTPROF, + show = show) + + nrow, ncol = self.getSubplots() + + self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1) + + + def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True', + xmin=None, xmax=None, ymin=None, ymax=None, + timerange=None, + save=False, figpath='./', figfile=None, ftp=False, wr_period=1, show=True, + server=None, folder=None, username=None, password=None): + + if channelList == None: + channelIndexList = dataOut.channelIndexList + channelList = dataOut.channelList + else: + channelIndexList = [] + for channel in channelList: + if channel not in dataOut.channelList: + raise ValueError, "Channel %d is not in dataOut.channelList" + channelIndexList.append(dataOut.channelList.index(channel)) + + if timerange != None: + self.timerange = timerange + + tmin = None + tmax = None + x = dataOut.getTimeRange() + y = dataOut.getHeiRange() + + + data = dataOut.data_spc + data = numpy.average(data,axis=1) + datadB = 10*numpy.log10(data) + +# factor = dataOut.normFactor +# noise = dataOut.getNoise()/factor +# noisedB = 10*numpy.log10(noise) + + #thisDatetime = dataOut.datatime + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " RTI: %s" %(thisDatetime.strftime("%d-%b-%Y")) + xlabel = "Local Time" + ylabel = "Intensity (dB)" + + if not self.isConfig: + + nplots = 1 + + self.setup(id=id, + nplots=nplots, + wintitle=wintitle, + showprofile=showprofile, + show=show) + + tmin, tmax = self.getTimeLim(x, xmin, xmax) + if ymin == None: ymin = numpy.nanmin(datadB) + if ymax == None: ymax = numpy.nanmax(datadB) + + self.name = thisDatetime.strftime("%Y%m%d_%H%M%S") + self.isConfig = True + + self.xdata = numpy.array([]) + self.ydata = numpy.array([]) + + self.setWinTitle(title) + + +# title = "RTI %s" %(thisDatetime.strftime("%d-%b-%Y")) + title = "RTI - %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + + legendlabels = ["channel %d"%idchannel for idchannel in channelList] + axes = self.axesList[0] + + self.xdata = numpy.hstack((self.xdata, x[0:1])) + + if len(self.ydata)==0: + self.ydata = datadB[channelIndexList].reshape(-1,1) + else: + self.ydata = numpy.hstack((self.ydata, datadB[channelIndexList].reshape(-1,1))) + + + axes.pmultilineyaxis(x=self.xdata, y=self.ydata, + xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='.', markersize=8, linestyle="solid", grid='both', + XAxisAsTime=True + ) + + self.draw() + + if save: + + if figfile == None: + figfile = self.getFilename(name = self.name) + + self.saveFigure(figpath, figfile) + + self.counter_imagwr += 1 + if (ftp and (self.counter_imagwr==wr_period)): + ftp_filename = os.path.join(figpath,figfile) + self.sendByFTP_Thread(ftp_filename, server, folder, username, password) + self.counter_imagwr = 0 + + if x[1] + (x[1]-x[0]) >= self.axesList[0].xmax: + self.isConfig = False + del self.xdata + del self.ydata diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py new file mode 100644 index 0000000..efdd9f6 --- /dev/null +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -0,0 +1,1196 @@ +''' +@author: Daniel Suarez +''' +import os +import datetime +import numpy + +from figure import Figure, isRealtime + +class SpectraPlot(Figure): + + isConfig = None + __nsubplots = None + + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'spc' + + def __init__(self): + + self.isConfig = False + self.__nsubplots = 1 + + self.WIDTH = 280 + self.HEIGHT = 250 + self.WIDTHPROF = 120 + self.HEIGHTPROF = 0 + self.counter_imagwr = 0 + + self.PLOT_CODE = 1 + self.FTP_WEI = None + self.EXP_CODE = None + self.SUB_EXP_CODE = None + self.PLOT_POS = None + + def getSubplots(self): + + ncol = int(numpy.sqrt(self.nplots)+0.9) + nrow = int(self.nplots*1./ncol + 0.9) + + return nrow, ncol + + def setup(self, id, nplots, wintitle, showprofile=True, show=True): + + self.__showprofile = showprofile + self.nplots = nplots + + ncolspan = 1 + colspan = 1 + if showprofile: + ncolspan = 3 + colspan = 2 + self.__nsubplots = 2 + + self.createFigure(id = id, + wintitle = wintitle, + widthplot = self.WIDTH + self.WIDTHPROF, + heightplot = self.HEIGHT + self.HEIGHTPROF, + show=show) + + nrow, ncol = self.getSubplots() + + counter = 0 + for y in range(nrow): + for x in range(ncol): + + if counter >= self.nplots: + break + + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1) + + if showprofile: + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1) + + counter += 1 + + def run(self, dataOut, id, wintitle="", channelList=None, showprofile=True, + xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, + save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1, + server=None, folder=None, username=None, password=None, + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + showProfile : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + zmin : None, + zmax : None + """ + + if dataOut.flagNoData: + return None + + if realtime: + if not(isRealtime(utcdatatime = dataOut.utctime)): + print 'Skipping this plot function' + return + + if channelList == None: + channelIndexList = dataOut.channelIndexList + else: + channelIndexList = [] + for channel in channelList: + if channel not in dataOut.channelList: + raise ValueError, "Channel %d is not in dataOut.channelList" + channelIndexList.append(dataOut.channelList.index(channel)) + + factor = dataOut.normFactor + + x = dataOut.getVelRange(1) + y = dataOut.getHeiRange() + + z = dataOut.data_spc[channelIndexList,:,:]/factor + z = numpy.where(numpy.isfinite(z), z, numpy.NAN) + avg = numpy.average(z, axis=1) + avg = numpy.nanmean(z, axis=1) + noise = dataOut.noise/factor + + zdB = 10*numpy.log10(z) + avgdB = 10*numpy.log10(avg) + noisedB = 10*numpy.log10(noise) + + #thisDatetime = dataOut.datatime + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " Spectra" + xlabel = "Velocity (m/s)" + ylabel = "Range (Km)" + + if not self.isConfig: + + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle=wintitle, + showprofile=showprofile, + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = numpy.nanmin(y) + if ymax == None: ymax = numpy.nanmax(y) + if zmin == None: zmin = numpy.nanmin(avgdB)*0.9 + if zmax == None: zmax = numpy.nanmax(avgdB)*0.9 + + self.FTP_WEI = ftp_wei + self.EXP_CODE = exp_code + self.SUB_EXP_CODE = sub_exp_code + self.PLOT_POS = plot_pos + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(self.nplots): + str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S")) + title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i]+1, noisedB[i], str_datetime) + axes = self.axesList[i*self.__nsubplots] + axes.pcolor(x, y, zdB[i,:,:], + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, + xlabel=xlabel, ylabel=ylabel, title=title, + ticksize=9, cblabel='') + + if self.__showprofile: + axes = self.axesList[i*self.__nsubplots +1] + axes.pline(avgdB[i], y, + xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax, + xlabel='dB', ylabel='', title='', + ytick_visible=False, + grid='x') + + noiseline = numpy.repeat(noisedB[i], len(y)) + axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2) + + self.draw() + + if save: + + self.counter_imagwr += 1 + if (self.counter_imagwr==wr_period): + if figfile == None: + str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + figfile = self.getFilename(name = str_datetime) + + self.saveFigure(figpath, figfile) + + if ftp: + #provisionalmente envia archivos en el formato de la web en tiempo real + name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS) + path = '%s%03d' %(self.PREFIX, self.id) + ftp_file = os.path.join(path,'ftp','%s.png'%name) + self.saveFigure(figpath, ftp_file) + ftp_filename = os.path.join(figpath,ftp_file) + self.sendByFTP_Thread(ftp_filename, server, folder, username, password) + self.counter_imagwr = 0 + + + self.counter_imagwr = 0 + +class CrossSpectraPlot(Figure): + + isConfig = None + __nsubplots = None + + WIDTH = None + HEIGHT = None + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'cspc' + + def __init__(self): + + self.isConfig = False + self.__nsubplots = 4 + self.counter_imagwr = 0 + self.WIDTH = 250 + self.HEIGHT = 250 + self.WIDTHPROF = 0 + self.HEIGHTPROF = 0 + + self.PLOT_CODE = 1 + self.FTP_WEI = None + self.EXP_CODE = None + self.SUB_EXP_CODE = None + self.PLOT_POS = None + + def getSubplots(self): + + ncol = 4 + nrow = self.nplots + + return nrow, ncol + + def setup(self, id, nplots, wintitle, showprofile=True, show=True): + + self.__showprofile = showprofile + self.nplots = nplots + + ncolspan = 1 + colspan = 1 + + self.createFigure(id = id, + wintitle = wintitle, + widthplot = self.WIDTH + self.WIDTHPROF, + heightplot = self.HEIGHT + self.HEIGHTPROF, + show=True) + + nrow, ncol = self.getSubplots() + + counter = 0 + for y in range(nrow): + for x in range(ncol): + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1) + + counter += 1 + + def run(self, dataOut, id, wintitle="", pairsList=None, + xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, + save=False, figpath='./', figfile=None, ftp=False, wr_period=1, + power_cmap='jet', coherence_cmap='jet', phase_cmap='RdBu_r', show=True, + server=None, folder=None, username=None, password=None, + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + showProfile : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + zmin : None, + zmax : None + """ + + if pairsList == None: + pairsIndexList = dataOut.pairsIndexList + else: + pairsIndexList = [] + for pair in pairsList: + if pair not in dataOut.pairsList: + raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair) + pairsIndexList.append(dataOut.pairsList.index(pair)) + + if pairsIndexList == []: + return + + if len(pairsIndexList) > 4: + pairsIndexList = pairsIndexList[0:4] + factor = dataOut.normFactor + x = dataOut.getVelRange(1) + y = dataOut.getHeiRange() + z = dataOut.data_spc[:,:,:]/factor +# z = numpy.where(numpy.isfinite(z), z, numpy.NAN) + avg = numpy.abs(numpy.average(z, axis=1)) + noise = dataOut.noise()/factor + + zdB = 10*numpy.log10(z) + avgdB = 10*numpy.log10(avg) + noisedB = 10*numpy.log10(noise) + + + #thisDatetime = dataOut.datatime + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " Cross-Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Velocity (m/s)" + ylabel = "Range (Km)" + + if not self.isConfig: + + nplots = len(pairsIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle=wintitle, + showprofile=False, + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = numpy.nanmin(y) + if ymax == None: ymax = numpy.nanmax(y) + if zmin == None: zmin = numpy.nanmin(avgdB)*0.9 + if zmax == None: zmax = numpy.nanmax(avgdB)*0.9 + + self.FTP_WEI = ftp_wei + self.EXP_CODE = exp_code + self.SUB_EXP_CODE = sub_exp_code + self.PLOT_POS = plot_pos + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(self.nplots): + pair = dataOut.pairsList[pairsIndexList[i]] + str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S")) + title = "Ch%d: %4.2fdB: %s" %(pair[0], noisedB[pair[0]], str_datetime) + zdB = 10.*numpy.log10(dataOut.data_spc[pair[0],:,:]/factor) + axes0 = self.axesList[i*self.__nsubplots] + axes0.pcolor(x, y, zdB, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, + xlabel=xlabel, ylabel=ylabel, title=title, + ticksize=9, colormap=power_cmap, cblabel='') + + title = "Ch%d: %4.2fdB: %s" %(pair[1], noisedB[pair[1]], str_datetime) + zdB = 10.*numpy.log10(dataOut.data_spc[pair[1],:,:]/factor) + axes0 = self.axesList[i*self.__nsubplots+1] + axes0.pcolor(x, y, zdB, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, + xlabel=xlabel, ylabel=ylabel, title=title, + ticksize=9, colormap=power_cmap, cblabel='') + + coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[pair[0],:,:]*dataOut.data_spc[pair[1],:,:]) + coherence = numpy.abs(coherenceComplex) +# phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi + phase = numpy.arctan2(coherenceComplex.imag, coherenceComplex.real)*180/numpy.pi + + title = "Coherence %d%d" %(pair[0], pair[1]) + axes0 = self.axesList[i*self.__nsubplots+2] + axes0.pcolor(x, y, coherence, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=0, zmax=1, + xlabel=xlabel, ylabel=ylabel, title=title, + ticksize=9, colormap=coherence_cmap, cblabel='') + + title = "Phase %d%d" %(pair[0], pair[1]) + axes0 = self.axesList[i*self.__nsubplots+3] + axes0.pcolor(x, y, phase, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=-180, zmax=180, + xlabel=xlabel, ylabel=ylabel, title=title, + ticksize=9, colormap=phase_cmap, cblabel='') + + + + self.draw() + + if save: + + self.counter_imagwr += 1 + if (self.counter_imagwr==wr_period): + if figfile == None: + str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + figfile = self.getFilename(name = str_datetime) + + self.saveFigure(figpath, figfile) + + if ftp: + #provisionalmente envia archivos en el formato de la web en tiempo real + name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS) + path = '%s%03d' %(self.PREFIX, self.id) + ftp_file = os.path.join(path,'ftp','%s.png'%name) + self.saveFigure(figpath, ftp_file) + ftp_filename = os.path.join(figpath,ftp_file) + + try: + self.sendByFTP(ftp_filename, server, folder, username, password) + except: + self.counter_imagwr = 0 + print ValueError, 'Error FTP' + + self.counter_imagwr = 0 + +class RTIPlot(Figure): + + isConfig = None + __nsubplots = None + + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'rti' + + def __init__(self): + + self.timerange = 2*60*60 + self.isConfig = False + self.__nsubplots = 1 + + self.WIDTH = 800 + self.HEIGHT = 150 + self.WIDTHPROF = 120 + self.HEIGHTPROF = 0 + self.counter_imagwr = 0 + + self.PLOT_CODE = 0 + self.FTP_WEI = None + self.EXP_CODE = None + self.SUB_EXP_CODE = None + self.PLOT_POS = None + self.tmin = None + self.tmax = None + + self.xmin = None + self.xmax = None + + def getSubplots(self): + + ncol = 1 + nrow = self.nplots + + return nrow, ncol + + def setup(self, id, nplots, wintitle, showprofile=True, show=True): + + self.__showprofile = showprofile + self.nplots = nplots + + ncolspan = 1 + colspan = 1 + if showprofile: + ncolspan = 7 + colspan = 6 + self.__nsubplots = 2 + + self.createFigure(id = id, + wintitle = wintitle, + widthplot = self.WIDTH + self.WIDTHPROF, + heightplot = self.HEIGHT + self.HEIGHTPROF, + show=show) + + nrow, ncol = self.getSubplots() + + counter = 0 + for y in range(nrow): + for x in range(ncol): + + if counter >= self.nplots: + break + + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1) + + if showprofile: + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1) + + counter += 1 + + def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True', + xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, + timerange=None, + save=False, figpath='./', lastone=0,figfile=None, ftp=False, wr_period=1, show=True, + server=None, folder=None, username=None, password=None, + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + showProfile : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + zmin : None, + zmax : None + """ + + if channelList == None: + channelIndexList = dataOut.channelIndexList + else: + channelIndexList = [] + for channel in channelList: + if channel not in dataOut.channelList: + raise ValueError, "Channel %d is not in dataOut.channelList" + channelIndexList.append(dataOut.channelList.index(channel)) + + if timerange != None: + self.timerange = timerange + + #tmin = None + #tmax = None + factor = dataOut.normFactor + x = dataOut.getTimeRange() + y = dataOut.getHeiRange() + + z = dataOut.data_spc[channelIndexList,:,:]/factor + z = numpy.where(numpy.isfinite(z), z, numpy.NAN) + avg = numpy.average(z, axis=1) + + avgdB = 10.*numpy.log10(avg) + + +# thisDatetime = dataOut.datatime + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " RTI" #: %s" %(thisDatetime.strftime("%d-%b-%Y")) + xlabel = "" + ylabel = "Range (Km)" + + if not self.isConfig: + + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle=wintitle, + showprofile=showprofile, + show=show) + + self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange) + +# if timerange != None: +# self.timerange = timerange +# self.xmin, self.tmax = self.getTimeLim(x, xmin, xmax, timerange) + + + + if ymin == None: ymin = numpy.nanmin(y) + if ymax == None: ymax = numpy.nanmax(y) + if zmin == None: zmin = numpy.nanmin(avgdB)*0.9 + if zmax == None: zmax = numpy.nanmax(avgdB)*0.9 + + self.FTP_WEI = ftp_wei + self.EXP_CODE = exp_code + self.SUB_EXP_CODE = sub_exp_code + self.PLOT_POS = plot_pos + + self.name = thisDatetime.strftime("%Y%m%d_%H%M%S") + self.isConfig = True + + + self.setWinTitle(title) + + if ((self.xmax - x[1]) < (x[1]-x[0])): + x[1] = self.xmax + + for i in range(self.nplots): + title = "Channel %d: %s" %(dataOut.channelList[i]+1, thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[i*self.__nsubplots] + zdB = avgdB[i].reshape((1,-1)) + axes.pcolorbuffer(x, y, zdB, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%") + + if self.__showprofile: + axes = self.axesList[i*self.__nsubplots +1] + axes.pline(avgdB[i], y, + xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax, + xlabel='dB', ylabel='', title='', + ytick_visible=False, + grid='x') + + self.draw() + +# if lastone: +# if dataOut.blocknow >= dataOut.last_block: +# if figfile == None: +# figfile = self.getFilename(name = self.name) +# self.saveFigure(figpath, figfile) +# +# if (save and not(lastone)): +# +# self.counter_imagwr += 1 +# if (self.counter_imagwr==wr_period): +# if figfile == None: +# figfile = self.getFilename(name = self.name) +# self.saveFigure(figpath, figfile) +# +# if ftp: +# #provisionalmente envia archivos en el formato de la web en tiempo real +# name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS) +# path = '%s%03d' %(self.PREFIX, self.id) +# ftp_file = os.path.join(path,'ftp','%s.png'%name) +# self.saveFigure(figpath, ftp_file) +# ftp_filename = os.path.join(figpath,ftp_file) +# self.sendByFTP_Thread(ftp_filename, server, folder, username, password) +# self.counter_imagwr = 0 +# +# self.counter_imagwr = 0 + + #if ((dataOut.utctime-time.timezone) >= self.axesList[0].xmax): + self.saveFigure(figpath, figfile) + if x[1] >= self.axesList[0].xmax: + self.saveFigure(figpath, figfile) + self.isConfig = False + +# if x[1] + (x[1]-x[0]) >= self.axesList[0].xmax: +# +# self.isConfig = False + +# if lastone: +# if figfile == None: +# figfile = self.getFilename(name = self.name) +# self.saveFigure(figpath, figfile) +# +# if ftp: +# #provisionalmente envia archivos en el formato de la web en tiempo real +# name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS) +# path = '%s%03d' %(self.PREFIX, self.id) +# ftp_file = os.path.join(path,'ftp','%s.png'%name) +# self.saveFigure(figpath, ftp_file) +# ftp_filename = os.path.join(figpath,ftp_file) +# self.sendByFTP_Thread(ftp_filename, server, folder, username, password) + +class CoherenceMap(Figure): + isConfig = None + __nsubplots = None + + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'cmap' + + def __init__(self): + self.timerange = 2*60*60 + self.isConfig = False + self.__nsubplots = 1 + + self.WIDTH = 800 + self.HEIGHT = 150 + self.WIDTHPROF = 120 + self.HEIGHTPROF = 0 + self.counter_imagwr = 0 + + self.PLOT_CODE = 3 + self.FTP_WEI = None + self.EXP_CODE = None + self.SUB_EXP_CODE = None + self.PLOT_POS = None + self.counter_imagwr = 0 + + self.xmin = None + self.xmax = None + + def getSubplots(self): + ncol = 1 + nrow = self.nplots*2 + + return nrow, ncol + + def setup(self, id, nplots, wintitle, showprofile=True, show=True): + self.__showprofile = showprofile + self.nplots = nplots + + ncolspan = 1 + colspan = 1 + if showprofile: + ncolspan = 7 + colspan = 6 + self.__nsubplots = 2 + + self.createFigure(id = id, + wintitle = wintitle, + widthplot = self.WIDTH + self.WIDTHPROF, + heightplot = self.HEIGHT + self.HEIGHTPROF, + show=True) + + nrow, ncol = self.getSubplots() + + for y in range(nrow): + for x in range(ncol): + + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1) + + if showprofile: + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1) + + def run(self, dataOut, id, wintitle="", pairsList=None, showprofile='True', + xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, + timerange=None, + save=False, figpath='./', figfile=None, ftp=False, wr_period=1, + coherence_cmap='jet', phase_cmap='RdBu_r', show=True, + server=None, folder=None, username=None, password=None, + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): + + if pairsList == None: + pairsIndexList = dataOut.pairsIndexList + else: + pairsIndexList = [] + for pair in pairsList: + if pair not in dataOut.pairsList: + raise ValueError, "Pair %s is not in dataOut.pairsList" %(pair) + pairsIndexList.append(dataOut.pairsList.index(pair)) + + if timerange != None: + self.timerange = timerange + + if pairsIndexList == []: + return + + if len(pairsIndexList) > 4: + pairsIndexList = pairsIndexList[0:4] + +# tmin = None +# tmax = None + x = dataOut.getTimeRange() + y = dataOut.getHeiRange() + + #thisDatetime = dataOut.datatime + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " CoherenceMap" #: %s" %(thisDatetime.strftime("%d-%b-%Y")) + xlabel = "" + ylabel = "Range (Km)" + + if not self.isConfig: + nplots = len(pairsIndexList) + self.setup(id=id, + nplots=nplots, + wintitle=wintitle, + showprofile=showprofile, + show=show) + + #tmin, tmax = self.getTimeLim(x, xmin, xmax) + + self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange) + + if ymin == None: ymin = numpy.nanmin(y) + if ymax == None: ymax = numpy.nanmax(y) + if zmin == None: zmin = 0. + if zmax == None: zmax = 1. + + self.FTP_WEI = ftp_wei + self.EXP_CODE = exp_code + self.SUB_EXP_CODE = sub_exp_code + self.PLOT_POS = plot_pos + + self.name = thisDatetime.strftime("%Y%m%d_%H%M%S") + + self.isConfig = True + + self.setWinTitle(title) + + if ((self.xmax - x[1]) < (x[1]-x[0])): + x[1] = self.xmax + + for i in range(self.nplots): + + pair = dataOut.pairsList[pairsIndexList[i]] +# coherenceComplex = dataOut.data_cspc[pairsIndexList[i],:,:]/numpy.sqrt(dataOut.data_spc[pair[0],:,:]*dataOut.data_spc[pair[1],:,:]) +# avgcoherenceComplex = numpy.average(coherenceComplex, axis=0) +# coherence = numpy.abs(avgcoherenceComplex) + +## coherence = numpy.abs(coherenceComplex) +## avg = numpy.average(coherence, axis=0) + + ccf = numpy.average(dataOut.data_cspc[pairsIndexList[i],:,:],axis=0) + powa = numpy.average(dataOut.data_spc[pair[0],:,:],axis=0) + powb = numpy.average(dataOut.data_spc[pair[1],:,:],axis=0) + + + avgcoherenceComplex = ccf/numpy.sqrt(powa*powb) + coherence = numpy.abs(avgcoherenceComplex) + + z = coherence.reshape((1,-1)) + + counter = 0 + + title = "Coherence %d%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + axes = self.axesList[i*self.__nsubplots*2] + axes.pcolorbuffer(x, y, z, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', colormap=coherence_cmap, cbsize="1%") + + if self.__showprofile: + counter += 1 + axes = self.axesList[i*self.__nsubplots*2 + counter] + axes.pline(coherence, y, + xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax, + xlabel='', ylabel='', title='', ticksize=7, + ytick_visible=False, nxticks=5, + grid='x') + + counter += 1 +# phase = numpy.arctan(-1*coherenceComplex.imag/coherenceComplex.real)*180/numpy.pi + phase = numpy.arctan2(avgcoherenceComplex.imag, avgcoherenceComplex.real)*180/numpy.pi +# avg = numpy.average(phase, axis=0) + z = phase.reshape((1,-1)) + + title = "Phase %d%d: %s" %(pair[0], pair[1], thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + axes = self.axesList[i*self.__nsubplots*2 + counter] + axes.pcolorbuffer(x, y, z, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=-180, zmax=180, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', colormap=phase_cmap, cbsize="1%") + + if self.__showprofile: + counter += 1 + axes = self.axesList[i*self.__nsubplots*2 + counter] + axes.pline(phase, y, + xmin=-180, xmax=180, ymin=ymin, ymax=ymax, + xlabel='', ylabel='', title='', ticksize=7, + ytick_visible=False, nxticks=4, + grid='x') + + self.draw() + + if x[1] >= self.axesList[0].xmax: + self.saveFigure(figpath, figfile) + self.isConfig = False + +# if save: +# +# self.counter_imagwr += 1 +# if (self.counter_imagwr==wr_period): +# if figfile == None: +# figfile = self.getFilename(name = self.name) +# self.saveFigure(figpath, figfile) +# +# if ftp: +# #provisionalmente envia archivos en el formato de la web en tiempo real +# name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS) +# path = '%s%03d' %(self.PREFIX, self.id) +# ftp_file = os.path.join(path,'ftp','%s.png'%name) +# self.saveFigure(figpath, ftp_file) +# ftp_filename = os.path.join(figpath,ftp_file) +# self.sendByFTP_Thread(ftp_filename, server, folder, username, password) +# self.counter_imagwr = 0 +# +# self.counter_imagwr = 0 +# +# +# if x[1] + (x[1]-x[0]) >= self.axesList[0].xmax: +# self.isConfig = False + +class PowerProfile(Figure): + isConfig = None + __nsubplots = None + + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'spcprofile' + + def __init__(self): + self.isConfig = False + self.__nsubplots = 1 + + self.WIDTH = 300 + self.HEIGHT = 500 + self.counter_imagwr = 0 + + def getSubplots(self): + ncol = 1 + nrow = 1 + + return nrow, ncol + + def setup(self, id, nplots, wintitle, show): + + self.nplots = nplots + + ncolspan = 1 + colspan = 1 + + self.createFigure(id = id, + wintitle = wintitle, + widthplot = self.WIDTH, + heightplot = self.HEIGHT, + show=show) + + nrow, ncol = self.getSubplots() + + counter = 0 + for y in range(nrow): + for x in range(ncol): + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1) + + def run(self, dataOut, id, wintitle="", channelList=None, + xmin=None, xmax=None, ymin=None, ymax=None, + save=False, figpath='./', figfile=None, show=True, wr_period=1, + server=None, folder=None, username=None, password=None,): + + if dataOut.flagNoData: + return None + + if channelList == None: + channelIndexList = dataOut.channelIndexList + channelList = dataOut.channelList + else: + channelIndexList = [] + for channel in channelList: + if channel not in dataOut.channelList: + raise ValueError, "Channel %d is not in dataOut.channelList" + channelIndexList.append(dataOut.channelList.index(channel)) + + try: + factor = dataOut.normFactor + except: + factor = 1 + + y = dataOut.getHeiRange() + + #for voltage + if dataOut.type == 'Voltage': + x = dataOut.data[channelIndexList,:] * numpy.conjugate(dataOut.data[channelIndexList,:]) + x = x.real + x = numpy.where(numpy.isfinite(x), x, numpy.NAN) + + #for spectra + if dataOut.type == 'Spectra': + x = dataOut.data_spc[channelIndexList,:,:]/factor + x = numpy.where(numpy.isfinite(x), x, numpy.NAN) + x = numpy.average(x, axis=1) + + + xdB = 10*numpy.log10(x) + + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " Power Profile %s" %(thisDatetime.strftime("%d-%b-%Y")) + xlabel = "dB" + ylabel = "Range (Km)" + + if not self.isConfig: + + nplots = 1 + + self.setup(id=id, + nplots=nplots, + wintitle=wintitle, + show=show) + + if ymin == None: ymin = numpy.nanmin(y) + if ymax == None: ymax = numpy.nanmax(y) + if xmin == None: xmin = numpy.nanmin(xdB)*0.9 + if xmax == None: xmax = numpy.nanmax(xdB)*0.9 + + self.__isConfig = True + + self.setWinTitle(title) + + title = "Power Profile: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + axes = self.axesList[0] + + legendlabels = ["channel %d"%x for x in channelList] + axes.pmultiline(xdB, y, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, + ytick_visible=True, nxticks=5, + grid='x') + + self.draw() + + if save: + date = thisDatetime.strftime("%Y%m%d") + if figfile == None: + figfile = self.getFilename(name = date) + + self.saveFigure(figpath, figfile) + + self.counter_imagwr += 1 + if (ftp and (self.counter_imagwr==wr_period)): + ftp_filename = os.path.join(figpath,figfile) + self.sendByFTP_Thread(ftp_filename, server, folder, username, password) + self.counter_imagwr = 0 + +class Noise(Figure): + + isConfig = None + __nsubplots = None + + PREFIX = 'noise' + + def __init__(self): + + self.timerange = 24*60*60 + self.isConfig = False + self.__nsubplots = 1 + self.counter_imagwr = 0 + self.WIDTH = 600 + self.HEIGHT = 300 + self.WIDTHPROF = 120 + self.HEIGHTPROF = 0 + self.xdata = None + self.ydata = None + + self.PLOT_CODE = 77 + self.FTP_WEI = None + self.EXP_CODE = None + self.SUB_EXP_CODE = None + self.PLOT_POS = None + + def getSubplots(self): + + ncol = 1 + nrow = 1 + + return nrow, ncol + + def openfile(self, filename): + f = open(filename,'w+') + f.write('\n\n') + f.write('JICAMARCA RADIO OBSERVATORY - Noise \n') + f.write('DD MM YYYY HH MM SS Channel0 Channel1 Channel2 Channel3\n\n' ) + f.close() + + def save_data(self, filename_phase, data, data_datetime): + f=open(filename_phase,'a') + timetuple_data = data_datetime.timetuple() + day = str(timetuple_data.tm_mday) + month = str(timetuple_data.tm_mon) + year = str(timetuple_data.tm_year) + hour = str(timetuple_data.tm_hour) + minute = str(timetuple_data.tm_min) + second = str(timetuple_data.tm_sec) + f.write(day+' '+month+' '+year+' '+hour+' '+minute+' '+second+' '+str(data[0])+' '+str(data[1])+' '+str(data[2])+' '+str(data[3])+'\n') + f.close() + + + def setup(self, id, nplots, wintitle, showprofile=True, show=True): + + self.__showprofile = showprofile + self.nplots = nplots + + ncolspan = 7 + colspan = 6 + self.__nsubplots = 2 + + self.createFigure(id = id, + wintitle = wintitle, + widthplot = self.WIDTH+self.WIDTHPROF, + heightplot = self.HEIGHT+self.HEIGHTPROF, + show=show) + + nrow, ncol = self.getSubplots() + + self.addAxes(nrow, ncol*ncolspan, 0, 0, colspan, 1) + + + def run(self, dataOut, id, wintitle="", channelList=None, showprofile='True', + xmin=None, xmax=None, ymin=None, ymax=None, + timerange=None, + save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1, + server=None, folder=None, username=None, password=None, + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): + + if channelList == None: + channelIndexList = dataOut.channelIndexList + channelList = dataOut.channelList + else: + channelIndexList = [] + for channel in channelList: + if channel not in dataOut.channelList: + raise ValueError, "Channel %d is not in dataOut.channelList" + channelIndexList.append(dataOut.channelList.index(channel)) + + if timerange != None: + self.timerange = timerange + + tmin = None + tmax = None + x = dataOut.getTimeRange() + y = dataOut.getHeiRange() + factor = dataOut.normFactor + noise = dataOut.noise()/factor + noisedB = 10*numpy.log10(noise) + + #thisDatetime = dataOut.datatime + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " Noise" # : %s" %(thisDatetime.strftime("%d-%b-%Y")) + xlabel = "" + ylabel = "Intensity (dB)" + + if not self.isConfig: + + nplots = 1 + + self.setup(id=id, + nplots=nplots, + wintitle=wintitle, + showprofile=showprofile, + show=show) + + tmin, tmax = self.getTimeLim(x, xmin, xmax) + if ymin == None: ymin = numpy.nanmin(noisedB) - 10.0 + if ymax == None: ymax = numpy.nanmax(noisedB) + 10.0 + + self.FTP_WEI = ftp_wei + self.EXP_CODE = exp_code + self.SUB_EXP_CODE = sub_exp_code + self.PLOT_POS = plot_pos + + + self.name = thisDatetime.strftime("%Y%m%d_%H%M%S") + self.isConfig = True + + self.xdata = numpy.array([]) + self.ydata = numpy.array([]) + + #open file beacon phase + path = '%s%03d' %(self.PREFIX, self.id) + noise_file = os.path.join(path,'%s.txt'%self.name) + self.filename_noise = os.path.join(figpath,noise_file) + self.openfile(self.filename_noise) + + + #store data beacon phase + self.save_data(self.filename_noise, noisedB, thisDatetime) + + + self.setWinTitle(title) + + + title = "Noise %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + + legendlabels = ["channel %d"%(idchannel+1) for idchannel in channelList] + axes = self.axesList[0] + + self.xdata = numpy.hstack((self.xdata, x[0:1])) + + if len(self.ydata)==0: + self.ydata = noisedB[channelIndexList].reshape(-1,1) + else: + self.ydata = numpy.hstack((self.ydata, noisedB[channelIndexList].reshape(-1,1))) + + + axes.pmultilineyaxis(x=self.xdata, y=self.ydata, + xmin=tmin, xmax=tmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title, legendlabels=legendlabels, marker='x', markersize=8, linestyle="solid", + XAxisAsTime=True, grid='both' + ) + + self.draw() + +# if save: +# +# if figfile == None: +# figfile = self.getFilename(name = self.name) +# +# self.saveFigure(figpath, figfile) + + if save: + + self.counter_imagwr += 1 + if (self.counter_imagwr==wr_period): + if figfile == None: + figfile = self.getFilename(name = self.name) + self.saveFigure(figpath, figfile) + + if ftp: + #provisionalmente envia archivos en el formato de la web en tiempo real + name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS) + path = '%s%03d' %(self.PREFIX, self.id) + ftp_file = os.path.join(path,'ftp','%s.png'%name) + self.saveFigure(figpath, ftp_file) + ftp_filename = os.path.join(figpath,ftp_file) + self.sendByFTP_Thread(ftp_filename, server, folder, username, password) + self.counter_imagwr = 0 + + self.counter_imagwr = 0 + + if x[1] + (x[1]-x[0]) >= self.axesList[0].xmax: + self.isConfig = False + del self.xdata + del self.ydata diff --git a/schainpy/model/graphics/jroplot_voltage.py b/schainpy/model/graphics/jroplot_voltage.py new file mode 100644 index 0000000..4860eb2 --- /dev/null +++ b/schainpy/model/graphics/jroplot_voltage.py @@ -0,0 +1,186 @@ +''' +@author: Daniel Suarez +''' + +import datetime +import numpy + +from figure import Figure + +class Scope(Figure): + + isConfig = None + + def __init__(self): + + self.isConfig = False + self.WIDTH = 300 + self.HEIGHT = 200 + self.counter_imagwr = 0 + + def getSubplots(self): + + nrow = self.nplots + ncol = 3 + return nrow, ncol + + def setup(self, id, nplots, wintitle, show): + + self.nplots = nplots + + self.createFigure(id=id, + wintitle=wintitle, + show=show) + + nrow,ncol = self.getSubplots() + colspan = 3 + rowspan = 1 + + for i in range(nplots): + self.addAxes(nrow, ncol, i, 0, colspan, rowspan) + + def plot_iq(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax): + yreal = y[channelIndexList,:].real + yimag = y[channelIndexList,:].imag + + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Range (Km)" + ylabel = "Intensity - IQ" + + if not self.isConfig: + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle='', + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = min(numpy.nanmin(yreal),numpy.nanmin(yimag)) + if ymax == None: ymax = max(numpy.nanmax(yreal),numpy.nanmax(yimag)) + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(len(self.axesList)): + title = "Channel %d" %(i) + axes = self.axesList[i] + + axes.pline(x, yreal[i,:], + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title) + + axes.addpline(x, yimag[i,:], idline=1, color="red", linestyle="solid", lw=2) + + def plot_power(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax): + y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:]) + yreal = y.real + + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Range (Km)" + ylabel = "Intensity" + + if not self.isConfig: + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle='', + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = numpy.nanmin(yreal) + if ymax == None: ymax = numpy.nanmax(yreal) + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(len(self.axesList)): + title = "Channel %d" %(i) + axes = self.axesList[i] + ychannel = yreal[i,:] + axes.pline(x, ychannel, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title) + + + def run(self, dataOut, id, wintitle="", channelList=None, + xmin=None, xmax=None, ymin=None, ymax=None, save=False, + figpath='./', figfile=None, show=True, wr_period=1, + server=None, folder=None, username=None, password=None, type='power'): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + """ + if dataOut.flagNoData: + return None + + if channelList == None: + channelIndexList = dataOut.channelIndexList + else: + channelIndexList = [] + for channel in channelList: + if channel not in dataOut.channelList: + raise ValueError, "Channel %d is not in dataOut.channelList" + channelIndexList.append(dataOut.channelList.index(channel)) + + x = dataOut.heightList + y = dataOut.data[channelIndexList,:] * numpy.conjugate(dataOut.data[channelIndexList,:]) + y = y.real + + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + + if type == "power": + self.plot_power(dataOut.heightList, + dataOut.data, + id, + channelIndexList, + thisDatetime, + wintitle, + show, + xmin, + xmax, + ymin, + ymax) + + if type == "iq": + self.plot_iq(dataOut.heightList, + dataOut.data, + id, + channelIndexList, + thisDatetime, + wintitle, + show, + xmin, + xmax, + ymin, + ymax) + + + self.draw() + + if save: + date = thisDatetime.strftime("%Y%m%d_%H%M%S") + if figfile == None: + figfile = self.getFilename(name = date) + + self.saveFigure(figpath, figfile) + + self.counter_imagwr += 1 + if (ftp and (self.counter_imagwr==wr_period)): + ftp_filename = os.path.join(figpath,figfile) + self.sendByFTP_Thread(ftp_filename, server, folder, username, password) + self.counter_imagwr = 0 diff --git a/schainpy/model/graphics/mpldriver.py b/schainpy/model/graphics/mpldriver.py index 8c9f971..691b4cb 100644 --- a/schainpy/model/graphics/mpldriver.py +++ b/schainpy/model/graphics/mpldriver.py @@ -7,7 +7,7 @@ if 'linux' in sys.platform: matplotlib.use("TKAgg") if 'darwin' in sys.platform: - matplotlib.use("GTKAgg") + matplotlib.use("TKAgg") import matplotlib.pyplot diff --git a/schainpy/model/io/__init__.py b/schainpy/model/io/__init__.py new file mode 100644 index 0000000..e69de29 --- /dev/null +++ b/schainpy/model/io/__init__.py diff --git a/schainpy/model/io/jroIO_base.py b/schainpy/model/io/jroIO_base.py new file mode 100644 index 0000000..fd15626 --- /dev/null +++ b/schainpy/model/io/jroIO_base.py @@ -0,0 +1,1334 @@ +''' + +''' +import os +import sys +import glob +import time +import numpy +import fnmatch +import time, datetime +import h5py +import traceback + +#try: +# import pyfits +#except: +# print "pyfits module has not been imported, it should be installed to save files in fits format" + +#from jrodata import * +#from jroheaderIO import * +#from jroprocessing import * + +#import re +#from xml.etree.ElementTree import Element, SubElement, ElementTree + + +LOCALTIME = True #-18000 + +from model.data.jroheaderIO import PROCFLAG, BasicHeader, SystemHeader, RadarControllerHeader, ProcessingHeader + +def isNumber(str): + """ + Chequea si el conjunto de caracteres que componen un string puede ser convertidos a un numero. + + Excepciones: + Si un determinado string no puede ser convertido a numero + Input: + str, string al cual se le analiza para determinar si convertible a un numero o no + + Return: + True : si el string es uno numerico + False : no es un string numerico + """ + try: + float( str ) + return True + except: + return False + +def isThisFileinRange(filename, startUTSeconds, endUTSeconds): + """ + Esta funcion determina si un archivo de datos se encuentra o no dentro del rango de fecha especificado. + + Inputs: + filename : nombre completo del archivo de datos en formato Jicamarca (.r) + + startUTSeconds : fecha inicial del rango seleccionado. La fecha esta dada en + segundos contados desde 01/01/1970. + endUTSeconds : fecha final del rango seleccionado. La fecha esta dada en + segundos contados desde 01/01/1970. + + Return: + Boolean : Retorna True si el archivo de datos contiene datos en el rango de + fecha especificado, de lo contrario retorna False. + + Excepciones: + Si el archivo no existe o no puede ser abierto + Si la cabecera no puede ser leida. + + """ + basicHeaderObj = BasicHeader(LOCALTIME) + + try: + fp = open(filename,'rb') + except IOError: + traceback.print_exc() + raise IOError, "The file %s can't be opened" %(filename) + + sts = basicHeaderObj.read(fp) + fp.close() + + if not(sts): + print "Skipping the file %s because it has not a valid header" %(filename) + return 0 + + if not ((startUTSeconds <= basicHeaderObj.utc) and (endUTSeconds > basicHeaderObj.utc)): + return 0 + + return 1 + +def isFileinThisTime(filename, startTime, endTime): + """ + Retorna 1 si el archivo de datos se encuentra dentro del rango de horas especificado. + + Inputs: + filename : nombre completo del archivo de datos en formato Jicamarca (.r) + + startTime : tiempo inicial del rango seleccionado en formato datetime.time + + endTime : tiempo final del rango seleccionado en formato datetime.time + + Return: + Boolean : Retorna True si el archivo de datos contiene datos en el rango de + fecha especificado, de lo contrario retorna False. + + Excepciones: + Si el archivo no existe o no puede ser abierto + Si la cabecera no puede ser leida. + + """ + + + try: + fp = open(filename,'rb') + except IOError: + traceback.print_exc() + raise IOError, "The file %s can't be opened" %(filename) + + basicHeaderObj = BasicHeader(LOCALTIME) + sts = basicHeaderObj.read(fp) + fp.close() + + thisDatetime = basicHeaderObj.datatime + thisTime = thisDatetime.time() + + if not(sts): + print "Skipping the file %s because it has not a valid header" %(filename) + return None + + if not ((startTime <= thisTime) and (endTime > thisTime)): + return None + + return thisDatetime + +def getFileFromSet(path, ext, set): + validFilelist = [] + fileList = os.listdir(path) + + # 0 1234 567 89A BCDE + # H YYYY DDD SSS .ext + + for thisFile in fileList: + try: + year = int(thisFile[1:5]) + doy = int(thisFile[5:8]) + except: + continue + + if (os.path.splitext(thisFile)[-1].lower() != ext.lower()): + continue + + validFilelist.append(thisFile) + + myfile = fnmatch.filter(validFilelist,'*%4.4d%3.3d%3.3d*'%(year,doy,set)) + + if len(myfile)!= 0: + return myfile[0] + else: + filename = '*%4.4d%3.3d%3.3d%s'%(year,doy,set,ext.lower()) + print 'the filename %s does not exist'%filename + print '...going to the last file: ' + + if validFilelist: + validFilelist = sorted( validFilelist, key=str.lower ) + return validFilelist[-1] + + return None + +def getlastFileFromPath(path, ext): + """ + Depura el fileList dejando solo los que cumplan el formato de "PYYYYDDDSSS.ext" + al final de la depuracion devuelve el ultimo file de la lista que quedo. + + Input: + fileList : lista conteniendo todos los files (sin path) que componen una determinada carpeta + ext : extension de los files contenidos en una carpeta + + Return: + El ultimo file de una determinada carpeta, no se considera el path. + """ + validFilelist = [] + fileList = os.listdir(path) + + # 0 1234 567 89A BCDE + # H YYYY DDD SSS .ext + + for thisFile in fileList: + + year = thisFile[1:5] + if not isNumber(year): + continue + + doy = thisFile[5:8] + if not isNumber(doy): + continue + + year = int(year) + doy = int(doy) + + if (os.path.splitext(thisFile)[-1].lower() != ext.lower()): + continue + + validFilelist.append(thisFile) + + if validFilelist: + validFilelist = sorted( validFilelist, key=str.lower ) + return validFilelist[-1] + + return None + +def checkForRealPath(path, foldercounter, year, doy, set, ext): + """ + Por ser Linux Case Sensitive entonces checkForRealPath encuentra el nombre correcto de un path, + Prueba por varias combinaciones de nombres entre mayusculas y minusculas para determinar + el path exacto de un determinado file. + + Example : + nombre correcto del file es .../.../D2009307/P2009307367.ext + + Entonces la funcion prueba con las siguientes combinaciones + .../.../y2009307367.ext + .../.../Y2009307367.ext + .../.../x2009307/y2009307367.ext + .../.../x2009307/Y2009307367.ext + .../.../X2009307/y2009307367.ext + .../.../X2009307/Y2009307367.ext + siendo para este caso, la ultima combinacion de letras, identica al file buscado + + Return: + Si encuentra la cobinacion adecuada devuelve el path completo y el nombre del file + caso contrario devuelve None como path y el la ultima combinacion de nombre en mayusculas + para el filename + """ + fullfilename = None + find_flag = False + filename = None + + prefixDirList = [None,'d','D'] + if ext.lower() == ".r": #voltage + prefixFileList = ['d','D'] + elif ext.lower() == ".pdata": #spectra + prefixFileList = ['p','P'] + else: + return None, filename + + #barrido por las combinaciones posibles + for prefixDir in prefixDirList: + thispath = path + if prefixDir != None: + #formo el nombre del directorio xYYYYDDD (x=d o x=D) + if foldercounter == 0: + thispath = os.path.join(path, "%s%04d%03d" % ( prefixDir, year, doy )) + else: + thispath = os.path.join(path, "%s%04d%03d_%02d" % ( prefixDir, year, doy , foldercounter)) + for prefixFile in prefixFileList: #barrido por las dos combinaciones posibles de "D" + filename = "%s%04d%03d%03d%s" % ( prefixFile, year, doy, set, ext ) #formo el nombre del file xYYYYDDDSSS.ext + fullfilename = os.path.join( thispath, filename ) #formo el path completo + + if os.path.exists( fullfilename ): #verifico que exista + find_flag = True + break + if find_flag: + break + + if not(find_flag): + return None, filename + + return fullfilename, filename + +def isDoyFolder(folder): + try: + year = int(folder[1:5]) + except: + return 0 + + try: + doy = int(folder[5:8]) + except: + return 0 + + return 1 + +class JRODataIO: + + c = 3E8 + + isConfig = False + + basicHeaderObj = None + + systemHeaderObj = None + + radarControllerHeaderObj = None + + processingHeaderObj = None + + online = 0 + + dtype = None + + pathList = [] + + filenameList = [] + + filename = None + + ext = None + + flagIsNewFile = 1 + + flagTimeBlock = 0 + + flagIsNewBlock = 0 + + fp = None + + firstHeaderSize = 0 + + basicHeaderSize = 24 + + versionFile = 1103 + + fileSize = None + +# ippSeconds = None + + fileSizeByHeader = None + + fileIndex = None + + profileIndex = None + + blockIndex = None + + nTotalBlocks = None + + maxTimeStep = 30 + + lastUTTime = None + + datablock = None + + dataOut = None + + blocksize = None + + def __init__(self): + + raise ValueError, "Not implemented" + + def run(self): + + raise ValueError, "Not implemented" + +class JRODataReader(JRODataIO): + + nReadBlocks = 0 + + delay = 10 #number of seconds waiting a new file + + nTries = 3 #quantity tries + + nFiles = 3 #number of files for searching + + path = None + + foldercounter = 0 + + flagNoMoreFiles = 0 + + datetimeList = [] + + __isFirstTimeOnline = 1 + + __printInfo = True + + profileIndex = None + + def __init__(self): + + """ + + """ + + raise ValueError, "This method has not been implemented" + + + def createObjByDefault(self): + """ + + """ + raise ValueError, "This method has not been implemented" + + def getBlockDimension(self): + + raise ValueError, "No implemented" + + def __searchFilesOffLine(self, + path, + startDate, + endDate, + startTime=datetime.time(0,0,0), + endTime=datetime.time(23,59,59), + set=None, + expLabel='', + ext='.r', + walk=True): + + pathList = [] + + if not walk: + #pathList.append(path) + multi_path = path.split(',') + for single_path in multi_path: + pathList.append(single_path) + + else: + #dirList = [] + multi_path = path.split(',') + for single_path in multi_path: + dirList = [] + for thisPath in os.listdir(single_path): + if not os.path.isdir(os.path.join(single_path,thisPath)): + continue + if not isDoyFolder(thisPath): + continue + + dirList.append(thisPath) + + if not(dirList): + return None, None + + thisDate = startDate + + while(thisDate <= endDate): + year = thisDate.timetuple().tm_year + doy = thisDate.timetuple().tm_yday + + matchlist = fnmatch.filter(dirList, '?' + '%4.4d%3.3d' % (year,doy) + '*') + if len(matchlist) == 0: + thisDate += datetime.timedelta(1) + continue + for match in matchlist: + pathList.append(os.path.join(single_path,match,expLabel)) + + thisDate += datetime.timedelta(1) + + if pathList == []: + print "Any folder was found for the date range: %s-%s" %(startDate, endDate) + return None, None + + print "%d folder(s) was(were) found for the date range: %s - %s" %(len(pathList), startDate, endDate) + + filenameList = [] + datetimeList = [] + pathDict = {} + filenameList_to_sort = [] + + for i in range(len(pathList)): + + thisPath = pathList[i] + + fileList = glob.glob1(thisPath, "*%s" %ext) + fileList.sort() + pathDict.setdefault(fileList[0]) + pathDict[fileList[0]] = i + filenameList_to_sort.append(fileList[0]) + + filenameList_to_sort.sort() + + for file in filenameList_to_sort: + thisPath = pathList[pathDict[file]] + + fileList = glob.glob1(thisPath, "*%s" %ext) + fileList.sort() + + for file in fileList: + + filename = os.path.join(thisPath,file) + thisDatetime = isFileinThisTime(filename, startTime, endTime) + + if not(thisDatetime): + continue + + filenameList.append(filename) + datetimeList.append(thisDatetime) + + if not(filenameList): + print "Any file was found for the time range %s - %s" %(startTime, endTime) + return None, None + + print "%d file(s) was(were) found for the time range: %s - %s" %(len(filenameList), startTime, endTime) + print + + for i in range(len(filenameList)): + print "%s -> [%s]" %(filenameList[i], datetimeList[i].ctime()) + + self.filenameList = filenameList + self.datetimeList = datetimeList + + return pathList, filenameList + + def __searchFilesOnLine(self, path, expLabel = "", ext = None, walk=True, set=None): + + """ + Busca el ultimo archivo de la ultima carpeta (determinada o no por startDateTime) y + devuelve el archivo encontrado ademas de otros datos. + + Input: + path : carpeta donde estan contenidos los files que contiene data + + expLabel : Nombre del subexperimento (subfolder) + + ext : extension de los files + + walk : Si es habilitado no realiza busquedas dentro de los ubdirectorios (doypath) + + Return: + directory : eL directorio donde esta el file encontrado + filename : el ultimo file de una determinada carpeta + year : el anho + doy : el numero de dia del anho + set : el set del archivo + + + """ + dirList = [] + + if not walk: + fullpath = path + foldercounter = 0 + else: + #Filtra solo los directorios + for thisPath in os.listdir(path): + if not os.path.isdir(os.path.join(path,thisPath)): + continue + if not isDoyFolder(thisPath): + continue + + dirList.append(thisPath) + + if not(dirList): + return None, None, None, None, None, None + + dirList = sorted( dirList, key=str.lower ) + + doypath = dirList[-1] + foldercounter = int(doypath.split('_')[1]) if len(doypath.split('_'))>1 else 0 + fullpath = os.path.join(path, doypath, expLabel) + + + print "%s folder was found: " %(fullpath ) + + if set == None: + filename = getlastFileFromPath(fullpath, ext) + else: + filename = getFileFromSet(fullpath, ext, set) + + if not(filename): + return None, None, None, None, None, None + + print "%s file was found" %(filename) + + if not(self.__verifyFile(os.path.join(fullpath, filename))): + return None, None, None, None, None, None + + year = int( filename[1:5] ) + doy = int( filename[5:8] ) + set = int( filename[8:11] ) + + return fullpath, foldercounter, filename, year, doy, set + + 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] + + if not(self.__verifyFile(filename)): + continue + + fileSize = os.path.getsize(filename) + fp = open(filename,'rb') + break + + self.flagIsNewFile = 1 + self.fileIndex = idFile + self.filename = filename + self.fileSize = fileSize + self.fp = fp + + print "Setting the file: %s"%self.filename + + return 1 + + def __setNextFileOnline(self): + """ + Busca el siguiente file que tenga suficiente data para ser leida, dentro de un folder especifico, si + no encuentra un file valido espera un tiempo determinado y luego busca en los posibles n files + siguientes. + + Affected: + self.flagIsNewFile + self.filename + self.fileSize + self.fp + self.set + self.flagNoMoreFiles + + Return: + 0 : si luego de una busqueda del siguiente file valido este no pudo ser encontrado + 1 : si el file fue abierto con exito y esta listo a ser leido + + Excepciones: + Si un determinado file no puede ser abierto + """ + nFiles = 0 + fileOk_flag = False + firstTime_flag = True + + self.set += 1 + + if self.set > 999: + self.set = 0 + self.foldercounter += 1 + + #busca el 1er file disponible + fullfilename, filename = checkForRealPath( self.path, self.foldercounter, self.year, self.doy, self.set, self.ext ) + if fullfilename: + if self.__verifyFile(fullfilename, False): + fileOk_flag = True + + #si no encuentra un file entonces espera y vuelve a buscar + if not(fileOk_flag): + for nFiles in range(self.nFiles+1): #busco en los siguientes self.nFiles+1 files posibles + + if firstTime_flag: #si es la 1era vez entonces hace el for self.nTries veces + tries = self.nTries + else: + tries = 1 #si no es la 1era vez entonces solo lo hace una vez + + for nTries in range( tries ): + if firstTime_flag: + print "\tWaiting %0.2f sec for the file \"%s\" , try %03d ..." % ( self.delay, filename, nTries+1 ) + time.sleep( self.delay ) + else: + print "\tSearching next \"%s%04d%03d%03d%s\" file ..." % (self.optchar, self.year, self.doy, self.set, self.ext) + + fullfilename, filename = checkForRealPath( self.path, self.foldercounter, self.year, self.doy, self.set, self.ext ) + if fullfilename: + if self.__verifyFile(fullfilename): + fileOk_flag = True + break + + if fileOk_flag: + break + + firstTime_flag = False + + print "\tSkipping the file \"%s\" due to this file doesn't exist" % filename + self.set += 1 + + if nFiles == (self.nFiles-1): #si no encuentro el file buscado cambio de carpeta y busco en la siguiente carpeta + self.set = 0 + self.doy += 1 + self.foldercounter = 0 + + if fileOk_flag: + self.fileSize = os.path.getsize( fullfilename ) + self.filename = fullfilename + self.flagIsNewFile = 1 + if self.fp != None: self.fp.close() + self.fp = open(fullfilename, 'rb') + self.flagNoMoreFiles = 0 + print 'Setting the file: %s' % fullfilename + else: + self.fileSize = 0 + self.filename = None + self.flagIsNewFile = 0 + self.fp = None + self.flagNoMoreFiles = 1 + print 'No more Files' + + return fileOk_flag + + def setNextFile(self): + if self.fp != None: + self.fp.close() + + if self.online: + newFile = self.__setNextFileOnline() + else: + newFile = self.__setNextFileOffline() + + if not(newFile): + return 0 + + self.__readFirstHeader() + self.nReadBlocks = 0 + return 1 + + def __waitNewBlock(self): + """ + Return 1 si se encontro un nuevo bloque de datos, 0 de otra forma. + + Si el modo de lectura es OffLine siempre retorn 0 + """ + if not self.online: + return 0 + + if (self.nReadBlocks >= self.processingHeaderObj.dataBlocksPerFile): + return 0 + + currentPointer = self.fp.tell() + + neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize + + for nTries in range( self.nTries ): + + self.fp.close() + self.fp = open( self.filename, 'rb' ) + self.fp.seek( currentPointer ) + + self.fileSize = os.path.getsize( self.filename ) + currentSize = self.fileSize - currentPointer + + if ( currentSize >= neededSize ): + self.basicHeaderObj.read(self.fp) + return 1 + + if self.fileSize == self.fileSizeByHeader: +# self.flagEoF = True + return 0 + + print "\tWaiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries+1) + time.sleep( self.delay ) + + + return 0 + + def waitDataBlock(self,pointer_location): + + currentPointer = pointer_location + + neededSize = self.processingHeaderObj.blockSize #+ self.basicHeaderSize + + for nTries in range( self.nTries ): + self.fp.close() + self.fp = open( self.filename, 'rb' ) + self.fp.seek( currentPointer ) + + self.fileSize = os.path.getsize( self.filename ) + currentSize = self.fileSize - currentPointer + + if ( currentSize >= neededSize ): + return 1 + + print "\tWaiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries+1) + time.sleep( self.delay ) + + return 0 + + def __jumpToLastBlock(self): + + if not(self.__isFirstTimeOnline): + return + + csize = self.fileSize - self.fp.tell() + blocksize = self.processingHeaderObj.blockSize + + #salta el primer bloque de datos + if csize > self.processingHeaderObj.blockSize: + self.fp.seek(self.fp.tell() + blocksize) + else: + return + + csize = self.fileSize - self.fp.tell() + neededsize = self.processingHeaderObj.blockSize + self.basicHeaderSize + while True: + + if self.fp.tell() 0: +# self.fp.seek(self.fp.tell() + factor*neededsize) + + self.flagIsNewFile = 0 + self.__isFirstTimeOnline = 0 + + def __setNewBlock(self): + + if self.fp == None: + return 0 + + if self.online: + self.__jumpToLastBlock() + + if self.flagIsNewFile: + return 1 + + self.lastUTTime = self.basicHeaderObj.utc + currentSize = self.fileSize - self.fp.tell() + neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize + + if (currentSize >= neededSize): + self.basicHeaderObj.read(self.fp) + return 1 + + if self.__waitNewBlock(): + return 1 + + if not(self.setNextFile()): + return 0 + + deltaTime = self.basicHeaderObj.utc - self.lastUTTime # + + self.flagTimeBlock = 0 + + if deltaTime > self.maxTimeStep: + self.flagTimeBlock = 1 + + return 1 + + def readNextBlock(self): + if not(self.__setNewBlock()): + return 0 + + if not(self.readBlock()): + return 0 + + return 1 + + def __readFirstHeader(self): + + self.basicHeaderObj.read(self.fp) + self.systemHeaderObj.read(self.fp) + self.radarControllerHeaderObj.read(self.fp) + self.processingHeaderObj.read(self.fp) + + self.firstHeaderSize = self.basicHeaderObj.size + + datatype = int(numpy.log2((self.processingHeaderObj.processFlags & PROCFLAG.DATATYPE_MASK))-numpy.log2(PROCFLAG.DATATYPE_CHAR)) + if datatype == 0: + datatype_str = numpy.dtype([('real',' %s" %(self.basicHeaderObj.dataBlock, self.nTotalBlocks, self.dataOut.datatime.ctime()) + self.dataOut.blocknow = self.basicHeaderObj.dataBlock + + def printInfo(self): + + if self.__printInfo == False: + return + + self.basicHeaderObj.printInfo() + self.systemHeaderObj.printInfo() + self.radarControllerHeaderObj.printInfo() + self.processingHeaderObj.printInfo() + + self.__printInfo = False + + + def run(self, **kwargs): + + if not(self.isConfig): + +# self.dataOut = dataOut + self.setup(**kwargs) + self.isConfig = True + + self.getData() + +class JRODataWriter(JRODataIO): + + """ + Esta clase permite escribir datos a archivos procesados (.r o ,pdata). La escritura + de los datos siempre se realiza por bloques. + """ + + blockIndex = 0 + + path = None + + setFile = None + + profilesPerBlock = None + + blocksPerFile = None + + nWriteBlocks = 0 + + def __init__(self, dataOut=None): + raise ValueError, "Not implemented" + + + def hasAllDataInBuffer(self): + raise ValueError, "Not implemented" + + + def setBlockDimension(self): + raise ValueError, "Not implemented" + + + def writeBlock(self): + raise ValueError, "No implemented" + + + def putData(self): + raise ValueError, "No implemented" + + + def setBasicHeader(self): + + self.basicHeaderObj.size = self.basicHeaderSize #bytes + self.basicHeaderObj.version = self.versionFile + self.basicHeaderObj.dataBlock = self.nTotalBlocks + + utc = numpy.floor(self.dataOut.utctime) + milisecond = (self.dataOut.utctime - utc)* 1000.0 + + self.basicHeaderObj.utc = utc + self.basicHeaderObj.miliSecond = milisecond + self.basicHeaderObj.timeZone = self.dataOut.timeZone + self.basicHeaderObj.dstFlag = self.dataOut.dstFlag + self.basicHeaderObj.errorCount = self.dataOut.errorCount + + def setFirstHeader(self): + """ + Obtiene una copia del First Header + + Affected: + + self.basicHeaderObj + self.systemHeaderObj + self.radarControllerHeaderObj + self.processingHeaderObj self. + + Return: + None + """ + + raise ValueError, "No implemented" + + def __writeFirstHeader(self): + """ + Escribe el primer header del file es decir el Basic header y el Long header (SystemHeader, RadarControllerHeader, ProcessingHeader) + + Affected: + __dataType + + Return: + None + """ + +# CALCULAR PARAMETROS + + sizeLongHeader = self.systemHeaderObj.size + self.radarControllerHeaderObj.size + self.processingHeaderObj.size + self.basicHeaderObj.size = self.basicHeaderSize + sizeLongHeader + + self.basicHeaderObj.write(self.fp) + self.systemHeaderObj.write(self.fp) + self.radarControllerHeaderObj.write(self.fp) + self.processingHeaderObj.write(self.fp) + + self.dtype = self.dataOut.dtype + + def __setNewBlock(self): + """ + Si es un nuevo file escribe el First Header caso contrario escribe solo el Basic Header + + Return: + 0 : si no pudo escribir nada + 1 : Si escribio el Basic el First Header + """ + if self.fp == None: + self.setNextFile() + + if self.flagIsNewFile: + return 1 + + if self.blockIndex < self.processingHeaderObj.dataBlocksPerFile: + self.basicHeaderObj.write(self.fp) + return 1 + + if not( self.setNextFile() ): + return 0 + + return 1 + + + def writeNextBlock(self): + """ + Selecciona el bloque siguiente de datos y los escribe en un file + + Return: + 0 : Si no hizo pudo escribir el bloque de datos + 1 : Si no pudo escribir el bloque de datos + """ + if not( self.__setNewBlock() ): + return 0 + + self.writeBlock() + + return 1 + + def setNextFile(self): + """ + Determina el siguiente file que sera escrito + + Affected: + self.filename + self.subfolder + self.fp + self.setFile + self.flagIsNewFile + + Return: + 0 : Si el archivo no puede ser escrito + 1 : Si el archivo esta listo para ser escrito + """ + ext = self.ext + path = self.path + + if self.fp != None: + self.fp.close() + + timeTuple = time.localtime( self.dataOut.utctime) + subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday) + + fullpath = os.path.join( path, subfolder ) + if not( os.path.exists(fullpath) ): + os.mkdir(fullpath) + self.setFile = -1 #inicializo mi contador de seteo + else: + filesList = os.listdir( fullpath ) + if len( filesList ) > 0: + filesList = sorted( filesList, key=str.lower ) + filen = filesList[-1] + # el filename debera tener el siguiente formato + # 0 1234 567 89A BCDE (hex) + # x YYYY DDD SSS .ext + if isNumber( filen[8:11] ): + self.setFile = int( filen[8:11] ) #inicializo mi contador de seteo al seteo del ultimo file + else: + self.setFile = -1 + else: + self.setFile = -1 #inicializo mi contador de seteo + + setFile = self.setFile + setFile += 1 + + file = '%s%4.4d%3.3d%3.3d%s' % (self.optchar, + timeTuple.tm_year, + timeTuple.tm_yday, + setFile, + ext ) + + filename = os.path.join( path, subfolder, file ) + + fp = open( filename,'wb' ) + + self.blockIndex = 0 + + #guardando atributos + self.filename = filename + self.subfolder = subfolder + self.fp = fp + self.setFile = setFile + self.flagIsNewFile = 1 + + self.setFirstHeader() + + print 'Writing the file: %s'%self.filename + + self.__writeFirstHeader() + + return 1 + + def setup(self, dataOut, path, blocksPerFile, profilesPerBlock=64, set=0, ext=None): + """ + Setea el tipo de formato en la cual sera guardada la data y escribe el First Header + + Inputs: + path : el path destino en el cual se escribiran los files a crear + format : formato en el cual sera salvado un file + set : el setebo del file + + Return: + 0 : Si no realizo un buen seteo + 1 : Si realizo un buen seteo + """ + + if ext == None: + ext = self.ext + + ext = ext.lower() + + self.ext = ext + + self.path = path + + self.setFile = set - 1 + + self.blocksPerFile = blocksPerFile + + self.profilesPerBlock = profilesPerBlock + + self.dataOut = dataOut + + if not(self.setNextFile()): + print "There isn't a next file" + return 0 + + self.setBlockDimension() + + return 1 + + def run(self, dataOut, **kwargs): + + if not(self.isConfig): + + self.setup(dataOut, **kwargs) + self.isConfig = True + + self.putData() + diff --git a/schainpy/model/io/jroIO_example.py b/schainpy/model/io/jroIO_example.py new file mode 100644 index 0000000..a26a393 --- /dev/null +++ b/schainpy/model/io/jroIO_example.py @@ -0,0 +1,100 @@ +''' +Created on Jul 3, 2014 + +@author: roj-idl71 +''' + +import os + +from schainpy.model.data.jrodata import Voltage +from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation + +class Reader(ProcessingUnit): + ''' + classdocs + ''' + + def __init__(self): + ''' + Constructor + ''' + + ProcessingUnit.__init__(self) + +# self.dataIn = None +# +# self.isConfig = False + + #Is really necessary create the output object in the initializer + self.dataOut = Voltage() + + def setup(self, path = None, + startDate = None, + endDate = None, + startTime = None, + endTime = None, + set = None, + expLabel = "", + ext = None, + online = False, + delay = 60, + walk = True): + ''' + In this method we should set all initial parameters. + + ''' + + self.isConfig = True + + def run(self, **kwargs): + ''' + This method will be called many times so here you should put all your code + ''' + + if not self.isConfig: + self.setup(**kwargs) + +class Writer(Operation): + ''' + classdocs + ''' + + def __init__(self): + ''' + Constructor + ''' + self.dataOut = None + + self.isConfig = False + + def setup(self, dataIn, path, blocksPerFile, set=0, ext=None): + ''' + In this method we should set all initial parameters. + + Input: + dataIn : Input data will also be outputa data + + ''' + self.dataOut = dataIn + + + + + + self.isConfig = True + + return + + def run(self, dataIn, **kwargs): + ''' + This method will be called many times so here you should put all your code + + Inputs: + + dataIn : object with the data + + ''' + + if not self.isConfig: + self.setup(dataIn, **kwargs) + \ No newline at end of file diff --git a/schainpy/model/io/jroIO_heispectra.py b/schainpy/model/io/jroIO_heispectra.py new file mode 100644 index 0000000..a3a7a04 --- /dev/null +++ b/schainpy/model/io/jroIO_heispectra.py @@ -0,0 +1,763 @@ +''' + +''' + +import os, sys +import time, datetime +import numpy +import fnmatch +import glob + +try: + import pyfits +except: + """ + """ + +from xml.etree.ElementTree import ElementTree + +from jroIO_base import isDoyFolder, isNumber +from model.proc.jroproc_base import Operation, ProcessingUnit + +class ParameterConf: + ELEMENTNAME = 'Parameter' + def __init__(self): + self.name = '' + self.value = '' + + def readXml(self, parmElement): + self.name = parmElement.get('name') + self.value = parmElement.get('value') + + def getElementName(self): + return self.ELEMENTNAME + +class Metadata: + + def __init__(self, filename): + self.parmConfObjList = [] + self.readXml(filename) + + def readXml(self, filename): + self.projectElement = None + self.procUnitConfObjDict = {} + self.projectElement = ElementTree().parse(filename) + self.project = self.projectElement.tag + + parmElementList = self.projectElement.getiterator(ParameterConf().getElementName()) + + for parmElement in parmElementList: + parmConfObj = ParameterConf() + parmConfObj.readXml(parmElement) + self.parmConfObjList.append(parmConfObj) + +class FitsWriter(Operation): + + def __init__(self): + self.isConfig = False + self.dataBlocksPerFile = None + self.blockIndex = 0 + self.flagIsNewFile = 1 + self.fitsObj = None + self.optchar = 'P' + self.ext = '.fits' + self.setFile = 0 + + def setFitsHeader(self, dataOut, metadatafile): + + header_data = pyfits.PrimaryHDU() + + metadata4fits = Metadata(metadatafile) + for parameter in metadata4fits.parmConfObjList: + parm_name = parameter.name + parm_value = parameter.value + +# if parm_value == 'fromdatadatetime': +# value = time.strftime("%b %d %Y %H:%M:%S", dataOut.datatime.timetuple()) +# elif parm_value == 'fromdataheights': +# value = dataOut.nHeights +# elif parm_value == 'fromdatachannel': +# value = dataOut.nChannels +# elif parm_value == 'fromdatasamples': +# value = dataOut.nFFTPoints +# else: +# value = parm_value + + header_data.header[parm_name] = parm_value + + + header_data.header['DATETIME'] = time.strftime("%b %d %Y %H:%M:%S", dataOut.datatime.timetuple()) + header_data.header['CHANNELLIST'] = str(dataOut.channelList) + header_data.header['NCHANNELS'] = dataOut.nChannels + #header_data.header['HEIGHTS'] = dataOut.heightList + header_data.header['NHEIGHTS'] = dataOut.nHeights + + header_data.header['IPPSECONDS'] = dataOut.ippSeconds + header_data.header['NCOHINT'] = dataOut.nCohInt + header_data.header['NINCOHINT'] = dataOut.nIncohInt + header_data.header['TIMEZONE'] = dataOut.timeZone + header_data.header['NBLOCK'] = self.blockIndex + + header_data.writeto(self.filename) + + self.addExtension(dataOut.heightList,'HEIGHTLIST') + + + def setup(self, dataOut, path, dataBlocksPerFile, metadatafile): + + self.path = path + self.dataOut = dataOut + self.metadatafile = metadatafile + self.dataBlocksPerFile = dataBlocksPerFile + + def open(self): + self.fitsObj = pyfits.open(self.filename, mode='update') + + + def addExtension(self, data, tagname): + self.open() + extension = pyfits.ImageHDU(data=data, name=tagname) + #extension.header['TAG'] = tagname + self.fitsObj.append(extension) + self.write() + + def addData(self, data): + self.open() + extension = pyfits.ImageHDU(data=data, name=self.fitsObj[0].header['DATATYPE']) + extension.header['UTCTIME'] = self.dataOut.utctime + self.fitsObj.append(extension) + self.blockIndex += 1 + self.fitsObj[0].header['NBLOCK'] = self.blockIndex + + self.write() + + def write(self): + + self.fitsObj.flush(verbose=True) + self.fitsObj.close() + + + def setNextFile(self): + + ext = self.ext + path = self.path + + timeTuple = time.localtime( self.dataOut.utctime) + subfolder = 'd%4.4d%3.3d' % (timeTuple.tm_year,timeTuple.tm_yday) + + fullpath = os.path.join( path, subfolder ) + if not( os.path.exists(fullpath) ): + os.mkdir(fullpath) + self.setFile = -1 #inicializo mi contador de seteo + else: + filesList = os.listdir( fullpath ) + if len( filesList ) > 0: + filesList = sorted( filesList, key=str.lower ) + filen = filesList[-1] + + if isNumber( filen[8:11] ): + self.setFile = int( filen[8:11] ) #inicializo mi contador de seteo al seteo del ultimo file + else: + self.setFile = -1 + else: + self.setFile = -1 #inicializo mi contador de seteo + + setFile = self.setFile + setFile += 1 + + file = '%s%4.4d%3.3d%3.3d%s' % (self.optchar, + timeTuple.tm_year, + timeTuple.tm_yday, + setFile, + ext ) + + filename = os.path.join( path, subfolder, file ) + + self.blockIndex = 0 + self.filename = filename + self.setFile = setFile + self.flagIsNewFile = 1 + + print 'Writing the file: %s'%self.filename + + self.setFitsHeader(self.dataOut, self.metadatafile) + + return 1 + + def writeBlock(self): + self.addData(self.dataOut.data_spc) + self.flagIsNewFile = 0 + + + def __setNewBlock(self): + + if self.flagIsNewFile: + return 1 + + if self.blockIndex < self.dataBlocksPerFile: + return 1 + + if not( self.setNextFile() ): + return 0 + + return 1 + + def writeNextBlock(self): + if not( self.__setNewBlock() ): + return 0 + self.writeBlock() + return 1 + + def putData(self): + if self.flagIsNewFile: + self.setNextFile() + self.writeNextBlock() + + def run(self, dataOut, **kwargs): + if not(self.isConfig): + self.setup(dataOut, **kwargs) + self.isConfig = True + self.putData() + + +class FitsReader(ProcessingUnit): + +# __TIMEZONE = time.timezone + + expName = None + datetimestr = None + utc = None + nChannels = None + nSamples = None + dataBlocksPerFile = None + comments = None + lastUTTime = None + header_dict = None + data = None + data_header_dict = None + + def __init__(self): + self.isConfig = False + self.ext = '.fits' + self.setFile = 0 + self.flagNoMoreFiles = 0 + self.flagIsNewFile = 1 + self.flagTimeBlock = None + self.fileIndex = None + self.filename = None + self.fileSize = None + self.fitsObj = None + self.timeZone = None + self.nReadBlocks = 0 + self.nTotalBlocks = 0 + self.dataOut = self.createObjByDefault() + self.maxTimeStep = 10# deberia ser definido por el usuario usando el metodo setup() + self.blockIndex = 1 + + def createObjByDefault(self): + + dataObj = Fits() + + return dataObj + + def isFileinThisTime(self, filename, startTime, endTime, useLocalTime=False): + try: + fitsObj = pyfits.open(filename,'readonly') + except: + raise IOError, "The file %s can't be opened" %(filename) + + header = fitsObj[0].header + struct_time = time.strptime(header['DATETIME'], "%b %d %Y %H:%M:%S") + utc = time.mktime(struct_time) - time.timezone #TIMEZONE debe ser un parametro del header FITS + + ltc = utc + if useLocalTime: + ltc -= time.timezone + thisDatetime = datetime.datetime.utcfromtimestamp(ltc) + thisTime = thisDatetime.time() + + if not ((startTime <= thisTime) and (endTime > thisTime)): + return None + + return thisDatetime + + def __setNextFileOnline(self): + raise ValueError, "No implemented" + + 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] + +# if not(self.__verifyFile(filename)): +# continue + + fileSize = os.path.getsize(filename) + fitsObj = pyfits.open(filename,'readonly') + break + + self.flagIsNewFile = 1 + self.fileIndex = idFile + self.filename = filename + self.fileSize = fileSize + self.fitsObj = fitsObj + self.blockIndex = 0 + print "Setting the file: %s"%self.filename + + return 1 + + def readHeader(self): + headerObj = self.fitsObj[0] + + self.header_dict = headerObj.header + if 'EXPNAME' in headerObj.header.keys(): + self.expName = headerObj.header['EXPNAME'] + + if 'DATATYPE' in headerObj.header.keys(): + self.dataType = headerObj.header['DATATYPE'] + + self.datetimestr = headerObj.header['DATETIME'] + channelList = headerObj.header['CHANNELLIST'] + channelList = channelList.split('[') + channelList = channelList[1].split(']') + channelList = channelList[0].split(',') + channelList = [int(ch) for ch in channelList] + self.channelList = channelList + self.nChannels = headerObj.header['NCHANNELS'] + self.nHeights = headerObj.header['NHEIGHTS'] + self.ippSeconds = headerObj.header['IPPSECONDS'] + self.nCohInt = headerObj.header['NCOHINT'] + self.nIncohInt = headerObj.header['NINCOHINT'] + self.dataBlocksPerFile = headerObj.header['NBLOCK'] + self.timeZone = headerObj.header['TIMEZONE'] + + self.timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt + + if 'COMMENT' in headerObj.header.keys(): + self.comments = headerObj.header['COMMENT'] + + self.readHeightList() + + def readHeightList(self): + self.blockIndex = self.blockIndex + 1 + obj = self.fitsObj[self.blockIndex] + self.heightList = obj.data + self.blockIndex = self.blockIndex + 1 + + def readExtension(self): + obj = self.fitsObj[self.blockIndex] + self.heightList = obj.data + self.blockIndex = self.blockIndex + 1 + + def setNextFile(self): + + if self.online: + newFile = self.__setNextFileOnline() + else: + newFile = self.__setNextFileOffline() + + if not(newFile): + return 0 + + self.readHeader() + + self.nReadBlocks = 0 +# self.blockIndex = 1 + return 1 + + def __searchFilesOffLine(self, + path, + startDate, + endDate, + startTime=datetime.time(0,0,0), + endTime=datetime.time(23,59,59), + set=None, + expLabel='', + ext='.fits', + walk=True): + + pathList = [] + + if not walk: + pathList.append(path) + + else: + dirList = [] + for thisPath in os.listdir(path): + if not os.path.isdir(os.path.join(path,thisPath)): + continue + if not isDoyFolder(thisPath): + continue + + dirList.append(thisPath) + + if not(dirList): + return None, None + + thisDate = startDate + + while(thisDate <= endDate): + year = thisDate.timetuple().tm_year + doy = thisDate.timetuple().tm_yday + + matchlist = fnmatch.filter(dirList, '?' + '%4.4d%3.3d' % (year,doy) + '*') + if len(matchlist) == 0: + thisDate += datetime.timedelta(1) + continue + for match in matchlist: + pathList.append(os.path.join(path,match,expLabel)) + + thisDate += datetime.timedelta(1) + + if pathList == []: + print "Any folder was found for the date range: %s-%s" %(startDate, endDate) + return None, None + + print "%d folder(s) was(were) found for the date range: %s - %s" %(len(pathList), startDate, endDate) + + filenameList = [] + datetimeList = [] + + for i in range(len(pathList)): + + thisPath = pathList[i] + + fileList = glob.glob1(thisPath, "*%s" %ext) + fileList.sort() + + for file in fileList: + + filename = os.path.join(thisPath,file) + thisDatetime = self.isFileinThisTime(filename, startTime, endTime) + + if not(thisDatetime): + continue + + filenameList.append(filename) + datetimeList.append(thisDatetime) + + if not(filenameList): + print "Any file was found for the time range %s - %s" %(startTime, endTime) + return None, None + + print "%d file(s) was(were) found for the time range: %s - %s" %(len(filenameList), startTime, endTime) + print + + for i in range(len(filenameList)): + print "%s -> [%s]" %(filenameList[i], datetimeList[i].ctime()) + + self.filenameList = filenameList + self.datetimeList = datetimeList + + return pathList, filenameList + + def setup(self, path=None, + startDate=None, + endDate=None, + startTime=datetime.time(0,0,0), + endTime=datetime.time(23,59,59), + set=0, + expLabel = "", + ext = None, + online = False, + delay = 60, + walk = True): + + if path == None: + raise ValueError, "The path is not valid" + + if ext == None: + ext = self.ext + + if not(online): + print "Searching files in offline mode ..." + pathList, filenameList = self.__searchFilesOffLine(path, startDate=startDate, endDate=endDate, + startTime=startTime, endTime=endTime, + set=set, expLabel=expLabel, ext=ext, + walk=walk) + + if not(pathList): + print "No *%s files into the folder %s \nfor the range: %s - %s"%(ext, path, + datetime.datetime.combine(startDate,startTime).ctime(), + datetime.datetime.combine(endDate,endTime).ctime()) + + sys.exit(-1) + + self.fileIndex = -1 + self.pathList = pathList + self.filenameList = filenameList + + self.online = online + self.delay = delay + ext = ext.lower() + self.ext = ext + + if not(self.setNextFile()): + if (startDate!=None) and (endDate!=None): + print "No files in range: %s - %s" %(datetime.datetime.combine(startDate,startTime).ctime(), datetime.datetime.combine(endDate,endTime).ctime()) + elif startDate != None: + print "No files in range: %s" %(datetime.datetime.combine(startDate,startTime).ctime()) + else: + print "No files" + + sys.exit(-1) + + + + def readBlock(self): + dataObj = self.fitsObj[self.blockIndex] + + self.data = dataObj.data + self.data_header_dict = dataObj.header + self.utc = self.data_header_dict['UTCTIME'] + + self.flagIsNewFile = 0 + self.blockIndex += 1 + self.nTotalBlocks += 1 + self.nReadBlocks += 1 + + return 1 + + def __jumpToLastBlock(self): + raise ValueError, "No implemented" + + def __waitNewBlock(self): + """ + Return 1 si se encontro un nuevo bloque de datos, 0 de otra forma. + + Si el modo de lectura es OffLine siempre retorn 0 + """ + if not self.online: + return 0 + + if (self.nReadBlocks >= self.dataBlocksPerFile): + return 0 + + currentPointer = self.fp.tell() + + neededSize = self.processingHeaderObj.blockSize + self.basicHeaderSize + + for nTries in range( self.nTries ): + + self.fp.close() + self.fp = open( self.filename, 'rb' ) + self.fp.seek( currentPointer ) + + self.fileSize = os.path.getsize( self.filename ) + currentSize = self.fileSize - currentPointer + + if ( currentSize >= neededSize ): + self.__rdBasicHeader() + return 1 + + print "\tWaiting %0.2f seconds for the next block, try %03d ..." % (self.delay, nTries+1) + time.sleep( self.delay ) + + + return 0 + + def __setNewBlock(self): + + if self.online: + self.__jumpToLastBlock() + + if self.flagIsNewFile: + return 1 + + self.lastUTTime = self.utc + + if self.online: + if self.__waitNewBlock(): + return 1 + + if self.nReadBlocks < self.dataBlocksPerFile: + return 1 + + if not(self.setNextFile()): + return 0 + + deltaTime = self.utc - self.lastUTTime + + self.flagTimeBlock = 0 + + if deltaTime > self.maxTimeStep: + self.flagTimeBlock = 1 + + return 1 + + + def readNextBlock(self): + if not(self.__setNewBlock()): + return 0 + + if not(self.readBlock()): + return 0 + + return 1 + + + def getData(self): + + if self.flagNoMoreFiles: + self.dataOut.flagNoData = True + print 'Process finished' + return 0 + + self.flagTimeBlock = 0 + self.flagIsNewBlock = 0 + + if not(self.readNextBlock()): + return 0 + + if self.data == None: + self.dataOut.flagNoData = True + return 0 + + self.dataOut.data = self.data + self.dataOut.data_header = self.data_header_dict + self.dataOut.utctime = self.utc + + self.dataOut.header = self.header_dict + self.dataOut.expName = self.expName + self.dataOut.nChannels = self.nChannels + self.dataOut.timeZone = self.timeZone + self.dataOut.dataBlocksPerFile = self.dataBlocksPerFile + self.dataOut.comments = self.comments + self.dataOut.timeInterval = self.timeInterval + self.dataOut.channelList = self.channelList + self.dataOut.heightList = self.heightList + self.dataOut.flagNoData = False + + return self.dataOut.data + + def run(self, **kwargs): + + if not(self.isConfig): + self.setup(**kwargs) + self.isConfig = True + + self.getData() + +class SpectraHeisWriter(Operation): +# set = None + setFile = None + idblock = None + doypath = None + subfolder = None + + def __init__(self): + self.wrObj = FITS() +# self.dataOut = dataOut + self.nTotalBlocks=0 +# self.set = None + self.setFile = None + self.idblock = 0 + self.wrpath = None + self.doypath = None + self.subfolder = None + self.isConfig = False + + def isNumber(str): + """ + Chequea si el conjunto de caracteres que componen un string puede ser convertidos a un numero. + + Excepciones: + Si un determinado string no puede ser convertido a numero + Input: + str, string al cual se le analiza para determinar si convertible a un numero o no + + Return: + True : si el string es uno numerico + False : no es un string numerico + """ + try: + float( str ) + return True + except: + return False + + def setup(self, dataOut, wrpath): + + if not(os.path.exists(wrpath)): + os.mkdir(wrpath) + + self.wrpath = wrpath +# self.setFile = 0 + self.dataOut = dataOut + + def putData(self): + name= time.localtime( self.dataOut.utctime) + ext=".fits" + + if self.doypath == None: + self.subfolder = 'F%4.4d%3.3d_%d' % (name.tm_year,name.tm_yday,time.mktime(datetime.datetime.now().timetuple())) + self.doypath = os.path.join( self.wrpath, self.subfolder ) + os.mkdir(self.doypath) + + if self.setFile == None: +# self.set = self.dataOut.set + self.setFile = 0 +# if self.set != self.dataOut.set: +## self.set = self.dataOut.set +# self.setFile = 0 + + #make the filename + file = 'D%4.4d%3.3d_%3.3d%s' % (name.tm_year,name.tm_yday,self.setFile,ext) + + filename = os.path.join(self.wrpath,self.subfolder, file) + + idblock = numpy.array([self.idblock],dtype="int64") + header=self.wrObj.cFImage(idblock=idblock, + year=time.gmtime(self.dataOut.utctime).tm_year, + month=time.gmtime(self.dataOut.utctime).tm_mon, + day=time.gmtime(self.dataOut.utctime).tm_mday, + hour=time.gmtime(self.dataOut.utctime).tm_hour, + minute=time.gmtime(self.dataOut.utctime).tm_min, + second=time.gmtime(self.dataOut.utctime).tm_sec) + + c=3E8 + deltaHeight = self.dataOut.heightList[1] - self.dataOut.heightList[0] + freq=numpy.arange(-1*self.dataOut.nHeights/2.,self.dataOut.nHeights/2.)*(c/(2*deltaHeight*1000)) + + colList = [] + + colFreq=self.wrObj.setColF(name="freq", format=str(self.dataOut.nFFTPoints)+'E', array=freq) + + colList.append(colFreq) + + nchannel=self.dataOut.nChannels + + for i in range(nchannel): + col = self.wrObj.writeData(name="PCh"+str(i+1), + format=str(self.dataOut.nFFTPoints)+'E', + data=10*numpy.log10(self.dataOut.data_spc[i,:])) + + colList.append(col) + + data=self.wrObj.Ctable(colList=colList) + + self.wrObj.CFile(header,data) + + self.wrObj.wFile(filename) + + #update the setFile + self.setFile += 1 + self.idblock += 1 + + return 1 + + def run(self, dataOut, **kwargs): + + if not(self.isConfig): + + self.setup(dataOut, **kwargs) + self.isConfig = True + + self.putData() \ No newline at end of file diff --git a/schainpy/model/io/jroIO_spectra.py b/schainpy/model/io/jroIO_spectra.py new file mode 100644 index 0000000..1fce650 --- /dev/null +++ b/schainpy/model/io/jroIO_spectra.py @@ -0,0 +1,761 @@ +''' +''' + +import numpy + +from jroIO_base import LOCALTIME, JRODataReader, JRODataWriter +from model.proc.jroproc_base import ProcessingUnit, Operation +from model.data.jroheaderIO import PROCFLAG, BasicHeader, SystemHeader, RadarControllerHeader, ProcessingHeader +from model.data.jrodata import Spectra + +class SpectraReader(JRODataReader, ProcessingUnit): + """ + Esta clase permite leer datos de espectros desde archivos procesados (.pdata). La lectura + de los datos siempre se realiza por bloques. Los datos leidos (array de 3 dimensiones) + son almacenados en tres buffer's para el Self Spectra, el Cross Spectra y el DC Channel. + + paresCanalesIguales * alturas * perfiles (Self Spectra) + paresCanalesDiferentes * alturas * perfiles (Cross Spectra) + canales * alturas (DC Channels) + + Esta clase contiene instancias (objetos) de las clases BasicHeader, SystemHeader, + RadarControllerHeader y Spectra. Los tres primeros se usan para almacenar informacion de la + cabecera de datos (metadata), y el cuarto (Spectra) para obtener y almacenar un bloque 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 = SpectraReader() + + readerObj.setup(dpath, startTime, endTime) + + while(True): + + readerObj.getData() + + print readerObj.data_spc + + print readerObj.data_cspc + + print readerObj.data_dc + + if readerObj.flagNoMoreFiles: + break + + """ + + pts2read_SelfSpectra = 0 + + pts2read_CrossSpectra = 0 + + pts2read_DCchannels = 0 + + ext = ".pdata" + + optchar = "P" + + dataOut = None + + nRdChannels = None + + nRdPairs = None + + rdPairList = [] + + def __init__(self): + """ + Inicializador de la clase SpectraReader para la lectura de datos de espectros. + + Inputs: + dataOut : Objeto de la clase Spectra. 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. + + Affected: + self.dataOut + + Return : None + """ + + #Eliminar de la base la herencia + ProcessingUnit.__init__(self) + +# self.isConfig = False + + self.pts2read_SelfSpectra = 0 + + self.pts2read_CrossSpectra = 0 + + self.pts2read_DCchannels = 0 + + self.datablock = None + + self.utc = None + + self.ext = ".pdata" + + self.optchar = "P" + + 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.lastUTTime = 0 + + self.maxTimeStep = 30 + + self.flagNoMoreFiles = 0 + + self.set = 0 + + self.path = None + + self.delay = 60 #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.flagTimeBlock = 0 + + self.flagIsNewBlock = 0 + + self.nTotalBlocks = 0 + + self.blocksize = 0 + + self.dataOut = self.createObjByDefault() + + self.profileIndex = 1 #Always + + + def createObjByDefault(self): + + dataObj = Spectra() + + return dataObj + + def __hasNotDataInBuffer(self): + return 1 + + + def getBlockDimension(self): + """ + Obtiene la cantidad de puntos a leer por cada bloque de datos + + Affected: + self.nRdChannels + self.nRdPairs + self.pts2read_SelfSpectra + self.pts2read_CrossSpectra + self.pts2read_DCchannels + self.blocksize + self.dataOut.nChannels + self.dataOut.nPairs + + Return: + None + """ + self.nRdChannels = 0 + self.nRdPairs = 0 + self.rdPairList = [] + + for i in range(0, self.processingHeaderObj.totalSpectra*2, 2): + if self.processingHeaderObj.spectraComb[i] == self.processingHeaderObj.spectraComb[i+1]: + self.nRdChannels = self.nRdChannels + 1 #par de canales iguales + else: + self.nRdPairs = self.nRdPairs + 1 #par de canales diferentes + self.rdPairList.append((self.processingHeaderObj.spectraComb[i], self.processingHeaderObj.spectraComb[i+1])) + + pts2read = self.processingHeaderObj.nHeights * self.processingHeaderObj.profilesPerBlock + + self.pts2read_SelfSpectra = int(self.nRdChannels * pts2read) + self.blocksize = self.pts2read_SelfSpectra + + if self.processingHeaderObj.flag_cspc: + self.pts2read_CrossSpectra = int(self.nRdPairs * pts2read) + self.blocksize += self.pts2read_CrossSpectra + + if self.processingHeaderObj.flag_dc: + self.pts2read_DCchannels = int(self.systemHeaderObj.nChannels * self.processingHeaderObj.nHeights) + self.blocksize += self.pts2read_DCchannels + +# self.blocksize = self.pts2read_SelfSpectra + self.pts2read_CrossSpectra + self.pts2read_DCchannels + + + def readBlock(self): + """ + 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 + + Return: None + + Variables afectadas: + + self.flagIsNewFile + self.flagIsNewBlock + self.nTotalBlocks + self.data_spc + self.data_cspc + self.data_dc + + Exceptions: + Si un bloque leido no es un bloque valido + """ + blockOk_flag = False + fpointer = self.fp.tell() + + spc = numpy.fromfile( self.fp, self.dtype[0], self.pts2read_SelfSpectra ) + spc = spc.reshape( (self.nRdChannels, self.processingHeaderObj.nHeights, self.processingHeaderObj.profilesPerBlock) ) #transforma a un arreglo 3D + + if self.processingHeaderObj.flag_cspc: + cspc = numpy.fromfile( self.fp, self.dtype, self.pts2read_CrossSpectra ) + cspc = cspc.reshape( (self.nRdPairs, self.processingHeaderObj.nHeights, self.processingHeaderObj.profilesPerBlock) ) #transforma a un arreglo 3D + + if self.processingHeaderObj.flag_dc: + dc = numpy.fromfile( self.fp, self.dtype, self.pts2read_DCchannels ) #int(self.processingHeaderObj.nHeights*self.systemHeaderObj.nChannels) ) + dc = dc.reshape( (self.systemHeaderObj.nChannels, self.processingHeaderObj.nHeights) ) #transforma a un arreglo 2D + + + if not(self.processingHeaderObj.shif_fft): + #desplaza a la derecha en el eje 2 determinadas posiciones + shift = int(self.processingHeaderObj.profilesPerBlock/2) + spc = numpy.roll( spc, shift , axis=2 ) + + if self.processingHeaderObj.flag_cspc: + #desplaza a la derecha en el eje 2 determinadas posiciones + cspc = numpy.roll( cspc, shift, axis=2 ) + +# self.processingHeaderObj.shif_fft = True + + spc = numpy.transpose( spc, (0,2,1) ) + self.data_spc = spc + + if self.processingHeaderObj.flag_cspc: + cspc = numpy.transpose( cspc, (0,2,1) ) + self.data_cspc = cspc['real'] + cspc['imag']*1j + else: + self.data_cspc = None + + if self.processingHeaderObj.flag_dc: + self.data_dc = dc['real'] + dc['imag']*1j + else: + self.data_dc = None + + 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() + +# self.dataOut.ippSeconds = self.ippSeconds + + self.dataOut.timeInterval = self.radarControllerHeaderObj.ippSeconds * self.processingHeaderObj.nCohInt * self.processingHeaderObj.nIncohInt * self.processingHeaderObj.profilesPerBlock + + self.dataOut.dtype = self.dtype + +# self.dataOut.nPairs = self.nPairs + + self.dataOut.pairsList = self.rdPairList + + self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock + + self.dataOut.nFFTPoints = self.processingHeaderObj.profilesPerBlock + + self.dataOut.nCohInt = self.processingHeaderObj.nCohInt + + self.dataOut.nIncohInt = self.processingHeaderObj.nIncohInt + + xf = self.processingHeaderObj.firstHeight + self.processingHeaderObj.nHeights*self.processingHeaderObj.deltaHeight + + self.dataOut.heightList = numpy.arange(self.processingHeaderObj.firstHeight, xf, self.processingHeaderObj.deltaHeight) + + self.dataOut.channelList = range(self.systemHeaderObj.nChannels) + + self.dataOut.flagShiftFFT = self.processingHeaderObj.shif_fft + + self.dataOut.flagDecodeData = False #asumo q la data no esta decodificada + + self.dataOut.flagDeflipData = True #asumo q la data no esta sin flip + + 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.flagDecodeData = True + + def getData(self): + """ + First method to execute before "RUN" is called. + + Copia el buffer de lectura a la clase "Spectra", + 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" + + Return: + 0 : Si no hay mas archivos disponibles + 1 : Si hizo una buena copia del buffer + + Affected: + self.dataOut + + self.flagTimeBlock + self.flagIsNewBlock + """ + + if self.flagNoMoreFiles: + self.dataOut.flagNoData = True + print 'Process finished' + return 0 + + self.flagTimeBlock = 0 + self.flagIsNewBlock = 0 + + if self.__hasNotDataInBuffer(): + + if not( self.readNextBlock() ): + self.dataOut.flagNoData = True + return 0 + + #data es un numpy array de 3 dmensiones (perfiles, alturas y canales) + + if self.data_dc == None: + self.dataOut.flagNoData = True + return 0 + + self.getBasicHeader() + + self.getFirstHeader() + + self.dataOut.data_spc = self.data_spc + + self.dataOut.data_cspc = self.data_cspc + + self.dataOut.data_dc = self.data_dc + + self.dataOut.flagNoData = False + + self.dataOut.realtime = self.online + + return self.dataOut.data_spc + +class SpectraWriter(JRODataWriter, Operation): + + """ + Esta clase permite escribir datos de espectros a archivos procesados (.pdata). La escritura + de los datos siempre se realiza por bloques. + """ + + ext = ".pdata" + + optchar = "P" + + shape_spc_Buffer = None + + shape_cspc_Buffer = None + + shape_dc_Buffer = None + + data_spc = None + + data_cspc = None + + data_dc = None + +# dataOut = None + + def __init__(self): + """ + Inicializador de la clase SpectraWriter para la escritura de datos de espectros. + + Affected: + self.dataOut + self.basicHeaderObj + self.systemHeaderObj + self.radarControllerHeaderObj + self.processingHeaderObj + + Return: None + """ + + Operation.__init__(self) + + self.isConfig = False + + self.nTotalBlocks = 0 + + self.data_spc = None + + self.data_cspc = None + + self.data_dc = None + + self.fp = None + + self.flagIsNewFile = 1 + + self.nTotalBlocks = 0 + + self.flagIsNewBlock = 0 + + self.setFile = None + + self.dtype = None + + self.path = None + + self.noMoreFiles = 0 + + self.filename = None + + self.basicHeaderObj = BasicHeader(LOCALTIME) + + self.systemHeaderObj = SystemHeader() + + self.radarControllerHeaderObj = RadarControllerHeader() + + self.processingHeaderObj = ProcessingHeader() + + + def hasAllDataInBuffer(self): + return 1 + + + 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.shape_spc_Buffer = (self.dataOut.nChannels, + self.processingHeaderObj.nHeights, + self.processingHeaderObj.profilesPerBlock) + + self.shape_cspc_Buffer = (self.dataOut.nPairs, + self.processingHeaderObj.nHeights, + self.processingHeaderObj.profilesPerBlock) + + self.shape_dc_Buffer = (self.dataOut.nChannels, + self.processingHeaderObj.nHeights) + + + def writeBlock(self): + """ + Escribe el buffer en el file designado + + Affected: + self.data_spc + self.data_cspc + self.data_dc + self.flagIsNewFile + self.flagIsNewBlock + self.nTotalBlocks + self.nWriteBlocks + + Return: None + """ + + spc = numpy.transpose( self.data_spc, (0,2,1) ) + if not( self.processingHeaderObj.shif_fft ): + spc = numpy.roll( spc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones + data = spc.reshape((-1)) + data = data.astype(self.dtype[0]) + data.tofile(self.fp) + + if self.data_cspc != None: + data = numpy.zeros( self.shape_cspc_Buffer, self.dtype ) + cspc = numpy.transpose( self.data_cspc, (0,2,1) ) + if not( self.processingHeaderObj.shif_fft ): + cspc = numpy.roll( cspc, self.processingHeaderObj.profilesPerBlock/2, axis=2 ) #desplaza a la derecha en el eje 2 determinadas posiciones + data['real'] = cspc.real + data['imag'] = cspc.imag + data = data.reshape((-1)) + data.tofile(self.fp) + + if self.data_dc != None: + data = numpy.zeros( self.shape_dc_Buffer, self.dtype ) + dc = self.data_dc + data['real'] = dc.real + data['imag'] = dc.imag + data = data.reshape((-1)) + data.tofile(self.fp) + + self.data_spc.fill(0) + + if self.data_dc != None: + self.data_dc.fill(0) + + if self.data_cspc != None: + self.data_cspc.fill(0) + + self.flagIsNewFile = 0 + self.flagIsNewBlock = 1 + self.nTotalBlocks += 1 + self.nWriteBlocks += 1 + self.blockIndex += 1 + + + def putData(self): + """ + Setea un bloque de datos y luego los escribe en un file + + Affected: + self.data_spc + self.data_cspc + self.data_dc + + Return: + 0 : Si no hay data o no hay mas files que puedan escribirse + 1 : Si se escribio la data de un bloque en un file + """ + + if self.dataOut.flagNoData: + return 0 + + self.flagIsNewBlock = 0 + + if self.dataOut.flagTimeBlock: + self.data_spc.fill(0) + self.data_cspc.fill(0) + self.data_dc.fill(0) + self.setNextFile() + + if self.flagIsNewFile == 0: + self.setBasicHeader() + + self.data_spc = self.dataOut.data_spc.copy() + if self.dataOut.data_cspc != None: + self.data_cspc = self.dataOut.data_cspc.copy() + self.data_dc = self.dataOut.data_dc.copy() + + # #self.processingHeaderObj.dataBlocksPerFile) + if self.hasAllDataInBuffer(): +# self.setFirstHeader() + self.writeNextBlock() + + return 1 + + + def __getProcessFlags(self): + + processFlags = 0 + + dtype0 = numpy.dtype([('real',' 1: + processFlags += PROCFLAG.INCOHERENT_INTEGRATION + + if self.dataOut.data_dc != None: + processFlags += PROCFLAG.SAVE_CHANNELS_DC + + return processFlags + + + def __getBlockSize(self): + ''' + Este metodos determina el cantidad de bytes para un bloque de datos de tipo Spectra + ''' + + dtype0 = numpy.dtype([('real',' 0: + channelList = [] + for channel in range(self.dataOut.nChannels): + channelList.append(channel) + channelList.append(channel) + + pairsList = [] + if self.dataOut.nPairs > 0: + for pair in self.dataOut.pairsList: + pairsList.append(pair[0]) + pairsList.append(pair[1]) + + spectraComb = channelList + pairsList + spectraComb = numpy.array(spectraComb,dtype="u1") + self.processingHeaderObj.spectraComb = spectraComb + sizeOfSpcComb = len(spectraComb) + processingHeaderSize += sizeOfSpcComb + +# The processing header should not have information about code +# if self.dataOut.code != None: +# self.processingHeaderObj.code = self.dataOut.code +# self.processingHeaderObj.nCode = self.dataOut.nCode +# self.processingHeaderObj.nBaud = self.dataOut.nBaud +# nCodeSize = 4 # bytes +# nBaudSize = 4 # bytes +# codeSize = 4 # bytes +# sizeOfCode = int(nCodeSize + nBaudSize + codeSize * self.dataOut.nCode * self.dataOut.nBaud) +# processingHeaderSize += sizeOfCode + + 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 + sizeOfFirstHeight = 4 + sizeOfdeltaHeight = 4 + sizeOfnHeights = 4 + sizeOfWindows = (sizeOfFirstHeight + sizeOfdeltaHeight + sizeOfnHeights)*self.processingHeaderObj.nWindows + processingHeaderSize += sizeOfWindows + + self.processingHeaderObj.size = processingHeaderSize + diff --git a/schainpy/model/io/jroIO_voltage.py b/schainpy/model/io/jroIO_voltage.py new file mode 100644 index 0000000..5d4a06c --- /dev/null +++ b/schainpy/model/io/jroIO_voltage.py @@ -0,0 +1,588 @@ +''' + +''' +import numpy + +from jroIO_base import LOCALTIME, JRODataReader, JRODataWriter +from model.proc.jroproc_base import ProcessingUnit, Operation +from model.data.jroheaderIO import PROCFLAG, BasicHeader, SystemHeader, RadarControllerHeader, ProcessingHeader +from 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.flagTimeBlock = 0 + + self.flagIsNewBlock = 0 + + self.nTotalBlocks = 0 + + self.blocksize = 0 + + self.dataOut = self.createObjByDefault() + + 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() + +# self.dataOut.ippSeconds = self.ippSeconds + + 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 + + xf = self.processingHeaderObj.firstHeight + self.processingHeaderObj.nHeights*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 y la copia a la clase "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 en 1. + + Return: + data : retorna un perfil de voltages (alturas * canales) copiados desde el + buffer. Si no hay mas archivos a leer retorna None. + + Variables afectadas: + self.dataOut + self.profileIndex + + Affected: + self.dataOut + self.profileIndex + self.flagTimeBlock + self.flagIsNewBlock + """ + + if self.flagNoMoreFiles: + self.dataOut.flagNoData = True + print 'Process finished' + return 0 + + self.flagTimeBlock = 0 + self.flagIsNewBlock = 0 + + if self.__hasNotDataInBuffer(): + + if not( self.readNextBlock() ): + return 0 + + self.getFirstHeader() + + if self.datablock == None: + self.dataOut.flagNoData = True + return 0 + + self.dataOut.data = self.datablock[:,self.profileIndex,:] + + self.dataOut.flagNoData = False + + self.getBasicHeader() + + self.profileIndex += 1 + + 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 + + def putData(self): + """ + Setea un bloque de datos y luego los escribe en un file + + Affected: + self.flagIsNewBlock + self.profileIndex + + Return: + 0 : Si no hay data o no hay mas files que puedan escribirse + 1 : Si se escribio la data de un bloque en un file + """ + if self.dataOut.flagNoData: + return 0 + + self.flagIsNewBlock = 0 + + if self.dataOut.flagTimeBlock: + + 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',' 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',' 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: + 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: + 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 \ No newline at end of file diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py new file mode 100644 index 0000000..813d9ce --- /dev/null +++ b/schainpy/model/proc/jroproc_spectra.py @@ -0,0 +1,1024 @@ +import numpy + +from jroproc_base import ProcessingUnit, Operation +from model.data.jrodata import Spectra + +class SpectraProc(ProcessingUnit): + + def __init__(self): + + ProcessingUnit.__init__(self) + + self.buffer = None + self.firstdatatime = None + self.profIndex = 0 + self.dataOut = Spectra() + + 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',' 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 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: + 1 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-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: + self.__buffer_cspc = None + else: + self.__buffer_cspc += data_cspc + + if data_dc == 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: + self.__buffer_spc = data_spc + + if data_cspc == None: + self.__buffer_cspc = None + else: + self.__buffer_cspc += data_cspc + + if data_dc == 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: + 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 + +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.profiles = data.shape[1] + 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.profiles + + def run(self, dataOut, m): + + dataOut.flagNoData = True + + if not self.isConfig: + self.setup(dataOut.data, m, 1) + self.isConfig = True + + 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 * 5 + dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight) + +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 isProfileInRange(self, minIndex, maxIndex): + + if self.profileIndex < minIndex: + return False + + if self.profileIndex > maxIndex: + return False + + return True + + def isProfileInList(self, profileList): + + if self.profileIndex not in profileList: + return False + + return True + + def run(self, dataOut, profileList=None, profileRangeList=None, beam=None): + + dataOut.flagNoData = True + self.nProfiles = dataOut.nProfiles + + if profileList != None: + if self.isProfileInList(profileList): + dataOut.flagNoData = False + + self.incIndex() + return 1 + + + elif profileRangeList != None: + minIndex = profileRangeList[0] + maxIndex = profileRangeList[1] + if self.isProfileInRange(minIndex, maxIndex): + dataOut.flagNoData = False + + self.incIndex() + return 1 + elif beam != None: + if self.isProfileInList(dataOut.beamRangeDict[beam]): + dataOut.flagNoData = False + + self.incIndex() + return 1 + + else: + raise ValueError, "ProfileSelector needs profileList or profileRangeList" + + return 0 diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py new file mode 100644 index 0000000..3c5853b --- /dev/null +++ b/schainpy/model/proc/jroproc_voltage.py @@ -0,0 +1,518 @@ +import numpy + +from jroproc_base import ProcessingUnit, Operation +from 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): + 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 +# +# 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: + 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 = 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: + 1 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: + 1 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-1 +# raise ValueError, "some value in (%d,%d) is not valid" % (minIndex, maxIndex) + +# nHeights = maxIndex - minIndex + 1 + + #voltage + data = self.dataOut.data[:,minIndex:maxIndex+1] + +# firstHeight = self.dataOut.heightList[minIndex] + + self.dataOut.data = data + self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex+1] + + 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.data.shape[1] % window + buffer = self.dataOut.data[:,0:self.dataOut.data.shape[1]-r] + buffer = buffer.reshape(self.dataOut.data.shape[0],self.dataOut.data.shape[1]/window,window) + buffer = numpy.sum(buffer,2) + self.dataOut.data = buffer + self.dataOut.heightList = numpy.arange(self.dataOut.heightList[0],newdelta*(self.dataOut.nHeights-r)/window,newdelta) + self.dataOut.windowOfFilter = window + + def deFlip(self): + self.dataOut.data *= self.flip + self.flip *= -1. + + 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): + """ + 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: + 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: + 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, 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.isConfig = False + + def setup(self, code, shape): + + self.__profIndex = 0 + + self.code = code + + self.nCode = len(code) + self.nBaud = len(code[0]) + + self.__nChannels, self.__nHeis = shape + + __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[:,:-self.nBaud+1] + + return datadec + + def convolutionInFreqOpt(self, data): + + fft_code = self.fft_code[self.__profIndex].reshape(1,-1) + + data = cfunctions.decoder(fft_code, data) + + datadec = data[:,:-self.nBaud+1] + + 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='valid') + + return self.datadecTime + + def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0): + + if code == None: + code = dataOut.code + else: + code = numpy.array(code).reshape(nCode,nBaud) + dataOut.code = code + dataOut.nCode = nCode + dataOut.nBaud = nBaud + dataOut.radarControllerHeaderObj.code = code + dataOut.radarControllerHeaderObj.nCode = nCode + dataOut.radarControllerHeaderObj.nBaud = nBaud + + + if not self.isConfig: + + self.setup(code, dataOut.data.shape) + self.isConfig = True + + 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.data = datadec + + dataOut.heightList = dataOut.heightList[0:self.ndatadec] + + dataOut.flagDecodeData = True #asumo q la data no 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 diff --git a/schainpy/model/proc/jroprocessing.py b/schainpy/model/proc/jroprocessing.py new file mode 100644 index 0000000..ed54ca9 --- /dev/null +++ b/schainpy/model/proc/jroprocessing.py @@ -0,0 +1,3 @@ +from jroproc_voltage import * +from jroproc_spectra import * +from jroproc_heispectra import * \ No newline at end of file