schain Commit - r1296:162740cfc930 · Jicamarca Repository

update for weather radar options

avaldez -

r1296:162740cfc930

parent child

Context file:

r1296:162740cfc930

schainpy/model/io/jroIO_simulator.py created 644 +473 0

@@ -0,0 +1,473
	1	import numpy,math,random,time
	2	#---------------1 Heredamos JRODatareader
	3	from schainpy.model.io.jroIO_base import *
	4	#---------------2 Heredamos las propiedades de ProcessingUnit
	5	from schainpy.model.proc.jroproc_base import ProcessingUnit,Operation,MPDecorator
	6	#---------------3 Importaremos las clases BascicHeader, SystemHeader, RadarControlHeader, ProcessingHeader
	7	from schainpy.model.data.jroheaderIO import PROCFLAG, BasicHeader,SystemHeader,RadarControllerHeader, ProcessingHeader
	8	#---------------4 Importaremos el objeto Voltge
	9	from schainpy.model.data.jrodata import Voltage
	10
	11	class SimulatorReader(JRODataReader, ProcessingUnit):
	12	incIntFactor = 1
	13	nFFTPoints = 0
	14	FixPP_IncInt = 1
	15	FixRCP_IPP = 1000
	16	FixPP_CohInt = 1
	17	Tau_0 = 250
	18	AcqH0_0 = 70
	19	H0 = AcqH0_0
	20	AcqDH_0 = 1.25
	21	DH0 = AcqDH_0
	22	Bauds = 32
	23	BaudWidth = None
	24	FixRCP_TXA = 40
	25	FixRCP_TXB = 70
	26	fAngle = 2.0math.pi(1/16)
	27	DC_level = 500
	28	stdev = 8
	29	Num_Codes = 2
	30	#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])
	31	#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])
	32	#Dyn_snCode = numpy.array([Num_Codes,Bauds])
	33	Dyn_snCode = None
	34	Samples = 200
	35	channels = 5
	36	pulses = None
	37	Reference = None
	38	pulse_size = None
	39	prof_gen = None
	40	Fdoppler = 100
	41	Hdoppler = 36
	42	Adoppler = 300
	43	frequency = 9345
	44	nTotalReadFiles = 1000
	45
	46	def __init__(self):
	47	"""
	48	Inicializador de la clases SimulatorReader para
	49	generar datos de voltage simulados.
	50	Input:
	51	dataOut: Objeto de la clase Voltage.
	52	Este Objeto sera utilizado apra almacenar
	53	un perfil de datos cada vez qe se haga
	54	un requerimiento (getData)
	55	"""
	56	ProcessingUnit.__init__(self)
	57	print(" [ START ] init - Metodo Simulator Reader")
	58
	59	self.isConfig = False
	60	self.basicHeaderObj = BasicHeader(LOCALTIME)
	61	self.systemHeaderObj = SystemHeader()
	62	self.radarControllerHeaderObj = RadarControllerHeader()
	63	self.processingHeaderObj = ProcessingHeader()
	64	self.profileIndex = 2**32-1
	65	self.dataOut = Voltage()
	66	#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])
	67	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])
	68	#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])
	69	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])
	70	#self.Dyn_snCode = numpy.array([code0,code1])
	71	self.Dyn_snCode = None
	72
	73	def set_kwargs(self, **kwargs):
	74	for key, value in kwargs.items():
	75	setattr(self, key, value)
	76
	77	def __hasNotDataInBuffer(self):
	78
	79	if self.profileIndex >= self.processingHeaderObj.profilesPerBlock* self.nTxs:
	80	if self.nReadBlocks>0:
	81	tmp = self.dataOut.utctime
	82	tmp_utc = int(self.dataOut.utctime)
	83	tmp_milisecond = int((tmp-tmp_utc)*1000)
	84	self.basicHeaderObj.utc = tmp_utc
	85	self.basicHeaderObj.miliSecond= tmp_milisecond
	86	return 1
	87	return 0
	88
	89	def setNextFile(self):
	90	"""Set the next file to be readed open it and parse de file header"""
	91
	92	if (self.nReadBlocks >= self.processingHeaderObj.dataBlocksPerFile):
	93	self.nReadFiles=self.nReadFiles+1
	94	if self.nReadFiles ==self.nTotalReadFiles:
	95	self.flagNoMoreFiles=1
	96	raise schainpy.admin.SchainWarning('No more files to read')
	97
	98	print('------------------- [Opening file] ------------------------------',self.nReadFiles)
	99	self.nReadBlocks = 0
	100
	101	def __setNewBlock(self):
	102	self.setNextFile()
	103	if self.flagIsNewFile:
	104	return 1
	105
	106	def readNextBlock(self):
	107	while True:
	108	self.__setNewBlock()
	109	if not(self.readBlock()):
	110	return 0
	111	self.getBasicHeader()
	112	break
	113	if self.verbose:
	114	print("[Reading] Block No. %d/%d -> %s" %(self.nReadBlocks,
	115	self.processingHeaderObj.dataBlocksPerFile,
	116	self.dataOut.datatime.ctime()) )
	117	return 1
	118
	119	def getFirstHeader(self):
	120	self.getBasicHeader()
	121	self.dataOut.processingHeaderObj = self.processingHeaderObj.copy()
	122	self.dataOut.systemHeaderObj = self.systemHeaderObj.copy()
	123	self.dataOut.radarControllerHeaderObj = self.radarControllerHeaderObj.copy()
	124
	125	self.dataOut.nProfiles = self.processingHeaderObj.profilesPerBlock
	126	self.dataOut.heightList = numpy.arange(self.processingHeaderObj.nHeights) * self.processingHeaderObj.deltaHeight + self.processingHeaderObj.firstHeight
	127	self.dataOut.channelList = list(range(self.systemHeaderObj.nChannels))
	128	self.dataOut.nCohInt = self.processingHeaderObj.nCohInt
	129	# asumo q la data no esta decodificada
	130	self.dataOut.flagDecodeData = self.processingHeaderObj.flag_decode
	131	# asumo q la data no esta sin flip
	132	self.dataOut.flagDeflipData = self.processingHeaderObj.flag_deflip
	133	self.dataOut.flagShiftFFT = self.processingHeaderObj.shif_fft
	134	self.dataOut.frequency = self.frequency
	135
	136	def getBasicHeader(self):
	137	self.dataOut.utctime = self.basicHeaderObj.utc + self.basicHeaderObj.miliSecond / \
	138	1000. + self.profileIndex * self.radarControllerHeaderObj.ippSeconds
	139
	140	self.dataOut.flagDiscontinuousBlock = self.flagDiscontinuousBlock
	141	self.dataOut.timeZone = self.basicHeaderObj.timeZone
	142	self.dataOut.dstFlag = self.basicHeaderObj.dstFlag
	143	self.dataOut.errorCount = self.basicHeaderObj.errorCount
	144	self.dataOut.useLocalTime = self.basicHeaderObj.useLocalTime
	145	self.dataOut.ippSeconds = self.radarControllerHeaderObj.ippSeconds / self.nTxs
	146
	147	def set_RCH(self, expType=2, nTx=1,ipp=None, txA=0, txB=0,
	148	nWindows=None, nHeights=None, firstHeight=None, deltaHeight=None,
	149	numTaus=0, line6Function=0, line5Function=0, fClock=None,
	150	prePulseBefore=0, prePulseAfter=0,
	151	codeType=0, nCode=0, nBaud=0, code=None,
	152	flip1=0, flip2=0):
	153	self.radarControllerHeaderObj.expType = expType
	154	self.radarControllerHeaderObj.nTx = nTx
	155	self.radarControllerHeaderObj.ipp = float(ipp)
	156	self.radarControllerHeaderObj.txA = float(txA)
	157	self.radarControllerHeaderObj.txB = float(txB)
	158	self.radarControllerHeaderObj.rangeIPP = ipp
	159	self.radarControllerHeaderObj.rangeTxA = txA
	160	self.radarControllerHeaderObj.rangeTxB = txB
	161
	162	self.radarControllerHeaderObj.nHeights = int(nHeights)
	163	self.radarControllerHeaderObj.firstHeight = numpy.array([firstHeight])
	164	self.radarControllerHeaderObj.deltaHeight = numpy.array([deltaHeight])
	165	self.radarControllerHeaderObj.samplesWin = numpy.array([nHeights])
	166
	167
	168	self.radarControllerHeaderObj.nWindows = nWindows
	169	self.radarControllerHeaderObj.numTaus = numTaus
	170	self.radarControllerHeaderObj.codeType = codeType
	171	self.radarControllerHeaderObj.line6Function = line6Function
	172	self.radarControllerHeaderObj.line5Function = line5Function
	173	self.radarControllerHeaderObj.fclock = fClock
	174	self.radarControllerHeaderObj.prePulseBefore= prePulseBefore
	175	self.radarControllerHeaderObj.prePulseAfter = prePulseAfter
	176
	177	self.radarControllerHeaderObj.nCode = nCode
	178	self.radarControllerHeaderObj.nBaud = nBaud
	179	self.radarControllerHeaderObj.code = code
	180	self.radarControllerHeaderObj.flip1 = flip1
	181	self.radarControllerHeaderObj.flip2 = flip2
	182
	183	self.radarControllerHeaderObj.code_size = int(numpy.ceil(nBaud / 32.)) * nCode * 4
	184
	185	if fClock is None and deltaHeight is not None:
	186	self.fClock = 0.15 / (deltaHeight * 1e-6)
	187
	188	def set_PH(self, dtype=0, blockSize=0, profilesPerBlock=0,
	189	dataBlocksPerFile=0, nWindows=0, processFlags=0, nCohInt=0,
	190	nIncohInt=0, totalSpectra=0, nHeights=0, firstHeight=0,
	191	deltaHeight=0, samplesWin=0, spectraComb=0, nCode=0,
	192	code=0, nBaud=None, shif_fft=False, flag_dc=False,
	193	flag_cspc=False, flag_decode=False, flag_deflip=False):
	194
	195	self.processingHeaderObj.profilesPerBlock = profilesPerBlock
	196	self.processingHeaderObj.dataBlocksPerFile = dataBlocksPerFile
	197	self.processingHeaderObj.nWindows = nWindows
	198	self.processingHeaderObj.nCohInt = nCohInt
	199	self.processingHeaderObj.nIncohInt = nIncohInt
	200	self.processingHeaderObj.totalSpectra = totalSpectra
	201	self.processingHeaderObj.nHeights = int(nHeights)
	202	self.processingHeaderObj.firstHeight = firstHeight
	203	self.processingHeaderObj.deltaHeight = deltaHeight
	204	self.processingHeaderObj.samplesWin = nHeights
	205
	206	def set_BH(self, utc = 0, miliSecond = 0, timeZone = 0):
	207	self.basicHeaderObj.utc = utc
	208	self.basicHeaderObj.miliSecond = miliSecond
	209	self.basicHeaderObj.timeZone = timeZone
	210
	211	def set_SH(self, nSamples=0, nProfiles=0, nChannels=0, adcResolution=14, pciDioBusWidth=0):
	212	self.systemHeaderObj.nSamples = nSamples
	213	self.systemHeaderObj.nProfiles = nProfiles
	214	self.systemHeaderObj.nChannels = nChannels
	215	self.systemHeaderObj.adcResolution = adcResolution
	216	self.systemHeaderObj.pciDioBusWidth = pciDioBusWidth
	217
	218	def init_acquisition(self):
	219
	220	if self.nFFTPoints != 0:
	221	self.incIntFactor = m_nProfilesperBlock/self.nFFTPoints
	222	if (self.FixPP_IncInt > self.incIntFactor):
	223	self.incIntFactor = self.FixPP_IncInt/ self.incIntFactor
	224	elif(self.FixPP_IncInt< self.incIntFactor):
	225	print("False alert...")
	226
	227	ProfilesperBlock = self.processingHeaderObj.profilesPerBlock
	228
	229	self.timeperblock =int(((self.FixRCP_IPP
	230	*ProfilesperBlock
	231	*self.FixPP_CohInt
	232	*self.incIntFactor)
	233	/150.0)
	234	*0.9
	235	+0.5)
	236	# para cada canal
	237	self.profiles = ProfilesperBlock*self.FixPP_CohInt
	238	self.profiles = ProfilesperBlock
	239	self.Reference = int((self.Tau_0-self.AcqH0_0)/(self.AcqDH_0)+0.5)
	240	self.BaudWidth = int((self.FixRCP_TXA/self.AcqDH_0)/self.Bauds + 0.5 )
	241
	242	if (self.BaudWidth==0):
	243	self.BaudWidth=1
	244
	245	def init_pulse(self,Num_Codes=Num_Codes,Bauds=Bauds,BaudWidth=BaudWidth,Dyn_snCode=Dyn_snCode):
	246
	247	Num_Codes = Num_Codes
	248	Bauds = Bauds
	249	BaudWidth = BaudWidth
	250	Dyn_snCode = Dyn_snCode
	251
	252	if Dyn_snCode:
	253	print("EXISTE")
	254	else:
	255	print("No existe")
	256
	257	if Dyn_snCode: # if Bauds:
	258	pulses = list(range(0,Num_Codes))
	259	num_codes = Num_Codes
	260	for i in range(num_codes):
	261	pulse_size = Bauds*BaudWidth
	262	pulses[i] = numpy.zeros(pulse_size)
	263	for j in range(Bauds):
	264	for k in range(BaudWidth):
	265	pulses[i][jBaudWidth+k] = int(Dyn_snCode[i][j]600)
	266	else:
	267	print("sin code")
	268	pulses = list(range(1))
	269	if self.AcqDH_0>0.149:
	270	pulse_size = int(self.FixRCP_TXB/0.15+0.5)
	271	else:
	272	pulse_size = int((self.FixRCP_TXB/self.AcqDH_0)+0.5) #0.0375
	273	pulses[0] = numpy.ones(pulse_size)
	274	pulses = 600*pulses[0]
	275
	276	return pulses,pulse_size
	277
	278	def jro_GenerateBlockOfData(self,Samples=Samples,DC_level= DC_level,stdev=stdev,
	279	Reference= Reference,pulses= pulses,
	280	Num_Codes= Num_Codes,pulse_size=pulse_size,
	281	prof_gen= prof_gen,H0 = H0,DH0=DH0,
	282	Adoppler=Adoppler,Fdoppler= Fdoppler,Hdoppler=Hdoppler):
	283	Samples = Samples
	284	DC_level = DC_level
	285	stdev = stdev
	286	m_nR = Reference
	287	pulses = pulses
	288	num_codes = Num_Codes
	289	ps = pulse_size
	290	prof_gen = prof_gen
	291	channels = self.channels
	292	H0 = H0
	293	DH0 = DH0
	294	ippSec = self.radarControllerHeaderObj.ippSeconds
	295	Fdoppler = self.Fdoppler
	296	Hdoppler = self.Hdoppler
	297	Adoppler = self.Adoppler
	298
	299	self.datablock = numpy.zeros([channels,prof_gen,Samples],dtype= numpy.complex64)
	300	for i in range(channels):
	301	for k in range(prof_gen):
	302	#·······················NOISE···············
	303	Noise_r = numpy.random.normal(DC_level,stdev,Samples)
	304	Noise_i = numpy.random.normal(DC_level,stdev,Samples)
	305	Noise = numpy.zeros(Samples,dtype=complex)
	306	Noise.real = Noise_r
	307	Noise.imag = Noise_i
	308	#·······················PULSOS··············
	309	Pulso = numpy.zeros(pulse_size,dtype=complex)
	310	Pulso.real = pulses[k%num_codes]
	311	Pulso.imag = pulses[k%num_codes]
	312	#····················· PULSES+NOISE··········
	313	InBuffer = numpy.zeros(Samples,dtype=complex)
	314	InBuffer[m_nR:m_nR+ps] = Pulso
	315	InBuffer = InBuffer+Noise
	316	#····················· ANGLE ·······························
	317	InBuffer.real[m_nR:m_nR+ps] = InBuffer.real[m_nR:m_nR+ps](math.cos( self.fAngle)5)
	318	InBuffer.imag[m_nR:m_nR+ps] = InBuffer.imag[m_nR:m_nR+ps](math.sin( self.fAngle)5)
	319	InBuffer=InBuffer
	320	self.datablock[i][k]= InBuffer
	321	#plot_cts(InBuffer,H0=H0,DH0=DH0)
	322	#wave_fft(x=InBuffer,plot_show=True)
	323	#time.sleep(1)
	324	#················DOPPLER SIGNAL...............................................
	325	time_vec = numpy.linspace(0,(prof_gen-1)ippSec,int(prof_gen))+self.nReadBlocksippSecprof_gen+(self.nReadFiles-1)ippSec*prof_gen
	326	fd = Fdoppler #+(600.0/120)*self.nReadBlocks
	327	d_signal = Adopplernumpy.array(numpy.exp(1.0j2.0math.pifd*time_vec),dtype=numpy.complex64)
	328	#·················· DATABLOCK + DOPPLER············...........................
	329	HD=int(Hdoppler/self.AcqDH_0)
	330	for i in range(12):
	331	self.datablock[:,:,HD+i]=self.datablock[:,:,HD+i]+ d_signal # RESULT
	332	'''
	333	a= numpy.zeros(10)
	334	for i in range(10):
	335	a[i]=i+self.nReadBlocks+20
	336	for i in a:
	337	self.datablock[0,:,int(i)]=self.datablock[0,:,int(i)]+ d_signal # RESULT
	338	'''
	339
	340	def readBlock(self):
	341
	342	self.jro_GenerateBlockOfData(Samples= self.samples,DC_level=self.DC_level,
	343	stdev=self.stdev,Reference= self.Reference,
	344	pulses = self.pulses,Num_Codes=self.Num_Codes,
	345	pulse_size=self.pulse_size,prof_gen=self.profiles,
	346	H0=self.H0,DH0=self.DH0)
	347
	348	self.profileIndex = 0
	349	self.flagIsNewFile = 0
	350	self.flagIsNewBlock = 1
	351	self.nTotalBlocks += 1
	352	self.nReadBlocks += 1
	353
	354	return 1
	355
	356
	357	def getData(self):
	358	if self.flagNoMoreFiles:
	359	self.dataOut.flagNodata = True
	360	return 0
	361	self.flagDiscontinuousBlock = 0
	362	self.flagIsNewBlock = 0
	363	if self.__hasNotDataInBuffer(): # aqui es verdad
	364	if not(self.readNextBlock()): # return 1 y por eso el if not salta a getBasic Header
	365	return 0
	366	self.getFirstHeader() # atributo
	367
	368	if not self.getByBlock:
	369	self.dataOut.flagDataAsBlock = False
	370	self.dataOut.data = self.datablock[:, self.profileIndex, :]
	371	self.dataOut.profileIndex = self.profileIndex
	372	self.profileIndex += 1
	373	else:
	374	pass
	375	self.dataOut.flagNoData = False
	376	self.getBasicHeader()
	377	self.dataOut.realtime = self.online
	378	return self.dataOut.data
	379
	380
	381	def setup(self,frequency=49.92e6,incIntFactor= 1, nFFTPoints = 0, FixPP_IncInt=1,FixRCP_IPP=1000,
	382	FixPP_CohInt= 1,Tau_0= 250,AcqH0_0 = 70 ,AcqDH_0=1.25, Bauds= 32,
	383	FixRCP_TXA = 40, FixRCP_TXB = 50, fAngle = 2.0math.pi(1/16),DC_level= 50,
	384	stdev= 8,Num_Codes = 1 , Dyn_snCode = None, samples=200,
	385	channels=2,Fdoppler=20,Hdoppler=36,Adoppler=500,nTotalReadFiles=10000,
	386	**kwargs):
	387
	388	self.set_kwargs(**kwargs)
	389	self.nReadBlocks = 0
	390	self.nReadFiles = 1
	391	print('------------------- [Opening file: ] ------------------------------',self.nReadFiles)
	392
	393	tmp = time.time()
	394	tmp_utc = int(tmp)
	395	tmp_milisecond = int((tmp-tmp_utc)*1000)
	396	print(" SETUP -basicHeaderObj.utc",datetime.datetime.utcfromtimestamp(tmp))
	397	if Dyn_snCode is None:
	398	Num_Codes=1
	399	Bauds =1
	400
	401
	402
	403	self.set_BH(utc= tmp_utc,miliSecond= tmp_milisecond,timeZone=300 )
	404
	405	self.set_RCH( expType=0, nTx=150,ipp=FixRCP_IPP, txA=FixRCP_TXA, txB= FixRCP_TXB,
	406	nWindows=1 , nHeights=samples, firstHeight=AcqH0_0, deltaHeight=AcqDH_0,
	407	numTaus=1, line6Function=0, line5Function=0, fClock=None,
	408	prePulseBefore=0, prePulseAfter=0,
	409	codeType=14, nCode=Num_Codes, nBaud=32, code=Dyn_snCode,
	410	flip1=0, flip2=0)
	411
	412	self.set_PH(dtype=0, blockSize=0, profilesPerBlock=300,
	413	dataBlocksPerFile=120, nWindows=1, processFlags=0, nCohInt=1,
	414	nIncohInt=1, totalSpectra=0, nHeights=samples, firstHeight=AcqH0_0,
	415	deltaHeight=AcqDH_0, samplesWin=samples, spectraComb=0, nCode=0,
	416	code=0, nBaud=None, shif_fft=False, flag_dc=False,
	417	flag_cspc=False, flag_decode=False, flag_deflip=False)
	418
	419	self.set_SH(nSamples=samples, nProfiles=300, nChannels=channels)
	420
	421
	422	self.frequency = frequency
	423	self.incIntFactor = incIntFactor
	424	self.nFFTPoints = nFFTPoints
	425	self.FixPP_IncInt = FixPP_IncInt
	426	self.FixRCP_IPP = FixRCP_IPP
	427	self.FixPP_CohInt = FixPP_CohInt
	428	self.Tau_0 = Tau_0
	429	self.AcqH0_0 = AcqH0_0
	430	self.H0 = AcqH0_0
	431	self.AcqDH_0 = AcqDH_0
	432	self.DH0 = AcqDH_0
	433	self.Bauds = Bauds
	434	self.FixRCP_TXA = FixRCP_TXA
	435	self.FixRCP_TXB = FixRCP_TXB
	436	self.fAngle = fAngle
	437	self.DC_level = DC_level
	438	self.stdev = stdev
	439	self.Num_Codes = Num_Codes
	440	self.Dyn_snCode = Dyn_snCode
	441	self.samples = samples
	442	self.channels = channels
	443	self.profiles = None
	444	self.m_nReference = None
	445	self.Baudwidth = None
	446	self.Fdoppler = Fdoppler
	447	self.Hdoppler = Hdoppler
	448	self.Adoppler = Adoppler
	449	self.nTotalReadFiles = int(nTotalReadFiles)
	450
	451	print("IPP ", self.FixRCP_IPP)
	452	print("Tau_0 ",self.Tau_0)
	453	print("AcqH0_0",self.AcqH0_0)
	454	print("samples,window ",self.samples)
	455	print("AcqDH_0",AcqDH_0)
	456	print("FixRCP_TXA",self.FixRCP_TXA)
	457	print("FixRCP_TXB",self.FixRCP_TXB)
	458	print("Dyn_snCode",Dyn_snCode)
	459	print("Fdoppler", Fdoppler)
	460	print("Hdoppler",Hdoppler)
	461	print("Vdopplermax",Fdoppler*(3.0e8/self.frequency)/2.0)
	462	print("nTotalReadFiles", nTotalReadFiles)
	463
	464	self.init_acquisition()
	465	self.pulses,self.pulse_size=self.init_pulse(Num_Codes=self.Num_Codes,Bauds=self.Bauds,BaudWidth=self.BaudWidth,Dyn_snCode=Dyn_snCode)
	466	print(" [ END ] - SETUP metodo")
	467	return
	468
	469	def run(self,**kwargs): # metodo propio
	470	if not(self.isConfig):
	471	self.setup(**kwargs)
	472	self.isConfig = True
	473	self.getData()

