schain Commit - r1323:133a5771abd8 · Jicamarca Repository

Merge branch 'v3.0-devel' of http://jro-dev.igp.gob.pe/rhodecode/schain into v3.0-devel

Juan C. Espinoza -

r1323:133a5771abd8

parent child

Context file:

r1323:133a5771abd8

Collapse all files

schainpy/scripts/pedestal_client.py created 644 +87 0

@@ -0,0 +1,87
	1	import numpy
	2	import sys
	3	import zmq
	4	import time
	5	import h5py
	6	import os
	7
	8	path="/home/alex/Downloads/pedestal"
	9	ext=".hdf5"
	10
	11	port ="5556"
	12	if len(sys.argv)>1:
	13	port = sys.argv[1]
	14	int(port)
	15
	16	if len(sys.argv)>2:
	17	port1 = sys.argv[2]
	18	int(port1)
	19
	20	#Socket to talk to server
	21	context = zmq.Context()
	22	socket = context.socket(zmq.SUB)
	23
	24	print("Collecting updates from weather server...")
	25	socket.connect("tcp://localhost:%s"%port)
	26
	27	if len(sys.argv)>2:
	28	socket.connect("tcp://localhost:%s"%port1)
	29
	30	#Subscribe to zipcode, default is NYC,10001
	31	topicfilter = "10001"
	32	socket.setsockopt_string(zmq.SUBSCRIBE,topicfilter)
	33	#Process 5 updates
	34	total_value=0
	35	count= -1
	36	azi= []
	37	elev=[]
	38	time0=[]
	39	#for update_nbr in range(250):
	40	while(True):
	41	string= socket.recv()
	42	topic,ang_elev,ang_elev_dec,ang_azi,ang_azi_dec,seconds,seconds_dec= string.split()
	43	ang_azi =float(ang_azi)+1e-3*float(ang_azi_dec)
	44	ang_elev =float(ang_elev)+1e-3*float(ang_elev_dec)
	45	seconds =float(seconds) +1e-6*float(seconds_dec)
	46	azi.append(ang_azi)
	47	elev.append(ang_elev)
	48	time0.append(seconds)
	49	count +=1
	50	if count == 100:
	51	timetuple=time.localtime()
	52	epoc = time.mktime(timetuple)
	53	#print(epoc)
	54	fullpath = path + ("/" if path[-1]!="/" else "")
	55
	56	if not os.path.exists(fullpath):
	57	os.mkdir(fullpath)
	58
	59	azi_array = numpy.array(azi)
	60	elev_array = numpy.array(elev)
	61	time0_array= numpy.array(time0)
	62	pedestal_array=numpy.array([azi,elev,time0])
	63	count=0
	64	azi= []
	65	elev=[]
	66	time0=[]
	67	#print(pedestal_array[0])
	68	#print(pedestal_array[1])
	69
	70	meta='PE'
	71	filex="%s%4.4d%3.3d%10.4d%s"%(meta,timetuple.tm_year,timetuple.tm_yday,epoc,ext)
	72	filename = os.path.join(fullpath,filex)
	73	fp = h5py.File(filename,'w')
	74	#print("Escribiendo HDF5...",epoc)
	75	#·················· Data·....······································
	76	grp = fp.create_group("Data")
	77	dset = grp.create_dataset("azimuth" , data=pedestal_array[0])
	78	dset = grp.create_dataset("elevacion", data=pedestal_array[1])
	79	dset = grp.create_dataset("utc" , data=pedestal_array[2])
	80	#·················· Metadata·······································
	81	grp = fp.create_group("Metadata")
	82	dset = grp.create_dataset("utctimeInit", data=pedestal_array[2][0])
	83	timeInterval = pedestal_array[2][1]-pedestal_array[2][0]
	84	dset = grp.create_dataset("timeInterval", data=timeInterval)
	85	fp.close()
	86
	87	#print ("Average messagedata value for topic '%s' was %dF" % ( topicfilter,total_value / update_nbr))

schainpy/scripts/pedestal_server.py created 644 +48 0

