diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index a811d99..f962b34 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -634,8 +634,9 @@ class Spectra(JROData): return freqrange def getAcfRange(self, extrapoints=0): - + #print "miay",self.ippFactor deltafreq = 10. / (self.getFmax() / (self.nFFTPoints * self.ippFactor)) + #print deltafreq freqrange = deltafreq * \ (numpy.arange(self.nFFTPoints + extrapoints) - self.nFFTPoints / 2.) - deltafreq / 2 diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index b4a8dfa..83ac624 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -234,10 +234,10 @@ class ACFPlot(Figure): self.isConfig = False self.__nsubplots = 1 - self.PLOT_CODE = POWER_CODE + self.PLOT_CODE = ACF_CODE - self.WIDTH = 700 - self.HEIGHT = 500 + self.WIDTH = 900 + self.HEIGHT = 700 self.counter_imagwr = 0 def getSubplots(self): @@ -266,7 +266,8 @@ class ACFPlot(Figure): for x in range(ncol): self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1) - def run(self, dataOut, id, wintitle="", channelList=None, + def run(self, dataOut, id, wintitle="", channelList=None,channel=None,nSamples=None, + nSampleList= None,resolutionFactor=None, xmin=None, xmax=None, ymin=None, ymax=None, save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1, server=None, @@ -274,9 +275,28 @@ class ACFPlot(Figure): xaxis="frequency"): - if channelList == None: + channel0 = channel + nSamples = nSamples + resFactor = resolutionFactor + + if nSamples == None: + nSamples = 20 + + if resFactor == None: + resFactor = 5 + #else: + # if nSamples not in dataOut.channelList: + # raise ValueError, "Channel %d is not in %s dataOut.channelList"%(channel0, dataOut.channelList) + + if channel0 == None: + channel0 = 0 + else: + if channel0 not in dataOut.channelList: + raise ValueError, "Channel %d is not in %s dataOut.channelList"%(channel0, dataOut.channelList) + + if channelList == None: channelIndexList = dataOut.channelIndexList - channelList = dataOut.channelList + channelList = dataOut.channelList else: channelIndexList = [] for channel in channelList: @@ -284,53 +304,70 @@ class ACFPlot(Figure): raise ValueError, "Channel %d is not in dataOut.channelList" channelIndexList.append(dataOut.channelList.index(channel)) - factor = dataOut.normFactor - - y = dataOut.getHeiRange() - #z = dataOut.data_spc/factor - deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] - print deltaHeight - z = dataOut.data_spc + factor = dataOut.normFactor + y = dataOut.getHeiRange() + deltaHeight = dataOut.heightList[1]-dataOut.heightList[0] + z = dataOut.data_acf + #z = numpy.where(numpy.isfinite(z), z, numpy.NAN) + shape = dataOut.data_acf.shape + hei_index = numpy.arange(shape[2]) + hei_plot = numpy.arange(nSamples)*resFactor + #print hei_plot + #import matplotlib.pyplot as plt + #c=z[0,:,0]*15+15 + #plt.plot(c) + #plt.show() + #print "HOLA# + + if nSampleList is not None: + for nsample in nSampleList: + if nsample not in dataOut.heightList/deltaHeight: + raise ValueError, "nsample %d is not in %s dataOut.heightList"%(nsample,dataOut.heightList) + + if nSampleList is not None: + hei_plot = numpy.array(nSampleList)*resFactor + + if hei_plot[-1] >= hei_index[-1]: + print ("La cantidad de puntos en altura es %d y la resolucion es %d Km"%(hei_plot.shape[0],deltaHeight*resFactor )) + raise ValueError, "resFactor %d multiplicado por el valor de %d nSamples es mayor a %d cantidad total de puntos"%(resFactor,nSamples,hei_index[-1]) + + #escalamiento -1 a 1 a resolucion (factor de resolucion en altura)* deltaHeight + min = numpy.min(z[0,:,0]) + max =numpy.max(z[0,:,0]) - #z = numpy.where(numpy.isfinite(z), z, numpy.NAN) - shape = dataOut.data_spc.shape for i in range(shape[0]): for j in range(shape[2]): - z[i,:,j]= (z[i,:,j]+1.0)*deltaHeight*5/2.0 + j*deltaHeight + z[i,:,j]= (((z[i,:,j]-min)/(max-min))*deltaHeight*resFactor + j*deltaHeight) #z[i,:,j]= (z[i,:,j]+1.0)*deltaHeight*dataOut.step/2.0 + j*deltaHeight*dataOut.step - hei_index = numpy.arange(shape[2]) - #print hei_index.shape - #b = [] - #for i in range(hei_index.shape[0]): - # if hei_index[i]%30 == 0: - # b.append(hei_index[i]) - - #hei_index= numpy.array(b) - hei_index = hei_index[300:320] - #hei_index = numpy.arange(20)*30+80 - hei_index= numpy.arange(20)*5 + #print deltaHeight + #print resFactor + #print numpy.max(z[0,:,0]) + #import matplotlib.pyplot as plt + #plt.plot((z[0,:,0])*deltaHeight) + #plt.show() + if xaxis == "frequency": x = dataOut.getFreqRange()/1000. - zdB = 10*numpy.log10(z[0,:,hei_index]) + zdB = 10*numpy.log10(z[channel0,:,hei_plot]) xlabel = "Frequency (kHz)" ylabel = "Power (dB)" elif xaxis == "time": x = dataOut.getAcfRange() - zdB = z[0,:,hei_index] + zdB = z[channel0,:,hei_plot] xlabel = "Time (ms)" ylabel = "ACF" else: x = dataOut.getVelRange() - zdB = 10*numpy.log10(z[0,:,hei_index]) + zdB = 10*numpy.log10(z[channel0,:,hei_plot]) xlabel = "Velocity (m/s)" ylabel = "Power (dB)" thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) - title = wintitle + " ACF Plot %s" %(thisDatetime.strftime("%d-%b-%Y")) + title = wintitle + " ACF Plot Ch %s %s" %(channel0,thisDatetime.strftime("%d-%b-%Y")) if not self.isConfig: @@ -346,14 +383,18 @@ class ACFPlot(Figure): if ymin == None: ymin = numpy.nanmin(zdB) if ymax == None: ymax = numpy.nanmax(zdB) + print ("El parametro resFactor es %d y la resolucion en altura es %d"%(resFactor,deltaHeight )) + print ("La cantidad de puntos en altura es %d y la nueva resolucion es %d Km"%(hei_plot.shape[0],deltaHeight*resFactor )) + print ("La altura maxima es %d Km"%(hei_plot[-1]*deltaHeight )) + self.isConfig = True self.setWinTitle(title) - title = "Spectra Cuts: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + title = "ACF Plot: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) axes = self.axesList[0] - legendlabels = ["Range = %dKm" %y[i] for i in hei_index] + legendlabels = ["Range = %dKm" %y[i] for i in hei_plot] axes.pmultilineyaxis( x, zdB, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, @@ -363,6 +404,13 @@ class ACFPlot(Figure): self.draw() + if figfile == None: + str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + name = str_datetime + if ((dataOut.azimuth!=None) and (dataOut.zenith!=None)): + name = name + '_az' + '_%2.2f'%(dataOut.azimuth) + '_zn' + '_%2.2f'%(dataOut.zenith) + figfile = self.getFilename(name) + self.save(figpath=figpath, figfile=figfile, save=save, diff --git a/schainpy/model/graphics/plotting_codes.py b/schainpy/model/graphics/plotting_codes.py index 5f6bc8a..dfcc9bd 100644 --- a/schainpy/model/graphics/plotting_codes.py +++ b/schainpy/model/graphics/plotting_codes.py @@ -12,7 +12,7 @@ TOTAL_CODE = 6 #Total Power. DRIFT_CODE = 7 #Drifts graphics. HEIGHT_CODE = 8 #Height profile. PHASE_CODE = 9 #Signal Phase. -AFC_CODE = 10 #Autocorrelation function. +ACF_CODE = 10 #Autocorrelation function. POWER_CODE = 16 NOISE_CODE = 17 diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index ad02f96..1b54d09 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -65,7 +65,7 @@ class SpectraProc(ProcessingUnit): self.dataOut.beam.azimuthList = self.dataIn.beam.azimuthList self.dataOut.beam.zenithList = self.dataIn.beam.zenithList - self.dataOut.step = self.dataIn.step + self.dataOut.step = self.dataIn.step # def __getFft(self): """ @@ -87,10 +87,8 @@ class SpectraProc(ProcessingUnit): # calculo de self-spectra fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,)) - #print "spec dtype 0",fft_volt.dtype spc = fft_volt * numpy.conjugate(fft_volt) spc = spc.real - #print "spec dtype 1",spc.dtype blocksize = 0 blocksize += dc.size @@ -127,10 +125,7 @@ class SpectraProc(ProcessingUnit): if self.dataIn.type == "Spectra": self.dataOut.copy(self.dataIn) - # if not pairsList: - # pairsList = itertools.combinations(self.dataOut.channelList, 2) - # if self.dataOut.data_cspc is not None: - # self.__selectPairs(pairsList) + print "hi",self.dataOut.ippSeconds if shift_fft: #desplaza a la derecha en el eje 2 determinadas posiciones shift = int(self.dataOut.nFFTPoints/2) @@ -164,16 +159,11 @@ class SpectraProc(ProcessingUnit): self.dataIn.heightList.shape[0]), dtype='complex') - #print self.buffer.shape,"spec2" - #print self.dataIn.heightList.shape[0],"spec3" + if self.dataIn.flagDataAsBlock: - # data dimension: [nChannels, nProfiles, nSamples] nVoltProfiles = self.dataIn.data.shape[1] - # nVoltProfiles = self.dataIn.nProfiles - #print nVoltProfiles,"spec1" - #print nProfiles if nVoltProfiles == nProfiles: self.buffer = self.dataIn.data.copy() self.profIndex = nVoltProfiles @@ -196,7 +186,6 @@ class SpectraProc(ProcessingUnit): else: self.buffer[:, self.profIndex, :] = self.dataIn.data.copy() self.profIndex += 1 - #print self.profIndex,"spectra D" if self.firstdatatime == None: self.firstdatatime = self.dataIn.utctime diff --git a/schainpy/model/proc/jroproc_spectra_acf.py b/schainpy/model/proc/jroproc_spectra_acf.py index 91ad0b0..5088a5c 100644 --- a/schainpy/model/proc/jroproc_spectra_acf.py +++ b/schainpy/model/proc/jroproc_spectra_acf.py @@ -99,16 +99,11 @@ class SpectraAFCProc(ProcessingUnit): dc = fft_volt[:,0,:] #calculo de self-spectra -# fft_volt = numpy.fft.fftshift(fft_volt, axes=(1,)) spc = fft_volt * numpy.conjugate(fft_volt) - - data = numpy.fft.ifft(spc, axis=1) data = numpy.fft.fftshift(data, axes=(1,)) - spc = data.real - - + spc = data blocksize = 0 blocksize += dc.size @@ -145,29 +140,31 @@ class SpectraAFCProc(ProcessingUnit): if self.dataIn.type == "Spectra": self.dataOut.copy(self.dataIn) - spc= self.dataOut.data_spc - data = numpy.fft.ifft(spc, axis=1) - data = numpy.fft.fftshift(data, axes=(1,)) - spc = data.real - shape = spc.shape #nchannels, nprofiles, nsamples + #print "hi",self.dataOut.ippSeconds + + spc = self.dataOut.data_spc + data = numpy.fft.ifft(spc, axis=1) + data = numpy.fft.fftshift(data, axes=(1,)) + acf = numpy.abs(data) # Autocorrelacion LLAMAR A ESTE VALOR ACF + #acf = data.imag + shape = acf.shape #nchannels, nprofiles, nsamples - #print spc.shape - for i in range(shape[0]): - for j in range(shape[2]): - spc[i,:,j]= spc[i,:,j] / numpy.max(numpy.abs(spc[i,:,j])) - #spc = spc[0,:,250] / numpy.max(numpy.abs(spc[0,:,250])) - #print spc.shape #import matplotlib.pyplot as plt - #print spc[0:10] - #plt.plot(spc[0,:,350]) + #plt.plot(acf[0,:,0] / numpy.max(numpy.abs(acf[0,:,0]))) #plt.show() + # Normalizando + for i in range(shape[0]): + for j in range(shape[2]): + acf[i,:,j]= acf[i,:,j] / numpy.max(numpy.abs(acf[i,:,j])) - self.dataOut.data_spc = spc + #import matplotlib.pyplot as plt + #plt.plot(acf[0,:,0]) + #plt.show() + self.dataOut.data_acf = acf return True - if code is not None: self.code = numpy.array(code).reshape(nCode,nBaud) else: diff --git a/schainpy/model/proc/jroproc_voltage.py b/schainpy/model/proc/jroproc_voltage.py index 75b061f..146bd5c 100644 --- a/schainpy/model/proc/jroproc_voltage.py +++ b/schainpy/model/proc/jroproc_voltage.py @@ -1248,8 +1248,8 @@ class SSheightProfiles(Operation): deltaHeight = dataOut.heightList[1] - dataOut.heightList[0] ippSeconds = (deltaHeight*1.0e-6)/(0.15) - - + #print "hi",dataOut.ippSeconds + #print ippSeconds dataOut.data = self.sshProfiles dataOut.flagNoData = False dataOut.heightList = numpy.arange(self.buffer.shape[1]) *self.step*deltaHeight + dataOut.heightList[0] @@ -1258,6 +1258,7 @@ class SSheightProfiles(Operation): dataOut.flagDataAsBlock = True dataOut.ippSeconds = ippSeconds dataOut.step = self.step + #print dataOut.ippSeconds # import collections