schain Commit - r1311:816ce18b0db1 · Jicamarca Repository

jrodata se cambio los atrituos de pulsepair,jroplot_voltage se cambio el ploteo de la potencial y la senal, jroIO_simulator se anadio la modificacion de bloquesporarchivo y perfilesporbloque,jro_proc_voltage se corrigio el calculo de noise con remocion DC y se asignaron nuevos nombres al dataOut, jroproc_parameters se anadio los nuevos atributos del dataOut jroproc_voltage pulsepair,test_sim0009.py es el nuevo test de escritura

avaldez -

r1311:816ce18b0db1 v3-devel-julio

parent child

Context file:

r1311:816ce18b0db1

Collapse all files

The requested changes are too big and content was truncated. Show full diff

schainpy/model/data/jrodata.py +42 -7

@@ -359,10 +359,12 class JROData(GenericData):
359	class Voltage(JROData):	359	class Voltage(JROData):
360		360
361	# data es un numpy array de 2 dmensiones (canales, alturas)	361	# data es un numpy array de 2 dmensiones (canales, alturas)
362	data = None	362	data = None
363	data~~_intensity~~ = None	363	dataPP_POW = None
364	data~~_velocity~~ = None	364	dataPP_DOP = None
365	data~~_specwidth~~ = None	365	dataPP_WIDTH = None
		366	dataPP_SNR = None
		367
366	def __init__(self):	368	def __init__(self):
367	'''	369	'''
368	Constructor	370	Constructor
@@ -421,6 +423,35 class Voltage(JROData):
421		423
422	return noise	424	return noise
423		425
		426	def getNoisebyHildebrandDC(self, channel=None,DC=0):
		427	"""
		428	Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
		429
		430	Return:
		431	noiselevel
		432	"""
		433
		434	if channel != None:
		435	data = self.data[channel]-DC
		436	nChannels = 1
		437	else:
		438	data = self.data
		439	nChannels = self.nChannels
		440
		441	noise = numpy.zeros(nChannels)
		442	power = data * numpy.conjugate(data)
		443
		444	for thisChannel in range(nChannels):
		445	if nChannels == 1:
		446	daux = power[:].real
		447	else:
		448	daux = power[thisChannel, :].real
		449	noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
		450
		451	return noise
		452
		453
		454
424	def getNoise(self, type=1, channel=None):	455	def getNoise(self, type=1, channel=None):
425		456
426	if type == 1:	457	if type == 1:
@@ -1262,15 +1293,19 class PlotterData(object):
1262	self.flagDataAsBlock = dataOut.flagDataAsBlock	1293	self.flagDataAsBlock = dataOut.flagDataAsBlock
1263	self.nProfiles = dataOut.nProfiles	1294	self.nProfiles = dataOut.nProfiles
1264	if plot == 'pp_power':	1295	if plot == 'pp_power':
1265	buffer = dataOut.data~~_intensity~~	1296	buffer = dataOut.dataPP_POWER
		1297	self.flagDataAsBlock = dataOut.flagDataAsBlock
		1298	self.nProfiles = dataOut.nProfiles
		1299	if plot == 'pp_signal':
		1300	buffer = dataOut.dataPP_POW
1266	self.flagDataAsBlock = dataOut.flagDataAsBlock	1301	self.flagDataAsBlock = dataOut.flagDataAsBlock
1267	self.nProfiles = dataOut.nProfiles	1302	self.nProfiles = dataOut.nProfiles
1268	if plot == 'pp_velocity':	1303	if plot == 'pp_velocity':
1269	buffer = dataOut.data~~_velocity~~	1304	buffer = dataOut.dataPP_DOP
1270	self.flagDataAsBlock = dataOut.flagDataAsBlock	1305	self.flagDataAsBlock = dataOut.flagDataAsBlock
1271	self.nProfiles = dataOut.nProfiles	1306	self.nProfiles = dataOut.nProfiles
1272	if plot == 'pp_specwidth':	1307	if plot == 'pp_specwidth':
1273	buffer = dataOut.data~~_specwidth~~	1308	buffer = dataOut.dataPP_WIDTH
1274	self.flagDataAsBlock = dataOut.flagDataAsBlock	1309	self.flagDataAsBlock = dataOut.flagDataAsBlock
1275	self.nProfiles = dataOut.nProfiles	1310	self.nProfiles = dataOut.nProfiles
1276		1311

schainpy/model/graphics/jroplot_voltage.py +29 -3

                     thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
                     if self.CODE == "pp_power":
                         scope = self.data['pp_power']
