diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index 6242474..37fa7d7 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -917,10 +917,18 @@ class Correlation(JROData): class Parameters(JROData): + #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 @@ -933,17 +941,25 @@ class Parameters(JROData): SNR = None #Signal to Noise Ratio - pairsList = None #List of Pairs for Cross correlations or Cross spectrum - initUtcTime = None #Initial UTC time paramInterval = None #Time interval to calculate Parameters in seconds - windsInterval = None #Time interval to calculate Winds in seconds + #Fitting + + constants = None + + error = None + + library = None + + #Output signal + + outputInterval = None #Time interval to calculate output signal in seconds + + data_output = None #Out signal - normFactor = None #Normalization Factor - winds = None #Wind estimations def __init__(self): ''' @@ -960,7 +976,7 @@ class Parameters(JROData): datatime = [] datatime.append(self.initUtcTime) - datatime.append(self.initUtcTime + self.windsInterval - 1) + datatime.append(self.initUtcTime + self.outputInterval - 1) datatime = numpy.array(datatime) diff --git a/schainpy/model/graphics/jroplot_parameters.py b/schainpy/model/graphics/jroplot_parameters.py index ada01c7..ce0b319 100644 --- a/schainpy/model/graphics/jroplot_parameters.py +++ b/schainpy/model/graphics/jroplot_parameters.py @@ -455,7 +455,7 @@ class WindProfilerPlot(Figure): # y = dataOut.heightRange y = dataOut.heightRange - z = dataOut.winds.copy() + z = dataOut.data_output.copy() nplots = z.shape[0] #Number of wind dimensions estimated nplotsw = nplots @@ -771,4 +771,408 @@ class ParametersPlot(Figure): if x[1] >= self.axesList[0].xmax: self.counter_imagwr = wr_period self.__isConfig = False + self.figfile = None + + +class SpectralFittingPlot(Figure): + + __isConfig = None + __nsubplots = None + + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'prm' + + + N = None + ippSeconds = None + + def __init__(self): + self.__isConfig = False + self.__nsubplots = 1 + + self.WIDTH = 450 + self.HEIGHT = 250 + self.WIDTHPROF = 0 + self.HEIGHTPROF = 0 + + def getSubplots(self): + + ncol = int(numpy.sqrt(self.nplots)+0.9) + nrow = int(self.nplots*1./ncol + 0.9) + + return nrow, ncol + + def setup(self, id, nplots, wintitle, showprofile=False, show=True): + + showprofile = False + self.__showprofile = showprofile + self.nplots = nplots + + ncolspan = 5 + colspan = 4 + if showprofile: + ncolspan = 5 + colspan = 4 + self.__nsubplots = 2 + + self.createFigure(id = id, + wintitle = wintitle, + widthplot = self.WIDTH + self.WIDTHPROF, + heightplot = self.HEIGHT + self.HEIGHTPROF, + show=show) + + nrow, ncol = self.getSubplots() + + counter = 0 + for y in range(nrow): + for x in range(ncol): + + if counter >= self.nplots: + break + + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan, colspan, 1) + + if showprofile: + self.addAxes(nrow, ncol*ncolspan, y, x*ncolspan+colspan, 1, 1) + + counter += 1 + + def run(self, dataOut, id, cutHeight=None, fit=False, wintitle="", channelList=None, showprofile=True, + xmin=None, xmax=None, ymin=None, ymax=None, + save=False, figpath='./', figfile=None, show=True): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + showProfile : + xmin : None, + xmax : None, + zmin : None, + zmax : None + """ + + if cutHeight==None: + h=270 + heightindex = numpy.abs(cutHeight - dataOut.heightList).argmin() + cutHeight = dataOut.heightList[heightindex] + + factor = dataOut.normFactor + x = dataOut.abscissaRange[:-1] + #y = dataOut.getHeiRange() + + z = dataOut.data_pre[:,:,heightindex]/factor + z = numpy.where(numpy.isfinite(z), z, numpy.NAN) + avg = numpy.average(z, axis=1) + listChannels = z.shape[0] + + #Reconstruct Function + if fit==True: + groupArray = dataOut.groupList + listChannels = groupArray.reshape((groupArray.size)) + listChannels.sort() + spcFitLine = numpy.zeros(z.shape) + constants = dataOut.constants + + nGroups = groupArray.shape[0] + nChannels = groupArray.shape[1] + nProfiles = z.shape[1] + + for f in range(nGroups): + groupChann = groupArray[f,:] + p = dataOut.data_param[f,:,heightindex] +# p = numpy.array([ 89.343967,0.14036615,0.17086219,18.89835291,1.58388365,1.55099167]) + fitLineAux = dataOut.library.modelFunction(p, constants)*nProfiles + fitLineAux = fitLineAux.reshape((nChannels,nProfiles)) + spcFitLine[groupChann,:] = fitLineAux +# spcFitLine = spcFitLine/factor + + z = z[listChannels,:] + spcFitLine = spcFitLine[listChannels,:] + spcFitLinedB = 10*numpy.log10(spcFitLine) + + zdB = 10*numpy.log10(z) + #thisDatetime = dataOut.datatime + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " Doppler Spectra: %s" %(thisDatetime.strftime("%d-%b-%Y %H:%M:%S")) + xlabel = "Velocity (m/s)" + ylabel = "Spectrum" + + if not self.__isConfig: + + nplots = listChannels.size + + self.setup(id=id, + nplots=nplots, + wintitle=wintitle, + showprofile=showprofile, + show=show) + + if xmin == None: xmin = numpy.nanmin(x) + if xmax == None: xmax = numpy.nanmax(x) + if ymin == None: ymin = numpy.nanmin(zdB) + if ymax == None: ymax = numpy.nanmax(zdB)+2 + + self.__isConfig = True + + self.setWinTitle(title) + for i in range(self.nplots): +# title = "Channel %d: %4.2fdB" %(dataOut.channelList[i]+1, noisedB[i]) + title = "Height %4.1f km\nChannel %d:" %(cutHeight, listChannels[i]+1) + axes = self.axesList[i*self.__nsubplots] + if fit == False: + axes.pline(x, zdB[i,:], + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title + ) + if fit == True: + fitline=spcFitLinedB[i,:] + y=numpy.vstack([zdB[i,:],fitline] ) + legendlabels=['Data','Fitting'] + axes.pmultilineyaxis(x, y, + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, + xlabel=xlabel, ylabel=ylabel, title=title, + legendlabels=legendlabels, marker=None, + linestyle='solid', grid='both') + + self.draw() + + if save: + date = thisDatetime.strftime("%Y%m%d_%H%M%S") + if figfile == None: + figfile = self.getFilename(name = date) + + self.saveFigure(figpath, figfile) + + +class EWDriftsPlot(Figure): + + __isConfig = None + __nsubplots = None + + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'drift' + + def __init__(self): + + self.timerange = 2*60*60 + self.isConfig = False + self.__nsubplots = 1 + + self.WIDTH = 800 + self.HEIGHT = 150 + self.WIDTHPROF = 120 + self.HEIGHTPROF = 0 + self.counter_imagwr = 0 + + self.PLOT_CODE = 0 + self.FTP_WEI = None + self.EXP_CODE = None + self.SUB_EXP_CODE = None + self.PLOT_POS = None + self.tmin = None + self.tmax = None + + self.xmin = None + self.xmax = None + + self.figfile = None + + def getSubplots(self): + + ncol = 1 + nrow = self.nplots + + return nrow, ncol + + def setup(self, id, nplots, wintitle, showprofile=True, show=True): + + self.__showprofile = showprofile + self.nplots = nplots + + ncolspan = 1 + colspan = 1 + + self.createFigure(id = id, + wintitle = wintitle, + widthplot = self.WIDTH + self.WIDTHPROF, + heightplot = self.HEIGHT + self.HEIGHTPROF, + show=show) + + nrow, ncol = self.getSubplots() + + counter = 0 + for y in range(nrow): + if counter >= self.nplots: + break + + self.addAxes(nrow, ncol*ncolspan, y, 0, colspan, 1) + counter += 1 + + def run(self, dataOut, id, wintitle="", channelList=None, + xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, + zmaxVertical = None, zminVertical = None, zmaxZonal = None, zminZonal = None, + timerange=None, SNRthresh = -numpy.inf, SNRmin = None, SNRmax = None, SNR_1 = False, + save=False, figpath='', lastone=0,figfile=None, ftp=False, wr_period=1, show=True, + server=None, folder=None, username=None, password=None, + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0): + """ + + Input: + dataOut : + id : + wintitle : + channelList : + showProfile : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + zmin : None, + zmax : 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)) + + if timerange != None: + self.timerange = timerange + + tmin = None + tmax = None + + x = dataOut.getTimeRange1() +# y = dataOut.heightRange + y = dataOut.heightList + + z = dataOut.data_output + nplots = z.shape[0] #Number of wind dimensions estimated + nplotsw = nplots + + #If there is a SNR function defined + if dataOut.SNR != None: + nplots += 1 + SNR = dataOut.SNR + + if SNR_1: + SNR += 1 + + SNRavg = numpy.average(SNR, axis=0) + + SNRdB = 10*numpy.log10(SNR) + SNRavgdB = 10*numpy.log10(SNRavg) + + ind = numpy.where(SNRavg < 10**(SNRthresh/10))[0] + + for i in range(nplotsw): + z[i,ind] = numpy.nan + + + showprofile = False +# thisDatetime = dataOut.datatime + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " EW Drifts" + xlabel = "" + ylabel = "Height (Km)" + + if not self.__isConfig: + + self.setup(id=id, + nplots=nplots, + wintitle=wintitle, + showprofile=showprofile, + show=show) + + self.xmin, self.xmax = self.getTimeLim(x, xmin, xmax, timerange) + + if ymin == None: ymin = numpy.nanmin(y) + if ymax == None: ymax = numpy.nanmax(y) + + if zmaxZonal == None: zmaxZonal = numpy.nanmax(abs(z[0,:])) + if zminZonal == None: zminZonal = -zmaxZonal + if zmaxVertical == None: zmaxVertical = numpy.nanmax(abs(z[1,:])) + if zminVertical == None: zminVertical = -zmaxVertical + + if dataOut.SNR != None: + if SNRmin == None: SNRmin = numpy.nanmin(SNRavgdB) + if SNRmax == None: SNRmax = numpy.nanmax(SNRavgdB) + + self.FTP_WEI = ftp_wei + self.EXP_CODE = exp_code + self.SUB_EXP_CODE = sub_exp_code + self.PLOT_POS = plot_pos + + self.name = thisDatetime.strftime("%Y%m%d_%H%M%S") + self.__isConfig = True + + + self.setWinTitle(title) + + if ((self.xmax - x[1]) < (x[1]-x[0])): + x[1] = self.xmax + + strWind = ['Zonal','Vertical'] + strCb = 'Velocity (m/s)' + zmaxVector = [zmaxZonal, zmaxVertical] + zminVector = [zminZonal, zminVertical] + + for i in range(nplotsw): + + title = "%s Drifts: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[i*self.__nsubplots] + + z1 = z[i,:].reshape((1,-1)) + + axes.pcolorbuffer(x, y, z1, + xmin=self.xmin, xmax=self.xmax, ymin=ymin, ymax=ymax, zmin=zminVector[i], zmax=zmaxVector[i], + xlabel=xlabel, ylabel=ylabel, title=title, rti=True, XAxisAsTime=True, + ticksize=9, cblabel=strCb, cbsize="1%", colormap="RdBu_r") + + if dataOut.SNR != None: + i += 1 + if SNR_1: + title = "Signal Noise Ratio + 1 (SNR+1): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + else: + title = "Signal Noise Ratio (SNR): %s" %(thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[i*self.__nsubplots] + SNRavgdB = SNRavgdB.reshape((1,-1)) + + axes.pcolorbuffer(x, y, SNRavgdB, + 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() + + if self.figfile == None: + str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + self.figfile = self.getFilename(name = str_datetime) + + if figpath != '': + + self.counter_imagwr += 1 + if (self.counter_imagwr>=wr_period): + # store png plot to local folder + self.saveFigure(figpath, self.figfile) + # store png plot to FTP server according to RT-Web format + name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS) + ftp_filename = os.path.join(figpath, name) + self.saveFigure(figpath, ftp_filename) + + self.counter_imagwr = 0 + + if x[1] >= self.axesList[0].xmax: + self.counter_imagwr = wr_period + self.__isConfig = False self.figfile = None \ No newline at end of file diff --git a/schainpy/model/graphics/mpldriver.py b/schainpy/model/graphics/mpldriver.py index 0c73a32..a8b1e3d 100644 --- a/schainpy/model/graphics/mpldriver.py +++ b/schainpy/model/graphics/mpldriver.py @@ -316,7 +316,7 @@ def createPmultilineYAxis(ax, x, y, xmin, xmax, ymin, ymax, xlabel='', ylabel='' matplotlib.pyplot.ioff() # lines = ax.plot(x, y.T, marker=marker,markersize=markersize,linestyle=linestyle) - lines = ax.plot(x, y.T, linestyle='None', marker='.', markersize=markersize) + lines = ax.plot(x, y.T, linestyle=linestyle, marker=marker, markersize=markersize) leg = ax.legend(lines, legendlabels, loc='upper left', bbox_to_anchor=(1.01, 1.00), numpoints=1, handlelength=1.5, \ handletextpad=0.5, borderpad=0.5, labelspacing=0.5, borderaxespad=0.) diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py index 3538b29..6e9f59d 100644 --- a/schainpy/model/proc/jroproc_parameters.py +++ b/schainpy/model/proc/jroproc_parameters.py @@ -7,7 +7,9 @@ from scipy import stats import re import datetime import copy - +import sys +import importlib +import itertools from jroproc_base import ProcessingUnit, Operation from model.data.jrodata import Parameters @@ -60,7 +62,7 @@ class ParametersProc(ProcessingUnit): def run(self, nSeconds = None, nProfiles = None): - self.dataOut.flagNoData = True + if self.firstdatatime == None: self.firstdatatime = self.dataIn.utctime @@ -68,6 +70,7 @@ class ParametersProc(ProcessingUnit): #---------------------- Voltage Data --------------------------- if self.dataIn.type == "Voltage": + self.dataOut.flagNoData = True if nSeconds != None: self.nSeconds = nSeconds self.nProfiles= int(numpy.floor(nSeconds/(self.dataIn.ippSeconds*self.dataIn.nCohInt))) @@ -100,6 +103,7 @@ class ParametersProc(ProcessingUnit): self.dataOut.abscissaRange = self.dataIn.getVelRange(1) self.dataOut.noise = self.dataIn.getNoise() self.dataOut.normFactor = self.dataIn.normFactor + self.dataOut.flagNoData = False #---------------------- Correlation Data --------------------------- @@ -112,14 +116,14 @@ class ParametersProc(ProcessingUnit): self.dataOut.noise = self.dataIn.noise self.dataOut.normFactor = self.dataIn.normFactor self.dataOut.SNR = self.dataIn.SNR - self.dataOut.pairsList = self.dataIn.pairsList + self.dataOut.groupList = self.dataIn.pairsList + self.dataOut.flagNoData = False self.__updateObjFromInput() - self.dataOut.flagNoData = False self.firstdatatime = None self.dataOut.initUtcTime = self.dataIn.ltctime - self.dataOut.windsInterval = self.dataIn.timeInterval + self.dataOut.outputInterval = self.dataIn.timeInterval #------------------- Get Moments ---------------------------------- def GetMoments(self, channelList = None): @@ -238,7 +242,7 @@ class ParametersProc(ProcessingUnit): self.dataOut.noise self.dataOut.normFactor self.dataOut.SNR - self.dataOut.pairsList + self.dataOut.groupList self.dataOut.nChannels Affected: @@ -251,7 +255,7 @@ class ParametersProc(ProcessingUnit): absc = self.dataOut.abscissaRange[:-1] noise = self.dataOut.noise SNR = self.dataOut.SNR - pairsList = self.dataOut.pairsList + pairsList = self.dataOut.groupList nChannels = self.dataOut.nChannels pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels) self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights)) @@ -1109,6 +1113,137 @@ class ParametersProc(ProcessingUnit): return heights, error + def SpectralFitting(self, getSNR = True, path=None, file=None, groupList=None): + + ''' + Function GetMoments() + + Input: + Output: + Variables modified: + ''' + if path != None: + sys.path.append(path) + self.dataOut.library = importlib.import_module(file) + + #To be inserted as a parameter + groupArray = numpy.array(groupList) +# groupArray = numpy.array([[0,1],[2,3]]) + self.dataOut.groupList = groupArray + + nGroups = groupArray.shape[0] + nChannels = self.dataIn.nChannels + nHeights=self.dataIn.heightList.size + + #Parameters Array + self.dataOut.data_param = None + + #Set constants + constants = self.dataOut.library.setConstants(self.dataIn) + self.dataOut.constants = constants + M = self.dataIn.normFactor + N = self.dataIn.nFFTPoints + ippSeconds = self.dataIn.ippSeconds + K = self.dataIn.nIncohInt + pairsArray = numpy.array(self.dataIn.pairsList) + + #List of possible combinations + listComb = itertools.combinations(numpy.arange(groupArray.shape[1]),2) + indCross = numpy.zeros(len(list(listComb)), dtype = 'int') + + if getSNR: + listChannels = groupArray.reshape((groupArray.size)) + listChannels.sort() + noise = self.dataIn.getNoise() + self.dataOut.SNR = self.__getSNR(self.dataIn.data_spc[listChannels,:,:], noise[listChannels]) + + for i in range(nGroups): + coord = groupArray[i,:] + + #Input data array + data = self.dataIn.data_spc[coord,:,:]/(M*N) + data = data.reshape((data.shape[0]*data.shape[1],data.shape[2])) + + #Cross Spectra data array for Covariance Matrixes + ind = 0 + for pairs in listComb: + pairsSel = numpy.array([coord[x],coord[y]]) + indCross[ind] = int(numpy.where(numpy.all(pairsArray == pairsSel, axis = 1))[0][0]) + ind += 1 + dataCross = self.dataIn.data_cspc[indCross,:,:]/(M*N) + dataCross = dataCross**2/K + + for h in range(nHeights): +# print self.dataOut.heightList[h] + + #Input + d = data[:,h] + + #Covariance Matrix + D = numpy.diag(d**2/K) + ind = 0 + for pairs in listComb: + #Coordinates in Covariance Matrix + x = pairs[0] + y = pairs[1] + #Channel Index + S12 = dataCross[ind,:,h] + D12 = numpy.diag(S12) + #Completing Covariance Matrix with Cross Spectras + D[x*N:(x+1)*N,y*N:(y+1)*N] = D12 + D[y*N:(y+1)*N,x*N:(x+1)*N] = D12 + ind += 1 + Dinv=numpy.linalg.inv(D) + L=numpy.linalg.cholesky(Dinv) + LT=L.T + + dp = numpy.dot(LT,d) + + #Initial values + data_spc = self.dataIn.data_spc[coord,:,h] + p0 = self.dataOut.library.initialValuesFunction(data_spc, constants) + + #Least Squares + minp,covp,infodict,mesg,ier = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants),full_output=True) +# minp,covp = optimize.leastsq(self.__residFunction,p0,args=(dp,LT,constants)) + #Chi square error + error0 = numpy.sum(infodict['fvec']**2)/(2*N) +# error0 = 0 + #Error with Jacobian + error1 = self.dataOut.library.errorFunction(minp,constants,LT) + #Save + if self.dataOut.data_param == None: + self.dataOut.data_param = numpy.zeros((nGroups, minp.size, nHeights))*numpy.nan + self.dataOut.error = numpy.zeros((nGroups, error1.size + 1, nHeights))*numpy.nan + + self.dataOut.error[i,:,h] = numpy.hstack((error0,error1)) + self.dataOut.data_param[i,:,h] = minp + return + + + def __residFunction(self, p, dp, LT, constants): + + fm = self.dataOut.library.modelFunction(p, constants) + fmp=numpy.dot(LT,fm) + + return dp-fmp + + def __getSNR(self, z, noise): + + avg = numpy.average(z, axis=1) + SNR = (avg.T-noise)/noise + SNR = SNR.T + return SNR + + def __chisq(p,chindex,hindex): + #similar to Resid but calculates CHI**2 + [LT,d,fm]=setupLTdfm(p,chindex,hindex) + dp=numpy.dot(LT,d) + fmp=numpy.dot(LT,fm) + chisq=numpy.dot((dp-fmp).T,(dp-fmp)) + return chisq + + class WindProfiler(Operation): @@ -1359,12 +1494,12 @@ class WindProfiler(Operation): winds = correctFactor*winds return winds - def __checkTime(self, currentTime, paramInterval, windsInterval): + def __checkTime(self, currentTime, paramInterval, outputInterval): dataTime = currentTime + paramInterval deltaTime = dataTime - self.__initime - if deltaTime >= windsInterval or deltaTime < 0: + if deltaTime >= outputInterval or deltaTime < 0: self.__dataReady = True return @@ -1462,7 +1597,7 @@ class WindProfiler(Operation): theta_y = theta_y[arrayChannel] velRadial0 = param[:,1,:] #Radial velocity - dataOut.winds, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function + dataOut.data_output, dataOut.heightRange, dataOut.SNR = self.techniqueDBS(velRadial0, theta_x, theta_y, azimuth, correctFactor, horizontalOnly, heightRange, SNR) #DBS Function elif technique == 'SA': @@ -1483,12 +1618,12 @@ class WindProfiler(Operation): tau = dataOut.data_param _lambda = dataOut.C/dataOut.frequency - pairsList = dataOut.pairsList + pairsList = dataOut.groupList nChannels = dataOut.nChannels - dataOut.winds = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor) + dataOut.data_output = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor) dataOut.initUtcTime = dataOut.ltctime - dataOut.windsInterval = dataOut.timeInterval + dataOut.outputInterval = dataOut.timeInterval elif technique == 'Meteors': dataOut.flagNoData = True @@ -1511,7 +1646,7 @@ class WindProfiler(Operation): hmax = kwargs['hmax'] else: hmax = 110 - dataOut.windsInterval = nHours*3600 + dataOut.outputInterval = nHours*3600 if self.__isConfig == False: # self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) @@ -1526,14 +1661,89 @@ class WindProfiler(Operation): else: self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) - self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.windsInterval) #Check if the buffer is ready + self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.outputInterval) #Check if the buffer is ready if self.__dataReady: dataOut.initUtcTime = self.__initime - self.__initime = self.__initime + dataOut.windsInterval #to erase time offset + self.__initime = self.__initime + dataOut.outputInterval #to erase time offset - dataOut.winds, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax) + dataOut.data_output, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax) dataOut.flagNoData = False self.__buffer = None - return \ No newline at end of file + return + +class EWDriftsEstimation(Operation): + + + def __init__(self): + Operation.__init__(self) + + def __correctValues(self, heiRang, phi, velRadial, SNR): + listPhi = phi.tolist() + maxid = listPhi.index(max(listPhi)) + minid = listPhi.index(min(listPhi)) + + rango = range(len(phi)) + # rango = numpy.delete(rango,maxid) + + heiRang1 = heiRang*math.cos(phi[maxid]) + heiRangAux = heiRang*math.cos(phi[minid]) + indOut = (heiRang1 < heiRangAux[0]).nonzero() + heiRang1 = numpy.delete(heiRang1,indOut) + + velRadial1 = numpy.zeros([len(phi),len(heiRang1)]) + SNR1 = numpy.zeros([len(phi),len(heiRang1)]) + + for i in rango: + x = heiRang*math.cos(phi[i]) + y1 = velRadial[i,:] + f1 = interpolate.interp1d(x,y1,kind = 'cubic') + + x1 = heiRang1 + y11 = f1(x1) + + y2 = SNR[i,:] + f2 = interpolate.interp1d(x,y2,kind = 'cubic') + y21 = f2(x1) + + velRadial1[i,:] = y11 + SNR1[i,:] = y21 + + return heiRang1, velRadial1, SNR1 + + def run(self, dataOut, zenith, zenithCorrection): + heiRang = dataOut.heightList + velRadial = dataOut.data_param[:,3,:] + SNR = dataOut.SNR + + zenith = numpy.array(zenith) + zenith -= zenithCorrection + zenith *= numpy.pi/180 + + heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, numpy.abs(zenith), velRadial, SNR) + + alp = zenith[0] + bet = zenith[1] + + w_w = velRadial1[0,:] + w_e = velRadial1[1,:] + + w = (w_w*numpy.sin(bet) - w_e*numpy.sin(alp))/(numpy.cos(alp)*numpy.sin(bet) - numpy.cos(bet)*numpy.sin(alp)) + u = (w_w*numpy.cos(bet) - w_e*numpy.cos(alp))/(numpy.sin(alp)*numpy.cos(bet) - numpy.sin(bet)*numpy.cos(alp)) + + winds = numpy.vstack((u,w)) + + dataOut.heightList = heiRang1 + dataOut.data_output = winds + dataOut.SNR = SNR1 + + dataOut.initUtcTime = dataOut.ltctime + dataOut.outputInterval = dataOut.timeInterval + return + + + + + + \ No newline at end of file diff --git a/schainpy/test/EWDrifts_estimation01.py b/schainpy/test/EWDrifts_estimation01.py new file mode 100644 index 0000000..b46e637 --- /dev/null +++ b/schainpy/test/EWDrifts_estimation01.py @@ -0,0 +1,135 @@ +# DIAS 19 Y 20 FEB 2014 +# Comprobacion de Resultados DBS con SA + +import os, sys + +path = os.path.split(os.getcwd())[0] +sys.path.append(path) + +from controller import * + +desc = "DBS Experiment Test" +filename = "DBStest.xml" + +controllerObj = Project() + +controllerObj.setup(id = '191', name='test01', description=desc) + +#Experimentos + +path = '/host/Jicamarca/EW_Drifts/d2012248' +pathFigure = '/home/propietario/workspace/Graficos/drifts' + + +path = "/home/soporte/Data/drifts" +pathFigure = '/home/soporte/workspace/Graficos/drifts/prueba' + +xmin = 11.75 +xmax = 14.75 +#------------------------------------------------------------------------------------------------ +readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader', + path=path, + startDate='2012/01/01', + endDate='2012/12/31', + startTime='00:00:00', + endTime='23:59:59', + online=0, + walk=1) + +opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock') + +#-------------------------------------------------------------------------------------------------- + +procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId()) + +opObj11 = procUnitConfObj0.addOperation(name='ProfileSelector', optype='other') +opObj11.addParameter(name='profileRangeList', value='0,127', format='intlist') + +opObj11 = procUnitConfObj0.addOperation(name='filterByHeights') +opObj11.addParameter(name='window', value='3', format='int') + +opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other') +# opObj11.addParameter(name='code', value='1,-1', format='floatlist') +# opObj11.addParameter(name='nCode', value='2', format='int') +# opObj11.addParameter(name='nBaud', value='1', format='int') + +procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId()) +procUnitConfObj1.addParameter(name='nFFTPoints', value='128', format='int') +procUnitConfObj1.addParameter(name='nProfiles', value='128', format='int') +procUnitConfObj1.addParameter(name='pairsList', value='(0,1),(2,3)', format='pairsList')#,(2,3) + +opObj11 = procUnitConfObj1.addOperation(name='selectHeights') +# # opObj11.addParameter(name='minHei', value='320.0', format='float') +# # opObj11.addParameter(name='maxHei', value='350.0', format='float') +opObj11.addParameter(name='minHei', value='200.0', format='float') +opObj11.addParameter(name='maxHei', value='600.0', format='float') + +opObj11 = procUnitConfObj1.addOperation(name='selectChannels') +opObj11.addParameter(name='channelList', value='0,1,2,3', format='intlist') + +opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other') +opObj11.addParameter(name='timeInterval', value='300.0', format='float') + +opObj13 = procUnitConfObj1.addOperation(name='removeDC') + +# opObj14 = procUnitConfObj1.addOperation(name='SpectraPlot', optype='other') +# opObj14.addParameter(name='id', value='1', format='int') +# # opObj14.addParameter(name='wintitle', value='Con interf', format='str') +# opObj14.addParameter(name='save', value='1', format='bool') +# opObj14.addParameter(name='figpath', value=pathFigure, format='str') +# # opObj14.addParameter(name='zmin', value='5', format='int') +# opObj14.addParameter(name='zmax', value='30', format='int') +# +# opObj12 = procUnitConfObj1.addOperation(name='RTIPlot', optype='other') +# opObj12.addParameter(name='id', value='2', format='int') +# opObj12.addParameter(name='wintitle', value='RTI Plot', format='str') +# opObj12.addParameter(name='save', value='1', format='bool') +# opObj12.addParameter(name='figpath', value = pathFigure, format='str') +# opObj12.addParameter(name='xmin', value=xmin, format='float') +# opObj12.addParameter(name='xmax', value=xmax, format='float') +# # opObj12.addParameter(name='zmin', value='5', format='int') +# opObj12.addParameter(name='zmax', value='30', format='int') + +#-------------------------------------------------------------------------------------------------- + +procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId()) +opObj20 = procUnitConfObj2.addOperation(name='SpectralFitting') +opObj20.addParameter(name='path', value='/home/soporte/workspace/RemoteSystemsTempFiles', format='str') +opObj20.addParameter(name='file', value='modelSpectralFitting', format='str') +opObj20.addParameter(name='groupList', value='(0,1),(2,3)',format='multiList') + +opObj11 = procUnitConfObj2.addOperation(name='SpectralFittingPlot', optype='other') +opObj11.addParameter(name='id', value='3', format='int') +opObj11.addParameter(name='wintitle', value='DopplerPlot', format='str') +opObj11.addParameter(name='cutHeight', value='350', format='int') +opObj11.addParameter(name='fit', value='1', format='int')#1--True/include fit +opObj11.addParameter(name='save', value='1', format='bool') +opObj11.addParameter(name='figpath', value = pathFigure, format='str') + +opObj11 = procUnitConfObj2.addOperation(name='EWDriftsEstimation', optype='other') +opObj11.addParameter(name='zenith', value='-3.80208,3.10658', format='floatlist') +opObj11.addParameter(name='zenithCorrection', value='0.183201', format='float') + +opObj23 = procUnitConfObj2.addOperation(name='EWDriftsPlot', optype='other') +opObj23.addParameter(name='id', value='4', format='int') +opObj23.addParameter(name='wintitle', value='EW Drifts', format='str') +opObj23.addParameter(name='save', value='1', format='bool') +opObj23.addParameter(name='figpath', value = pathFigure, format='str') +opObj23.addParameter(name='zminZonal', value='-150', format='int') +opObj23.addParameter(name='zmaxZonal', value='150', format='int') +opObj23.addParameter(name='zminVertical', value='-30', format='float') +opObj23.addParameter(name='zmaxVertical', value='30', format='float') +opObj23.addParameter(name='SNR_1', value='1', format='bool') +opObj23.addParameter(name='SNRmax', value='5', format='int') +# opObj23.addParameter(name='SNRthresh', value='-50', format='float') +opObj23.addParameter(name='xmin', value=xmin, format='float') +opObj23.addParameter(name='xmax', value=xmax, format='float') +#-------------------------------------------------------------------------------------------------- +print "Escribiendo el archivo XML" +controllerObj.writeXml(filename) +print "Leyendo el archivo XML" +controllerObj.readXml(filename) + +controllerObj.createObjects() +controllerObj.connectObjects() +controllerObj.run() \ No newline at end of file