schainpy/model/data/jrodata.py +23 -16

             '''
             $Author: murco $
             $Id: JROData.py 173 2012-11-20 15:06:21Z murco $
             '''
             import copy
             import numpy
             import datetime
             import json
             import schainpy.admin
             from schainpy.utils import log
             from .jroheaderIO import SystemHeader, RadarControllerHeader
             from schainpy.model.data import _noise
             def getNumpyDtype(dataTypeCode):
                 if dataTypeCode == 0:
                     numpyDtype = numpy.dtype([('real', '<i1'), ('imag', '<i1')])
                 elif dataTypeCode == 1:
                     numpyDtype = numpy.dtype([('real', '<i2'), ('imag', '<i2')])
                 elif dataTypeCode == 2:
                     numpyDtype = numpy.dtype([('real', '<i4'), ('imag', '<i4')])
                 elif dataTypeCode == 3:
                     numpyDtype = numpy.dtype([('real', '<i8'), ('imag', '<i8')])
                 elif dataTypeCode == 4:
                     numpyDtype = numpy.dtype([('real', '<f4'), ('imag', '<f4')])
                 elif dataTypeCode == 5:
                     numpyDtype = numpy.dtype([('real', '<f8'), ('imag', '<f8')])
                 else:
                     raise ValueError('dataTypeCode was not defined')
                 return numpyDtype
             def getDataTypeCode(numpyDtype):
                 if numpyDtype == numpy.dtype([('real', '<i1'), ('imag', '<i1')]):
                     datatype = 0
                 elif numpyDtype == numpy.dtype([('real', '<i2'), ('imag', '<i2')]):
                     datatype = 1
                 elif numpyDtype == numpy.dtype([('real', '<i4'), ('imag', '<i4')]):
                     datatype = 2
                 elif numpyDtype == numpy.dtype([('real', '<i8'), ('imag', '<i8')]):
                     datatype = 3
                 elif numpyDtype == numpy.dtype([('real', '<f4'), ('imag', '<f4')]):
                     datatype = 4
                 elif numpyDtype == numpy.dtype([('real', '<f8'), ('imag', '<f8')]):
                     datatype = 5
                 else:
                     datatype = None
                 return datatype
             def hildebrand_sekhon(data, navg):
                 """
                 This method is for the objective determination of the noise level in Doppler spectra. This
                 implementation technique is based on the fact that the standard deviation of the spectral
                 densities is equal to the mean spectral density for white Gaussian noise
                 Inputs:
                     Data    :    heights
                     navg    :    numbers of averages
                 Return:
                     mean    :    noise's level
                 """
                 sortdata = numpy.sort(data, axis=None)
                 '''
                 lenOfData = len(sortdata)
                 nums_min = lenOfData*0.2
                 if nums_min <= 5:
                     nums_min = 5
                 sump = 0.
                 sumq = 0.
                 j = 0
                 cont = 1
                 while((cont == 1)and(j < lenOfData)):
                     sump += sortdata[j]
                     sumq += sortdata[j]**2
                     if j > nums_min:
                         rtest = float(j)/(j-1) + 1.0/navg
                         if ((sumq*j) > (rtest*sump**2)):
                             j = j - 1
                             sump = sump - sortdata[j]
                             sumq = sumq - sortdata[j]**2
                             cont = 0
                     j += 1
                 lnoise = sump / j
                 '''
                 return _noise.hildebrand_sekhon(sortdata, navg)
             class Beam:
                 def __init__(self):
                     self.codeList = []
                     self.azimuthList = []
                     self.zenithList = []
             class GenericData(object):
                 flagNoData = True
                 def copy(self, inputObj=None):
                     if inputObj == None:
                         return copy.deepcopy(self)
                     for key in list(inputObj.__dict__.keys()):
                         attribute = inputObj.__dict__[key]
                         # If this attribute is a tuple or list
                         if type(inputObj.__dict__[key]) in (tuple, list):
                             self.__dict__[key] = attribute[:]
                             continue
                         # If this attribute is another object or instance
                         if hasattr(attribute, '__dict__'):
                             self.__dict__[key] = attribute.copy()
                             continue
                         self.__dict__[key] = inputObj.__dict__[key]
                 def deepcopy(self):
                     return copy.deepcopy(self)
                 def isEmpty(self):
                     return self.flagNoData
                 def isReady(self):
                     return not self.flagNoData
             class JROData(GenericData):
                 #    m_BasicHeader = BasicHeader()
                 #    m_ProcessingHeader = ProcessingHeader()
                 systemHeaderObj = SystemHeader()
                 radarControllerHeaderObj = RadarControllerHeader()
             #    data = None
                 type = None
                 datatype = None  # dtype but in string
             #     dtype = None
             #    nChannels = None
             #    nHeights = None
                 nProfiles = None
                 heightList = None
                 channelList = None
                 flagDiscontinuousBlock = False
                 useLocalTime = False
                 utctime = None
                 timeZone = None
                 dstFlag = None
                 errorCount = None
                 blocksize = None
             #     nCode = None
             #     nBaud = None
             #     code = None
                 flagDecodeData = False  # asumo q la data no esta decodificada
                 flagDeflipData = False  # asumo q la data no esta sin flip
                 flagShiftFFT = False
             #     ippSeconds = None
             #     timeInterval = None
                 nCohInt = None
             #     noise = None
                 windowOfFilter = 1
                 # Speed of ligth
                 C = 3e8
                 frequency = 49.92e6
                 realtime = False
                 beacon_heiIndexList = None
                 last_block = None
                 blocknow = None
                 azimuth = None
                 zenith = None
                 beam = Beam()
                 profileIndex = None
                 error = None
                 data = None
                 nmodes = None
                 def __str__(self):
                     return '{} - {}'.format(self.type, self.getDatatime())
                 def getNoise(self):
                     raise NotImplementedError
                 def getNChannels(self):
                     return len(self.channelList)
                 def getChannelIndexList(self):
                     return list(range(self.nChannels))
                 def getNHeights(self):
                     return len(self.heightList)
                 def getHeiRange(self, extrapoints=0):
                     heis = self.heightList
             #        deltah = self.heightList[1] - self.heightList[0]
             #
             #        heis.append(self.heightList[-1])
                     return heis
                 def getDeltaH(self):
                     delta = self.heightList[1] - self.heightList[0]
                     return delta
                 def getltctime(self):
                     if self.useLocalTime:
                         return self.utctime - self.timeZone * 60
                     return self.utctime
                 def getDatatime(self):
                     datatimeValue = datetime.datetime.utcfromtimestamp(self.ltctime)
                     return datatimeValue
                 def getTimeRange(self):
                     datatime = []
                     datatime.append(self.ltctime)
                     datatime.append(self.ltctime + self.timeInterval + 1)
                     datatime = numpy.array(datatime)
                     return datatime
                 def getFmaxTimeResponse(self):
                     period = (10**-6) * self.getDeltaH() / (0.15)
                     PRF = 1. / (period * self.nCohInt)
                     fmax = PRF
                     return fmax
                 def getFmax(self):
                     PRF = 1. / (self.ippSeconds * self.nCohInt)
                     fmax = PRF
                     return fmax
                 def getVmax(self):
                     _lambda = self.C / self.frequency
                     vmax = self.getFmax() * _lambda / 2
                     return vmax
                 def get_ippSeconds(self):
                     '''
                     '''
                     return self.radarControllerHeaderObj.ippSeconds
                 def set_ippSeconds(self, ippSeconds):
                     '''
                     '''
                     self.radarControllerHeaderObj.ippSeconds = ippSeconds
                     return
                 def get_dtype(self):
                     '''
                     '''
                     return getNumpyDtype(self.datatype)
                 def set_dtype(self, numpyDtype):
                     '''
                     '''
                     self.datatype = getDataTypeCode(numpyDtype)
                 def get_code(self):
                     '''
                     '''
                     return self.radarControllerHeaderObj.code
                 def set_code(self, code):
                     '''
                     '''
                     self.radarControllerHeaderObj.code = code
                     return
                 def get_ncode(self):
                     '''
                     '''
                     return self.radarControllerHeaderObj.nCode
                 def set_ncode(self, nCode):
                     '''
                     '''
                     self.radarControllerHeaderObj.nCode = nCode
                     return
                 def get_nbaud(self):
                     '''
                     '''
                     return self.radarControllerHeaderObj.nBaud
                 def set_nbaud(self, nBaud):
                     '''
                     '''
                     self.radarControllerHeaderObj.nBaud = nBaud
                     return
                 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
                 channelIndexList = property(
                     getChannelIndexList, "I'm the 'channelIndexList' property.")
                 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
                 #noise = property(getNoise, "I'm the 'nHeights' property.")
                 datatime = property(getDatatime, "I'm the 'datatime' property")
                 ltctime = property(getltctime, "I'm the 'ltctime' property")
                 ippSeconds = property(get_ippSeconds, set_ippSeconds)
                 dtype = property(get_dtype, set_dtype)
             #     timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
                 code = property(get_code, set_code)
                 nCode = property(get_ncode, set_ncode)
                 nBaud = property(get_nbaud, set_nbaud)
             class Voltage(JROData):
                 # data es un numpy array de 2 dmensiones (canales, alturas)
                 data = None