+                    elif self.CODE == "pp_signal":
+                        scope = self.data["pp_signal"]
                     elif self.CODE == "pp_velocity":
                         scope = self.data["pp_velocity"]
                     elif self.CODE == "pp_specwidth":
                                            thisDatetime,
                                            wintitle
                                            )
+                            if self.CODE=="pp_signal":
+                                self.plot_weatherpower(self.data.heights,
+                                           scope[:,i,:],
+                                           channels,
+                                           thisDatetime,
+                                           wintitle
+                                           )
                             if self.CODE=="pp_velocity":
                                 self.plot_weathervelocity(scope[:,i,:],
                                            self.data.heights,
                                                 thisDatetime,
                                                 wintitle
                                                    )
+                        if self.CODE=="pp_signal":
+                            self.plot_weatherpower(self.data.heights,
+                                                scope,
+                                                channels,
+                                                thisDatetime,
+                                                wintitle
+                                                   )
                         if self.CODE=="pp_velocity":
                             self.plot_weathervelocity(scope,
                                                    self.data.heights,
             class PulsepairPowerPlot(ScopePlot):
                 '''
-                Plot for
+                Plot for  P= S+N
                 '''
                 CODE = 'pp_power'
             class PulsepairVelocityPlot(ScopePlot):
                 '''
-                Plot for
+                Plot for VELOCITY
                 '''
                 CODE = 'pp_velocity'
                 plot_name = 'PulsepairVelocity'
             class PulsepairSpecwidthPlot(ScopePlot):
                 '''
-                Plot for
+                Plot for WIDTH
                 '''
                 CODE = 'pp_specwidth'
                 plot_name = 'PulsepairSpecwidth'
                 plot_type = 'scatter'
                 buffering = False
+            class PulsepairSignalPlot(ScopePlot):
+                '''
+                Plot for S
+                '''
+                CODE = 'pp_signal'
+                plot_name = 'PulsepairSignal'
+                plot_type = 'scatter'
+                buffering = False

schainpy/model/io/jroIO_simulator.py +21 -14

                     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)
                     #·············Señal con ancho espectral····················
-                    #specw_sig  = numpy.linspace(-149,150,300)
+                    if prof_gen%2==0:
-                    #w          = 8
+                        min = int(prof_gen/2.0-1.0)
-                    #A          = 20
+                        max = int(prof_gen/2.0)
-                    #specw_sig   = specw_sig/w
+                    else:
-                    #specw_sig   = numpy.sinc(specw_sig)
+                        min = int(prof_gen/2.0)
-                    #specw_sig   =  A*numpy.array(specw_sig,dtype=numpy.complex64)
+                        max = int(prof_gen/2.0)
+                    specw_sig  = numpy.linspace(-min,max,prof_gen)
+                    w          = 4
+                    A          = 20
+                    specw_sig   = specw_sig/w
+                    specw_sig   = numpy.sinc(specw_sig)
+                    specw_sig   =  A*numpy.array(specw_sig,dtype=numpy.complex64)
                     #·················· DATABLOCK + DOPPLER····················
                     HD=int(Hdoppler/self.AcqDH_0)
                     for  i in range(12):
                         self.datablock[0,:,HD+i]=self.datablock[0,:,HD+i]+ d_signal# RESULT
                     #·················· DATABLOCK + DOPPLER*Sinc(x)····················
-                    #HD=int(Hdoppler/self.AcqDH_0)
+                    HD=int(Hdoppler/self.AcqDH_0)
-                    #HD=int(HD/2)
+                    HD=int(HD/2)
-                    #for  i in range(12):
+                    for  i in range(12):
-                    #    self.datablock[0,:,HD+i]=self.datablock[0,:,HD+i]+ specw_sig*d_signal# RESULT
+                        self.datablock[0,:,HD+i]=self.datablock[0,:,HD+i]+ specw_sig*d_signal# RESULT
                 def readBlock(self):
                                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,
+                               channels=2,Fdoppler=20,Hdoppler=36,Adoppler=500,
+                               profilesPerBlock=300,dataBlocksPerFile=120,nTotalReadFiles=10000,
                                **kwargs):
                     self.set_kwargs(**kwargs)
                              codeType=0, nCode=Num_Codes, nBaud=32, code=Dyn_snCode,
                              flip1=0, flip2=0,Taus=Tau_0)