@@ -0,0 +1,48
	1	###########################################################################
	2	############################### SERVIDOR###################################
	3	######################### SIMULADOR DE PEDESTAL############################
	4	###########################################################################
	5	import time
	6	import math
	7	import numpy
	8	import struct
	9	from time import sleep
	10	import zmq
	11	import pickle
	12	port="5556"
	13	context = zmq.Context()
	14	socket = context.socket(zmq.PUB)
	15	socket.bind("tcp://*:%s"%port)
	16	###### PARAMETROS DE ENTRADA################################
	17	print("PEDESTAL RESOLUCION 0.01")
	18	print("MAXIMA VELOCIDAD DEL PEDESTAL")
	19	ang_elev = 4.12
	20	ang_azi = 30
	21	velocidad= input ("Ingresa velocidad:")
	22	velocidad= float(velocidad)
	23	print (velocidad)
	24	############################################################
	25	sleep(3)
	26	print("Start program")
	27	t1 = time.time()
	28	count=0
	29	while(True):
	30	tmp_vuelta = int(360/velocidad)
	31	t1=t1+tmp_vuelta*count
	32	count= count+1
	33	muestras_seg = 100
	34	t2 = time.time()
	35	for i in range(tmp_vuelta):
	36	for j in range(muestras_seg):
	37	tmp_variable = (i+j/100.0)
	38	ang_azi = (tmp_variable)*float(velocidad)
	39	seconds = t1+ tmp_variable
	40	topic=10001
	41	print ("Azim°: ","%.4f"%ang_azi,"Time:" ,"%.5f"%seconds)
	42	seconds_dec=(seconds-int(seconds))*1e6
	43	ang_azi_dec= (ang_azi-int(ang_azi))*1e3
	44	ang_elev_dec=(ang_elev-int(ang_elev))*1e3
	45	sleep(0.0088)
	46	socket.send_string("%d %d %d %d %d %d %d"%(topic,ang_elev,ang_elev_dec,ang_azi,ang_azi_dec,seconds,seconds_dec))
	47	t3 = time.time()
	48	print ("Total time for 1 vuelta in Seconds",t3-t2)

schainpy/scripts/test_sim00010.py created 644 +82 0

@@ -0,0 +1,82
	1	import os, sys
	2	import datetime
	3	import time
	4	from schainpy.controller import Project
	5
	6	desc = "USRP_test"
	7	filename = "USRP_processing.xml"
	8	controllerObj = Project()
	9	controllerObj.setup(id = '191', name='Test_USRP', description=desc)
	10
	11	############## USED TO PLOT IQ VOLTAGE, POWER AND SPECTRA #############
	12	######PATH DE LECTURA, ESCRITURA, GRAFICOS Y ENVIO WEB#################
	13	path = '/home/alex/Downloads/test_rawdata'
	14	figpath = '/home/alex/Downloads/hdf5_test'
	15	######################## UNIDAD DE LECTURA#############################
	16	'''
	17	readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
	18	path=path,
	19	startDate="2020/01/01", #"2020/01/01",#today,
	20	endDate= "2020/12/01", #"2020/12/30",#today,
	21	startTime='00:00:00',
	22	endTime='23:59:59',
	23	delay=0,
	24	#set=0,
	25	online=0,
	26	walk=1)
	27
	28	'''
	29	readUnitConfObj = controllerObj.addReadUnit(datatype='SimulatorReader',
	30	frequency=9.345e9,
	31	FixRCP_IPP= 60,
	32	Tau_0 = 30,
	33	AcqH0_0=0,
	34	samples=330,
	35	AcqDH_0=0.15,
	36	FixRCP_TXA=0.15,
	37	FixRCP_TXB=0.15,
	38	Fdoppler=600.0,
	39	Hdoppler=36,
	40	Adoppler=300,#300
	41	delay=0,
	42	online=0,
	43	walk=0,
	44	profilesPerBlock=625,
	45	dataBlocksPerFile=100)
	46	#nTotalReadFiles=2)
	47
	48
	49	#opObj11 = readUnitConfObj.addOperation(name='printInfo')
	50
	51	procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
	52
	53	procUnitConfObjB = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjA.getId())
	54	procUnitConfObjB.addParameter(name='nFFTPoints', value=625, format='int')
	55	procUnitConfObjB.addParameter(name='nProfiles', value=625, format='int')
	56
	57	opObj11 = procUnitConfObjB.addOperation(name='removeDC')
	58	opObj11.addParameter(name='mode', value=2)
	59	#opObj11 = procUnitConfObjB.addOperation(name='SpectraPlot')
	60	#opObj11 = procUnitConfObjB.addOperation(name='PowerProfilePlot')
	61
	62	procUnitConfObjC= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjB.getId())
	63	procUnitConfObjC.addOperation(name='SpectralMoments')
	64	#opObj11 = procUnitConfObjC.addOperation(name='PowerPlot')
	65
	66	'''
	67	opObj11 = procUnitConfObjC.addOperation(name='SpectralMomentsPlot')
	68	#opObj11.addParameter(name='xmin', value=14)
	69	#opObj11.addParameter(name='xmax', value=15)
	70	#opObj11.addParameter(name='save', value=figpath)
	71	opObj11.addParameter(name='showprofile', value=1)
	72	#opObj11.addParameter(name='save_period', value=10)
	73	'''
	74
	75	opObj10 = procUnitConfObjC.addOperation(name='ParameterWriter')
	76	opObj10.addParameter(name='path',value=figpath)
	77	#opObj10.addParameter(name='mode',value=0)
	78	opObj10.addParameter(name='blocksPerFile',value='100',format='int')
	79	opObj10.addParameter(name='metadataList',value='utctimeInit,timeInterval',format='list')
	80	opObj10.addParameter(name='dataList',value='data_POW,data_DOP,data_WIDTH,data_SNR')#,format='list'
	81
	82	controllerObj.start()

