diff --git a/schainpy/model/data/jrodata.py b/schainpy/model/data/jrodata.py index 1491ba6..6242474 100644 --- a/schainpy/model/data/jrodata.py +++ b/schainpy/model/data/jrodata.py @@ -743,4 +743,225 @@ class Fits: noise = property(getNoise, "I'm the 'nHeights' property.") datatime = property(getDatatime, "I'm the 'datatime' property") ltctime = property(getltctime, "I'm the 'ltctime' property") + +class Correlation(JROData): + + noise = None + + SNR = None + + pairsAutoCorr = None #Pairs of Autocorrelation + + #-------------------------------------------------- + + data_corr = None + + data_volt = None + + lagT = None # each element value is a profileIndex + + lagR = None # each element value is in km + + pairsList = None + + calculateVelocity = None + + nPoints = None + + nAvg = None + + bufferSize = 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.flagTimeBlock = 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 getLagTRange(self, extrapoints=0): + + lagTRange = self.lagT + diff = lagTRange[1] - lagTRange[0] + extra = numpy.arange(1,extrapoints + 1)*diff + lagTRange[-1] + lagTRange = numpy.hstack((lagTRange, extra)) + + return lagTRange + + def getLagRRange(self, extrapoints=0): + + return self.lagR + + def getPairsList(self): + + return self.pairsList + + def getCalculateVelocity(self): + + return self.calculateVelocity + + def getNPoints(self): + + return self.nPoints + + def getNAvg(self): + + return self.nAvg + + def getBufferSize(self): + + return self.bufferSize + + def getPairsAutoCorr(self): + pairsList = self.pairsList + pairsAutoCorr = numpy.zeros(self.nChannels, dtype = 'int')*numpy.nan + + for l in range(len(pairsList)): + firstChannel = pairsList[l][0] + secondChannel = pairsList[l][1] + + #Obteniendo pares de Autocorrelacion + if firstChannel == secondChannel: + pairsAutoCorr[firstChannel] = int(l) + + pairsAutoCorr = pairsAutoCorr.astype(int) + + return pairsAutoCorr + + 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[range(0,1)] = ind_vel[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(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 + +# pairsList = property(getPairsList, "I'm the 'pairsList' property.") +# nPoints = property(getNPoints, "I'm the 'nPoints' property.") + calculateVelocity = property(getCalculateVelocity, "I'm the 'calculateVelocity' property.") + nAvg = property(getNAvg, "I'm the 'nAvg' property.") + bufferSize = property(getBufferSize, "I'm the 'bufferSize' property.") + + +class Parameters(JROData): + + inputUnit = None #Type of data to be processed + + operation = None #Type of operation to parametrize + + data_param = None #Parameters obtained + + data_pre = None #Data Pre Parametrization + heightRange = None #Heights + + abscissaRange = None #Abscissa, can be velocities, lags or time + + noise = None #Noise Potency + + 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 + + normFactor = None #Normalization Factor + + winds = None #Wind estimations + + def __init__(self): + ''' + Constructor + ''' + self.radarControllerHeaderObj = RadarControllerHeader() + + self.systemHeaderObj = SystemHeader() + + self.type = "Parameters" + + def getTimeRange1(self): + + datatime = [] + + datatime.append(self.initUtcTime) + datatime.append(self.initUtcTime + self.windsInterval - 1) + + datatime = numpy.array(datatime) + + return datatime diff --git a/schainpy/model/graphics/figure.py b/schainpy/model/graphics/figure.py index dfc239a..f60e49b 100644 --- a/schainpy/model/graphics/figure.py +++ b/schainpy/model/graphics/figure.py @@ -572,6 +572,23 @@ class Axes: ylabel=ylabel, title=title, colormap=colormap) + + def polar(self, x, y, + title='', xlabel='',ylabel='',**kwargs): + + if self.__firsttime: + self.plot = self.__driver.createPolar(self.ax, x, y, title = title, xlabel = xlabel, ylabel = ylabel) + self.__firsttime = False + self.x_buffer = x + self.y_buffer = y + return + + self.x_buffer = numpy.hstack((self.x_buffer,x)) + self.y_buffer = numpy.hstack((self.y_buffer,y)) + self.__driver.polar(self.plot, self.x_buffer, self.y_buffer, xlabel=xlabel, + ylabel=ylabel, + title=title) + def __fillGaps(self, x_buffer, y_buffer, z_buffer): deltas = x_buffer[1:] - x_buffer[0:-1] diff --git a/schainpy/model/graphics/jroplot.py b/schainpy/model/graphics/jroplot.py index 7a53039..b298167 100644 --- a/schainpy/model/graphics/jroplot.py +++ b/schainpy/model/graphics/jroplot.py @@ -1,3 +1,5 @@ from jroplot_voltage import * from jroplot_spectra import * from jroplot_heispectra import * +from jroplot_correlation import * +from jroplot_parameters import * \ No newline at end of file diff --git a/schainpy/model/graphics/jroplot_correlation.py b/schainpy/model/graphics/jroplot_correlation.py new file mode 100644 index 0000000..08c845e --- /dev/null +++ b/schainpy/model/graphics/jroplot_correlation.py @@ -0,0 +1,196 @@ +import os +import datetime +import numpy +import copy + +from figure import Figure, isRealtime + +class CorrelationPlot(Figure): + + isConfig = None + __nsubplots = None + + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'corr' + + def __init__(self): + + self.isConfig = False + self.__nsubplots = 1 + + self.WIDTH = 280 + self.HEIGHT = 250 + self.WIDTHPROF = 120 + self.HEIGHTPROF = 0 + self.counter_imagwr = 0 + + self.PLOT_CODE = 1 + self.FTP_WEI = None + self.EXP_CODE = None + self.SUB_EXP_CODE = None + self.PLOT_POS = None + + 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 = 1 + colspan = 1 + if showprofile: + ncolspan = 3 + colspan = 2 + 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, wintitle="", channelList=None, showprofile=False, + xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, + save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1, + server=None, folder=None, username=None, password=None, + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + showProfile : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + zmin : None, + zmax : None + """ + + if dataOut.flagNoData: + return None + + if realtime: + if not(isRealtime(utcdatatime = dataOut.utctime)): + print 'Skipping this plot function' + return + + 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)) + + factor = dataOut.normFactor + lenfactor = factor.shape[1] + x = dataOut.getLagTRange(1) + y = dataOut.getHeiRange() + + z = copy.copy(dataOut.data_corr[:,:,0,:]) + for i in range(dataOut.data_corr.shape[0]): + z[i,:,:] = z[i,:,:]/factor[i,:] + zdB = numpy.abs(z) + + avg = numpy.average(z, axis=1) +# avg = numpy.nanmean(z, axis=1) +# noise = dataOut.noise/factor + + #thisDatetime = dataOut.datatime + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " Correlation" + xlabel = "Lag T (s)" + ylabel = "Range (Km)" + + if not self.isConfig: + + nplots = dataOut.data_corr.shape[0] + + 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(y) + if ymax == None: ymax = numpy.nanmax(y) + if zmin == None: zmin = 0 + if zmax == None: zmax = 1 + + self.FTP_WEI = ftp_wei + self.EXP_CODE = exp_code + self.SUB_EXP_CODE = sub_exp_code + self.PLOT_POS = plot_pos + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(self.nplots): + str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S")) + title = "Channel %d and %d: : %s" %(dataOut.pairsList[i][0]+1,dataOut.pairsList[i][1]+1 , str_datetime) + axes = self.axesList[i*self.__nsubplots] + axes.pcolor(x, y, zdB[i,:,:], + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, + xlabel=xlabel, ylabel=ylabel, title=title, + ticksize=9, cblabel='') + +# if self.__showprofile: +# axes = self.axesList[i*self.__nsubplots +1] +# axes.pline(avgdB[i], y, +# xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax, +# xlabel='dB', ylabel='', title='', +# ytick_visible=False, +# grid='x') +# +# noiseline = numpy.repeat(noisedB[i], len(y)) +# axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2) + + self.draw() + + if figfile == None: + str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + 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, 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 diff --git a/schainpy/model/graphics/jroplot_parameters.py b/schainpy/model/graphics/jroplot_parameters.py new file mode 100644 index 0000000..c27ce59 --- /dev/null +++ b/schainpy/model/graphics/jroplot_parameters.py @@ -0,0 +1,579 @@ +import os +import datetime +import numpy + +from figure import Figure, isRealtime + +class MomentsPlot(Figure): + + isConfig = None + __nsubplots = None + + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'prm' + + def __init__(self): + + self.isConfig = False + self.__nsubplots = 1 + + self.WIDTH = 280 + self.HEIGHT = 250 + self.WIDTHPROF = 120 + self.HEIGHTPROF = 0 + self.counter_imagwr = 0 + + self.PLOT_CODE = 1 + self.FTP_WEI = None + self.EXP_CODE = None + self.SUB_EXP_CODE = None + self.PLOT_POS = None + + 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=True, show=True): + + self.__showprofile = showprofile + self.nplots = nplots + + ncolspan = 1 + colspan = 1 + if showprofile: + ncolspan = 3 + colspan = 2 + 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, wintitle="", channelList=None, showprofile=True, + xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, + save=False, figpath='', figfile=None, show=True, ftp=False, wr_period=1, + server=None, folder=None, username=None, password=None, + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + showProfile : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + zmin : None, + zmax : None + """ + + if dataOut.flagNoData: + return None + + if realtime: + if not(isRealtime(utcdatatime = dataOut.utctime)): + print 'Skipping this plot function' + return + + 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)) + + factor = dataOut.normFactor + x = dataOut.abscissaRange + y = dataOut.heightRange + + z = dataOut.data_pre[channelIndexList,:,:]/factor + z = numpy.where(numpy.isfinite(z), z, numpy.NAN) + avg = numpy.average(z, axis=1) + noise = dataOut.noise/factor + + zdB = 10*numpy.log10(z) + avgdB = 10*numpy.log10(avg) + noisedB = 10*numpy.log10(noise) + + #thisDatetime = dataOut.datatime + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " Parameters" + xlabel = "Velocity (m/s)" + ylabel = "Range (Km)" + + if not self.isConfig: + + nplots = len(channelIndexList) + + 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(y) + if ymax == None: ymax = numpy.nanmax(y) + if zmin == None: zmin = numpy.nanmin(avgdB)*0.9 + if zmax == None: zmax = numpy.nanmax(avgdB)*0.9 + + self.FTP_WEI = ftp_wei + self.EXP_CODE = exp_code + self.SUB_EXP_CODE = sub_exp_code + self.PLOT_POS = plot_pos + + self.isConfig = True + + self.setWinTitle(title) + + for i in range(self.nplots): + str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S")) + title = "Channel %d: %4.2fdB: %s" %(dataOut.channelList[i]+1, noisedB[i], str_datetime) + axes = self.axesList[i*self.__nsubplots] + axes.pcolor(x, y, zdB[i,:,:], + xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, zmin=zmin, zmax=zmax, + xlabel=xlabel, ylabel=ylabel, title=title, + ticksize=9, cblabel='') + #Mean Line + mean = dataOut.data_param[i, 1, :] + axes.addpline(mean, y, idline=0, color="black", linestyle="solid", lw=1) + + if self.__showprofile: + axes = self.axesList[i*self.__nsubplots +1] + axes.pline(avgdB[i], y, + xmin=zmin, xmax=zmax, ymin=ymin, ymax=ymax, + xlabel='dB', ylabel='', title='', + ytick_visible=False, + grid='x') + + noiseline = numpy.repeat(noisedB[i], len(y)) + axes.addpline(noiseline, y, idline=1, color="black", linestyle="dashed", lw=2) + + self.draw() + + if figfile == None: + str_datetime = thisDatetime.strftime("%Y%m%d_%H%M%S") + 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, 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 + +class SkyMapPlot(Figure): + + __isConfig = None + __nsubplots = None + + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'prm' + + def __init__(self): + + self.__isConfig = False + self.__nsubplots = 1 + +# self.WIDTH = 280 +# self.HEIGHT = 250 + self.WIDTH = 600 + self.HEIGHT = 600 + self.WIDTHPROF = 120 + self.HEIGHTPROF = 0 + self.counter_imagwr = 0 + + self.PLOT_CODE = 1 + self.FTP_WEI = None + self.EXP_CODE = None + self.SUB_EXP_CODE = None + self.PLOT_POS = None + + 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): + + 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 = 1,1 + counter = 0 + x = 0 + y = 0 + self.addAxes(1, 1, 0, 0, 1, 1, True) + + def run(self, dataOut, id, wintitle="", channelList=None, showprofile=False, + xmin=None, xmax=None, ymin=None, ymax=None, zmin=None, zmax=None, + save=False, figpath='./', figfile=None, show=True, ftp=False, wr_period=1, + server=None, folder=None, username=None, password=None, + ftp_wei=0, exp_code=0, sub_exp_code=0, plot_pos=0, realtime=False): + + """ + + Input: + dataOut : + id : + wintitle : + channelList : + showProfile : + xmin : None, + xmax : None, + ymin : None, + ymax : None, + zmin : None, + zmax : None + """ + + arrayParameters = dataOut.data_param + error = arrayParameters[:,-1] + indValid = numpy.where(error == 0)[0] + finalMeteor = arrayParameters[indValid,:] + finalAzimuth = finalMeteor[:,4] + finalZenith = finalMeteor[:,5] + + x = finalAzimuth*numpy.pi/180 + y = finalZenith + + + #thisDatetime = dataOut.datatime + thisDatetime = datetime.datetime.utcfromtimestamp(dataOut.getTimeRange()[1]) + title = wintitle + " Parameters" + xlabel = "Zonal Zenith Angle (deg) " + ylabel = "Meridional Zenith Angle (deg)" + + if not self.__isConfig: + + nplots = 1 + + self.setup(id=id, + nplots=nplots, + wintitle=wintitle, + showprofile=showprofile, + show=show) + + 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.firstdate = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S")) + self.__isConfig = True + + self.setWinTitle(title) + + i = 0 + str_datetime = '%s %s'%(thisDatetime.strftime("%Y/%m/%d"),thisDatetime.strftime("%H:%M:%S")) + + axes = self.axesList[i*self.__nsubplots] + nevents = axes.x_buffer.shape[0] + x.shape[0] + title = "Meteor Detection Sky Map\n %s - %s \n Number of events: %5.0f\n" %(self.firstdate,str_datetime,nevents) + axes.polar(x, y, + title=title, xlabel=xlabel, ylabel=ylabel, + ticksize=9, cblabel='') + + self.draw() + + if save: + + self.counter_imagwr += 1 + if (self.counter_imagwr==wr_period): + + if figfile == None: + figfile = self.getFilename(name = self.name) + self.saveFigure(figpath, figfile) + + if ftp: + #provisionalmente envia archivos en el formato de la web en tiempo real + name = self.getNameToFtp(thisDatetime, self.FTP_WEI, self.EXP_CODE, self.SUB_EXP_CODE, self.PLOT_CODE, self.PLOT_POS) + path = '%s%03d' %(self.PREFIX, self.id) + ftp_file = os.path.join(path,'ftp','%s.png'%name) + self.saveFigure(figpath, ftp_file) + ftp_filename = os.path.join(figpath,ftp_file) + + + try: + self.sendByFTP(ftp_filename, server, folder, username, password) + except: + self.counter_imagwr = 0 + raise ValueError, 'Error FTP' + + self.counter_imagwr = 0 + + +class WindProfilerPlot(Figure): + + __isConfig = None + __nsubplots = None + + WIDTHPROF = None + HEIGHTPROF = None + PREFIX = 'wind' + + 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, + zmax_ver = None, zmin_ver = None, SNRmin = None, SNRmax = None, + timerange=None, SNRthresh = None, + 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 + + z = dataOut.winds + 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 + SNRavg = numpy.average(SNR, axis=0) + + SNRdB = 10*numpy.log10(SNR) + SNRavgdB = 10*numpy.log10(SNRavg) + + if SNRthresh == None: SNRthresh = -5.0 + 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 + "Wind" + xlabel = "" + ylabel = "Range (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 zmax == None: zmax = numpy.nanmax(abs(z[range(2),:])) + #if numpy.isnan(zmax): zmax = 50 + if zmin == None: zmin = -zmax + + if nplotsw == 3: + if zmax_ver == None: zmax_ver = numpy.nanmax(abs(z[2,:])) + if zmin_ver == None: zmin_ver = -zmax_ver + + 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', 'Meridional', 'Vertical'] + strCb = ['Velocity (m/s)','Velocity (m/s)','Velocity (cm/s)'] + zmaxVector = [zmax, zmax, zmax_ver] + zminVector = [zmin, zmin, zmin_ver] + windFactor = [1,1,100] + + for i in range(nplotsw): + + title = "%s Wind: %s" %(strWind[i], thisDatetime.strftime("%Y/%m/%d %H:%M:%S")) + axes = self.axesList[i*self.__nsubplots] + + z1 = z[i,:].reshape((1,-1))*windFactor[i] + + 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[i], cbsize="1%", colormap="RdBu_r" ) + + if dataOut.SNR != None: + i += 1 + 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 x[1] >= self.axesList[0].xmax: + self.counter_imagwr = wr_period + self.__isConfig = False + + + 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 + + \ No newline at end of file diff --git a/schainpy/model/graphics/jroplot_spectra.py b/schainpy/model/graphics/jroplot_spectra.py index 1abfc73..fd65d65 100644 --- a/schainpy/model/graphics/jroplot_spectra.py +++ b/schainpy/model/graphics/jroplot_spectra.py @@ -121,7 +121,7 @@ class SpectraPlot(Figure): z = dataOut.data_spc[channelIndexList,:,:]/factor z = numpy.where(numpy.isfinite(z), z, numpy.NAN) avg = numpy.average(z, axis=1) - avg = numpy.nanmean(z, axis=1) + #avg = numpy.nanmean(z, axis=1) noise = dataOut.noise/factor zdB = 10*numpy.log10(z) diff --git a/schainpy/model/graphics/mpldriver.py b/schainpy/model/graphics/mpldriver.py index 7ac5267..0c73a32 100644 --- a/schainpy/model/graphics/mpldriver.py +++ b/schainpy/model/graphics/mpldriver.py @@ -52,14 +52,15 @@ def setTitle(fig, title): fig.suptitle(title) -def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan): +def createAxes(fig, nrow, ncol, xpos, ypos, colspan, rowspan, polar=False): matplotlib.pyplot.ioff() matplotlib.pyplot.figure(fig.number) axes = matplotlib.pyplot.subplot2grid((nrow, ncol), (xpos, ypos), colspan=colspan, - rowspan=rowspan) + rowspan=rowspan, + polar=polar) matplotlib.pyplot.ion() return axes @@ -374,6 +375,49 @@ def pmultilineyaxis(iplot, x, y, xlabel='', ylabel='', title=''): for i in range(len(ax.lines)): line = ax.lines[i] line.set_data(x,y[i,:]) + +def createPolar(ax, x, y, + xlabel='', ylabel='', title='', ticksize = 9, + colormap='jet',cblabel='', cbsize="5%", + XAxisAsTime=False): + + matplotlib.pyplot.ioff() + + ax.plot(x,y,'bo', markersize=5) +# ax.set_rmax(90) + ax.set_ylim(0,90) + ax.set_yticks(numpy.arange(0,90,20)) + ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center' ,size='11') +# ax.text(100, 100, 'example', ha='left', va='center', rotation='vertical') + printLabels(ax, xlabel, '', title) + iplot = ax.lines[-1] + + if '0.' in matplotlib.__version__[0:2]: + print "The matplotlib version has to be updated to 1.1 or newer" + return iplot + + if '1.0.' in matplotlib.__version__[0:4]: + print "The matplotlib version has to be updated to 1.1 or newer" + return iplot + +# if grid != None: +# ax.grid(b=True, which='major', axis=grid) + + matplotlib.pyplot.tight_layout() + + matplotlib.pyplot.ion() + + + return iplot + +def polar(iplot, x, y, xlabel='', ylabel='', title=''): + + ax = iplot.get_axes() + +# ax.text(0, -110, ylabel, rotation='vertical', va ='center', ha = 'center',size='11') + printLabels(ax, xlabel, '', title) + + set_linedata(ax, x, y, idline=0) def draw(fig): diff --git a/schainpy/model/proc/jroproc_correlation.py b/schainpy/model/proc/jroproc_correlation.py new file mode 100644 index 0000000..b608adb --- /dev/null +++ b/schainpy/model/proc/jroproc_correlation.py @@ -0,0 +1,246 @@ +import numpy + +from jroproc_base import ProcessingUnit, Operation +from model.data.jrodata import Correlation + +class CorrelationProc(ProcessingUnit): + + def __init__(self): + + ProcessingUnit.__init__(self) + + self.objectDict = {} + self.buffer = None + self.firstdatatime = None + self.profIndex = 0 + self.dataOut = Correlation() + + def __updateObjFromInput(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.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy() + self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy() + self.dataOut.channelList = self.dataIn.channelList + self.dataOut.heightList = self.dataIn.heightList + self.dataOut.dtype = numpy.dtype([('real','= 0 and AC2 >= 0): + + data1 = numpy.abs(self.dataOut.data_corr[AC1,:,indR,:]) + data2 = numpy.abs(self.dataOut.data_corr[AC2,:,indR,:]) + maxim1 = data1.max(axis = 0) + maxim2 = data1.max(axis = 0) + maxim = numpy.sqrt(maxim1*maxim2) + else: + #In case there is no autocorrelation for the pair + data = numpy.abs(self.dataOut.data_corr[l,:,indR,:]) + maxim = numpy.max(data, axis = 0) + + normFactor[l,:] = maxim + + self.dataOut.normFactor = normFactor + + return 1 + + def run(self, lagT=None, lagR=None, pairsList=None, + nPoints=None, nAvg=None, bufferSize=None, + fullT = False, fullR = False, removeDC = False): + + self.dataOut.flagNoData = True + + if self.dataIn.type == "Correlation": + + self.dataOut.copy(self.dataIn) + + return + + if self.dataIn.type == "Voltage": + + if pairsList == None: + pairsList = [numpy.array([0,0])] + + if nPoints == None: + nPoints = 128 + #------------------------------------------------------------ + #Condicionales para calcular Correlaciones en Tiempo y Rango + if fullT: + lagT = numpy.arange(nPoints*2 - 1) - nPoints + 1 + elif lagT == None: + lagT = numpy.array([0]) + else: + lagT = numpy.array(lagT) + + if fullR: + lagR = numpy.arange(self.dataOut.nHeights) + elif lagR == None: + lagR = numpy.array([0]) + #------------------------------------------------------------- + + if nAvg == None: + nAvg = 1 + + if bufferSize == None: + bufferSize = 0 + + deltaH = self.dataIn.heightList[1] - self.dataIn.heightList[0] + self.dataOut.lagR = numpy.round(numpy.array(lagR)/deltaH) + self.dataOut.pairsList = pairsList + self.dataOut.nPoints = nPoints +# channels = numpy.sort(list(set(list(itertools.chain.from_iterable(pairsList))))) + + if self.buffer == None: + + self.buffer = numpy.zeros((self.dataIn.nChannels,self.dataIn.nProfiles,self.dataIn.nHeights),dtype='complex') + + + self.buffer[:,self.profIndex,:] = self.dataIn.data.copy()[:,:] + + self.profIndex += 1 + + if self.firstdatatime == None: + + self.firstdatatime = self.dataIn.utctime + + if self.profIndex == nPoints: + + tmp = self.buffer[:,0:nPoints,:] + self.buffer = None + self.buffer = tmp + + #--------------- Remover DC ------------ + if removeDC: + self.buffer = self.removeDC(self.buffer) + #--------------------------------------------- + self.dataOut.data_volts = self.buffer + self.__updateObjFromInput() + self.dataOut.data_corr = numpy.zeros((len(pairsList), + len(lagT),len(lagR), + self.dataIn.nHeights), + dtype='complex') + + for l in range(len(pairsList)): + + firstChannel = pairsList[l][0] + secondChannel = pairsList[l][1] + + tmp = None + tmp = numpy.zeros((len(lagT),len(lagR),self.dataIn.nHeights),dtype='complex') + + for t in range(len(lagT)): + + for r in range(len(lagR)): + + idxT = lagT[t] + idxR = lagR[r] + + if idxT >= 0: + vStacked = numpy.vstack((self.buffer[secondChannel,idxT:,:], + numpy.zeros((idxT,self.dataIn.nHeights),dtype='complex'))) + else: + vStacked = numpy.vstack((numpy.zeros((-idxT,self.dataIn.nHeights),dtype='complex'), + self.buffer[secondChannel,:(nPoints + idxT),:])) + + if idxR >= 0: + hStacked = numpy.hstack((vStacked[:,idxR:],numpy.zeros((nPoints,idxR),dtype='complex'))) + else: + hStacked = numpy.hstack((numpy.zeros((nPoints,-idxR),dtype='complex'),vStacked[:,(self.dataOut.nHeights + idxR)])) + + + tmp[t,r,:] = numpy.sum((numpy.conjugate(self.buffer[firstChannel,:,:])*hStacked),axis=0) + + + hStacked = None + vStacked = None + + self.dataOut.data_corr[l,:,:,:] = tmp[:,:,:] + + #Se Calcula los factores de Normalizacion + self.dataOut.pairsAutoCorr = self.dataOut.getPairsAutoCorr() + self.dataOut.lagT = lagT*self.dataIn.ippSeconds*self.dataIn.nCohInt + self.dataOut.lagR = lagR + + self.calculateNormFactor() + + self.dataOut.flagNoData = False + self.buffer = None + self.firstdatatime = None + self.profIndex = 0 + + return \ No newline at end of file diff --git a/schainpy/model/proc/jroproc_parameters.py b/schainpy/model/proc/jroproc_parameters.py new file mode 100644 index 0000000..d0ec1b8 --- /dev/null +++ b/schainpy/model/proc/jroproc_parameters.py @@ -0,0 +1,1521 @@ +import numpy +import math +from scipy import optimize +from scipy import interpolate +from scipy import signal +from scipy import stats +import re +import datetime +import copy + + +from jroproc_base import ProcessingUnit, Operation +from model.data.jrodata import Parameters + + +class ParametersProc(ProcessingUnit): + + nSeconds = None + + def __init__(self): + ProcessingUnit.__init__(self) + + self.objectDict = {} + self.buffer = None + self.firstdatatime = None + self.profIndex = 0 + self.dataOut = Parameters() + + def __updateObjFromInput(self): + + self.dataOut.inputUnit = self.dataIn.type + + 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.radarControllerHeaderObj = self.dataIn.radarControllerHeaderObj.copy() + self.dataOut.systemHeaderObj = self.dataIn.systemHeaderObj.copy() + self.dataOut.channelList = self.dataIn.channelList + self.dataOut.heightList = self.dataIn.heightList + self.dataOut.dtype = numpy.dtype([('real',' m): ss1 = m + + valid = numpy.asarray(range(int(m + bb0 - ss1 + 1))) + ss1 + power = ((spec2[valid] - n0)*fwindow[valid]).sum() + fd = ((spec2[valid]- n0)*freq[valid]*fwindow[valid]).sum()/power + w = math.sqrt(((spec2[valid] - n0)*fwindow[valid]*(freq[valid]- fd)**2).sum()/power) + snr = (spec2.mean()-n0)/n0 + + if (snr < 1.e-20) : + snr = 1.e-20 + + vec_power[ind] = power + vec_fd[ind] = fd + vec_w[ind] = w + vec_snr[ind] = snr + + moments = numpy.vstack((vec_snr, vec_power, vec_fd, vec_w)) + return moments + + #------------------- Get Lags ---------------------------------- + + def GetLags(self): + ''' + Function GetMoments() + + Input: + self.dataOut.data_pre + self.dataOut.abscissaRange + self.dataOut.noise + self.dataOut.normFactor + self.dataOut.SNR + self.dataOut.pairsList + self.dataOut.nChannels + + Affected: + self.dataOut.data_param + + ''' + data = self.dataOut.data_pre + normFactor = self.dataOut.normFactor + nHeights = self.dataOut.nHeights + absc = self.dataOut.abscissaRange[:-1] + noise = self.dataOut.noise + SNR = self.dataOut.SNR + pairsList = self.dataOut.pairsList + nChannels = self.dataOut.nChannels + pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels) + self.dataOut.data_param = numpy.zeros((len(pairsCrossCorr)*2 + 1, nHeights)) + + dataNorm = numpy.abs(data) + for l in range(len(pairsList)): + dataNorm[l,:,:] = dataNorm[l,:,:]/normFactor[l,:] + + self.dataOut.data_param[:-1,:] = self.__calculateTaus(dataNorm, pairsCrossCorr, pairsAutoCorr, absc) + self.dataOut.data_param[-1,:] = self.__calculateLag1Phase(data, pairsAutoCorr, absc) + return + + def __getPairsAutoCorr(self, pairsList, nChannels): + + pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan + + for l in range(len(pairsList)): + firstChannel = pairsList[l][0] + secondChannel = pairsList[l][1] + + #Obteniendo pares de Autocorrelacion + if firstChannel == secondChannel: + pairsAutoCorr[firstChannel] = int(l) + + pairsAutoCorr = pairsAutoCorr.astype(int) + + pairsCrossCorr = range(len(pairsList)) + pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr) + + return pairsAutoCorr, pairsCrossCorr + + def __calculateTaus(self, data, pairsCrossCorr, pairsAutoCorr, lagTRange): + + Pt0 = data.shape[1]/2 + #Funcion de Autocorrelacion + dataAutoCorr = stats.nanmean(data[pairsAutoCorr,:,:], axis = 0) + + #Obtencion Indice de TauCross + indCross = data[pairsCrossCorr,:,:].argmax(axis = 1) + #Obtencion Indice de TauAuto + indAuto = numpy.zeros(indCross.shape,dtype = 'int') + CCValue = data[pairsCrossCorr,Pt0,:] + for i in range(pairsCrossCorr.size): + indAuto[i,:] = numpy.abs(dataAutoCorr - CCValue[i,:]).argmin(axis = 0) + + #Obtencion de TauCross y TauAuto + tauCross = lagTRange[indCross] + tauAuto = lagTRange[indAuto] + + Nan1, Nan2 = numpy.where(tauCross == lagTRange[0]) + + tauCross[Nan1,Nan2] = numpy.nan + tauAuto[Nan1,Nan2] = numpy.nan + tau = numpy.vstack((tauCross,tauAuto)) + + return tau + + def __calculateLag1Phase(self, data, pairs, lagTRange): + data1 = stats.nanmean(data[pairs,:,:], axis = 0) + lag1 = numpy.where(lagTRange == 0)[0][0] + 1 + + phase = numpy.angle(data1[lag1,:]) + + return phase + #------------------- Detect Meteors ------------------------------ + + def DetectMeteors(self, hei_ref = None, tauindex = 0, + predefinedPhaseShifts = None, centerReceiverIndex = 2, + cohDetection = False, cohDet_timeStep = 1, cohDet_thresh = 25, + noise_timeStep = 4, noise_multiple = 4, + multDet_timeLimit = 1, multDet_rangeLimit = 3, + phaseThresh = 20, SNRThresh = 8, + hmin = 70, hmax=110, azimuth = 0) : + + ''' + Function DetectMeteors() + Project developed with paper: + HOLDSWORTH ET AL. 2004 + + Input: + self.dataOut.data_pre + + centerReceiverIndex: From the channels, which is the center receiver + + hei_ref: Height reference for the Beacon signal extraction + tauindex: + predefinedPhaseShifts: Predefined phase offset for the voltge signals + + cohDetection: Whether to user Coherent detection or not + cohDet_timeStep: Coherent Detection calculation time step + cohDet_thresh: Coherent Detection phase threshold to correct phases + + noise_timeStep: Noise calculation time step + noise_multiple: Noise multiple to define signal threshold + + multDet_timeLimit: Multiple Detection Removal time limit in seconds + multDet_rangeLimit: Multiple Detection Removal range limit in km + + phaseThresh: Maximum phase difference between receiver to be consider a meteor + SNRThresh: Minimum SNR threshold of the meteor signal to be consider a meteor + + hmin: Minimum Height of the meteor to use it in the further wind estimations + hmax: Maximum Height of the meteor to use it in the further wind estimations + azimuth: Azimuth angle correction + + Affected: + self.dataOut.data_param + + Rejection Criteria (Errors): + 0: No error; analysis OK + 1: SNR < SNR threshold + 2: angle of arrival (AOA) ambiguously determined + 3: AOA estimate not feasible + 4: Large difference in AOAs obtained from different antenna baselines + 5: echo at start or end of time series + 6: echo less than 5 examples long; too short for analysis + 7: echo rise exceeds 0.3s + 8: echo decay time less than twice rise time + 9: large power level before echo + 10: large power level after echo + 11: poor fit to amplitude for estimation of decay time + 12: poor fit to CCF phase variation for estimation of radial drift velocity + 13: height unresolvable echo: not valid height within 70 to 110 km + 14: height ambiguous echo: more then one possible height within 70 to 110 km + 15: radial drift velocity or projected horizontal velocity exceeds 200 m/s + 16: oscilatory echo, indicating event most likely not an underdense echo + + 17: phase difference in meteor Reestimation + + Data Storage: + Meteors for Wind Estimation (8): + Day Hour | Range Height + Azimuth Zenith errorCosDir + VelRad errorVelRad + TypeError + + ''' + #Get Beacon signal + newheis = numpy.where(self.dataOut.heightList>self.dataOut.radarControllerHeaderObj.Taus[tauindex]) + + if hei_ref != None: + newheis = numpy.where(self.dataOut.heightList>hei_ref) + + heiRang = self.dataOut.getHeiRange() + #Pairs List + pairslist = [] + nChannel = self.dataOut.nChannels + for i in range(nChannel): + if i != centerReceiverIndex: + pairslist.append((centerReceiverIndex,i)) + + #****************REMOVING HARDWARE PHASE DIFFERENCES*************** + # see if the user put in pre defined phase shifts + voltsPShift = self.dataOut.data_pre.copy() + + if predefinedPhaseShifts != None: + hardwarePhaseShifts = numpy.array(predefinedPhaseShifts)*numpy.pi/180 + else: + #get hardware phase shifts using beacon signal + hardwarePhaseShifts = self.__getHardwarePhaseDiff(self.dataOut.data_pre, pairslist, newheis, 10) + hardwarePhaseShifts = numpy.insert(hardwarePhaseShifts,centerReceiverIndex,0) + + voltsPShift = numpy.zeros((self.dataOut.data_pre.shape[0],self.dataOut.data_pre.shape[1],self.dataOut.data_pre.shape[2]), dtype = 'complex') + for i in range(self.dataOut.data_pre.shape[0]): + voltsPShift[i,:,:] = self.__shiftPhase(self.dataOut.data_pre[i,:,:], hardwarePhaseShifts[i]) + #******************END OF REMOVING HARDWARE PHASE DIFFERENCES********* + + #Remove DC + voltsDC = numpy.mean(voltsPShift,1) + voltsDC = numpy.mean(voltsDC,1) + for i in range(voltsDC.shape[0]): + voltsPShift[i] = voltsPShift[i] - voltsDC[i] + + #Don't considerate last heights, theyre used to calculate Hardware Phase Shift + voltsPShift = voltsPShift[:,:,:newheis[0][0]] + + #************ FIND POWER OF DATA W/COH OR NON COH DETECTION (3.4) ********** + #Coherent Detection + if cohDetection: + #use coherent detection to get the net power + cohDet_thresh = cohDet_thresh*numpy.pi/180 + voltsPShift = self.__coherentDetection(voltsPShift, cohDet_timeStep, self.dataOut.timeInterval, pairslist, cohDet_thresh) + + #Non-coherent detection! + powerNet = numpy.nansum(numpy.abs(voltsPShift[:,:,:])**2,0) + #********** END OF COH/NON-COH POWER CALCULATION********************** + + #********** FIND THE NOISE LEVEL AND POSSIBLE METEORS **************** + #Get noise + noise, noise1 = self.__getNoise(powerNet, noise_timeStep, self.dataOut.timeInterval) +# noise = self.getNoise1(powerNet, noise_timeStep, self.dataOut.timeInterval) + #Get signal threshold + signalThresh = noise_multiple*noise + #Meteor echoes detection + listMeteors = self.__findMeteors(powerNet, signalThresh) + #******* END OF NOISE LEVEL AND POSSIBLE METEORS CACULATION ********** + + #************** REMOVE MULTIPLE DETECTIONS (3.5) *************************** + #Parameters + heiRange = self.dataOut.getHeiRange() + rangeInterval = heiRange[1] - heiRange[0] + rangeLimit = multDet_rangeLimit/rangeInterval + timeLimit = multDet_timeLimit/self.dataOut.timeInterval + #Multiple detection removals + listMeteors1 = self.__removeMultipleDetections(listMeteors, rangeLimit, timeLimit) + #************ END OF REMOVE MULTIPLE DETECTIONS ********************** + + #********************* METEOR REESTIMATION (3.7, 3.8, 3.9, 3.10) ******************** + #Parameters + phaseThresh = phaseThresh*numpy.pi/180 + thresh = [phaseThresh, noise_multiple, SNRThresh] + #Meteor reestimation (Errors N 1, 6, 12, 17) + listMeteors2, listMeteorsPower, listMeteorsVolts = self.__meteorReestimation(listMeteors1, voltsPShift, pairslist, thresh, noise, self.dataOut.timeInterval, self.dataOut.frequency) +# listMeteors2, listMeteorsPower, listMeteorsVolts = self.meteorReestimation3(listMeteors2, listMeteorsPower, listMeteorsVolts, voltsPShift, pairslist, thresh, noise) + #Estimation of decay times (Errors N 7, 8, 11) + listMeteors3 = self.__estimateDecayTime(listMeteors2, listMeteorsPower, self.dataOut.timeInterval, self.dataOut.frequency) + #******************* END OF METEOR REESTIMATION ******************* + + #********************* METEOR PARAMETERS CALCULATION (3.11, 3.12, 3.13) ************************** + #Calculating Radial Velocity (Error N 15) + radialStdThresh = 10 + listMeteors4 = self.__getRadialVelocity(listMeteors3, listMeteorsVolts, radialStdThresh, pairslist, self.dataOut.timeInterval) + + if len(listMeteors4) > 0: + #Setting New Array + date = repr(self.dataOut.datatime) + arrayMeteors4, arrayParameters = self.__setNewArrays(listMeteors4, date, heiRang) + + #Calculate AOA (Error N 3, 4) + #JONES ET AL. 1998 + AOAthresh = numpy.pi/8 + error = arrayParameters[:,-1] + phases = -arrayMeteors4[:,9:13] + pairsList = [] + pairsList.append((0,3)) + pairsList.append((1,2)) + arrayParameters[:,4:7], arrayParameters[:,-1] = self.__getAOA(phases, pairsList, error, AOAthresh, azimuth) + + #Calculate Heights (Error N 13 and 14) + error = arrayParameters[:,-1] + Ranges = arrayParameters[:,2] + zenith = arrayParameters[:,5] + arrayParameters[:,3], arrayParameters[:,-1] = self.__getHeights(Ranges, zenith, error, hmin, hmax) + #********************* END OF PARAMETERS CALCULATION ************************** + + #***************************+ SAVE DATA IN HDF5 FORMAT ********************** + self.dataOut.data_param = arrayParameters + + return + + def __getHardwarePhaseDiff(self, voltage0, pairslist, newheis, n): + + minIndex = min(newheis[0]) + maxIndex = max(newheis[0]) + + voltage = voltage0[:,:,minIndex:maxIndex+1] + nLength = voltage.shape[1]/n + nMin = 0 + nMax = 0 + phaseOffset = numpy.zeros((len(pairslist),n)) + + for i in range(n): + nMax += nLength + phaseCCF = -numpy.angle(self.__calculateCCF(voltage[:,nMin:nMax,:], pairslist, [0])) + phaseCCF = numpy.mean(phaseCCF, axis = 2) + phaseOffset[:,i] = phaseCCF.transpose() + nMin = nMax +# phaseDiff, phaseArrival = self.estimatePhaseDifference(voltage, pairslist) + + #Remove Outliers + factor = 2 + wt = phaseOffset - signal.medfilt(phaseOffset,(1,5)) + dw = numpy.std(wt,axis = 1) + dw = dw.reshape((dw.size,1)) + ind = numpy.where(numpy.logical_or(wt>dw*factor,wt<-dw*factor)) + phaseOffset[ind] = numpy.nan + phaseOffset = stats.nanmean(phaseOffset, axis=1) + + return phaseOffset + + def __shiftPhase(self, data, phaseShift): + #this will shift the phase of a complex number + dataShifted = numpy.abs(data) * numpy.exp((numpy.angle(data)+phaseShift)*1j) + return dataShifted + + def __estimatePhaseDifference(self, array, pairslist): + nChannel = array.shape[0] + nHeights = array.shape[2] + numPairs = len(pairslist) +# phaseCCF = numpy.zeros((nChannel, 5, nHeights)) + phaseCCF = numpy.angle(self.__calculateCCF(array, pairslist, [-2,-1,0,1,2])) + + #Correct phases + derPhaseCCF = phaseCCF[:,1:,:] - phaseCCF[:,0:-1,:] + indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi) + + if indDer[0].shape[0] > 0: + for i in range(indDer[0].shape[0]): + signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i],indDer[2][i]]) + phaseCCF[indDer[0][i],indDer[1][i]+1:,:] += signo*2*numpy.pi + +# for j in range(numSides): +# phaseCCFAux = self.calculateCCF(arrayCenter, arraySides[j,:,:], [-2,1,0,1,2]) +# phaseCCF[j,:,:] = numpy.angle(phaseCCFAux) +# + #Linear + phaseInt = numpy.zeros((numPairs,1)) + angAllCCF = phaseCCF[:,[0,1,3,4],0] + for j in range(numPairs): + fit = stats.linregress([-2,-1,1,2],angAllCCF[j,:]) + phaseInt[j] = fit[1] + #Phase Differences + phaseDiff = phaseInt - phaseCCF[:,2,:] + phaseArrival = phaseInt.reshape(phaseInt.size) + + #Dealias + indAlias = numpy.where(phaseArrival > numpy.pi) + phaseArrival[indAlias] -= 2*numpy.pi + indAlias = numpy.where(phaseArrival < -numpy.pi) + phaseArrival[indAlias] += 2*numpy.pi + + return phaseDiff, phaseArrival + + def __coherentDetection(self, volts, timeSegment, timeInterval, pairslist, thresh): + #this function will run the coherent detection used in Holdworth et al. 2004 and return the net power + #find the phase shifts of each channel over 1 second intervals + #only look at ranges below the beacon signal + numProfPerBlock = numpy.ceil(timeSegment/timeInterval) + numBlocks = int(volts.shape[1]/numProfPerBlock) + numHeights = volts.shape[2] + nChannel = volts.shape[0] + voltsCohDet = volts.copy() + + pairsarray = numpy.array(pairslist) + indSides = pairsarray[:,1] +# indSides = numpy.array(range(nChannel)) +# indSides = numpy.delete(indSides, indCenter) +# +# listCenter = numpy.array_split(volts[indCenter,:,:], numBlocks, 0) + listBlocks = numpy.array_split(volts, numBlocks, 1) + + startInd = 0 + endInd = 0 + + for i in range(numBlocks): + startInd = endInd + endInd = endInd + listBlocks[i].shape[1] + + arrayBlock = listBlocks[i] +# arrayBlockCenter = listCenter[i] + + #Estimate the Phase Difference + phaseDiff, aux = self.__estimatePhaseDifference(arrayBlock, pairslist) + #Phase Difference RMS + arrayPhaseRMS = numpy.abs(phaseDiff) + phaseRMSaux = numpy.sum(arrayPhaseRMS < thresh,0) + indPhase = numpy.where(phaseRMSaux==4) + #Shifting + if indPhase[0].shape[0] > 0: + for j in range(indSides.size): + arrayBlock[indSides[j],:,indPhase] = self.__shiftPhase(arrayBlock[indSides[j],:,indPhase], phaseDiff[j,indPhase].transpose()) + voltsCohDet[:,startInd:endInd,:] = arrayBlock + + return voltsCohDet + + def __calculateCCF(self, volts, pairslist ,laglist): + + nHeights = volts.shape[2] + nPoints = volts.shape[1] + voltsCCF = numpy.zeros((len(pairslist), len(laglist), nHeights),dtype = 'complex') + + for i in range(len(pairslist)): + volts1 = volts[pairslist[i][0]] + volts2 = volts[pairslist[i][1]] + + for t in range(len(laglist)): + idxT = laglist[t] + if idxT >= 0: + vStacked = numpy.vstack((volts2[idxT:,:], + numpy.zeros((idxT, nHeights),dtype='complex'))) + else: + vStacked = numpy.vstack((numpy.zeros((-idxT, nHeights),dtype='complex'), + volts2[:(nPoints + idxT),:])) + voltsCCF[i,t,:] = numpy.sum((numpy.conjugate(volts1)*vStacked),axis=0) + + vStacked = None + return voltsCCF + + def __getNoise(self, power, timeSegment, timeInterval): + numProfPerBlock = numpy.ceil(timeSegment/timeInterval) + numBlocks = int(power.shape[0]/numProfPerBlock) + numHeights = power.shape[1] + + listPower = numpy.array_split(power, numBlocks, 0) + noise = numpy.zeros((power.shape[0], power.shape[1])) + noise1 = numpy.zeros((power.shape[0], power.shape[1])) + + startInd = 0 + endInd = 0 + + for i in range(numBlocks): #split por canal + startInd = endInd + endInd = endInd + listPower[i].shape[0] + + arrayBlock = listPower[i] + noiseAux = numpy.mean(arrayBlock, 0) +# noiseAux = numpy.median(noiseAux) +# noiseAux = numpy.mean(arrayBlock) + noise[startInd:endInd,:] = noise[startInd:endInd,:] + noiseAux + + noiseAux1 = numpy.mean(arrayBlock) + noise1[startInd:endInd,:] = noise1[startInd:endInd,:] + noiseAux1 + + return noise, noise1 + + def __findMeteors(self, power, thresh): + nProf = power.shape[0] + nHeights = power.shape[1] + listMeteors = [] + + for i in range(nHeights): + powerAux = power[:,i] + threshAux = thresh[:,i] + + indUPthresh = numpy.where(powerAux > threshAux)[0] + indDNthresh = numpy.where(powerAux <= threshAux)[0] + + j = 0 + + while (j < indUPthresh.size - 2): + if (indUPthresh[j + 2] == indUPthresh[j] + 2): + indDNAux = numpy.where(indDNthresh > indUPthresh[j]) + indDNthresh = indDNthresh[indDNAux] + + if (indDNthresh.size > 0): + indEnd = indDNthresh[0] - 1 + indInit = indUPthresh[j] + + meteor = powerAux[indInit:indEnd + 1] + indPeak = meteor.argmax() + indInit + FLA = sum(numpy.conj(meteor)*numpy.hstack((meteor[1:],0))) + + listMeteors.append(numpy.array([i,indInit,indPeak,indEnd,FLA])) #CHEQUEAR!!!!! + j = numpy.where(indUPthresh == indEnd)[0] + 1 + else: j+=1 + else: j+=1 + + return listMeteors + + def __removeMultipleDetections(self,listMeteors, rangeLimit, timeLimit): + + arrayMeteors = numpy.asarray(listMeteors) + listMeteors1 = [] + + while arrayMeteors.shape[0] > 0: + FLAs = arrayMeteors[:,4] + maxFLA = FLAs.argmax() + listMeteors1.append(arrayMeteors[maxFLA,:]) + + MeteorInitTime = arrayMeteors[maxFLA,1] + MeteorEndTime = arrayMeteors[maxFLA,3] + MeteorHeight = arrayMeteors[maxFLA,0] + + #Check neighborhood + maxHeightIndex = MeteorHeight + rangeLimit + minHeightIndex = MeteorHeight - rangeLimit + minTimeIndex = MeteorInitTime - timeLimit + maxTimeIndex = MeteorEndTime + timeLimit + + #Check Heights + indHeight = numpy.logical_and(arrayMeteors[:,0] >= minHeightIndex, arrayMeteors[:,0] <= maxHeightIndex) + indTime = numpy.logical_and(arrayMeteors[:,3] >= minTimeIndex, arrayMeteors[:,1] <= maxTimeIndex) + indBoth = numpy.where(numpy.logical_and(indTime,indHeight)) + + arrayMeteors = numpy.delete(arrayMeteors, indBoth, axis = 0) + + return listMeteors1 + + def __meteorReestimation(self, listMeteors, volts, pairslist, thresh, noise, timeInterval,frequency): + numHeights = volts.shape[2] + nChannel = volts.shape[0] + + thresholdPhase = thresh[0] + thresholdNoise = thresh[1] + thresholdDB = float(thresh[2]) + + thresholdDB1 = 10**(thresholdDB/10) + pairsarray = numpy.array(pairslist) + indSides = pairsarray[:,1] + + pairslist1 = list(pairslist) + pairslist1.append((0,1)) + pairslist1.append((3,4)) + + listMeteors1 = [] + listPowerSeries = [] + listVoltageSeries = [] + #volts has the war data + + if frequency == 30e6: + timeLag = 45*10**-3 + else: + timeLag = 15*10**-3 + lag = numpy.ceil(timeLag/timeInterval) + + for i in range(len(listMeteors)): + + ###################### 3.6 - 3.7 PARAMETERS REESTIMATION ######################### + meteorAux = numpy.zeros(16) + + #Loading meteor Data (mHeight, mStart, mPeak, mEnd) + mHeight = listMeteors[i][0] + mStart = listMeteors[i][1] + mPeak = listMeteors[i][2] + mEnd = listMeteors[i][3] + + #get the volt data between the start and end times of the meteor + meteorVolts = volts[:,mStart:mEnd+1,mHeight] + meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1) + + #3.6. Phase Difference estimation + phaseDiff, aux = self.__estimatePhaseDifference(meteorVolts, pairslist) + + #3.7. Phase difference removal & meteor start, peak and end times reestimated + #meteorVolts0.- all Channels, all Profiles + meteorVolts0 = volts[:,:,mHeight] + meteorThresh = noise[:,mHeight]*thresholdNoise + meteorNoise = noise[:,mHeight] + meteorVolts0[indSides,:] = self.__shiftPhase(meteorVolts0[indSides,:], phaseDiff) #Phase Shifting + powerNet0 = numpy.nansum(numpy.abs(meteorVolts0)**2, axis = 0) #Power + + #Times reestimation + mStart1 = numpy.where(powerNet0[:mPeak] < meteorThresh[:mPeak])[0] + if mStart1.size > 0: + mStart1 = mStart1[-1] + 1 + + else: + mStart1 = mPeak + + mEnd1 = numpy.where(powerNet0[mPeak:] < meteorThresh[mPeak:])[0][0] + mPeak - 1 + mEndDecayTime1 = numpy.where(powerNet0[mPeak:] < meteorNoise[mPeak:])[0] + if mEndDecayTime1.size == 0: + mEndDecayTime1 = powerNet0.size + else: + mEndDecayTime1 = mEndDecayTime1[0] + mPeak - 1 +# mPeak1 = meteorVolts0[mStart1:mEnd1 + 1].argmax() + + #meteorVolts1.- all Channels, from start to end + meteorVolts1 = meteorVolts0[:,mStart1:mEnd1 + 1] + meteorVolts2 = meteorVolts0[:,mPeak + lag:mEnd1 + 1] + if meteorVolts2.shape[1] == 0: + meteorVolts2 = meteorVolts0[:,mPeak:mEnd1 + 1] + meteorVolts1 = meteorVolts1.reshape(meteorVolts1.shape[0], meteorVolts1.shape[1], 1) + meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1], 1) + ##################### END PARAMETERS REESTIMATION ######################### + + ##################### 3.8 PHASE DIFFERENCE REESTIMATION ######################## +# if mEnd1 - mStart1 > 4: #Error Number 6: echo less than 5 samples long; too short for analysis + if meteorVolts2.shape[1] > 0: + #Phase Difference re-estimation + phaseDiff1, phaseDiffint = self.__estimatePhaseDifference(meteorVolts2, pairslist1) #Phase Difference Estimation +# phaseDiff1, phaseDiffint = self.estimatePhaseDifference(meteorVolts2, pairslist) + meteorVolts2 = meteorVolts2.reshape(meteorVolts2.shape[0], meteorVolts2.shape[1]) + phaseDiff11 = numpy.reshape(phaseDiff1, (phaseDiff1.shape[0],1)) + meteorVolts2[indSides,:] = self.__shiftPhase(meteorVolts2[indSides,:], phaseDiff11[0:4]) #Phase Shifting + + #Phase Difference RMS + phaseRMS1 = numpy.sqrt(numpy.mean(numpy.square(phaseDiff1))) + powerNet1 = numpy.nansum(numpy.abs(meteorVolts1[:,:])**2,0) + #Data from Meteor + mPeak1 = powerNet1.argmax() + mStart1 + mPeakPower1 = powerNet1.max() + noiseAux = sum(noise[mStart1:mEnd1 + 1,mHeight]) + mSNR1 = (sum(powerNet1)-noiseAux)/noiseAux + Meteor1 = numpy.array([mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1]) + Meteor1 = numpy.hstack((Meteor1,phaseDiffint)) + PowerSeries = powerNet0[mStart1:mEndDecayTime1 + 1] + #Vectorize + meteorAux[0:7] = [mHeight, mStart1, mPeak1, mEnd1, mPeakPower1, mSNR1, phaseRMS1] + meteorAux[7:11] = phaseDiffint[0:4] + + #Rejection Criterions + if phaseRMS1 > thresholdPhase: #Error Number 17: Phase variation + meteorAux[-1] = 17 + elif mSNR1 < thresholdDB1: #Error Number 1: SNR < threshold dB + meteorAux[-1] = 1 + + + else: + meteorAux[0:4] = [mHeight, mStart, mPeak, mEnd] + meteorAux[-1] = 6 #Error Number 6: echo less than 5 samples long; too short for analysis + PowerSeries = 0 + + listMeteors1.append(meteorAux) + listPowerSeries.append(PowerSeries) + listVoltageSeries.append(meteorVolts1) + + return listMeteors1, listPowerSeries, listVoltageSeries + + def __estimateDecayTime(self, listMeteors, listPower, timeInterval, frequency): + + threshError = 10 + #Depending if it is 30 or 50 MHz + if frequency == 30e6: + timeLag = 45*10**-3 + else: + timeLag = 15*10**-3 + lag = numpy.ceil(timeLag/timeInterval) + + listMeteors1 = [] + + for i in range(len(listMeteors)): + meteorPower = listPower[i] + meteorAux = listMeteors[i] + + if meteorAux[-1] == 0: + + try: + indmax = meteorPower.argmax() + indlag = indmax + lag + + y = meteorPower[indlag:] + x = numpy.arange(0, y.size)*timeLag + + #first guess + a = y[0] + tau = timeLag + #exponential fit + popt, pcov = optimize.curve_fit(self.__exponential_function, x, y, p0 = [a, tau]) + y1 = self.__exponential_function(x, *popt) + #error estimation + error = sum((y - y1)**2)/(numpy.var(y)*(y.size - popt.size)) + + decayTime = popt[1] + riseTime = indmax*timeInterval + meteorAux[11:13] = [decayTime, error] + + #Table items 7, 8 and 11 + if (riseTime > 0.3): #Number 7: Echo rise exceeds 0.3s + meteorAux[-1] = 7 + elif (decayTime < 2*riseTime) : #Number 8: Echo decay time less than than twice rise time + meteorAux[-1] = 8 + if (error > threshError): #Number 11: Poor fit to amplitude for estimation of decay time + meteorAux[-1] = 11 + + + except: + meteorAux[-1] = 11 + + + listMeteors1.append(meteorAux) + + return listMeteors1 + + #Exponential Function + + def __exponential_function(self, x, a, tau): + y = a*numpy.exp(-x/tau) + return y + + def __getRadialVelocity(self, listMeteors, listVolts, radialStdThresh, pairslist, timeInterval): + + pairslist1 = list(pairslist) + pairslist1.append((0,1)) + pairslist1.append((3,4)) + numPairs = len(pairslist1) + #Time Lag + timeLag = 45*10**-3 + c = 3e8 + lag = numpy.ceil(timeLag/timeInterval) + freq = 30e6 + + listMeteors1 = [] + + for i in range(len(listMeteors)): + meteor = listMeteors[i] + meteorAux = numpy.hstack((meteor[:-1], 0, 0, meteor[-1])) + if meteor[-1] == 0: + mStart = listMeteors[i][1] + mPeak = listMeteors[i][2] + mLag = mPeak - mStart + lag + + #get the volt data between the start and end times of the meteor + meteorVolts = listVolts[i] + meteorVolts = meteorVolts.reshape(meteorVolts.shape[0], meteorVolts.shape[1], 1) + + #Get CCF + allCCFs = self.__calculateCCF(meteorVolts, pairslist1, [-2,-1,0,1,2]) + + #Method 2 + slopes = numpy.zeros(numPairs) + time = numpy.array([-2,-1,1,2])*timeInterval + angAllCCF = numpy.angle(allCCFs[:,[0,1,3,4],0]) + + #Correct phases + derPhaseCCF = angAllCCF[:,1:] - angAllCCF[:,0:-1] + indDer = numpy.where(numpy.abs(derPhaseCCF) > numpy.pi) + + if indDer[0].shape[0] > 0: + for i in range(indDer[0].shape[0]): + signo = -numpy.sign(derPhaseCCF[indDer[0][i],indDer[1][i]]) + angAllCCF[indDer[0][i],indDer[1][i]+1:] += signo*2*numpy.pi + +# fit = scipy.stats.linregress(numpy.array([-2,-1,1,2])*timeInterval, numpy.array([phaseLagN2s[i],phaseLagN1s[i],phaseLag1s[i],phaseLag2s[i]])) + for j in range(numPairs): + fit = stats.linregress(time, angAllCCF[j,:]) + slopes[j] = fit[0] + + #Remove Outlier +# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes))) +# slopes = numpy.delete(slopes,indOut) +# indOut = numpy.argmax(numpy.abs(slopes - numpy.mean(slopes))) +# slopes = numpy.delete(slopes,indOut) + + radialVelocity = -numpy.mean(slopes)*(0.25/numpy.pi)*(c/freq) + radialError = numpy.std(slopes)*(0.25/numpy.pi)*(c/freq) + meteorAux[-2] = radialError + meteorAux[-3] = radialVelocity + + #Setting Error + #Number 15: Radial Drift velocity or projected horizontal velocity exceeds 200 m/s + if numpy.abs(radialVelocity) > 200: + meteorAux[-1] = 15 + #Number 12: Poor fit to CCF variation for estimation of radial drift velocity + elif radialError > radialStdThresh: + meteorAux[-1] = 12 + + listMeteors1.append(meteorAux) + return listMeteors1 + + def __setNewArrays(self, listMeteors, date, heiRang): + + #New arrays + arrayMeteors = numpy.array(listMeteors) + arrayParameters = numpy.zeros((len(listMeteors),10)) + + #Date inclusion + date = re.findall(r'\((.*?)\)', date) + date = date[0].split(',') + date = map(int, date) + date = [date[0]*10000 + date[1]*100 + date[2], date[3]*10000 + date[4]*100 + date[5]] + arrayDate = numpy.tile(date, (len(listMeteors), 1)) + + #Meteor array + arrayMeteors[:,0] = heiRang[arrayMeteors[:,0].astype(int)] + arrayMeteors = numpy.hstack((arrayDate, arrayMeteors)) + + #Parameters Array + arrayParameters[:,0:3] = arrayMeteors[:,0:3] + arrayParameters[:,-3:] = arrayMeteors[:,-3:] + + return arrayMeteors, arrayParameters + + def __getAOA(self, phases, pairsList, error, AOAthresh, azimuth): + + arrayAOA = numpy.zeros((phases.shape[0],3)) + cosdir0, cosdir = self.__getDirectionCosines(phases, pairsList) + + arrayAOA[:,:2] = self.__calculateAOA(cosdir, azimuth) + cosDirError = numpy.sum(numpy.abs(cosdir0 - cosdir), axis = 1) + arrayAOA[:,2] = cosDirError + + azimuthAngle = arrayAOA[:,0] + zenithAngle = arrayAOA[:,1] + + #Setting Error + #Number 3: AOA not fesible + indInvalid = numpy.where(numpy.logical_and((numpy.logical_or(numpy.isnan(zenithAngle), numpy.isnan(azimuthAngle))),error == 0))[0] + error[indInvalid] = 3 + #Number 4: Large difference in AOAs obtained from different antenna baselines + indInvalid = numpy.where(numpy.logical_and(cosDirError > AOAthresh,error == 0))[0] + error[indInvalid] = 4 + return arrayAOA, error + + def __getDirectionCosines(self, arrayPhase, pairsList): + + #Initializing some variables + ang_aux = numpy.array([-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8])*2*numpy.pi + ang_aux = ang_aux.reshape(1,ang_aux.size) + + cosdir = numpy.zeros((arrayPhase.shape[0],2)) + cosdir0 = numpy.zeros((arrayPhase.shape[0],2)) + + + for i in range(2): + #First Estimation + phi0_aux = arrayPhase[:,pairsList[i][0]] + arrayPhase[:,pairsList[i][1]] + #Dealias + indcsi = numpy.where(phi0_aux > numpy.pi) + phi0_aux[indcsi] -= 2*numpy.pi + indcsi = numpy.where(phi0_aux < -numpy.pi) + phi0_aux[indcsi] += 2*numpy.pi + #Direction Cosine 0 + cosdir0[:,i] = -(phi0_aux)/(2*numpy.pi*0.5) + + #Most-Accurate Second Estimation + phi1_aux = arrayPhase[:,pairsList[i][0]] - arrayPhase[:,pairsList[i][1]] + phi1_aux = phi1_aux.reshape(phi1_aux.size,1) + #Direction Cosine 1 + cosdir1 = -(phi1_aux + ang_aux)/(2*numpy.pi*4.5) + + #Searching the correct Direction Cosine + cosdir0_aux = cosdir0[:,i] + cosdir0_aux = cosdir0_aux.reshape(cosdir0_aux.size,1) + #Minimum Distance + cosDiff = (cosdir1 - cosdir0_aux)**2 + indcos = cosDiff.argmin(axis = 1) + #Saving Value obtained + cosdir[:,i] = cosdir1[numpy.arange(len(indcos)),indcos] + + return cosdir0, cosdir + + def __calculateAOA(self, cosdir, azimuth): + cosdirX = cosdir[:,0] + cosdirY = cosdir[:,1] + + zenithAngle = numpy.arccos(numpy.sqrt(1 - cosdirX**2 - cosdirY**2))*180/numpy.pi + azimuthAngle = numpy.arctan2(cosdirX,cosdirY)*180/numpy.pi + azimuth #0 deg north, 90 deg east + angles = numpy.vstack((azimuthAngle, zenithAngle)).transpose() + + return angles + + def __getHeights(self, Ranges, zenith, error, minHeight, maxHeight): + + Ramb = 375 #Ramb = c/(2*PRF) + Re = 6371 #Earth Radius + heights = numpy.zeros(Ranges.shape) + + R_aux = numpy.array([0,1,2])*Ramb + R_aux = R_aux.reshape(1,R_aux.size) + + Ranges = Ranges.reshape(Ranges.size,1) + + Ri = Ranges + R_aux + hi = numpy.sqrt(Re**2 + Ri**2 + (2*Re*numpy.cos(zenith*numpy.pi/180)*Ri.transpose()).transpose()) - Re + + #Check if there is a height between 70 and 110 km + h_bool = numpy.sum(numpy.logical_and(hi > minHeight, hi < maxHeight), axis = 1) + ind_h = numpy.where(h_bool == 1)[0] + + hCorr = hi[ind_h, :] + ind_hCorr = numpy.where(numpy.logical_and(hi > minHeight, hi < maxHeight)) + + hCorr = hi[ind_hCorr] + heights[ind_h] = hCorr + + #Setting Error + #Number 13: Height unresolvable echo: not valid height within 70 to 110 km + #Number 14: Height ambiguous echo: more than one possible height within 70 to 110 km + + indInvalid2 = numpy.where(numpy.logical_and(h_bool > 1, error == 0))[0] + error[indInvalid2] = 14 + indInvalid1 = numpy.where(numpy.logical_and(h_bool == 0, error == 0))[0] + error[indInvalid1] = 13 + + return heights, error + + +class WindProfiler(Operation): + + __isConfig = False + + __initime = None + __lastdatatime = None + __integrationtime = None + + __buffer = None + + __dataReady = False + + __firstdata = None + + n = None + + def __init__(self): + Operation.__init__(self) + + def __calculateAngles(self, theta_x, theta_y, azimuth): + + dir_cosw = numpy.sqrt(1-theta_x**2-theta_y**2) + zenith_arr = numpy.arccos(dir_cosw) + azimuth_arr = numpy.arctan2(theta_x,theta_y) + azimuth*math.pi/180 + + dir_cosu = numpy.sin(azimuth_arr)*numpy.sin(zenith_arr) + dir_cosv = numpy.cos(azimuth_arr)*numpy.sin(zenith_arr) + + return azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw + + def __calculateMatA(self, dir_cosu, dir_cosv, dir_cosw, horOnly): + +# + if horOnly: + A = numpy.c_[dir_cosu,dir_cosv] + else: + A = numpy.c_[dir_cosu,dir_cosv,dir_cosw] + A = numpy.asmatrix(A) + A1 = numpy.linalg.inv(A.transpose()*A)*A.transpose() + + return A1 + + 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 __calculateVelUVW(self, A, velRadial): + + #Operacion Matricial +# velUVW = numpy.zeros((velRadial.shape[1],3)) +# for ind in range(velRadial.shape[1]): +# velUVW[ind,:] = numpy.dot(A,velRadial[:,ind]) +# velUVW = velUVW.transpose() + velUVW = numpy.zeros((A.shape[0],velRadial.shape[1])) + velUVW[:,:] = numpy.dot(A,velRadial) + + + return velUVW + + def techniqueDBS(self, velRadial0, dirCosx, disrCosy, azimuth, correct, horizontalOnly, heiRang, SNR0): + """ + Function that implements Doppler Beam Swinging (DBS) technique. + + Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth, + Direction correction (if necessary), Ranges and SNR + + Output: Winds estimation (Zonal, Meridional and Vertical) + + Parameters affected: Winds, height range, SNR + """ + azimuth_arr, zenith_arr, dir_cosu, dir_cosv, dir_cosw = self.__calculateAngles(dirCosx, disrCosy, azimuth) + heiRang1, velRadial1, SNR1 = self.__correctValues(heiRang, zenith_arr, correct*velRadial0, SNR0) + A = self.__calculateMatA(dir_cosu, dir_cosv, dir_cosw, horizontalOnly) + + #Calculo de Componentes de la velocidad con DBS + winds = self.__calculateVelUVW(A,velRadial1) + + return winds, heiRang1, SNR1 + + def __calculateDistance(self, posx, posy, pairsCrossCorr, pairsList, pairs, azimuth = None): + + posx = numpy.asarray(posx) + posy = numpy.asarray(posy) + + #Rotacion Inversa para alinear con el azimuth + if azimuth!= None: + azimuth = azimuth*math.pi/180 + posx1 = posx*math.cos(azimuth) + posy*math.sin(azimuth) + posy1 = -posx*math.sin(azimuth) + posy*math.cos(azimuth) + else: + posx1 = posx + posy1 = posy + + #Calculo de Distancias + distx = numpy.zeros(pairsCrossCorr.size) + disty = numpy.zeros(pairsCrossCorr.size) + dist = numpy.zeros(pairsCrossCorr.size) + ang = numpy.zeros(pairsCrossCorr.size) + + for i in range(pairsCrossCorr.size): + distx[i] = posx1[pairsList[pairsCrossCorr[i]][1]] - posx1[pairsList[pairsCrossCorr[i]][0]] + disty[i] = posy1[pairsList[pairsCrossCorr[i]][1]] - posy1[pairsList[pairsCrossCorr[i]][0]] + dist[i] = numpy.sqrt(distx[i]**2 + disty[i]**2) + ang[i] = numpy.arctan2(disty[i],distx[i]) + #Calculo de Matrices + nPairs = len(pairs) + ang1 = numpy.zeros((nPairs, 2, 1)) + dist1 = numpy.zeros((nPairs, 2, 1)) + + for j in range(nPairs): + dist1[j,0,0] = dist[pairs[j][0]] + dist1[j,1,0] = dist[pairs[j][1]] + ang1[j,0,0] = ang[pairs[j][0]] + ang1[j,1,0] = ang[pairs[j][1]] + + return distx,disty, dist1,ang1 + + def __calculateVelVer(self, phase, lagTRange, _lambda): + + Ts = lagTRange[1] - lagTRange[0] + velW = -_lambda*phase/(4*math.pi*Ts) + + return velW + + def __calculateVelHorDir(self, dist, tau1, tau2, ang): + nPairs = tau1.shape[0] + vel = numpy.zeros((nPairs,3,tau1.shape[2])) + + angCos = numpy.cos(ang) + angSin = numpy.sin(ang) + + vel0 = dist*tau1/(2*tau2**2) + vel[:,0,:] = (vel0*angCos).sum(axis = 1) + vel[:,1,:] = (vel0*angSin).sum(axis = 1) + + ind = numpy.where(numpy.isinf(vel)) + vel[ind] = numpy.nan + + return vel + + def __getPairsAutoCorr(self, pairsList, nChannels): + + pairsAutoCorr = numpy.zeros(nChannels, dtype = 'int')*numpy.nan + + for l in range(len(pairsList)): + firstChannel = pairsList[l][0] + secondChannel = pairsList[l][1] + + #Obteniendo pares de Autocorrelacion + if firstChannel == secondChannel: + pairsAutoCorr[firstChannel] = int(l) + + pairsAutoCorr = pairsAutoCorr.astype(int) + + pairsCrossCorr = range(len(pairsList)) + pairsCrossCorr = numpy.delete(pairsCrossCorr,pairsAutoCorr) + + return pairsAutoCorr, pairsCrossCorr + + def techniqueSA(self, pairsSelected, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, lagTRange, correctFactor): + """ + Function that implements Spaced Antenna (SA) technique. + + Input: Radial velocities, Direction cosines (x and y) of the Beam, Antenna azimuth, + Direction correction (if necessary), Ranges and SNR + + Output: Winds estimation (Zonal, Meridional and Vertical) + + Parameters affected: Winds + """ + #Cross Correlation pairs obtained + pairsAutoCorr, pairsCrossCorr = self.__getPairsAutoCorr(pairsList, nChannels) + pairsArray = numpy.array(pairsList)[pairsCrossCorr] + pairsSelArray = numpy.array(pairsSelected) + pairs = [] + + #Wind estimation pairs obtained + for i in range(pairsSelArray.shape[0]/2): + ind1 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i], axis = 1))[0][0] + ind2 = numpy.where(numpy.all(pairsArray == pairsSelArray[2*i + 1], axis = 1))[0][0] + pairs.append((ind1,ind2)) + + indtau = tau.shape[0]/2 + tau1 = tau[:indtau,:] + tau2 = tau[indtau:-1,:] + tau1 = tau1[pairs,:] + tau2 = tau2[pairs,:] + phase1 = tau[-1,:] + + #--------------------------------------------------------------------- + #Metodo Directo + distx, disty, dist, ang = self.__calculateDistance(position_x, position_y, pairsCrossCorr, pairsList, pairs,azimuth) + winds = self.__calculateVelHorDir(dist, tau1, tau2, ang) + winds = stats.nanmean(winds, axis=0) + #--------------------------------------------------------------------- + #Metodo General +# distx, disty, dist = self.calculateDistance(position_x,position_y,pairsCrossCorr, pairsList, azimuth) +# #Calculo Coeficientes de Funcion de Correlacion +# F,G,A,B,H = self.calculateCoef(tau1,tau2,distx,disty,n) +# #Calculo de Velocidades +# winds = self.calculateVelUV(F,G,A,B,H) + + #--------------------------------------------------------------------- + winds[2,:] = self.__calculateVelVer(phase1, lagTRange, _lambda) + winds = correctFactor*winds + return winds + + def __checkTime(self, currentTime, paramInterval, windsInterval): + + dataTime = currentTime + paramInterval + deltaTime = dataTime - self.__initime + + if deltaTime >= windsInterval or deltaTime < 0: + self.__dataReady = True + return + + def techniqueMeteors(self, arrayMeteor, meteorThresh, heightMin, heightMax): + ''' + Function that implements winds estimation technique with detected meteors. + + Input: Detected meteors, Minimum meteor quantity to wind estimation + + Output: Winds estimation (Zonal and Meridional) + + Parameters affected: Winds + ''' + #Settings + nInt = (heightMax - heightMin)/2 + winds = numpy.zeros((2,nInt))*numpy.nan + + #Filter errors + error = numpy.where(arrayMeteor[:,-1] == 0)[0] + finalMeteor = arrayMeteor[error,:] + + #Meteor Histogram + finalHeights = finalMeteor[:,3] + hist = numpy.histogram(finalHeights, bins = nInt, range = (heightMin,heightMax)) + nMeteorsPerI = hist[0] + heightPerI = hist[1] + + #Sort of meteors + indSort = finalHeights.argsort() + finalMeteor2 = finalMeteor[indSort,:] + + # Calculating winds + ind1 = 0 + ind2 = 0 + + for i in range(nInt): + nMet = nMeteorsPerI[i] + ind1 = ind2 + ind2 = ind1 + nMet + + meteorAux = finalMeteor2[ind1:ind2,:] + + if meteorAux.shape[0] >= meteorThresh: + vel = meteorAux[:, 7] + zen = meteorAux[:, 5]*numpy.pi/180 + azim = meteorAux[:, 4]*numpy.pi/180 + + n = numpy.cos(zen) + # m = (1 - n**2)/(1 - numpy.tan(azim)**2) + # l = m*numpy.tan(azim) + l = numpy.sin(zen)*numpy.sin(azim) + m = numpy.sin(zen)*numpy.cos(azim) + + A = numpy.vstack((l, m)).transpose() + A1 = numpy.dot(numpy.linalg.inv( numpy.dot(A.transpose(),A) ),A.transpose()) + windsAux = numpy.dot(A1, vel) + + winds[0,i] = windsAux[0] + winds[1,i] = windsAux[1] + + return winds, heightPerI[:-1] + + def run(self, dataOut, technique, **kwargs): + + param = dataOut.data_param + if dataOut.abscissaRange != None: + absc = dataOut.abscissaRange[:-1] + noise = dataOut.noise + heightRange = dataOut.getHeiRange() + SNR = dataOut.SNR + + if technique == 'DBS': + + theta_x = numpy.array(kwargs['dirCosx']) + theta_y = numpy.array(kwargs['dirCosy']) + azimuth = kwargs['azimuth'] + if kwargs.has_key('horizontalOnly'): + horizontalOnly = kwargs['horizontalOnly'] + else: horizontalOnly = False + if kwargs.has_key('correctFactor'): + correctFactor = kwargs['correctFactor'] + else: correctFactor = 1 + if kwargs.has_key('channelList'): + channelList = kwargs['channelList'] + if len(channelList) == 2: + horizontalOnly = True + arrayChannel = numpy.array(channelList) + param = param[arrayChannel,:,:] + theta_x = theta_x[arrayChannel] + 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.initUtcTime = dataOut.ltctime + dataOut.windsInterval = dataOut.timeInterval + + elif technique == 'SA': + + #Parameters + position_x = kwargs['positionX'] + position_y = kwargs['positionY'] + azimuth = kwargs['azimuth'] + + if kwargs.has_key('crosspairsList'): + pairs = kwargs['crosspairsList'] + else: + pairs = None + + if kwargs.has_key('correctFactor'): + correctFactor = kwargs['correctFactor'] + else: + correctFactor = 1 + + tau = dataOut.data_param + _lambda = dataOut.C/dataOut.frequency + pairsList = dataOut.pairsList + nChannels = dataOut.nChannels + + dataOut.winds = self.techniqueSA(pairs, pairsList, nChannels, tau, azimuth, _lambda, position_x, position_y, absc, correctFactor) + dataOut.initUtcTime = dataOut.ltctime + dataOut.windsInterval = dataOut.timeInterval + + elif technique == 'Meteors': + dataOut.flagNoData = True + self.__dataReady = False + + if kwargs.has_key('nHours'): + nHours = kwargs['nHours'] + else: + nHours = 1 + + if kwargs.has_key('meteorsPerBin'): + meteorThresh = kwargs['meteorsPerBin'] + else: + meteorThresh = 6 + + if kwargs.has_key('hmin'): + hmin = kwargs['hmin'] + else: hmin = 70 + if kwargs.has_key('hmax'): + hmax = kwargs['hmax'] + else: hmax = 110 + + dataOut.windsInterval = nHours*3600 + + if self.__isConfig == False: +# self.__initime = dataOut.datatime.replace(minute = 0, second = 0, microsecond = 03) + #Get Initial LTC time + self.__initime = (dataOut.datatime.replace(minute = 0, second = 0, microsecond = 0) - datetime.datetime(1970, 1, 1)).total_seconds() + self.__isConfig = True + + if self.__buffer == None: + self.__buffer = dataOut.data_param + self.__firstdata = copy.copy(dataOut) + + else: + self.__buffer = numpy.vstack((self.__buffer, dataOut.data_param)) + + self.__checkTime(dataOut.ltctime, dataOut.paramInterval, dataOut.windsInterval) #Check if the buffer is ready + + if self.__dataReady: + dataOut.initUtcTime = self.__initime + self.__initime = self.__initime + dataOut.windsInterval #to erase time offset + + dataOut.winds, dataOut.heightRange = self.techniqueMeteors(self.__buffer, meteorThresh, hmin, hmax) + dataOut.flagNoData = False + self.__buffer = None + + return \ No newline at end of file diff --git a/schainpy/model/proc/jroproc_spectra.py b/schainpy/model/proc/jroproc_spectra.py index 58e9ea4..4384a60 100644 --- a/schainpy/model/proc/jroproc_spectra.py +++ b/schainpy/model/proc/jroproc_spectra.py @@ -1,4 +1,5 @@ import numpy +import math from jroproc_base import ProcessingUnit, Operation from model.data.jrodata import Spectra diff --git a/schainpy/model/proc/jroprocessing.py b/schainpy/model/proc/jroprocessing.py index ab197bb..00333fc 100644 --- a/schainpy/model/proc/jroprocessing.py +++ b/schainpy/model/proc/jroprocessing.py @@ -1,4 +1,6 @@ from jroproc_voltage import * from jroproc_spectra import * from jroproc_heispectra import * -from jroproc_amisr import * \ No newline at end of file +from jroproc_amisr import * +from jroproc_correlation import * +from jroproc_parameters import * \ No newline at end of file diff --git a/schainpy/test/Meteor_JASMET30MHz_Winds.py b/schainpy/test/Meteor_JASMET30MHz_Winds.py new file mode 100644 index 0000000..334aa1c --- /dev/null +++ b/schainpy/test/Meteor_JASMET30MHz_Winds.py @@ -0,0 +1,103 @@ +# 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 = "JASMET Experiment Test" +filename = "JASMETtest.xml" + +controllerObj = Project() + +controllerObj.setup(id = '191', name='test01', description=desc) + +#Experimentos + +#2014051 20 Feb 2014 +path = '/home/soporte/Data/JASMET/JASMET_30/2014106' +pathFigure = '/home/soporte/workspace/Graficos/JASMET/prueba1' + +startTime = '00:00:00' +endTime = '23:59:59' +xmin ='19.0' +xmax = '34.0' + +#------------------------------------------------------------------------------------------------ +readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader', + path=path, + startDate='2014/04/15', + endDate='2014/04/16', + startTime=startTime, + endTime=endTime, + online=0, + delay=5, + walk=0) + +opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock') + + +#-------------------------------------------------------------------------------------------------- + +procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId()) + +opObj00 = procUnitConfObj0.addOperation(name='selectChannels') +opObj00.addParameter(name='channelList', value='0, 1, 2, 3, 4', format='intlist') + +opObj01 = procUnitConfObj0.addOperation(name='setRadarFrequency') +opObj01.addParameter(name='frequency', value='30.e6', format='float') + +#opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other') +#-------------------------------------------------------------------------------------------------- + +procUnitConfObj1 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj0.getId()) +procUnitConfObj1.addParameter(name='nSeconds', value='100', format='int') + +opObj10 = procUnitConfObj1.addOperation(name='DetectMeteors') +opObj10.addParameter(name='predefinedPhaseShifts', value='-89.5, 41.5, 0.0, -138.0, -85.5', format='floatlist') +opObj10.addParameter(name='cohDetection', value='0', format='bool') +opObj10.addParameter(name='noise_multiple', value='4', format='int') +opObj10.addParameter(name='SNRThresh', value='5', format='float') +opObj10.addParameter(name='phaseThresh', value='20', format='float') +opObj10.addParameter(name='azimuth', value='45', format='float') +opObj10.addParameter(name='hmin', value='68', format='float') +opObj10.addParameter(name='hmax', value='112', format='float') + +opObj13 = procUnitConfObj1.addOperation(name='SkyMapPlot', optype='other') +opObj13.addParameter(name='id', value='1', format='int') +opObj13.addParameter(name='wintitle', value='Sky Map', format='str') +opObj13.addParameter(name='save', value='1', format='bool') +opObj13.addParameter(name='figpath', value=pathFigure, format='str') + +opObj11 = procUnitConfObj1.addOperation(name='WindProfiler', optype='other') +opObj11.addParameter(name='technique', value='Meteors', format='str') +opObj11.addParameter(name='nHours', value='1.0', format='float') +opObj11.addParameter(name='SNRThresh', value='12.0', format='float') + +opObj12 = procUnitConfObj1.addOperation(name='WindProfilerPlot', optype='other') +opObj12.addParameter(name='id', value='2', format='int') +opObj12.addParameter(name='wintitle', value='Wind Profiler', format='str') +opObj12.addParameter(name='save', value='1', format='bool') +opObj12.addParameter(name='figpath', value = pathFigure, format='str') +opObj12.addParameter(name='zmin', value='-120', format='int') +opObj12.addParameter(name='zmax', value='120', format='int') +# opObj12.addParameter(name='zmin_ver', value='-0.8', format='float') +# opObj12.addParameter(name='zmax_ver', value='0.8', format='float') +# opObj23.addParameter(name='SNRmin', value='-10', format='int') +# opObj23.addParameter(name='SNRmax', value='60', format='int') +# opObj23.addParameter(name='SNRthresh', value='0', format='float') +opObj12.addParameter(name='xmin', value=xmin, format='float') +opObj12.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 diff --git a/schainpy/test/Meteor_JASMET_30mhz.py b/schainpy/test/Meteor_JASMET30Mhz_Beacon.py similarity index 67% rename from schainpy/test/Meteor_JASMET_30mhz.py rename to schainpy/test/Meteor_JASMET30Mhz_Beacon.py index b6eb250..a56aa48 100644 --- a/schainpy/test/Meteor_JASMET_30mhz.py +++ b/schainpy/test/Meteor_JASMET30Mhz_Beacon.py @@ -15,39 +15,42 @@ controllerObj = Project() controllerObj.setup(id = '191', name='meteor_test01', description=desc) path = '/home/dsuarez/.gvfs/data on 10.10.20.13/Jasmet50' +pathFigure = '/home/jasmet/jasmet30_phase' +path = '/home/soporte/Data/JASMET/JASMET_30/2014106' +pathFigure = '/home/soporte/workspace/Graficos/JASMET/prueba1' -readUnitConfObj = controllerObj.addReadUnit(datatype='Voltage', +readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader', path=path, startDate='2014/04/15', endDate='2014/04/15', - startTime='17:00:00', + startTime='20:00:00', endTime='23:59:59', online=0, - walk=1) + walk=0) opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock') -procUnitConfObj0 = controllerObj.addProcUnit(datatype='Voltage', inputId=readUnitConfObj.getId()) +procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId()) -opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other') +# opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other') +# +# +# opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other') +# opObj11.addParameter(name='n', value='2', format='int') - -opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other') -opObj11.addParameter(name='n', value='2', format='int') - -opObj11 = procUnitConfObj0.addOperation(name='VoltageWriter', optype='other') -opObj11.addParameter(name='path', value='/home/jasmet/jasmet30_abril') -opObj11.addParameter(name='blocksPerFile', value='100', format='int') -opObj11.addParameter(name='profilesPerBlock', value='200', format='int') +# opObj11 = procUnitConfObj0.addOperation(name='VoltageWriter', optype='other') +# opObj11.addParameter(name='path', value='/home/jasmet/jasmet30_abril') +# opObj11.addParameter(name='blocksPerFile', value='100', format='int') +# opObj11.addParameter(name='profilesPerBlock', value='200', format='int') """ ########################################### BEACON ########################################## """ -procUnitConfObjBeacon = controllerObj.addProcUnit(datatype='Spectra', inputId=procUnitConfObj0.getId()) +procUnitConfObjBeacon = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId()) procUnitConfObjBeacon.addParameter(name='nProfiles', value='200', format='int') procUnitConfObjBeacon.addParameter(name='nFFTPoints', value='200', format='int') procUnitConfObjBeacon.addParameter(name='pairsList', value='(2,0),(2,1),(2,3),(2,4)', format='pairsList') @@ -61,11 +64,11 @@ opObj11 = procUnitConfObjBeacon.addOperation(name='BeaconPhase', optype='other') opObj11.addParameter(name='id', value='201', format='int') opObj11.addParameter(name='wintitle', value='Beacon Phase', format='str') opObj11.addParameter(name='timerange', value='300', format='int') -opObj11.addParameter(name='xmin', value='0', format='float') -opObj11.addParameter(name='xmax', value='24', format='float') +opObj11.addParameter(name='xmin', value='20', format='float') +opObj11.addParameter(name='xmax', value='22', format='float') opObj11.addParameter(name='ymin', value='-180', format='float') opObj11.addParameter(name='ymax', value='180', format='float') -opObj11.addParameter(name='figpath', value='/home/jasmet/jasmet30_phase', format='str') +opObj11.addParameter(name='figpath', value=pathFigure, format='str') print "Escribiendo el archivo XML" diff --git a/schainpy/test/Meteor_JASMET50MHz_Winds.py b/schainpy/test/Meteor_JASMET50MHz_Winds.py new file mode 100644 index 0000000..6dc9ffc --- /dev/null +++ b/schainpy/test/Meteor_JASMET50MHz_Winds.py @@ -0,0 +1,105 @@ +# 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 = "JASMET Experiment Test" +filename = "JASMETtest.xml" + +controllerObj = Project() + +controllerObj.setup(id = '191', name='test01', description=desc) + +#Experimentos + +#2014051 20 Feb 2014 +path = '/home/soporte/Data/JASMET/JASMET_50/2014106' +pathFigure = '/home/soporte/workspace/Graficos/JASMET' + +startTime = '00:00:00' +endTime = '23:59:59' +xmin ='19.0' +xmax = '33.0' + + +#------------------------------------------------------------------------------------------------ +readUnitConfObj = controllerObj.addReadUnit(datatype='Voltage', + path=path, + startDate='2014/04/15', + endDate='2014/04/16', + startTime=startTime, + endTime=endTime, + online=0, + delay=5, + walk=0) + +opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock') + + +#-------------------------------------------------------------------------------------------------- + +procUnitConfObj0 = controllerObj.addProcUnit(datatype='Voltage', inputId=readUnitConfObj.getId()) + +opObj00 = procUnitConfObj0.addOperation(name='selectChannels') +opObj00.addParameter(name='channelList', value='0, 1, 2, 3, 4', format='intlist') + +opObj01 = procUnitConfObj0.addOperation(name='setRadarFrequency') +opObj01.addParameter(name='frequency', value='50.e6', format='float') + +#opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other') +#-------------------------------------------------------------------------------------------------- + +procUnitConfObj1 = controllerObj.addProcUnit(datatype='Parameters', inputId=procUnitConfObj0.getId()) +procUnitConfObj1.addParameter(name='nSeconds', value='100', format='int') + +opObj10 = procUnitConfObj1.addOperation(name='DetectMeteors') +opObj10.addParameter(name='predefinedPhaseShifts', value='-17.0, -139.0, 0.0, 48.0, -130.0', format='floatlist') +opObj10.addParameter(name='cohDetection', value='0', format='bool') +opObj10.addParameter(name='noise_multiple', value='4', format='int') +opObj10.addParameter(name='SNRThresh', value='5', format='float') +opObj10.addParameter(name='phaseThresh', value='20', format='float') +opObj10.addParameter(name='azimuth', value='45', format='float') +opObj10.addParameter(name='hmin', value='68', format='float') +opObj10.addParameter(name='hmax', value='112', format='float') +opObj10.addParameter(name='savefile', value='1', format='bool') +opObj10.addParameter(name='filename', value=filehdf5, format='str') + +opObj13 = procUnitConfObj1.addOperation(name='SkyMapPlot', optype='other') +opObj13.addParameter(name='id', value='1', format='int') +opObj13.addParameter(name='wintitle', value='Sky Map', format='str') +opObj13.addParameter(name='save', value='1', format='bool') +opObj13.addParameter(name='figpath', value=pathFigure, format='str') + +opObj11 = procUnitConfObj1.addOperation(name='WindProfiler', optype='other') +opObj11.addParameter(name='technique', value='Meteors', format='str') +opObj11.addParameter(name='nHours', value='1.0', format='float') + +opObj12 = procUnitConfObj1.addOperation(name='WindProfilerPlot', optype='other') +opObj12.addParameter(name='id', value='2', format='int') +opObj12.addParameter(name='wintitle', value='Wind Profiler', format='str') +opObj12.addParameter(name='save', value='1', format='bool') +opObj12.addParameter(name='figpath', value = pathFigure, format='str') +opObj12.addParameter(name='zmin', value='-120', format='int') +opObj12.addParameter(name='zmax', value='120', format='int') +# opObj12.addParameter(name='zmin_ver', value='-0.8', format='float') +# opObj12.addParameter(name='zmax_ver', value='0.8', format='float') +# opObj23.addParameter(name='SNRmin', value='-10', format='int') +# opObj23.addParameter(name='SNRmax', value='60', format='int') +# opObj23.addParameter(name='SNRthresh', value='0', format='float') +opObj12.addParameter(name='xmin', value=xmin, format='float') +opObj12.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 diff --git a/schainpy/test/Meteor_JASMET_50mhz.py b/schainpy/test/Meteor_JASMET50Mhz_Beacon.py similarity index 100% rename from schainpy/test/Meteor_JASMET_50mhz.py rename to schainpy/test/Meteor_JASMET50Mhz_Beacon.py diff --git a/schainpy/test/WindProfiler_DBS01.py b/schainpy/test/WindProfiler_DBS01.py new file mode 100644 index 0000000..d007ae4 --- /dev/null +++ b/schainpy/test/WindProfiler_DBS01.py @@ -0,0 +1,144 @@ +# 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 + +#2014050 19 Feb 2014 +# path = '/home/soporte/Documents/MST_Data/DBS/d2014050' +# pathFigure = '/home/soporte/workspace/Graficos/DBS/d2014050p/' +# xmin = '15.5' +# xmax = '23.99999999' +# startTime = '17:25:00' +# filehdf5 = "DBS_2014050.hdf5" + +#2014051 20 Feb 2014 +path = '/home/soporte/Data/MST/DBS/d2014051' +pathFigure = '/home/soporte/workspace/Graficos/DBS/prueba1/' +xmin = '0.0' +xmax = '8.0' +startTime = '00:00:00' +filehdf5 = "DBS_2014051.hdf5" + + + +#------------------------------------------------------------------------------------------------ +readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader', + path=path, + startDate='2014/01/31', + endDate='2014/03/31', + startTime=startTime, + endTime='23:59:59', + online=0, + delay=5, + walk=0) + +opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock') + + +#-------------------------------------------------------------------------------------------------- + +procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId()) + +opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other') + +opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other') +opObj11.addParameter(name='n', value='256', format='int') +# opObj11.addParameter(name='n', value='16', format='int') + +opObj11 = procUnitConfObj0.addOperation(name='selectHeightsByIndex') +opObj11.addParameter(name='minIndex', value='10', format='float') +opObj11.addParameter(name='maxIndex', value='60', format='float') + +#--------------------------------------------------------------------------------------------------- + +procUnitConfObj1 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObj0.getId()) +procUnitConfObj1.addParameter(name='nFFTPoints', value='64', format='int') +procUnitConfObj1.addParameter(name='nProfiles', value='64', format='int') +# procUnitConfObj1.addParameter(name='ippFactor', value='2', format='int') +procUnitConfObj1.addParameter(name='pairsList', value='(0,0),(0,1),(2,1)', format='pairsList') + +opObj11 = procUnitConfObj1.addOperation(name='IncohInt', optype='other') +opObj11.addParameter(name='n', value='5', format='int') + +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='90', format='int') + +opObj12 = procUnitConfObj1.addOperation(name='removeInterference') +opObj13 = procUnitConfObj1.addOperation(name='removeDC') +opObj13.addParameter(name='mode', value='1', 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='90', format='int') + +#-------------------------------------------------------------------------------------------------- + +procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId()) +opObj20 = procUnitConfObj2.addOperation(name='GetMoments') + +opObj21 = procUnitConfObj2.addOperation(name='MomentsPlot', optype='other') +opObj21.addParameter(name='id', value='3', format='int') +opObj21.addParameter(name='wintitle', value='Moments Plot', format='str') +opObj21.addParameter(name='save', value='1', format='bool') +opObj21.addParameter(name='figpath', value=pathFigure, format='str') +opObj21.addParameter(name='zmin', value='5', format='int') +opObj21.addParameter(name='zmax', value='90', format='int') + +opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other') +opObj22.addParameter(name='technique', value='DBS', format='str') +opObj22.addParameter(name='azimuth', value='51.06', format='float') +opObj22.addParameter(name='correctFactor', value='-1', format='float') +opObj22.addParameter(name='dirCosx', value='0.041016, 0, -0.054688', format='floatlist') +opObj22.addParameter(name='dirCosy', value='-0.041016, 0.025391, -0.023438', format='floatlist') +# opObj22.addParameter(name='horizontalOnly', value='1', format='bool') +# opObj22.addParameter(name='channelList', value='1,2', format='intlist') + +opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other') +opObj23.addParameter(name='id', value='4', format='int') +opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str') +opObj23.addParameter(name='save', value='1', format='bool') +opObj23.addParameter(name='figpath', value = pathFigure, format='str') +opObj23.addParameter(name='zmin', value='-10', format='int') +opObj23.addParameter(name='zmax', value='10', format='int') +opObj23.addParameter(name='zmin_ver', value='-80', format='float') +opObj23.addParameter(name='zmax_ver', value='80', format='float') +opObj23.addParameter(name='SNRmin', value='-10', format='int') +opObj23.addParameter(name='SNRmax', value='60', format='int') +opObj23.addParameter(name='SNRthresh', value='0', 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 diff --git a/schainpy/test/WindProfiler_SA01.py b/schainpy/test/WindProfiler_SA01.py new file mode 100644 index 0000000..15babd8 --- /dev/null +++ b/schainpy/test/WindProfiler_SA01.py @@ -0,0 +1,139 @@ +# 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 = "SA Experiment Test" +filename = "SA2014050.xml" + +controllerObj = Project() + +controllerObj.setup(id = '191', name='test01', description=desc) + + +#Experimentos + +#2014050 19 Feb 2014 +# path = '/home/soporte/Documents/MST_Data/SA/d2014050' +# pathFigure = '/home/soporte/workspace/Graficos/SA/d2014050_prueba/' +# xmin = '15.5' +# xmax = '23.99999999' +# startTime = '15:30:00' +# filehdf5 = "SA_2014050.hdf5" + +#2014051 20 Feb 2014 +path = '/home/soporte/Data/MST/SA/d2014051' +pathFigure = '/home/soporte/workspace/Graficos/SA/prueba1/' +xmin = '0.0' +xmax = '8.0' +startTime = '06:00:00' +filehdf5 = "SA_2014051.hdf5" + +readUnitConfObj = controllerObj.addReadUnit(datatype='VoltageReader', + path=path, + startDate='2014/01/01', + endDate='2014/03/31', + startTime=startTime, + endTime='23:59:59', + online=0, + delay=5, + walk=0) + +opObj11 = readUnitConfObj.addOperation(name='printNumberOfBlock') + + +#-------------------------------------------------------------------------------------------------- + +procUnitConfObj0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=readUnitConfObj.getId()) + +opObj11 = procUnitConfObj0.addOperation(name='Decoder', optype='other') + +opObj11 = procUnitConfObj0.addOperation(name='CohInt', optype='other') +opObj11.addParameter(name='n', value='600', format='int') +# opObj11.addParameter(name='n', value='10', format='int') + +opObj11 = procUnitConfObj0.addOperation(name='selectHeightsByIndex') +opObj11.addParameter(name='minIndex', value='10', format='float') +opObj11.addParameter(name='maxIndex', value='60', format='float') +#--------------------------------------------------------------------------------------------------- +procUnitConfObj1 = controllerObj.addProcUnit(datatype='CorrelationProc', inputId=procUnitConfObj0.getId()) +# procUnitConfObj1.addParameter(name='pairsList', value='(0,0),(1,1),(2,2),(3,3),(1,0),(2,3)', format='pairsList') +procUnitConfObj1.addParameter(name='pairsList', value='(0,0),(1,1),(2,2),(3,3),(0,3),(0,2),(1,3),(1,2)', format='pairsList') +procUnitConfObj1.addParameter(name='fullT', value='1', format='bool') +procUnitConfObj1.addParameter(name='removeDC', value='1', format='bool') +#procUnitConfObj1.addParameter(name='lagT', value='0,1,2,3', format='intlist') + +opObj12 = procUnitConfObj1.addOperation(name='CorrelationPlot', optype='other') +opObj12.addParameter(name='id', value='1', format='int') +opObj12.addParameter(name='wintitle', value='CrossCorrelation Plot', format='str') +opObj12.addParameter(name='save', value='1', format='bool') +opObj12.addParameter(name='zmin', value='0', format='int') +opObj12.addParameter(name='zmax', value='1', format='int') +opObj12.addParameter(name='figpath', value = pathFigure, format='str') + +opObj12 = procUnitConfObj1.addOperation(name='removeNoise') +opObj12.addParameter(name='mode', value='2', format='int') +opObj12 = procUnitConfObj1.addOperation(name='calculateNormFactor') + +opObj12 = procUnitConfObj1.addOperation(name='CorrelationPlot', optype='other') +opObj12.addParameter(name='id', value='2', format='int') +opObj12.addParameter(name='wintitle', value='CrossCorrelation Plot', format='str') +opObj12.addParameter(name='save', value='1', format='bool') +opObj12.addParameter(name='zmin', value='0', format='int') +opObj12.addParameter(name='zmax', value='1', format='int') +opObj12.addParameter(name='figpath', value = pathFigure, format='str') + +#--------------------------------------------------------------------------------------------------- +procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObj1.getId()) +opObj20 = procUnitConfObj2.addOperation(name='GetLags') + +opObj21 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other') +opObj21.addParameter(name='technique', value='SA', format='str') +# opObj21.addParameter(name='correctFactor', value='-1', format='float') +opObj21.addParameter(name='positionX', value='36,0,36,0', format='floatlist') +opObj21.addParameter(name='positionY', value='36,0,0,36', format='floatlist') +opObj21.addParameter(name='azimuth', value='51.06', format='float') +opObj21.addParameter(name='crosspairsList', value='(0,3),(0,2),(1,3),(1,2)', format='pairsList')#COrregir +# +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='figpath', value = pathFigure, format='str') +opObj22.addParameter(name='zmin', value='-15', format='int') +opObj22.addParameter(name='zmax', value='15', format='int') +opObj22.addParameter(name='zmin_ver', value='-80', format='float') +opObj22.addParameter(name='zmax_ver', value='80', format='float') +opObj22.addParameter(name='SNRmin', value='-20', format='int') +opObj22.addParameter(name='SNRmax', value='40', format='int') +opObj22.addParameter(name='SNRthresh', value='-3.5', format='float') +opObj22.addParameter(name='xmin', value=xmin, format='float') +opObj22.addParameter(name='xmax', value=xmax, format='float') +# #----------------------------------------------------------------------------------- +# +# procUnitConfObj2 = controllerObj.addProcUnit(datatype='Spectra', inputId=procUnitConfObj0.getId()) +# procUnitConfObj2.addParameter(name='nFFTPoints', value='128', format='int') +# procUnitConfObj2.addParameter(name='nProfiles', value='128', format='int') +# procUnitConfObj2.addParameter(name='pairsList', value='(0,0),(0,1),(2,1)', format='pairsList') +# +# opObj22 = procUnitConfObj2.addOperation(name='SpectraPlot', optype='other') +# opObj22.addParameter(name='id', value='5', format='int') +# opObj22.addParameter(name='wintitle', value='Spectra Plot', format='str') +# opObj22.addParameter(name='save', value='1', format='bool') +# opObj22.addParameter(name='figpath', value = pathFigure, format='str') + +#----------------------------------------------------------------------------------- + +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 diff --git a/schainpy/test/amisr_windEstimation.py b/schainpy/test/amisr_windEstimation.py new file mode 100644 index 0000000..ec94d46 --- /dev/null +++ b/schainpy/test/amisr_windEstimation.py @@ -0,0 +1,131 @@ +import os, sys + +path = os.path.split(os.getcwd())[0] +sys.path.append(path) + +from controller import * + +desc = "AMISR Experiment" + +filename = "amisr_reader.xml" + +controllerObj = Project() + +controllerObj.setup(id = '191', name='test01', description=desc) + + +path = os.path.join(os.environ['HOME'],'Development/amisr/data') +path = '/home/soporte/Data/AMISR' +figpath = os.path.join(os.environ['HOME'],'Pictures/amisr') + +pathFigure = '/home/soporte/workspace/Graficos/DBS/amisr/' +xmin = '12.0' +xmax = '14.0' + +readUnitConfObj = controllerObj.addReadUnit(datatype='AMISRReader', + path=path, + startDate='2014/10/21', + endDate='2014/10/21', + startTime='00:00:00', + endTime='23:59:59', + walk=1, + timezone='lt') + +#AMISR Processing Unit +procUnitAMISRBeam0 = controllerObj.addProcUnit(datatype='AMISRProc', inputId=readUnitConfObj.getId()) + +opObj11 = procUnitAMISRBeam0.addOperation(name='PrintInfo', optype='other') + +#Reshaper +opObj11 = procUnitAMISRBeam0.addOperation(name='ProfileToChannels', optype='other') + + +#Beam Selector +#opObj11 = procUnitAMISRBeam0.addOperation(name='BeamSelector', optype='other') +#opObj11.addParameter(name='beam', value='0', format='int') + +#Voltage Processing Unit +procUnitConfObjBeam0 = controllerObj.addProcUnit(datatype='VoltageProc', inputId=procUnitAMISRBeam0.getId()) +#Coherent Integration +opObj11 = procUnitConfObjBeam0.addOperation(name='CohInt', optype='other') +opObj11.addParameter(name='n', value='8', format='int') +#Spectra Unit Processing, getting spectras with nProfiles and nFFTPoints +procUnitConfObjSpectraBeam0 = controllerObj.addProcUnit(datatype='SpectraProc', inputId=procUnitConfObjBeam0.getId()) +procUnitConfObjSpectraBeam0.addParameter(name='nFFTPoints', value=32, format='int') +procUnitConfObjSpectraBeam0.addParameter(name='nProfiles', value=32, format='int') +#RemoveDc +opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='removeDC') + +#Noise Estimation +opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='getNoise') +opObj11.addParameter(name='minHei', value='5', format='float') +opObj11.addParameter(name='maxHei', value='20', format='float') + +#SpectraPlot +opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='SpectraPlot', optype='other') +opObj11.addParameter(name='id', value='100', format='int') +opObj11.addParameter(name='wintitle', value='AMISR Beam 0', format='str') +opObj11.addParameter(name='zmin', value='30', format='int') +opObj11.addParameter(name='zmax', value='80', format='int') + +#RTIPlot +#title0 = 'RTI AMISR Beam 0' +#opObj11 = procUnitConfObjSpectraBeam0.addOperation(name='RTIPlot', optype='other') +#opObj11.addParameter(name='id', value='200', format='int') +#opObj11.addParameter(name='wintitle', value=title0, format='str') +#opObj11.addParameter(name='showprofile', value='0', format='int') +##Setting RTI time using xmin,xmax +#opObj11.addParameter(name='xmin', value='15', format='int') +#opObj11.addParameter(name='xmax', value='23', format='int') +#Setting dB range with zmin, zmax +#opObj11.addParameter(name='zmin', value='45', format='int') +#opObj11.addParameter(name='zmax', value='70', format='int') +#Save RTI +#figfile0 = 'amisr_rti_beam0.png' +#opObj11.addParameter(name='figpath', value=figpath, format='str') +#opObj11.addParameter(name='figfile', value=figfile0, format='str') + +#----------------------------------------------------------------------------------------------- +procUnitConfObj2 = controllerObj.addProcUnit(datatype='ParametersProc', inputId=procUnitConfObjSpectraBeam0 .getId()) +opObj20 = procUnitConfObj2.addOperation(name='GetMoments') + +# opObj21 = procUnitConfObj2.addOperation(name='MomentsPlot', optype='other') +# opObj21.addParameter(name='id', value='3', format='int') +# opObj21.addParameter(name='wintitle', value='Moments Plot', format='str') +# opObj21.addParameter(name='save', value='1', format='bool') +# opObj21.addParameter(name='figpath', value=pathFigure, format='str') +# opObj21.addParameter(name='zmin', value='5', format='int') +# opObj21.addParameter(name='zmax', value='90', format='int') + +opObj22 = procUnitConfObj2.addOperation(name='WindProfiler', optype='other') +opObj22.addParameter(name='technique', value='DBS', format='str') +opObj22.addParameter(name='azimuth', value='51.06', format='float') +opObj22.addParameter(name='correctFactor', value='-1', format='float') +opObj22.addParameter(name='dirCosx', value='9.63770247e-01, 5.93066547e-17, 1, 5.93066547e-17,-9.68583161e-01', format='floatlist') +opObj22.addParameter(name='dirCosy', value=' 0,-0.96858316,0,0.96858316,0', format='floatlist') +# opObj22.addParameter(name='horizontalOnly', value='1', format='bool') +# opObj22.addParameter(name='channelList', value='1,2', format='intlist') + +opObj23 = procUnitConfObj2.addOperation(name='WindProfilerPlot', optype='other') +opObj23.addParameter(name='id', value='4', format='int') +opObj23.addParameter(name='wintitle', value='Wind Profiler', format='str') +opObj23.addParameter(name='save', value='1', format='bool') +opObj23.addParameter(name='figpath', value = pathFigure, format='str') +opObj23.addParameter(name='zmin', value='-10', format='int') +opObj23.addParameter(name='zmax', value='10', format='int') +opObj23.addParameter(name='zmin_ver', value='-80', format='float') +opObj23.addParameter(name='zmax_ver', value='80', format='float') +opObj23.addParameter(name='SNRmin', value='-10', format='int') +opObj23.addParameter(name='SNRmax', value='60', format='int') +opObj23.addParameter(name='SNRthresh', value='0', 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()