-                    self.set_PH(dtype=0, blockSize=0, profilesPerBlock=300,
+                    self.set_PH(dtype=0, blockSize=0, profilesPerBlock=profilesPerBlock,
-                                  dataBlocksPerFile=120, nWindows=1, processFlags=numpy.array([1024]), nCohInt=1,
+                                  dataBlocksPerFile=dataBlocksPerFile, nWindows=1, processFlags=numpy.array([1024]), 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.set_SH(nSamples=samples, nProfiles=profilesPerBlock, nChannels=channels)
                     self.readFirstHeader()

schainpy/model/proc/jroproc_parameters.py 0 0

	1	NO CONTENT: modified file			NO CONTENT: modified file
The requested commit or file is too big and content was truncated. Show full diff

schainpy/model/proc/jroproc_voltage.py +57 -36

             import numpy,math
             from scipy import interpolate
             from schainpy.model.proc.jroproc_base import ProcessingUnit, Operation, MPDecorator
-            from schainpy.model.data.jrodata import Voltage
+            from schainpy.model.data.jrodata import Voltage,hildebrand_sekhon
+            from schainpy.model.data import _noise
             from schainpy.utils import log
             from time import time
                                                        n,
                                                        dataOut.nHeights),
                                                       dtype='complex')
-                    #self.noise  = numpy.zeros([self.__nch,self.__nHeis])
-                    #for i in range(self.__nch):
-                    #    self.noise[i]=dataOut.getNoise(channel=i)
                 def putData(self,data):
                     '''
                     Return the PULSEPAIR and the profiles used in the operation
                     Affected :  self.__profileIndex
                     '''
+                    #················· Remove DC···································
                     if self.removeDC==True:
                         mean    = numpy.mean(self.__buffer,1)
                         tmp     = mean.reshape(self.__nch,1,self.__nHeis)
                         dc= numpy.tile(tmp,[1,self.__nProf,1])
                         self.__buffer = self.__buffer -  dc
+                    #··················Calculo de Potencia ························
+                    pair0       = self.__buffer*numpy.conj(self.__buffer)
+                    pair0       = pair0.real
+                    lag_0       = numpy.sum(pair0,1)
+                    #··················Calculo de Ruido x canal····················
+                    self.noise  = numpy.zeros(self.__nch)
+                    for i in range(self.__nch):
+                        daux         = numpy.sort(pair0[i,:,:],axis= None)
+                        self.noise[i]=hildebrand_sekhon( daux ,self.nCohInt)
+                    self.noise       = self.noise.reshape(self.__nch,1)
+                    self.noise       = numpy.tile(self.noise,[1,self.__nHeis])
+                    noise_buffer     = self.noise.reshape(self.__nch,1,self.__nHeis)
+                    noise_buffer     = numpy.tile(noise_buffer,[1,self.__nProf,1])
+                    #·················· Potencia recibida= P , Potencia senal = S , Ruido= N··
+                    #··················   P= S+N  ,P=lag_0/N ·································
+                    #···················· Power ··················································
+                    data_power       = lag_0/(self.n*self.nCohInt)
+                    #------------------  Senal  ···················································
+                    data_intensity   = pair0 - noise_buffer
+                    data_intensity   = numpy.sum(data_intensity,axis=1)*(self.n*self.nCohInt)#*self.nCohInt)
+                    #data_intensity   = (lag_0-self.noise*self.n)*(self.n*self.nCohInt)
+                    for i in range(self.__nch):
+                        for j in range(self.__nHeis):
+                            if data_intensity[i][j]  < 0:
+                                data_intensity[i][j] = numpy.min(numpy.absolute(data_intensity[i][j]))
-                    lag_0            = numpy.sum(self.__buffer*numpy.conj(self.__buffer),1)
+                    #················· Calculo de Frecuencia y Velocidad doppler········
-                    data_intensity   = lag_0/(self.n*self.nCohInt)#*self.nCohInt)
                     pair1            = self.__buffer[:,:-1,:]*numpy.conjugate(self.__buffer[:,1:,:])
                     lag_1            = numpy.sum(pair1,1)
-                    #angle            = numpy.angle(numpy.sum(pair1,1))*180/(math.pi)
+                    data_freq        = (-1/(2.0*math.pi*self.ippSec*self.nCohInt))*numpy.angle(lag_1)
-                    data_velocity    = (-1.0*self.lambda_/(4*math.pi*self.ippSec))*numpy.angle(lag_1)#self.ippSec*self.nCohInt
+                    data_velocity    = (self.lambda_/2.0)*data_freq
-                    self.noise  = numpy.zeros([self.__nch,self.__nHeis])
+                    #················ Potencia promedio estimada de la Senal···········
-                    for i in range(self.__nch):
+                    lag_0            = lag_0/self.n
-                        self.noise[i]=dataOut.getNoise(channel=i)
+                    S                = lag_0-self.noise
-                    lag_0            = lag_0.real/(self.n)
+                    #················ Frecuencia Doppler promedio ·····················
                     lag_1            = lag_1/(self.n-1)
                     R1               = numpy.abs(lag_1)