+                data_intensity = None
+                data_velocity = None
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.useLocalTime = True
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Voltage"
                     self.data = None
             #         self.dtype = None
             #        self.nChannels = 0
             #        self.nHeights = 0
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #        self.channelIndexList = None
                     self.flagNoData = True
                     self.flagDiscontinuousBlock = False
                     self.utctime = None
                     self.timeZone = None
                     self.dstFlag = None
                     self.errorCount = None
                     self.nCohInt = None
                     self.blocksize = None
                     self.flagDecodeData = False  # asumo q la data no esta decodificada
                     self.flagDeflipData = False  # asumo q la data no esta sin flip
                     self.flagShiftFFT = False
                     self.flagDataAsBlock = False  # Asumo que la data es leida perfil a perfil
                     self.profileIndex = 0
                 def getNoisebyHildebrand(self, channel=None):
                     """
                     Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
                     Return:
                         noiselevel
                     """
                     if channel != None:
                         data = self.data[channel]
                         nChannels = 1
                     else:
                         data = self.data
                         nChannels = self.nChannels
                     noise = numpy.zeros(nChannels)
                     power = data * numpy.conjugate(data)
                     for thisChannel in range(nChannels):
                         if nChannels == 1:
                             daux = power[:].real
                         else:
                             daux = power[thisChannel, :].real
                         noise[thisChannel] = hildebrand_sekhon(daux, self.nCohInt)
                     return noise
                 def getNoise(self, type=1, channel=None):
                     if type == 1:
                         noise = self.getNoisebyHildebrand(channel)
                     return noise
                 def getPower(self, channel=None):
                     if channel != None:
                         data = self.data[channel]
                     else:
                         data = self.data
                     power = data * numpy.conjugate(data)
                     powerdB = 10 * numpy.log10(power.real)
                     powerdB = numpy.squeeze(powerdB)
                     return powerdB
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt
                     return timeInterval
                 noise = property(getNoise, "I'm the 'nHeights' property.")
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
             class Spectra(JROData):
                 # data spc es un numpy array de 2 dmensiones (canales, perfiles, alturas)
                 data_spc = None
                 # data cspc es un numpy array de 2 dmensiones (canales, pares, alturas)
                 data_cspc = None
                 # data dc es un numpy array de 2 dmensiones (canales, alturas)
                 data_dc = None
                 # data power
                 data_pwr = None
                 nFFTPoints = None
             #     nPairs = None
                 pairsList = None
                 nIncohInt = None
                 wavelength = None  # Necesario para cacular el rango de velocidad desde la frecuencia
                 nCohInt = None  # se requiere para determinar el valor de timeInterval
                 ippFactor = None
                 profileIndex = 0
                 plotting = "spectra"
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.useLocalTime = True
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Spectra"
             #        self.data = None
             #         self.dtype = None
             #        self.nChannels = 0
             #        self.nHeights = 0
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #        self.channelIndexList = None
                     self.pairsList = None
                     self.flagNoData = True
                     self.flagDiscontinuousBlock = False
                     self.utctime = None
                     self.nCohInt = None
                     self.nIncohInt = None
                     self.blocksize = None
                     self.nFFTPoints = None
                     self.wavelength = None
                     self.flagDecodeData = False  # asumo q la data no esta decodificada
                     self.flagDeflipData = False  # asumo q la data no esta sin flip
                     self.flagShiftFFT = False
                     self.ippFactor = 1
                     #self.noise = None
                     self.beacon_heiIndexList = []
                     self.noise_estimation = None
                 def getNoisebyHildebrand(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
                     """
                     Determino el nivel de ruido usando el metodo Hildebrand-Sekhon
                     Return:
                         noiselevel
                     """
                     noise = numpy.zeros(self.nChannels)
                     for channel in range(self.nChannels):
                         daux = self.data_spc[channel,
                                              xmin_index:xmax_index, ymin_index:ymax_index]
                         noise[channel] = hildebrand_sekhon(daux, self.nIncohInt)
                     return noise
                 def getNoise(self, xmin_index=None, xmax_index=None, ymin_index=None, ymax_index=None):
                     if self.noise_estimation is not None:
                         # this was estimated by getNoise Operation defined in jroproc_spectra.py
                         return self.noise_estimation
                     else:
                         noise = self.getNoisebyHildebrand(
                             xmin_index, xmax_index, ymin_index, ymax_index)
                         return noise
                 def getFreqRangeTimeResponse(self, extrapoints=0):
                     deltafreq = self.getFmaxTimeResponse() / (self.nFFTPoints * self.ippFactor)
                     freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2
                     return freqrange
                 def getAcfRange(self, extrapoints=0):
                     deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor))
                     freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
                     return freqrange
                 def getFreqRange(self, extrapoints=0):
                     deltafreq = self.getFmax() / (self.nFFTPoints * self.ippFactor)
                     freqrange = deltafreq * (numpy.arange(self.nFFTPoints + extrapoints) -self.nFFTPoints / 2.) - deltafreq / 2
                     return freqrange
                 def getVelRange(self, extrapoints=0):
                     deltav = self.getVmax() / (self.nFFTPoints * self.ippFactor)
                     velrange = deltav * (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.)
                     if self.nmodes:
                         return velrange/self.nmodes
                     else:
                         return velrange
                 def getNPairs(self):
                     return len(self.pairsList)
                 def getPairsIndexList(self):
                     return list(range(self.nPairs))
                 def getNormFactor(self):
                     pwcode = 1
                     if self.flagDecodeData:
                         pwcode = numpy.sum(self.code[0]**2)
                     #normFactor = min(self.nFFTPoints,self.nProfiles)*self.nIncohInt*self.nCohInt*pwcode*self.windowOfFilter
                     normFactor = self.nProfiles * self.nIncohInt * self.nCohInt * pwcode * self.windowOfFilter
                     return normFactor
                 def getFlagCspc(self):
                     if self.data_cspc is None:
                         return True
                     return False
                 def getFlagDc(self):
                     if self.data_dc is None:
                         return True
                     return False
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt * self.nProfiles * self.ippFactor
                     if self.nmodes:
                         return self.nmodes*timeInterval
                     else:
                         return timeInterval
                 def getPower(self):
                     factor = self.normFactor
                     z = self.data_spc / factor
                     z = numpy.where(numpy.isfinite(z), z, numpy.NAN)
                     avg = numpy.average(z, axis=1)
                     return 10 * numpy.log10(avg)
                 def getCoherence(self, pairsList=None, phase=False):
                     z = []
                     if pairsList is None:
                         pairsIndexList = self.pairsIndexList
                     else:
                         pairsIndexList = []
                         for pair in pairsList:
                             if pair not in self.pairsList:
                                 raise ValueError("Pair %s is not in dataOut.pairsList" % (
                                     pair))
                             pairsIndexList.append(self.pairsList.index(pair))
                     for i in range(len(pairsIndexList)):
                         pair = self.pairsList[pairsIndexList[i]]
                         ccf = numpy.average(self.data_cspc[pairsIndexList[i], :, :], axis=0)
                         powa = numpy.average(self.data_spc[pair[0], :, :], axis=0)
                         powb = numpy.average(self.data_spc[pair[1], :, :], axis=0)
                         avgcoherenceComplex = ccf / numpy.sqrt(powa * powb)
                         if phase:
                             data = numpy.arctan2(avgcoherenceComplex.imag,
                                                  avgcoherenceComplex.real) * 180 / numpy.pi
                         else:
                             data = numpy.abs(avgcoherenceComplex)
                         z.append(data)
                     return numpy.array(z)
                 def setValue(self, value):
                     print("This property should not be initialized")
                     return
                 nPairs = property(getNPairs, setValue, "I'm the 'nPairs' property.")
                 pairsIndexList = property(
                     getPairsIndexList, setValue, "I'm the 'pairsIndexList' property.")
                 normFactor = property(getNormFactor, setValue,
                                       "I'm the 'getNormFactor' property.")
                 flag_cspc = property(getFlagCspc, setValue)
                 flag_dc = property(getFlagDc, setValue)
                 noise = property(getNoise, setValue, "I'm the 'nHeights' property.")
                 timeInterval = property(getTimeInterval, setValue,
                                         "I'm the 'timeInterval' property")
             class SpectraHeis(Spectra):
                 data_spc = None
                 data_cspc = None
                 data_dc = None
                 nFFTPoints = None
             #     nPairs = None
                 pairsList = None
                 nCohInt = None
                 nIncohInt = None
                 def __init__(self):
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "SpectraHeis"
             #         self.dtype = None
             #        self.nChannels = 0
             #        self.nHeights = 0
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #        self.channelIndexList = None
                     self.flagNoData = True
                     self.flagDiscontinuousBlock = False
             #         self.nPairs = 0
                     self.utctime = None
                     self.blocksize = None
                     self.profileIndex = 0
                     self.nCohInt = 1
                     self.nIncohInt = 1
                 def getNormFactor(self):
                     pwcode = 1
                     if self.flagDecodeData:
                         pwcode = numpy.sum(self.code[0]**2)
                     normFactor = self.nIncohInt * self.nCohInt * pwcode
                     return normFactor
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
                     return timeInterval
                 normFactor = property(getNormFactor, "I'm the 'getNormFactor' property.")
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
             class Fits(JROData):
                 heightList = None
                 channelList = None
                 flagNoData = True
                 flagDiscontinuousBlock = False
                 useLocalTime = False
                 utctime = None
                 timeZone = None
             #     ippSeconds = None
             #     timeInterval = None
                 nCohInt = None
                 nIncohInt = None
                 noise = None
                 windowOfFilter = 1
                 # Speed of ligth
                 C = 3e8
                 frequency = 49.92e6
                 realtime = False
                 def __init__(self):
                     self.type = "Fits"
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
             #         self.channelIndexList = None
                     self.flagNoData = True
                     self.utctime = None
                     self.nCohInt = 1
                     self.nIncohInt = 1
                     self.useLocalTime = True
                     self.profileIndex = 0
             #         self.utctime = None
             #         self.timeZone = None
             #         self.ltctime = None
             #         self.timeInterval = None
             #         self.header = None
             #         self.data_header = None
             #         self.data = None
             #         self.datatime = None
             #         self.flagNoData = False
             #         self.expName = ''
             #         self.nChannels = None
             #         self.nSamples = None
             #         self.dataBlocksPerFile = None
             #         self.comments = ''
             #
                 def getltctime(self):
                     if self.useLocalTime:
                         return self.utctime - self.timeZone * 60
                     return self.utctime
                 def getDatatime(self):
                     datatime = datetime.datetime.utcfromtimestamp(self.ltctime)
                     return datatime
                 def getTimeRange(self):
                     datatime = []
                     datatime.append(self.ltctime)
                     datatime.append(self.ltctime + self.timeInterval)
                     datatime = numpy.array(datatime)
                     return datatime
                 def getHeiRange(self):
                     heis = self.heightList
                     return heis
                 def getNHeights(self):
                     return len(self.heightList)
                 def getNChannels(self):
                     return len(self.channelList)
                 def getChannelIndexList(self):
                     return list(range(self.nChannels))
                 def getNoise(self, type=1):
                     #noise = numpy.zeros(self.nChannels)
                     if type == 1:
                         noise = self.getNoisebyHildebrand()
                     if type == 2:
                         noise = self.getNoisebySort()
                     if type == 3:
                         noise = self.getNoisebyWindow()
                     return noise
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt * self.nIncohInt
                     return timeInterval
                 def get_ippSeconds(self):
                     '''
                     '''
                     return self.ipp_sec
                 datatime = property(getDatatime, "I'm the 'datatime' property")
                 nHeights = property(getNHeights, "I'm the 'nHeights' property.")
                 nChannels = property(getNChannels, "I'm the 'nChannel' property.")
                 channelIndexList = property(
                     getChannelIndexList, "I'm the 'channelIndexList' property.")
                 noise = property(getNoise, "I'm the 'nHeights' property.")
                 ltctime = property(getltctime, "I'm the 'ltctime' property")
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
                 ippSeconds = property(get_ippSeconds, '')
             class Correlation(JROData):
                 noise = None
                 SNR = None
                 #--------------------------------------------------
                 mode = None
                 split = False
                 data_cf = None
                 lags = None
                 lagRange = None
                 pairsList = None
                 normFactor = None
                 #--------------------------------------------------
             #     calculateVelocity = None
                 nLags = None
                 nPairs = None
                 nAvg = None
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Correlation"
                     self.data = None
                     self.dtype = None
                     self.nProfiles = None
                     self.heightList = None
                     self.channelList = None
                     self.flagNoData = True
                     self.flagDiscontinuousBlock = False
                     self.utctime = None
                     self.timeZone = None
                     self.dstFlag = None
                     self.errorCount = None
                     self.blocksize = None
                     self.flagDecodeData = False  # asumo q la data no esta decodificada
                     self.flagDeflipData = False  # asumo q la data no esta sin flip
                     self.pairsList = None
                     self.nPoints = None
                 def getPairsList(self):
                     return self.pairsList
                 def getNoise(self, mode=2):
                     indR = numpy.where(self.lagR == 0)[0][0]
                     indT = numpy.where(self.lagT == 0)[0][0]
                     jspectra0 = self.data_corr[:, :, indR, :]
                     jspectra = copy.copy(jspectra0)
                     num_chan = jspectra.shape[0]
                     num_hei = jspectra.shape[2]
                     freq_dc = jspectra.shape[1] / 2
                     ind_vel = numpy.array([-2, -1, 1, 2]) + freq_dc
                     if ind_vel[0] < 0:
                         ind_vel[list(range(0, 1))] = ind_vel[list(
                             range(0, 1))] + self.num_prof
                     if mode == 1:
                         jspectra[:, freq_dc, :] = (
                             jspectra[:, ind_vel[1], :] + jspectra[:, ind_vel[2], :]) / 2  # CORRECCION
                     if mode == 2:
                         vel = numpy.array([-2, -1, 1, 2])
                         xx = numpy.zeros([4, 4])
                         for fil in range(4):
                             xx[fil, :] = vel[fil]**numpy.asarray(list(range(4)))
                         xx_inv = numpy.linalg.inv(xx)
                         xx_aux = xx_inv[0, :]
                         for ich in range(num_chan):
                             yy = jspectra[ich, ind_vel, :]
                             jspectra[ich, freq_dc, :] = numpy.dot(xx_aux, yy)
                             junkid = jspectra[ich, freq_dc, :] <= 0
                             cjunkid = sum(junkid)
                             if cjunkid.any():
                                 jspectra[ich, freq_dc, junkid.nonzero()] = (
                                     jspectra[ich, ind_vel[1], junkid] + jspectra[ich, ind_vel[2], junkid]) / 2
                     noise = jspectra0[:, freq_dc, :] - jspectra[:, freq_dc, :]
                     return noise
                 def getTimeInterval(self):
                     timeInterval = self.ippSeconds * self.nCohInt * self.nProfiles
                     return timeInterval
                 def splitFunctions(self):
                     pairsList = self.pairsList
                     ccf_pairs = []
                     acf_pairs = []
                     ccf_ind = []
                     acf_ind = []
                     for l in range(len(pairsList)):
                         chan0 = pairsList[l][0]
                         chan1 = pairsList[l][1]
                         # Obteniendo pares de Autocorrelacion
                         if chan0 == chan1:
                             acf_pairs.append(chan0)
                             acf_ind.append(l)
                         else:
                             ccf_pairs.append(pairsList[l])
                             ccf_ind.append(l)
                     data_acf = self.data_cf[acf_ind]
                     data_ccf = self.data_cf[ccf_ind]
                     return acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf
                 def getNormFactor(self):
                     acf_ind, ccf_ind, acf_pairs, ccf_pairs, data_acf, data_ccf = self.splitFunctions()
                     acf_pairs = numpy.array(acf_pairs)
                     normFactor = numpy.zeros((self.nPairs, self.nHeights))
                     for p in range(self.nPairs):
                         pair = self.pairsList[p]
                         ch0 = pair[0]
                         ch1 = pair[1]
                         ch0_max = numpy.max(data_acf[acf_pairs == ch0, :, :], axis=1)
                         ch1_max = numpy.max(data_acf[acf_pairs == ch1, :, :], axis=1)
                         normFactor[p, :] = numpy.sqrt(ch0_max * ch1_max)
                     return normFactor
                 timeInterval = property(getTimeInterval, "I'm the 'timeInterval' property")
                 normFactor = property(getNormFactor, "I'm the 'normFactor property'")
             class Parameters(Spectra):
                 experimentInfo = None  # Information about the experiment
                 # Information from previous data
                 inputUnit = None  # Type of data to be processed
                 operation = None  # Type of operation to parametrize
                 # normFactor = None       #Normalization Factor
                 groupList = None  # List of Pairs, Groups, etc
                 # Parameters
                 data_param = None  # Parameters obtained
                 data_pre = None  # Data Pre Parametrization
                 data_SNR = None  # Signal to Noise Ratio
             #     heightRange = None      #Heights
                 abscissaList = None  # Abscissa, can be velocities, lags or time
             #    noise = None            #Noise Potency
                 utctimeInit = None  # Initial UTC time
                 paramInterval = None  # Time interval to calculate Parameters in seconds
                 useLocalTime = True
                 # Fitting
                 data_error = None  # Error of the estimation
                 constants = None
                 library = None
                 # Output signal
                 outputInterval = None  # Time interval to calculate output signal in seconds
                 data_output = None  # Out signal
                 nAvg = None
                 noise_estimation = None
                 GauSPC = None  # Fit gaussian SPC
                 def __init__(self):
                     '''
                     Constructor
                     '''
                     self.radarControllerHeaderObj = RadarControllerHeader()
                     self.systemHeaderObj = SystemHeader()
                     self.type = "Parameters"
                 def getTimeRange1(self, interval):
                     datatime = []
                     if self.useLocalTime:
                         time1 = self.utctimeInit - self.timeZone * 60
                     else:
                         time1 = self.utctimeInit
                     datatime.append(time1)
                     datatime.append(time1 + interval)
                     datatime = numpy.array(datatime)
                     return datatime
                 def getTimeInterval(self):
                     if hasattr(self, 'timeInterval1'):
                         return self.timeInterval1
                     else:
                         return self.paramInterval
                 def setValue(self, value):
                     print("This property should not be initialized")
                     return
                 def getNoise(self):
                     return self.spc_noise
                 timeInterval = property(getTimeInterval)
                 noise = property(getNoise, setValue, "I'm the 'Noise' property.")
             class PlotterData(object):
                 '''
                 Object to hold data to be plotted
                 '''
                 MAXNUMX = 100
                 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
                     self.buffering = buffering
                     self.ready = False
                     self.flagNoData = False
                     self.localtime = False
                     self.data = {}
                     self.meta = {}
                     self.__times = []
                     self.__heights = []
                     if 'snr' in code:
                         self.plottypes = ['snr']
                     elif code == 'spc':
                         self.plottypes = ['spc', 'noise', 'rti']
                     elif code == 'rti':
                         self.plottypes = ['noise', 'rti']
                     else:
                         self.plottypes = [code]
                     if 'snr' not in self.plottypes and snr:
                         self.plottypes.append('snr')
                     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))
                 def __len__(self):
                     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:
                         ret = self.data[key]
                     elif 'scope' in key:
                         ret = numpy.array(self.data[key][float(self.tm)])
                     else:
                         ret = numpy.array([self.data[key][x] for x in self.times])
                         if ret.ndim > 1:
                             ret = numpy.swapaxes(ret, 0, 1)
                     return ret
                 def __contains__(self, key):
                     return key in self.data
                 def setup(self):
                     '''
                     Configure object
                     '''
                     self.type = ''
                     self.ready = False
                     self.data = {}
                     self.__times = []
                     self.__heights = []
                     self.__all_heights = set()
                     for plot in self.plottypes:
                         if 'snr' in plot:
                             plot = 'snr'
                         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'] = {}
                         if 'noise' not in self.plottypes:
                             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
                     '''
                     if len(self.data[key]):
                         if 'spc' in key or not self.buffering:
                             return self.data[key].shape
                         return self.data[key][self.__times[0]].shape
                     return (0,)
                 def update(self, dataOut, tm):
                     '''
                     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]:
                         self.xrange = (dataOut.getFreqRange(1)/1000.,
                                        dataOut.getAcfRange(1), dataOut.getVelRange(1))
                         self.factor = dataOut.normFactor
                     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
                             buffer = 10*numpy.log10(z)
                         if plot == 'cspc':
                             z = dataOut.data_spc/dataOut.normFactor
                             buffer = (dataOut.data_spc, dataOut.data_cspc)
                         if plot == 'noise':
                             buffer = 10*numpy.log10(dataOut.getNoise()/dataOut.normFactor)
                         if plot in ('rti', 'spcprofile'):
                             buffer = dataOut.getPower()
                         if plot == 'snr_db':
                             buffer = dataOut.data_SNR
                         if plot == 'snr':
                             buffer = 10*numpy.log10(dataOut.data_SNR)
                         if plot == 'dop':
                             buffer = dataOut.data_DOP
                         if plot == 'pow':
                             buffer = 10*numpy.log10(dataOut.data_POW)
                         if plot == 'width':
                             buffer = dataOut.data_WIDTH
                         if plot == 'coh':
                             buffer = dataOut.getCoherence()
                         if plot == 'phase':
                             buffer = dataOut.getCoherence(phase=True)
                         if plot == 'output':
                             buffer = dataOut.data_output
                         if plot == 'param':
                             buffer = dataOut.data_param
                         if plot == 'scope':
                             buffer = dataOut.data
                             self.flagDataAsBlock = dataOut.flagDataAsBlock
                             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':
                             self.data['spc'] = buffer[0]
                             self.data['cspc'] = buffer[1]
                         elif plot == 'spc_moments':
                             self.data['spc'] = buffer
                             self.data['moments'][tm] = dataOut.moments
                         else:
                             if self.buffering:
                                 self.data[plot][tm] = buffer
                             else:
                                 self.data[plot] = buffer
                     if dataOut.channelList is None:
                         self.channels = range(buffer.shape[0])
                     else:
                         self.channels = dataOut.channelList
                     if buffer is None:
                         self.flagNoData = True
                         raise schainpy.admin.SchainWarning('Attribute data_{} is empty'.format(self.key))
                 def normalize_heights(self):
                     '''
                     Ensure same-dimension of the data for different heighList
                     '''
                     H = numpy.array(list(self.__all_heights))
                     H.sort()
                     for key in self.data:
                         shape = self.shape(key)[:-1] + H.shape
                         for tm, obj in list(self.data[key].items()):
                             h = self.__heights[self.__times.index(tm)]
                             if H.size == h.size:
                                 continue
                             index = numpy.where(numpy.in1d(H, h))[0]
                             dummy = numpy.zeros(shape) + numpy.nan
                             if len(shape) == 2:
                                 dummy[:, index] = obj
                             else:
                                 dummy[index] = obj
                             self.data[key][tm] = dummy
                     self.__heights = [H for tm in self.__times]
                 def jsonify(self, plot_name, plot_type, decimate=False):
                     '''
                     Convert data to json
                     '''
                     tm = self.times[-1]
                     dy = int(self.heights.size/self.MAXNUMY) + 1
                     if self.key in ('spc', 'cspc') or not self.buffering:
                         dx = int(self.data[self.key].shape[1]/self.MAXNUMX) + 1
                         data = self.roundFloats(
                             self.data[self.key][::, ::dx, ::dy].tolist())
                     else:
                         data = self.roundFloats(self.data[self.key][tm].tolist())
                         if self.key is 'noise':
                             data = [[x] for x in data]
                     meta = {}
                     ret = {
                         'plot': plot_name,
                         'code': self.exp_code,
                         'time': float(tm),
                         'data': data,
                         }
                     meta['type'] = plot_type
                     meta['interval'] = float(self.interval)
                     meta['localtime'] = self.localtime
                     meta['yrange'] = self.roundFloats(self.heights[::dy].tolist())
                     if 'spc' in self.data or 'cspc' in self.data:
                         meta['xrange'] = self.roundFloats(self.xrange[2][::dx].tolist())
                     else:
                         meta['xrange'] = []
                     meta.update(self.meta)
                     ret['metadata'] = meta
                     return json.dumps(ret)
                 @property
                 def times(self):
                     '''
                     Return the list of times of the current data
                     '''
                     ret = numpy.array(self.__times)
                     ret.sort()
                     return ret
                 @property
                 def min_time(self):
                     '''
                     Return the minimun time value
                     '''
                     return self.times[0]
                 @property
                 def max_time(self):
                     '''
                     Return the maximun time value
                     '''
                     return self.times[-1]
                 @property
                 def heights(self):
                     '''
                     Return the list of heights of the current data
                     '''
                     return numpy.array(self.__heights[-1])
                 @staticmethod
                 def roundFloats(obj):
                     if isinstance(obj, list):
                         return list(map(PlotterData.roundFloats, obj))
                     elif isinstance(obj, float):
                         return round(obj, 2)