schainpy/scripts/wr_integrador.py created 644 +162 0

@@ -0,0 +1,162
	1	import os,numpy,h5py
	2	from shutil import copyfile
	3
	4	def isNumber(str):
	5	try:
	6	float(str)
	7	return True
	8	except:
	9	return False
	10
	11	def getfirstFilefromPath(path,meta,ext):
	12	validFilelist = []
	13	fileList = os.listdir(path)
	14	if len(fileList)<1:
	15	return None
	16	# meta 1234 567 8-18 BCDE
	17	# H,D,PE YYYY DDD EPOC .ext
	18
	19	for thisFile in fileList:
	20	if meta =="PE":
	21	try:
	22	number= int(thisFile[len(meta)+7:len(meta)+17])
	23	except:
	24	print("There is a file or folder with different format")
	25	if meta == "D":
	26	try:
	27	number= int(thisFile[8:11])
	28	except:
	29	print("There is a file or folder with different format")
	30
	31	if not isNumber(str=number):
	32	continue
	33	if (os.path.splitext(thisFile)[-1].lower() != ext.lower()):
	34	continue
	35	validFilelist.sort()
	36	validFilelist.append(thisFile)
	37	if len(validFilelist)>0:
	38	validFilelist = sorted(validFilelist,key=str.lower)
	39	return validFilelist
	40	return None
	41
	42	def gettimeutcfromDirFilename(path,file):
	43	dir_file= path+"/"+file
	44	fp = h5py.File(dir_file,'r')
	45	epoc = fp['Metadata'].get('utctimeInit')[()]
	46	fp.close()
	47	return epoc
	48
	49	def getDatavaluefromDirFilename(path,file,value):
	50	dir_file= path+"/"+file
	51	fp = h5py.File(dir_file,'r')
	52	array = fp['Data'].get(value)[()]
	53	fp.close()
	54	return array
	55
	56
	57	#·········· Velocidad de Pedestal·················
	58	w = input ("Ingresa velocidad de Pedestal: ")
	59	w = 4
	60	w = float(w)
	61	#·········· Resolucion minimo en grados···········
	62	alfa = input ("Ingresa resolucion minima en grados: ")
	63	alfa = 1
	64	alfa = float(alfa)
	65	#·········· IPP del Experimento ··················
	66	IPP = input ("Ingresa el IPP del experimento: ")
	67	IPP = 0.0004
	68	IPP = float(IPP)
	69	#·········· MODE ··················
	70	mode = input ("Ingresa el MODO del experimento T or F: ")
	71	mode = "T"
	72	mode = str(mode)
	73
	74	#·········· Tiempo en generar la resolucion min···
	75	#············ MCU ·· var_ang = w * (var_tiempo)···
	76	var_tiempo = alfa/w
	77	#·········· Tiempo Equivalente en perfiles········
	78	#·········· var_tiempo = IPP * ( num_perfiles )·
	79	num_perfiles = int(var_tiempo/IPP)
	80
	81	#··········DATA PEDESTAL··························
	82	dir_pedestal = "/home/alex/Downloads/pedestal"
	83	#·········· DATA ADQ······························
	84	if mode=="T":
	85	dir_adq = "/home/alex/Downloads/hdf5_testPP/d2020194" # Time domain
	86	else:
	87	dir_adq = "/home/alex/Downloads/hdf5_test/d2020194" # Frequency domain
	88
	89	print( "Velocidad angular :", w)
	90	print( "Resolucion minima en grados :", alfa)
	91	print( "Numero de perfiles equivalente:", num_perfiles)
	92	print( "Mode :", mode)
	93
	94	#············ First File·············
	95	list_pedestal = getfirstFilefromPath(path=dir_pedestal,meta="PE",ext=".hdf5")
	96	list_adq = getfirstFilefromPath(path=dir_adq ,meta="D",ext=".hdf5")
	97
	98	#············ utc time ··············
	99	utc_pedestal= gettimeutcfromDirFilename(path=dir_pedestal,file=list_pedestal[0])
	100	utc_adq = gettimeutcfromDirFilename(path=dir_adq ,file=list_adq[0])
	101
	102	print("utc_pedestal :",utc_pedestal)
	103	print("utc_adq :",utc_adq)
	104	#·············Relacion: utc_adq (+/-) var_tiempo*nro_file= utc_pedestal
	105	time_Interval_p = 0.01
	106	n_perfiles_p = 100
	107	if utc_adq>utc_pedestal:
	108	nro_file = int((int(utc_adq) - int(utc_pedestal))/(time_Interval_p*n_perfiles_p))
	109	ff_pedestal = list_pedestal[nro_file]
	110	utc_pedestal = gettimeutcfromDirFilename(path=dir_pedestal,file=ff_pedestal)
	111	nro_key_p = int((utc_adq-utc_pedestal)/time_Interval_p)
	112	if utc_adq >utc_pedestal:
	113	ff_pedestal = ff_pedestal
	114	else:
	115	nro_file = nro_file-1
	116	ff_pedestal = list_pedestal[nro_file]
	117	angulo = getDatavaluefromDirFilename(path=dir_pedestal,file=ff_pedestal,value="azimuth")
	118	nro_key_p = int((utc_adq-utc_pedestal)/time_Interval_p)
	119	print("nro_file :",nro_file)
	120	print("name_file :",ff_pedestal)
	121	print("utc_pedestal_file :",utc_pedestal)
	122	print("nro_key_p :",nro_key_p)
	123	print("utc_pedestal_init :",utc_pedestal+nro_key_p*time_Interval_p)
	124	print("angulo_array :",angulo[nro_key_p])
	125	#4+25+25+25+21
	126	#while True:
	127	list_pedestal = getfirstFilefromPath(path=dir_pedestal,meta="PE",ext=".hdf5")
	128	list_adq = getfirstFilefromPath(path=dir_adq ,meta="D",ext=".hdf5")
	129
	130	nro_file = nro_file #10
	131	nro_key_perfil = nro_key_p
	132	blocksPerFile = 100
	133	wr_path = "/home/alex/Downloads/hdf5_wr/"
	134	# Lectura de archivos de adquisicion para adicion de azimuth
	135	for thisFile in range(len(list_adq)):
	136	print("thisFileAdq",thisFile)
	137	angulo_adq = numpy.zeros(blocksPerFile)
	138	tmp = 0
	139	for j in range(blocksPerFile):
	140	iterador = nro_key_perfil + 25*(j-tmp)
	141	if iterador < n_perfiles_p:
	142	nro_file = nro_file
	143	else:
	144	nro_file = nro_file+1
	145	tmp = j
	146	iterador = nro_key_perfil
	147	ff_pedestal = list_pedestal[nro_file]
	148	angulo = getDatavaluefromDirFilename(path=dir_pedestal,file=ff_pedestal,value="azimuth")
	149	angulo_adq[j]= angulo[iterador]
	150	copyfile(dir_adq+"/"+list_adq[thisFile],wr_path+list_adq[thisFile])
	151	fp = h5py.File(wr_path+list_adq[thisFile],'a')
	152	grp = fp.create_group("Pedestal")
	153	dset = grp.create_dataset("azimuth" , data=angulo_adq)
	154	fp.close()
	155	print("Angulo",angulo_adq)
	156	print("Angulo",len(angulo_adq))
	157	nro_key_perfil=iterador + 25
	158	if nro_key_perfil< n_perfiles_p:
	159	nro_file = nro_file
	160	else:
	161	nro_file = nro_file+1
	162	nro_key_perfil= nro_key_p