-                    S                = (lag_0-self.noise)
+                    #················ Calculo del SNR··································
                     data_snrPP       = S/self.noise
-                    data_snrPP       = numpy.where(data_snrPP<0,1,data_snrPP)
+                    for i in range(self.__nch):
+                        for j in range(self.__nHeis):
+                            if data_snrPP[i][j]  < 1.e-20:
+                                data_snrPP[i][j] = 1.e-20
+                    #················· Calculo del ancho espectral ······················
                     L                = S/R1
                     L                = numpy.where(L<0,1,L)
                     L                = numpy.log(L)
                     tmp              = numpy.sqrt(numpy.absolute(L))
+                    data_specwidth   = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec*self.nCohInt))*tmp*numpy.sign(L)
-                    data_specwidth   = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec))*tmp*numpy.sign(L)
-                    #data_specwidth   = (self.lambda_/(2*math.sqrt(2)*math.pi*self.ippSec))*k
                     n                = self.__profIndex
                     self.__buffer    = numpy.zeros((self.__nch, self.__nProf,self.__nHeis),  dtype='complex')
                     self.__profIndex = 0
-                    return data_intensity,data_velocity,data_snrPP,data_specwidth,n
+                    return data_power,data_intensity,data_velocity,data_snrPP,data_specwidth,n
                 def pulsePairbyProfiles(self,dataOut):
                     self.__dataReady     =  False
+                    data_power           =  None
                     data_intensity       =  None
                     data_velocity        =  None
                     data_specwidth       =  None
                     data_snrPP           =  None
                     self.putData(data=dataOut.data)
                     if self.__profIndex  == self.n:
-                        #self.noise  = numpy.zeros([self.__nch,self.__nHeis])
+                        data_power,data_intensity, data_velocity,data_snrPP,data_specwidth, n   = self.pushData(dataOut=dataOut)
-                        #for i in range(self.__nch):
-                        #    self.noise[i]=data.getNoise(channel=i)
-                        #print(self.noise.shape)
-                        data_intensity, data_velocity,data_snrPP,data_specwidth, n   = self.pushData(dataOut=dataOut)
                         self.__dataReady                   = True
-                    return data_intensity, data_velocity,data_snrPP,data_specwidth
+                    return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth
                 def pulsePairOp(self, dataOut, datatime= None):
                     if self.__initime == None:
                         self.__initime = datatime
-                    #print("hola")
+                    data_power, data_intensity, data_velocity, data_snrPP, data_specwidth = self.pulsePairbyProfiles(dataOut)
-                    data_intensity, data_velocity,data_snrPP,data_specwidth = self.pulsePairbyProfiles(dataOut)
                     self.__lastdatatime           = datatime
-                    if data_intensity is None:
+                    if data_power is None:
-                        return None, None,None,None,None
+                        return None, None, None,None,None,None
                     avgdatatime    = self.__initime
                     deltatime      = datatime - self.__lastdatatime
                     self.__initime = datatime
-                    return data_intensity, data_velocity,data_snrPP,data_specwidth,avgdatatime
+                    return data_power, data_intensity, data_velocity, data_snrPP, data_specwidth, avgdatatime
                 def run(self, dataOut,n = None,removeDC= False, overlapping= False,**kwargs):
                     if not self.isConfig:
                         self.setup(dataOut = dataOut, n    = n , removeDC=removeDC , **kwargs)
                         self.isConfig   = True
-                    data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
+                    data_power, data_intensity, data_velocity,data_snrPP,data_specwidth, avgdatatime = self.pulsePairOp(dataOut, dataOut.utctime)
                     dataOut.flagNoData                         = True
                     if self.__dataReady:
                         dataOut.nCohInt        *= self.n
-                        dataOut.data_intensity  = data_intensity #valor para intensidad
+                        dataOut.dataPP_POW      = data_intensity # S
-                        dataOut.data_velocity   = data_velocity  #valor para velocidad
+                        dataOut.dataPP_POWER    = data_power     # P
-                        dataOut.data_snrPP      = data_snrPP     # valor para snr
+                        dataOut.dataPP_DOP      = data_velocity
-                        dataOut.data_specwidth  = data_specwidth
+                        dataOut.dataPP_SNR      = data_snrPP
+                        dataOut.dataPP_WIDTH    = data_specwidth
                         dataOut.PRFbyAngle      = self.n         #numero de PRF*cada angulo rotado que equivale a un tiempo.
                         dataOut.utctime         = avgdatatime
                         dataOut.flagNoData      = False