schainpy/model/graphics/jroplot_voltage.py +147 -48

@@ -1,134 +1,233
1	'''	1	'''
2	Created on Jul 9, 2014	2	Created on Jul 9, 2014
3		3
4	@author: roj-idl71	4	@author: roj-idl71
5	'''	5	'''
6	import os	6	import os
7	import datetime	7	import datetime
8	import numpy	8	import numpy
9		9
10	from schainpy.model.graphics.jroplot_base import Plot, plt	10	from schainpy.model.graphics.jroplot_base import Plot, plt
11		11
12		12
13	class ScopePlot(Plot):	13	class ScopePlot(Plot):
14		14
15	'''	15	'''
16	Plot for Scope	16	Plot for Scope
17	'''	17	'''
18		18
19	CODE = 'scope'	19	CODE = 'scope'
20	plot_name = 'Scope'	20	plot_name = 'Scope'
21	plot_type = 'scatter'	21	plot_type = 'scatter'
22		22
23	def setup(self):	23	def setup(self):
24		24
25	self.xaxis = 'Range (Km)'	25	self.xaxis = 'Range (Km)'
26	self.ncols = 1	26	self.ncols = 1
27	self.nrows = 1	27	self.nrows = 1
28	self.nplots = 1	28	self.nplots = 1
29	self.ylabel = 'Intensity [dB]'	29	self.ylabel = 'Intensity [dB]'
30	self.titles = ['Scope']	30	self.titles = ['Scope']
31	self.colorbar = False	31	self.colorbar = False
32	self.width = 6	32	self.width = 6
33	self.height = 4	33	self.height = 4
34		34
35	def plot_iq(self, x, y, channelIndexList, thisDatetime, wintitle):	35	def plot_iq(self, x, y, channelIndexList, thisDatetime, wintitle):
36		36
37	yreal = y[channelIndexList,:].real	37	yreal = y[channelIndexList,:].real
38	yimag = y[channelIndexList,:].imag	38	yimag = y[channelIndexList,:].imag
39	title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))	39	title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
40	self.xlabel = "Range (Km)"	40	self.xlabel = "Range (Km)"
41	self.ylabel = "Intensity - IQ"	41	self.ylabel = "Intensity - IQ"
42		42
43	self.y = yreal	43	self.y = yreal
44	self.x = x	44	self.x = x
45	self.xmin = min(x)	45	self.xmin = min(x)
46	self.xmax = max(x)	46	self.xmax = max(x)
47
48		47
49	self.titles[0] = title	48
		49	self.titles[0] = title
