From 823d4012bd344813ebfe0337bc963f53233874b0 2018-05-10 16:53:43 From: ebocanegra Date: 2018-05-10 16:53:43 Subject: [PATCH] Last Commit foreva! --- diff --git a/schainpy/model/graphics/jroplot_parameters.py b/schainpy/model/graphics/jroplot_parameters.py index 89ee920..efe6882 100644 --- a/schainpy/model/graphics/jroplot_parameters.py +++ b/schainpy/model/graphics/jroplot_parameters.py @@ -5,6 +5,181 @@ import inspect from figure import Figure, isRealtime, isTimeInHourRange from plotting_codes import * +class ParamLine(Figure): + + isConfig = None + + def __init__(self): + + self.isConfig = False + self.WIDTH = 300 + self.HEIGHT = 200 + self.counter_imagwr = 0 + + def getSubplots(self): + + nrow = self.nplots + ncol = 3 + return nrow, ncol + + def setup(self, id, nplots, wintitle, show): + + self.nplots = nplots + + self.createFigure(id=id, + wintitle=wintitle, + show=show) + + nrow,ncol = self.getSubplots() + colspan = 3 + rowspan = 1 + + for i in range(nplots): + self.addAxes(nrow, ncol, i, 0, colspan, rowspan) + + def plot_iq(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax): + yreal = y[channelIndexList,:].real + yimag = y[channelIndexList,:].imag + + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Range (Km)" + ylabel = "Intensity - IQ" + + if not self.isConfig: + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle='', + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = min(numpy.nanmin(yreal),numpy.nanmin(yimag)) + if ymax == None: ymax = max(numpy.nanmax(yreal),numpy.nanmax(yimag)) + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(len(self.axesList)): + title = "Channel %d" %(i) + axes = self.axesList[i] + + axes.pline(x, yreal[i,:], + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title) + + axes.addpline(x, yimag[i,:], idline=1, color="red", linestyle="solid", lw=2) + + def plot_power(self, x, y, id, channelIndexList, thisDatetime, wintitle, show, xmin, xmax, ymin, ymax): + y = y[channelIndexList,:] * numpy.conjugate(y[channelIndexList,:]) + yreal = y.real + + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Range (Km)" + ylabel = "Intensity" + + if not self.isConfig: + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle='', + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = numpy.nanmin(yreal) + if ymax == None: ymax = numpy.nanmax(yreal) + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(len(self.axesList)): + title = "Channel %d" %(i) + axes = self.axesList[i] + ychannel = yreal[i,:] + axes.pline(x, ychannel, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title) + + + def run(self, dataOut, id, wintitle="", channelList=None, + xmin=None, xmax=None, ymin=None, ymax=None, save=False, + figpath='./', figfile=None, show=True, wr_period=1, + ftp=False, server=None, folder=None, username=None, password=None): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + """ + + if channelList == None: + channelIndexList = dataOut.channelIndexList + else: + channelIndexList = [] + for channel in channelList: + if channel not in dataOut.channelList: + raise ValueError, "Channel %d is not in dataOut.channelList" + channelIndexList.append(dataOut.channelList.index(channel)) + + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) + + y = dataOut.RR + + title = wintitle + " Scope: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Range (Km)" + ylabel = "Intensity" + + if not self.isConfig: + nplots = len(channelIndexList) + + self.setup(id=id, + nplots=nplots, + wintitle='', + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = numpy.nanmin(y) + if ymax == None: ymax = numpy.nanmax(y) + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(len(self.axesList)): + title = "Channel %d" %(i) + axes = self.axesList[i] + ychannel = y[i,:] + axes.pline(x, ychannel, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title) + + + self.draw() + + str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + "_" + str(dataOut.profileIndex) + figfile = self.getFilename(name = str_datetime) + + self.save(figpath=figpath, + figfile=figfile, + save=save, + ftp=ftp, + wr_period=wr_period, + thisDatetime=thisDatetime) + + class SpcParamPlot(Figure): @@ -800,8 +975,8 @@ class ParametersPlot(Figure): self.isConfig = False self.__nsubplots = 1 - self.WIDTH = 800 - self.HEIGHT = 250 + self.WIDTH = 300 + self.HEIGHT = 550 self.WIDTHPROF = 120 self.HEIGHTPROF = 0 self.counter_imagwr = 0 @@ -910,7 +1085,7 @@ class ParametersPlot(Figure): # thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[0]) title = wintitle + " Parameters Plot" #: %s" %(thisDatetime.strftime("%d-%b-%Y")) xlabel = "" - ylabel = "Range (Km)" + ylabel = "Range (km)" update_figfile = False @@ -954,24 +1129,81 @@ class ParametersPlot(Figure): self.setWinTitle(title) - for i in range(self.nchan): - index = channelIndexList[i] - title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) - axes = self.axesList[i*self.plotFact] - z1 = z[i,:].reshape((1,-1)) - axes.pcolorbuffer(x, y, z1, - xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, - xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, - ticksize=9, cblabel='', cbsize="1%",colormap=colormap) - - if showSNR: - title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) - axes = self.axesList[i*self.plotFact + 1] - SNRdB1 = SNRdB[i,:].reshape((1,-1)) - axes.pcolorbuffer(x, y, SNRdB1, - xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, - xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, - ticksize=9, cblabel='', cbsize="1%",colormap='jet') +# for i in range(self.nchan): +# index = channelIndexList[i] +# title = "Channel %d: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) +# axes = self.axesList[i*self.plotFact] +# z1 = z[i,:].reshape((1,-1)) +# axes.pcolorbuffer(x, y, z1, +# xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, +# xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, +# ticksize=9, cblabel='', cbsize="1%",colormap=colormap) +# +# if showSNR: +# title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) +# axes = self.axesList[i*self.plotFact + 1] +# SNRdB1 = SNRdB[i,:].reshape((1,-1)) +# axes.pcolorbuffer(x, y, SNRdB1, +# xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, +# xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, +# ticksize=9, cblabel='', cbsize="1%",colormap='jet') + + i=0 + index = channelIndexList[i] + title = "Factor de reflectividad Z [dBZ]" + axes = self.axesList[i*self.plotFact] + z1 = z[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, z1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap=colormap) + + if showSNR: + title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[i*self.plotFact + 1] + SNRdB1 = SNRdB[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, SNRdB1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='jet') + + i=1 + index = channelIndexList[i] + title = "Velocidad vertical Doppler [m/s]" + axes = self.axesList[i*self.plotFact] + z1 = z[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, z1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=-10, zmax=10, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='seismic_r') + + if showSNR: + title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[i*self.plotFact + 1] + SNRdB1 = SNRdB[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, SNRdB1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='jet') + + i=2 + index = channelIndexList[i] + title = "Intensidad de lluvia [mm/h]" + axes = self.axesList[i*self.plotFact] + z1 = z[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, z1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=0, zmax=40, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='ocean_r') + + if showSNR: + title = "Channel %d SNR: %s" %(dataOut.channelList[index], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[i*self.plotFact + 1] + SNRdB1 = SNRdB[i,:].reshape((1,-1)) + axes.pcolorbuffer(x, y, SNRdB1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=SNRmin, zmax=SNRmax, + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel='', cbsize="1%",colormap='jet') self.draw() diff --git a/schainpy/model/io/jroIO_bltr.py b/schainpy/model/io/jroIO_bltr.py index 89e3915..70ea6b7 100644 --- a/schainpy/model/io/jroIO_bltr.py +++ b/schainpy/model/io/jroIO_bltr.py @@ -568,7 +568,7 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa if self.flagNoMoreFiles: self.dataOut.flagNoData = True - print 'NoData se vuelve true' + #print 'NoData se vuelve true' return 0 self.fp=self.path @@ -578,7 +578,7 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa self.dataOut.data_cspc =self.data_cspc self.dataOut.data_output=self.data_output - print 'self.dataOut.data_output', shape(self.dataOut.data_output) + #print 'self.dataOut.data_output', shape(self.dataOut.data_output) #self.removeDC() return self.dataOut.data_spc @@ -596,7 +596,7 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa ''' #The address of the folder is generated the name of the .fdt file that will be read - print "File: ",self.fileSelector+1 + #print "File: ",self.fileSelector+1 if self.fileSelector < len(self.filenameList): @@ -608,7 +608,7 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa self.readBlock() #Block reading else: - print 'readFile FlagNoData becomes true' + #print 'readFile FlagNoData becomes true' self.flagNoMoreFiles=True self.dataOut.flagNoData = True return 0 @@ -634,8 +634,8 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa ''' - if self.BlockCounter < self.nFDTdataRecors-2: - print self.nFDTdataRecors, 'CONDICION!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + if self.BlockCounter < self.nFDTdataRecors-1: + #print self.nFDTdataRecors, 'CONDICION' if self.ReadMode==1: rheader = RecordHeaderBLTR(RecCounter=self.BlockCounter+1) elif self.ReadMode==0: @@ -675,7 +675,7 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa self.dataOut.outputInterval= self.dataOut.ippSeconds * self.dataOut.nCohInt * self.dataOut.nIncohInt * self.nProfiles self.data_output=numpy.ones([3,rheader.nHeights])*numpy.NaN - print 'self.data_output', shape(self.data_output) + #print 'self.data_output', shape(self.data_output) self.dataOut.velocityX=[] self.dataOut.velocityY=[] self.dataOut.velocityV=[] @@ -711,9 +711,17 @@ class BLTRSpectraReader (ProcessingUnit, FileHeaderBLTR, RecordHeaderBLTR, JRODa if self.DualModeIndex==self.ReadMode: self.data_fft = numpy.fromfile( startDATA, [('complex',' 0.0001) : -# -# try: -# popt,pcov = curve_fit(gaus,xSamples,yMean,p0=[1,meanGauss,sigma]) -# -# if numpy.amax(popt)>numpy.amax(yMean)*0.3: -# FitGauss=gaus(xSamples,*popt) -# -# else: -# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) -# print 'Verificador: Dentro', Height -# except RuntimeError: -# -# try: -# for j in range(len(ySamples[1])): -# yMean2=numpy.append(yMean2,numpy.average([ySamples[1,j],ySamples[2,j]])) -# popt,pcov = curve_fit(gaus,xSamples,yMean2,p0=[1,meanGauss,sigma]) -# FitGauss=gaus(xSamples,*popt) -# print 'Verificador: Exepcion1', Height -# except RuntimeError: -# -# try: -# popt,pcov = curve_fit(gaus,xSamples,ySamples[1],p0=[1,meanGauss,sigma]) -# FitGauss=gaus(xSamples,*popt) -# print 'Verificador: Exepcion2', Height -# except RuntimeError: -# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) -# print 'Verificador: Exepcion3', Height -# else: -# FitGauss=numpy.ones(len(xSamples))*numpy.mean(yMean) -# #print 'Verificador: Fuera', Height -# -# -# -# Maximun=numpy.amax(yMean) -# eMinus1=Maximun*numpy.exp(-1) -# -# HWpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-eMinus1))) -# HalfWidth= xFrec[HWpos] -# GCpos=Find(FitGauss, numpy.amax(FitGauss)) -# Vpos=Find(FactNorm, numpy.amax(FactNorm)) -# #Vpos=numpy.sum(FactNorm)/len(FactNorm) -# #Vpos=Find(FactNorm, min(FactNorm, key=lambda value:abs(value- numpy.mean(FactNorm) ))) -# #print 'GCpos',GCpos, numpy.amax(FitGauss), 'HWpos',HWpos -# '''****** Getting Fij ******''' -# -# GaussCenter=xFrec[GCpos] -# if (GaussCenter<0 and HalfWidth>0) or (GaussCenter>0 and HalfWidth<0): -# Fij=abs(GaussCenter)+abs(HalfWidth)+0.0000001 -# else: -# Fij=abs(GaussCenter-HalfWidth)+0.0000001 -# -# '''****** Getting Frecuency range of significant data ******''' -# -# Rangpos=Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.10))) -# -# if Rangpos5 and len(FrecRange) 0.: -# self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, Vzon) #Vmag -# #print 'Vmag',Vmag -# else: -# self.dataOut.velocityX=numpy.append(self.dataOut.velocityX, NaN) -# -# if abs(Vx)<100 and abs(Vx) > 0.: -# self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, Vmer) #Vang -# #print 'Vang',Vang -# else: -# self.dataOut.velocityY=numpy.append(self.dataOut.velocityY, NaN) -# -# if abs(GaussCenter)<2: -# self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, xFrec[Vpos]) -# -# else: -# self.dataOut.velocityV=numpy.append(self.dataOut.velocityV, NaN) -# -# -# # print '********************************************' -# # print 'HalfWidth ', HalfWidth -# # print 'Maximun ', Maximun -# # print 'eMinus1 ', eMinus1 -# # print 'Rangpos ', Rangpos -# # print 'GaussCenter ',GaussCenter -# # print 'E01 ',E01 -# # print 'N01 ',N01 -# # print 'E02 ',E02 -# # print 'N02 ',N02 -# # print 'E12 ',E12 -# # print 'N12 ',N12 -# #print 'self.dataOut.velocityX ', self.dataOut.velocityX -# # print 'Fij ', Fij -# # print 'cC ', cC -# # print 'cF ', cF -# # print 'cG ', cG -# # print 'cA ', cA -# # print 'cB ', cB -# # print 'cH ', cH -# # print 'Vx ', Vx -# # print 'Vy ', Vy -# # print 'Vmag ', Vmag -# # print 'Vang ', Vang*180/numpy.pi -# # print 'PhaseSlope ',PhaseSlope[0] -# # print 'PhaseSlope ',PhaseSlope[1] -# # print 'PhaseSlope ',PhaseSlope[2] -# # print '********************************************' -# #print 'data_output',shape(self.dataOut.velocityX), shape(self.dataOut.velocityY) -# -# #print 'self.dataOut.velocityX', len(self.dataOut.velocityX) -# #print 'self.dataOut.velocityY', len(self.dataOut.velocityY) -# #print 'self.dataOut.velocityV', self.dataOut.velocityV -# -# self.data_output[0]=numpy.array(self.dataOut.velocityX) -# self.data_output[1]=numpy.array(self.dataOut.velocityY) -# self.data_output[2]=numpy.array(self.dataOut.velocityV) -# -# prin= self.data_output[0][~numpy.isnan(self.data_output[0])] -# print ' ' -# print 'VmagAverage',numpy.mean(prin) -# print ' ' -# # plt.figure(5) -# # plt.subplot(211) -# # plt.plot(self.dataOut.velocityX,'yo:') -# # plt.subplot(212) -# # plt.plot(self.dataOut.velocityY,'yo:') -# -# # plt.figure(1) -# # # plt.subplot(121) -# # # plt.plot(xFrec,ySamples[0],'k',label='Ch0') -# # # plt.plot(xFrec,ySamples[1],'g',label='Ch1') -# # # plt.plot(xFrec,ySamples[2],'r',label='Ch2') -# # # plt.plot(xFrec,FitGauss,'yo:',label='fit') -# # # plt.legend() -# # plt.title('DATOS A ALTURA DE 2850 METROS') -# # -# # plt.xlabel('Frecuencia (KHz)') -# # plt.ylabel('Magnitud') -# # # plt.subplot(122) -# # # plt.title('Fit for Time Constant') -# # #plt.plot(xFrec,zline) -# # #plt.plot(xFrec,SmoothSPC,'g') -# # plt.plot(xFrec,FactNorm) -# # plt.axis([-4, 4, 0, 0.15]) -# # # plt.xlabel('SelfSpectra KHz') -# # -# # plt.figure(10) -# # # plt.subplot(121) -# # plt.plot(xFrec,ySamples[0],'b',label='Ch0') -# # plt.plot(xFrec,ySamples[1],'y',label='Ch1') -# # plt.plot(xFrec,ySamples[2],'r',label='Ch2') -# # # plt.plot(xFrec,FitGauss,'yo:',label='fit') -# # plt.legend() -# # plt.title('SELFSPECTRA EN CANALES') -# # -# # plt.xlabel('Frecuencia (KHz)') -# # plt.ylabel('Magnitud') -# # # plt.subplot(122) -# # # plt.title('Fit for Time Constant') -# # #plt.plot(xFrec,zline) -# # #plt.plot(xFrec,SmoothSPC,'g') -# # # plt.plot(xFrec,FactNorm) -# # # plt.axis([-4, 4, 0, 0.15]) -# # # plt.xlabel('SelfSpectra KHz') -# # -# # plt.figure(9) -# # -# # -# # plt.title('DATOS SUAVIZADOS') -# # plt.xlabel('Frecuencia (KHz)') -# # plt.ylabel('Magnitud') -# # plt.plot(xFrec,SmoothSPC,'g') -# # -# # #plt.plot(xFrec,FactNorm) -# # plt.axis([-4, 4, 0, 0.15]) -# # # plt.xlabel('SelfSpectra KHz') -# # # -# # plt.figure(2) -# # # #plt.subplot(121) -# # plt.plot(xFrec,yMean,'r',label='Mean SelfSpectra') -# # plt.plot(xFrec,FitGauss,'yo:',label='Ajuste Gaussiano') -# # # plt.plot(xFrec[Rangpos],FitGauss[Find(FitGauss,min(FitGauss, key=lambda value:abs(value-Maximun*0.1)))],'bo') -# # # #plt.plot(xFrec,phase) -# # # plt.xlabel('Suavizado, promediado KHz') -# # plt.title('SELFSPECTRA PROMEDIADO') -# # # #plt.subplot(122) -# # # #plt.plot(xSamples,zline) -# # plt.xlabel('Frecuencia (KHz)') -# # plt.ylabel('Magnitud') -# # plt.legend() -# # # -# # # plt.figure(3) -# # # plt.subplot(311) -# # # #plt.plot(xFrec,phase[0]) -# # # plt.plot(xFrec,phase[0],'g') -# # # plt.subplot(312) -# # # plt.plot(xFrec,phase[1],'g') -# # # plt.subplot(313) -# # # plt.plot(xFrec,phase[2],'g') -# # # #plt.plot(xFrec,phase[2]) -# # # -# # # plt.figure(4) -# # # -# # # plt.plot(xSamples,coherence[0],'b') -# # # plt.plot(xSamples,coherence[1],'r') -# # # plt.plot(xSamples,coherence[2],'g') -# # plt.show() -# # # -# # # plt.clf() -# # # plt.cla() -# # # plt.close() -# -# print ' ' - - - - self.BlockCounter+=2 - - else: - self.fileSelector+=1 - self.BlockCounter=0 - print "Next File" - - - - - diff --git a/schainpy/model/io/jroIO_param.py b/schainpy/model/io/jroIO_param.py index 862be2c..a661d0f 100644 --- a/schainpy/model/io/jroIO_param.py +++ b/schainpy/model/io/jroIO_param.py @@ -645,7 +645,10 @@ class ParamWriter(Operation): dsDict['variable'] = self.dataList[i] #--------------------- Conditionals ------------------------ #There is no data + + if dataAux is None: + return 0 #Not array, just a number @@ -820,7 +823,7 @@ class ParamWriter(Operation): return False def setNextFile(self): - + ext = self.ext path = self.path setFile = self.setFile @@ -1076,10 +1079,11 @@ class ParamWriter(Operation): return def run(self, dataOut, **kwargs): - + if not(self.isConfig): + flagdata = self.setup(dataOut, **kwargs) - + if not(flagdata): return diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index f0a93af..747f3e2 100755 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -182,7 +182,7 @@ class ParametersProc(ProcessingUnit): def target(tups): obj, args = tups - #print 'TARGETTT', obj, args + return obj.FitGau(args) @@ -238,12 +238,6 @@ class SpectralFilters(Operation): Delta = self.Num_Bin/2 - Breaker1R[0] - #Breaker1W=VelRange[numpy.abs(VelRange-(-LimitW)).argmin()] - #Breaker1W=numpy.where(VelRange == Breaker1W) - - #Breaker2W=VelRange[numpy.abs(VelRange-(LimitW)).argmin()] - #Breaker2W=numpy.where(VelRange == Breaker2W) - '''Reacomodando SPCrange''' @@ -279,16 +273,6 @@ class SpectralFilters(Operation): SPCroll[i]=SPCroll[i]-dataOut.noise[i] SPCroll[ numpy.where( SPCroll<0 ) ] = 1e-20 - #self.spc[i, 0:int(Breaker1W[0]) ,:] = dataOut.noise[i] - #self.spc[i, int(Breaker2W[0]):self.Num_Bin ,:] = dataOut.noise[i] - - #self.cspc[i, 0:int(Breaker1W[0]) ,:] = dataOut.noise[i] - #self.cspc[i, int(Breaker2W[0]):self.Num_Bin ,:] = dataOut.noise[i] - - - - - SPC_ch1 = SPCroll SPC_ch2 = SPCcut @@ -296,9 +280,6 @@ class SpectralFilters(Operation): SPCparam = (SPC_ch1, SPC_ch2, self.spc) dataOut.SPCparam = numpy.asarray(SPCparam) - #dataOut.data_pre= (self.spc , self.cspc) - - #dataOut.data_preParam = (self.spc , self.cspc) dataOut.spcparam_range=numpy.zeros([self.Num_Chn,self.Num_Bin+1]) @@ -306,8 +287,6 @@ class SpectralFilters(Operation): dataOut.spcparam_range[1]=TimeRange dataOut.spcparam_range[0]=FrecRange - - class GaussianFit(Operation): @@ -338,21 +317,7 @@ class GaussianFit(Operation): self.spc = dataOut.data_pre[0].copy() - print 'SelfSpectra Shape', numpy.asarray(self.spc).shape - - - #plt.figure(50) - #plt.subplot(121) - #plt.plot(self.spc,'k',label='spc(66)') - #plt.plot(xFrec,ySamples[1],'g',label='Ch1') - #plt.plot(xFrec,ySamples[2],'r',label='Ch2') - #plt.plot(xFrec,FitGauss,'yo:',label='fit') - #plt.legend() - #plt.title('DATOS A ALTURA DE 7500 METROS') - #plt.show() - self.Num_Hei = self.spc.shape[2] - #self.Num_Bin = len(self.spc) self.Num_Bin = self.spc.shape[1] self.Num_Chn = self.spc.shape[0] Vrange = dataOut.abscissaList @@ -376,15 +341,8 @@ class GaussianFit(Operation): gauSPC = pool.map(target, attrs) dataOut.SPCparam = numpy.asarray(SPCparam) - - - print '========================================================' - print 'total_time: ', time.time()-start_time - # re-normalizing spc and noise # This part differs from gg1 - - ''' Parameters: 1. Amplitude @@ -398,10 +356,7 @@ class GaussianFit(Operation): def FitGau(self, X): Vrange, ch, pnoise, noise_, num_intg, SNRlimit = X - #print 'VARSSSS', ch, pnoise, noise, num_intg - - #print 'HEIGHTS', self.Num_Hei - + SPCparam = [] SPC_ch1 = numpy.empty([self.Num_Bin,self.Num_Hei]) SPC_ch2 = numpy.empty([self.Num_Bin,self.Num_Hei]) @@ -411,10 +366,6 @@ class GaussianFit(Operation): for ht in range(self.Num_Hei): - #print (numpy.asarray(self.spc).shape) - - #print 'TTTTT', ch , ht - #print self.spc.shape spc = numpy.asarray(self.spc)[ch,:,ht] @@ -436,16 +387,13 @@ class GaussianFit(Operation): noisebl=wnoise*0.9; noisebh=wnoise*1.1 spc=spc-wnoise - # print 'wnoise', noise_[0], spc_norm_max, wnoise + minx=numpy.argmin(spc) #spcs=spc.copy() spcs=numpy.roll(spc,-minx) cum=numpy.cumsum(spcs) tot_noise=wnoise * self.Num_Bin #64; - #print 'spc' , spcs[5:8] , 'tot_noise', tot_noise - #tot_signal=sum(cum[-5:])/5.; ''' How does this line work? ''' - #snr=tot_signal/tot_noise - #snr=cum[-1]/tot_noise + snr = sum(spcs)/tot_noise snrdB=10.*numpy.log10(snr) @@ -455,8 +403,6 @@ class GaussianFit(Operation): SPC_ch1[:,ht] = 0#numpy.NaN SPCparam = (SPC_ch1,SPC_ch2) continue - #print 'snr',snrdB #, sum(spcs) , tot_noise - #if snrdB<-18 or numpy.isnan(snrdB) or num_intg<4: @@ -602,57 +548,10 @@ class GaussianFit(Operation): shift1,width1,Amplitude1,p1 = [0,0,0,0]#4*[numpy.NaN] -# if choice==0: # pick the single gaussian fit -# Amplitude0=lsq1[0][2] -# shift0=lsq1[0][0] -# width0=lsq1[0][1] -# p0=lsq1[0][3] -# Amplitude1=0. -# shift1=0. -# width1=0. -# p1=0. -# noise=lsq1[0][4] -# elif choice==1: # take the first one of the 2 gaussians fitted -# Amplitude0 = lsq2[0][2] -# shift0 = lsq2[0][0] -# width0 = lsq2[0][1] -# p0 = lsq2[0][3] -# Amplitude1 = lsq2[0][6] # This is 0 in gg1 -# shift1 = lsq2[0][4] # This is 0 in gg1 -# width1 = lsq2[0][5] # This is 0 in gg1 -# p1 = lsq2[0][7] # This is 0 in gg1 -# noise = lsq2[0][8] -# else: # the second one -# Amplitude0 = lsq2[0][6] -# shift0 = lsq2[0][4] -# width0 = lsq2[0][5] -# p0 = lsq2[0][7] -# Amplitude1 = lsq2[0][2] # This is 0 in gg1 -# shift1 = lsq2[0][0] # This is 0 in gg1 -# width1 = lsq2[0][1] # This is 0 in gg1 -# p1 = lsq2[0][3] # This is 0 in gg1 -# noise = lsq2[0][8] - - #print len(noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0) SPC_ch1[:,ht] = noise + Amplitude0*numpy.exp(-0.5*(abs(x-shift0))/width0)**p0 SPC_ch2[:,ht] = noise + Amplitude1*numpy.exp(-0.5*(abs(x-shift1))/width1)**p1 - #print 'SPC_ch1.shape',SPC_ch1.shape - #print 'SPC_ch2.shape',SPC_ch2.shape - #dataOut.data_param = SPC_ch1 SPCparam = (SPC_ch1,SPC_ch2) - #GauSPC[1] = SPC_ch2 - -# print 'shift0', shift0 -# print 'Amplitude0', Amplitude0 -# print 'width0', width0 -# print 'p0', p0 -# print '========================' -# print 'shift1', shift1 -# print 'Amplitude1', Amplitude1 -# print 'width1', width1 -# print 'p1', p1 -# print 'noise', noise -# print 's_noise', wnoise + return GauSPC @@ -686,7 +585,7 @@ class GaussianFit(Operation): def misfit2(self,state,y_data,x,num_intg): return num_intg*sum((numpy.log(y_data)-numpy.log(self.y_model2(x,state)))**2)#/(64-9.) - + class PrecipitationProc(Operation): @@ -703,6 +602,12 @@ class PrecipitationProc(Operation): Parameters affected: ''' + + def __init__(self, **kwargs): + Operation.__init__(self, **kwargs) + self.i=0 + + def gaus(self,xSamples,Amp,Mu,Sigma): return ( Amp / ((2*numpy.pi)**0.5 * Sigma) ) * numpy.exp( -( xSamples - Mu )**2 / ( 2 * (Sigma**2) )) @@ -739,11 +644,16 @@ class PrecipitationProc(Operation): else: self.spc = dataOut.SPCparam[1].copy() #dataOut.data_pre[0].copy() # + + """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX""" + + self.spc[:,:,0:7]= numpy.NaN + + """##########################################""" + self.Num_Hei = self.spc.shape[2] self.Num_Bin = self.spc.shape[1] self.Num_Chn = self.spc.shape[0] - print '==================== SPC SHAPE',numpy.shape(self.spc) - ''' Se obtiene la constante del RADAR ''' @@ -758,10 +668,8 @@ class PrecipitationProc(Operation): Numerator = ( (4*numpy.pi)**3 * aL**2 * 16 * numpy.log(2) ) Denominator = ( Pt * Gt * Gr * Lambda**2 * SPEED_OF_LIGHT * tauW * numpy.pi * ThetaT * ThetaR) - RadarConstant = 5e-27 * Numerator / Denominator # - print '***' - print '*** RadarConstant' , RadarConstant, '****' - print '***' + RadarConstant = 5e-26 * Numerator / Denominator # + ''' ============================= ''' self.spc[0] = (self.spc[0]-dataOut.noise[0]) @@ -781,8 +689,6 @@ class PrecipitationProc(Operation): VelMeteoro = numpy.mean(SPCmean,axis=0) - #print '==================== Vel SHAPE',VelMeteoro - D_range = numpy.zeros([self.Num_Bin,self.Num_Hei]) SIGMA = numpy.zeros([self.Num_Bin,self.Num_Hei]) N_dist = numpy.zeros([self.Num_Bin,self.Num_Hei]) @@ -847,7 +753,7 @@ class PrecipitationProc(Operation): dBZe = 10*numpy.log10(Ze) dBZ = 10*numpy.log10(Z) - dataOut.data_output = Z + dataOut.data_output = RR[8] dataOut.data_param = numpy.ones([3,self.Num_Hei]) dataOut.channelList = [0,1,2] @@ -855,31 +761,6 @@ class PrecipitationProc(Operation): dataOut.data_param[1]=V_mean dataOut.data_param[2]=RR - #print 'VELRANGE', Velrange - #print 'Range', len(Range) - #print 'delv',del_V -# print 'DRANGE', D_range[:,56] -# #print 'NOISE', dataOut.noise[0], 10*numpy.log10(dataOut.noise[0]) -# print 'radarconstant', RadarConstant -# #print 'ETAn SHAPE', ETAn.shape -# # print 'ETAn ', numpy.nansum(ETAn, axis=0) -# # print 'ETAd ', numpy.nansum(ETAd, axis=0) -# print 'Pr ', numpy.nansum(Pr, axis=0) -# print 'dataOut.SPCparam[1]', numpy.nansum(dataOut.SPCparam[1][0], axis=0) -# print 'Ze ', dBZe - print 'Z ', dBZ -# print 'Ndist',N_dist[:,56] - print 'RR2 ', RR2 - print 'RR ', RR - print 'Vr', V_mean - #print 'RR2 ', dBRR2 - #print 'D_mean', D_mean - #print 'del_V', del_V - #print 'D_range',D_range.shape, D_range[:,30] - #print 'Velrange', Velrange - #print 'numpy.nansum( N_dist[:,R]', numpy.nansum( N_dist, axis=0) - #print 'dataOut.data_param SHAPE', dataOut.data_param.shape - def dBZeMODE2(self, dataOut): # Processing for MIRA35C @@ -894,7 +775,7 @@ class PrecipitationProc(Operation): data_output = numpy.ones([self.Num_Chn , self.Num_Hei])*numpy.NaN ETA = numpy.sum(SNR,1) - print 'ETA' , ETA + ETA = numpy.where(ETA is not 0. , ETA, numpy.NaN) Ze = numpy.ones([self.Num_Chn, self.Num_Hei] ) @@ -956,6 +837,14 @@ class FullSpectralAnalysis(Operation): spc = dataOut.data_pre[0].copy() cspc = dataOut.data_pre[1] + """NOTA SE DEBE REMOVER EL RANGO DEL PULSO TX""" + + SNRspc = spc.copy() + SNRspc[:,:,0:7]= numpy.NaN + + """##########################################""" + + nChannel = spc.shape[0] nProfiles = spc.shape[1] nHeights = spc.shape[2] @@ -979,15 +868,10 @@ class FullSpectralAnalysis(Operation): data = dataOut.data_pre noise = dataOut.noise - dataOut.data_SNR = (numpy.mean(spc,axis=1)- noise[0]) / noise[0] + dataOut.data_SNR = (numpy.mean(SNRspc,axis=1)- noise[0]) / noise[0] dataOut.data_SNR[numpy.where( dataOut.data_SNR <0 )] = 1e-20 - - #FirstMoment = dataOut.moments[0,1,:]#numpy.average(dataOut.data_param[:,1,:],0) - #SecondMoment = numpy.average(dataOut.moments[:,2,:],0) - - #SNRdBMean = [] data_output=numpy.ones([spc.shape[0],spc.shape[2]])*numpy.NaN @@ -1000,7 +884,7 @@ class FullSpectralAnalysis(Operation): dbSNR = numpy.average(dbSNR,0) for Height in range(nHeights): - + [Vzon,Vmer,Vver, GaussCenter, PhaseSlope, FitGaussCSPC]= self.WindEstimation(spc, cspc, pairsList, ChanDist, Height, noise, dataOut.spc_range.copy(), dbSNR[Height], SNRlimit) PhaseLine = numpy.append(PhaseLine, PhaseSlope) @@ -1008,25 +892,19 @@ class FullSpectralAnalysis(Operation): velocityX=numpy.append(velocityX, Vzon)#Vmag else: - #print 'Vzon',Vzon velocityX=numpy.append(velocityX, numpy.NaN) if abs(Vmer)<100. and abs(Vmer) > 0.: velocityY=numpy.append(velocityY, -Vmer)#Vang else: - #print 'Vmer',Vmer velocityY=numpy.append(velocityY, numpy.NaN) if dbSNR[Height] > SNRlimit: velocityV=numpy.append(velocityV, -Vver)#FirstMoment[Height]) else: velocityV=numpy.append(velocityV, numpy.NaN) - #FirstMoment[Height]= numpy.NaN -# if SNRdBMean[Height] <12: -# FirstMoment[Height] = numpy.NaN -# velocityX[Height] = numpy.NaN -# velocityY[Height] = numpy.NaN + @@ -1034,21 +912,6 @@ class FullSpectralAnalysis(Operation): data_output[1] = numpy.array(velocityY) #self.moving_average(numpy.array(velocityY) , N=1) data_output[2] = velocityV#FirstMoment - print 'data_output', data_output.shape - #print FirstMoment -# print 'velocityX',numpy.shape(data_output[0]) -# print 'velocityX',data_output[0] -# print ' ' -# print 'velocityY',numpy.shape(data_output[1]) -# print 'velocityY',data_output[1] -# print 'velocityV',data_output[2] -# print 'PhaseLine',PhaseLine - #print numpy.array(velocityY) - #print 'SNR' - #print 10*numpy.log10(dataOut.data_SNR) - #print numpy.shape(10*numpy.log10(dataOut.data_SNR)) - print ' ' - xFrec=FrecRange[0:spc.shape[1]] dataOut.data_output=data_output @@ -1144,16 +1007,7 @@ class FullSpectralAnalysis(Operation): self.Moments(numpy.abs(CSPCSamples[1]), xSamples), self.Moments(numpy.abs(CSPCSamples[2]), xSamples)]) - #print '##### SUMA de SPC #####', len(ySamples) - #print numpy.sum(ySamples[0]) - #print '##### SUMA de CSPC #####', len(coherence) - #print numpy.sum(numpy.abs(CSPCNorm)) - #print numpy.sum(coherence[0]) -# print 'len',len(xSamples) -# print 'CSPCmoments', numpy.shape(CSPCmoments) -# print CSPCmoments -# print '#######################' - + popt=[1e-10,0,1e-10] popt01, popt02, popt12 = [1e-10,1e-10,1e-10], [1e-10,1e-10,1e-10] ,[1e-10,1e-10,1e-10] FitGauss01, FitGauss02, FitGauss12 = numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0, numpy.empty(len(xSamples))*0 @@ -1235,24 +1089,6 @@ class FullSpectralAnalysis(Operation): else: Fij = xSamples[PointGauCenter] - xSamples[PointFij] -# print 'CSPCopt' -# print CSPCopt -# print 'popt' -# print popt -# print '#######################################' - #print 'dataOut.data_param', numpy.shape(data_param) - #print 'dataOut.data_param0', data_param[0,0,Height] - #print 'dataOut.data_param1', data_param[0,1,Height] - #print 'dataOut.data_param2', data_param[0,2,Height] - - -# print 'yMoments', yMoments -# print 'Moments', SPCmoments -# print 'Fij2 Moment', Fij -# #print 'Fij', Fij, 'popt[2]/2',popt[2]/2 -# print 'Fijcspc',Fijcspc -# print '#######################################' - '''****** Taking frequency ranges from SPCs ******''' @@ -1275,24 +1111,13 @@ class FullSpectralAnalysis(Operation): VelRange = xVel[ Range[0] : Range[1] ] - #print 'RANGE: ', Range - #print 'FrecRange', numpy.shape(FrecRange)#,FrecRange - #print 'len: ', len(FrecRange) - '''****** Getting SCPC Slope ******''' for i in range(spc.shape[0]): if len(FrecRange)>5 and len(FrecRange) \ No newline at end of file + \ No newline at end of file diff --git a/schainpy/scripts/testProcBLTR.py b/schainpy/scripts/testProcBLTR.py index 34ffffb..0da67aa 100644 --- a/schainpy/scripts/testProcBLTR.py +++ b/schainpy/scripts/testProcBLTR.py @@ -13,7 +13,7 @@ xmax = '24' desc = "ProcBLTR Test" filename = "ProcBLTR.xml" -figpath = '/media/erick/6F60F7113095A154/BLTR' +figpath = '/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/2018' controllerObj = Project() @@ -22,10 +22,11 @@ controllerObj.setup(id='191', name='test01', description=desc) readUnitConfObj = controllerObj.addReadUnit(datatype='BLTRSpectraReader', #path='/media/erick/6F60F7113095A154/BLTR/', + #path='/data/BLTR', path='/home/erick/Documents/Data/BLTR_Data/fdt/', - endDate='2017/10/19', - startTime='13:00:00', - startDate='2016/11/8', + endDate='2016/11/01', + startTime='0:00:00', + startDate='2016/11/01', endTime='23:59:59', @@ -82,12 +83,12 @@ opObj10 = procUnitConfObj1.addOperation(name='removeDC') # opObj11.addParameter(name='showprofile', value='1', format='int') # opObj11.addParameter(name='timerange', value=str(2*60*60), format='int') -opObj11 = procUnitConfObj1.addOperation(name='CrossSpectraPlot', optype='other') -procUnitConfObj1.addParameter(name='pairsList', value='(0,1),(0,2),(1,2)', format='pairsList') -opObj11.addParameter(name='id', value='2005', format='int') -opObj11.addParameter(name='wintitle', value='CrossSpectraPlot_ShortPulse', format='str') -# opObj11.addParameter(name='exp_code', value='13', format='int') -opObj11.addParameter(name='xaxis', value='Velocity', format='str') +# opObj11 = procUnitConfObj1.addOperation(name='CrossSpectraPlot', optype='other') +# procUnitConfObj1.addParameter(name='pairsList', value='(0,1),(0,2),(1,2)', format='pairsList') +# opObj11.addParameter(name='id', value='2005', format='int') +# opObj11.addParameter(name='wintitle', value='CrossSpectraPlot_ShortPulse', format='str') +# # opObj11.addParameter(name='exp_code', value='13', format='int') +# opObj11.addParameter(name='xaxis', value='Velocity', format='str') #opObj11.addParameter(name='xmin', value='-10', format='float') #opObj11.addParameter(name='xmax', value='10', format='float') #opObj11.addParameter(name='ymin', value='225', format='float') @@ -99,14 +100,14 @@ opObj11.addParameter(name='xaxis', value='Velocity', format='str') # procUnitConfObj2.addParameter(name='pairsList', value='(0,1),(0,2),(1,2)', format='pairsList') procUnitConfObj2 = controllerObj.addProcUnit(datatype='Parameters', inputId=procUnitConfObj1.getId()) -opObj11 = procUnitConfObj2.addOperation(name='SpectralMoments', optype='other') +#opObj11 = procUnitConfObj2.addOperation(name='SpectralMoments', optype='other') opObj22 = procUnitConfObj2.addOperation(name='FullSpectralAnalysis', optype='other') -opObj22.addParameter(name='SNRlimit', value='7', format='float') +opObj22.addParameter(name='SNRlimit', value='2', format='float') opObj22 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other') opObj22.addParameter(name='id', value='4', format='int') opObj22.addParameter(name='wintitle', value='Wind Profiler', format='str') -opObj22.addParameter(name='save', value='1', format='bool') +#opObj22.addParameter(name='save', value='1', format='bool') # opObj22.addParameter(name='figpath', value = '/home/erick/Pictures', format='str') opObj22.addParameter(name='zmin', value='-20', format='int') diff --git a/schainpy/scripts/testRawData.py b/schainpy/scripts/testRawData.py index e2286c0..0085cd7 100644 --- a/schainpy/scripts/testRawData.py +++ b/schainpy/scripts/testRawData.py @@ -7,20 +7,20 @@ if __name__ == '__main__': desc = "Segundo Test" filename = "schain.xml" - pathW='/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra' - figpath = '/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra' + pathW='/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra' + figpath = '/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/pdataCLAIRE/Extra' controllerObj = Project() controllerObj.setup(id='191', name='test01', description=desc) readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader', - path='/media/erick/6F60F7113095A154/CLAIRE/CLAIRE_WINDS_2MHZ/DATA/', + path='/data/CLAIRE/CLAIRE_WINDS_2MHZ/DATA', #path='/home/erick/Documents/Data/Claire_Data/raw', - startDate='2018/02/01', - endDate='2018/02/01', - startTime='12:00:00', - endTime='20:00:00', + startDate='2018/02/22', + endDate='2018/02/24', + startTime='00:00:00', + endTime='23:59:00', online=0, walk=1)