schainpy/scripts/test_sim0009.py +36 -9

@@ -2,8 +2,11 import os,sys
2	import datetime	2	import datetime
3	import time	3	import time
4	from schainpy.controller import Project	4	from schainpy.controller import Project
5	path = '/home/alex/Downloads/NEW_WR2/spc16removeDC'	5	#path = '/home/alex/Downloads/NEW_WR2/spc16removeDC'
6	figpath = path	6	#figpath = path
		7
		8	path = '/home/alex/Downloads/test_rawdata'
		9	figpath = '/home/alex/Downloads/hdf5_testPP'
7	desc = "Simulator Test"	10	desc = "Simulator Test"
8		11
9	controllerObj = Project()	12	controllerObj = Project()
@@ -25,25 +28,49 readUnitConfObj = controllerObj.addReadUnit(datatype='SimulatorReader',
25	delay=0,	28	delay=0,
26	online=0,	29	online=0,
27	walk=0,	30	walk=0,
28	~~nTotalReadFiles~~=3)	31	profilesPerBlock=625,
29		32	dataBlocksPerFile=100)#,#nTotalReadFiles=2)
30	opObj11 = readUnitConfObj.addOperation(name='printInfo')
31		33
		34	'''
		35	readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
		36	path=path,
		37	startDate="2020/01/01", #"2020/01/01",#today,
		38	endDate= "2020/12/01", #"2020/12/30",#today,
		39	startTime='00:00:00',
		40	endTime='23:59:59',
		41	delay=0,
		42	#set=0,
		43	online=0,
		44	walk=1)
		45	'''
		46	#opObj11 = readUnitConfObj.addOperation(name='printInfo')
32	procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())	47	procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
33	#opObj11 = procUnitConfObjA.addOperation(name='CohInt', optype='other')	48	#opObj11 = procUnitConfObjA.addOperation(name='CohInt', optype='other')
34	#opObj11.addParameter(name='n', value='10', format='int')	49	#opObj11.addParameter(name='n', value='4', format='int')
35		50
36	#opObj10 = procUnitConfObjA.addOperation(name='selectChannels')	51	#opObj10 = procUnitConfObjA.addOperation(name='selectChannels')
37	#opObj10.addParameter(name='channelList', value=[0])	52	#opObj10.addParameter(name='channelList', value=[0])
38	opObj11 = procUnitConfObjA.addOperation(name='PulsePairVoltage', optype='other')	53	opObj11 = procUnitConfObjA.addOperation(name='PulsePairVoltage', optype='other')
39	opObj11.addParameter(name='n', value='~~300~~', format='int')#10	54	opObj11.addParameter(name='n', value='625', format='int')#10
40	opObj11.addParameter(name='removeDC', value=1, format='int')	55	opObj11.addParameter(name='removeDC', value=1, format='int')
41		56
42	#opObj11 = procUnitConfObjA.addOperation(name='PulsepairPowerPlot', optype='other')	57	#opObj11 = procUnitConfObjA.addOperation(name='PulsepairPowerPlot', optype='other')
		58	#opObj11 = procUnitConfObjA.addOperation(name='PulsepairSignalPlot', optype='other')
43		59
44	opObj11 = procUnitConfObjA.addOperation(name='PulsepairVelocityPlot', optype='other')	60
		61	#opObj11 = procUnitConfObjA.addOperation(name='PulsepairVelocityPlot', optype='other')
45	#opObj11.addParameter(name='xmax', value=8)	62	#opObj11.addParameter(name='xmax', value=8)
46		63
47	opObj11 = procUnitConfObjA.addOperation(name='PulsepairSpecwidthPlot', optype='other')	64	#opObj11 = procUnitConfObjA.addOperation(name='PulsepairSpecwidthPlot', optype='other')
		65
		66	procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId())
		67
		68
		69	opObj10 = procUnitConfObjB.addOperation(name='ParameterWriter')
		70	opObj10.addParameter(name='path',value=figpath)
		71	#opObj10.addParameter(name='mode',value=0)
		72	opObj10.addParameter(name='blocksPerFile',value='100',format='int')
		73	opObj10.addParameter(name='metadataList',value='utctimeInit,timeInterval',format='list')
		74	opObj10.addParameter(name='dataList',value='dataPP_POW,dataPP_DOP,dataPP_SNR,dataPP_WIDTH')#,format='list'
48		75
49	controllerObj.start()	76	controllerObj.start()

General Comments 0

Write
Preview

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

No TODOs yet

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

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