50		50
51	for i,ax in enumerate(self.axes):	51	for i,ax in enumerate(self.axes):
52	title = "Channel %d" %(i)	52	title = "Channel %d" %(i)
53	if ax.firsttime:	53	if ax.firsttime:
54	ax.plt_r = ax.plot(x, yreal[i,:], color='b')[0]	54	ax.plt_r = ax.plot(x, yreal[i,:], color='b')[0]
55	ax.plt_i = ax.plot(x, yimag[i,:], color='r')[0]	55	ax.plt_i = ax.plot(x, yimag[i,:], color='r')[0]
56	else:	56	else:
57	ax.plt_r.set_data(x, yreal[i,:])	57	ax.plt_r.set_data(x, yreal[i,:])
58	ax.plt_i.set_data(x, yimag[i,:])	58	ax.plt_i.set_data(x, yimag[i,:])
59		59
60	def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle):	60	def plot_power(self, x, y, channelIndexList, thisDatetime, wintitle):
61	y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])	61	y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:])
62	yreal = y.real	62	yreal = y.real
		63	yreal = 10*numpy.log10(yreal)
63	self.y = yreal	64	self.y = yreal
64	title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))	65	title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y"))
65	self.xlabel = "Range (Km)"	66	self.xlabel = "Range (Km)"
66	self.ylabel = "Intensity"	67	self.ylabel = "Intensity"
67	self.xmin = min(x)	68	self.xmin = min(x)
68	self.xmax = max(x)	69	self.xmax = max(x)
69		70
70		71
71	self.titles[0] = title	72	self.titles[0] = title
72		73
73	for i,ax in enumerate(self.axes):	74	for i,ax in enumerate(self.axes):
74	title = "Channel %d" %(i)	75	title = "Channel %d" %(i)
75		76
76	ychannel = yreal[i,:]	77	ychannel = yreal[i,:]
77		78
78	if ax.firsttime:	79	if ax.firsttime:
79	ax.plt_r = ax.plot(x, ychannel)[0]	80	ax.plt_r = ax.plot(x, ychannel)[0]
80	else:	81	else:
81	#pass	82	#pass
82	ax.plt_r.set_data(x, ychannel)	83	ax.plt_r.set_data(x, ychannel)
83		84
		85	def plot_weatherpower(self, x, y, channelIndexList, thisDatetime, wintitle):
		86
		87
		88	y = y[channelIndexList,:]
		89	yreal = y.real
		90	yreal = 10*numpy.log10(yreal)
		91	self.y = yreal
		92	title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
		93	self.xlabel = "Range (Km)"
		94	self.ylabel = "Intensity"
		95	self.xmin = min(x)
		96	self.xmax = max(x)
		97
		98	self.titles[0] =title
		99	for i,ax in enumerate(self.axes):
		100	title = "Channel %d" %(i)
		101
		102	ychannel = yreal[i,:]
		103
		104	if ax.firsttime:
		105	ax.plt_r = ax.plot(x, ychannel)[0]
		106	else:
		107	#pass
		108	ax.plt_r.set_data(x, ychannel)
		109
		110	def plot_weathervelocity(self, x, y, channelIndexList, thisDatetime, wintitle):
		111
		112	x = x[channelIndexList,:]
		113	yreal = y
		114	self.y = yreal
		115	title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S"))
		116	self.xlabel = "Velocity (m/s)"
		117	self.ylabel = "Range (Km)"
		118	self.xmin = numpy.min(x)
		119	self.xmax = numpy.max(x)
		120	self.titles[0] =title
		121	for i,ax in enumerate(self.axes):
		122	title = "Channel %d" %(i)
		123	xchannel = x[i,:]
		124	if ax.firsttime:
		125	ax.plt_r = ax.plot(xchannel, yreal)[0]
		126	else:
		127	#pass
		128	ax.plt_r.set_data(xchannel, yreal)
84		129
85	def plot(self):	130	def plot(self):
86
87	if self.channels:	131	if self.channels:
88	channels = self.channels	132	channels = self.channels
89	else:	133	else:
90	channels = self.data.channels	134	channels = self.data.channels
91		135
92	thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])	136	thisDatetime = datetime.datetime.utcfromtimestamp(self.data.times[-1])
93		137	if self.CODE == "pp_power":
94	scope = self.data['~~scope~~']	138	scope = self.data['pp_power']
95		139	elif self.CODE == "pp_velocity":
		140	scope = self.data["pp_velocity"]
		141	else:
		142	scope =self.data["scope"]
		143
96	if self.data.flagDataAsBlock:	144	if self.data.flagDataAsBlock:
97		145
98	for i in range(self.data.nProfiles):	146	for i in range(self.data.nProfiles):
99		147
100	wintitle1 = " [Profile = %d] " %i	148	wintitle1 = " [Profile = %d] " %i
101		149	if self.CODE =="scope":
102	if self.type == "power":	150	if self.type == "power":
103	self.plot_power(self.data.heights,	151	self.plot_power(self.data.heights,
104	scope[:,i,:],	152	scope[:,i,:],
105	channels,	153	channels,
106	thisDatetime,	154	thisDatetime,
107	wintitle1	155	wintitle1
		156	)
		157
		158	if self.type == "iq":
		159	self.plot_iq(self.data.heights,
		160	scope[:,i,:],
		161	channels,
		162	thisDatetime,
		163	wintitle1
108	)	164	)
		165	if self.CODE=="pp_power":
		166	self.plot_weatherpower(self.data.heights,
		167	scope[:,i,:],
		168	channels,
		169	thisDatetime,
		170	wintitle
		171	)
		172	if self.CODE=="pp_velocity":
		173	self.plot_weathervelocity(scope[:,i,:],
		174	self.data.heights,
		175	channels,
		176	thisDatetime,
		177	wintitle
		178	)
		179	else:
		180	wintitle = " [Profile = %d] " %self.data.profileIndex
		181	if self.CODE== "scope":
		182	if self.type == "power":
		183	self.plot_power(self.data.heights,
		184	scope,
		185	channels,
		186	thisDatetime,
		187	wintitle
		188	)
109		189
110	if self.type == "iq":	190	if self.type == "iq":
111	self.plot_iq(self.data.heights,	191	self.plot_iq(self.data.heights,
112	scope~~[:,~~i~~,:]~~,	192	scope,
113	channels,	193	channels,
114	thisDatetime,	194	thisDatetime,
115	wintitle1	195	wintitle
116	)	196	)
117	else:	197	if self.CODE=="pp_power":
118	wintitle = " [Profile = %d] " %self.data.profileIndex	198	self.plot_weatherpower(self.data.heights,
119		199	scope,
120	if self.type == "power":	200	channels,
121	self.plot_power(self.data.heights,	201	thisDatetime,
122	~~scope~~,	202	wintitle
123	~~channels~~,	203	)
124	thisDatetime,	204	if self.CODE=="pp_velocity":
125	wintitle	205	self.plot_weathervelocity(scope,
126	)	206	self.data.heights,
127		207	channels,
128	if self.type == "iq":	208	thisDatetime,
129	self.plot_iq(self.data.heights,	209	wintitle
130	~~scope~~,	210	)
131	channels,	211
132	thisDatetime,	212
133	wintitle	213
134	)	214	class PulsepairPowerPlot(ScopePlot):
		215	'''
		216	Plot for
		217	'''
		218
		219	CODE = 'pp_power'
		220	plot_name = 'PulsepairPower'
		221	plot_type = 'scatter'
		222	buffering = False
		223
		224
		225
		226	class PulsepairVelocityPlot(ScopePlot):
		227	'''
		228	Plot for
		229	'''
		230	CODE = 'pp_velocity'
		231	plot_name = 'PulsepairVelocity'
		232	plot_type = 'scatter'
		233	buffering = False

