From 162740cfc9304f3e53f4893cf52237c6ddcc91d7 2020-06-17 10:26:10 From: avaldezp22 Date: 2020-06-17 10:26:10 Subject: [PATCH] update for weather radar options --- diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index 042ae22..39a5287 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -360,7 +360,8 @@ class Voltage(JROData): # data es un numpy array de 2 dmensiones (canales, alturas) data = None - + data_intensity = None + data_velocity = None def __init__(self): ''' Constructor @@ -555,7 +556,7 @@ class Spectra(JROData): deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor) velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - + if self.nmodes: return velrange/self.nmodes else: @@ -1111,7 +1112,7 @@ class PlotterData(object): MAXNUMY = 100 def __init__(self, code, throttle_value, exp_code, buffering=True, snr=False): - + self.key = code self.throttle = throttle_value self.exp_code = exp_code @@ -1139,6 +1140,7 @@ class PlotterData(object): for plot in self.plottypes: self.data[plot] = {} + def __str__(self): dum = ['{}{}'.format(key, self.shape(key)) for key in self.data] return 'Data[{}][{}]'.format(';'.join(dum), len(self.__times)) @@ -1147,7 +1149,7 @@ class PlotterData(object): return len(self.__times) def __getitem__(self, key): - + if key not in self.data: raise KeyError(log.error('Missing key: {}'.format(key))) if 'spc' in key or not self.buffering: @@ -1167,7 +1169,6 @@ class PlotterData(object): ''' Configure object ''' - self.type = '' self.ready = False self.data = {} @@ -1180,7 +1181,7 @@ class PlotterData(object): elif 'spc_moments' == plot: plot = 'moments' self.data[plot] = {} - + if 'spc' in self.data or 'rti' in self.data or 'cspc' in self.data or 'moments' in self.data: self.data['noise'] = {} self.data['rti'] = {} @@ -1188,7 +1189,7 @@ class PlotterData(object): self.plottypes.append('noise') if 'rti' not in self.plottypes: self.plottypes.append('rti') - + def shape(self, key): ''' Get the shape of the one-element data for the given key @@ -1204,20 +1205,19 @@ class PlotterData(object): ''' Update data object with new dataOut ''' - if tm in self.__times: return self.profileIndex = dataOut.profileIndex self.tm = tm self.type = dataOut.type self.parameters = getattr(dataOut, 'parameters', []) - + if hasattr(dataOut, 'meta'): self.meta.update(dataOut.meta) - + if hasattr(dataOut, 'pairsList'): self.pairs = dataOut.pairsList - + self.interval = dataOut.getTimeInterval() self.localtime = dataOut.useLocalTime if True in ['spc' in ptype for ptype in self.plottypes]: @@ -1227,7 +1227,6 @@ class PlotterData(object): self.__heights.append(dataOut.heightList) self.__all_heights.update(dataOut.heightList) self.__times.append(tm) - for plot in self.plottypes: if plot in ('spc', 'spc_moments', 'spc_cut'): z = dataOut.data_spc/dataOut.normFactor @@ -1260,8 +1259,16 @@ class PlotterData(object): if plot == 'scope': buffer = dataOut.data self.flagDataAsBlock = dataOut.flagDataAsBlock - self.nProfiles = dataOut.nProfiles - + self.nProfiles = dataOut.nProfiles + if plot == 'pp_power': + buffer = dataOut.data_intensity + self.flagDataAsBlock = dataOut.flagDataAsBlock + self.nProfiles = dataOut.nProfiles + if plot == 'pp_velocity': + buffer = dataOut.data_velocity + self.flagDataAsBlock = dataOut.flagDataAsBlock + self.nProfiles = dataOut.nProfiles + if plot == 'spc': self.data['spc'] = buffer elif plot == 'cspc': @@ -1283,7 +1290,7 @@ class PlotterData(object): if buffer is None: self.flagNoData = True - raise schainpy.admin.SchainWarning('Attribute data_{} is empty'.format(self.key)) + raise schainpy.admin.SchainWarning('Attribute data_{} is empty'.format(self.key)) def normalize_heights(self): ''' @@ -1340,7 +1347,7 @@ class PlotterData(object): else: meta['xrange'] = [] - meta.update(self.meta) + meta.update(self.meta) ret['metadata'] = meta return json.dumps(ret) diff --git a/schainpy/model/graphics/jroplot_voltage.py b/schainpy/model/graphics/jroplot_voltage.py index 7755fe0..0ffc027 100644 --- a/schainpy/model/graphics/jroplot_voltage.py +++ b/schainpy/model/graphics/jroplot_voltage.py @@ -14,12 +14,12 @@ class ScopePlot(Plot): ''' Plot for Scope - ''' + ''' CODE = 'scope' plot_name = 'Scope' plot_type = 'scatter' - + def setup(self): self.xaxis = 'Range (Km)' @@ -33,20 +33,20 @@ class ScopePlot(Plot): self.height = 4 def plot_iq(self, x, y, channelIndexList, thisDatetime, wintitle): - + yreal = y[channelIndexList,:].real yimag = y[channelIndexList,:].imag title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y")) self.xlabel = "Range (Km)" self.ylabel = "Intensity - IQ" - + self.y = yreal self.x = x self.xmin = min(x) self.xmax = max(x) - - self.titles[0] = title + + self.titles[0] = title for i,ax in enumerate(self.axes): title = "Channel %d" %(i) @@ -56,79 +56,178 @@ class ScopePlot(Plot): else: ax.plt_r.set_data(x, yreal[i,:]) ax.plt_i.set_data(x, yimag[i,:]) - + def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle): y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:]) yreal = y.real + yreal = 10*numpy.log10(yreal) self.y = yreal title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y")) self.xlabel = "Range (Km)" self.ylabel = "Intensity" self.xmin = min(x) self.xmax = max(x) - - + + self.titles[0] = title for i,ax in enumerate(self.axes): title = "Channel %d" %(i) - + ychannel = yreal[i,:] - - if ax.firsttime: + + if ax.firsttime: ax.plt_r = ax.plot(x, ychannel)[0] else: #pass ax.plt_r.set_data(x, ychannel) - + + def plot_weatherpower(self, x, y, channelIndexList, thisDatetime, wintitle): + + + y = y[channelIndexList,:] + yreal = y.real + yreal = 10*numpy.log10(yreal) + self.y = yreal + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + self.xlabel = "Range (Km)" + self.ylabel = "Intensity" + self.xmin = min(x) + self.xmax = max(x) + + self.titles[0] =title + for i,ax in enumerate(self.axes): + title = "Channel %d" %(i) + + ychannel = yreal[i,:] + + if ax.firsttime: + ax.plt_r = ax.plot(x, ychannel)[0] + else: + #pass + ax.plt_r.set_data(x, ychannel) + + def plot_weathervelocity(self, x, y, channelIndexList, thisDatetime, wintitle): + + x = x[channelIndexList,:] + yreal = y + self.y = yreal + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + self.xlabel = "Velocity (m/s)" + self.ylabel = "Range (Km)" + self.xmin = numpy.min(x) + self.xmax = numpy.max(x) + self.titles[0] =title + for i,ax in enumerate(self.axes): + title = "Channel %d" %(i) + xchannel = x[i,:] + if ax.firsttime: + ax.plt_r = ax.plot(xchannel, yreal)[0] + else: + #pass + ax.plt_r.set_data(xchannel, yreal) def plot(self): - if self.channels: channels = self.channels else: channels = self.data.channels thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1]) - - scope = self.data['scope'] - + if self.CODE == "pp_power": + scope = self.data['pp_power'] + elif self.CODE == "pp_velocity": + scope = self.data["pp_velocity"] + else: + scope =self.data["scope"] + if self.data.flagDataAsBlock: - + for i in range(self.data.nProfiles): wintitle1 = " [Profile = %d] " %i - - if self.type == "power": - self.plot_power(self.data.heights, - scope[:,i,:], - channels, - thisDatetime, - wintitle1 + if self.CODE =="scope": + if self.type == "power": + self.plot_power(self.data.heights, + scope[:,i,:], + channels, + thisDatetime, + wintitle1 + ) + + if self.type == "iq": + self.plot_iq(self.data.heights, + scope[:,i,:], + channels, + thisDatetime, + wintitle1 ) + if self.CODE=="pp_power": + self.plot_weatherpower(self.data.heights, + scope[:,i,:], + channels, + thisDatetime, + wintitle + ) + if self.CODE=="pp_velocity": + self.plot_weathervelocity(scope[:,i,:], + self.data.heights, + channels, + thisDatetime, + wintitle + ) + else: + wintitle = " [Profile = %d] " %self.data.profileIndex + if self.CODE== "scope": + if self.type == "power": + self.plot_power(self.data.heights, + scope, + channels, + thisDatetime, + wintitle + ) if self.type == "iq": - self.plot_iq(self.data.heights, - scope[:,i,:], - channels, - thisDatetime, - wintitle1 + self.plot_iq(self.data.heights, + scope, + channels, + thisDatetime, + wintitle ) - else: - wintitle = " [Profile = %d] " %self.data.profileIndex - - if self.type == "power": - self.plot_power(self.data.heights, - scope, - channels, - thisDatetime, - wintitle - ) - - if self.type == "iq": - self.plot_iq(self.data.heights, - scope, - channels, - thisDatetime, - wintitle - ) + if self.CODE=="pp_power": + self.plot_weatherpower(self.data.heights, + scope, + channels, + thisDatetime, + wintitle + ) + if self.CODE=="pp_velocity": + self.plot_weathervelocity(scope, + self.data.heights, + channels, + thisDatetime, + wintitle + ) + + + +class PulsepairPowerPlot(ScopePlot): + ''' + Plot for + ''' + + CODE = 'pp_power' + plot_name = 'PulsepairPower' + plot_type = 'scatter' + buffering = False + + + +class PulsepairVelocityPlot(ScopePlot): + ''' + Plot for + ''' + CODE = 'pp_velocity' + plot_name = 'PulsepairVelocity' + plot_type = 'scatter' + buffering = False diff --git a/schainpy/model/io/jroIO_simulator.py b/schainpy/model/io/jroIO_simulator.py new file mode 100644 index 0000000..bc2bbe1 --- /dev/null +++ b/schainpy/model/io/jroIO_simulator.py @@ -0,0 +1,473 @@ +import numpy,math,random,time +#---------------1 Heredamos JRODatareader +from schainpy.model.io.jroIO_base import * +#---------------2 Heredamos las propiedades de ProcessingUnit +from schainpy.model.proc.jroproc_base import ProcessingUnit,Operation,MPDecorator +#---------------3 Importaremos las clases BascicHeader, SystemHeader, RadarControlHeader, ProcessingHeader +from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader,SystemHeader,RadarControllerHeader, ProcessingHeader +#---------------4 Importaremos el objeto Voltge +from schainpy.model.data.jrodata import Voltage + +class SimulatorReader(JRODataReader, ProcessingUnit): + incIntFactor = 1 + nFFTPoints = 0 + FixPP_IncInt = 1 + FixRCP_IPP = 1000 + FixPP_CohInt = 1 + Tau_0 = 250 + AcqH0_0 = 70 + H0 = AcqH0_0 + AcqDH_0 = 1.25 + DH0 = AcqDH_0 + Bauds = 32 + BaudWidth = None + FixRCP_TXA = 40 + FixRCP_TXB = 70 + fAngle = 2.0*math.pi*(1/16) + DC_level = 500 + stdev = 8 + Num_Codes = 2 + #code0 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1]) + #code1 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0]) + #Dyn_snCode = numpy.array([Num_Codes,Bauds]) + Dyn_snCode = None + Samples = 200 + channels = 5 + pulses = None + Reference = None + pulse_size = None + prof_gen = None + Fdoppler = 100 + Hdoppler = 36 + Adoppler = 300 + frequency = 9345 + nTotalReadFiles = 1000 + + def __init__(self): + """ + Inicializador de la clases SimulatorReader para + generar datos de voltage simulados. + Input: + dataOut: Objeto de la clase Voltage. + Este Objeto sera utilizado apra almacenar + un perfil de datos cada vez qe se haga + un requerimiento (getData) + """ + ProcessingUnit.__init__(self) + print(" [ START ] init - Metodo Simulator Reader") + + self.isConfig = False + self.basicHeaderObj = BasicHeader(LOCALTIME) + self.systemHeaderObj = SystemHeader() + self.radarControllerHeaderObj = RadarControllerHeader() + self.processingHeaderObj = ProcessingHeader() + self.profileIndex = 2**32-1 + self.dataOut = Voltage() + #code0 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,1,1,1,0,1,1,0,1,0,0,0,1,1,1,0,1]) + code0 = numpy.array([1,1,1,-1,1,1,-1,1,1,1,1,-1,-1,-1,1,-1,1,1,1,-1,1,1,-1,1,-1,-1,-1,1,1,1,-1,1]) + #code1 = numpy.array([1,1,1,0,1,1,0,1,1,1,1,0,0,0,1,0,0,0,0,1,0,0,1,0,1,1,1,0,0,0,1,0]) + code1 = numpy.array([1,1,1,-1,1,1,-1,1,1,1,1,-1,-1,-1,1,-1,-1,-1,-1,1,-1,-1,1,-1,1,1,1,-1,-1,-1,1,-1]) + #self.Dyn_snCode = numpy.array([code0,code1]) + self.Dyn_snCode = None + + def set_kwargs(self, **kwargs): + for key, value in kwargs.items(): + setattr(self, key, value) + + def __hasNotDataInBuffer(self): + + if self.profileIndex >= self.processingHeaderObj.profilesPerBlock* self.nTxs: + if self.nReadBlocks>0: + tmp = self.dataOut.utctime + tmp_utc = int(self.dataOut.utctime) + tmp_milisecond = int((tmp-tmp_utc)*1000) + self.basicHeaderObj.utc = tmp_utc + self.basicHeaderObj.miliSecond= tmp_milisecond + return 1 + return 0 + + def setNextFile(self): + """Set the next file to be readed open it and parse de file header""" + + if (self.nReadBlocks >= self.processingHeaderObj.dataBlocksPerFile): + self.nReadFiles=self.nReadFiles+1 + if self.nReadFiles ==self.nTotalReadFiles: + self.flagNoMoreFiles=1 + raise schainpy.admin.SchainWarning('No more files to read') + + print('------------------- [Opening file] ------------------------------',self.nReadFiles) + self.nReadBlocks = 0 + + def __setNewBlock(self): + self.setNextFile() + if self.flagIsNewFile: + return 1 + + def readNextBlock(self): + while True: + self.__setNewBlock() + if not(self.readBlock()): + return 0 + self.getBasicHeader() + break + if self.verbose: + print("[Reading] Block No. %d/%d -> %s" %(self.nReadBlocks, + self.processingHeaderObj.dataBlocksPerFile, + self.dataOut.datatime.ctime()) ) + return 1 + + def getFirstHeader(self): + self.getBasicHeader() + self.dataOut.processingHeaderObj = self.processingHeaderObj.copy() + self.dataOut.systemHeaderObj = self.systemHeaderObj.copy() + self.dataOut.radarControllerHeaderObj = self.radarControllerHeaderObj.copy() + + self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock + self.dataOut.heightList = numpy.arange(self.processingHeaderObj.nHeights) * self.processingHeaderObj.deltaHeight + self.processingHeaderObj.firstHeight + self.dataOut.channelList = list(range(self.systemHeaderObj.nChannels)) + self.dataOut.nCohInt = self.processingHeaderObj.nCohInt + # asumo q la data no esta decodificada + self.dataOut.flagDecodeData = self.processingHeaderObj.flag_decode + # asumo q la data no esta sin flip + self.dataOut.flagDeflipData = self.processingHeaderObj.flag_deflip + self.dataOut.flagShiftFFT = self.processingHeaderObj.shif_fft + self.dataOut.frequency = self.frequency + + def getBasicHeader(self): + self.dataOut.utctime = self.basicHeaderObj.utc + self.basicHeaderObj.miliSecond / \ + 1000. + self.profileIndex * self.radarControllerHeaderObj.ippSeconds + + self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock + self.dataOut.timeZone = self.basicHeaderObj.timeZone + self.dataOut.dstFlag = self.basicHeaderObj.dstFlag + self.dataOut.errorCount = self.basicHeaderObj.errorCount + self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime + self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs + + def set_RCH(self, expType=2, nTx=1,ipp=None, txA=0, txB=0, + nWindows=None, nHeights=None, firstHeight=None, deltaHeight=None, + numTaus=0, line6Function=0, line5Function=0, fClock=None, + prePulseBefore=0, prePulseAfter=0, + codeType=0, nCode=0, nBaud=0, code=None, + flip1=0, flip2=0): + self.radarControllerHeaderObj.expType = expType + self.radarControllerHeaderObj.nTx = nTx + self.radarControllerHeaderObj.ipp = float(ipp) + self.radarControllerHeaderObj.txA = float(txA) + self.radarControllerHeaderObj.txB = float(txB) + self.radarControllerHeaderObj.rangeIPP = ipp + self.radarControllerHeaderObj.rangeTxA = txA + self.radarControllerHeaderObj.rangeTxB = txB + + self.radarControllerHeaderObj.nHeights = int(nHeights) + self.radarControllerHeaderObj.firstHeight = numpy.array([firstHeight]) + self.radarControllerHeaderObj.deltaHeight = numpy.array([deltaHeight]) + self.radarControllerHeaderObj.samplesWin = numpy.array([nHeights]) + + + self.radarControllerHeaderObj.nWindows = nWindows + self.radarControllerHeaderObj.numTaus = numTaus + self.radarControllerHeaderObj.codeType = codeType + self.radarControllerHeaderObj.line6Function = line6Function + self.radarControllerHeaderObj.line5Function = line5Function + self.radarControllerHeaderObj.fclock = fClock + self.radarControllerHeaderObj.prePulseBefore= prePulseBefore + self.radarControllerHeaderObj.prePulseAfter = prePulseAfter + + self.radarControllerHeaderObj.nCode = nCode + self.radarControllerHeaderObj.nBaud = nBaud + self.radarControllerHeaderObj.code = code + self.radarControllerHeaderObj.flip1 = flip1 + self.radarControllerHeaderObj.flip2 = flip2 + + self.radarControllerHeaderObj.code_size = int(numpy.ceil(nBaud / 32.)) * nCode * 4 + + if fClock is None and deltaHeight is not None: + self.fClock = 0.15 / (deltaHeight * 1e-6) + + def set_PH(self, dtype=0, blockSize=0, profilesPerBlock=0, + dataBlocksPerFile=0, nWindows=0, processFlags=0, nCohInt=0, + nIncohInt=0, totalSpectra=0, nHeights=0, firstHeight=0, + deltaHeight=0, samplesWin=0, spectraComb=0, nCode=0, + code=0, nBaud=None, shif_fft=False, flag_dc=False, + flag_cspc=False, flag_decode=False, flag_deflip=False): + + self.processingHeaderObj.profilesPerBlock = profilesPerBlock + self.processingHeaderObj.dataBlocksPerFile = dataBlocksPerFile + self.processingHeaderObj.nWindows = nWindows + self.processingHeaderObj.nCohInt = nCohInt + self.processingHeaderObj.nIncohInt = nIncohInt + self.processingHeaderObj.totalSpectra = totalSpectra + self.processingHeaderObj.nHeights = int(nHeights) + self.processingHeaderObj.firstHeight = firstHeight + self.processingHeaderObj.deltaHeight = deltaHeight + self.processingHeaderObj.samplesWin = nHeights + + def set_BH(self, utc = 0, miliSecond = 0, timeZone = 0): + self.basicHeaderObj.utc = utc + self.basicHeaderObj.miliSecond = miliSecond + self.basicHeaderObj.timeZone = timeZone + + def set_SH(self, nSamples=0, nProfiles=0, nChannels=0, adcResolution=14, pciDioBusWidth=0): + self.systemHeaderObj.nSamples = nSamples + self.systemHeaderObj.nProfiles = nProfiles + self.systemHeaderObj.nChannels = nChannels + self.systemHeaderObj.adcResolution = adcResolution + self.systemHeaderObj.pciDioBusWidth = pciDioBusWidth + + def init_acquisition(self): + + if self.nFFTPoints != 0: + self.incIntFactor = m_nProfilesperBlock/self.nFFTPoints + if (self.FixPP_IncInt > self.incIntFactor): + self.incIntFactor = self.FixPP_IncInt/ self.incIntFactor + elif(self.FixPP_IncInt< self.incIntFactor): + print("False alert...") + + ProfilesperBlock = self.processingHeaderObj.profilesPerBlock + + self.timeperblock =int(((self.FixRCP_IPP + *ProfilesperBlock + *self.FixPP_CohInt + *self.incIntFactor) + /150.0) + *0.9 + +0.5) + # para cada canal + self.profiles = ProfilesperBlock*self.FixPP_CohInt + self.profiles = ProfilesperBlock + self.Reference = int((self.Tau_0-self.AcqH0_0)/(self.AcqDH_0)+0.5) + self.BaudWidth = int((self.FixRCP_TXA/self.AcqDH_0)/self.Bauds + 0.5 ) + + if (self.BaudWidth==0): + self.BaudWidth=1 + + def init_pulse(self,Num_Codes=Num_Codes,Bauds=Bauds,BaudWidth=BaudWidth,Dyn_snCode=Dyn_snCode): + + Num_Codes = Num_Codes + Bauds = Bauds + BaudWidth = BaudWidth + Dyn_snCode = Dyn_snCode + + if Dyn_snCode: + print("EXISTE") + else: + print("No existe") + + if Dyn_snCode: # if Bauds: + pulses = list(range(0,Num_Codes)) + num_codes = Num_Codes + for i in range(num_codes): + pulse_size = Bauds*BaudWidth + pulses[i] = numpy.zeros(pulse_size) + for j in range(Bauds): + for k in range(BaudWidth): + pulses[i][j*BaudWidth+k] = int(Dyn_snCode[i][j]*600) + else: + print("sin code") + pulses = list(range(1)) + if self.AcqDH_0>0.149: + pulse_size = int(self.FixRCP_TXB/0.15+0.5) + else: + pulse_size = int((self.FixRCP_TXB/self.AcqDH_0)+0.5) #0.0375 + pulses[0] = numpy.ones(pulse_size) + pulses = 600*pulses[0] + + return pulses,pulse_size + + def jro_GenerateBlockOfData(self,Samples=Samples,DC_level= DC_level,stdev=stdev, + Reference= Reference,pulses= pulses, + Num_Codes= Num_Codes,pulse_size=pulse_size, + prof_gen= prof_gen,H0 = H0,DH0=DH0, + Adoppler=Adoppler,Fdoppler= Fdoppler,Hdoppler=Hdoppler): + Samples = Samples + DC_level = DC_level + stdev = stdev + m_nR = Reference + pulses = pulses + num_codes = Num_Codes + ps = pulse_size + prof_gen = prof_gen + channels = self.channels + H0 = H0 + DH0 = DH0 + ippSec = self.radarControllerHeaderObj.ippSeconds + Fdoppler = self.Fdoppler + Hdoppler = self.Hdoppler + Adoppler = self.Adoppler + + self.datablock = numpy.zeros([channels,prof_gen,Samples],dtype= numpy.complex64) + for i in range(channels): + for k in range(prof_gen): + #·······················NOISE··············· + Noise_r = numpy.random.normal(DC_level,stdev,Samples) + Noise_i = numpy.random.normal(DC_level,stdev,Samples) + Noise = numpy.zeros(Samples,dtype=complex) + Noise.real = Noise_r + Noise.imag = Noise_i + #·······················PULSOS·············· + Pulso = numpy.zeros(pulse_size,dtype=complex) + Pulso.real = pulses[k%num_codes] + Pulso.imag = pulses[k%num_codes] + #····················· PULSES+NOISE·········· + InBuffer = numpy.zeros(Samples,dtype=complex) + InBuffer[m_nR:m_nR+ps] = Pulso + InBuffer = InBuffer+Noise + #····················· ANGLE ······························· + InBuffer.real[m_nR:m_nR+ps] = InBuffer.real[m_nR:m_nR+ps]*(math.cos( self.fAngle)*5) + InBuffer.imag[m_nR:m_nR+ps] = InBuffer.imag[m_nR:m_nR+ps]*(math.sin( self.fAngle)*5) + InBuffer=InBuffer + self.datablock[i][k]= InBuffer + #plot_cts(InBuffer,H0=H0,DH0=DH0) + #wave_fft(x=InBuffer,plot_show=True) + #time.sleep(1) + #················DOPPLER SIGNAL............................................... + time_vec = numpy.linspace(0,(prof_gen-1)*ippSec,int(prof_gen))+self.nReadBlocks*ippSec*prof_gen+(self.nReadFiles-1)*ippSec*prof_gen + fd = Fdoppler #+(600.0/120)*self.nReadBlocks + d_signal = Adoppler*numpy.array(numpy.exp(1.0j*2.0*math.pi*fd*time_vec),dtype=numpy.complex64) + #·················· DATABLOCK + DOPPLER············........................... + HD=int(Hdoppler/self.AcqDH_0) + for i in range(12): + self.datablock[:,:,HD+i]=self.datablock[:,:,HD+i]+ d_signal # RESULT + ''' + a= numpy.zeros(10) + for i in range(10): + a[i]=i+self.nReadBlocks+20 + for i in a: + self.datablock[0,:,int(i)]=self.datablock[0,:,int(i)]+ d_signal # RESULT + ''' + + def readBlock(self): + + self.jro_GenerateBlockOfData(Samples= self.samples,DC_level=self.DC_level, + stdev=self.stdev,Reference= self.Reference, + pulses = self.pulses,Num_Codes=self.Num_Codes, + pulse_size=self.pulse_size,prof_gen=self.profiles, + H0=self.H0,DH0=self.DH0) + + self.profileIndex = 0 + self.flagIsNewFile = 0 + self.flagIsNewBlock = 1 + self.nTotalBlocks += 1 + self.nReadBlocks += 1 + + return 1 + + + def getData(self): + if self.flagNoMoreFiles: + self.dataOut.flagNodata = True + return 0 + self.flagDiscontinuousBlock = 0 + self.flagIsNewBlock = 0 + if self.__hasNotDataInBuffer(): # aqui es verdad + if not(self.readNextBlock()): # return 1 y por eso el if not salta a getBasic Header + return 0 + self.getFirstHeader() # atributo + + if not self.getByBlock: + self.dataOut.flagDataAsBlock = False + self.dataOut.data = self.datablock[:, self.profileIndex, :] + self.dataOut.profileIndex = self.profileIndex + self.profileIndex += 1 + else: + pass + self.dataOut.flagNoData = False + self.getBasicHeader() + self.dataOut.realtime = self.online + return self.dataOut.data + + + def setup(self,frequency=49.92e6,incIntFactor= 1, nFFTPoints = 0, FixPP_IncInt=1,FixRCP_IPP=1000, + FixPP_CohInt= 1,Tau_0= 250,AcqH0_0 = 70 ,AcqDH_0=1.25, Bauds= 32, + FixRCP_TXA = 40, FixRCP_TXB = 50, fAngle = 2.0*math.pi*(1/16),DC_level= 50, + stdev= 8,Num_Codes = 1 , Dyn_snCode = None, samples=200, + channels=2,Fdoppler=20,Hdoppler=36,Adoppler=500,nTotalReadFiles=10000, + **kwargs): + + self.set_kwargs(**kwargs) + self.nReadBlocks = 0 + self.nReadFiles = 1 + print('------------------- [Opening file: ] ------------------------------',self.nReadFiles) + + tmp = time.time() + tmp_utc = int(tmp) + tmp_milisecond = int((tmp-tmp_utc)*1000) + print(" SETUP -basicHeaderObj.utc",datetime.datetime.utcfromtimestamp(tmp)) + if Dyn_snCode is None: + Num_Codes=1 + Bauds =1 + + + + self.set_BH(utc= tmp_utc,miliSecond= tmp_milisecond,timeZone=300 ) + + self.set_RCH( expType=0, nTx=150,ipp=FixRCP_IPP, txA=FixRCP_TXA, txB= FixRCP_TXB, + nWindows=1 , nHeights=samples, firstHeight=AcqH0_0, deltaHeight=AcqDH_0, + numTaus=1, line6Function=0, line5Function=0, fClock=None, + prePulseBefore=0, prePulseAfter=0, + codeType=14, nCode=Num_Codes, nBaud=32, code=Dyn_snCode, + flip1=0, flip2=0) + + self.set_PH(dtype=0, blockSize=0, profilesPerBlock=300, + dataBlocksPerFile=120, nWindows=1, processFlags=0, nCohInt=1, + nIncohInt=1, totalSpectra=0, nHeights=samples, firstHeight=AcqH0_0, + deltaHeight=AcqDH_0, samplesWin=samples, spectraComb=0, nCode=0, + code=0, nBaud=None, shif_fft=False, flag_dc=False, + flag_cspc=False, flag_decode=False, flag_deflip=False) + + self.set_SH(nSamples=samples, nProfiles=300, nChannels=channels) + + + self.frequency = frequency + self.incIntFactor = incIntFactor + self.nFFTPoints = nFFTPoints + self.FixPP_IncInt = FixPP_IncInt + self.FixRCP_IPP = FixRCP_IPP + self.FixPP_CohInt = FixPP_CohInt + self.Tau_0 = Tau_0 + self.AcqH0_0 = AcqH0_0 + self.H0 = AcqH0_0 + self.AcqDH_0 = AcqDH_0 + self.DH0 = AcqDH_0 + self.Bauds = Bauds + self.FixRCP_TXA = FixRCP_TXA + self.FixRCP_TXB = FixRCP_TXB + self.fAngle = fAngle + self.DC_level = DC_level + self.stdev = stdev + self.Num_Codes = Num_Codes + self.Dyn_snCode = Dyn_snCode + self.samples = samples + self.channels = channels + self.profiles = None + self.m_nReference = None + self.Baudwidth = None + self.Fdoppler = Fdoppler + self.Hdoppler = Hdoppler + self.Adoppler = Adoppler + self.nTotalReadFiles = int(nTotalReadFiles) + + print("IPP ", self.FixRCP_IPP) + print("Tau_0 ",self.Tau_0) + print("AcqH0_0",self.AcqH0_0) + print("samples,window ",self.samples) + print("AcqDH_0",AcqDH_0) + print("FixRCP_TXA",self.FixRCP_TXA) + print("FixRCP_TXB",self.FixRCP_TXB) + print("Dyn_snCode",Dyn_snCode) + print("Fdoppler", Fdoppler) + print("Hdoppler",Hdoppler) + print("Vdopplermax",Fdoppler*(3.0e8/self.frequency)/2.0) + print("nTotalReadFiles", nTotalReadFiles) + + self.init_acquisition() + self.pulses,self.pulse_size=self.init_pulse(Num_Codes=self.Num_Codes,Bauds=self.Bauds,BaudWidth=self.BaudWidth,Dyn_snCode=Dyn_snCode) + print(" [ END ] - SETUP metodo") + return + + def run(self,**kwargs): # metodo propio + if not(self.isConfig): + self.setup(**kwargs) + self.isConfig = True + self.getData() diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index 6ce7724..8e99b76 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -1,5 +1,5 @@ import sys -import numpy +import numpy,math from scipy import interpolate from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator from schainpy.model.data.jrodata import Voltage @@ -8,8 +8,8 @@ from time import time -class VoltageProc(ProcessingUnit): - +class VoltageProc(ProcessingUnit): + def __init__(self): ProcessingUnit.__init__(self) @@ -51,22 +51,20 @@ class VoltageProc(ProcessingUnit): self.dataOut.beam.codeList = self.dataIn.beam.codeList self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList self.dataOut.beam.zenithList = self.dataIn.beam.zenithList - + class selectChannels(Operation): - + def run(self, dataOut, channelList): channelIndexList = [] self.dataOut = dataOut - for channel in channelList: if channel not in self.dataOut.channelList: raise ValueError("Channel %d is not in %s" %(channel, str(self.dataOut.channelList))) index = self.dataOut.channelList.index(channel) channelIndexList.append(index) - self.selectChannelsByIndex(channelIndexList) return self.dataOut @@ -105,6 +103,7 @@ class selectChannels(Operation): self.dataOut.data = data # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList] self.dataOut.channelList = range(len(channelIndexList)) + elif self.dataOut.type == 'Spectra': data_spc = self.dataOut.data_spc[channelIndexList, :] data_dc = self.dataOut.data_dc[channelIndexList, :] @@ -146,7 +145,7 @@ class selectChannels(Operation): return class selectHeights(Operation): - + def run(self, dataOut, minHei=None, maxHei=None): """ Selecciona un bloque de datos en base a un grupo de valores de alturas segun el rango @@ -262,7 +261,7 @@ class selectHeights(Operation): self.dataOut.data_dc = data_dc self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1] - + return 1 @@ -286,7 +285,7 @@ class filterByHeights(Operation): """ Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis] """ - buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)] + buffer = dataOut.data[:, :, 0:int(dataOut.nHeights-r)] buffer = buffer.reshape(dataOut.nChannels, dataOut.nProfiles, int(dataOut.nHeights/window), window) buffer = numpy.sum(buffer,3) @@ -580,8 +579,8 @@ class CohInt(Operation): # print self.__bufferStride[self.__profIndexStride - 1] # raise return self.__bufferStride[self.__profIndexStride - 1] - - + + return None, None def integrate(self, data, datatime=None): @@ -603,7 +602,7 @@ class CohInt(Operation): avgdatatime = self.__initime deltatime = datatime - self.__lastdatatime - + if not self.__withOverlapping: self.__initime = datatime else: @@ -629,7 +628,7 @@ class CohInt(Operation): avgdatatime = (times - 1) * timeInterval + dataOut.utctime self.__dataReady = True return avgdata, avgdatatime - + def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs): if not self.isConfig: @@ -643,12 +642,12 @@ class CohInt(Operation): avgdata, avgdatatime = self.integrateByBlock(dataOut) dataOut.nProfiles /= self.n else: - if stride is None: + if stride is None: avgdata, avgdatatime = self.integrate(dataOut.data, dataOut.utctime) else: avgdata, avgdatatime = self.integrateByStride(dataOut.data, dataOut.utctime) - + # dataOut.timeInterval *= n dataOut.flagNoData = True @@ -753,11 +752,11 @@ class Decoder(Operation): junk = junk.flatten() code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud)) profilesList = range(self.__nProfiles) - - for i in range(self.__nChannels): - for j in profilesList: - self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:] - return self.datadecTime + + for i in range(self.__nChannels): + for j in profilesList: + self.datadecTime[i,j,:] = numpy.correlate(data[i,j,:], code_block[j,:], mode='full')[self.nBaud-1:] + return self.datadecTime def __convolutionByBlockInFreq(self, data): @@ -774,7 +773,7 @@ class Decoder(Operation): return data - + def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None): if dataOut.flagDecodeData: @@ -805,7 +804,7 @@ class Decoder(Operation): self.__nProfiles = dataOut.nProfiles datadec = None - + if mode == 3: mode = 0 @@ -1186,9 +1185,9 @@ class SplitProfiles(Operation): if shape[2] % n != 0: raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[2])) - + new_shape = shape[0], shape[1]*n, int(shape[2]/n) - + dataOut.data = numpy.reshape(dataOut.data, new_shape) dataOut.flagNoData = False @@ -1272,6 +1271,155 @@ class CombineProfiles(Operation): dataOut.ippSeconds *= n return dataOut + +class PulsePairVoltage(Operation): + ''' + Function PulsePair(Signal Power, Velocity) + The real component of Lag[0] provides Intensity Information + The imag component of Lag[1] Phase provides Velocity Information + + Configuration Parameters: + nPRF = Number of Several PRF + theta = Degree Azimuth angel Boundaries + + Input: + self.dataOut + lag[N] + Affected: + self.dataOut.spc + ''' + isConfig = False + __profIndex = 0 + __initime = None + __lastdatatime = None + __buffer = None + __buffer2 = [] + __buffer3 = None + __dataReady = False + n = None + __nch = 0 + __nHeis = 0 + removeDC = False + ipp = None + lambda_ = 0 + + def __init__(self,**kwargs): + Operation.__init__(self,**kwargs) + + def setup(self, dataOut, n = None, removeDC=False): + ''' + n= Numero de PRF's de entrada + ''' + self.__initime = None + self.__lastdatatime = 0 + self.__dataReady = False + self.__buffer = 0 + self.__buffer2 = [] + self.__buffer3 = 0 + self.__profIndex = 0 + + self.__nch = dataOut.nChannels + self.__nHeis = dataOut.nHeights + self.removeDC = removeDC + self.lambda_ = 3.0e8/(9345.0e6) + self.ippSec = dataOut.ippSeconds + self.nCohInt = dataOut.nCohInt + print("IPPseconds",dataOut.ippSeconds) + + print("ELVALOR DE n es:", n) + if n == None: + raise ValueError("n should be specified.") + + if n != None: + if n<2: + raise ValueError("n should be greater than 2") + + self.n = n + self.__nProf = n + + self.__buffer = numpy.zeros((dataOut.nChannels, + n, + dataOut.nHeights), + dtype='complex') + + def putData(self,data): + ''' + Add a profile to he __buffer and increase in one the __profiel Index + ''' + self.__buffer[:,self.__profIndex,:]= data + self.__profIndex += 1 + return + + def pushData(self): + ''' + Return the PULSEPAIR and the profiles used in the operation + Affected : self.__profileIndex + ''' + + if self.removeDC==True: + mean = numpy.mean(self.__buffer,1) + tmp = mean.reshape(self.__nch,1,self.__nHeis) + dc= numpy.tile(tmp,[1,self.__nProf,1]) + self.__buffer = self.__buffer - dc + + data_intensity = numpy.sum(self.__buffer*numpy.conj(self.__buffer),1)/(self.n*self.nCohInt)#*self.nCohInt) + pair1 = self.__buffer[:,1:,:]*numpy.conjugate(self.__buffer[:,:-1,:]) + angle = numpy.angle(numpy.sum(pair1,1))*180/(math.pi) + #print(angle.shape)#print("__ANGLE__") #print("angle",angle[:,:10]) + data_velocity = (self.lambda_/(4*math.pi*self.ippSec))*numpy.angle(numpy.sum(pair1,1))#self.ippSec*self.nCohInt + n = self.__profIndex + + self.__buffer = numpy.zeros((self.__nch, self.__nProf,self.__nHeis), dtype='complex') + self.__profIndex = 0 + return data_intensity,data_velocity,n + + def pulsePairbyProfiles(self,data): + + self.__dataReady = False + data_intensity = None + data_velocity = None + self.putData(data) + if self.__profIndex == self.n: + data_intensity, data_velocity, n = self.pushData() + self.__dataReady = True + + return data_intensity, data_velocity + + def pulsePairOp(self, data, datatime= None): + + if self.__initime == None: + self.__initime = datatime + + data_intensity, data_velocity = self.pulsePairbyProfiles(data) + self.__lastdatatime = datatime + + if data_intensity is None: + return None, None, None + + avgdatatime = self.__initime + deltatime = datatime - self.__lastdatatime + self.__initime = datatime + + return data_intensity, data_velocity, avgdatatime + + def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs): + + if not self.isConfig: + self.setup(dataOut = dataOut, n = n , removeDC=removeDC , **kwargs) + self.isConfig = True + data_intensity, data_velocity, avgdatatime = self.pulsePairOp(dataOut.data, dataOut.utctime) + dataOut.flagNoData = True + + if self.__dataReady: + dataOut.nCohInt *= self.n + dataOut.data_intensity = data_intensity #valor para intensidad + dataOut.data_velocity = data_velocity #valor para velocidad + dataOut.PRFbyAngle = self.n #numero de PRF*cada angulo rotado que equivale a un tiempo. + dataOut.utctime = avgdatatime + dataOut.flagNoData = False + return dataOut + + # import collections # from scipy.stats import mode #