schainpy/model/data/jrodata.py +12 -6

                 # data es un numpy array de 2 dmensiones (canales, alturas)
                 data = None
-                data_intensity = None
+                dataPP_POW   = None
-                data_velocity = None
+                dataPP_DOP   = None
-                data_specwidth = None
+                dataPP_WIDTH = None
+                dataPP_SNR   = None
                 def __init__(self):
                     '''
                     Constructor
                             self.flagDataAsBlock = dataOut.flagDataAsBlock
                             self.nProfiles = dataOut.nProfiles
                         if plot == 'pp_power':
-                            buffer = dataOut.data_intensity
+                            buffer = dataOut.dataPP_POWER
+                            self.flagDataAsBlock = dataOut.flagDataAsBlock
+                            self.nProfiles = dataOut.nProfiles
+                        if plot == 'pp_signal':
+                            buffer = dataOut.dataPP_POW
                             self.flagDataAsBlock = dataOut.flagDataAsBlock
                             self.nProfiles = dataOut.nProfiles
                         if plot == 'pp_velocity':
-                            buffer = dataOut.data_velocity
+                            buffer = dataOut.dataPP_DOP
                             self.flagDataAsBlock = dataOut.flagDataAsBlock
                             self.nProfiles = dataOut.nProfiles
                         if plot == 'pp_specwidth':
-                            buffer = dataOut.data_specwidth
+                            buffer = dataOut.dataPP_WIDTH
                             self.flagDataAsBlock = dataOut.flagDataAsBlock
                             self.nProfiles = dataOut.nProfiles

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 +14 -1

@@ -109,6 +109,20 class ParametersProc(ProcessingUnit):
109	self.dataOut.flagNoData = False	109	self.dataOut.flagNoData = False
110	self.dataOut.utctimeInit = self.dataIn.utctime	110	self.dataOut.utctimeInit = self.dataIn.utctime
111	self.dataOut.paramInterval = self.dataIn.nProfilesself.dataIn.nCohIntself.dataIn.ippSeconds	111	self.dataOut.paramInterval = self.dataIn.nProfilesself.dataIn.nCohIntself.dataIn.ippSeconds
		112	if hasattr(self.dataIn, 'dataPP_POW'):
		113	self.dataOut.dataPP_POW = self.dataIn.dataPP_POW
		114
		115	if hasattr(self.dataIn, 'dataPP_POWER'):
		116	self.dataOut.dataPP_POWER = self.dataIn.dataPP_POWER
		117
		118	if hasattr(self.dataIn, 'dataPP_DOP'):
		119	self.dataOut.dataPP_DOP = self.dataIn.dataPP_DOP
		120
		121	if hasattr(self.dataIn, 'dataPP_SNR'):
		122	self.dataOut.dataPP_SNR = self.dataIn.dataPP_SNR
		123
		124	if hasattr(self.dataIn, 'dataPP_WIDTH'):
		125	self.dataOut.dataPP_WIDTH = self.dataIn.dataPP_WIDTH
112	return	126	return
113		127
114	#---------------------- Spectra Data ---------------------------	128	#---------------------- Spectra Data ---------------------------
@@ -3984,4 +3998,3 class SMOperations():
3984	# error[indInvalid1] = 13	3998	# error[indInvalid1] = 13
3985	#	3999	#
3986	# return heights, error	4000	# return heights, error
3987	No newline at end of file

schainpy/model/proc/jroproc_voltage.py +59 -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.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
                     return dataOut
             # import collections
             # from scipy.stats import mode
             #

schainpy/scripts/test_sim0009.py +29 -5

@@ -25,8 +25,20 readUnitConfObj = controllerObj.addReadUnit(datatype='SimulatorReader',
25	delay=0,	25	delay=0,
26	online=0,	26	online=0,
27	walk=0,	27	walk=0,
28	~~nTotalReadFiles~~=3)	28	profilesPerBlock=625,
29		29	dataBlocksPerFile=100)#,#nTotalReadFiles=2)
		30	'''
		31	readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader',
		32	path=path,
		33	startDate="2020/01/01", #"2020/01/01",#today,
		34	endDate= "2020/12/01", #"2020/12/30",#today,
		35	startTime='00:00:00',
		36	endTime='23:59:59',
		37	delay=0,
		38	#set=0,
		39	online=0,
		40	walk=1)
		41	'''
30	opObj11 = readUnitConfObj.addOperation(name='printInfo')	42	opObj11 = readUnitConfObj.addOperation(name='printInfo')
31		43
32	procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())	44	procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId())
@@ -36,14 +48,26 procUnitConfObjA = controllerObj.addProcUnit(datatype='VoltageProc', inputId=rea
36	#opObj10 = procUnitConfObjA.addOperation(name='selectChannels')	48	#opObj10 = procUnitConfObjA.addOperation(name='selectChannels')
37	#opObj10.addParameter(name='channelList', value=[0])	49	#opObj10.addParameter(name='channelList', value=[0])
38	opObj11 = procUnitConfObjA.addOperation(name='PulsePairVoltage', optype='other')	50	opObj11 = procUnitConfObjA.addOperation(name='PulsePairVoltage', optype='other')
39	opObj11.addParameter(name='n', value='~~300~~', format='int')#10	51	opObj11.addParameter(name='n', value='625', format='int')#10
40	opObj11.addParameter(name='removeDC', value=1, format='int')	52	opObj11.addParameter(name='removeDC', value=1, format='int')
41		53
42	#opObj11 = procUnitConfObjA.addOperation(name='PulsepairPowerPlot', optype='other')	54	#opObj11 = procUnitConfObjA.addOperation(name='PulsepairPowerPlot', optype='other')
		55	#opObj11 = procUnitConfObjA.addOperation(name='PulsepairSignalPlot', optype='other')
		56
43		57
44	opObj11 = procUnitConfObjA.addOperation(name='PulsepairVelocityPlot', optype='other')	58	#opObj11 = procUnitConfObjA.addOperation(name='PulsepairVelocityPlot', optype='other')
45	#opObj11.addParameter(name='xmax', value=8)	59	#opObj11.addParameter(name='xmax', value=8)
46		60
47	opObj11 = procUnitConfObjA.addOperation(name='PulsepairSpecwidthPlot', optype='other')	61	#opObj11 = procUnitConfObjA.addOperation(name='PulsepairSpecwidthPlot', optype='other')
		62
		63	procUnitConfObjB= controllerObj.addProcUnit(datatype='ParametersProc',inputId=procUnitConfObjA.getId())
		64
		65
		66	opObj10 = procUnitConfObjB.addOperation(name='ParameterWriter')
		67	opObj10.addParameter(name='path',value=figpath)
		68	#opObj10.addParameter(name='mode',value=0)
		69	opObj10.addParameter(name='blocksPerFile',value='100',format='int')
		70	opObj10.addParameter(name='metadataList',value='utctimeInit,timeInterval',format='list')
		71	opObj10.addParameter(name='dataList',value='dataPP_POW,dataPP_DOP,dataPP_SNR,dataPP_WIDTH')#,format='list'
48		72
49	controllerObj.start()	73	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