schainpy/model/proc/jroproc_voltage.py +173 -25

             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
             from schainpy.utils import log
             from time import time
             class VoltageProc(ProcessingUnit):
                 def __init__(self):
                     ProcessingUnit.__init__(self)
                     self.dataOut = Voltage()
                     self.flip = 1
                     self.setupReq = False
                 def run(self):
                     if self.dataIn.type == 'AMISR':
                         self.__updateObjFromAmisrInput()
                     if self.dataIn.type == 'Voltage':
                         self.dataOut.copy(self.dataIn)
                 def __updateObjFromAmisrInput(self):
                     self.dataOut.timeZone = self.dataIn.timeZone
                     self.dataOut.dstFlag = self.dataIn.dstFlag
                     self.dataOut.errorCount = self.dataIn.errorCount
                     self.dataOut.useLocalTime = self.dataIn.useLocalTime
                     self.dataOut.flagNoData = self.dataIn.flagNoData
                     self.dataOut.data = self.dataIn.data
                     self.dataOut.utctime = self.dataIn.utctime
                     self.dataOut.channelList = self.dataIn.channelList
                     #self.dataOut.timeInterval = self.dataIn.timeInterval
                     self.dataOut.heightList = self.dataIn.heightList
                     self.dataOut.nProfiles = self.dataIn.nProfiles
                     self.dataOut.nCohInt = self.dataIn.nCohInt
                     self.dataOut.ippSeconds = self.dataIn.ippSeconds
                     self.dataOut.frequency = self.dataIn.frequency
                     self.dataOut.azimuth = self.dataIn.azimuth
                     self.dataOut.zenith = self.dataIn.zenith
                     self.dataOut.beam.codeList = self.dataIn.beam.codeList
                     self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList
                     self.dataOut.beam.zenithList = self.dataIn.beam.zenithList
             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
                 def selectChannelsByIndex(self, channelIndexList):
                     """
                     Selecciona un bloque de datos en base a canales segun el channelIndexList
                     Input:
                         channelIndexList    :    lista sencilla de canales a seleccionar por ej. [2,3,7]
                     Affected:
                         self.dataOut.data
                         self.dataOut.channelIndexList
                         self.dataOut.nChannels
                         self.dataOut.m_ProcessingHeader.totalSpectra
                         self.dataOut.systemHeaderObj.numChannels
                         self.dataOut.m_ProcessingHeader.blockSize
                     Return:
                         None
                     """
                     for channelIndex in channelIndexList:
                         if channelIndex not in self.dataOut.channelIndexList:
                             raise ValueError("The value %d in channelIndexList is not valid" %channelIndex)
                     if self.dataOut.type == 'Voltage':
                         if self.dataOut.flagDataAsBlock:
                             """
                             Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
                             """
                             data = self.dataOut.data[channelIndexList,:,:]
                         else:
                             data = self.dataOut.data[channelIndexList,:]
                         self.dataOut.data = data
                         # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
                         self.dataOut.channelList = range(len(channelIndexList))
                     elif self.dataOut.type == 'Spectra':
                         data_spc = self.dataOut.data_spc[channelIndexList, :]
                         data_dc = self.dataOut.data_dc[channelIndexList, :]
                         self.dataOut.data_spc = data_spc
                         self.dataOut.data_dc = data_dc
                         # self.dataOut.channelList = [self.dataOut.channelList[i] for i in channelIndexList]
                         self.dataOut.channelList = range(len(channelIndexList))
                         self.__selectPairsByChannel(channelIndexList)
                     return 1
                 def __selectPairsByChannel(self, channelList=None):
                     if channelList == None:
                         return
                     pairsIndexListSelected = []
                     for pairIndex in self.dataOut.pairsIndexList:
                         # First pair
                         if self.dataOut.pairsList[pairIndex][0] not in channelList:
                             continue
                         # Second pair
                         if self.dataOut.pairsList[pairIndex][1] not in channelList:
                             continue
                         pairsIndexListSelected.append(pairIndex)
                     if not pairsIndexListSelected:
                         self.dataOut.data_cspc = None
                         self.dataOut.pairsList = []
                         return
                     self.dataOut.data_cspc = self.dataOut.data_cspc[pairsIndexListSelected]
                     self.dataOut.pairsList = [self.dataOut.pairsList[i]
                                               for i in pairsIndexListSelected]
                     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
                     minHei <= height <= maxHei
                     Input:
                         minHei    :    valor minimo de altura a considerar
                         maxHei    :    valor maximo de altura a considerar
                     Affected:
                         Indirectamente son cambiados varios valores a travez del metodo selectHeightsByIndex
                     Return:
 si el metodo se ejecuto con exito caso contrario devuelve 0
                     """
                     self.dataOut = dataOut
                     if minHei == None:
                         minHei = self.dataOut.heightList[0]
                     if maxHei == None:
                         maxHei = self.dataOut.heightList[-1]
                     if (minHei < self.dataOut.heightList[0]):
                         minHei = self.dataOut.heightList[0]
                     if (maxHei > self.dataOut.heightList[-1]):
                         maxHei = self.dataOut.heightList[-1]
                     minIndex = 0
                     maxIndex = 0
                     heights = self.dataOut.heightList
                     inda = numpy.where(heights >= minHei)
                     indb = numpy.where(heights <= maxHei)
                     try:
                         minIndex = inda[0][0]
                     except:
                         minIndex = 0
                     try:
                         maxIndex = indb[0][-1]
                     except:
                         maxIndex = len(heights)
                     self.selectHeightsByIndex(minIndex, maxIndex)
                     return self.dataOut
                 def selectHeightsByIndex(self, minIndex, maxIndex):
                     """
                     Selecciona un bloque de datos en base a un grupo indices de alturas segun el rango
                     minIndex <= index <= maxIndex
                     Input:
                         minIndex    :    valor de indice minimo de altura a considerar
                         maxIndex    :    valor de indice maximo de altura a considerar
                     Affected:
                         self.dataOut.data
                         self.dataOut.heightList
                     Return:
 si el metodo se ejecuto con exito caso contrario devuelve 0
                     """
                     if self.dataOut.type == 'Voltage':
                         if (minIndex < 0) or (minIndex > maxIndex):
                             raise ValueError("Height index range (%d,%d) is not valid" % (minIndex, maxIndex))
                         if (maxIndex >= self.dataOut.nHeights):
                             maxIndex = self.dataOut.nHeights
                         #voltage
                         if self.dataOut.flagDataAsBlock:
                             """
                             Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
                             """
                             data = self.dataOut.data[:,:, minIndex:maxIndex]
                         else:
                             data = self.dataOut.data[:, minIndex:maxIndex]
                         #         firstHeight = self.dataOut.heightList[minIndex]
                         self.dataOut.data = data
                         self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex]
                         if self.dataOut.nHeights <= 1:
                             raise ValueError("selectHeights: Too few heights. Current number of heights is %d" %(self.dataOut.nHeights))
                     elif self.dataOut.type == 'Spectra':
                         if (minIndex < 0) or (minIndex > maxIndex):
                             raise ValueError("Error selecting heights: Index range (%d,%d) is not valid" % (
                                 minIndex, maxIndex))
                         if (maxIndex >= self.dataOut.nHeights):
                             maxIndex = self.dataOut.nHeights - 1
                         # Spectra
                         data_spc = self.dataOut.data_spc[:, :, minIndex:maxIndex + 1]
                         data_cspc = None
                         if self.dataOut.data_cspc is not None:
                             data_cspc = self.dataOut.data_cspc[:, :, minIndex:maxIndex + 1]
                         data_dc = None
                         if self.dataOut.data_dc is not None:
                             data_dc = self.dataOut.data_dc[:, minIndex:maxIndex + 1]
                         self.dataOut.data_spc = data_spc
                         self.dataOut.data_cspc = data_cspc
                         self.dataOut.data_dc = data_dc
                         self.dataOut.heightList = self.dataOut.heightList[minIndex:maxIndex + 1]
                     return 1
             class filterByHeights(Operation):
                 def run(self, dataOut, window):
                     deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                     if window == None:
                         window = (dataOut.radarControllerHeaderObj.txA/dataOut.radarControllerHeaderObj.nBaud) / deltaHeight
                     newdelta = deltaHeight * window
                     r = dataOut.nHeights % window
                     newheights = (dataOut.nHeights-r)/window
                     if newheights <= 1:
                         raise ValueError("filterByHeights: Too few heights. Current number of heights is %d and window is %d" %(dataOut.nHeights, window))
                     if dataOut.flagDataAsBlock:
                         """
                         Si la data es obtenida por bloques, dimension = [nChannels, nProfiles, nHeis]
                         """
                         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)
                     else:
                         buffer = dataOut.data[:,0:int(dataOut.nHeights-r)]
                         buffer = buffer.reshape(dataOut.nChannels,int(dataOut.nHeights/window),int(window))
                         buffer = numpy.sum(buffer,2)
                     dataOut.data = buffer
                     dataOut.heightList = dataOut.heightList[0] + numpy.arange( newheights )*newdelta
                     dataOut.windowOfFilter = window
                     return dataOut
             class setH0(Operation):
                 def run(self, dataOut, h0, deltaHeight = None):
                     if not deltaHeight:
                         deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                     nHeights = dataOut.nHeights
                     newHeiRange = h0 + numpy.arange(nHeights)*deltaHeight
                     dataOut.heightList = newHeiRange
                     return dataOut
             class deFlip(Operation):
                 def run(self, dataOut, channelList = []):
                     data = dataOut.data.copy()
                     if dataOut.flagDataAsBlock:
                         flip = self.flip
                         profileList = list(range(dataOut.nProfiles))
                         if not channelList:
                             for thisProfile in profileList:
                                 data[:,thisProfile,:] = data[:,thisProfile,:]*flip
                                 flip *= -1.0
                         else:
                             for thisChannel in channelList:
                                 if thisChannel not in dataOut.channelList:
                                     continue
                                 for thisProfile in profileList:
                                     data[thisChannel,thisProfile,:] = data[thisChannel,thisProfile,:]*flip
                                     flip *= -1.0
                         self.flip = flip
                     else:
                         if not channelList:
                             data[:,:] = data[:,:]*self.flip
                         else:
                             for thisChannel in channelList:
                                 if thisChannel not in dataOut.channelList:
                                     continue
                                 data[thisChannel,:] = data[thisChannel,:]*self.flip
                         self.flip *= -1.
                     dataOut.data = data
                     return dataOut
             class setAttribute(Operation):
                 '''
                 Set an arbitrary attribute(s) to dataOut
                 '''
                 def __init__(self):
                     Operation.__init__(self)
                     self._ready = False
                 def run(self, dataOut, **kwargs):
                     for key, value in kwargs.items():
                         setattr(dataOut, key, value)
                     return dataOut
             class interpolateHeights(Operation):
                 def run(self, dataOut, topLim, botLim):
                     #69 al 72 para julia
                     #82-84 para meteoros
                     if len(numpy.shape(dataOut.data))==2:
                         sampInterp = (dataOut.data[:,botLim-1] + dataOut.data[:,topLim+1])/2
                         sampInterp = numpy.transpose(numpy.tile(sampInterp,(topLim-botLim + 1,1)))
                         #dataOut.data[:,botLim:limSup+1] = sampInterp
                         dataOut.data[:,botLim:topLim+1] = sampInterp
                     else:
                         nHeights = dataOut.data.shape[2]
                         x = numpy.hstack((numpy.arange(botLim),numpy.arange(topLim+1,nHeights)))
                         y = dataOut.data[:,:,list(range(botLim))+list(range(topLim+1,nHeights))]
                         f = interpolate.interp1d(x, y, axis = 2)
                         xnew = numpy.arange(botLim,topLim+1)
                         ynew = f(xnew)
                         dataOut.data[:,:,botLim:topLim+1]  = ynew
                     return dataOut
             class CohInt(Operation):
                 isConfig = False
                 __profIndex = 0
                 __byTime = False
                 __initime = None
                 __lastdatatime = None
                 __integrationtime = None
                 __buffer = None
                 __bufferStride = []
                 __dataReady = False
                 __profIndexStride = 0
                 __dataToPutStride = False
                 n = None
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     #   self.isConfig = False
                 def setup(self, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False):
                     """
                     Set the parameters of the integration class.
                     Inputs:
                         n               :    Number of coherent integrations
                         timeInterval    :    Time of integration. If the parameter "n" is selected this one does not work
                         overlapping     :
                     """
                     self.__initime = None
                     self.__lastdatatime = 0
                     self.__buffer = None
                     self.__dataReady = False
                     self.byblock = byblock
                     self.stride = stride
                     if n == None and timeInterval == None:
                         raise ValueError("n or timeInterval should be specified ...")
                     if n != None:
                         self.n = n
                         self.__byTime = False
                     else:
                         self.__integrationtime = timeInterval #* 60. #if (type(timeInterval)!=integer) -> change this line
                         self.n = 9999
                         self.__byTime = True
                     if overlapping:
                         self.__withOverlapping = True
                         self.__buffer = None
                     else:
                         self.__withOverlapping = False
                         self.__buffer = 0
                     self.__profIndex = 0
                 def putData(self, data):
                     """
                     Add a profile to the __buffer and increase in one the __profileIndex
                     """
                     if not self.__withOverlapping:
                         self.__buffer += data.copy()
                         self.__profIndex += 1
                         return
                     #Overlapping data
                     nChannels, nHeis = data.shape
                     data = numpy.reshape(data, (1, nChannels, nHeis))
                     #If the buffer is empty then it takes the data value
                     if self.__buffer is None:
                         self.__buffer = data
                         self.__profIndex += 1
                         return
                     #If the buffer length is lower than n then stakcing the data value
                     if self.__profIndex < self.n:
                         self.__buffer = numpy.vstack((self.__buffer, data))
                         self.__profIndex += 1
                         return
                     #If the buffer length is equal to n then replacing the last buffer value with the data value
                     self.__buffer = numpy.roll(self.__buffer, -1, axis=0)
                     self.__buffer[self.n-1] = data
                     self.__profIndex = self.n
                     return
                 def pushData(self):
                     """
                     Return the sum of the last profiles and the profiles used in the sum.
                     Affected:
                     self.__profileIndex
                     """
                     if not self.__withOverlapping:
                         data = self.__buffer
                         n = self.__profIndex
                         self.__buffer = 0
                         self.__profIndex = 0
                         return data, n
                     #Integration with Overlapping
                     data = numpy.sum(self.__buffer, axis=0)
                     # print data
                     # raise
                     n = self.__profIndex
                     return data, n
                 def byProfiles(self, data):
                     self.__dataReady = False
                     avgdata = None
                     #         n = None
                     # print data
                     # raise
                     self.putData(data)
                     if self.__profIndex == self.n:
                         avgdata, n = self.pushData()
                         self.__dataReady = True
                     return avgdata
                 def byTime(self, data, datatime):
                     self.__dataReady = False
                     avgdata = None
                     n = None
                     self.putData(data)
                     if (datatime - self.__initime) >= self.__integrationtime:
                         avgdata, n = self.pushData()
                         self.n = n
                         self.__dataReady = True
                     return avgdata
                 def integrateByStride(self, data, datatime):
                     # print data
                     if self.__profIndex == 0:
                         self.__buffer = [[data.copy(), datatime]]
                     else:
                         self.__buffer.append([data.copy(),datatime])
                     self.__profIndex += 1
                     self.__dataReady = False
                     if self.__profIndex == self.n * self.stride :
                         self.__dataToPutStride = True
                         self.__profIndexStride = 0
                         self.__profIndex = 0
                         self.__bufferStride = []
                         for i in range(self.stride):
                             current = self.__buffer[i::self.stride]
                             data = numpy.sum([t[0] for t in current], axis=0)
                             avgdatatime = numpy.average([t[1] for t in current])
                             # print data
                             self.__bufferStride.append((data, avgdatatime))
                     if self.__dataToPutStride:
                         self.__dataReady = True
                         self.__profIndexStride += 1
                         if self.__profIndexStride == self.stride:
                             self.__dataToPutStride = False
                         # print self.__bufferStride[self.__profIndexStride - 1]
                         # raise
                         return self.__bufferStride[self.__profIndexStride - 1]
                     return None, None
                 def integrate(self, data, datatime=None):
                     if self.__initime == None:
                         self.__initime = datatime
                     if self.__byTime:
                         avgdata = self.byTime(data, datatime)
                     else:
                         avgdata = self.byProfiles(data)
                     self.__lastdatatime = datatime
                     if avgdata is None:
                         return None, None
                     avgdatatime = self.__initime
                     deltatime = datatime - self.__lastdatatime
                     if not self.__withOverlapping:
                         self.__initime = datatime
                     else:
                         self.__initime += deltatime
                     return avgdata, avgdatatime
                 def integrateByBlock(self, dataOut):
                     times = int(dataOut.data.shape[1]/self.n)
                     avgdata = numpy.zeros((dataOut.nChannels, times, dataOut.nHeights), dtype=numpy.complex)
                     id_min = 0
                     id_max = self.n
                     for i in range(times):
                         junk = dataOut.data[:,id_min:id_max,:]
                         avgdata[:,i,:] = junk.sum(axis=1)
                         id_min += self.n
                         id_max += self.n
                     timeInterval = dataOut.ippSeconds*self.n
                     avgdatatime = (times - 1) * timeInterval + dataOut.utctime
                     self.__dataReady = True
                     return avgdata, avgdatatime
                 def run(self, dataOut, n=None, timeInterval=None, stride=None, overlapping=False, byblock=False, **kwargs):
                     if not self.isConfig:
                         self.setup(n=n, stride=stride, timeInterval=timeInterval, overlapping=overlapping, byblock=byblock, **kwargs)
                         self.isConfig = True
                     if dataOut.flagDataAsBlock:
                         """
                         Si la data es leida por bloques, dimension = [nChannels, nProfiles, nHeis]
                         """
                         avgdata, avgdatatime = self.integrateByBlock(dataOut)
                         dataOut.nProfiles /= self.n
                     else:
                         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
                     if self.__dataReady:
                         dataOut.data = avgdata
                         dataOut.nCohInt *= self.n
                         dataOut.utctime = avgdatatime
                         # print avgdata, avgdatatime
                         # raise
                         #   dataOut.timeInterval = dataOut.ippSeconds * dataOut.nCohInt
                         dataOut.flagNoData = False
                     return dataOut
             class Decoder(Operation):
                 isConfig = False
                 __profIndex = 0
                 code = None
                 nCode = None
                 nBaud = None
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.times = None
                     self.osamp = None
                 #         self.__setValues = False
                     self.isConfig = False
                     self.setupReq = False
                 def setup(self, code, osamp, dataOut):
                     self.__profIndex = 0
                     self.code = code
                     self.nCode = len(code)
                     self.nBaud = len(code[0])
                     if (osamp != None) and (osamp >1):
                         self.osamp = osamp
                         self.code = numpy.repeat(code, repeats=self.osamp, axis=1)
                         self.nBaud = self.nBaud*self.osamp
                     self.__nChannels = dataOut.nChannels
                     self.__nProfiles = dataOut.nProfiles
                     self.__nHeis = dataOut.nHeights
                     if self.__nHeis < self.nBaud:
                         raise ValueError('Number of heights (%d) should be greater than number of bauds (%d)' %(self.__nHeis, self.nBaud))
                     #Frequency
                     __codeBuffer = numpy.zeros((self.nCode, self.__nHeis), dtype=numpy.complex)
                     __codeBuffer[:,0:self.nBaud] = self.code
                     self.fft_code = numpy.conj(numpy.fft.fft(__codeBuffer, axis=1))
                     if dataOut.flagDataAsBlock:
                         self.ndatadec = self.__nHeis #- self.nBaud + 1
                         self.datadecTime = numpy.zeros((self.__nChannels, self.__nProfiles, self.ndatadec), dtype=numpy.complex)
                     else:
                         #Time
                         self.ndatadec = self.__nHeis #- self.nBaud + 1
                         self.datadecTime = numpy.zeros((self.__nChannels, self.ndatadec), dtype=numpy.complex)
                 def __convolutionInFreq(self, data):
                     fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
                     fft_data = numpy.fft.fft(data, axis=1)
                     conv = fft_data*fft_code
                     data = numpy.fft.ifft(conv,axis=1)
                     return data
                 def __convolutionInFreqOpt(self, data):
                     raise NotImplementedError
                 def __convolutionInTime(self, data):
                     code = self.code[self.__profIndex]
                     for i in range(self.__nChannels):
                         self.datadecTime[i,:] = numpy.correlate(data[i,:], code, mode='full')[self.nBaud-1:]
                     return self.datadecTime
                 def __convolutionByBlockInTime(self, data):
                     repetitions = int(self.__nProfiles / self.nCode)
                     junk = numpy.lib.stride_tricks.as_strided(self.code, (repetitions, self.code.size), (0, self.code.itemsize))
                     junk = junk.flatten()
                     code_block = numpy.reshape(junk, (self.nCode*repetitions, self.nBaud))
                     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
                 def __convolutionByBlockInFreq(self, data):
                     raise NotImplementedError("Decoder by frequency fro Blocks not implemented")
                     fft_code = self.fft_code[self.__profIndex].reshape(1,-1)
                     fft_data = numpy.fft.fft(data, axis=2)
                     conv = fft_data*fft_code
                     data = numpy.fft.ifft(conv,axis=2)
                     return data
                 def run(self, dataOut, code=None, nCode=None, nBaud=None, mode = 0, osamp=None, times=None):
                     if dataOut.flagDecodeData:
                         print("This data is already decoded, recoding again ...")
                     if not self.isConfig:
                         if code is None:
                             if dataOut.code is None:
                                 raise ValueError("Code could not be read from %s instance. Enter a value in Code parameter" %dataOut.type)
                             code = dataOut.code
                         else:
                             code = numpy.array(code).reshape(nCode,nBaud)
                         self.setup(code, osamp, dataOut)
                         self.isConfig = True
                         if mode == 3:
                             sys.stderr.write("Decoder Warning: mode=%d is not valid, using mode=0\n" %mode)
                         if times != None:
                             sys.stderr.write("Decoder Warning: Argument 'times' in not used anymore\n")
                     if self.code is None:
                         print("Fail decoding: Code is not defined.")
                         return
                     self.__nProfiles = dataOut.nProfiles
                     datadec = None
                     if mode == 3:
                         mode = 0
                     if dataOut.flagDataAsBlock:
                         """
                         Decoding when data have been read as block,
                         """
                         if mode == 0:
                             datadec = self.__convolutionByBlockInTime(dataOut.data)
                         if mode == 1:
                             datadec = self.__convolutionByBlockInFreq(dataOut.data)
                     else:
                         """
                         Decoding when data have been read profile by profile
                         """
                         if mode == 0:
                             datadec = self.__convolutionInTime(dataOut.data)
                         if mode == 1:
                             datadec = self.__convolutionInFreq(dataOut.data)
                         if mode == 2:
                             datadec = self.__convolutionInFreqOpt(dataOut.data)
                     if datadec is None:
                         raise ValueError("Codification mode selected is not valid: mode=%d. Try selecting 0 or 1" %mode)
                     dataOut.code = self.code
                     dataOut.nCode = self.nCode
                     dataOut.nBaud = self.nBaud
                     dataOut.data = datadec
                     dataOut.heightList = dataOut.heightList[0:datadec.shape[-1]]
                     dataOut.flagDecodeData = True #asumo q la data esta decodificada
                     if self.__profIndex == self.nCode-1:
                         self.__profIndex = 0
                         return dataOut
                     self.__profIndex += 1
                     return dataOut
                 #        dataOut.flagDeflipData = True #asumo q la data no esta sin flip
             class ProfileConcat(Operation):
                 isConfig = False
                 buffer = None
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.profileIndex = 0
                 def reset(self):
                     self.buffer = numpy.zeros_like(self.buffer)
                     self.start_index = 0
                     self.times = 1
                 def setup(self, data, m, n=1):
                     self.buffer = numpy.zeros((data.shape[0],data.shape[1]*m),dtype=type(data[0,0]))
                     self.nHeights = data.shape[1]#.nHeights
                     self.start_index = 0
                     self.times = 1
                 def concat(self, data):
                     self.buffer[:,self.start_index:self.nHeights*self.times] = data.copy()
                     self.start_index = self.start_index + self.nHeights
                 def run(self, dataOut, m):
                     dataOut.flagNoData = True
                     if not self.isConfig:
                         self.setup(dataOut.data, m, 1)
                         self.isConfig = True
                     if dataOut.flagDataAsBlock:
                         raise ValueError("ProfileConcat can only be used when voltage have been read profile by profile, getBlock = False")
                     else:
                         self.concat(dataOut.data)
                         self.times += 1
                         if self.times > m:
                             dataOut.data = self.buffer
                             self.reset()
                             dataOut.flagNoData = False
                             # se deben actualizar mas propiedades del header y del objeto dataOut, por ejemplo, las alturas
                             deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                             xf = dataOut.heightList[0] + dataOut.nHeights * deltaHeight * m
                             dataOut.heightList = numpy.arange(dataOut.heightList[0], xf, deltaHeight)
                             dataOut.ippSeconds *= m
                     return dataOut
             class ProfileSelector(Operation):
                 profileIndex = None
                 # Tamanho total de los perfiles
                 nProfiles = None
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.profileIndex = 0
                 def incProfileIndex(self):
                     self.profileIndex += 1
                     if self.profileIndex >= self.nProfiles:
                         self.profileIndex = 0
                 def isThisProfileInRange(self, profileIndex, minIndex, maxIndex):
                     if profileIndex < minIndex:
                         return False
                     if profileIndex > maxIndex:
                         return False
                     return True
                 def isThisProfileInList(self, profileIndex, profileList):
                     if profileIndex not in profileList:
                         return False
                     return True
                 def run(self, dataOut, profileList=None, profileRangeList=None, beam=None, byblock=False, rangeList = None, nProfiles=None):
                     """
                     ProfileSelector:
                     Inputs:
                         profileList        :    Index of profiles selected. Example: profileList = (0,1,2,7,8)
                         profileRangeList    :    Minimum and maximum profile indexes. Example: profileRangeList = (4, 30)
                         rangeList            :    List of profile ranges. Example: rangeList = ((4, 30), (32, 64), (128, 256))
                     """
                     if rangeList is not None:
                         if type(rangeList[0]) not in (tuple, list):
                             rangeList = [rangeList]
                     dataOut.flagNoData = True
                     if dataOut.flagDataAsBlock:
                         """
                         data dimension  = [nChannels, nProfiles, nHeis]
                         """
                         if profileList != None:
                             dataOut.data = dataOut.data[:,profileList,:]
                         if profileRangeList != None:
                             minIndex = profileRangeList[0]
                             maxIndex = profileRangeList[1]
                             profileList = list(range(minIndex, maxIndex+1))
                             dataOut.data = dataOut.data[:,minIndex:maxIndex+1,:]
                         if rangeList != None:
                             profileList = []
                             for thisRange in rangeList:
                                 minIndex = thisRange[0]
                                 maxIndex = thisRange[1]
                                 profileList.extend(list(range(minIndex, maxIndex+1)))
                             dataOut.data = dataOut.data[:,profileList,:]
                         dataOut.nProfiles = len(profileList)
                         dataOut.profileIndex = dataOut.nProfiles - 1
                         dataOut.flagNoData = False
                         return dataOut
                     """
                     data dimension  = [nChannels, nHeis]
                     """
                     if profileList != None:
                         if self.isThisProfileInList(dataOut.profileIndex, profileList):
                             self.nProfiles = len(profileList)
                             dataOut.nProfiles = self.nProfiles
                             dataOut.profileIndex = self.profileIndex
                             dataOut.flagNoData = False
                             self.incProfileIndex()
                         return dataOut
                     if profileRangeList != None:
                         minIndex = profileRangeList[0]
                         maxIndex = profileRangeList[1]
                         if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
                             self.nProfiles = maxIndex - minIndex + 1
                             dataOut.nProfiles = self.nProfiles
                             dataOut.profileIndex = self.profileIndex
                             dataOut.flagNoData = False
                             self.incProfileIndex()
                         return dataOut
                     if rangeList != None:
                         nProfiles = 0
                         for thisRange in rangeList:
                             minIndex = thisRange[0]
                             maxIndex = thisRange[1]
                             nProfiles += maxIndex - minIndex + 1
                         for thisRange in rangeList:
                             minIndex = thisRange[0]
                             maxIndex = thisRange[1]
                             if self.isThisProfileInRange(dataOut.profileIndex, minIndex, maxIndex):
                                 self.nProfiles = nProfiles
                                 dataOut.nProfiles = self.nProfiles
                                 dataOut.profileIndex = self.profileIndex
                                 dataOut.flagNoData = False
                                 self.incProfileIndex()
                                 break
                         return dataOut
                     if beam != None: #beam is only for AMISR data
                         if self.isThisProfileInList(dataOut.profileIndex, dataOut.beamRangeDict[beam]):
                             dataOut.flagNoData = False
                             dataOut.profileIndex = self.profileIndex
                             self.incProfileIndex()
                         return dataOut
                     raise ValueError("ProfileSelector needs profileList, profileRangeList or rangeList parameter")
             class Reshaper(Operation):
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.__buffer = None
                     self.__nitems = 0
                 def __appendProfile(self, dataOut, nTxs):
                     if self.__buffer is None:
                         shape = (dataOut.nChannels, int(dataOut.nHeights/nTxs) )
                         self.__buffer = numpy.empty(shape, dtype = dataOut.data.dtype)
                     ini = dataOut.nHeights * self.__nitems
                     end = ini + dataOut.nHeights
                     self.__buffer[:, ini:end] = dataOut.data
                     self.__nitems += 1
                     return int(self.__nitems*nTxs)
                 def __getBuffer(self):
                     if self.__nitems == int(1./self.__nTxs):
                         self.__nitems = 0
                         return self.__buffer.copy()
                     return None
                 def __checkInputs(self, dataOut, shape, nTxs):
                     if shape is None and nTxs is None:
                         raise ValueError("Reshaper: shape of factor should be defined")
                     if nTxs:
                         if nTxs < 0:
                             raise ValueError("nTxs should be greater than 0")
                         if nTxs < 1 and dataOut.nProfiles % (1./nTxs) != 0:
                             raise ValueError("nProfiles= %d is not divisibled by (1./nTxs) = %f" %(dataOut.nProfiles, (1./nTxs)))
                         shape = [dataOut.nChannels, dataOut.nProfiles*nTxs, dataOut.nHeights/nTxs]
                         return shape, nTxs
                     if len(shape) != 2 and len(shape) !=  3:
                         raise ValueError("shape dimension should be equal to 2 or 3. shape = (nProfiles, nHeis) or (nChannels, nProfiles, nHeis). Actually shape = (%d, %d, %d)" %(dataOut.nChannels, dataOut.nProfiles, dataOut.nHeights))
                     if len(shape) == 2:
                         shape_tuple = [dataOut.nChannels]
                         shape_tuple.extend(shape)
                     else:
                         shape_tuple = list(shape)
                     nTxs = 1.0*shape_tuple[1]/dataOut.nProfiles
                     return shape_tuple, nTxs
                 def run(self, dataOut, shape=None, nTxs=None):
                     shape_tuple, self.__nTxs = self.__checkInputs(dataOut, shape, nTxs)
                     dataOut.flagNoData = True
                     profileIndex = None
                     if dataOut.flagDataAsBlock:
                         dataOut.data = numpy.reshape(dataOut.data, shape_tuple)
                         dataOut.flagNoData = False
                         profileIndex = int(dataOut.nProfiles*self.__nTxs) - 1
                     else:
                         if self.__nTxs < 1:
                             self.__appendProfile(dataOut, self.__nTxs)
                             new_data = self.__getBuffer()
                             if new_data is not None:
                                 dataOut.data = new_data
                                 dataOut.flagNoData = False
                                 profileIndex = dataOut.profileIndex*nTxs
                         else:
                             raise ValueError("nTxs should be greater than 0 and lower than 1, or use VoltageReader(..., getblock=True)")
                     deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                     dataOut.heightList = numpy.arange(dataOut.nHeights/self.__nTxs) * deltaHeight + dataOut.heightList[0]
                     dataOut.nProfiles = int(dataOut.nProfiles*self.__nTxs)
                     dataOut.profileIndex = profileIndex
                     dataOut.ippSeconds /= self.__nTxs
                     return dataOut
             class SplitProfiles(Operation):
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                 def run(self, dataOut, n):
                     dataOut.flagNoData = True
                     profileIndex = None
                     if dataOut.flagDataAsBlock:
                         #nchannels, nprofiles, nsamples
                         shape = dataOut.data.shape
                         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
                         profileIndex = int(dataOut.nProfiles/n) - 1
                     else:
                         raise ValueError("Could not split the data when is read Profile by Profile. Use VoltageReader(..., getblock=True)")
                     deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                     dataOut.heightList = numpy.arange(dataOut.nHeights/n) * deltaHeight + dataOut.heightList[0]
                     dataOut.nProfiles = int(dataOut.nProfiles*n)
                     dataOut.profileIndex = profileIndex
                     dataOut.ippSeconds /= n
                     return dataOut
             class CombineProfiles(Operation):
                 def __init__(self, **kwargs):
                     Operation.__init__(self, **kwargs)
                     self.__remData = None
                     self.__profileIndex = 0
                 def run(self, dataOut, n):
                     dataOut.flagNoData = True
                     profileIndex = None
                     if dataOut.flagDataAsBlock:
                         #nchannels, nprofiles, nsamples
                         shape = dataOut.data.shape
                         new_shape = shape[0], shape[1]/n, shape[2]*n
                         if shape[1] % n != 0:
                             raise ValueError("Could not split the data, n=%d has to be multiple of %d" %(n, shape[1]))
                         dataOut.data = numpy.reshape(dataOut.data, new_shape)
                         dataOut.flagNoData = False
                         profileIndex = int(dataOut.nProfiles*n) - 1
                     else:
                         #nchannels, nsamples
                         if self.__remData is None:
                             newData = dataOut.data
                         else:
                             newData = numpy.concatenate((self.__remData, dataOut.data), axis=1)
                         self.__profileIndex += 1
                         if self.__profileIndex < n:
                             self.__remData = newData
                             #continue
                             return
                         self.__profileIndex = 0
                         self.__remData = None
                         dataOut.data = newData
                         dataOut.flagNoData = False
                         profileIndex = dataOut.profileIndex/n
                     deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
                     dataOut.heightList = numpy.arange(dataOut.nHeights*n) * deltaHeight + dataOut.heightList[0]
                     dataOut.nProfiles = int(dataOut.nProfiles/n)
                     dataOut.profileIndex = profileIndex
                     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
             #
             # class Synchronize(Operation):
             #
             #     isConfig = False
             #     __profIndex = 0
             #
             #     def __init__(self, **kwargs):
             #
             #         Operation.__init__(self, **kwargs)
             # #         self.isConfig = False
             #         self.__powBuffer = None
             #         self.__startIndex = 0
             #         self.__pulseFound = False
             #
             #     def __findTxPulse(self, dataOut, channel=0, pulse_with = None):
             #
             #         #Read data
             #
             #         powerdB = dataOut.getPower(channel = channel)
             #         noisedB = dataOut.getNoise(channel = channel)[0]
             #
             #         self.__powBuffer.extend(powerdB.flatten())
             #
             #         dataArray = numpy.array(self.__powBuffer)
             #
             #         filteredPower = numpy.correlate(dataArray, dataArray[0:self.__nSamples], "same")
             #
             #         maxValue = numpy.nanmax(filteredPower)
             #
             #         if maxValue < noisedB + 10:
             #             #No se encuentra ningun pulso de transmision
             #             return None
             #
             #         maxValuesIndex = numpy.where(filteredPower > maxValue - 0.1*abs(maxValue))[0]
             #
             #         if len(maxValuesIndex) < 2:
             #             #Solo se encontro un solo pulso de transmision de un baudio, esperando por el siguiente TX
             #             return None
             #
             #         phasedMaxValuesIndex = maxValuesIndex - self.__nSamples
             #
             #         #Seleccionar solo valores con un espaciamiento de nSamples
             #         pulseIndex = numpy.intersect1d(maxValuesIndex, phasedMaxValuesIndex)
             #
             #         if len(pulseIndex) < 2:
             #             #Solo se encontro un pulso de transmision con ancho mayor a 1
             #             return None
             #
             #         spacing = pulseIndex[1:] - pulseIndex[:-1]
             #
             #         #remover senales que se distancien menos de 10 unidades o muestras
             #         #(No deberian existir IPP menor a 10 unidades)
             #
             #         realIndex = numpy.where(spacing > 10 )[0]
             #
             #         if len(realIndex) < 2:
             #             #Solo se encontro un pulso de transmision con ancho mayor a 1
             #             return None
             #
             #         #Eliminar pulsos anchos (deja solo la diferencia entre IPPs)
             #         realPulseIndex = pulseIndex[realIndex]
             #
             #         period = mode(realPulseIndex[1:] - realPulseIndex[:-1])[0][0]
             #
             #         print "IPP = %d samples" %period
             #
             #         self.__newNSamples = dataOut.nHeights #int(period)
             #         self.__startIndex = int(realPulseIndex[0])
             #
             #         return 1
             #
             #
             #     def setup(self, nSamples, nChannels, buffer_size = 4):
             #
             #         self.__powBuffer = collections.deque(numpy.zeros( buffer_size*nSamples,dtype=numpy.float),
             #                                           maxlen = buffer_size*nSamples)
             #
             #         bufferList = []
             #
             #         for i in range(nChannels):
             #             bufferByChannel = collections.deque(numpy.zeros( buffer_size*nSamples, dtype=numpy.complex) +  numpy.NAN,
             #                                           maxlen = buffer_size*nSamples)
             #
             #             bufferList.append(bufferByChannel)
             #
             #         self.__nSamples = nSamples
             #         self.__nChannels = nChannels
             #         self.__bufferList = bufferList
             #
             #     def run(self, dataOut, channel = 0):
             #
             #         if not self.isConfig:
             #             nSamples = dataOut.nHeights
             #             nChannels = dataOut.nChannels
             #             self.setup(nSamples, nChannels)
             #             self.isConfig = True
             #
             #         #Append new data to internal buffer
             #         for thisChannel in range(self.__nChannels):
             #             bufferByChannel = self.__bufferList[thisChannel]
             #             bufferByChannel.extend(dataOut.data[thisChannel])
             #
             #         if self.__pulseFound:
             #             self.__startIndex -= self.__nSamples
             #
             #         #Finding Tx Pulse
             #         if not self.__pulseFound:
             #             indexFound = self.__findTxPulse(dataOut, channel)
             #
             #             if indexFound == None:
             #                 dataOut.flagNoData = True
             #                 return
             #
             #             self.__arrayBuffer = numpy.zeros((self.__nChannels, self.__newNSamples), dtype = numpy.complex)
             #             self.__pulseFound = True
             #             self.__startIndex = indexFound
             #
             #         #If pulse was found ...
             #         for thisChannel in range(self.__nChannels):
             #             bufferByChannel = self.__bufferList[thisChannel]
             #             #print self.__startIndex
             #             x = numpy.array(bufferByChannel)
             #             self.__arrayBuffer[thisChannel] = x[self.__startIndex:self.__startIndex+self.__newNSamples]
             #
             #         deltaHeight = dataOut.heightList[1] - dataOut.heightList[0]
             #         dataOut.heightList = numpy.arange(self.__newNSamples)*deltaHeight
             # #             dataOut.ippSeconds = (self.__newNSamples / deltaHeight)/1e6
             #
             #         dataOut.data = self.__arrayBuffer
             #
             #         self.__startIndex += self.__newNSamples
             #
             #         return

General Comments 0

Write
Preview

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

No TODOs yet

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

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